Chapter 10 Analytical background

In this chapter I give a brief sketch of the computational approach, which is discussed in detail in Kirkilionis et al. (2001), Diekmann, Gyllenberg, and Metz (2003) and De Roos (2008). The description is far from complete, but only captures the basic idea of the computational machinery implemented in the package.

Consider the following generic model for the interaction of a size-structured consumer population foraging on an unstructured resource:

\[\begin{align*} &\dfrac{\partial n(t, s)}{\partial t}\:+\:\dfrac{\partial\left(g(s, R) n(t, s)\right)}{\partial s}\;=\;-\mu(s, R)\,n(t,s)\qquad\qquad\qquad\\[2ex] &g(s_b, R)\,n(t, s_b)\;=\;\int_{s_b}^{s_m} \beta(s, R)\,n(t,s)\,ds\\[2ex] &\dfrac{dR}{dt}\;=\;G(R)\;-\;\int_{s_b}^{s_m} \gamma(s, R)\,n(t, s)\,ds \end{align*}\]

In this model \(n(t,s)\) represents the size distribution of the consumer population at time \(t\) and \(R(t)\) is the resource density. The functions \(g(s,R)\), \(\beta(s,R)\) and \(\mu(s,R)\) represent the growth rate in size of an individual with size \(s\), its fecundity and its mortality rate, respectively. The function \(G(R)\) described the autonomous dynamics of the resource \(R\) in the absence of consumers.

The computational approach is based on the idea that this model can also be expressed as a system of integro-differential equations of the following form:

\[\begin{equation} \begin{aligned} b(t)&=\;\int_0^\infty\beta(s(t, a, R_t), R(t))\,\mathcal{F}(t, a, R_t)\,b(t-a)\,da\qquad\qquad\qquad\\[2ex] \dfrac{dR}{dt}&=\;G(R(t))\;-\;\int_0^\infty\gamma(s(t, a, R_t), R(t))\,\mathcal{F}(t, a, R_t)\,b(t-a)\,da \end{aligned}\tag{10.1} \end{equation}\]

in which \(b(t)\) is the population birth rate of the consumer population at time \(t\) and \(R_t\) represents the history of the resource density prior to time \(t\), i.e. the function \(R(\xi)\) with \(\xi\in(-\infty,t]\).

The function \(s(t, a, R_t)\) represents the body size of an individual consumer that is of age \(a\) at time \(t\) and has been exposed to the resource densities \(R_t\) since its birth. This body size is the integrated result of the growth rate \(g(s,R)\) that the individual has experienced since birth:

\[s(t, a, R_t)\;=\;s^0\;+\;\int_0^a g(s(t-\alpha, \alpha, R_{t-\alpha}), R(t-\alpha))\,d\alpha\]

The function \(\mathcal{F}(t, a, R_t)\) represents the probability that an individual that is of age \(a\) at time \(t\) and has been exposed to the resource densities \(R_t\) since its birth is still alive. \(\mathcal{F}(t, a, R_t)\) is related to the mortality rate \(\mu(s,R)\) following:

\[\mathcal{F}(t, a, R_t)\;=\;\exp\left(-\int_0^a \mu(s(t-\alpha, \alpha, R_{t-\alpha}), R(t-\alpha))\,d\alpha\right)\]

Figure 10.1 below illustrates how the integro-differential equation system relates the birth rate in the past to the birth rate at time \(t\) through the intervening history of the resource density and the development of the consumers that have experienced this resource history.

Figure 10.1: Schematic representation of the integro-differential equation system for the size-structured consumer-resource model, showing how the population birth rate at time \(t-a\) contributes to the birth rate at time \(t\) through consumers of age \(a\) that have grown during their life from their size at birth \(s_b\) till their current body size \(s(t,a,R_t)\) and have survived with a probability \(\mathcal{F}(t, a, R_t)\), which both depend on the history of the resource \(R_t\) that these consumers have experienced.

10.1 The system of equations determining the population growth rate

