Spatial Fast X

Abrupt_reversal

###################################################################################
###################################################################################
###################		DETERMINISTIC MODEL 		###################
###################################################################################
###################################################################################


# This is a one-dimensional stepping stone model considering
# - an abrupt environmental transition (Gradual environmental transition model available in Table 2)
# - a dominance reversal model (alternative model: parallel dominance)


############ PARAMETERS and CONSTANTS ############
# Even integer of demes
H = 100
# Number of generations per run
G = 10000
# Sex-specific migration coefficients between two adjacent patches
mm = 0.1
mf = 0.1
# Sex-specific selection coefficient against the deleterious allele
sm = 0.05
sf = 0.05
# Dominance coefficient of the deleterious allele
h = 0.5



############ GENOTYPE FREQUENCIES ############
## Vectors to hold genotype frequencies in FEMALES
	# autosomal locus
fAA_a = rep(0,H)
fAB_a = rep(0,H)
fBB_a = rep(0,H)
	# X-linked locus
fAA_x = rep(0,H)
fAB_x = rep(0,H)
fBB_x = rep(0,H)

## Vectors to hold genotype frequencies in MALES
	# autosomal locus
gAA_a = rep(0,H)
gAB_a = rep(0,H)
gBB_a = rep(0,H)
	# X-linked locus
gA_x = rep(0,H)
gB_x = rep(0,H)


## INITIAL GENOTYPE FREQUENCIES at the autosomal and and X-linked loci, in females and males
# from the first deme to the point of environmental transition (set at mid-range in this model)
for (i in 1:(H/2) ){
### IN FEMALES ###
  fAA_a[i] = 0.25
  fAB_a[i] = 0.5
  fBB_a[i] = 0.25
  fAA_x[i] = 0.25
  fAB_x[i] = 0.5
  fBB_x[i] = 0.25
### IN MALES ###
  gAA_a[i] = 0.25
  gAB_a[i] = 0.5
  gBB_a[i] = 0.25
  gA_x[i] = 0.5
  gB_x[i] = 0.5
}

# from mid-range to deme H
for (i in (H/2+1):H){
### IN FEMALES ###
  fAA_a[i] = 0.25
  fAB_a[i] = 0.5
  fBB_a[i] = 0.25
  fAA_x[i] = 0.25
  fAB_x[i] = 0.5
  fBB_x[i] = 0.25
### IN MALES ###
  gAA_a[i] = 0.25
  gAB_a[i] = 0.5
  gBB_a[i] = 0.25
  gA_x[i] = 0.5
  gB_x[i] = 0.5
}



######### 	1. MIGRATION	#########
## Vectors to hold genotype frequencies AFTER SEX-SPECIFIC MIGRATION

### IN FEMALES ###
	# autosomal locus
fAA_a_mig = rep(0,H)
fAB_a_mig = rep(0,H)
fBB_a_mig = rep(0,H)
	# X-linked locus
fAA_x_mig = rep(0,H)
fAB_x_mig = rep(0,H)
fBB_x_mig = rep(0,H)

### IN MALES ###
	# autosomal locus
gAA_a_mig = rep(0,H)
gAB_a_mig = rep(0,H)
gBB_a_mig = rep(0,H)
	# X-linked locus
gA_x_mig = rep(0,H)
gB_x_mig = rep(0,H)


######### 	2. SELECTION	##########
## Vectors to hold genotype frequencies AFTER SEX-SPECIFIC SELECTION

### IN FEMALES ###
	# autosomal locus
fAA_a_sel = rep(0,H)
fAB_a_sel = rep(0,H)
fBB_a_sel = rep(0,H)
	# X-linked locus
fAA_x_sel = rep(0,H)
fAB_x_sel = rep(0,H)
fBB_x_sel = rep(0,H)

### IN MALES ###
	# autosomal locus
gAA_a_sel = rep(0,H)
gAB_a_sel = rep(0,H)
gBB_a_sel = rep(0,H)
	# X-linked locus
gA_x_sel = rep(0,H)
gB_x_sel = rep(0,H)



#### FITNESS ####
## Vectors to hold AVERAGE POPULATIONS FITNESS

### IN FEMALES ###
	# autosomal locus
