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