For demographic analysis of a linear PSPM only the integral equation (10.1) is relevant. In linear PSPMs the individual life history is not influenced by any density dependence or by any dependence on environment variables. We can hence drop the dependence of the development rate, fecundity and mortality rate on environment variables and generalize the integral equation (10.1) for an arbitrary choice of the individual state to:

\[b(t)\;=\;\int_0^\infty \beta(\boldsymbol{\chi}(a))\,\mathcal{F}(a)\,b(t-a)\,da\]

in which \(\boldsymbol{\chi}(a)\) is the state that individuals reach at age \(a\) provided they were born with state \(\boldsymbol{\chi}_b\). \(\boldsymbol{\chi}(a)\) is formally given by:

\[\boldsymbol{\chi}(a)\;=\;\boldsymbol{\chi}_b\;+\;\int_0^a g(\boldsymbol{\chi}(\alpha))\,d\alpha\]

and \(\mathcal{F}(a)\) is the probability of survival up to age \(a\):

\[\mathcal{F}(a)\;=\;\exp\left(\displaystyle -\int_0^a\mu(\boldsymbol{\chi}(\alpha))\,d\alpha\right)\]

Assuming exponential growth of the population birth rate:

\[b(t)=e^{ra}b(t-a)\]

leads to Lotka’s integral equation for the population growth rate \(r\):

\[\begin{equation} \int_0^\infty e^{-ra}\beta(\boldsymbol{\chi}(a))\,\mathcal{F}(a)\,da\;=\;1\tag{10.2} \end{equation}\]

Define the function \(H(a, r)\) as the value of Lotka’s integral up to age \(a\):

\[H(a,r)\;=\;\int_0^a e^{-r\alpha}\beta(\boldsymbol{\chi}(\alpha))\,\mathcal{F}(\alpha)\,d\alpha\]

Equation (10.2) can then be expressed as:

\[\begin{equation} H(\infty, r)\;=\;1\tag{10.3} \end{equation}\]

which is a non-linear equation for the population growth rate \(r\). This is the equation that is solved by the software package for the unknown quantity \(r\) using an iterative approach based on the Newton-Chord method. For more details of the Newton-Chord method I refer to Kuznetsov (1995), which source I have used to a large extent for the iterative calculation of the solution \(\tilde{r}\).

The central idea of the computational approach relates to the evaluation of the function \(H(\infty,r)\), which is computed by solving an ordinary differential equation. To derive this ODE differentiate \(\mathcal{F}(a)\) with respect to \(a\) using the chain rule:

\[\dfrac{d}{da}\mathcal{F}(a)\,=\,-\exp\left(\displaystyle -\int_0^a\mu(\boldsymbol{\chi}(\alpha))\,d\alpha\right)\dfrac{d}{da}\left(\int_0^a\mu(\boldsymbol{\chi}(\alpha))\,d\alpha\right)\]

Applying Leibniz rule for differentiation of an integral:

\[\dfrac{d}{d\theta}\left(\int_{a(\theta)}^{b(\theta)}\!\!f(\boldsymbol{\chi},\theta)d\boldsymbol{\chi}\right)= \int_{a(\theta)}^{b(\theta)}\!\!f_\theta(\boldsymbol{\chi},\theta)d\boldsymbol{\chi}+f(b(\theta),\theta)\,b^\prime(\theta)-f(a(\theta),\theta)\,a^\prime(\theta)\]

then leads to:

\[\dfrac{d\mathcal{F}}{da}\,=\,-\mu(\boldsymbol{\chi}(a))\,\mathcal{F}(a),\qquad\quad\mathcal{F}(0)\;=\;1\]

Similarly, differentiate \(H(a,r)\) with respect to \(a\) and applying Leibniz rule yields:

\[\dfrac{dH}{da}\,=\,e^{-ra}\,\beta(\boldsymbol{\chi}(a))\,\mathcal{F}(a),\qquad\quad H(0)\;=\;0\]

The value of \(H(\infty, r)\) can hence be calculated by (numerical) integration of the ODEs:

\[\left\{ \begin{aligned} \dfrac{d\boldsymbol{\chi}}{da}&=\;g(\boldsymbol{\chi}), &\boldsymbol{\chi}(0)\,=\,\boldsymbol{\chi}_b\\[1ex] \dfrac{d\mathcal{F}}{da}&=\;-\mu(\boldsymbol{\chi}(a))\,\mathcal{F}(a), &\mathcal{F}(0)\,=\,1\\[1ex] \dfrac{dH}{da}&=\;e^{-ra}\,\beta(\boldsymbol{\chi}(a))\,\mathcal{F}(a), &H(0)\,=\,0 \end{aligned}\right.\]

In practice numerical integration of these ODEs is carried out up to \(a=A_{max}\) with \(A_{max}\) either a fixed value or by \(\mathcal{F}(A_{max})=\epsilon\) with \(\epsilon\) a very small value (i.e. \(\epsilon=10^{-9}\)). Whenever an evaluation of the function \(H(\infty, r)\) is required in the Newton iterations of equation (10.3) this system of ODEs has to be integrated numerically. Once, a solution \(\tilde{r}\) has been found, the sensitivities of this solution with respect to the model parameters are calculated using numerical differentiation.

10.2 The system of equations determining an equilibrium

The idea discussed above for the demographic analysis of a linear PSPM extends to the computation of an equilibrium of a nonlinear PSPM. In such a nonlinear PSPM the fecundity and the development and mortality rates of individuals does depend on their environment, but in equilibrium this environment is necessarily constant: \(\tilde{E}\). Therefore, Lotka’s integral equation should determine as before the population growth rate \(r\):

\[\int_0^\infty e^{-ra}\beta(\boldsymbol{\chi}(a, \tilde{E}), \tilde{E})\,\mathcal{F}(a, \tilde{E})\,da\;=\;1\]

Note that all parts of life history now depend on \(E\). \(\boldsymbol{\chi}\) and \(\mathcal{F}\) do because of \(g(\boldsymbol{\chi},E)\) and \(\mu(\boldsymbol{\chi},E)\). However, \(r\) should equal \(0\) for equilibrium of the structured population:

\[\int_0^\infty \beta(\boldsymbol{\chi}(a, \tilde{E}), \tilde{E})\,\mathcal{F}(a, \tilde{E})\,da\;=\;1\]

In addition, the autonomous dynamics of the environment should be balanced by the impact of the population:

\[G(\tilde{E})\;=\;\tilde{b}\int_0^\infty \gamma(\boldsymbol{\chi}(a, \tilde{E}), \tilde{E})\,\mathcal{F}(a, \tilde{E})\,da\]

The survival rate \(\mathcal{F}(a, \tilde{E})\) and the value of the cumulative reproduction integral:

\[H(a, \tilde{E})\;=\;\int_0^a \beta(\boldsymbol{\chi}(\alpha,\tilde{E}), \tilde{E})\,\mathcal{F}(\alpha, \tilde{E})\,d\alpha\]

can be computed as before by solving the corresponding ODEs. To compute the impact of the population on the environment, define the function \(I(a, \tilde{E})\) as:

\[I(a, \tilde{E})\;=\;\int_0^a \gamma(\boldsymbol{\chi}(\alpha,\tilde{E}), \tilde{E})\,\mathcal{F}(\alpha,\tilde{E})\,d\alpha\]

\(I(a, \tilde{E})\) represents the cumulative, expected impact that a single individual exerts on its environment until age \(a\). Differentiating \(I(a, \tilde{E})\) with respect to \(a\) yields after applying Leibniz rule:

\[\dfrac{dI}{da}\,=\,\gamma(\boldsymbol{\chi}(a,\tilde{E}), \tilde{E})\,\mathcal{F}(a, \tilde{E}),\qquad\quad I(0)\;=\;0\]

The equilibrium of a nonlinear structured population model is therefore determined by the system of equations:

\[\begin{align*} &H(\infty, \tilde{E})\;=\;1\\[2pt] &\tilde{b}I(\infty, \tilde{E})\;=\;G(\tilde{E}) \end{align*}\]

which has to be solved numerically and iteratively for the unknowns \(\tilde{E}\) and \(\tilde{b}\). These equations are solved by the software package for the unknown quantities using the Newton-Chord method as discussed before. Whenever the functions \(H(\infty, \tilde{E})\) and \(I(\infty, \tilde{E})\) have to be evaluated in this iterative procedure, the following system of ODEs is integrated numerically:

\[\left\{ \begin{aligned} \dfrac{d\boldsymbol{\chi}}{da}&=\;g(\boldsymbol{\chi}(a, \tilde{E}), \tilde{E}), &\boldsymbol{\chi}(0, \tilde{E})\,=\,\boldsymbol{\chi}_b\\[0.5ex] \dfrac{d\mathcal{F}}{da}&=\;-\mu(\boldsymbol{\chi}(a, \tilde{E}), \tilde{E})\,\mathcal{F}(a, \tilde{E}),\qquad &\mathcal{F}(0, \tilde{E})\,=\,1\\[0.5ex] \dfrac{dH}{da}&=\;\beta(\boldsymbol{\chi}(a, \tilde{E}), \tilde{E})\,\mathcal{F}(a, \tilde{E}), &H(0, \tilde{E})\,=\,0\\[0.5ex] \dfrac{dI}{da}&=\;\gamma(\boldsymbol{\chi}(a, \tilde{E}), \tilde{E})\,\mathcal{F}(a, \tilde{E}), &I(0, \tilde{E})\,=\,0 \end{aligned} \right.\qquad\qquad\qquad\qquad\]

10.3 Curve continuation and detection of bifurcation points

The package uses the Newton-Chord method with Broyden updating of the Jacobian matrix to solve for the root of the nonlinear system of equations that determines the population growth rate of linear PSPMs or the equilibrium of nonlinear PSPMs. In addition, pseudo-arclength continuation is used to compute a curve of either the population growth rate or the equilibrium as a function of a single parameter. The numerical details about the Newton-Chord method as well as the pseudo-arclength continuation will not be discussed here. For details I refer to the appropriate sections in Kuznetsov (1995), which has been used as the basis for the implementations in the package. Both the Newton-Chord method as well as the pseudo-arclength continuation method make extensive use of partial derivatives of the system of equations with respect to variables and parameters. These partial derivatives, which for example make up the Jacobian matrix of the system of equations, are always computed numerically using a central-differencing approach.

The partial derivatives also play a role in the detection of bifurcation points, as explained in section 5.1. For example, the evolutionary analysis of PSPM using Adaptive Dynamics (AD) centers around the analysis of \(s_x(y)\), which is the population growth rate of a mutant with trait \(y\) in an environment that is completely determined by a resident population with trait \(x\) (Geritz et al. 1998). An evolutionary fixed point occurs at \(x=x^*\) where

\[\dfrac{\partial s_{x}(y)|_{x,y=x^{*}}}{\partial y}\;=\;0\]

The evolutionary fixed point can be classified as a convergent stable strategy (CSS), an evolutionary repellor (ERP) or evolutionary branching point (EBP) based on the value of Geritz et al. (1998):

\[\dfrac{\partial^2 s_{x}(y)|_{x,y=x^{*}}}{\partial y^2}\qquad\mathrm{and}\qquad \dfrac{\partial^2 s_{x}(y)|_{x,y=x^{*}}}{\partial x^2}\]

Because the equilibrium conditions for a structured model

\[\int_0^\infty \beta(\boldsymbol{\chi}(a, \tilde{E}), \tilde{E})\,\mathcal{F}(a, \tilde{E})\,da\,-\,1\;=\;R_0\,-\,1=0\]

is sign-equivalent with \(s_x(y)\) AD analysis can be performed using \(R_0\) and its (partial) derivatives with respect to resident and mutant traits \(x\) and \(y\), respectively (Geritz et al. 1998). Hence, the detection of evolutionary fixed points, their classification as convergent stable strategies, repellors or branching points, as well as the continuation of these evolutionary singularities as a function of two parameters, which is discussed in sections 5.1 to 5.3 relies on the computation of these partial derivatives, which are computed numerically as pointed out above.