Wavg_a = rep(0,H)
	# X-linked locus
Wavg_x = rep(0,H)

### IN MALES ###
	# autosomal locus
Vavg_a = rep(0,H)
	# X-linked locus
Vavg_x = rep(0,H)


## Vectors for GENOTYPE FITNESS 

### IN FEMALES ###
	# autosomal locus
WAA_a = rep(0,H)
WAB_a = rep(0,H)
WBB_a = rep(0,H)
	# X-linked locus
WAA_x = rep(0,H)
WAB_x = rep(0,H)
WBB_x = rep(0,H)

### IN MALES ###
	# autosomal locus
VAA_a = rep(0,H)
VAB_a = rep(0,H)
VBB_a = rep(0,H)
	# X-linked locus
VA_x = rep(0,H)
VB_x = rep(0,H) 


## Fitness loops
for (i in 1:(H/2) ) {
### IN FEMALES ###
  WAA_a[i] = 1-sf
  WAB_a[i] = 1-h*sf
  WBB_a[i] = 1
  WAA_x[i] = 1-sf
  WAB_x[i] = 1-h*sf
  WBB_x[i] = 1
### IN MALES ###
  VAA_a[i] = 1-sm
  VAB_a[i] = 1-h*sm
  VBB_a[i] = 1
  VA_x[i] = 1-sm
  VB_x[i] = 1
}

for (i in (H/2+1):H){
### IN FEMALES ###
  WAA_a[i] = 1
  WAB_a[i] = 1-h*sf
  WBB_a[i] = 1-sf
  WAA_x[i] = 1
  WAB_x[i] = 1-h*sf
  WBB_x[i] = 1-sf
### IN MALES ### 
  VAA_a[i] = 1
  VAB_a[i] = 1-h*sm
  VBB_a[i] = 1-sm
  VA_x[i] = 1
  VB_x[i] = 1-sm
}

######### CLINE SLOPE EXTRACTION #########

	# autosomal locus
slope_Auto=rep(0,H)
	# X-linked locus
slope_X=rep(0,H)



##################################################################### 
###################### 	   GENERATION LOOP 	##################### 
#####################################################################
  

