In ecology, the Euler-Lotka equation describes the growth of a population in terms of $P(x)$, the fraction of individuals alive at age $x$ and $m(x)$, the mean number of live females born per time period per female alive during that time period: $$ \sum_{x=\alpha}^\beta P(x)m(x)e^{-rx} = 1, $$ where $\alpha$ and $\beta$ are the boundary ages for reproduction defining the discrete growth rate, $\lambda = e^r$. $r = \ln\lambda$ is known as Lotka's intrinsic rate of natural increase.
In a paper by Leslie and Ranson [1], $P(x)$ and $m(x)$ were measured for a population of voles (\emph{Microtus agrestis}) using a time period of 8 weeks. The data are given in the table below.
| $x$ /weeks | $m(x)$ | $P(x)$ |
|---|---|---|
| 8 | 0.6504 | 0.83349 |
| 16 | 2.3939 | 0.73132 |
| 24 | 2.9727 | 0.58809 |
| 32 | 2.4662 | 0.43343 |
| 40 | 1.7043 | 0.29277 |
| 48 | 1.0815 | 0.18126 |
| 56 | 0.6683 | 0.10285 |
| 64 | 0.4286 | 0.05348 |
| 72 | 0.3000 | 0.02549 |
The sum $R_0 = \sum_{x=\alpha}^\beta P(x)m(x)$ gives the ratio between the total number of female births in successive generations; a population grows if $R_0 > 1$ and $r$ determines how fast this growth is. In order to find $r$, Leslie and Ranson used an approximate numerical method; the code below determines $r$ by finding the real root of the Lotka-Euler equation directly (it can be shown that there is only one).
import numpy as np
from scipy.optimize import brentq
# The data, from Table 6 of:
# P. H. Leslie and R. M. Ranson, J. Anim. Ecol. 9, 27 (1940)
x = np.linspace(8, 72, 9)
m = np.array( [0.6504, 2.3939, 2.9727, 2.4662, 1.7043,
1.0815, 0.6683, 0.4286, 0.3000] )
P = np.array( [0.83349, 0.73132, 0.58809, 0.43343, 0.29277,
0.18126, 0.10285, 0.05348, 0.02549] )
# Calculate the product sequence f and R0, the ratio between the number of
# female births in successive generations.
f = P * m
R0 = np.sum(f)
if R0 > 1:
msg = 'R0 > 1: population grows'
else:
msg = 'Population does not grow'
# The Euler-Lotka equation: we seek the one real root in r
def func(r):
return np.sum(f * np.exp(-r * x)) - 1
# Bracket the root and solve with scipy.optimize.brentq
a, b = 0, 10
r = brentq(func, a, b)
print('R0 = {:.3f} ({})'.format(R0, msg))
print('r = {:.5f} (lambda = {:.5f})'.format(r, np.exp(r)))
The output of this program is as follows:
R0 = 5.904 (R0 > 1: population grows)
r = 0.08742 (lambda = 1.09135)
This value of $r$ may be compared with the approximate value obtained by Leslie and Ranson, who comment:
The required root is 0.087703 which slightly overestimates the value of $r$, to which the series is approaching. This lies between 0.0861 (the third degree approximation) and 0.0877, but nearer the latter than the former, the error being probably in the last decimal place.