The size-structured PDE
theory
This page outlines the demographic machinery used to model size-structured populations in the plant package. We outline the system of equations being solved for individual plants, patches of competing plants, and the age structure of patches across a metacommunity. These equations were primarily developed elsewhere, in particular in Falster, Brännström, Dieckmann, & Westoby (2011) and Falster, Brännström, Westoby, & Dieckmann (2015), following general principles laid out in De Roos, Taljapurkar, & Caswell (1997), Kohyama (1993), and Moorcroft, Hurtt, & Pacala (2001). They are presented here so that users can understand the full system of equations being solved within plant.
The companion page on adaptive dynamics builds on this machinery to define invasion fitness and evolutionary endpoints.
Table 1 provides a list of names and definitions used throughout this page.
plant package.
| Symbol | Unit | Description |
|---|---|---|
| Patch state variables | ||
| \(K\) | Number of species in the metacommunity | |
| \(a, a^{\prime}\) | yr | Patch age (time since disturbance) |
| \(a_0\) | yr | Patch age when plant germinates |
| \(E_a\) | Profile of canopy openness within a patch of age \(a\) | |
| \(E_a(z)\) | Canopy openness at height \(z\) within a patch of age \(a\) | |
| Plant state variables | ||
| \(x\) | Vector of traits for a species | |
| \(H\) | m | Height of a plant |
| \(H_0(x)\) | m | Height of a seedling with traits \(x\) after germination |
| \(z\) | m | Height in canopy |
| \(A_{\rm l}(H)\) | m\(^2\) | Leaf area of a plant with height \(H\) |
| \(Q(z, H)\) | Fraction of leaf area above height \(z\) for a plant with height \(H\) | |
| Abundance measures | ||
| \(P(a)\) | yr\(^{-1}\) | Frequency-density of patches of age \(a\) |
| \(Y_{x}\) | m\(^{-2}\) yr\(^{-1}\) | Global birth rate for a species with traits \(x\) |
| \(N(H, x, a)\) | m\(^{-1}\) m\(^{-2}\) | Density of plants at height \(H\) per unit ground area for given traits \(x\) and patch age \(a\) |
| Demographic rates | ||
| \(g(x, H, E_a)\) | m yr\(^{-1}\) | Height growth rate of a plant with traits \(x\) and height \(H\) in the light environment \(E_a\) in a patch of age \(a\) |
| \(f(x, H, E_a)\) | yr\(^{-1}\) | Seed production rate of a plant with traits \(x\) and height \(H\) in the light environment \(E_a\) in a patch of age \(a\) |
| \(d(x, H, E_a)\) | yr\(^{-1}\) | Instantaneous mortality rate of a plant with traits \(x\) and height \(H\) in the light environment \(E_a\) in a patch of age \(a\) |
| \(d_{\rm P}(a)\) | yr\(^{-1}\) | Instantaneous disturbance rate of a patch of age \(a\) |
| Demographic outcomes | ||
| \(H(x, a_0, a)\) | m | Height of a plant with traits \(x\) in a patch of age \(a\) having germinated at patch age \(a_0\) |
| \(S_{\rm D}\) | Probability a seed survives dispersal | |
| \(S_{\rm G} (x, H_0, E_{{\rm a}0})\) | Probability the seed of a plant with traits \(x\) germinates successfully at height \(H_0\) in the light environment \(E_{{\rm a}0}\) in a patch of age \(a_0\) | |
| \(S_{\rm I} (x, a_0, a)\) | Probability a plant survives from patch age \(a_0\) to patch age \(a\) | |
| \(S_{\rm P} ( a_0, a)\) | Probability a patch remains undisturbed from patch age \(a_0\) to patch age \(a\) | |
| \(\tilde{R}(x, a_0, a)\) | Cumulative seed output for plants with traits \(x\) from patch age \(a_0\) to patch age \(a\) | |
| \(R\left(x^\prime, x\right)\) | Basic reproduction ratio for a mutant plant with traits \(x^\prime\) growing in a metapopulation of resident plants with traits \(x\) | |
| Miscellaneous constants | ||
| \(k_{\rm I}\) | Light extinction coefficient | |
| \(\bar{a}\) | yr | Mean interval between patch disturbances |
Metacommunity dynamics
The following material is written assuming there exists a physiological sub-model that takes as inputs a plant’s traits \(x\), height \(H\), and light environment \(E_a\), and returns rates of growth, mortality, and fecundity. The physiological model must describe how much leaf area the plant contributes to shading other plants in the patch. Specifically, the physiological model is responsible for calculating the following variables from Table 1: \(A_{\rm l}(H)\), \(Q(z, H)\), \(H_0(x)\), \(g(x, H, E_a)\), \(f(x, H, E_a)\), \(d(x, H, E_a)\), and \(S_{\rm G} (x, H_0, E_{{\rm a}0})\). All are other variables in Table 1 are calculated in the demographic model using the equations provided below.
Individual plants
We first consider the dynamics of an individual plant. Throughout, we refer to a plant as having traits \(x\) and size \(H\), with the latter given by its height. The plant grows in a light environment \(E\), a function describing the distribution of light with respect to height within a patch. Ultimately, \(E\) depends on the composition of plants in a patch, and thus on the patch age \(a\). To indicate this dependence, we write \(E_a\). Further, let \(a_0\)be the age in which the plant germinated, and \(H_0(x)\) its initial height. The functions \(g(x, H, E_a)\), \(f(x, H, E_a)\), and \(d(x, H, E_a)\) denote the growth, death, and fecundity rates of the plant. Then, \[ H(x, a_0, a) = H_0(x) + \int_{a_0}^{a} g(x, H(x, a_0, a^\prime), E_{a^\prime}) \, {\rm d}a^\prime \tag{1}\] is the trajectory of plant height, \[ S_{\rm I} (x, a_0, a) = S_{\rm G} (x, H_0, E_{{\rm a}0}) \, \exp\left(- \int_{a_0}^{a} d(x, H(x, a_0, a^\prime), E_{a^\prime}) \, {\rm d} a^\prime \right) \tag{2}\] is the probability of plant survival within the patch, where \(S_{\rm G} (x, H_0, E_{{\rm a}0})\) is the probability a seed germinates successfully, and \[ \tilde{R}(x, a_0, a) = \int_{a_0}^{a} f(x, H(x, a_0, a^\prime), E_{a^\prime}) \, S_{\rm I} (x, a_0, a^\prime) \, {\rm d} a^\prime \tag{3}\] is the plant’s cumulative seed output, from its birth at patch age \(a = a_0\) until patch age \(a\).
The notational complexity required in Equation 13 potentially obscures an important point: these equations are simply the general, non-linear solutions to integrating growth, mortality, and fecundity over time.
Patches of competing plants
We now consider a patch of competing plants. At any age \(a\), the patch is described by the size-density distribution \(N(H | x, a)\) of plants with traits \(x\) and height \(H\). In a finite-sized patch, \(N\) is described by a collection of points, each indicating the height of one individual, whereas in a very (infinitely) large patch, \(N\) is a continuous function. In either case, the demographic behaviour of the plants within the patch is given by Equation 13. Integrating the dynamics over time is complicated by- two other factors: (i) plants interact, thereby altering \(E_a\) with age, and (ii) new individuals may establish, expanding the system of equations.
In the current version of plant, plants interact by shading one another. Following standards biophysical principles, we let the canopy openness \(E_a(z)\) at height \(z\) in a patch of age \(a\) decline exponentially with the total amount of leaf area above \(z\), \[
E_a(z) = \exp \left(-k_{\rm I} \sum_{i = 1}^{K} \int_{0}^{\infty} A_{\rm l}(H) \, Q(z, H) \, N(H | x_i, a) \, {\rm d}H \right),
\tag{4}\] where \(A_{\rm l}(H)\) is the total leaf area of a plant with height \(H\), \(Q(z, H)\) is the fraction of this leaf area situated above height \(z\) for plants of height \(H\), \(k_{\rm I}\) is the light extinction coefficient, and \(K\) is the number of species in the metacommunity.
Assuming patches are sufficiently large, the dynamics of \(N\) can be modelled deterministically via the following partial differential equation (PDE) (Kohyama, 1993; De Roos et al., 1997; Moorcroft et al., 2001), \[ \frac{\partial}{\partial a} N(H | x, a) = - d(x, H, E_a) \, N(H | x, a) - \frac{\partial}{\partial H} \left[g(x, H, E_a) \, N(H | x, a)\right]. \tag{5}\] See Section 2.2 for the derivation of this PDE.
Equation 5 has two boundary conditions. The first links the flux of individuals across the lower bound \((H_0(x))\) of the size-density distribution to the rate \(Y_{x}\) at which seeds arrive in the patch, \[ N(H_0 | x, a_0) = \left\{ \begin{array}{ll} \frac{Y_{x} \, S_{\rm G} (x, H_0, E_{{\rm a}0}) }{ g(x, H_0, E_{{\rm a}0}) } & \textrm{if } g(x, H_0, E_{{\rm a}0}) > 0 \\ 0 & \textrm{otherwise.} \end{array} \right. \tag{6}\] The function \(S_{\rm G} (x, H_0, E_{{\rm a}0})\) denotes survival through germination and must be chosen such that \(S_{\rm G} (x, H_0, E_{{\rm a}0}) / g(x, H_0, E_{{\rm a}0}) \rightarrow 0\) as \(g(x, H_0, E_{{\rm a}0}) \rightarrow 0\), to ensure a smooth decline in initial density as conditions deteriorate (Falster et al., 2011).
The second boundary condition of Equation 5 specifies the initial size-density distribution in patches right after a disturbance, i.e., for \(a = 0\). Throughout, we consider only situations starting with an empty patch, \[ N\left(H|x,0\right) = 0, \tag{7}\] although non–zero distributions could alternatively be specified (e.g. Kohyama, 1993; Moorcroft et al., 2001).
Age structure of patches
We now consider the distribution of patch age \(a\) in the metacommunity. With \(a\) denoting the time since last disturbance, \(P(a)\) denotes the frequency-density of patch age \(a\) and \(d_{\rm P}(a)\) is the age-dependent probability that a patch of age \(a\) is transformed into a patch of age 0 by a disturbance event. Here we focus on situations where the age structure has reached an equilibrium. See Section 2.1 for the derivation and non-equilibrium situations. The dynamics of \(P\) are given by (McKendrick, 1926; Foerster, 1959) \[ \frac{{\rm d}}{{\rm d} a} P(a) = - d_{\rm P}(a) \, P(a) , \tag{8}\] with the boundary condition \[ P(0) = \int_0^\infty d_{\rm P}(a) \, P(a) \, {\rm d} a. \tag{9}\] The probability that a patch remains undisturbed from patch age \(a_0\) to patch age \(a\) is then given by \[ S_{\rm P} (a_0, a) = \exp\left( - \int_{a_0}^{a} d_{\rm P}(a^\prime) \, {\rm d} a^\prime \right). \tag{10}\] These equations lead to an equilibrium distribution of patch ages, \[ P(a) = P(0) S_{\rm P} (0, a), \tag{11}\] where \[ P(0) = \frac1{\int_0^\infty S_{\rm P} (0, a) {\rm d}a} \tag{12}\] is the average disturbance frequency of a patch and, at the same time, the frequency-density of patches of age \(0\).
The default approach in plant is to assume that \(d_{\rm P}(a)\) increases linearly with patch age, which leads to a Weibull distribution for \(P\), specified by a single parameter \(\bar{a}\) measuring the mean interval between disturbances (see Section 2.1.2 for details).
Trait-, size-, and patch-structured metacommunities
We consider a large area of habitat where: (i) disturbances (such as fires, storms, landslides, floods, or disease outbreaks) strike patches of the habitat on a stochastic basis, killing all individuals within the affected patches; (ii) individuals compete for resources within patches, with a spatial scale of competitive interactions that renders negligible such interactions between individuals in adjacent patches; and (iii) there is high connectivity via dispersal between all patches in the habitat, allowing empty patches to be quickly re-colonised. Such a system can be modelled as a metapopulation (for a single species) or metacommunity (for multiple species). The dynamics of such a metacommunity are described by the PDEs in Equation 512.
The offspring production of each species in the metacommunity is given by the rate at which seeds are produced across all patches, \[ Y_x = \int_0^{\infty} \int_0^{\infty} P(a) \, S_{\rm D} \, f(x, H, E_a) \, N(H | x, a) \, {\rm d} H \, {\rm d} a, \tag{13}\] where \(S_{\rm D}\) is the survival probability of seeds during dispersal.
A convenient feature of Equation 5 - Equation 12 is that the dynamics of a single patch scale up to give the dynamics of the entire metacommunity. Note that the rate \(Y_x\) at which offspring of individuals with traits \(x\) arrive from the disperser pool is constant when the metacommunity is at equilibrium. Combined with the assumption that all patches have the same initial (empty) size-density distribution after a disturbance, the assumption of constant birth rates ensures that all patches show the same temporal behaviour, the only difference between them being their ages.
To model the temporal dynamics of an archetypal patch, we need only a value for \(Y_x\). The numerical challenge is therefore to find the right value for \(Y_x\), by solving Equation 613 simultaneously, for all species with traits \(x\) in the metacommunity.
Emergent properties of metacommunities
Summary statistics of a metacommunity are obtained by integrating over the size-density distribution, weighting by the frequency-density \(P(a)\). In particular, the average density of individuals per unit ground area across the metacommunity is and the average size-density of plants with height \(H\) is \[ \bar{N}(x, H) = \int_{0}^{\infty}P(a) \, N(H | x, a) \, {\rm d}a. \tag{15}\]
Averages for other individual-level quantities can also be calculated. Denoting by \(w(x, H, E_a)\) a quantity of interest, either a demographic rate (growth, mortality) or a state (plant height, leaf area, light environment), the average of \(w\) for plants with height \(H\) and traits \(x\) is \[ \bar{w}(x, H) = \frac1{\bar{N}(x, H)}\int_{0}^{\infty}P(a) \, N(H | x, a) \, w(x, H, E_a) \, {\rm d}a. \tag{16}\] The average of \(w\) across all individuals of the species is \[ \hat{w}(x) = \frac1{\hat{N}(x) }\int_{0}^{\infty} \int_{0}^{\infty}P(a) \, N(H | x, a) \, w(x, H, E_a) \, {\rm d}a \, {\rm d}H. \tag{17}\] When calculating the average mortality rate, one must decide whether mortality due to patch disturbance is included. Non-disturbance mortality is obtained by setting \(w(x, H, E_a) = d(x, H, E_a)\), while the total mortality due to growth processes and disturbance is obtained by setting \(w(x, H, E_a) = d(x, H, E_a) + d_{\rm P}(a) S_{\rm P}(0, a).\)
Similarly, we can integrate over the size-density distribution to extract aggregate features of the vegetation within a patch, \[ W(a) = \sum_{i = 1}^{N} \int_{0}^{\infty} N(H | x_i, a)\, w(x_i, H, E_a) \, {\rm d}H, \tag{18}\] and across the entire metacommunity, \[ \hat{W} = \int_{0}^{\infty} P(a) \, W(a) \, {\rm d}a. \tag{19}\]
For the analysis of invasion fitness and evolutionary endpoints that builds on these dynamics, see the adaptive dynamics page.
Appendices
Derivation of PDE describing age-structured dynamics
We consider patches of habitat that are subject to an intermittent disturbance, with the age of a patch measuring the time since the last disturbance. Denoting by \(P(a, t)\) the frequency-density of patch age \(a\) at time \(t\) and by \(d_{\rm P}(a)\) the age-dependent probability that a patch of age \(a\) is transformed into a patch of age \(0\) through disturbance, the dynamics of \(P(a, t)\) are given by \[ \frac{\partial}{\partial t} P(a, t) = -\frac{\partial}{\partial a} P(a, t)-d_{\rm P}(a, t) P(a, t), \tag{20}\] with boundary condition \[ P(0, t) = \int^{\infty}_{0}d_{\rm P}(a, t) P(a, t) \, {\rm d}a. \tag{21}\]
The frequency-density of patches of age \(a < a^{\prime}\) is given by \(\int_{0}^{a^{\prime}}P(a, t) \, {\rm d}a\), with \(\int_{0}^{\infty} P(a, t) \, {\rm d}a = 1\). If \(\frac{\partial}{\partial t}d_{\rm P}(a, t) = 0\), then \(P(a)\) will approach an equilibrium solution given by \[ P(a) = P(0) \, S_{\rm P}(0, a), \tag{22}\] where \[ S_{\rm P}(0, a) = \exp \left( - \int_{0}^{a} d_{\rm P}(a^\prime) \, {\rm d}a^\prime\right) \tag{23}\] is the probability that a patch remains undisturbed for duration \(a\), and \[ P(0) = \frac1{ \int_{0}^{\infty}S_{\rm P}(0, a) \, {\rm d}a} \tag{24}\] is the frequency-density of patches of age \(0\). The rate of disturbance for patches of age \(a\) is given by \(\frac{\partial (1-S_{\rm P}(0, a))}{\partial a} = - \frac{\partial S_{\rm P}(0, a)}{\partial a}\), while the expected lifetime of patches is \(- \int_0^\infty a \frac{\partial}{\partial a} S_{\rm P}(0, a) \, {\rm d} a = \int_0^\infty S_{\rm P}(0, a) \, {\rm d} a = \frac1{P(0)}\) (first step made using integration by parts).
An equilibrium distribution of patch ages may be achieved under a variety of conditions, for example, if \(d_{\rm P}(a, t)\) depends on patch age \(a\) but not on time \(t\). The rate of disturbance may also depend on features of the vegetation in the patch, rather than on patch age directly, in which case an equilibrium distribution of patch ages can still arise, provided the vegetation is also assumed to be at equilibrium.
Exponential distribution
If the rate \(d_{\rm P}\) of patch disturbance is constant with respect to patch age, then the rates at which patches of age \(a\) are disturbed follow an exponential distribution, \(-\partial S_{\rm P}(0, a)/ \partial a = d_{\rm P} \, \exp(-d_{\rm P} a)\). The distribution of patch ages is then given by \[ S_{\rm P}(0, a) = \exp\left(-d_{\rm P} a\right), \, P(0) = d_{\rm P}. \tag{25}\]
Weibull distribution
If the rate of patch disturbance changes with patch age according to \(d_{\rm P}(a) = \lambda \psi a^{\psi-1}\), then the rates at which patches of age \(a\) are disturbed follow a Weibull distribution, \(-\partial S_{\rm P}(0, a)/ \partial a = \lambda \psi a^{\psi -1}e^{-\lambda a^\psi}\). \(\psi>1\) implies that the probability of disturbance increases with patch age, while \(\psi<1\) implies that it decreases with patch age. For \(\psi = 1\), we obtain the exponential distribution, a special case of the Weibull distribution. The Weibull distribution results in the following distribution of patch ages, \[ S_{\rm P}(0, a) = \exp(-\lambda a^\psi), \, P(0) = \frac{\psi \lambda^{\frac1{\psi}}}{\Gamma\left(\frac1{\psi}\right)}, \tag{26}\] where \(\Gamma(x) = \int_{0}^{\infty}e^{-t}t^{x-1} \, dt\) is the gamma function. We can also specify this distribution by the mean disturbance interval \(\bar{a} = \frac1{P(0)}\). From this, we can calculate the relevant value for \(\lambda = \left(\frac{\Gamma\left(\frac1{\psi}\right)}{\psi \bar{a}}\right)^{\psi}\).
The default in plant is to assume \(\psi=2\), such that \(d_{\rm P}\) increases as a linear function of patch age. The distribution of patch ages is then specified by a single parameter, \(\bar{a}\).
Derivation of PDE describing size-structured dynamics
To model the metacommunity, we use a PDE describing the dynamics for a thin slice \(\Delta H\). This PDE can be derived as follows, following De Roos et al. (1997). Note that in this sub-section the dependencies on the traits \(x\) are deliberately not made explicit. Assuming that all rates are constant within the interval \(\Delta H\), the total number of individuals within the interval spanned by \([H - 0.5\Delta H, H + 0.5\Delta H)\) is \(N(H, a)\Delta H\). The flux of individuals in and out of this size interval can be expressed as \[ \begin{array}{ll} &g(H - 0.5 \Delta H, a) \, N(H - 0.5 \Delta H, a) - g(H + 0.5 \Delta H, a) \, N(H + 0.5 \Delta H, a) \\ & - d (H, a) \, N(H, a)\Delta H\\ \end{array}. \tag{27}\] The first two terms describe the flux in and out of the considered size interval through growth, while the last term describes losses through mortality. The change in the density of individuals within this size interval during a time step $\(a is thus\)$ \[\begin{array}{ll} N(H, a + \Delta a)\Delta H - N(H, a)\Delta H = &g(H - 0.5 \Delta H, a) \, N(H - 0.5 \Delta H, a)\Delta a \\ & - g(H + 0.5 \Delta H, a) \, N(H + 0.5 \Delta H, a)\Delta a\\& - d (H, a) \, N(H, a)\Delta H\Delta a. \end{array}\] \[ {#eq-flux-2} By rearranging, we obtain \] \[\begin{array}{ll} \frac{N(H, a + \Delta a) - N(H, a)}{\Delta a} = & - d (H, a) \, N(H, a) \\ & - \frac{g(H + 0.5 \Delta H, a) \, N(H + 0.5 \Delta H, a) - g(H - 0.5 \Delta H, a) \, N(H - 0.5 \Delta H, a)}{\Delta H}. \end{array}\]\[ {#eq-flux-3} The left-hand side above corresponds to the derivative of $N$ as $\Delta a\to 0$. For thin slices, $\Delta H \to 0$, this yields \] N(H, a) = - d (H, a) , N(H, a) - (g(H, a) , N(H, a)). $$ {#eq-pde-app}
To complete the model, this PDE must be supplemented with boundary conditions that specify the density at the lower end of heights for all \(a\), as well as the initial distribution \(N(H,0)\). The former is derived by integrating the PDE with respect to \(H\) over the interval \((H_{0}, H_{\infty} )\), yielding \[ \frac{\partial}{\partial a} \int _{H_{0} }^{H_{\infty}}N(H, a) \, {\rm d} H = g(H_{0} , a) \, N(H_{0} , a) - g(H_{\infty}, a) \, N(H_{\infty}, a) - \int _{H_{0} }^{H_{\infty}}d (H, a) \, N(H, a) \, {\rm d} H. \tag{28}\] The left-hand size above is evidently the rate of change of the total density of individuals in the population, while the right-hand side is the population’s total death rate. Further, we can assume that \(N(H_{\infty}, a) = 0\). Thus, to balance total births and deaths, \(g(H_{0} , a) \, N(H_{0} , a)\) must equal the population’s total birth rate \(B\), yielding the boundary condition \[ g(H_{0} , a) \, N(H_{0} , a) = B. \tag{29}\]
Converting density from one size unit to another
The quantification of size-density depends on the chosen size unit (Equation 5). Thus, what if we want to express size-density with respect to another size unit? A relation between the two corresponding values of size-density can be derived by noting that the total density of individuals within a given size range must be equal. We consider a situation in which size-density is expressed in units of size \(M\), but we want it in units of size \(H\). First, we require a one-to-one function that translates the first size unit into the second, \(H = \hat{H}(M)\). Then the following must hold \[ \int_{M_1}^{M_2} N(M | x, a) \, {\rm d}M = \int_{\hat{H}(M_1)}^{\hat{H}(M_2)} N^\prime(H | x, a) \, {\rm d}H . \tag{30}\] For very small size intervals, this equation is equivalent to \[ \left(M_2 - M_1 \right) \, N(M_1 | x, a) = \left( \hat{H}(M_2) - \hat{H}(M_1)\right) \, N^\prime(\hat{H}(M_1) | x, a). \tag{31}\] Rearranging gives \[ N^\prime(\hat{H}(M_1) | x, a) = N(M_1 | x, a) \, \frac{M_2 - M_1}{\hat{H}(M_2) - \hat{H}(M_1)}. \tag{32}\] Noting that the second term on the right-hand side is simply the definition of \(\frac{{\rm d} M}{{\rm d} H}\) evaluated at \(M_1\), we have \[ N^\prime(H | x, a) = N(M | x, a) \, \frac{{\rm d} M}{{\rm d} H}. \tag{33}\]