Preditor-prey modeling

This could be fun. I wonder if the Monte Carlo code could be modified to make a preditor-prey model. Two populations, each subject to survival advantage forces. But here the survival advantage for one population/species is not fixed, it is a function of it's own population, the population of the other (preditor) species, it's own available food supply, and other environmental forces.

The basic equations are the Lotka–Volterra equations.

Random variable for each population, reproducing as in function evolve(4), with same assumptions about random mating, etc.  But also with preditor rate term in the binomial step.

Of course this is a very well-studied problem, a standard student exercise. I've come across it in a variety of places, including numerous lovely interactive graphic simulations on the web.   There are even Monte Carlo solutions online, so I'm certainly talking about a well-known problem. Still, as an exercise, it will be fun and instructive to do it myself!!

function mutate1()

mutate1= function( n.gen=50000, n.variables=2, A.init=0.00, mu.init=0.0001, nu.init = 0.00001, title = "Mutation: equilibrating allele frequencies")
{
    # n.samp = 1 means there is no sample population, we're computing the mean.
 
    # mu = mutation rate per generation A to a.
    # nu = mutation rate per generation a to A.
 
    # p = allele frequency of A(for all individuals and for all generations)
 
    # Parameters for environment
    mu = mu.init     # mutation rate per generation A to a.
    nu = nu.init     # mutation rate per generation a to A.
 
    # Initialize matrix
    k = matrix(-999, nr=n.variables, nc=n.gen) #matrix to store values of p
 
    for (i in 1:n.variables)
    {
        k[1,1] = A.init    #  initial freq of A
        k[2,1] = 1- A.init    #  initial freq of a
     
        for (j in 2:(n.gen) )   # n.gen generations
        {
            k[1, j] = (1-mu) *  k[1, j-1] + nu * (1 - k[1, j-1]) # freq of A
            k[2, j] = (1-mu) *  k[2, j-1] + nu * (1 - k[2, j-1]) # freq of a
        }
    }
    #######
    #
    # Plot allele frequency vs. generation
    graphics.off()
    plot(1,1, xlim=c(1,n.gen), ylim=c(0,1), ty='n', xlab='Generation', ylab='Allele Frequency' )
    mtext(line=0, cex=0.7, paste('gen=', n.gen,  'A.init=', A.init,  'a.init=', 1-A.init, 'mu.init=', mu.init, 'nu.init=', nu.init ))
    mtext(line=1, cex=1.0, title)
    #
    lines(1:n.gen, k[1, ], lwd=2, col=1, lty=1)  # allele A, black
    lines(1:n.gen, k[2, ], lwd=2, col=4, lty=2)  # allele a, blue
    legend(45000, 0.9, c("A","a"), cex=0.8,
    col=c(1,4), lty=c(1,1), lwd=2)
 
} ###### END

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()

#

function evolve4()

evolve4= function( n.pop1=1000, n.pop2=1000, n.gen=60, n.samp=100, f.init=0.05, s.adv=0.2, dom = FALSE, histogram=FALSE, plot.mean=TRUE, color=FALSE, 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.

#  Usage: 

#     Recreate Fig 2.3 (upper panel) in F&M

#        evolve4(f.init=0.95, s.adv= -0.2, title="No dominance, neg. advantage", histogram=TRUE, color=TRUE) 

#       evolve4(f.init=0.05, s.adv=  0.2, title="No dominance, pos. advantage", histogram=TRUE, color=FALSE)  

#       evolve4(f.init=0.005, s.adv=  0.2, n.gen=150, title="No dominance, pos. advantage", histogram=TRUE, color=TRUE) 

#       Drift in 5 generations#     evolve4(n.pop1=100, n.pop2=100, n.gen=5, n.samp=1000, f.init=0.5, s.adv=0.0)

#       evolve4(f.init=0.5, s.adv=  0.1, n.pop1=200, n.pop2=200, n.gen=200, title="No dominance, pos. advantage", histogram=TRUE, color=TRUE, dom=FALSE) 

#      evolve4(f.init=0.5, s.adv=  -0.1, n.pop1=1000, n.pop2=1000, n.gen=200, title="No dominance, neg. advantage", histogram=TRUE, color=TRUE, dom=FALSE) 

#      evolve4(f.init=0.1, s.adv=  0.005, n.pop1=5000, n.pop2=5000, n.gen=200, title="", histogram=FALSE, color=TRUE, dom=FALSE) 

# n.pop1 = starting population size

# n.pop2 = ending population size

# 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

#  graphics.off()

#  setwd("/Users/paulreasenberg/Desktop/evolution/plots")

#

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 = as.integer( 2* (n.pop1 + ( (n.pop2-n.pop1) *j / (n.gen-1) ) ) )  

      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()

#