### START LOOP ###
for(j in 1:G){
  
## 1. SEX-SPECIFIC MIGRATION
# 1.1 genotype frequencies after migration in patches 2 to H-1
  
  	for(i in 2:(H-1)){
    
  ### IN FEMALES ###
    # autosomal locus
  fAA_a_mig[i] = fAA_a[i]*(1-mf) + (mf/2)*(fAA_a[i-1] + fAA_a[i+1])
  fAB_a_mig[i] = fAB_a[i]*(1-mf) + (mf/2)*(fAB_a[i-1] + fAB_a[i+1])
  fBB_a_mig[i] = fBB_a[i]*(1-mf) + (mf/2)*(fBB_a[i-1] + fBB_a[i+1])
    # X-linked locus
  fAA_x_mig[i] = fAA_x[i]*(1-mf) + (fAA_x[i-1] + fAA_x[i+1])*(mf/2)
  fAB_x_mig[i] = fAB_x[i]*(1-mf) + (fAB_x[i-1] + fAB_x[i+1])*(mf/2) 
  fBB_x_mig[i] = fBB_x[i]*(1-mf) + (fBB_x[i-1] + fBB_x[i+1])*(mf/2)
    
  ### IN MALES ###
    # autosomal locus
  gAA_a_mig[i] = gAA_a[i]*(1-mm) + (gAA_a[i-1] + gAA_a[i+1])*(mm/2)
  gAB_a_mig[i] = gAB_a[i]*(1-mm) + (gAB_a[i-1] + gAB_a[i+1])*(mm/2)
  gBB_a_mig[i] = gBB_a[i]*(1-mm) + (gBB_a[i-1] + gBB_a[i+1])*(mm/2)
    # X-linked locus
  gA_x_mig[i] = gA_x[i]*(1-mm) + (gA_x[i-1] + gA_x[i+1])*(mm/2)
  gB_x_mig[i] = gB_x[i]*(1-mm) + (gB_x[i-1] + gB_x[i+1])*(mm/2) 
  }
 
  
# 1.2 genotype frequencies after migration in patch 1
  
  ### IN FEMALES ###
  	# autosomal locus
  fAA_a_mig[1] = fAA_a[2]*(mf/2) + fAA_a[1]*(1-(mf/2))
  fAB_a_mig[1] = fAB_a[2]*(mf/2) + fAB_a[1]*(1-(mf/2))
  fBB_a_mig[1] = fBB_a[2]*(mf/2) + fBB_a[1]*(1-(mf/2))
  	# X-linked locus
  fAA_x_mig[1] = fAA_x[2]*(mf/2) + fAA_x[1]*(1-(mf/2))
  fAB_x_mig[1] = fAB_x[2]*(mf/2) + fAB_x[1]*(1-(mf/2))
  fBB_x_mig[1] = fBB_x[2]*(mf/2) + fBB_x[1]*(1-(mf/2))
  
  ### IN MALES ###
  	# autosomal locus
  gAA_a_mig[1] = gAA_a[2]*(mm/2) + gAA_a[1]*(1-(mm/2))
  gAB_a_mig[1] = gAB_a[2]*(mm/2) + gAB_a[1]*(1-(mm/2))
  gBB_a_mig[1] = gBB_a[2]*(mm/2) + gBB_a[1]*(1-(mm/2))
  	# X-linked locus
  gA_x_mig[1] = gA_x[2]*(mm/2) + gA_x[1]*(1-(mm/2))
  gB_x_mig[1] = gB_x[2]*(mm/2) + gB_x[1]*(1-(mm/2))
 
  
# 1.3 genotype frequencies after migration in patch H
  
  ### IN FEMALES ###
  	# autosomal locus
  fAA_a_mig[H] = fAA_a[H-1]*(mf/2) + fAA_a[H]*(1-(mf/2))
  fAB_a_mig[H] = fAB_a[H-1]*(mf/2) + fAB_a[H]*(1-(mf/2))
  fBB_a_mig[H] = fBB_a[H-1]*(mf/2) + fBB_a[H]*(1-(mf/2))
  	# X-linked locus
  fAA_x_mig[H] = fAA_x[H-1]*(mf/2) + fAA_x[H]*(1-(mf/2))
  fAB_x_mig[H] = fAB_x[H-1]*(mf/2) + fAB_x[H]*(1-(mf/2))
  fBB_x_mig[H] = fBB_x[H-1]*(mf/2) + fBB_x[H]*(1-(mf/2)) 
  
  ### IN MALES ###
  	# autosomal locus
  gAA_a_mig[H] = gAA_a[H-1]*(mm/2) + gAA_a[H]*(1-(mm/2))
  gAB_a_mig[H] = gAB_a[H-1]*(mm/2) + gAB_a[H]*(1-(mm/2))
  gBB_a_mig[H] = gBB_a[H-1]*(mm/2) + gBB_a[H]*(1-(mm/2))
 	# X-linked locus 
  gA_x_mig[H] = gA_x[H-1]*(mm/2) + gA_x[H]*(1-(mm/2))
  gB_x_mig[H] = gB_x[H-1]*(mm/2) + gB_x[H]*(1-(mm/2))
  
  
  
## 2. SEX-SPECIFIC SELECTION
# 2.1. average population fitness
  
  for(i in 1:H){
  
  ### FEMALE POPULATION ###
    # autosomal locus
  Wavg_a[i] = WAA_a[i]*fAA_a_mig[i] + WAB_a[i]*fAB_a_mig[i] + WBB_a[i]*fBB_a_mig[i]
    # X-linked locus
  Wavg_x[i] = WAA_x[i]*fAA_x_mig[i] + WAB_x[i]*fAB_x_mig[i] + WBB_x[i]*fBB_x_mig[i]
    
  ### MALE POPULATION ###
    # autosomal locus
  Vavg_a[i] = VAA_a[i]*gAA_a_mig[i] + VAB_a[i]*gAB_a_mig[i] + VBB_a[i]*gBB_a_mig[i]
    # X-linked locus
  Vavg_x[i] = VA_x[i]*gA_x_mig[i] + VB_x[i]*gB_x_mig[i]
  }
  
# 2.2. genotype frequencies after selection
  
  for (i in 1:H){  
  
  ### IN FEMALES ###
    # autosomal locus
  fAA_a_sel[i] = (WAA_a[i]*fAA_a_mig[i]) / Wavg_a[i]
  fAB_a_sel[i] = (WAB_a[i]*fAB_a_mig[i]) / Wavg_a[i]
  fBB_a_sel[i] = (WBB_a[i]*fBB_a_mig[i]) / Wavg_a[i]
    # X-linked locus
  fAA_x_sel[i] = (WAA_x[i]*fAA_x_mig[i]) / Wavg_x[i]
  fAB_x_sel[i] = (WAB_x[i]*fAB_x_mig[i]) / Wavg_x[i]
  fBB_x_sel[i] = (WBB_x[i]*fBB_x_mig[i]) / Wavg_x[i]
    
  ### IN MALES ###
    # autosomal locus
  gAA_a_sel[i] = (VAA_a[i]*gAA_a_mig[i]) / Vavg_a[i]
  gAB_a_sel[i] = (VAB_a[i]*gAB_a_mig[i]) / Vavg_a[i]
  gBB_a_sel[i] = (VBB_a[i]*gBB_a_mig[i]) / Vavg_a[i] 
    # X-linked locus 
  gA_x_sel[i] = (VA_x[i]*gA_x_mig[i]) / Vavg_x[i]
  gB_x_sel[i] = (VB_x[i]*gB_x_mig[i]) / Vavg_x[i]
  }
  
# 3. Random mating (no drift)
  
  for(i in 1:H) {
    
  ### IN FEMALES ###
    # autosomal locus
  fAA_a[i] = fAA_a_sel[i]*gAA_a_sel[i] + fAA_a_sel[i]*gAB_a_sel[i]/2 + fAB_a_sel[i]*gAA_a_sel[i]/2 + fAB_a_sel[i]*gAB_a_sel[i]/4
  fBB_a[i] = fBB_a_sel[i]*gBB_a_sel[i] + fBB_a_sel[i]*gAB_a_sel[i]/2 + fAB_a_sel[i]*gBB_a_sel[i]/2 + fAB_a_sel[i]*gAB_a_sel[i]/4
  fAB_a[i] = 1 - fAA_a[i] - fBB_a[i] 
    # X-linked locus
  fAA_x[i] = fAA_x_sel[i]*gA_x_sel[i] + fAB_x_sel[i]*gA_x_sel[i]/2
  fBB_x[i] = fBB_x_sel[i]*gB_x_sel[i] + fAB_x_sel[i]*gB_x_sel[i]/2
  fAB_x[i] = 1 - fAA_x[i] - fBB_x[i]
    
  ### IN MALES ###
    # autosomal locus
  gAA_a[i] = fAA_a[i]
  gBB_a[i] = fBB_a[i]
  gAB_a[i] = fAB_a[i]
    # X-linked locus
  gA_x[i] = fAA_x_sel[i] + fAB_x_sel[i]/2
  gB_x[i] = fBB_x_sel[i] + fAB_x_sel[i]/2
  }

# calculation of autosomal and X-linked cline slopes using changes in frequency of allele A between adjacent patches
  
  for(i in 2:H){
    slope_Auto[i-1] = ((fAA_a[i] + fAB_a[i]/2 + gAA_a[i] +gAB_a[i]/2)/2) - ((fAA_a[i-1] + fAB_a[i-1]/2 + gAA_a[i-1] +gAB_a[i-1]/2)/2)
    slope_X[i-1]= (2*(fAA_x[i] + fAB_x[i]/2)/3 + gA_x[i]/3) - (2*(fAA_x[i-1] + fAB_x[i-1]/2)/3 + gA_x[i-1]/3)
  }
  
} #### END OF THE LOOP




