PSPManalysis: Steady‐state and bifurcation analysis of physiologically structured population models

How environmental conditions affect the life history of individuals and how these effects shape population and community dynamics on ecological and evolutionary time‐scales is a central question in many eco‐evolutionary studies. Physiologically structured population models (PSPMs) allow to address this question theoretically as PSPMs are built on a function‐based life‐history model, which explicitly describes how life history depends on individual traits and environmental factors. PSPMs furthermore explicitly account for population feedback on these environmental factors, which translates into density‐dependent effects on life history. PSPManalysis is an r package that allows to simulate ecological dynamics of PSPMs, compute their ecological steady states as a function of model parameters and detect bifurcation points in the computed curves where dynamics change drastically. It furthermore allows for analysing evolutionary dynamics and evolutionary singular states of PSPMs based on Adaptive Dynamics theory. The package only requires a relatively straightforward specification of the life‐history functions as input. Compared to dynamic simulations alone, PSPManalysis uses methods from bifurcation analysis to gain a more complete and comprehensive understanding of model behaviour which is much less dependent on particular parameter values or initial model conditions. Given the central role of the individual life history in many studies, there is substantial scope for using PSPManalysis in fields as diverse as ecology, ecotoxicology, conservation biology and evolutionary biology.


| INTRODUC TI ON
The individual life history plays a central role in ecology and evolution, determining demography and persistence of populations and, together with interactions with other species, shaping the dynamics of interacting populations and communities. An individual's life history furthermore influences its fitness and thereby governs the evolutionary change in species traits. Methodologies to assess how life-history characteristics translate into consequences at the population level are hence a core part of many eco-evolutionary studies.

de ROOS
Physiologically structured population models (PSPMs, de Roos, 1997;Metz & Diekmann, 1986) provide a theoretical approach to analyse links between individual life history and population dynamics, as PSPMs describe population dynamics on the basis of a function-based model of the life history, which also accounts for effects of environmental variables, such as food availability and predator density (de Roos, 2020). PSPMs furthermore account for how these environmental variables are impacted by the population as a whole. While matrix (Caswell, 2001) and integral projection models (Ellner et al., 2016) are better suited to analyse life-history observations and infer their population dynamical consequences (de Roos, 2020), PSPMs capture with considerable mechanistic detail how individual-level processes, like energetics, together with interactions of the individual with its environment shape its life history and how feedback of the entire population on this environment has density-dependent impacts on that life history. PSPMs are therefore particularly useful to analyse how particular mechanisms or aspects of the life history or ecology of an individual would affect the population and community dynamics. The downside of PSPMs, however, is their mathematical tractability. The most simple PSPMs that only account for population size-structure can be formulated using partial differential equations (de Roos, 1997;Metz & Diekmann, 1986), but in general PSPMs are more appropriately couched in terms of coupled systems of nonlinear renewal equations and differential delay equations (Diekmann et al., 2007). Recently, numerical methodology (Kirkilionis et al., 2001;Sánchez-Sanz & Getto, 2016) was developed that allows for analysing ecological and evolutionary dynamics of even fairly complicated PSPMs (ten Brink et al., 2019;Chaparro-Pedraza & de Roos, 2020;Hin & de Roos, 2019) without bothering about these complicated population-level model formulations. This paper introduces the package PSPManalysis implementing this methodology. Because of their formulation in terms of partial differential equations PSPMs are often analysed using numerical simulations of the dynamics at particular parameter values for specific initial conditions. PSPManalysis implements functions for such simulations of ecological dynamics using the Escalator Boxcar Train (EBT) method (de Roos, 1997;de Roos et al., 1992). More importantly, however, PSPManalysis uses the theory on bifurcations in nonlinear dynamical systems (Kuznetsov, 1998) for model analysis. Bifurcation analysis provides more complete and comprehensive understanding of model dynamics because the results are much less dependent on particular parameter values or initial conditions (Kuznetsov, 1998).

