Demographic stochasticity in Malthusian populations

We denote by $N\left(t\right)$ the population abundance (more precisely, the number of females) at time $t$ keeping in mind that, as the population is small, $N$ must be considered as an integer, not a real number, because actually each organism is a discrete entity. It should also be noted that, although the population initially consists of exactly $N_0$ individuals, the abundance $N\left(t\right)$ at time $t>0$ is a random variable. The simplest question we can ask about this random variable is how its average changes over time. The average is defined as

   E$$\left[N(t)\right]=\overline N (t) = \sum\limits_{i = 0}^\infty {i\, {p_i}(t)} = \sum\limits_{i = 1}^\infty {i\, {p_i}(t)}$$

where $p_i\left(t\right)$ is the probability that at time $t$ the population consists of an integer number $i$ of individuals. The result that stems from the theory of stochastic processes (see e.g. Pielou, 1977) is extremely simple, because the average population size satisfies the equation

$$\frac{{d\overline N }}{{dt}} = \left(\nu - \mu\right)\overline N = r\overline N$$

which is the same equation that governs the dynamics of the deterministic Malthusian model. Therefore, it turns out

$$\overline N (t) = \overline N (0)\exp (rt) = {N_0}\exp (rt)$$

which implies that the average population size increases or decreases over time in accordance with the instantaneous rate of growth $r = \nu - \mu$ being positive or negative. Note that the average must be calculated from the infinite possible time evolutions (the realizations of the stochastic process) that can be followed by the population starting from the initial condition $N_0$.

However, one must not think that nothing changes with respect to the deterministic model. In particular, this result does not imply that the population can become extinct only if the growth rate $r$ is negative. After some non-elementary calculations starting from an infinite system of nonlinear differential equations (the so-called Kolmogorov equations) in the variables $p_i(t)$ - i.e. the probabilities that the population be composed of exactly $i=0, 1, 2, \ldots$ individuals - one gets the most important result of the theory of demographic stochasticity. This result concerns the dynamics over time $t$ of the probability $p_0 (t)$ that the population be composed of zero individuals, or in other words of the extinction risk. One can in fact prove Iannelli and Pugliese (2014) that the time dynamics of $p_0$ is given by the formula