###############################
#######  CALCULATIONS	#######
###############################

## CLINE SLOPES calculated with the REACTION DIFFUSION APPROXIMATION 
  # autosomal locus
approx_auto = sqrt((sm + sf)*(6*h + 5)/(48*(mf + mm)))
approx_auto
  # X-linked locus
approx.X = sqrt((sf*(6*h + 5) + 8*sm)/(24*(2*mf + mm)))
approx.X
  # calculation of the X/Autosomal ratio
approx.X/approx_auto


## CLINE SLOPES calculated with the STEPPING STONE MODEL
# We use the frequency of the allele A at the autosomal locus (and assumed equal sex ratio) of the offspring generation (after random mating)
  # autosomal locus
fA_a = (fAA_a + fAB_a/2 + gAA_a +gAB_a/2) / 2
  # X-linked locus
fA_x = 2*(fAA_x + fAB_x/2)/3 + gA_x/3

# And calculate the difference in allele frequency between the patches H/2 and H/2+1 (patches on either sides of the environmental transition)
  # autosomal locus
slope_a = fA_a[H/2+1] - fA_a[H/2]
slope_a
  # X-linked locus
slope_x = fA_x[H/2+1] - fA_x[H/2]
slope_x
  # calculation of the X/Autosomal ratio
slope_x/slope_a


