https://github.com/Schumerlab/hybridization_review
Raw File
Tip revision: 69a398b89365cc069c6856d990c2b74293b52486 authored by schumerm on 02 June 2021, 22:34:41 UTC
add R models
Tip revision: 69a398b
two_locus_DMI_model.R

#####DMI

#Simulates change in allele frequency as a result of two locus selection based on Karlin 1975, Lewontin & Kojima

# % parent 1 ancestry
M=0.25
B=1-M

#set number of generations to simulate
numreps=2000

# initial haplotype frequencies
x<-c((1-B)*M, (1-B)*(1-M), B*M, B*(1-M))
ancestry_track1<-{}
ancestry_track2<-{}

# Set fitness and dominance
s1<-0.1
s2<-0
hbc=0
ha=0.5

#fitness function for the diagonal is multiplicative here for coevolving incompatibilities but can be set differently
diagonal=(1-((1-ha)*hbc*s1))*(1-(ha*(1-hbc)*s2))

# fitness table
# w11 and w44 are conspecific homozygotes; w22 and w33 are heterospecific homozygotes, and so on

w<-function(s){matrix(c(
    1,       1-(1-hbc)*s2,         1-(1-ha)*s1,     diagonal,
    1-(1-hbc)*s2,       1-s2,      diagonal,     1-ha*s2,
    1-(1-ha)*s1,  diagonal,    1-s1,     1-hbc*s1,
    diagonal,  1-ha*s2,    1-hbc*s1,     1
    ), nrow=4,ncol=4)}
    
    
# intitialise and store haplotype frequency over generations
X<-x
d <- x[1]*x[4] - x[2]*x[3]
w_mean<-0
for (j in 1:4){
	for (k in 1:4){
		w_mean = w_mean + as.numeric(x[j]*x[k]*w(s)[j,k])
	}
}
w_bar<-w_mean


# start simulation here
for (z in 1:numreps){ # number of generations
    # disequilibrium
	D = x[1]*x[4] - x[2]*x[3]; D
	d<-c(d,as.numeric(D)) # tracks D over time
    
    #from allele frequencies assuming random mating
    f_M1=x[1]+x[2]
    f_M2=x[1]+x[3]
    
    #store allele frequencies
    ancestry_track1<-c(ancestry_track1,f_M1)
    ancestry_track2<-c(ancestry_track2,f_M2)
    
    # marginal fitnesses
	wstar<-{}
	for (i in 1:4){
        wstar<-c(wstar, x[1]*w(s)[i,1] + x[2]*w(s)[i,2] + x[3]*w(s)[i,3] + x[4]*w(s)[i,4])
	}
    # mean fitness
	w_mean<-0
	for (j in 1:4){
		for (k in 1:4){
			w_mean = w_mean + as.numeric(x[j]*x[k]*w(s)[j,k])
		}
	}
	w_bar<-c(w_bar, w_mean)
    # new haplotype frequencies
	y<-{}
	y<-c(y,as.numeric(round(x[1]*(wstar[1]/w_mean)-(0.5*D*w(s)[1,4]/w_mean),10)))
	y<-c(y,as.numeric(round(x[2]*(wstar[2]/w_mean)+(0.5*D*w(s)[1,4]/w_mean),10)))
	y<-c(y,as.numeric(round(x[3]*(wstar[3]/w_mean)+(0.5*D*w(s)[1,4]/w_mean),10)))
	y<-c(y,as.numeric(round(x[4]*(wstar[4]/w_mean)-(0.5*D*w(s)[1,4]/w_mean),10)))
	x<-y
	X<-rbind(X,as.numeric(y)) # store information
} # for numreps


plot(ancestry_track1,ylim=c(0,0.6),col="darkblue",pch=20,cex.lab=1.4,cex.axis=1.2,ylab="Minor parent ancestry",xlab="Generations since admixture",type="l",lwd=3,xlim=c(0,2000))

