library(normalp) #parameters loi_mig<-"norm" #choose your kernel among neg_bino (=negative binomial) pois (=poisson) bino (=binomial) norm (=gaussian) mixt_norm (=mixture of gaussian) or exp_power (=exponential power) S<-700 # number of demes N1<-c(1000,200,400,600,800) #number of individuals per deme on the right of the contact point N2<-1000 #number of individuals per dem on the left of the contact point r<-0.5 #recomnination rate t<-5 #number of simulated generations+1 mu<-c(400,20,1,3,10,30,100) #mean number of steps for discrete dispersal kernels (bino, pois, neg_bino) (1-10) beta<-c(2,6,2.8,1.4,1) #shape parameter of the exponential power kernel poids<-c(1,0.843,0.8,0.856,0.929) # weight of the first Gaussian in the Gaussian mixture mu1<-c(400,300,200,100,50) # variance of the first Gaussian law in the Gaussian mixture mu2<-c(0,936.94,1200,2183.3,4979.6) # variance of the second Gaussian law in the Gaussian mixture k<-c(1000,3,1,0.3,0.1) # shape parameter of the negative binomial (0.2-5) #Initialisation gentot<-array(0,c(4,4,S,t)) # array with gamete number in each deme at t (4 dimensions) envtot<-array(0,S) freqA1<-array(0,c(t,S)) # array with frequency of A1 in each deme at t freqB1<-array(0,c(t,S)) # array with frequency of B1 in each deme at t DL<-array(0,c(t,S)) # array of linkage disequilibrium in each deme at t for (mui in 1:1) { #beginning loop on mu for (betai in 1:1) { #beginning loop on beta for (z in 1:1) { #beginning loop on N1 N<-c(rep(N1[z],S/2),rep(N2,S/2)) for (simul in 1:5) { #beginning loop on number of simulations #p1<-0.6+0.4*runif(1,0,1) #frequency of A1 on the left before contact #q1<-0.6+0.4*runif(1,0,1) #frequency of B1 on the left before contact #p2<-0.4-0.4*runif(1,0,1) #frequency of A1 on the right before contact #q2<-0.4-0.4*runif(1,0,1) #frequency of B1 on the right before contact p1<-1 #frequency of A1 on the left before contact q1<-1 #frequency of B1 on the left before contact p2<-0 #frequency of A1 on the right before contact q2<-0 #frequency of B1 on the right before contact #initial step on each side onf the contact point = generation 1 env<-array(0,c(4,4,S)) for (i in 1:(S/2)) { env[1,1,i]<-round(N1[z]*p1*q1 * p1*q1) env[1,2,i]<-round(N1[z]*p1*q1 * p1*(1-q1)) env[1,3,i]<-round(N1[z]*p1*q1 * (1-p1)*q1) env[1,4,i]<-round(N1[z]*p1*q1 * (1-p1)*(1-q1)) env[2,1,i]<-round(N1[z]*p1*(1-q1) * p1*q1) env[2,2,i]<-round(N1[z]*p1*(1-q1) * p1*(1-q1)) env[2,3,i]<-round(N1[z]*p1*(1-q1) * (1-p1)*q1) env[2,4,i]<-round(N1[z]*p1*(1-q1) * (1-p1)*(1-q1)) env[3,1,i]<-round(N1[z]*(1-p1)*q1 * p1*q1) env[3,2,i]<-round(N1[z]*(1-p1)*q1 * p1*(1-q1)) env[3,3,i]<-round(N1[z]*(1-p1)*q1 * (1-p1)*q1) env[3,4,i]<-round(N1[z]*(1-p1)*q1 * (1-p1)*(1-q1)) env[4,1,i]<-round(N1[z]*(1-p1)*(1-q1) * p1*q1) env[4,2,i]<-round(N1[z]*(1-p1)*(1-q1) * p1*(1-q1)) env[4,3,i]<-round(N1[z]*(1-p1)*(1-q1) * (1-p1)*q1) env[4,4,i]<-round(N1[z]*(1-p1)*(1-q1) * (1-p1)*(1-q1)) } for (i in ((S/2)+1):S) { env[1,1,i]<-round(N2*p2*q2 * p2*q2) env[1,2,i]<-round(N2*p2*q2 * p2*(1-q2)) env[1,3,i]<-round(N2*p2*q2 * (1-p2)*q2) env[1,4,i]<-round(N2*p2*q2 * (1-p2)*(1-q2)) env[2,1,i]<-round(N2*p2*(1-q2) * p2*q2) env[2,2,i]<-round(N2*p2*(1-q2) * p2*(1-q2)) env[2,3,i]<-round(N2*p2*(1-q2) * (1-p2)*q2) env[2,4,i]<-round(N2*p2*(1-q2) * (1-p2)*(1-q2)) env[3,1,i]<-round(N2*(1-p2)*q2 * p2*q2) env[3,2,i]<-round(N2*(1-p2)*q2 * p2*(1-q2)) env[3,3,i]<-round(N2*(1-p2)*q2 * (1-p2)*q2) env[3,4,i]<-round(N2*(1-p2)*q2 * (1-p2)*(1-q2)) env[4,1,i]<-round(N2*(1-p2)*(1-q2) * p2*q2) env[4,2,i]<-round(N2*(1-p2)*(1-q2) * p2*(1-q2)) env[4,3,i]<-round(N2*(1-p2)*(1-q2) * (1-p2)*q2) env[4,4,i]<-round(N2*(1-p2)*(1-q2) * (1-p2)*(1-q2)) } generation<-1 for (i in 1:S) {freqA1[generation,i]<-((sum(env[1:2,,i])+sum(env[,1:2,i]))/2)/sum(env[,,i]) #freq A1 freqB1[generation,i]<-((sum(env[1,,i])+sum(env[3,,i])+sum(env[,1,i])+sum(env[,3,i]))/2)/sum(env[,,i]) #freq B1 DL[generation,i]<-((sum(env[1,,i])+sum(env[,1,i]))/2)/sum(env[,,i])-freqA1[generation,i]*freqB1[generation,i]} #freq (A1B1) - freqA1*freqB1 #Next generations while (generation0){ freq<-env[,,i]/sum(env[,,i]) freq_gam_parents<-rep(NA,4) freq_gam_parents[1]<-(1-r)/2*(sum(freq[1,])+sum(freq[,1]))+r/2*(sum(freq[1,1:3])+sum(freq[1:3,1])+freq[2,3]+freq[3,2]) #calcul pour chaque type de gamète les différentes combinaisons dont il est issu (cas par cas) freq_gam_parents[2]<-(1-r)/2*(sum(freq[2,])+sum(freq[,2]))+r/2*(freq[1,2]+freq[1,4]+freq[2,1]+2*freq[2,2]+freq[2,4]+freq[4,1]+freq[4,2]) freq_gam_parents[3]<-(1-r)/2*(sum(freq[3,])+sum(freq[,3]))+r/2*(sum(freq[1,3:4])+freq[3,1]+2*freq[3,3]+freq[3,4]+freq[4,1]+freq[4,3]) freq_gam_parents[4]<-(1-r)/2*(sum(freq[4,])+sum(freq[,4]))+r/2*(freq[2,3]+freq[3,2]+sum(freq[2:4,4])+sum(freq[4,2:4])) for (j in 1:4) { offspring[j,j]<-round(freq_gam_parents[j]*freq_gam_parents[j]*N[i])} for (j in 1:3){ for (l in (j+1):4) {offspring[j,l]<-round(2*freq_gam_parents[j]*freq_gam_parents[l]*N[i])}} #DISPERSAL for (j in 1:4) { for (l in j:4) { if (offspring[j,l]>0) { if (loi_mig=="neg_bino") {Nsteps<-rnbinom(offspring[j,l],k[ki],(k[ki]/(k[ki]+mu[mui])))} if (loi_mig=="bino") {Nsteps<-rep(mu[mui],offspring[j,l])} if (loi_mig=="pois") {Nsteps<-rpois(offspring[j,l],mu[mui])} if (loi_mig=="norm") {Nsteps<-rep(0,offspring[j,l])} if (loi_mig=="mixt_norm") {Nsteps<-rep(0,offspring[j,l])} if (loi_mig=="exp_power") {Nsteps<-rep(0,offspring[j,l])} for (m in 1:length(Nsteps)) { if ( (loi_mig=='neg_bino')|(loi_mig=='bino')|(loi_mig=='pois') ) {steps<-rbinom(Nsteps[m],1,0.5) steps_pos<-sum(steps) steps_neg<-Nsteps[m]-steps_pos move<-steps_pos-steps_neg } if (loi_mig=="norm") {move<-ifelse(runif(1,0,1)<=0.5,floor(rnorm(1,0,mu[mui]^0.5)),ceiling(rnorm(1,0,mu[mui]^0.5)))} if (loi_mig=="exp_power") {move<-ifelse(runif(1,0,1)<=0.5,floor(rnormp(1,0,mu[mui]^0.5,beta[betai])),ceiling(rnormp(1,0,mu[mui]^0.5,beta[betai])))} if (loi_mig=='mixt_norm') {move<-ifelse(runif(1,0,1)<=poids[mui],ifelse(runif(1,0,1)<=0.5,floor(rnorm(1,0,mu1[mui]^0.5)),ceiling(rnorm(1,0,mu1[mui]^0.5))),ifelse(runif(1,0,1)<=0.5,floor(rnorm(1,0,mu2[mui]^0.5)),ceiling(rnorm(1,0,mu2[mui]^0.5))))} if ((i+move)>0 && (i+move)<(S+1)) {gen[j,l,(i+move)]<-gen[j,l,(i+move)]+1} if ((i+move)<1) {gen[j,l,max(1,-(i+move))]<-gen[j,l,max(1,-(i+move))]+1} if ((i+move)>S) {gen[j,l,(S-((i+move)-S))]<-gen[j,l,(S-((i+move)-S))]+1} }} #fin boucles if (offspring[j,l]>0) et for sur m }} #fin boucles for sur j,l } #fin boucle if (sum(env[,,i]>0)) envtot[i]<-sum(env[,,i]) } #fin boucle for sur i gentot[,,,generation]<-gen generation<-generation+1 env<-gen for (i in 1:S) {freqA1[generation,i]<-((sum(gen[1:2,,i])+sum(gen[,1:2,i]))/2)/sum(gen[,,i]) #freq A1 freqB1[generation,i]<-((sum(gen[1,,i])+sum(gen[3,,i])+sum(gen[,1,i])+sum(gen[,3,i]))/2)/sum(gen[,,i]) #freq B1 DL[generation,i]<-((sum(gen[1,,i])+sum(gen[,1,i]))/2)/sum(gen[,,i])-freqA1[generation,i]*freqB1[generation,i]} #freq (A1B1) - freqA1*freqB1 } #end of loop on the generation time t #Plot of the frequency clines and of the linkage disequilibrium during five generations x11() plot(seq(1,S,1),seq(1,S,1)/S,col="white") lines(seq(1,S,1),freqA1[1,],col="blue") lines(seq(1,S,1),freqA1[3,],col="blue") lines(seq(1,S,1),freqA1[5,],col="blue") lines(seq(1,S,1),freqB1[1,],col="green") lines(seq(1,S,1),freqB1[3,],col="green") lines(seq(1,S,1),freqB1[5,],col="green") lines(seq(1,S,1),DL[1,],col="green") lines(seq(1,S,1),DL[3,],col="green") lines(seq(1,S,1),DL[5,],col="green") #writing of the results in 3 files write.table(freqA1,paste("freqA1_","file_out.csv"),sep=";") write.table(freqB1,paste("freqB1_","file_out.csv"),sep=";") write.table(DL,paste("DL_","file_out.csv"),sep=";") } #end of the loop on the number of simulations }}} #ends of the loops on mui, betai et N1z