## We check that the maximum cline slope using the slope_Auto and slope_X vectors match the above results
  # autosomal locus
max.slope_Auto = max(slope_Auto)
max.slope_Auto
  # X-linked locus	
max.slope_X = max(slope_X)
max.slope_X
  # calculation of the X/Autosomal ratio MAXI X/Auto
max.slope_X / max.slope_Auto
	

###########################################
#######    VISUALISING THE CLINES   #######
###########################################
par(mfrow=c(1,1))
  # AUTOSOMAL frequency of the A allele (red line)
plot(1:H,(fAA_a+fAB_a/2+gAA_a+gAB_a/2)/2 ,ylim=c(0,1),type="l",col="red")
par(new=TRUE)
  # X-LINKED frequency of the A allele (blue line)
plot(1:H, 2*(fAA_x + fAB_x/2)/3 + gA_x/3 ,ylim=c(0,1),type="l",col="blue")









###################################################################################
###################################################################################
###################		STOCHASTIC MODEL 		###################
###################################################################################
###################################################################################

# This is a one-dimensional stepping stone model considering
# - an abrupt environmental transition (Gradual environmental transition model available in Table 2)
# - a dominance reversal model (alternative model: parallel dominance)


############ DEPENDENCIES and FUNCTIONS ############

#  Calculating sex-specific migration
migr  <-  function(genotype, H, resMigRate, immMigRate) {
  nonBoundary  <-  genotype[-c(1,H)]
  lPatchIndex  <-  (seq_along(genotype) - 1)[-c(1,H)]
  rPatchIndex  <-  (seq_along(genotype) + 1)[-c(1,H)]
  
  newNonBoundaryFreq  <-  nonBoundary*resMigRate + (genotype[lPatchIndex] + genotype[rPatchIndex])*(immMigRate/2)
  newLBoundaryFreq    <-  genotype[1]*(resMigRate/2) + genotype[2]*(immMigRate/2) 
  newRBoundaryFreq    <-  genotype[H]*(resMigRate/2) + genotype[H-1]*(immMigRate/2) 
  
  newFreq  <-  c(newLBoundaryFreq, newNonBoundaryFreq, newRBoundaryFreq)
  newFreq
}

# drift technique for autosomal genes
pAfterDrift  <-  function(N, homoA, het, homoB) {
  a_drift  <-  rmultinom(1,N,c( homoA, het, homoB )) / N
  a_drift[1]+a_drift[2]/2
}

# drift technique for X-linked genes
pmxAfterDrift  <-  function(N, A, B) {
  M_x_drift  <-  rmultinom(1,N,c( A, B )) / N
  M_x_drift[1]
}

# Calculating average allele frequency between 2 adjacent patches
## at an autosomal locus
pAverageA <- function(pF, pM){
  nonBoundPf   <-  pF[-1]
  nonBoundPm   <-  pM[-1]
  lPatchIndex  <- (seq_along(pF) - 1)[-1] #??
  
  (((nonBoundPf+nonBoundPm)/2 + (pF[lPatchIndex]+pM[lPatchIndex])/2) /2)
}

## at an X-linked locus
pAverageX <- function(pF, pM){
  nonBoundPf   <-  pF[-1]
  nonBoundPm   <-  pM[-1]
  lPatchIndex  <- (seq_along(pF) - 1)[-1] #??
  
  (((2*nonBoundPf+nonBoundPm)/3 + (2*pF[lPatchIndex]+pM[lPatchIndex])/3) /2)
}

