BIOS13:Demographic stochasticity

Demographic stochasticity (3p)

Small populations are vulnerable to extinction simply because they are small. The fate of single individuals may have large effect on the growth of the population. A few failed breeding attempts or fatal infections can cause a dramatic population decline. For large populations these effects tend to average out, such that the population as a whole has more predictable growth. This source of uncertainty is called demographic stochasticity. A model of population growth taking demographic stochasticity is described below.

Assume a (asexual) population of N individuals. Each year, each individual produces a few offspring and then dies. The offspring constitute the next generation, and so on. The number of offspring produced by a single individual follow a Poisson distribution with a mean of 1. Thus, each individual replaces itself on average, but there may be considerable variation between individuals. Some are lucky and get many offspring, some get none. Total population size N will in this manner vary over time and may even reach extinction, N = 0.

In R, the function rpois(N,lambda) generates N random numbers with mean lambda.

• Write a function in R that takes parameters N0 and tmax as input and returns the number of individuals in the population after tmax generations, starting with N0 individuals. (1p)

  generation = function(N0, Tmax){
    N = N0
    for(i in 1:Tmax){
      offspring = rpois(N, 1)
      N = sum(offspring)
    }
    return(N)
  }
  

• Write another function that takes N0 and tmax as input parameters and calls the function from a) repeatedly (1000 times) to calculate the risk of extinction within tmax generations, starting with N0 individuals. (1p)

  extinctionRisk = function(Ts, N0, Tmax){
    extinction = 0
    for(i in 1:Ts){
      N = generation(N0, Tmax)
      
      if(N == 0){
        extinction = extinction + 1
      }
    }
    return(extinction/Ts)
  }
  

• Finally, write a script that uses the function from b) to plot the risk of extinction within 100 generations as a function of N0, with N0 taking the values 10, 20, 30, … 100. (1p)

  risk = function(N){
    Ns = seq(10, N, by = 10)
    r = rep(0, length(Ns))
    for(i in 1:length(Ns)){
      r[i] = extinctionRisk(1000, Ns[i], 100)
    }
    return(r)
  }
  plot(seq(10, 100, by = 10), risk(100), type = 'b', xlim= c(10, 100), ylim = c(0, 1), xlab = 'Initial Population Size', ylab = 'Extinction Risk')