# This simualates the Wright-Fisher Model of selection and random genetic drift
# for an asexual organism with two alleles and no mutation
# By R. Gomulkiewicz (5 September 2013)
# Updated 8 Sep 2015
s <- 0.01 # selection coefficient of the mutant type
init.num <- 1 # initial number of the mutant type
N <-10 # population size N
reps <- 15 # number of replicate simulation runs
gens <- 40 # generations of evolution to simulate for each replicate run
#create a matrix with gens+1 rows and reps columns whose every entry is the initial number of mutants
#below,we will replace entries 2 through gens+1 of each column with simulated numbers of mutants in each replicate run
j=matrix(init.num,gens+1,reps)
#loops that execute the replicate simulations
#uses simplified form of the post-selection expected frequency pstar = j(1+s)/[j(1+s)+(N-j)] = j(1+s)/(js + N)
for(k in 1:reps){ #this loops over the number of replicate runs
for(i in 2:gens+1) { #simulate W-F model starting with init.num mutants
p.star <- j[i-1,k]*(1+s)/(j[i-1,k]*s+N) #compute the post-selection expected frequency, given j
j[i,k]=rbinom(1,N,p.star) #generate the next j as a single binomial random variable with parameters N and p.star
}
}
#Show the results as mutant frequencies (i.e., plot each mutant count divided by 2N)
#Display all the replicate runs on a single plot
matplot(0:gens,j/N,type="l",ylab="freq(A)",xlab="generation")