# Calculating FST between 2 adjacent patches
## at an autosomal locus
FstA <- function(pF,pM,pAVG){
  nonBoundPf   <-  pF[-1]
  nonBoundPm   <-  pM[-1]
  lPatchIndex  <- (seq_along(pF) - 1)[-1]
  FST<- (((pF[lPatchIndex]+pM[lPatchIndex])/2 - (nonBoundPf+nonBoundPm)/2)^2 / (4*pAVG*(1-pAVG)))
  
  c(FST, max(FST[is.nan(FST) == FALSE]))
}

## at an X-linked locus
FstX <- function(pF,pM,pAVG){
  nonBoundPf   <-  pF[-1]
  nonBoundPm   <-  pM[-1]
  lPatchIndex  <- (seq_along(pF) - 1)[-1]
  FST<- (((2*pF[lPatchIndex]+pM[lPatchIndex])/3 - (2*nonBoundPf+nonBoundPm)/3)^2 / (4*pAVG*(1-pAVG)))
  
  c(FST, max(FST[is.nan(FST) == FALSE]))
}


############ PARAMETERS and CONSTANTS ############

# Number of replicate simulations
runs = 1000
# Even integer of demes
H = 100
# Sex-specific migration coefficients between two adjacent patches
mm = 0.1
mf = 0.1
# Sex-specific selection coefficient against the deleterious allele
sm = 0.05
sf = 0.05
# Dominance coefficient of the deleterious allele
h = 0.5
## sex-specific population size
Nf = 1400
Nm = 9800
## Number of generation G = 2*Ne with Ne = 4*Nf*Nm/(Nf+Nm)
G = 2*(4*Nf*Nm/(Nf+Nm))

## Storage matrix for results
matrixFstA  <-  matrix(0, nrow=runs, ncol=H) #this goes at the top of your code
matrixFstX  <-  matrix(0, nrow=runs, ncol=H) #this goes at the top of your code


###########################################################################
######################### Replicate Simulation Loop #######################
###########################################################################

