#Project frequencies of a new mutation with selection coefficient s #Show DETERMINISTIC projections for 3 values of s: s < 0, s = 0, and s > 0 #By R. Gomulkiewicz (29 Aug 2013) #Updated 8 Sep 2015 #create a single list called "s" that holds all three values of s s <- c(0.01, 0, -0.01) cases <- length(s) #length(x) returns the number of list items #size of population when mutation first appears initial.size <- 20 #starting population size steps <- 200 #number of generations (iterations) #create a table to hold the mutant frequencies at all the steps for each value of s q <- matrix(0, steps+1, cases) #store the initial mutant frequency for each case q[1,] <- 1/initial.size #this function computes the mutant frequency at the next time step next.q <- function(q, s) {q*(1+s)/(q*(1+s)+(1-q))} #loop through the cases of s #for each case, run a loop that projects the mutant at all steps for(j in 1:cases){ for(i in 1:steps){ q[i+1,j] <- next.q(q[i,j],s[j]) } } #plot the results for all cases together matplot(0:steps,q,type="l",xlab="steps",ylab="mutant frequency") #log plots of the ratio q/(1-q) for all three cases matplot(0:steps,q/(1-q),type="l",xlab="steps",ylab="mutant:non-mutant ratio",log="y") #display frequency dynamics and ratio dynamics together par(mfrow=c(2,1)) #par(mfrow=c(a,b) displays a*b figures as a table with a rows and b columns matplot(0:steps,q,type="l",xlab="steps",ylab="mutant frequency") matplot(0:steps,q/(1-q),type="l",xlab="time step",ylab="mutant:non-mutant ratio",log="y")