Image credit: James Maino

Unconstrained population growth

Let’s derive some population growth functions!

How to grow

Populations grow. They grow positively, if rates of reproduction > mortality, or negatively, if reproduction < mortality. For an unconstrained population of size \(N\) this diffence in per capita reproduction and mortality is referred to as the intrinsic growth rate, \(r\) and has the units individuals per individual per time \(N.N^{-1}.t^{-1}\) . The value of \(r\) is a constant if the age-distribution is constant (e.g. the proportion of reproductively active individuals in the population does not change). Taking \(r\) as constant, the change in population size with time is simply:

\[\frac{dN}{dt} = rN\]

This is a rudimentary first-order differential equation which can be solved in a few simple steps.

\[\int 1/N dN = \int r dt\] \[ln|N| = rt + C_0\] \[N = C_1 e^{rt}\]

Solving \(C_1\) for \(N = N_0\) at \(t = t_0\):

\[N = N_0 e^{rt}\]

While the parameter \(r\) hopefully makes conceptual sense, it is not straightforward to measure for a given species, as it will depend on the age structure (all juveniles means no reproduction), resource availability (food for growth), and environmental temperature (slow cold ectotherms). The age structure is assumed to be in a steady state, i.e. the relative proportion of different age classes does not change. Likewise, resource availability and other environmental conditions such as temperature are also taken as fixed. Steady state can be illustrated with a simple transition matrix model. Although there are initially 1000 individuals in the first age class (e.g. a mass hatching event triggered by rain), the age classes eventually stabilise.

trans = matrix(
  c(0.00, 0.50, 0.00,
    0.00, 0.00, 0.50,
    4.00, 0.00, 0.10),
  ncol = 3,byrow = T
  )
rownames(trans) = c('V1','V2','V3')
colnames(trans) = c('V1','V2','V3')
trans
##    V1  V2  V3
## V1  0 0.5 0.0
## V2  0 0.0 0.5
## V3  4 0.0 0.1
Tmax = 200
N = matrix(rep(0, Tmax*3), nrow = 3)
N[1,1] = 1000
for(t in 2:Tmax) N[, t] = N[, t-1] %*% trans

Measuring intrinsic growth rate

There are a number of ways that \(r\) can be measured, but they all rely on the above assumptions.

Regressing population vs. time

The most direct way to access \(r\) is to look at how an unconstrained growing population changes through time. Let’s use the numbers from our previous matrix model simulation.

df = data.frame(
      t = 0:(Tmax-1),
      N = colSums(N)
     )
knitr::kable(df[seq(0,200, by = 20),],digits = 0, align = c('r','r','r'))
t N
20 19 1177
40 39 2137
60 59 4041
80 79 8262
100 99 15973
120 119 31683
140 139 62267
160 159 122686
180 179 241627
200 199 475851

To estimate \(r\) through regression, notice that a log tranformation of the growth function results in:

\[\ln N = \ln (N_0 e^{rt})\] \[\ln N = \ln N_0 + rt\] This is a linear regression with an intercept \(\ln N_0\) and slope \(r\) which we can recover through fitting a simple linear model.

lm1 = lm(log(N) ~ t , data = df)
coef(lm1)
## (Intercept)           t 
##  6.30299077  0.03407767

Notice that the recovered intercept exp(coef(lm1)[1]) = 546 was actually lower than the real values of 1000. This was due to the age distribution not being in steady state, and the population not growing properly exponential.

Now let’s try and run the same matrix model simulation but with the 1000 individuals distributed across age classes at the steady state proportions.

N1_stable = 1000*N[,Tmax]/sum(N[,Tmax])
N = matrix(rep(0, Tmax*3), nrow = 3)
N[,1] = N1_stable 

for(t in 2:Tmax) N[, t] = N[, t-1] %*% trans

df = data.frame(
      t = 0:(Tmax-1),
      N = colSums(N)
     )

lm1 = lm(log(N) ~ t , data = df)
coef(lm1)
## (Intercept)           t 
##  6.90777973  0.03388836

After exponentiating the intercept we are able to recover the initial population size.

exp(coef(lm1)[1])
## (Intercept) 
##    1000.024

Measuring vital rates and building a life table

Old school population demographers also built life tables for a cohort of individuals followed through life, which include direct measurements of mortality, and reproduction to calculate \(r\).

Let’s use an example published by Carey and Bradley1 of the Atlantic Spider mite. The below table shows for each day \(x\) the proportion surviving \(P(x)\) and mean daily reproduction rate per individual for a synchronised cohort through time starting at egg lay.

x P(x) m(x)
0 1.00 0.00
1 1.00 0.00
2 1.00 0.00
3 1.00 0.00
4 1.00 0.00
5 1.00 0.00
6 1.00 0.00
7 1.00 0.37
8 1.00 3.53
9 1.00 7.43
10 0.96 7.24
11 0.85 4.31
12 0.73 6.57
13 0.68 5.59
14 0.50 4.86
15 0.46 4.70
16 0.38 3.67
17 0.33 4.26
18 0.29 3.93
19 0.25 2.48
20 0.23 3.83
21 0.21 5.09
22 0.19 3.30
23 0.15 2.34
24 0.09 1.65
25 0.04 1.10

Firstly, we estimate net generational reproduction \(R_0 = \sum_{x} P(x)m(x) = 42.8\) . As not all individuals will reproduce at \(P(x)\), i.e. some will die, this value is scaled by \(m(x)\). Now let’s plot mortality and reproduction data. But how fast does a population grow by this factor?

An analytical estimate of r from this data relies on the Lotka equation.

\[\sum_{x=\alpha}^\beta P(x)m(x)e^{-rx} = 1\]

Here we can set \(x = T\) where

\[T = \frac{\sum_{x}xP(x)m(x)}{\sum_{x}P(x)m(x)}\] Substituting into the Lotka equation and recognisign \(\sum_{x}P(x)m(x) = R_0\) we get :

\[e^{rT} \approx R_0\] \[r \approx \ln R_0/ T\]

Alternatively, we can solve this equation numerically.


  1. Carey, J.R. and Bradley, J.W., 1982. Developmental rates, vital schedules, sex ratios and life tables for Tetranychus urticae, T. turkestani and T. pacificus (Acarina: Tetranychidae) on cotton. Acarologia, 23(4), pp.333-345.

Related

comments powered by Disqus