$$p_0(t)= {\left( {\frac{{\mu \exp (rt) - \mu }}{{\nu \exp (rt) - \... ...{\frac{{\mu \exp (rt) - \mu }}{{(r + \mu )\exp (rt) - \mu }}} \right)^{{N_0}}},$$ (3.2)

where $N_0$ is the initial number of adult females (for simplicity we can assume that the sex ratio remains constant over time and does not affect the birth rate).

First of all, we note that the extinction probability is positive even for $r >$ 0. This is a fundamental result: demographic stochasticity can lead to extinction even Malthusian populations with a positive growth rate. Secondly, we can calculate the extinction probability of the population as $t$ tends to infinity. This calculation allows us to understand how likely the population is to die out in the long run. Obviously, the asymptotic extinction risk depends on the value of parameter $r$. Suppose first that $r <0$; one obtains

$${\lim _{t \to \infty }}{p_0}(t) = {\left( {\frac{{ - \mu }}{{ - \mu }}} \right)^{{N_0}}} =1 .$$

So in the long run extinction is certain for populations with negative growth rate. The result is indeed a foregone conclusion!

Suppose now that $r >$ 0. It is easy to get

$${\lim _{t \to \infty }}{p_0}(t) = {\lim _{t \to \infty }}{\left( ... ...{\nu }} \right)^{{N_0}}} = {\left( {\frac{\mu }{{r + \mu }}} \right)^{{N_0}}.}$$

Therefore, in the long run, there exists a non-vanishing extinction probability even for populations that on the average are growing, as shown in Fig. 9.

Figure 9: Two out of five stochastic simulations of the Malthusian process shown in figure terminate in extinction, although the growth rate is $r =$ 0.02 time-unit$^{-1}$. In bold, the trend of the theoretical average $\overline N(t)$

However, it is worth noting that a long time might be necessary before the extinction probability becomes close to the long-term value. Recalling that $\frac {1}{\mu}$ is the average lifetime of an organism belonging to the population - and is also the average generation length in the case of semelparous populations - we can wonder what the risk of extinction is after a reasonable time, for instance 10 or 50 generations. Calculations are easy thanks to eq. 7. Table 2 reports the results of computations that show that, for high mortality and small (albeit positive) growth rate, the probability of extinction after 50 generations is not very different from the asymptotic one. It should be noted that the extinction probability is anyway very small for populations with more than 50 individuals.

Table 2: Extinction probability for a Malthusian population characterized by $r = 0.05$ time-unit$^{-1}$ and $\mu = 0.5$ time-unit$^{-1}$ and initially composed by $N_0$ individuals: asymptotic probability (second column), probability after 10 and 50 generations (third and fourth column).

$N_0$	Asymptotic Prob.	10 gener. Prob.	50 gener. Prob.
1	0.909	0.863	0.908
5	0.621	0.480	0.619
10	0.386	0.230	0.383
50	0.00852	0.000647	0.00826
100	7.26x$10^{-5}$	4.18x$10^{-7}$	6.82x$10^{-5}$
500	2.01x$10^{-21}$	1.28x$10^{-32}$	1.48x$10^{-21}$

Finally, we can calculate the extinction probability for stationary populations ($r = 0$, namely they are stationary from a deterministic viewpoint). This is a bit critical, because we must resort to de L'H$\hat{o}$pital rule for indeterminate forms:

$${\lim_{r \to 0}{p_0}(t)} = {\lim _{r \to 0}}{\left( {\frac{{\mu t... ...}}} \right)^{{N_0}}} = {\left( {\frac{{\mu t}}{{1 + \mu t}}} \right)^{{N_0}}}.$$

We can then easily calculate the asymptotic extinction probability

$${\lim _{t \to \infty }}{p_0}(t) = {\lim _{t \to \infty }}{\left( {\frac{{\mu t}}{{1 + \mu t}}} \right)^{{N_0}}} = 1.$$

Therefore, long-term-extinction is certain in stationary populations, even if the average value of $N(t)$ is always constant and equal to the initial population $N_0$. This may sound like a paradox, but basically means this: if we replicated the demographic process very many times with initial population size equal to $N_0$, in almost all cases the population would become extinct sooner or later, but in some very rare cases the population would succeed in avoiding small numbers and grow exponentially thus keeping the average size constant. Note that, if $N_0$ is large, $p_0 (t)$ can tend to 1 very slowly (see Table 3). Anyway, after 1000 generations the extinction probability of a stationary population starting from 500 individuals is about 60%.

Table 3: Extinction probability at different times $t$ in stationary populations ($r = 0$, death-rate = $\mu$) for increasing initial population size $N_0$.

$N_0$	$t=1/\mu$	$t=10/\mu$	$t=100/\mu$	$t=1000/\mu$
1	0.5	0.909	0.990	0.999
5	0.031	0.621	0.951	0.995
10	0.00097	0.386	0.905	0.990
50	8.88x$10^{-6}$	0.0085	0.608	0.951
100	7.89x$10^{-31}$	0.000072	0.370	0.905
500	$\sim0$	2.01x$10^{-21}$	0.0069	0.607

It is also interesting to see how long it takes on average for a population to become extinct. A very simple index is the median time to extinction. As the time to extinction is a stochastic variable, then $p_0 (t)$ is also the probability that the time to extinction be $\leq t$. Therefore the median time to extinction $t_{\text{med}}$ is the one at which $p_0(t_{\text{med}}) = 0.5$. With easy calculations one gets for stationary populations ($r = 0$)

$${t_{{\rm {med}}}} = \frac{1}{\mu }\frac{1}{{{2^{\frac{1}{{{N_0}}}... ...{1}{{\frac{{\ln 2}}{{{N_0}}} + 1 - 1}} = \frac{1}{\mu }\frac{{{N_0}}}{{\ln 2}}.$$

So with $N_0=5$, the median extinction time of a stationary population is about 7 generations, with $N_0=50$ it is approximately 70 generations and with $N_0=500$ about 700 generations. If the population is non-stationary, one must distinguish between decreasing ($r <0$) and increasing ($r > 0$) populations. In decreasing populations, the asymptotic probability of extinction is one. Therefore, using the formula given by eq. 7 and equating the left-hand-side to 0.5, one can solve with respect to time and obtain a general expression for the median time to extinction for decreasing populations (the derivation is left to the reader as an exercise):

$${{t}_{med}}=\frac{1}{r}\ln \left( \frac{{{2}^{\frac{1}{{{N}_{0}}}}}-1}{{{2}^{\frac{1}{{{N}_{0}}}}}-1-\frac{r}{\mu }} \right).$$

With $r=-0.05$ time-unit$^{-1}$, the time to extinction is 6 generations for $N_0=5$, 30 generations for $N_0=50$, and 42 generations for $N_0 = 100$.

In increasing population, only the fraction $\left( {\frac{\mu }{\nu }} \right)^{{N_0}}$ of long-term possible time evolutions will become extinct. The remaining fraction will increase indefinitely. So the median time to extinction must be calculated conditional on population extinction. The probability for the population to become extinct within time $t$ conditional on extinction is ${{{p}_{0}}(t)}/{{{\left( {\mu }/{\nu }\; \right)}^{{{N}_{0}}}}}\;$. By equating this expression to $0.5$, one can obtain, after boring calculations, the following formula

$${{t}_{med}}=\frac{1}{r}\ln \left( \frac{{{2}^{\frac{1}{{{N}_{0}}}... ... {{2}^{\frac{1}{{{N}_{0}}}}}-1 \right)\left( 1+\frac{r}{\mu } \right)} \right).$$

With $r=0.02$ time-unit$^{-1}$, the conditional median time to extinction is 6 generations fo $N_0=5$, 44 generations for $N_0=50$, and 67 generations for $N_0 = 100$.