| OVERVIE W OF SOF T WARE
For that purpose, PSPManalysis includes functions to compute demographic quantities, such as population growth and stable population states (de Roos, 2008), as well as equilibrium population states (Diekmann et al., 2003) of PSPMs. These models can be of arbitrary complexity with individuals characterized by multiple variables (traits) and multiple environmental variables, such as resource and predator densities, influencing their life history. PSPManalysis also uses curve continuation techniques (Sánchez-Sanz & Getto, 2016) to compute these equilibrium states over ranges of model parameters and detects the so-called bifurcation points, at which a qualitative change in model dynamics occurs (Kuznetsov, 1998). In addition, PSPManalysis allows for an evolutionary analysis of the computed ecological steady states based on the framework of 'Adaptive Dynamics' (Dieckmann & Law, 1996;Metz et al., 1996).
PSPManalysis thus allows for more extensive analysis of a substantially larger class of PSPMs than other packages geared at structured population models such as Mizer (Scott et al., 2014) and plant (Falster et al., 2015) for dynamic simulations of size-spectrum models and size-structured plant populations, respectively, or IPMpack (Metcalf et al., 2013) for demographic analysis of integral projection models.
The main computational routines of the PSPManalysis package are implemented in C for performance reasons with R functions (R Core Team, 2020) providing the interface to these routines. Finally, the package includes a very detailed manual, which discusses the full functionality of the package and illustrates its use with step-by-step instructions.

| ANALYS IS OF AN E X AMPLE MODEL
The life-history model described in Chaparro-Pedraza and de Roos (2020) is used here to illustrate the functionality of the PSPManalysis package. This model is based on the life history of salmon with individuals starting life in a nursery habitat, where they are protected from predation but compete for a shared resource X. Individuals subsequently shift to a growing habitat with negligible and ad libitum food but where they experience predation mortality. All model equations are presented in Table 1, while Table S1 lists default parameters values.
Individuals are characterized by their age a as well as their length l and the focal population (from here on called the 'consumer') is therefore age-and length-structured. Migration to the growing habitat and maturation occur on reaching length l = l s and l = l m respectively. Resource feeding, growth in body size and reproduction are linked through a dynamic energy budget (DEB) model (Chaparro-Pedraza & de Roos, 2020). The DEB model predicts individuals to grow following a von Bertalanffy growth curve with ultimate body length equal to l ∞ X/(K + X) and l ∞ in the nursery and growth habitat respectively (Table 1). Reproduction occurs after migration to the growth habitat since l s < l m , while fecundity is proportional to squared individual length.
Mortality in the nursery habitat is constant, but in the growth habitat decreases with body size (Table 1). In contrast to Chaparro-Pedraza and de Roos (2020), size-dependent mortality is assumed proportional to the density of predators that prey on consumers following a linear functional response and experience a mortality rate μ p . The scaled predator density P incorporates the conversion efficiency between ingested consumer biomass and the predator's numerical response B (Table 1). The contribution of consumer individuals to predator intake equals the product of their vulnerability to predation (l −d ) and their biomass (l 3 ). Resource turnover in the nursery habitat follows semi-chemostat growth dynamics.
The interaction between resource, structured consumer and predator is fully determined by seven functions (Table 1): three life-history functions, describing development, reproduction and mortality, two functions describing consumers' impact on their environment through resource foraging and contribution to predator intake and two functions determining resource and predator dynamics in the nursery and growth habitat respectively. The Supporting Information (Section 2) shows how to implement these functions in an R script for analysis with PSPManalysis (implementation in C is possible and in fact preferable as it speeds up computations by a factor 40-50).
The function PSPMequi, which is the main R function in the PSPManalysis package, uses the methodology of Kirkilionis et al. (2001) and Sánchez-Sanz and Getto (2016) (Supporting Information, Section 5) to compute steady states of PSPMs as a function of a model parameter, called the bifurcation parameter. Figure 1, which was produced from output of PSPMequi using basic plotting commands in R, shows the results of such computations as a function of maximum resource density X max for the example model in Table 1. Three computational steps generated the data for Figure 1.
The R commands for these computational steps are discussed below (see Supporting Information, Section 3, for a more detailed presentation). After installing and loading the PSPManalysis package (see Data availability statement below) running the command: demo("Salmon", package = "PSPManalysis", echo = FALSE) illustrates the generation of all figures in a step-by-step manner and thus supports the following presentation.
For low maximum resource densities, the example model only allows a steady state with zero consumer and predator density and the resource equal to its maximum density (X = X max ), because consumers do not encounter sufficient food to mature (maturation length equals l ∞ X/(K + X)). Starting in this resource-only equilibrium for X max = 0.1, the following call to PSPMequi was used to compute resource-only equilibria for increasing values of X max : EqR <-PSPMequi( modelname = "Salmon.R", biftype = "EQ", startpoint = c(0.1, 0.1), stepsize = 0.5, parbnds = c(1, 0.01, 10), options = c("popZE", "0", "envZE", "1")) This call resulted in the thin curve section with increasing densities X at low values of X max in Figure  The two elements "popZE" and "0" in the options argument of the first command instruct the function to assume a zero equilibrium density for the structured population with index 0 TA B L E 1 Life-history functions a of the model of Chaparro-Pedraza and de Roos (2020)

