M7116 Interactively
Mathematical background
More details
The matrix population model
\[ n(t + 1) = A n(t) \]
is essentially a vector linear recurrence relation, or a system of linear autonomous difference equations.. An analytical solution is given by
\[ n(t) = A^t n(0), \]
where \(n(0)\) is the vector of initial conditions. Rewriting this using the Jordan canonical form yields
\[ n(t) = W J^t W^{-1} n(0). \]
From Perron–Frobenius theory, we can derive the asymptotic behavior for certain matrices.
See Caswell (1989, 57)
Primitive Matrices
A primitive matrix \(A\) has a positive eigenvalue \(\lambda_1\) that is strictly greater in modulus than any other eigenvalue. This \(\lambda_1\) is called the dominant eigenvalue of \(A\). The associated right eigenvector \(w_1\) is positive is referred to as dominant. Similarly, \(v_1\) is the positive left eigenvector associated with \(\lambda_1\).
We can write:
\[ n(t) = \begin{pmatrix} w_1 & \widetilde{W} \end{pmatrix} \begin{pmatrix} \lambda_1^t & \mathbf{0}^T \\ \mathbf{0} & \widetilde{J}^t \end{pmatrix} \begin{pmatrix} v_1^T \\ \widetilde{V} \end{pmatrix} n(0) = \lambda_1^t w_1 v_1^T n(0) + \widetilde{W} \widetilde{J}^t \widetilde{V} n(0). \]
For large \(t\), we have:
\[ n(t) \sim c \lambda_1^t w_1, \]
where \(c = v_1^T n(0) > 0\).
The total population size at time \(t\) is given by the 1-norm of the population vector:
\[ N(t) = \| n(t) \|_1 = \sum_{i=1}^k n_i(t) \sim \lambda_1^t c \| w_1 \|_1. \]
Regardless of the initial population structure \(n(0)\), the population size behaves asymptotically like a geometric sequence with ratio \(\lambda_1\), and the population composition is proportional to the components of the eigenvector \(w_1\).
See Caswell (1989, 61)
Imprimitive and Irreducible Matrices
In this case, the matrix \(A\) has eigenvalues:
\[ \lambda_j = \lambda_1 e^{\frac{2(j-1)}{d} \pi i}, \quad j = 1, 2, \dots, d, \]
where \(d\) is the index of imprimitivity of the matrix \(A\), and \(\lambda_1\) is a positive real number. Other eigenvalues satisfy \(\lambda_1 > |\lambda|\).
We derive:
\[ n(t) = \lambda_1^t \sum_{j=1}^d c_j e^{\frac{2(j-1)}{d} t \pi i} w_j + \widetilde{W} \widetilde{J}^t \widetilde{V} n(0), \]
where \(c_j = v_j^T n(0)\). Thus, the solution of the projection equation with an irreducible and imprimitive matrix \(A\) is asymptotically equivalent to the vector sequence: \[ n(t) \sim \lambda_1^t \sum_{j=1}^d c_j e^{\frac{2(j-1)}{d} t \pi i} w_j. \]
Defining
\[ s(t) = \sum_{j=1}^d c_j e^{\frac{2(j-1)}{d} t \pi i} w_j, \]
yields:
\[ n(t) \sim \lambda_1^t s(t). \]
Also,
\[ s(t + d) = s(t), \quad \frac{1}{d} \sum_{l=0}^{d-1} s(t + l) = \frac{c_1}{d} w_1. \]
This result implies that the solution behaves like a geometric sequence with ratio \(\lambda_1\), modulated by a \(d\)-periodic structure.
If all eigenvectors \(w_j\) are normalized, the total population size satisfies:
\[ N(t) \sim \lambda_1^t \left\| \sum_{j=1}^d c_j e^{\frac{2(j-1)}{d} t \pi i} w_j \right\|_1. \]
See Caswell (1989, 63)
Damping ratio
Let \(\mu\) be the largest modulus of the eigenvalues of the \(k \times k\) matrix \(A\) that are smaller than \(\lambda_1\); in the case \(d = k\), where \(d\) is the imprimitivity index, we define \(\mu = 0\).
Then the damping ratio is defined as:
\[ \rho = \frac{\lambda_1}{\mu}, \]
and for \(\mu = 0\) we set \(\rho = \infty\).
The damping ratio expresses the rate of convergence to the stabilized structure. See Caswell (1989, 69)
Sensitivity of the growth rate
The sensitivity of the characteristic \(\chi = \chi(A)\) to the component \(a_{ij}\) of the projection matrix \(A\) is defined as
\[ \frac{\partial \chi}{\partial a_{ij}}. \]
In particular, for the growth rate \(\lambda_1\), if matrix \(A\) is ireducible, we have
\[ \frac{\partial \lambda_1}{\partial a_{ij}} = \frac{v_i w_j}{v^T w}. \]
The sensitivity matrix \(S\) is then defined by
\[ S = (s_{ij}) = \left( \frac{\partial \lambda_1}{\partial a_{ij}} \right) = \left( \frac{v_i w_j}{v^T w} \right) = \frac{1}{v^T w} v w^T. \]
See Caswell (1989, 118)
Elasticity of the growth rate
The elasticity of the characteristic \(\chi = \chi(A)\) with respect to the component \(a_{ij}\) of the projection matrix \(A\) is defined as
\[ \frac{a_{ij}}{\chi} \cdot \frac{\partial \chi}{\partial a_{ij}} = \frac{\partial \ln \chi}{\partial \ln a_{ij}}. \]
In particular, the elasticity \(e_{ij}\) of the growth rate \(\lambda_1\), if matrix \(A\) is ireducible, with respect to the entry \(a_{ij}\) is
\[ e_{ij} = \frac{1}{\lambda_1 v^T w} a_{ij} v_i w_j. \]
The elasticity matrix is given by
\[ E = (e_{ij}) = \frac{1}{\lambda_1 v^T w} A \circ (v w^T), \]
where \(\circ\) denotes the Hadamard (element-wise) product of matrices.
See Caswell (1989, 132)
Long-term behavior of eigenvalues for Caswell’s \(2 \times 2\) matrix model
The long-term behavior of the population vector \(n(t)\), given by the equation
\[ n(t) = A^t n(0), \]
depends on the eigenvalues \(\lambda_i\) of the projection matrix \(A\), as they are raised to higher and higher powers. For this particular model, we have two real eigenvalues. In this situation, we can in some way analyze the behavior of this system given by matrix \(A\) using all of its eigenvalues.
(In more general cases, the contribution of all the smaller eigenvalues to the behavior tends to be less relevant, and the deciding factors are the dominant eigenvalue and whether the matrix is primitive or not, and at least irreducible.)
If \(\lambda_i \in \mathbb{R}\), then:
- If \(\lambda_1 > 1\): the population grows exponentially.
- If \(\lambda_1 = 1\): the population stabilizes over time.
- If \(0 < \lambda_2 < 1\): the rate of convergence to the stable structure increases as \(\lambda_2\) tends to zero.
- If \(\lambda_2 = 0\): the population stabilizes after just one step.
- If \(-1 < \lambda_2 < 0\): the population stabilizes with damped oscillations.
- If \(\lambda_2 = -1\): the ratio of the components of the solution \(n(t)\) changes periodically, with a period of 2.
- If \(0 < \lambda_1 < 1\): the population decays exponentially.
Projecton matrix of Caswell two stage model:
\[ \begin{pmatrix} \sigma_1 (1-\gamma) & \varphi \\ \sigma_1 \gamma & \sigma_2 \\ \end{pmatrix} \]
Caswell two stage model
Here is the table for the reader to try some of the aforementioned situations. The first four matrices are primitive, and the last two are imprimitive and irreducible. For all examples, the fertility \(\varphi = 1\) and the initial conditions vector is \(n(0) = (1, 0)^T\).
| \(\sigma_1\) | \(\sigma_2\) | \(\gamma\) | Matrix \(A\) | \(\lambda_1\) | \(\lambda_2\) | Behavior of population \(||n(t)||_1\) | Damping ratio |
|---|---|---|---|---|---|---|---|
| \(1\) | \(1\) | \(1\) | \(\pmatrix{0 & 1 \\ 1 & 1}\) | \(\frac {1 + \sqrt{5}} 2\) | \(\frac{1 - \sqrt{5}} 2\) | \(\lambda_1 > 1\) : population grows exponentially | \(\frac{\sqrt{5} - 1} 2\) |
| \(\frac 2 3\) | \(\frac 1 3\) | \(1\) | \(\pmatrix{\frac 1 3 & 1 \\ \frac 2 3 & 1}\) | \(1\) | \(\frac 1 3\) | \(\lambda_1 = 1\), \(0 < \lambda_2 < 1\): population stabilizes over time | \(3\) |
| \(\frac 2 3\) | \(\frac 1 3\) | \(\frac 1 2\) | \(\pmatrix{\frac 1 2 & 1 \\ \frac 1 2 & 1}\) | \(1\) | \(0\) | \(\lambda_1 = 1\), \(\lambda_2 = 0\): population stabilizes after one step | \(\infty\) |
| \(\frac 5 9\) | \(1\) | \(1\) | \(\pmatrix{0 & 1 \\ \frac 5 9 & 1}\) | \(1\) | \(-\frac 5 9\) | \(\lambda_1 = 1\), \(\lambda_2 < 0\): population stabilizes, but exhibits damped oscillations | \(\frac 9 5\) |
| \(1\) | \(0\) | \(1\) | \(\pmatrix{0 & 1 \\ 1 & 0}\) | \(1\) | \(-1\) | \(\lambda_1 = 1\), \(\lambda_2 = -1\): population is stable, but its structure oscillates with period \(2\) | \(\infty\) |
| \(\frac 2 3\) | \(0\) | \(1\) | \(\pmatrix{0 & 1 \\ \frac 2 3 & 0}\) | \(\frac {\sqrt 2} 3\) | \(-\frac {\sqrt 2} 3\) | \(|\lambda_1| < 1\), \(\lambda_2 = \overline{\lambda_1}\): population decreases, but its structure oscillates with period \(2\) | \(\infty\) |
Model with constant Matrix
| Model | Matrix |
|---|---|
| Leslie model | \(\pmatrix{ 0 & 1 & 5 \\ 0.3 & 0 & 0 \\ 0 & 0.5 & 0}\) |
| Reducible-Improper model | \(\pmatrix{ 0 & 0 & 0 & 0 \\ 0.1 & 0 & 0 & 25 \\ 0 & 0.3 & 0 & 0 \\ 0 & 0 & 0.2 & 0}\) |
| Killer Whale model | \(\pmatrix{ 0 & 0.0043 & 0.1132 & 0 \\ 0.9775 & 0.9111 & 0 & 0 \\ 0 & 0.0736 & 0.9534 & 0 \\ 0 & 0 & 0.0452 & 0.9804}\) |
| Teasel model | \(\pmatrix{ 0 & 0 & 0 & 0 & 0 & 322.38 \\ 0.966 & 0 & 0 & 0 & 0 & 0 \\ 0.013 & 0.010 & 0.125 & 0 & 0 & 3.448 \\ 0.007 & 0 & 0.125 & 0.238 & 0 & 30.170 \\ 0.001 & 0 & 0.036 & 0.245 & 0.167 & 0.863 \\ 0 & 0 & 0 & 0.023 & 0.750 & 0}\) |
Interactive matrix dynamics: general (n x n) projection
Select the size of the projection matrix, define its entries, and explore the resulting dynamics of the structured population. For visualisation of the elasticity and sensitivity matrices, we chose two colour scales. Grey for \(i,j\) , where \(A_{ij} = 0\) and green or orange for \(A_{ij} > 0\), where \(A\) is the projection matrix. More vibrant color indicates higher value. Keep in mind, vibrant of the colours are scaled according to the maximum values in the plotted senistivity/elasticity matrix, so be careful in comparing these visualisation across different models.
Select matrix size (n × n)
Define the entries of matrix A
Length of projection (T)
Models on a continuous spatial domain
Fisher equation
In mathematics, Fisher equation (named after Ronald Fisher) also known as the KPP equation (named after Andrey Kolmogorov, Ivan Petrovsky, and Nikolai Piskunov) or Fisher-KPP equation (Fisher 1937) is the partial differential equation: \[ \frac{\partial u}{\partial t} = D\frac{\partial^2 u}{\partial x^2} + ru(1-u), \qquad x\in\mathbb{R},\quad t>0. \] It is a kind of reaction-diffusion system that can be used to model population growth and wave propagation.
If we suppose the following initial conditions: \[ u(0,x) = 0.01 \delta(x), \quad x\in\mathbb{R}. \]
We can find a solution in a form of travelling wave.
The travelling wave solution can be written as \[ u(t,x) = v(x - ct), \] for all \(x\in\mathbb{R}\), \(t>0\) and for some \(c>0\). The wave propagates with a constant velocity \(c\). However, the wave speed must be greater than \(2\) (or rather \(c\geq2\sqrt{rD}\)). Furthermore, for this solution, it holds \[ \lim_{t\to\infty} u(t,0) = 1, \qquad \lim_{x\to\infty} u(t,x) = 0,\quad \forall t>0. \] That is, the solution switches from the equilibrium state \(u = 0\) to the equilibrium state \(u = 1\). No such solution exists for \(c < 2\). For further reading we refer the reader to (Canosa 1973).
Reaction-diffusion equation
The problem above can be generalized to any reaction function, yielding the reaction-diffusion equation: \[ \frac{\partial u}{\partial t} = \frac{\partial}{\partial x}D\frac{\partial u}{\partial x} + f(u), \qquad x\in\mathbb{R},\quad t>0. \] In this equation \(u = u(t,x)\) represents the population density at location \(x\) at time \(t\), i.e. population size on interval \((a,b)\) at time \(t\) equals \(\int_a^bu(t,x)dx\); \(D = D(x)\) is the diffusion coefficient and \(f\) represents density dependent intensity of “demographic events”, i.e. intensity of “the new individuals production” (birth) and of “the old ones removal” (death).
Here the individuals constituting the modeled population can appear, die and move at any time. However, demographic events do not always happen on a continuous scale, as can be demonstrated, for example, by mating seasons in nature.
Time-discrete reaction-dispersion equation
In time-discrete analogy, we need to modify our assumptions. The life cycle of the modeled population consists of two phases \(1.\) demographic phase and \(2.\) dispersion phase. Demographic events occur at time instants \(t=t_1, t_2,\dots\) and individuals disperse at the remaining time. Therefore, our model satisfies \[ \begin{aligned} \frac{\partial u}{\partial t} &= \frac{\partial}{\partial x}D\frac{\partial u}{\partial x}, \qquad x\in\mathbb{R},\ t\in\left(t_i, t_{i+1}\right),\ i=1,2,\dots\ ,\\ \lim_{t\to i^+}u(t,x) &= f(\lim_{t\to i^-}u(t,x)), \qquad x\in\mathbb{R},\ i=1,2,\dots\ . \end{aligned} \] This is a well defined system, that for a given initial condition and boundary conditions has a solution. However, if one wants, this can be solved with theory from calculus, namely, with the Green function. It then can be derived that this equality holds \[ u(t+1, x) = \int_{-\infty}^\infty k(x-y)f(u(t,y))dy = k(x) \ast f(u(t,x)), \qquad t=1,2,\dots\ , \] where \(\ast\) denotes convolution and \(k(x) = \frac{1}{2\sqrt{\pi D}}e^{-\frac{x^2}{4D}}\) or a kernel in general.
We could take this further and compute the minimal wave speed for given kernel and reaction term, however, that would take some fun from the simulations, and we will just state the result: \[ c_{\min} = \min_{\mu>0}\left\{ \frac{1}{\mu}\ln\left( f'(0)\int_{-\infty}^\infty k(s)e^{\mu s} ds \right) \right\}. \] For further reading we refer reader to (Kot 1992).
Let us now consider some specific examples.
- Beverton-Holt reaction term: \[ f(u) = \frac{ru}{1 + (r - 1)u/K}. \]
- Logistic reaction term (adjusted for \(K\) being the carrying capacity): \[ f(u) = (1 + r)u - \frac{r}{K}u^2. \]
- Ricker reaction term \[ f(u) = ue^{r(1 - \frac{u}{K})}. \]
Interactive simulation
Considering everything what was mentioned above, the reader may make their own simulation.
We suggest to start with lower growth rate (\(r\)) and then proceed to higher values. One might notice regular and irregular patterns, presumably, induced by chaotic behavior of the reaction term.
Parameters of the simulation
Output of the simulation
Bibliography