Mutation<-function(haplotype, U, add.effect) { Nb_mut1<-rpois(1,U) Nb_mut2<-rpois(1,U) Nb_mut3<-rpois(1,U) Nb_mut4<-rpois(1,U) Nb_mut5<-rpois(1,U) Nb_mut6<-rpois(1,U) Nb_mut7<-rpois(1,U) Nb_mut8<-rpois(1,U) Nb_mut9<-rpois(1,U) Nb_mut10<-rpois(1,U) Nb_mut = Nb_mut1 + Nb_mut2 + Nb_mut3 + Nb_mut4 + Nb_mut5 + Nb_mut6 + Nb_mut7 + Nb_mut8 + Nb_mut9 + Nb_mut10 position_mut<-c(sample(1:50, Nb_mut1, replace = F),sample(51:100, Nb_mut2, replace = F),sample(101:150, Nb_mut3, replace = F),sample(151:200, Nb_mut4, replace = F),sample(201:250, Nb_mut5, replace = F),sample(251:300, Nb_mut6, replace = F),sample(301:350, Nb_mut7, replace = F),sample(351:400, Nb_mut8, replace = F),sample(401:450, Nb_mut9, replace = F),sample(451:500, Nb_mut10, replace = F)) if(Nb_mut != 0){ for(i in 1:length(position_mut)){ if(haplotype[position_mut[i]] != 0){haplotype[position_mut[i]] = 0} else{haplotype[position_mut[i]] = add.effect[position_mut[i]]} } } return(haplotype) } Selection <-function(Npop, dist.phenotype,om_2, selfing.rate) { selected.parent = 0 list.parent = c(NULL) repeat{ fit.max = max(exp(-((dist.phenotype)^2)/(2*om_2))) rand.repro1 = sample(1:Npop, 1) fit.samp = exp(-((dist.phenotype[rand.repro1])^2)/(2*om_2)) if( (fit.samp/fit.max) >= runif(1, 0, 1) ){ selected.parent = selected.parent + 1 list.parent[selected.parent]=rand.repro1 # Selfing or outcrossing if( selfing.rate < runif(1, 0, 1) ){ # Outcrossing repeat{ rand.repro2 = sample(1:Npop, 1) fit.samp = exp(-((dist.phenotype[rand.repro2])^2)/(2*om_2)) if( rand.repro1 != rand.repro2 && (fit.samp/fit.max) >= runif(1, 0, 1) ){ selected.parent = selected.parent + 1 list.parent[selected.parent]=rand.repro2 break } } } else{ # Selfing selected.parent = selected.parent + 1 list.parent[selected.parent]=rand.repro1 } } if(selected.parent == 2*Npop) break } return(list.parent) } Genome.offspring <-function(list.parent, L, Npop, U, genome, add.effect) { for(n in 1:Npop){ haplo.temp1 = c(NULL) haplo.temp2 = c(NULL) for(i in 1:L){ add_loc1 = c(genome[2*list.parent[2*n], i], genome[2*list.parent[2*n] - 1, i]) haplo.temp1[i] = sample(add_loc1, 1) add_loc2 = c(genome[2*list.parent[2*n - 1], i], genome[2*list.parent[2*n - 1] - 1, i]) haplo.temp2[i] = sample(add_loc2, 1) } haplo.temp.mut1 = Mutation(haplo.temp1, U, add.effect) haplo.temp.mut2 = Mutation(haplo.temp2, U, add.effect) if(n == 1){genotype.off = haplo.temp.mut1 genotype.off = rbind(genotype.off, haplo.temp.mut2) } else{genotype.off = rbind(genotype.off, haplo.temp.mut1, haplo.temp.mut2)} } return(genotype.off) } dist.phenotype.offspring <-function(Genome,add.effect, dom.effect,opt, opt.foc, Npop) { geno.off = c(NULL) for(i in 1:Npop){ dist_opt = 0 pheno.temp = 0 for(j in 1:50){ if(Genome[2*i,j] == Genome[2*i - 1,j]){pheno.temp = pheno.temp + 2*Genome[2*i,j]} else{pheno.temp = pheno.temp + dom.effect[j]} } dist_opt = dist_opt + ((pheno.temp + rnorm(1, 0, 0.1) - opt.foc)^2) pheno.temp = 0 for(j in 51:100){ if(Genome[2*i,j] == Genome[2*i - 1,j]){pheno.temp = pheno.temp + 2*Genome[2*i,j]} else{pheno.temp = pheno.temp + dom.effect[j]} } dist_opt = dist_opt + ((pheno.temp + rnorm(1, 0, 0.1) - opt)^2) pheno.temp = 0 for(j in 101:150){ if(Genome[2*i,j] == Genome[2*i - 1,j]){pheno.temp = pheno.temp + 2*Genome[2*i,j]} else{pheno.temp = pheno.temp + dom.effect[j]} } dist_opt = dist_opt + ((pheno.temp + rnorm(1, 0, 0.1) - opt)^2) pheno.temp = 0 for(j in 151:200){ if(Genome[2*i,j] == Genome[2*i - 1,j]){pheno.temp = pheno.temp + 2*Genome[2*i,j]} else{pheno.temp = pheno.temp + dom.effect[j]} } dist_opt = dist_opt + ((pheno.temp + rnorm(1, 0, 0.1) - opt)^2) pheno.temp = 0 for(j in 201:250){ if(Genome[2*i,j] == Genome[2*i - 1,j]){pheno.temp = pheno.temp + 2*Genome[2*i,j]} else{pheno.temp = pheno.temp + dom.effect[j]} } dist_opt = dist_opt + ((pheno.temp + rnorm(1, 0, 0.1) - opt)^2) pheno.temp = 0 for(j in 251:300){ if(Genome[2*i,j] == Genome[2*i - 1,j]){pheno.temp = pheno.temp + 2*Genome[2*i,j]} else{pheno.temp = pheno.temp + dom.effect[j]} } dist_opt = dist_opt + ((pheno.temp + rnorm(1, 0, 0.1) - opt)^2) pheno.temp = 0 for(j in 301:350){ if(Genome[2*i,j] == Genome[2*i - 1,j]){pheno.temp = pheno.temp + 2*Genome[2*i,j]} else{pheno.temp = pheno.temp + dom.effect[j]} } dist_opt = dist_opt + ((pheno.temp + rnorm(1, 0, 0.1) - opt)^2) pheno.temp = 0 for(j in 351:400){ if(Genome[2*i,j] == Genome[2*i - 1,j]){pheno.temp = pheno.temp + 2*Genome[2*i,j]} else{pheno.temp = pheno.temp + dom.effect[j]} } dist_opt = dist_opt + ((pheno.temp + rnorm(1, 0, 0.1) - opt)^2) pheno.temp = 0 for(j in 401:450){ if(Genome[2*i,j] == Genome[2*i - 1,j]){pheno.temp = pheno.temp + 2*Genome[2*i,j]} else{pheno.temp = pheno.temp + dom.effect[j]} } dist_opt = dist_opt + ((pheno.temp + rnorm(1, 0, 0.1) - opt)^2) pheno.temp = 0 for(j in 451:500){ if(Genome[2*i,j] == Genome[2*i - 1,j]){pheno.temp = pheno.temp + 2*Genome[2*i,j]} else{pheno.temp = pheno.temp + dom.effect[j]} } dist_opt = dist_opt + ((pheno.temp + rnorm(1, 0, 0.1) - opt)^2) geno.off[i] = dist_opt } return(geno.off) } dist.genotype.offspring <-function(Genome,add.effect, dom.effect,opt, opt.foc, Npop) { geno.off = c(NULL) for(i in 1:Npop){ dist_opt = 0 pheno.temp = 0 for(j in 1:50){ if(Genome[2*i,j] == Genome[2*i - 1,j]){pheno.temp = pheno.temp + 2*Genome[2*i,j]} else{pheno.temp = pheno.temp + dom.effect[j]} } dist_opt = dist_opt + ((pheno.temp + rnorm(1, 0, 0) - opt.foc)^2) pheno.temp = 0 for(j in 51:100){ if(Genome[2*i,j] == Genome[2*i - 1,j]){pheno.temp = pheno.temp + 2*Genome[2*i,j]} else{pheno.temp = pheno.temp + dom.effect[j]} } dist_opt = dist_opt + ((pheno.temp + rnorm(1, 0, 0) - opt)^2) pheno.temp = 0 for(j in 101:150){ if(Genome[2*i,j] == Genome[2*i - 1,j]){pheno.temp = pheno.temp + 2*Genome[2*i,j]} else{pheno.temp = pheno.temp + dom.effect[j]} } dist_opt = dist_opt + ((pheno.temp + rnorm(1, 0, 0) - opt)^2) pheno.temp = 0 for(j in 151:200){ if(Genome[2*i,j] == Genome[2*i - 1,j]){pheno.temp = pheno.temp + 2*Genome[2*i,j]} else{pheno.temp = pheno.temp + dom.effect[j]} } dist_opt = dist_opt + ((pheno.temp + rnorm(1, 0, 0) - opt)^2) pheno.temp = 0 for(j in 201:250){ if(Genome[2*i,j] == Genome[2*i - 1,j]){pheno.temp = pheno.temp + 2*Genome[2*i,j]} else{pheno.temp = pheno.temp + dom.effect[j]} } dist_opt = dist_opt + ((pheno.temp + rnorm(1, 0, 0) - opt)^2) pheno.temp = 0 for(j in 251:300){ if(Genome[2*i,j] == Genome[2*i - 1,j]){pheno.temp = pheno.temp + 2*Genome[2*i,j]} else{pheno.temp = pheno.temp + dom.effect[j]} } dist_opt = dist_opt + ((pheno.temp + rnorm(1, 0, 0) - opt)^2) pheno.temp = 0 for(j in 301:350){ if(Genome[2*i,j] == Genome[2*i - 1,j]){pheno.temp = pheno.temp + 2*Genome[2*i,j]} else{pheno.temp = pheno.temp + dom.effect[j]} } dist_opt = dist_opt + ((pheno.temp + rnorm(1, 0, 0) - opt)^2) pheno.temp = 0 for(j in 351:400){ if(Genome[2*i,j] == Genome[2*i - 1,j]){pheno.temp = pheno.temp + 2*Genome[2*i,j]} else{pheno.temp = pheno.temp + dom.effect[j]} } dist_opt = dist_opt + ((pheno.temp + rnorm(1, 0, 0) - opt)^2) pheno.temp = 0 for(j in 401:450){ if(Genome[2*i,j] == Genome[2*i - 1,j]){pheno.temp = pheno.temp + 2*Genome[2*i,j]} else{pheno.temp = pheno.temp + dom.effect[j]} } dist_opt = dist_opt + ((pheno.temp + rnorm(1, 0, 0) - opt)^2) pheno.temp = 0 for(j in 451:500){ if(Genome[2*i,j] == Genome[2*i - 1,j]){pheno.temp = pheno.temp + 2*Genome[2*i,j]} else{pheno.temp = pheno.temp + dom.effect[j]} } dist_opt = dist_opt + ((pheno.temp + rnorm(1, 0, 0) - opt)^2) geno.off[i] = dist_opt } return(geno.off) } Phenotype.offspring <-function(Genotype, Npop) { pheno.off = c(NULL) for(i in 1:Npop){ pheno.off[i] = Genotype[i] + rnorm(1, 0, 0.1) # genotype + environmental effects } return(pheno.off) } genotype.offspring <-function(Genome,add.effect, dom.effect, Npop) { geno.off = c(NULL) for(i in 1:Npop){ pheno.temp = 0 for(j in 1:50){ if(Genome[2*i,j] == Genome[2*i - 1,j]){pheno.temp = pheno.temp + 2*Genome[2*i,j]} else{pheno.temp = pheno.temp + dom.effect[j]} } geno.off[i] = pheno.temp } return(geno.off) } is.wholenumber <- function(x, tol = .Machine$double.eps^0.5) abs(x - round(x)) < tol variances <- function(phenotype, genotype, genome,add.effect, Npop,dom.effect){ mean.phenotype = mean(phenotype) phenotypic.variance = (sum((phenotype - mean.phenotype)^2))/Npop mean.genotype = mean(genotype) genetic.variance = (sum((genotype - mean.genotype)^2))/Npop par_locus<-0 for(i in 1:50){ phenotype_loc = c(NULL) for(j in 1:Npop){ if(genome[2*j, i] == genome[2*j-1, i]){phenotype_loc[j] = 2*genome[2*j, i]} else{phenotype_loc[j] = dom.effect[i]} } mean.phenotype = mean(phenotype_loc) par_locus = par_locus + (sum((phenotype_loc - mean.phenotype)^2))/Npop } genic.variance = 0 cov.additive = 0 cov.dominance = 0 cov.add.dom = 0 cov.dom.add = 0 cov.genetic = 0 dominance.variance = 0 autozygite.dominance.variance = 0 cov.breed.autozygote.dev = 0 var.inbreeding.depression = 0 hetero.cor = 0 ### Within-loci components for(i in 1:50){ if(var(genome[,i]) > 0){ Homo.1 = 0 Homo.2 = 2*add.effect[i] Hetero = dom.effect[i] ### Frequencies of alleles nb.all1 = 0 for(j in 1:(2*Npop)){ if(genome[j,i] == 0){nb.all1 = nb.all1 + 1} } freq.all1 = nb.all1/(2*Npop) freq.all2 = 1 - freq.all1 ### Genotypique values moyenne.loc = ((freq.all1^2)*Homo.1) + (2*freq.all1*freq.all2*Hetero) + ((freq.all2^2)*Homo.2) Geno.val.Homo1 = Homo.1 - moyenne.loc Geno.val.Homo2 = Homo.2 - moyenne.loc Geno.val.Hetero = Hetero - moyenne.loc ### Breeding values alpha.all1 = freq.all1*Geno.val.Homo1 + freq.all2*Geno.val.Hetero alpha.all2 = freq.all1*Geno.val.Hetero + freq.all2*Geno.val.Homo2 Breeding.Homo1 = 2*alpha.all1 Breeding.Homo2 = 2*alpha.all2 Breeding.Hetero = alpha.all1 + alpha.all2 ### dominance deviations delta.11 = Geno.val.Homo1 - Breeding.Homo1 delta.22 = Geno.val.Homo2 - Breeding.Homo2 delta.12 = Geno.val.Hetero - Breeding.Hetero genic.variance = genic.variance + (2*freq.all1*(alpha.all1^2)) + (2*freq.all2*(alpha.all2^2)) dominance.variance = dominance.variance + (2*freq.all1*freq.all2*(Hetero-((Homo.1+Homo.2)/2)))^2 autozygite.dominance.variance = autozygite.dominance.variance + (freq.all1*(delta.11^2)) + (freq.all2*(delta.22^2)) - (((freq.all1*delta.11) + (freq.all2*delta.22))^2) cov.breed.autozygote.dev = cov.breed.autozygote.dev + (freq.all1*alpha.all1*delta.11) + (freq.all2*alpha.all2*delta.22) var.inbreeding.depression = var.inbreeding.depression + ((freq.all1*delta.11) + (freq.all2*delta.22))^2 hetero.cor = hetero.cor + (freq.all1*delta.11) + (freq.all2*delta.22) } } ### additive covariances for(i in 1:50){ Homo.1 = 0 Homo.2 = 2*add.effect[i] Hetero = dom.effect[i] ### Frequencies of alleles nb.all1 = 0 for(h in 1:(2*Npop)){ if(genome[h,i] == 0){nb.all1 = nb.all1 + 1} } freq.all1 = nb.all1/(2*Npop) freq.all2 = 1 - freq.all1 ### Genotypique values moyenne.loc = ((freq.all1^2)*Homo.1) + (2*freq.all1*freq.all2*Hetero) + ((freq.all2^2)*Homo.2) Geno.val.Homo1 = Homo.1 - moyenne.loc Geno.val.Homo2 = Homo.2 - moyenne.loc Geno.val.Hetero = Hetero - moyenne.loc ### Breeding values alpha.all1 = freq.all1*Geno.val.Homo1 + freq.all2*Geno.val.Hetero alpha.all2 = freq.all1*Geno.val.Hetero + freq.all2*Geno.val.Homo2 Breeding.Homo1 = 2*alpha.all1 Breeding.Homo2 = 2*alpha.all2 Breeding.Hetero = alpha.all1 + alpha.all2 ### dominance deviations delta.11 = Geno.val.Homo1 - Breeding.Homo1 delta.22 = Geno.val.Homo2 - Breeding.Homo2 delta.12 = Geno.val.Hetero - Breeding.Hetero product.add1 = c(NULL) product.dom1 = c(NULL) product.geno1 = c(NULL) for(h in 1:Npop){ if(genome[2*h,i] == 0 && genome[2*h-1,i] == 0){product.add1[h] = Breeding.Homo1 product.dom1[h] = delta.11 product.geno1[h] = Homo.1} else if(genome[2*h,i] == add.effect[i] && genome[2*h-1,i] == add.effect[i]){product.add1[h] = Breeding.Homo2 product.dom1[h] = delta.22 product.geno1[h] = Homo.2} else{product.add1[h] = Breeding.Hetero product.dom1[h] = delta.12 product.geno1[h] = Hetero} } for(j in (i+1):50){ if(i != 50){ Homo.1 = 0 Homo.2 = 2*add.effect[j] Hetero = dom.effect[j] ### Frequencies of alleles nb.all1 = 0 for(h in 1:(2*Npop)){ if(genome[h,j] == 0){nb.all1 = nb.all1 + 1} } freq.all1 = nb.all1/(2*Npop) freq.all2 = 1 - freq.all1 ### Genotypique values moyenne.loc = ((freq.all1^2)*Homo.1) + (2*freq.all1*freq.all2*Hetero) + ((freq.all2^2)*Homo.2) Geno.val.Homo1 = Homo.1 - moyenne.loc Geno.val.Homo2 = Homo.2 - moyenne.loc Geno.val.Hetero = Hetero - moyenne.loc ### Breeding values alpha.all1 = freq.all1*Geno.val.Homo1 + freq.all2*Geno.val.Hetero alpha.all2 = freq.all1*Geno.val.Hetero + freq.all2*Geno.val.Homo2 Breeding.Homo1 = 2*alpha.all1 Breeding.Homo2 = 2*alpha.all2 Breeding.Hetero = alpha.all1 + alpha.all2 ### dominance deviations delta.11 = Geno.val.Homo1 - Breeding.Homo1 delta.22 = Geno.val.Homo2 - Breeding.Homo2 delta.12 = Geno.val.Hetero - Breeding.Hetero product.add2 = c(NULL) product.dom2 = c(NULL) product.geno2 = c(NULL) for(h in 1:Npop){ if(genome[2*h,j] == 0 && genome[2*h-1,j] == 0){product.add2[h] = Breeding.Homo1 product.dom2[h] = delta.11 product.geno2[h] = Homo.1} else if(genome[2*h,j] == add.effect[j] && genome[2*h-1,j] == add.effect[j]){product.add2[h] = Breeding.Homo2 product.dom2[h] = delta.22 product.geno2[h] = Homo.2} else{product.add2[h] = Breeding.Hetero product.dom2[h] = delta.12 product.geno2[h] = Hetero} } cov.additive = cov.additive + sum((product.add1-mean(product.add1))*(product.add2-mean(product.add2)))*(1/Npop) cov.dominance = cov.dominance + sum((product.dom1-mean(product.dom1))*(product.dom2-mean(product.dom2)))*(1/Npop) cov.add.dom = cov.add.dom + sum((product.add1-mean(product.add1))*(product.dom2-mean(product.dom2)))*(1/Npop) cov.dom.add = cov.dom.add + sum((product.dom1-mean(product.dom1))*(product.add2-mean(product.add2)))*(1/Npop) cov.genetic = cov.genetic + sum((product.geno1-mean(product.geno1))*(product.geno2-mean(product.geno2)))*(1/Npop) } } } cov.additive = 2*cov.additive cov.dominance = 2*cov.dominance cov.add.dom = 2*cov.add.dom cov.dom.add = 2*cov.dom.add cov.genetic = 2*cov.genetic hetero.cor = (hetero.cor^2) - var.inbreeding.depression summary.variance = c(phenotypic.variance, genetic.variance,par_locus, genic.variance,dominance.variance,autozygite.dominance.variance,cov.breed.autozygote.dev,var.inbreeding.depression,hetero.cor,cov.additive,cov.dominance,cov.add.dom,cov.dom.add,cov.genetic) return(summary.variance) } variances.fitness <- function(dist.phenotype, dist.genotype,genome,add.effect, Npop, L,dom.effect,opt, opt.foc,om.2){ phenotype1 = exp(-(dist.phenotype)/(2*(om.2))) mean.phenotype = mean(phenotype1) phenotypic.variance = (sum((phenotype1 - mean.phenotype)^2))/Npop genotype1 = exp(-(dist.genotype)/(2*(om.2))) mean.genotype = mean(genotype1) genetic.variance = (sum((genotype1 - mean.genotype)^2))/Npop par_locus<-0 for(i in 1:50){ phenotype_loc = c(NULL) for(j in 1:Npop){ if(genome[2*j, i] == genome[2*j-1, i]){phenotype_loc[j] = exp(-((2*genome[2*j, i]-opt.foc)^2)/(2*om.2))} else{phenotype_loc[j] = exp(-((dom.effect[i]-opt.foc)^2)/(2*om.2))} } mean.phenotype2 = mean(phenotype_loc) par_locus = par_locus + (sum((phenotype_loc - mean.phenotype2)^2))/Npop } for(i in 51:500){ phenotype_loc = c(NULL) for(j in 1:Npop){ if(genome[2*j, i] == genome[2*j-1, i]){phenotype_loc[j] = exp(-((2*genome[2*j, i]-opt)^2)/(2*om.2))} else{phenotype_loc[j] = exp(-((dom.effect[i]-opt)^2)/(2*om.2))} } mean.phenotype2 = mean(phenotype_loc) par_locus = par_locus + (sum((phenotype_loc - mean.phenotype2)^2))/Npop } genic.variance = 0 dominance.variance = 0 autozygite.dominance.variance = 0 cov.breed.autozygote.dev = 0 var.inbreeding.depression = 0 hetero.cor = 0 ### Within-loci components for(i in 1:50){ if(var(genome[,i]) > 0){ Homo.1 = exp(-((0-opt.foc)^2)/(2*(om.2))) Homo.2 = exp(-((2*add.effect[i]-opt.foc)^2)/(2*(om.2))) Hetero = exp(-((dom.effect[i]-opt.foc)^2)/(2*(om.2))) ### Frequencies of alleles nb.all1 = 0 for(j in 1:(2*Npop)){ if(genome[j,i] == 0){nb.all1 = nb.all1 + 1} } freq.all1 = nb.all1/(2*Npop) freq.all2 = 1 - freq.all1 ### Genotypique values moyenne.loc = ((freq.all1^2)*Homo.1) + (2*freq.all1*freq.all2*Hetero) + ((freq.all2^2)*Homo.2) Geno.val.Homo1 = Homo.1 - moyenne.loc Geno.val.Homo2 = Homo.2 - moyenne.loc Geno.val.Hetero = Hetero - moyenne.loc ### Breeding values alpha.all1 = freq.all1*Geno.val.Homo1 + freq.all2*Geno.val.Hetero alpha.all2 = freq.all1*Geno.val.Hetero + freq.all2*Geno.val.Homo2 Breeding.Homo1 = 2*alpha.all1 Breeding.Homo2 = 2*alpha.all2 Breeding.Hetero = alpha.all1 + alpha.all2 ### dominance deviations delta.11 = Geno.val.Homo1 - Breeding.Homo1 delta.22 = Geno.val.Homo2 - Breeding.Homo2 delta.12 = Geno.val.Hetero - Breeding.Hetero genic.variance = genic.variance + (2*freq.all1*(alpha.all1^2)) + (2*freq.all2*(alpha.all2^2)) dominance.variance = dominance.variance + (2*freq.all1*freq.all2*(Hetero-((Homo.1+Homo.2)/2)))^2 autozygite.dominance.variance = autozygite.dominance.variance + (freq.all1*(delta.11^2)) + (freq.all2*(delta.22^2)) - (((freq.all1*delta.11) + (freq.all2*delta.22))^2) cov.breed.autozygote.dev = cov.breed.autozygote.dev + (freq.all1*alpha.all1*delta.11) + (freq.all2*alpha.all2*delta.22) var.inbreeding.depression = var.inbreeding.depression + ((freq.all1*delta.11) + (freq.all2*delta.22))^2 hetero.cor = hetero.cor + (freq.all1*delta.11) + (freq.all2*delta.22) } } for(i in 51:500){ if(var(genome[,i]) > 0){ Homo.1 = exp(-((0-opt)^2)/(2*(om.2))) Homo.2 = exp(-((2*add.effect[i]-opt)^2)/(2*(om.2))) Hetero = exp(-((dom.effect[i]-opt)^2)/(2*(om.2))) ### Frequencies of alleles nb.all1 = 0 for(j in 1:(2*Npop)){ if(genome[j,i] == 0){nb.all1 = nb.all1 + 1} } freq.all1 = nb.all1/(2*Npop) freq.all2 = 1 - freq.all1 ### Genotypique values moyenne.loc = ((freq.all1^2)*Homo.1) + (2*freq.all1*freq.all2*Hetero) + ((freq.all2^2)*Homo.2) Geno.val.Homo1 = Homo.1 - moyenne.loc Geno.val.Homo2 = Homo.2 - moyenne.loc Geno.val.Hetero = Hetero - moyenne.loc ### Breeding values alpha.all1 = freq.all1*Geno.val.Homo1 + freq.all2*Geno.val.Hetero alpha.all2 = freq.all1*Geno.val.Hetero + freq.all2*Geno.val.Homo2 Breeding.Homo1 = 2*alpha.all1 Breeding.Homo2 = 2*alpha.all2 Breeding.Hetero = alpha.all1 + alpha.all2 ### dominance deviations delta.11 = Geno.val.Homo1 - Breeding.Homo1 delta.22 = Geno.val.Homo2 - Breeding.Homo2 delta.12 = Geno.val.Hetero - Breeding.Hetero genic.variance = genic.variance + (2*freq.all1*(alpha.all1^2)) + (2*freq.all2*(alpha.all2^2)) dominance.variance = dominance.variance + (2*freq.all1*freq.all2*(Hetero-((Homo.1+Homo.2)/2)))^2 autozygite.dominance.variance = autozygite.dominance.variance + (freq.all1*(delta.11^2)) + (freq.all2*(delta.22^2)) - (((freq.all1*delta.11) + (freq.all2*delta.22))^2) cov.breed.autozygote.dev = cov.breed.autozygote.dev + (freq.all1*alpha.all1*delta.11) + (freq.all2*alpha.all2*delta.22) var.inbreeding.depression = var.inbreeding.depression + ((freq.all1*delta.11) + (freq.all2*delta.22))^2 hetero.cor = hetero.cor + (freq.all1*delta.11) + (freq.all2*delta.22) } } hetero.cor = (hetero.cor^2) - var.inbreeding.depression summary.variance = c(phenotypic.variance, genetic.variance,par_locus, genic.variance,dominance.variance,autozygite.dominance.variance,cov.breed.autozygote.dev,var.inbreeding.depression,hetero.cor) return(summary.variance) } frequence <- function(genome,L,Npop){ freq.all0=c(NULL) for(i in 1:L){ nb.all1 = 0 for(j in 1:(2*Npop)){ if(genome[j,i] == 0){nb.all1 = nb.all1 + 1} } freq.all1 = nb.all1/(2*Npop) freq.all0[i]=freq.all1 } return(freq.all0) } L = 500 Npop = 100 var.add.eff = 0.05 dom.val = 1 var.dom.eff = 0.05 U = 0.1 om_2 = 1 selfing.rate = 1 Simulation.programme <- function(Number.sim, L, var.add.eff, dom.val,var.dom.eff, Npop, U, om_2, selfing.rate){ Nsimul = 0 repeat{ Nsimul = Nsimul + 1 mean.fitness = c(NULL) opt = 0 opt.foc = 0 om.2 = om_2 ### Initialisation add.effect = c(rnorm(n = L, mean = 0, sd = var.add.eff^0.5)) dom.effect = c(NULL) for(i in 1:L){ if(add.effect[i] >= 0){dom.effect[i] = rnorm(1, dom.val*2*add.effect[i], var.dom.eff)} else{dom.effect[i] = dom.effect[i] = rnorm(1, (1-dom.val)*2*add.effect[i], var.dom.eff)} } geno.init = c(rep(0, times=L)) for(i in 1:(2*Npop)){ if(i == 1)(genome = geno.init) else{genome = rbind(genome, geno.init)} } dist.phenotype = dist.phenotype.offspring(genome, add.effect,dom.effect,opt, opt.foc, Npop) generation = 0 repeat{ generation = generation + 1 fitness = exp(-(dist.phenotype)/(2*om_2)) mean.fitness[generation] = mean(fitness) Parent.sel = Selection(Npop, dist.phenotype, om_2, selfing.rate) genome = Genome.offspring(Parent.sel, L, Npop, U, genome, add.effect) genotype = genotype.offspring(genome, add.effect,dom.effect, Npop) phenotype = Phenotype.offspring(genotype, Npop) dist.phenotype = dist.phenotype.offspring(genome, add.effect,dom.effect,opt, opt.foc, Npop) dist.genotype = dist.genotype.offspring(genome, add.effect,dom.effect,opt,opt.foc, Npop) if( is.wholenumber(generation/1000) == TRUE && (generation/1000) > 1 ){ mean.1<-mean(mean.fitness[(generation-999):generation]) mean.2<-mean(mean.fitness[(generation-1999):(generation-1000)]) if( abs(1 - (mean.1/mean.2)) <= 0.01){break} } } fitness.eq = mean.fitness[generation] variances.genetic = variances(phenotype,genotype, genome, add.effect, Npop,dom.effect) var.fit = variances.fitness(dist.phenotype, dist.genotype, genome, add.effect, Npop, L,dom.effect,opt,opt.foc,om.2) tout.fit = c(var.fit,fitness.eq) #freq_cumule = frequence(genome, L, Npop) ## Step 2 opt.foc = 2.5 mut=0 fis = selfing.rate/(2-selfing.rate) dyna_obs = c(NULL) dyna_kelly = c(NULL) dyna_bias = c(NULL) for(z in 1:21){ if(z == 1){ generation = generation + 1 dist.phenotype = dist.phenotype.offspring(genome, add.effect,dom.effect,opt, opt.foc, Npop) dist.genotype = dist.genotype.offspring(genome, add.effect,dom.effect,opt,opt.foc, Npop) fitness = exp(-((dist.phenotype)/(2*om_2))) mean.fitness[generation] = mean(fitness) vari.time = variances.fitness(dist.phenotype,dist.genotype, genome, add.effect, Npop, L,dom.effect,opt,opt.foc,om.2) dyna_kelly[z] = ((1+fis)*vari.time[4] + ((3*fis)+(fis^2))*vari.time[7]+(fis^2)*vari.time[6] + (vari.time[2]-vari.time[3]))/mean.fitness[generation] dyna_bias[z] = ((1+fis)*vari.time[4] + (vari.time[2]-vari.time[3]))/mean.fitness[generation] Parent.sel = Selection(Npop, dist.phenotype, om_2, selfing.rate) genome = Genome.offspring(Parent.sel, L, Npop, U, genome, add.effect) genotype = genotype.offspring(genome, add.effect,dom.effect, Npop) phenotype = Phenotype.offspring(genotype, Npop) dist.phenotype = dist.phenotype.offspring(genome, add.effect,dom.effect,opt, opt.foc, Npop) dist.genotype = dist.genotype.offspring(genome, add.effect,dom.effect,opt,opt.foc, Npop) } else if(z == 21){ generation = generation + 1 fitness = exp(-(dist.phenotype)/(2*om_2)) mean.fitness[generation] = mean(fitness) dyna_obs[z-1] = mean.fitness[generation] - mean.fitness[generation-1] } else{ generation = generation + 1 fitness = exp(-(dist.phenotype)/(2*om_2)) mean.fitness[generation] = mean(fitness) vari.time = variances.fitness(dist.phenotype,dist.genotype, genome, add.effect, Npop, L,dom.effect,opt, opt.foc,om.2) dyna_kelly[z] = ((1+fis)*vari.time[4] + ((3*fis)+(fis^2))*vari.time[7]+(fis^2)*vari.time[6] + (vari.time[2]-vari.time[3]))/mean.fitness[generation] dyna_bias[z] = ((1+fis)*vari.time[4] + (vari.time[2]-vari.time[3]))/mean.fitness[generation] dyna_obs[z-1] = mean.fitness[generation] - mean.fitness[generation-1] Parent.sel = Selection(Npop, dist.phenotype, om_2, selfing.rate) genome = Genome.offspring(Parent.sel, L, Npop, U, genome, add.effect) genotype = genotype.offspring(genome, add.effect,dom.effect, Npop) phenotype = Phenotype.offspring(genotype, Npop) dist.phenotype = dist.phenotype.offspring(genome, add.effect,dom.effect,opt, opt.foc, Npop) dist.genotype = dist.genotype.offspring(genome, add.effect,dom.effect,opt,opt.foc, Npop) } } if(Nsimul == 1){recap = variances.genetic recap.fit = tout.fit recap.dyfit = dyna_obs recap.kelly = dyna_kelly recap.bias = dyna_bias #freq_recap = freq_cumule } else{recap = rbind(recap, variances.genetic) recap.fit = rbind(recap.fit, tout.fit) recap.dyfit = rbind(recap.dyfit, dyna_obs) recap.kelly = rbind(recap.kelly, dyna_kelly) recap.bias = rbind(recap.bias, dyna_bias) #freq_recap = c(freq_recap,freq_cumule) } if(Nsimul == Number.sim) break } write.table(recap, file=paste("variance_s",selfing.rate,"_d",dom.val,"_om",om_2,"_N",Npop,".txt",sep=""),row.names = F, dec = ".") write.table(recap.fit, file=paste("variance_fit_s",selfing.rate,"_d",dom.val,"_om",om_2,"_N",Npop,".txt",sep=""),row.names = F, dec = ".") write.table(recap.dyfit, file=paste("dyfit_s",selfing.rate,"_d",dom.val,"_om",om_2,"_N",Npop,".txt",sep=""),row.names = F, dec = ".") write.table(recap.kelly, file=paste("kelly_s",selfing.rate,"_d",dom.val,"_om",om_2,"_N",Npop,".txt",sep=""),row.names = F, dec = ".") write.table(recap.bias, file=paste("bias_s",selfing.rate,"_d",dom.val,"_om",om_2,"_N",Npop,".txt",sep=""),row.names = F, dec = ".") #write.table(freq_recap, file=paste("freq_s",selfing.rate,"_d",dom.val,"_om",om_2,"_N",Npop,".txt",sep=""),row.names = F, dec = ".") } Simulation.programme(10,500,0.05,1, 0.05,1000,0.05, 1, 0)