Function Description
Life-history functions Mortality Impacts on environment Resource forageing in nursery habitat Contribution to predator numerical response Functions related to the environment (scaled) predator density; B (day −1 ): predator numerical response; ξ (day −1 ): von Bertalanffy growth rate; l s and l m (cm): body size at habitat shift and maturation; l ∞ (cm): maximum body size at maximum feeding; K (g/m 3 ): half-saturation resource density; B max (cm −2 /day) and I max (g cm −2 day −1 ): fecundity and ingestion proportionality constants; μ 1 and μ 2 (day −1 ): mortality rate in habitats 1 and 2; d (−): exponent in sizedependent predation; ϕ (cm d m 3 /day) and μ p (day −1 ): predator attack and mortality rate; ρ (day −1 ) and X max (g/m 3 ): turnover rate and maximum density of resource; See Table S1 for full details.  (Table S1) Age (days) Body length (cm) (the consumer; because PSPManalysis is written in C, the first vector element has index 0 rather than 1) and the two elements "envZE" and "1" do the same for the environmental variable with index 1 (= 2nd vector element, the unstructured predator).
Indicating that these two populations have zero equilibrium density simplifies and speeds up computations. Because these two populations have zero density, the argument startpoint contains only two values: the value of the bifurcation parameter X max at which to start the computations and an estimate for the equilibrium resource density in this state. This computational step is only useful because PSPMequi detects a bifurcation point along the curve, labelling it with the string 'BP #0' (Figure 1) to indicate it represents a branching (or transcritical bifurcation) point (Kuznetsov, 1998) (Figure 1), which is returned in the output list element EqCR$bifpoints of the previous step, as starting point to compute steady states with positive predator density as a function of X max : EqPCR <-PSPMequi( modelname = "Salmon.R", biftype = "EQ", startpoint = EqCR$bifpoints[1,1:4], stepsize = -0.5, parbnds = c(1, 0.0, 10)) The computation of this curve starts off to lower values of X max (notice the negative stepsize argument) as otherwise negative predator densities would result. The result of this computation is a folded curve that extends to a minimum at X max ≈ 4. The function PSPMequi labels this minimum as 'LP' (Figure 1), thus classifying it as a limit (or saddle-node bifurcation) point (Kuznetsov, 1998).
Ecologically, this minimum X max value represents the persistence boundary of the predator, whereas the branching point labelled 'BPE #1' represents its invasion boundary. General bifurcation theory (Kuznetsov, 1998) stipulates that the curve section connecting these two bifurcation points represents unstable steady states or saddle points. PSPMequi also allows for computing the location of the predator's invasion and persistence boundary (labelled 'BPE #1' and 'LP', respectively, in Figure 1) dependent on a second model parameter. Figure 2 illustrates this for the maximum resource density X max and predator mortality rate μ p . For parameter combinations between the invasion ('BPE #1') and persistence boundary ('LP') in Figure 2, two potentially stable steady states occur, a tritrophic steady state with predators and a consumer-resource steady state that predators cannot invade.
The PSPManalysis package also allows for calculating evolutionary singular strategies using the framework of Adaptive Dynamics (Dieckmann & Law, 1996;Metz et al., 1996) and thus permits studying the evolution of life-history traits in a population and community context. While computing equilibrium curves F I G U R E 2 Location of the bifurcation points labelled 'BPE #1' and 'LP' in Figure 1 dependent on maximum resource density X max and predator mortality rate μ p . Otherwise default parameter values (Table S1). See Supporting Information (Section 3) for a detailed discussion of the R commands to compute the data for this figure Maximum resource density (g/m 3 ) Predator mortality (day -1 ) de ROOS as function of a life-history parameter the function PSPMequi can produce as output the selection gradient on this parameter, which equals the derivative of the lifetime reproductive output R 0 with respect to the life-history parameter (Diekmann et al., 2003). Figure 3 shows the equilibrium curves of the example model dependent on the length at habitat shift l s , including its selection gradient (see Supporting Information, Section 3, for details about the R commands for these computations). In the absence of predators, smaller sizes at habitat shift are selected, but with predators present PSPMequi detects an evolutionarily singular state (Brännström et al., 2013), which it labels as 'CSS #0' on the basis of second-order derivatives of R 0 with respect to the length at habitat shift l s (Geritz et al., 1998). The label indicates that this singular state is convergent stable, such that the value of l s will evolve towards this CSS value, while after fixation mutants with slightly different values of l s cannot invade. Starting from this CSS, the function PSPMequi can also compute the curve with zero mutant fitness dependent on both resident and mutant life-history trait value, which corresponds to the boundary separating regions with positive and negative mutant fitness in the pairwise invasibility plot (or PIP, van Tienderen & de Jong, 1986) shown in Figure 4 (left).
Finally, the evolutionary dynamics of life-history parameters, as described by the canonical equation of Adaptive Dynamics (Dieckmann & Law, 1996), can be simulated with the function PSPMevodyn in the PSPManalysis package (Figure 4, right). The trajectory of the length at habitat shift l s over evolutionary time shown F I G U R E 3 Equilibrium densities of predator (top) and consumers (middle) in the nursery (blue) and growth habitat (red) dependent on the length at habitat shift, plus the selection gradient (bottom) on this parameter. Solid and dashed lines represent possibly stable equilibria and saddle points respectively. Curve sections representing consumer-resource steady states that can be invaded by predators are omitted for clarity. See text for details about the bifurcation points labelled 'BPE #1', 'LP' and 'CSS #0'. Otherwise default parameter values (Table S1). See Supporting Information (Section 3) for a detailed discussion of the R commands to compute the data for this figure Methods in Ecology and Evoluঞon de ROOS in Figure 4 (right) confirms that length at habitat shift evolves to its CSS value, but since the function PSPMevodyn cannot simulate combined mutant and resident dynamics, it cannot verify whether evolutionary branching occurs at this evolutionary singular state.

