#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")