function evolve5()

evolve5= function( n.pop1=5000, n.pop2=100, a.pop=0.4, n.gen=400, n.samp=100, f.init=0.2, s.adv=0.01, dom = TRUE, histogram=TRUE, plot.mean=FALSE, color=TRUE, title = "") 

{

#Based on code by Lex Flagel in his blog A Hopeful Monster

#

#      -allows for a linear time-varying population

#      -uses simplified forms of the Hardy-Weinberg expressions for calculating the change of allele frequency  (see Falconer and Mackay 4th ed., p.29, equations 2.7 and 2.8)

#      -case of either no dominance or complete dominance

#      -assumes: s.adv is small, or allele frequency is small.  Usually I use 0 < s.adv < 0.2, as in F&M page 29.

#   evolve5() employs the Heaviside function to create a step in population at a specified time. The population size starts out at n.pop1, and then suddenly changes to n.pop2 at the time in the n.gen generations defined by a.pop

#

#  Usage: 

#   Example of bottleneck (sudden drop in population size):

#   evolve5( n.pop1=5000, n.pop2=100, a.pop=0.4, n.gen=400, n.samp=100, f.init=0.2, s.adv=0.015, dom = TRUE, histogram=TRUE, plot.mean=FALSE, color=TRUE, title = "Example of bottleneck")

#  Counter-example - no bottleneck (population remains large throughout):

#   evolve5( n.pop1=5000, n.pop2=5000, a.pop=0.4, n.gen=400, n.samp=100, f.init=0.2, s.adv=0.015, dom = TRUE, histogram=TRUE, plot.mean=FALSE, color=TRUE, "Population remains large") 

# n.pop1 = starting population size

# n.pop2 = ending population size

# a.pop = normalized time (value between 0 and 1) of population change from n.pop1 to n.pop2 

# n.gen = number of generations to run

# n.samp = number of individual lines (monte carlo realizations) to calculate

# s.adv = reproductive advantage (=0.0 for just random drift)

# f.init = initial allele frequency in population

#

    k = matrix(-999, nr=n.samp, nc=n.gen-1) #matrix to store results

    sz.vector = n.pop1 + (n.pop2-n.pop1)*Heaviside(1:n.gen, as.integer(a.pop*n.gen))

    graphics.off()

#

for (i in 1:n.samp)

 {

#  Falcon & Mackay parameters

      q = f.init  # allele frequency

      p = 1-q

      s = s.adv

#

   for (j in 1:(n.gen-1) )   # n.gen generations require n.gen - 1 samples

{     

      dq = 0.5 * s * q * p                  #  No dominance

      if (dom) {dq = s * q * q * p}      #  Complete dominance

      prob= min((q + dq), 1.0)    #truncate prob from above

      prob= max(prob, 0.0)         #truncate prob from below

      sz = sz.vector[j]                    # get population size from pre-calculated size vector

      q = rbinom(1, size=sz, prob) / sz

      p = 1 - q

#      

      k[i, j] = q  #update matrix

}

}

#

k = cbind(rep(f.init, n.samp), k) #bind first gen. first column of matrix

#

#

#######

##       that’s it for sampling, code below is plotting

#

# Plot all n.samp random walks

dev.new(width=5,height=4)

#jpeg(“randomwalks.jpg”)

plot(1,1, xlim=c(1,n.gen), ylim=c(0,1), ty='n', xlab='generation', ylab='allele freq')

abline(h=0.5) #

mtext(line=0, cex=0.7, paste('pop=', n.pop1, n.pop2, 'gen=', n.gen, 'samp=', n.samp, 'f.init=', f.init, 'adv=', s.adv, 'dom=', dom ))

mtext(line=1, cex=1.0, title)

#

for (i in 1:n.samp){

ic=1

if (color) {ic = round(502 * i/n.samp, digits=0)}

  lines(1:n.gen, k[i, ], lwd=.3, col=ic)  #draw n.samp random walks as lines in color col=ic

  points(n.gen+0.1, k[i, n.gen-1], pch='-'    )  #add red rug on right

}

#

if (histogram) {

#  Histogram of final allele frequencies

dev.new(width=5,height=4)

hist (k[, n.gen], xlim=c(0,1),  main='', xlab="allele freq") 

mtext(line=0, cex=0.7, paste('Gen', n.gen, 'mean=', round(mean(k[, n.gen]),3), 'sigma=', round(sd(k[, n.gen]),3) ) )

mtext(line=1, cex=1, "Final allele frequency distribution" )

}

if (plot.mean)  {

# Plot the mean and std. dev. as function of generation

dev.new(width=5,height=4)

q.mean =  rep(-999, n.gen)

q.std =   rep(-999, n.gen)

q.sqmean = rep(-999, n.gen)

for (j in 1:(n.gen))

   {    # accumulate mean and sd of freq at each generation

          q.mean[j] = mean(k[,j])

          q.sqmean[j] = mean( (k[,j])*(k[,j]) )

          q.std[j] = sd(k[,j])

   }

plot(1,1, xlim=c(1,n.gen-1), ylim=c(0, 1 ), ty='n'  , xlab='Generation', ylab='Means of q and q^2')

lines (1:n.gen,  q.mean,  lwd=2)

lines (1:n.gen,  q.mean-q.std, lwd=0.5, lty=3)

lines (1:n.gen,  q.mean+q.std, lwd=0.5, lty=3)

lines (1:n.gen,  q.sqmean, lwd=2)

}

} ###### END 

###############################

#  Turn off all graphics windows

#  graphics.off()

#