for(k in 1:runs){
  ######################################
  ######### INITIAL CONDITIONS #########
  ######################################
  
  # AUTOSOMAL Initial allele frequencies in females and males
  pf_a  <- rep(0.5, times=H)
  pm_a  <- rep(0.5, times=H)
  
  
  # X-LINKED Initial allele frequencies in females and males
  pf_x  <-  rep(0.5, times=H)
  pm_x  <-  rep(0.5, times=H)
  
  
  #################################
  ######### RANDOM MATING #########
  #################################
  
  ##  genotype frequencies in offspring
  
  # AUTOSOMAL LOCUS in females
  fAA_a = rep(0,H)
  fAB_a = rep(0,H)
  fBB_a = rep(0,H)
  # AUTOSOMAL LOCUS in males
  gAA_a = rep(0,H)
  gAB_a = rep(0,H)
  gBB_a = rep(0,H)
  
  # X-LINKED LOCUS in females
  fAA_x = rep(0,H)
  fAB_x = rep(0,H)
  fBB_x = rep(0,H)
  # X-LINKED LOCUS in males
  gA_x = rep(0,H)
  gB_x = rep(0,H)
  
  #############################
  ######### MIGRATION #########
  #############################
  
  ## Genotype frequencies after migration
  
  # AUTOSOMAL LOCUS in females
  fAA_a_mig = rep(0,H)
  fAB_a_mig = rep(0,H)
  fBB_a_mig = rep(0,H)
  # AUTOSOMAL LOCUS in males
  gAA_a_mig = rep(0,H)
  gAB_a_mig = rep(0,H)
  gBB_a_mig = rep(0,H)
  
  # X-LINKED LOCUS in females
  fAA_x_mig = rep(0,H)
  fAB_x_mig = rep(0,H)
  fBB_x_mig = rep(0,H)
  # X-LINKED LOCUS  in males
  gA_x_mig = rep(0,H)
  gB_x_mig = rep(0,H)
  
  #############################
  ######### SELECTION #########
  #############################
  
  ## Vectors for genotype frequencies after selection
  
  # AUTOSOMAL LOCUS in females
  fAA_a_sel = rep(0,H)
  fAB_a_sel = rep(0,H)
  fBB_a_sel = rep(0,H)
  # AUTOSOMAL LOCUS in males
  gAA_a_sel = rep(0,H)
  gAB_a_sel = rep(0,H)
  gBB_a_sel = rep(0,H)
  
  # X-LINKED LOCUS in females
  fAA_x_sel = rep(0,H)
  fAB_x_sel = rep(0,H)
  fBB_x_sel = rep(0,H)
  # X-LINKED LOCUS in males
  gA_x_sel = rep(0,H)
  gB_x_sel = rep(0,H)
  
  ## Vectors for average population fitness per patch
  # AUTOSOMAL LOCUS in females
  W_a = rep (0,H)
  # AUTOSOMAL LOCUS in males
  V_a = rep (0,H)
  
  # X-LINKED LOCUS in females
  W_x = rep (0,H)
  # X-LINKED LOCUS in males
  V_x = rep (0,H)
  
  ## fitness vectors for each genotype at autosomal and X-linked loci
  
  # AUTOSOMAL LOCUS in females
  wAA_a  <-  c(rep(1-sf, times=H/2), rep(1,times=H/2))
  wAB_a  <-  c(rep(1-h*sf, times=H))
  wBB_a  <-  c(rep(1, times=H/2), rep(1-sf,times=H/2))
  # AUTOSOMAL LOCUS in males 
  vAA_a  <-  c(rep(1-sm, times=H/2), rep(1,times=H/2))
  vAB_a  <-  c(rep(1-h*sm, times=H))
  vBB_a  <-  c(rep(1, times=H/2), rep(1-sm,times=H/2))
  
  # X-LINKED LOCUS IN FEMALES
  wAA_x  <-  c(rep(1-sf, times=H/2), rep(1,times=H/2))
  wAB_x  <-  c(rep(1-h*sf, times=H))
  wBB_x  <-  c(rep(1, times=H/2), rep(1-sf,times=H/2))
  # X-LINKED LOCUS IN MALES
  vA_x  <-  c(rep(1-sm, times=H/2), rep(1,times=H/2))
  vB_x  <-  c(rep(1, times=H/2), rep(1-sm,times=H/2))  
  
  
  ##################################################################
  ####################### GENERATION LOOP ##########################
  ##################################################################
  
  
  ### START LOOP ###
  for(j in 1:G){
    
    # Random mating #
    
    # AUTOSOMAL LOCUS - genotype frequencies in female offsprings
    fAA_a = pf_a*pm_a
    fAB_a = pf_a*(1-pm_a) + pm_a*(1-pf_a)
    fBB_a = (1-pf_a)*(1-pm_a)
    # AUTOSOMAL LOCUS - genotype frequencies in male offsprings
    gAA_a = pf_a*pm_a
    gAB_a = pf_a*(1-pm_a)+pm_a*(1-pf_a)
    gBB_a = (1-pf_a)*(1-pm_a)
    
    # X-LINKED LOCUS - genotype frequencies in female offsprings
    fAA_x = pf_x*pm_x
    fAB_x = pf_x*(1-pm_x) + pm_x*(1-pf_x)
    fBB_x = (1-pf_x)*(1-pm_x)
    # X-LINKED LOCUS - genotype frequencies in male offsprings
    gA_x = pf_x
    gB_x = (1-pf_x)
    
    
    # Sex-specific migration #
    ## genotype frequencies after migration in patches 2 to H-1
    
    # AUTOSOMAL LOCUS in females
    fAA_a_mig = migr(genotype=fAA_a, H=H, resMigRate=(1-mf), immMigRate=mf)
    fAB_a_mig = migr(genotype=fAB_a, H=H, resMigRate=(1-mf), immMigRate=mf)      
    fBB_a_mig = migr(genotype=fBB_a, H=H, resMigRate=(1-mf), immMigRate=mf)
    gAA_a_mig = migr(genotype=gAA_a, H=H, resMigRate=(1-mm), immMigRate=mm)
    gAB_a_mig = migr(genotype=gAB_a, H=H, resMigRate=(1-mm), immMigRate=mm)      
    gBB_a_mig = migr(genotype=gBB_a, H=H, resMigRate=(1-mm), immMigRate=mm)
    
    # X-LINKED LOCUS in females
    fAA_x_mig = migr(genotype=fAA_x, H=H, resMigRate=(1-mf), immMigRate=mf)
    fAB_x_mig = migr(genotype=fAB_x, H=H, resMigRate=(1-mf), immMigRate=mf)      
    fBB_x_mig = migr(genotype=fBB_x, H=H, resMigRate=(1-mf), immMigRate=mf)
    
    # X-LINKED LOCUS in males
    gA_x_mig = migr(genotype=gA_x, H=H, resMigRate=(1-mm), immMigRate=mm)
    gB_x_mig = migr(genotype=gB_x, H=H, resMigRate=(1-mm), immMigRate=mm)
    
    
  
    # Sex-specific selection #
    
    ## Average fitness
    # AUTOSOMAL LOCUS in females
    W_a = fAA_a_mig*wAA_a + fAB_a_mig*wAB_a + fBB_a_mig*wBB_a
    # AUTOSOMAL LOCUS in males
    V_a = gAA_a_mig*vAA_a + gAB_a_mig*vAB_a + gBB_a_mig*vBB_a
    
    # X-LINKED LOCUS in females
    W_x = fAA_x_mig*wAA_x + fAB_x_mig*wAB_x + fBB_x_mig*wBB_x
    # X-LINKED LOCUS in males
    V_x= gA_x_mig*vA_x + gB_x_mig*vB_x
    
    
    ## genotype frequencies after selection      
    # AUTOSOMAL LOCUS in females
    fAA_a_sel = fAA_a_mig*wAA_a / W_a 
    fAB_a_sel = fAB_a_mig*wAB_a / W_a
    fBB_a_sel = fBB_a_mig*wBB_a / W_a 
    # AUTOSOMAL LOCUS in males
    gAA_a_sel = gAA_a_mig*vAA_a / V_a
    gAB_a_sel = gAB_a_mig*vAB_a / V_a
    gBB_a_sel = gBB_a_mig*vBB_a / V_a
    
    # X-LINKED LOCUS in females
    fAA_x_sel = fAA_x_mig*wAA_x / W_x 
    fAB_x_sel = fAB_x_mig*wAB_x / W_x
    fBB_x_sel = fBB_x_mig*wBB_x / W_x 
    # X-LINKED LOCUS in males
    gA_x_sel = gA_x_mig*vA_x / V_x
    gB_x_sel = gB_x_mig*vB_x / V_x
    
    
    # Drift in adults #
    
    for(i in 1:H) {
      
      # AUTOSOMAL LOCUS in females
      pf_a[i]  <-  pAfterDrift(Nf, fAA_a_sel[i], fAB_a_sel[i], fBB_a_sel[i])
      # AUTOSOMAL LOCUS in males
      pm_a[i]  <-  pAfterDrift(Nm, gAA_a_sel[i], gAB_a_sel[i], gBB_a_sel[i])
      # X-LINKED LOCUS in females
      pf_x[i]  <-  pAfterDrift(Nf, fAA_x_sel[i], fAB_x_sel[i], fBB_x_sel[i])
      # X-LINKED LOCUS in males
      pm_x[i]  <-  pmxAfterDrift(Nm, gA_x_sel[i], gB_x_sel[i])
    }
    
    
    
  } #### END GENERATION LOOP ####
  
  # Calculate average allele frequency differences between adjacent patches 
  # at the end of the run
  
  # AUTOSOMAL LOCUS
  p_a.avg <- pAverageA (pF=pf_a, pM=pm_a)
  # X-LINKED LOCUS
  p_x.avg <- pAverageX (pF=pf_x, pM=pm_x)

  
  
  
  # Calculate Fst between adjacent patches #
  # AUTOSOMAL LOCUS
  matrixFstA[k,]  <-  FstA(pF=pf_a, pM=pm_a, pAVG=p_a.avg)
  # X-LINKED LOCUS
  matrixFstX[k,]  <-  FstX(pF=pf_x, pM=pm_x, pAVG=p_x.avg)
  
  
} #### END REPLICATE SIMULATION LOOP ####

matrixFstA
matrixFstX