#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)
s <- 0.01 # selection coefficient of the mutant type
init.num <- 1 # initial number of the mutant type
two.N <-10 # population size 2N
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)+(2N-j)] = j(1+s)/(js + 2N)
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+two.N) #compute the post-selection expected frequency
j[i,k]=rbinom(1,two.N,p.star) #generate single binomial random variable with parameters two.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/two.N,type="l",ylab="freq(A)",xlab="generation")