points(ancestry_track2,col="darkblue",pch=20,cex.lab=1.3,cex.axis=1.2,ylab="Proportion Parent1",xlab="Generations",type="l",lwd=3,lty=2)


M=0.4
B=1-M

#set number of generations to simulate
numreps=500

# initial haplotype frequencies
x<-c((1-B)*M, (1-B)*(1-M), B*M, B*(1-M))
ancestry_track1<-{}
ancestry_track2<-{}

# Set fitness and dominance
s1<-0.1
s2<-0.1
hbc=0.5
ha=0.5

#fitness function for the diagonal is multiplicative here for coevolving incompatibilities but can be set differently
diagonal=(1-((1-ha)*hbc*s1))*(1-(ha*(1-hbc)*s2))

# fitness table
# w11 and w44 are conspecific homozygotes; w22 and w33 are heterospecific homozygotes, and so on

w<-function(s){matrix(c(
  1,       1-(1-hbc)*s2,         1-(1-ha)*s1,     diagonal,
  1-(1-hbc)*s2,       1-s2,      diagonal,     1-ha*s2,
  1-(1-ha)*s1,  diagonal,    1-s1,     1-hbc*s1,
  diagonal,  1-ha*s2,    1-hbc*s1,     1
), nrow=4,ncol=4)}


# intitialise and store haplotype frequency over generations
X<-x
d <- x[1]*x[4] - x[2]*x[3]
w_mean<-0
for (j in 1:4){
  for (k in 1:4){
    w_mean = w_mean + as.numeric(x[j]*x[k]*w(s)[j,k])
  }
}
w_bar<-w_mean


# start simulation here
for (z in 1:numreps){ # number of generations
  # disequilibrium
  D = x[1]*x[4] - x[2]*x[3]; D
  d<-c(d,as.numeric(D)) # tracks D over time
  
  #from allele frequencies assuming random mating
  f_M1=x[1]+x[2]
  f_M2=x[1]+x[3]
  
  #store allele frequencies
  ancestry_track1<-c(ancestry_track1,f_M1)
  ancestry_track2<-c(ancestry_track2,f_M2)
  
  # marginal fitnesses
  wstar<-{}
  for (i in 1:4){
    wstar<-c(wstar, x[1]*w(s)[i,1] + x[2]*w(s)[i,2] + x[3]*w(s)[i,3] + x[4]*w(s)[i,4])
  }
  # mean fitness
  w_mean<-0
  for (j in 1:4){
    for (k in 1:4){
      w_mean = w_mean + as.numeric(x[j]*x[k]*w(s)[j,k])
    }
  }
  w_bar<-c(w_bar, w_mean)
  # new haplotype frequencies
  y<-{}
  y<-c(y,as.numeric(round(x[1]*(wstar[1]/w_mean)-(0.5*D*w(s)[1,4]/w_mean),10)))
  y<-c(y,as.numeric(round(x[2]*(wstar[2]/w_mean)+(0.5*D*w(s)[1,4]/w_mean),10)))
  y<-c(y,as.numeric(round(x[3]*(wstar[3]/w_mean)+(0.5*D*w(s)[1,4]/w_mean),10)))
  y<-c(y,as.numeric(round(x[4]*(wstar[4]/w_mean)-(0.5*D*w(s)[1,4]/w_mean),10)))
  x<-y
  X<-rbind(X,as.numeric(y)) # store information
} # for numreps

plot(ancestry_track1,ylim=c(0,0.6),col="darkblue",pch=20,cex.lab=1.4,cex.axis=1.2,ylab="Minor parent ancestry",xlab="Generations since admixture",type="l",lwd=3,xlim=c(0,500))

points(ancestry_track2,col="darkblue",pch=20,cex.lab=1.3,cex.axis=1.2,ylab="Proportion Parent1",xlab="Generations",type="l",lwd=3,lty=2)


back to top