This document outlines methods used to model demography in the plant package. We first outline the system of equations being solved. 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. We then describe the numerical techniques used to solve the equations. Table 1 provides a list of names and definitions used throughout this document.
| 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 seed rain 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 |
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.
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 Eqs. (1)-(3) potentially obscures an important point: these equations are simply the general, non-linear solutions to integrating growth, mortality, and fecundity over time.
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 Eqs. (1)-(3). 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 ofplant, 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\),
\[\begin{equation}
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}
\end{equation}\]
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)), \[\begin{equation} \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} \end{equation}\]See section derivation-of-pde-describing-size-structured-dynamics for the derivation of this PDE.
Eq. (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, \[\begin{equation} 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} \end{equation}\]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 Eq. (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, \[\begin{equation} N\left(H|x,0\right) = 0, \tag{7} \end{equation}\]although non–zero distributions could alternatively be specified (e.g. Kohyama, 1993, Moorcroft et al. (2001)).
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 derivation-of-pde-describing-age-structured-dynamics for details).
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 Eqs. (5)-(12).
The seed rain of each species in the metacommunity is given by the rate at which seeds are produced across all patches, \[\begin{equation} 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} \end{equation}\]where \(S_{\rm D}\) is the survival probability of seeds during dispersal.
A convenient feature of Eqs. (5) - (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 seed rains 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 Eqs. (6)-(13) simultaneously, for all species with traits \(x\) in the metacommunity.
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, \[\begin{equation} 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} \end{equation}\] and across the entire metacommunity, \[\begin{equation} \hat{W} = \int_{0}^{\infty} P(a) \, W(a) \, {\rm d}a. \tag{19} \end{equation}\]We now consider how to estimate the fitness of a rare mutant individual with traits \(x^\prime\) growing in the light environment of a resident community with traits \(x\). We focus on phenotype-dependent components of fitness – describing the aggregate consequences of a given set of traits for growth, fecundity, and mortality – taking into account the non-linear effects of competition, but ignoring the underlying genetic basis for trait inheritance and expression. We also adhere to standard conventions in such analyses by assuming that the mutant phenotype is sufficiently rare to have a negligible effect on the light environment where it is growing (Geritz, Kisdi, Meszéna, & Metz, 1998). The approach for calculating fitness implemented here follows the approach described by Falster et al. (2015).
In general, invasion fitness is defined as the long-term per capita growth rate of a rare mutant phenotype in the environment determined by a resident phenotype (Metz, Nisbet, & Geritz, 1992). Calculating per capita growth rates, however, is particularly challenging in a structured metacommunity model (Gyllenberg & Metz, 2001, Metz & Gyllenberg (2001)). As an alternative measure of invasion fitness in metacommunities, we can use the basic reproduction ratio measuring the expected number of new dispersers arising from a single dispersal event. For metacommunities at demographic equilibrium, evolutionary inferences made using basic reproduction ratios are equivalent to those made using per capita growth rates (Gyllenberg & Metz, 2001, Metz & Gyllenberg (2001)).
We denote by \(R\left(x^\prime, x\right)\) the basic reproduction ratio of individuals with traits \(x^\prime\) growing in the competitive light environment of the resident traits \(x\). Recalling that patches of age \(a\) have frequency-density \(P(a)\) in the landscape, it follows that any seed with traits \(x^\prime\) has a probability of \(P(a)\) of landing in a patch of age \(a\). The basic reproduction ratio for individuals with traits \(x^\prime\) is then \[\begin{equation} R\left(x^\prime,x\right) = \int _0^{\infty} P\left(a\right) \, \tilde{R}\left(x^\prime, a, \infty \right) \, {\rm d}a , \tag{20} \end{equation}\] where \(\tilde{R}\left(x^\prime, a_0, a \right)\) is the expected number of dispersing offspring produced by a single dispersing seed with traits \(x^\prime\) arriving in a patch of age \(a_0\) up until age \(a\) (Gyllenberg & Metz, 2001, Metz & Gyllenberg (2001)). In turn, \(\tilde{R}\left(x^\prime, a,\infty\right)\) is calculated by integrating an individual’s fecundity over the expected lifetime of a patch, taking into account competitive shading from residents with traits \(x\), the individual’s probability of surviving, and its traits, \[\begin{equation} \tilde{R}(x^\prime, a_0, \infty) = \int_{a_0}^{\infty} S_{\rm D} \, f(x^\prime, H(x^\prime, a_0, a), E_{a}) \, S_{\rm I} (x^\prime, a_0, a) \, S_{\rm P} (a_0, a) \, {\rm d} a. \tag{21} \end{equation}\]The equations below describe how we solve the trait-, size-, and patch-structured population dynamics specified in the previous section. Our approach to solving for the size-density distribution is based on the characteristic method (Angulo & López-Marcos, 2004, Angulo, López-Marcos, & López-Marcos (2014), Angulo, López-Marcos, & López-Marcos (2016)). The characteristic method is related to the Escalator Boxcar Train technique (De Roos, 1988, De Roos et al. (1997), De Roos, Diekmann, & Metz (1992)), but not identical to it.
When simulating an individual plant, or the development of a patch, we need to solve for the size, survival, and seed output of individual plants. When solving for the size-density distribution in a large patch, we also need to estimate the average abundance of individuals. Each of these problems is formulated as an initial-value ODE problem (IVP), which can be solved using an ODE stepper.
where \(S_{\rm I}\) is calculated as described above and \(S_{\rm P}\) is calculated as in Eq. (10)$.
Eq. (29) states that the size-density \(N\) at a specific patch age \(a\) is the product of the size-density at patch age \(a_{0}\) adjusted for changes through growth and mortality. Size-density decreases through time because of mortality, as in a typical survival equation, but also changes because of growth. If growth is slowing with size, (i.e., \(\partial g / \partial H < 0\)), size-density will increase since the characteristics compress. Conversely, size-density will increases if \(\partial g / \partial H > 0\).
Denoting by \(\left[H_0, H_{ + } \right)\) the range of heights attainable by any individual, our algorithm for solving metacommunity dynamics proceeds by sub-dividing this interval into a series of cohorts with heights \(H_0 < H_1 < \ldots < H_k\) at the initial points of the characteristic curves. These cohorts are then transported along the characteristics of Eq. (5). The placement of cohorts is controlled indirectly, via the schedule of patch ages at which new cohorts are introduced into the metacommunity. We then track the demography of each such cohort.
The integral in Eq. (29) is solved via the IVP \[\begin{equation} \frac{dy}{dt} = \frac{\partial g(x, H_i(t), E_t)}{\partial H} + d(x, H_i(t), E_t), \tag{30} \end{equation}\] \[\begin{equation} y(0) = - \ln\left(N(H_0 | x, a_0) \, S_{\rm G} (x, H_0, E_{{\rm a}0}) \right), \tag{31} \end{equation}\] from which we obtain the size-density \[\begin{equation} N(H_0 | x, a_0) = \exp( - y(a)). \tag{32} \end{equation}\]We now outline how to control the error of the approximate solutions to the system of equations described above. In our algorithm, numerical solutions are required to address a variety of problems:
As with all numerical techniques, solutions to each of these problems are accurate only up to a specified level. These levels are controlled via parameters in the plant code. Below, we provide a brief overview of the different numerical techniques being applied and outline how error tolerance can be increased or decreased. We refer to various control parameters that can be found within the control object. For a worked example illustrating how to modify these control parameters, see the section parameters of Appendix S3.
When a seed germinates, it produces a seedling of given height. The height of these seedlings is assumed to vary with the seed mass. Because there is no analytical solution relating seedling height to seed mass – at least when using the default FF16 physiological model – we must solve for this height numerically. The calculation is performed by the function height_seed within the physiological model, using the Boost library’s one-dimensional bisect routine (Schäling, 2014, Eddelbuettel, Emerson, & Kane (2015)). The accuracy of the solution is controlled by the parameter plant_seed_tol.
All of the IVPs outlined above (in section approach) must be stepped through time. For this, plant uses the embedded Runge-Kutta Cash-Karp 4-5 algorithm (Cash & Karp, 1990), with code ported directly from the GNU Scientific Library (Galassi, 2009). The accuracy of the solver is controlled by two control parameters for relative and absolute accuracy, ode_tol_rel and ode_tol_abs.
Errors in the approximation of the size-density distribution arise from two sources: (i) coarse stepping of cohorts through time and (ii) poor spacing of cohorts across the attainable size range.
As described above, the stepping of the ODE solver is controlled by two control parameters for relative and absolute accuracy, ode_tol_rel and ode_tol_abs.
A second factor controlling the accuracy with which cohorts are stepped through time is the accuracy of the derivative calculation according to Eq. (29), calculated via standard finite differencing (Abramowitz & Stegun, 2012). When the parameter cohort_gradient_richardson is TRUE a Richardson extrapolation (Stoer & Bulirsch, 2002) is used to refine the estimate, up to depth cohort_gradient_richardson. The overall accuracy of the derivative is controlled by cohort_gradient_eps.
The primary factor controlling the spacing of cohorts is the schedule of cohort introduction times. Because the system of equations to be integrated is deterministic, the schedule of cohort introduction times determines the spacing of cohorts throughout the entire development of a patch. Poor cohort spacing introduces error because various emergent properties – such as total leaf area, biomass, or seed output – are estimated by integrating over the size-density distribution. The accuracy of these integrations declines directly with the inappropriate spacing of cohorts. Thus, our algorithm aims to build an appropriately refined schedule, which allows the required integrations to be performed with the desired accuracy at every time point. At the same time, for reasons of computational efficiency, we want as few cohorts as possible.
A suitable schedule is found using the function build_schedule. This function takes an initial vector of introduction times and considers for each cohort whether removing that cohort causes the error introduced when integrating two specified functions over the size-density distribution to jump over the allowable error threshold schedule_eps. This calculation is repeated for every time step in the development of the patch. A new cohort is introduced immediately prior to any cohort failing these tests. The dynamics of the patch are then simulated again and the process is repeated, until all integrations at all time points have an error below the tolerable limit schedule_eps. Decreasing schedule_eps demands higher accuracy from the solver, and thus increases the number of cohorts being introduced. Note that we are assessing whether removing an existing cohort causes the error to jump above the threshold limit, and using this to decide whether an extra cohort – in addition to the one used in the test – should be introduced. Thus, the actual error is likely to be lower than, but at most equal to, schedule_eps. The general idea of adaptively refining the vector of introduction times was first applied by Falster et al. (2011), and was described further by Falster et al. (2015).
To determine the error associated with a given cohort, we integrate two different functions over the size-density distribution, within the function run_scm_error. We then assess how much removing the focal cohort increases the error in these two integrations. The first integration, performed by the function area_leaf_error, determines how much the removal of the focal cohort increases the error in the estimate of the total leaf area in the patch. The second integration, performed by the function seed_rain_error, determines how much the removal of the focal cohort increases the error in the estimate of the total seed production from the patch. The relative error in each integration is then calculated using the local_error_integration function.
For a worked example illustrating the build_schedule function, see the section cohort_spacing of Appendix S3.
To progress with solving the system of ODEs requires that we calculate the amount of shading on each of the cohorts, from all other plants in the patch.
Calculating the canopy openness \(E_a(z)\) at a given height \(z\) in a patch of age \(a\) requires that we integrate over the size-density distribution (Eq. (4). This integration is performed using the trapezium rule, within the function area_leaf_above in species.h. The main factor controlling the accuracy of the integration is the spacing of cohorts. The cohort introduction times determining the spacing of cohorts are adaptively refined as described above. This implies that also the trapezium integration within the area_leaf_above function is adaptively refined via the build_schedule function.
The cost of calculating \(E_a(z)\) linearly increases with the number of cohorts in the metacommunity. Since the same calculation must be repeated for every cohort, the overall computational cost of a step increases as \(O(k^2 )\), where \(k\) is the total number of cohorts across all species. This disproportionate increase in computational cost with the number of cohorts is highly undesirable.
We reduce the computational cost from \(O(k^2)\) to \(O(k)\) by approximating \(E_a(z)\) with a spline. Eq. (4) describes a function monotonically increasing with size. This function is easily approximated using a piecewise continuous spline fitted to a limited number of points. Once fitted, the spline can be used to estimate any additional evaluations of competitive effect. Since spline evaluations are computationally cheaper than integrating over the size-density distribution, this approach reduces the overall cost of stepping the resident population. A new spline is constructed for each time step.
The accuracy of the spline interpolation depends on the number of points used in its construction and on their placement along the size axis. We select the number and locations of points via an adaptive algorithm. Starting with an initial set of 33 points, we assess how much each point contributes to the accuracy of the spline fit at the location of each cohort, first via exact calculation, and second by linearly interpolating from adjacent cohorts. The absolute difference in these values is compared to the control parameter environment_light_tol. If the error is greater than this tolerance, the interval is bisected and the is process repeated (see adaptive_interpolator.h for details).
Plants have leaf area distributed over a range of heights. Estimating a plant’s assimilation at each time step thus requires integrating leaf-level rates over the plant. The integration is performed using Gaussian quadrature, using the QUADPACK routines (Piessens, Doncker-Kapenga, Überhuber, & Kahaner, 1983) adapted from the GNU Scientific Library (Galassi, 2009); see qag.h for details. If the control parameter plant_assimilation_adaptive is TRUE, the integration is performed using adaptive refinement with an accuracy controlled by the parameter plant_assimilation_tol.
For a single species, solving for \(Y_x\) is a straightforward one-dimensional root-finding problem, which can be solved with accuracy equilibrium_eps via a simple bisection algorithm (see equlibrium.R for details).
Solving for seed rains in metacommunities with multiple species is significantly harder, because there is no generally applicable method for multi-dimensional root finding. In plant, we have therefore implemented several different approaches (see equlibrium.R for details).
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.
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}\).
Abramowitz, M., & Stegun, I. A. (2012). Handbook of mathematical functions: With formulas, graphs, and mathematical tables. Courier Corporation.
Angulo, & López-Marcos. (2004). Numerical integration of fully nonlinear size-structured population models. Applied Numerical Mathematics, 50(3-4), 291–327. doi:10.1016/j.apnum.2004.01.007
Angulo, O., López-Marcos, J. C., & López-Marcos, M. A. (2014). Analysis of an efficient integrator for a size-structured population model with a dynamical resource. Computers & Mathematics with Applications, 68(9), 941–961. doi:10.1016/j.camwa.2014.04.009
Angulo, O., López-Marcos, J. C., & López-Marcos, M. A. (2016). Study on the efficiency in the numerical integration of size-structured population models: Error and computational cost. Journal of Computational and Applied Mathematics, 291, 391–401. doi:10.1016/j.cam.2015.03.022
Cash, J. R., & Karp, A. H. (1990). A variable order Runge-Kutta method for initial value problems with rapidly varying right-hand sides. ACM Transactions on Mathematical Software, 16(3), 201–222. doi:10.1145/79505.79507
De Roos, A. M. (1988). Numerical methods for structured population models: The Escalator Boxcar Train. Numerical Methods for Partial Differential Equations, 4(3), 173–195. doi:10.1002/num.1690040303
De Roos, A. M., Diekmann, O., & Metz, J. A. J. (1992). Studying the dynamics of structured population models: A versatile technique and its application to daphnia. American Naturalist, 139(1), 123. doi:abs/10.1086/285316
De Roos, A. M., Taljapurkar, S., & Caswell, H. (1997). A gentle introduction to physiologically structured population models. In Structured population models in marine, terrestrial and fresh-water systems (pp. 119–204). New York: Chapman & Hall.
Eddelbuettel, D., Emerson, J. W., & Kane, M. J. (2015). BH: Boost c++ header files. Retrieved from http://CRAN.R-project.org/package=BH
Falster, D. S., Brännström, Å., Dieckmann, U., & Westoby, M. (2011). Influence of four major plant traits on average height, leaf-area cover, net primary productivity, and biomass density in single-species forests: A theoretical investigation. Journal of Ecology, 99(1), 148–164. doi:10.1111/j.1365-2745.2010.01735.x
Falster, D. S., Brännström, Westoby, M., & Dieckmann, U. (2015). Multi-trait eco-evolutionary dynamics explain niche diversity and evolved neutrality in forests. bioRxiv, 014605. doi:10.1101/014605
Foerster, H. V. (1959). Some remarks on changing populations. In F. Stohlman (Ed.), The kinetics of cellular proliferation (pp. 382–407). Grune et Stratton, New York.
Galassi, M. (Ed.). (2009). GNU scientific library: Reference manual (3. ed., for GSL version 1.12). s.l.: Network Theory.
Geritz, S. A. H., Kisdi, Meszéna, G., & Metz, J. A. J. (1998). Evolutionarily singular strategies and the adaptive growth and branching of the evolutionary tree. Evolutionary Ecology, 12(1), 35–57. doi:10.1023/A:1006554906681
Gyllenberg, M., & Metz, J. A. J. (2001). On fitness in structured metapopulations. Journal of Mathematical Biology, 43(6), 545–560.
Kohyama, T. (1993). Size-structured tree populations in gap-dynamic forest: The forest architecture hypothesis for the stable coexistence of species. Journal of Ecology, 81(1), 131–143.
McKendrick, A. G. (1926). Applications of mathematics to medical problems. Proceedings of the Edinburgh Mathematical Society, 44, 1–34.
Metz, J. A. J., & Gyllenberg, M. (2001). How should we define fitness in structured metapopulation models? Including an application to the calculation of evolutionarily stable dispersal strategies. Proceedings of the Royal Society of London - Series B: Biological Sciences, 268(1466), 499–508.
Metz, J. A. J., Nisbet, R. M., & Geritz, S. A. H. (1992). How should we define ’fitness’ for general ecological scenarios? Trends in Ecology and Evolution, 7(6), 198–202. doi:10.1016/0169-5347(92)90073-K
Moorcroft, P. R., Hurtt, G. C., & Pacala, S. W. (2001). A method for scaling vegetation dynamics: The Ecosystem Demography model (ED). Ecological Monographs, 71(4), 557–586. doi:10.1890/0012-9615(2001)071[0557:AMFSVD]2.0.CO;2
Piessens, R., Doncker-Kapenga, E. D., Überhuber, C., & Kahaner, D. (1983). QUADPACK, a subroutine package for automatic integration. Berlin Heidelberg: Springer.
Schäling, B. (2014). The Boost C++ libraries (2nd ed.). XML Press.
Stoer, J., & Bulirsch, R. (2002). Introduction to numerical analysis. Springer Science & Business Media.