| D ISCUSS I ON
The methodology for analysing PSPMs provided by PSPManalysis has previously been used to investigate how ontogeny affects ecological dynamics of size-structured communities (de Roos & Persson, 2013) and, more recently, to study evolution of metamorphosis (ten Brink et al., 2019), cannibalism (Hin & de Roos, 2019) and timing of habitat shifts (Chaparro-Pedraza & de Roos, 2020).
Given the importance of individual life history for eco-evolutionary dynamics and of environmental feedback on this life history, the methodology is, however, applicable to a wide range of eco-evolutionary questions. This includes questions in ecotoxicology and conservation biology, as PSPManalysis is especially suited to investigate links between the dynamic energy budget (DEB) of individuals and its population consequences (de Roos & Persson, 2013) and DEB models are widely used to assess the consequences of toxicants and changing temperature (Nisbet et al., 2000).
The function PSPMequi is the main component of PSPManalysis and its functionality is more extensive than discussed here. Next The author thanks two anonymous reviewers for their helpful comments that greatly improved the paper.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns. com/publo n/10.1111/2041-210X.13527.

DATA AVA I L A B I L I T Y S TAT E M E N T
The PSPManalysis package is available on CRAN (https://CRAN.Rproje ct.org/packa ge=PSPMa nalysis) and can be installed and subsequently loaded using the commands: install.packages("PSPManalysis") library(PSPManalysis) After loading the package, all computations necessary to reproduce the figures presented in this paper can be executed step-by-step by executing the following demo() command: demo("Salmon", package = "PSPManalysis", echo = FALSE) No data have been used in this paper, other than the contents of the PSPManalysis package available on CRAN.