Predicting evolution over multiple generations in deteriorating environments using evolutionarily explicit Integral Projection Models

Abstract Human impacts on the natural world often generate environmental trends that can have detrimental effects on distributions of phenotypic traits. We do not have a good understanding of how deteriorating environments might impact evolutionary trajectories across multiple generations, even though effects of environmental trends are often significant in the statistical quantitative genetic analyses of phenotypic trait data that are used to estimate additive genetic (co)variances. These environmental trends capture reaction norms, where the same (average) genotype expresses different phenotypic trait values in different environments. Not incorporated into the predictive models typically parameterised from statistical analyses to predict evolution, such as the breeder's equation. We describe how these environmental effects can be incorporated into multi‐generational, evolutionarily explicit, structured population models before exploring how these effects can influence evolutionary dynamics. The paper is primarily a description of the modelling approach, but we also show how incorporation into models of the types of environmental trends that human activity has generated can have considerable impacts on the evolutionary dynamics that are predicted.

or abiotic environmental driver that alters survival and recruitment consequently has the potential to alter selection pressures and the rate, and potentially direction, of evolution. Evolution of these fitness-related traits in response to human-induced environmental change is a type of biotic change and can in turn influence survival and recruitment rates, and consequently the population dynamics, generating eco-evolutionary feedbacks (Hendry, 2016).
Human-induced biotic and abiotic environmental change can also impact phenotypic traits via phenotypic plasticity and nongenetic inheritance (Reed et al., 2011;Salinas et al., 2013;Via & Lande, 1985). These processes alter the map between genotype and phenotypic trait value such that the same genotype may generate different phenotypic trait values in different environments. The difference between the two is that phenotypic plasticity leads to environment-induced phenotypic changes in the individual experiencing environmental change, while nongenetic inheritance causes phenotypic changes in its offspring (Pigliucci, 2001). If the phenotypic trait an individual expresses is assumed to consist of a breeding value, determined by its genotype potentially at very many loci, and an environmental component (Falconer, 1960), phenotypic plasticity and nongenetic inheritance occur when environmental change alters the value of the environmental component of the phenotype (Via & Lande, 1985). Such dynamics are captured by reaction norms that describe how environmental variation influences phenotypic trait expression within a genotype (Falconer, 1990;Lande, 2009). If phenotypic plasticity or nongenetic inheritance change the distribution of phenotypic traits, this can alter survival and recruitment rates -and consequently population dynamics and selection. Such dynamics can occur even if the fitness functions themselves are not altered by environmental variation . Phenotypic plasticity and genetic inheritance consequently have the potential to generate eco-evolutionary feedbacks in the presence of unchanging fitness functions, with human-induced environmental change having considerable potential to be a major driver of such feedbacks. The question we ask here is how does the impact of human-induced environmental change on the environmental component of the phenotype influence eco-evolutionary dynamics, and the way populations respond to environmental change? Our results extend to any type of environmental change, but we couch this paper in terms of humaninduced change, and in particular in the impacts of a deteriorating environment.
Environmental variation can have substantial effects on phenotypic trait values as is widely appreciated in statistical quantitative genetics -a powerful framework for studying evolution (Falconer, 1960;Lynch & Walsh, 1998). In quantitative genetic analyses, it is often essential to fit variables into statistical models to correct for environmental influences on phenotypic trait expression (Kruuk, 2004;Merilä et al., 2001). For example, variables such as population density or weather attributes -that often show temporal trends as a result of human activity -are sometimes fitted as fixed effects into animal models of free-living populations (Fletcher et al., 2015;Kruuk et al., 2002;Potter et al., 2021), or year is fitted as a random effect (Kruuk, 2004;Lynch & Walsh, 1998). These environmental variables statistically adjust for reaction norms, allowing more robust estimates of additive genetic (co)variances by comparing phenotypic trait values amongst individuals of known relatedness once the effect of environmental variation on phenotypic trait values has been accounted for.
In predictive models widely applied to empirical systems, such as the breeder's equation, the additive genetic (co)variances are used to make evolutionary predictions but the evolutionary effects of environmental variation on phenotypic trait distributions are usually not incorporated (Chevin et al., 2010). We know that models such as the breeder's equation can provide accurate estimates of evolution over a single generation, but given the potential effects of environmental variation on selection via phenotypic plasticity and nongenetic inheritance, and via impacts on the fitness functions directly, these approaches may fail for the longer-term predictions required to understand how anthropogenic environmental change will impact populations (Morrissey et al., 2010). This leads us to pose the following hypothesis: to make multi-generational predictions of evolutionary change for populations in human-induced deteriorating environments it is necessary to model the effects of the environment on the dynamics of both the breeding value (via selection) and environmental component of the phenotype (via selection, phenotypic plasticity and nongenetic inheritance). We test this prediction by constructing simple evolutionarily explicit Integral Projection Models (hereafter called EE-IPMs) Coulson et al., 2017;Rees & Ellner, 2019).
Integral projection models (IPMs) are discrete-time population models structured by one or more continuous traits (Ellner et al., 2016). In addition, they can be structured by discrete characteristics such as age or sex (Ellner & Rees, 2006;Schindler et al., 2015).
The models are constructed from mathematical functions typically identified from statistical analyses. These functions describe (i) associations between the values of one or several phenotypic traits measured at time t and per-time step fitness (the fitness function) and (ii) phenotypic transitions between time t and t + 1 (transition functions) (Ellner et al., 2016). Most applications of IPMs assume time steps that are shorter than the generation length of the species being modelled. The fitness functions are typically divided into (i) the expected survival from t to t + 1 (the survival function), and (ii) the expected number of offspring produced between t and t + 1 that survive to recruit to the population at t + 1 (the recruitment function), while the transition functions are split into (iii) the values of trait(s) measured between t and t + 1 amongst surviving individuals (the development function) and (iv) the values of trait(s) measured in offspring when they recruit to the population at time t + 1.
This last function has been referred to as the inheritance function by some authors (Coulson et al., 2010), and this has caused confusion (Chevin, 2015). We refer to it here as the parent-offspring phenotypic similarity function. IPMs can also be constructed on a per-generation time step, where the recruitment function describes the association between a phenotypic trait and lifetime reproductive success, and the parent-offspring phenotypic difference function describes phenotypic trait similarity between parents and their offspring (Coulson et al., 2018). We use per-generation time step models here. Regardless of the approach, functions can be statistically estimated from individual-based phenotypic trait and demographic data that are used in statistical quantitative genetics and can include fixed and random effects describing how elements of the biotic or abiotic environment affect associations between trait values and each response variable (Coulson, 2012;Ellner et al., 2016).
The statistical functions are then combined to produce a projection model that iterates forward the distribution of phenotypic trait values from time t to time t + 1 (Easterling et al., 2000). At each time t, the projection model is usually approximated as a Lefkovich stagestructured matrix. The random and fixed effects identified in the statistical analyses of each function can be included in the projection model if the modeller desires so or else can be ignored (Coulson, 2012;Ellner et al., 2016). If the model includes elements of the biotic and abiotic environment, including human-induced environmental trends, the values in each matrix at each time t may vary between successive time steps. This generates a series of time-varying matrices that can often be analysed using approaches from random matrix theory (Tuljapurkar, 2013;Tuljapurkar & Caswell, 2012).
In evolutionarily explicit IPMs the phenotypic trait distribution is described as a multivariate distribution of components of the phenotypic trait(s) involved -for example each trait is decomposed into a bivariate distribution of the breeding values and the environmental components Coulson et al., 2017;Rees & Ellner, 2019). Selection operates on the phenotypic trait(s) under study, and this selection is then transmitted to each component of the phenotype. In models where the time step is shorter than the generation length, the breeding values remain fixed within individuals as they age, while, if desired, the environmental component of the phenotype may vary with the environment, generating phenotypic plasticity. The breeding values are genetically inherited , and assumptions about the effects of selection and inheritance on the additive genetic variance need to be explicitly specified . If desired, the environmental component of the phenotypic trait(s) in offspring can be either random, a function of that of their parents (nongenetic inheritance), or dependent upon the abiotic or biotic environment experienced by the offspring. So far, the first option has been most commonly used Coulson et al., 2017;Simmonds et al., 2020). Routinely, ran- Evolutionarily explicit Integral Projection Models are quite complex to construct and analyse , and there is a gap in the literature describing how simple versions of these models can be quite easily constructed. We attempt to fill this gap here. In doing this, we provide some novel biological insight by demonstrating how evolution will be fastest when it is cryptic, and slowest when phenotypic plasticity and nongenetic inheritance are adaptive.
Statistical quantitative genetics and structured population modelling both use similar data, and both have achieved considerable success in shining light on the complex patterns of phenotypic trait evolution, life-history evolution and demographic changes observed in the wild. EE-IPMs have been parameterized with estimates obtained through application of the animal model Simmonds et al., 2020). Despite this, the two approaches have largely independent histories of development, different lexicons, cite different literatures, and are used by different communities. As a consequence, crosstalk between advocates of the two approaches is not as frequent or constructive as it could be. We do not claim that our approach is the only way to link structured population modelling and quantitative genetics nor that our integration is complete.
We also do not generate new theory. Our aims, instead, are (i) to illustrate connections between the two approaches and thereby to, hopefully, encourage crosstalk, and (ii) explore how a humaninduced deteriorating environments might be expected to impact evolutionary trajectories.

| Modelling approach
In general, EE-IPMs assume that Coulson et al., 2017;Rees & Ellner, 2019): 1. An individual i's phenotypic trait value z i is the sum of a breeding value A i and an environmental component E i . The bivariate distribution of the components of a hypothetical phenotypic trait is given in Figure 1a.
2. The environmental component of the phenotype can be determined by random developmental noise and aspects of the external abiotic or biotic environment θ. Note that such effects are frequently corrected for in quantitative genetic statistical analyses, but are rarely incorporated into predictive models (Chevin et al., 2010). A temporally deteriorating environment, caused, for Note that such effects are sometimes corrected for in quantitative genetic statistical analyses as parental environmental effects (Lynch & Walsh, 1998) but are rarely incorporated into predictive models. In the models we construct here, we include nongenetic inheritance yet are silent on its underlying mechanistic causes.
We also make more assumptions specific to the models we re- The integral limits are taken to be below and above all possible values of A and E but are not displayed to simplify notation (Ellner et al., 2016). In addition, from now on we simply use a single, rather than a double, integral sign, with the reader determining the variables over which the integral is taken by the infinitesimals on the far-right hand side (dA and dE).

The kernel
, t is constructed to include the biological processes that determine the dynamics of the bivariate distribution of A and E and their drivers. The 'trick' in formulating a model is to specify the functions and the rules they encode. In our simple model, we will consider a fitness func- selection, followed by a parent-offspring phenotypic difference noting the implications of the notation change described above. We next turn to the second function in the kernel  (Easterling et al., 2000).
We now need to specify this parent-offspring phenotypic dif-  (Falconer, 1990). Next, because we want our model to be dynamic and to make predictions over multiple generations, we need to make assumptions about the dynamics of the variance of breeding values (the additive genetic variance), and how that changes (or not) from one generation to the next (Lande, 2009;Turelli & Barton, 1994). There are four ways in which the dynamics of the variance have been treated in structured models (Table 1), and the choice will depend upon the assumptions the researcher wishes to make.
Arguably the most intuitive way to generate the offspring distribution from the quantitative genetic perspective is to work with an algorithm as follows (approach 1 in Table 1 An alternative approach is to simply assume that the distribution of breeding values amongst offspring is always Gaussian and has a distribution with a mean breeding value A s , and a constant variance that does not change with time (Lande, 1982(Lande, , 2009). The two approaches differ and produce slightly different dynamics, because the approach based on convolutions does not necessarily produce a Gaussian distribution of offspring breeding values.
Next, we turn to rules for the environmental component of the phenotype. There may be three aspects we wish to incorporate into the postselection dynamics of E. First, random developmental noise ( Figure 2b) with a mean of zero and a fixed variance; second, the effects of the abiotic or biotic value in year t on the mean of the environmental component of the phenotype (Figure 2c) -that is, the processes that generate reaction norms (Chevin et al., 2010;Reed et al., 2011); and third, nongenetic inheritance -that is, a correlation between parental and offspring environmental components of the (1) phenotype caused by nongenetic inheritance (Figure 2d) (Salinas et al., 2013). If we can incorporate these processes into our model, we can examine the effects the environmental effects often included in statistical quantitative genetic analyses on evolutionary dynamics. In this paper, to keep things simple, we focus on the first two processes only.
Let us start with the assumption that the environmental component of the phenotype is determined solely by random developmen- types (Lynch & Walsh, 1998). The similarity between parental and offspring environmental components of the phenotype can generate the types of dynamics depicted in Figure 2d.
We will give examples of the form and parameterization of these functions below. But before we do, we describe the steps required to implement an EE-IPM.

| Model implementation
We approximate Equation 1 into matrix form: N (t + 1) = D (t) R (t) N (t) (Easterling et al., 2000). First, we need to approximate the bivariate distribution N (A, E, t) by categorizing it into many small bins to generate the column vector N (t). In Table 2, below we use 10,000 bins.  We specify the vector N (t = 1) as bivariate normal with two means (one each for A and E) and a variance-covariance matrix that can be estimated from statistical analyses used in quantitative genetics (Falconer, 1960;Lynch & Walsh, 1998). The fit-  Table 2) are then used to construct a diagonal matrix R (t) describing the fitness of each z = A + E. The matrix R (t) is square and with the same length and width dimensions as the length of N (t) (e.g., 10,000 in our example). 1 Note, also, that the value of R (t) for a value of A = 2 and E = 3 would be the same as the value of R (t) for A = 2.5 and E = 2.5 as both give the same phenotypic trait value of 1 It is computationally faster to work with R (t) as a column vector and using element-wise multiplication but we don't dwell on this here. TA B L E 1 Four approaches have been used in evolutionarily structured models to determine the dynamics of the additive genetic variance. Approach 1 and 2 give indistinguishable results in cases where they have been compared, and these are similar to those obtained in approach 3   Additive genetic variance in each generation can be non-Gaussian Childs et al. (2016) and Coulson et al. (2011) 2. Construct a linear and Gaussian probability density function (typical of a standard IPM) that passes through the point A s , A s , has a slope of 0.5, and generates a constant variance.
Additive genetic variance in each generation can be non-Gaussian Coulson et al. (2017) 3. Generate a Gaussian distribution of breeding values in offspring with a constant variance and a mean equal to the mean of the breeding value distribution postselection (A s ).

Additive genetic variance in each generation is Gaussian
Lande (1982) and Simmonds et al. (2020) 4. Allow selection to erode the additive genetic variance.
Additive genetic variance in each generation can be non-Gaussian Coulson et al. (2017)  The kernel D(A ′ , E ′ | A, E, , t) can also be approximated as a ma- It is not always computationally necessary to construct the matrix D (t) and it can be significantly computationally faster to construct the vector N (t + 1) directly from N s (t) = R (t) N (t). We illustrate this here, because for readers who are unfamiliar with Lefkovitch matrices their construction can be opaque, and we do not have space to elaborate here.

TA B L E 2 Describing a bivariate distribution of and as a column vector
The following algorithm could be used to construct the vector N (t + 1), removing the need to populate the matrix D (t): 1. Calculate the sum n (t + 1) = ∑ N s (t) where n (t + 1) is the population size at birth of the next generation.

Calculate the mean of
column vector of the values of A in the third column of Table 1. 3. Choose your assumption about the effect of selection on the additive genetic variance. We will assume that it remains Gaussian with a constant variance σ AA . (3), generate a Gaussian probability distribution with a mean of A s and a variance σ AA discretized into the number of unique bins used to categorize the distribution of A. Call this N short A (t + 1).

A
(t + 1) by the number of unique bins used to categorize E (4th column, Table 1) to generate N long A (t + 1). This will generate the 'blocks' of breeding values with the same value depicted in Table 2 (third column?).
6. Standardize the vector produced in (5) to sum to unity.

E
(t + 1) to sum to unity.
The code to implement these algorithms is provided on Zenodo.

| Choice of parameters
In this section, we describe the two models we use in this paper. The parameters we choose are simply for illustration, and the general results we report are not specific to the parameterization although the specific rates are. We do not identify parameter values from statistical analyses. Simmonds et al. (2019), Simmonds et al. (2020) show how parameters from the analysis of data can be used to parameterize EE-IPMs.
We use 500 bins for both A and for E and upper and lower integration limits for each of 0 and 40. We define the initial bivariate distribution of N (A, E, t = 1) as Gaussian with A = E = 18 and In both models, we specify a constant fitness function R

| RE SULTS
When the mean environment deteriorates with time -perhaps due to human-induced environmental change -and this trend influences the value of E (model 2), evolution (defined as change in A ) occurs at a faster rate than when E remains constant with time (model 1; Figure 3a). In our example, evolution is consequently faster in a human-induced deteriorating environment than in a constant one. The reason for this is that the mean phenotype z = A + E changes more slowly when E trends downwards than when it does not ( Figure 3b). In both our models, the additive genetic variance, the variance in the environmental component of the phenotype, and the phenotypic variance remain constant with time ( Figure 3c).
However, the population growth rate (which equals mean lifetime reproductive success) evolves much more slowly when the environment deteriorates over time compared to when it does not. It is the contrasting dynamics of the population growth rate that generates the difference in rates of evolution between the two models, with the difference in the population dynamics driven by the temporal dynamics of E. All other aspects of the two models are identical.

| DISCUSS ION
Our aim here is to make evolutionarily explicit IPMs accessible to readers who do not have a background in structured population modelling and to explore how a deteriorating environment that mimics the effect of human-induced biotic or abiotic change influences evolutionary dynamics. We have done this by (i) providing background that has not been explicitly described in previous papers using EE-IPMs and (ii) introducing very simple models, one of which includes the effects of a deteriorating on the environmental component of the phenotype. Nonetheless, even these simplified models provide interesting insight.
In both our models, we have a constant, linear, fitness function.
The only difference between our two models is that one contains a deteriorating environment designed to mimic human-induced environmental change that impacts the mean of the environmental component of the phenotype, while the other does not. Such effects of the mean environment on the mean value of phenotypic traits are well-documented in statistical analyses used by statistical quantitative geneticists (Kruuk, 2004;Kruuk et al., 2002;Wilson et al., 2006) but are rarely incorporated into predictive models (Morrissey et al., 2010). In this case, Δz ≠ ΔA due to the trend in E. However, in addition, the trend in E generates divergence in the dynamics of selection between the two models, and this leads to a difference in the rate of evolution.
A selection differential on a phenotypic trait can be written as cov(z,w) w where w, in our model, is absolute lifetime reproductive success and w is mean lifetime reproductive success. For an annual life history, w is the population growth rate (Fisher, 1930). In both our models cov (z, w) remains constant with time -it is determined by the slope of our linear fitness function. However, w changes with time at different rates between our two models. The reason for this is we have a constant fitness function: mean fitness consequently increases with the mean of the phenotype. When there is a trend in E with time, there is consequently a trend in the mean phenotype with time, and hence, mean fitness will change at a different rate compared to when there is no trend in E. When the direction of selection is positive, a negative trend in E will accelerate evolution, while a positive trend will slow it via its effect on mean fitness. For example, when nongenetic inheritance is adaptive it will accelerate the rate of change in mean fitness, and consequently decrease the selection differential, thus slowing the rate of evolution. In contrast, when nongenetic inheritance is maladaptive, it will decrease the rate of change in mean fitness and consequently accelerate the rate of evolution ( Figure 4, and see also Coulson et al., 2017).
The dynamics of selection are rarely decomposed into the dynamics of the covariance between phenotypic traits and absolute fitness and the dynamics of mean fitness. However, doing this does allow some useful insight. For example, because evolution simultaneously alters the mean value of phenotypic traits and mean fitness (Fisher, 1930), it should not be assumed that selection is constant with time when making evolutionary predictions with a constant fitness function. Our results show that the dynamics of E -typically ignored in traditional quantitative genetic approaches, but potentially important when investigating human-induced evolution -can change the denominator of the selection differential, by modifying the rate of change of mean fitness. Selection differentials can vary with time due solely to evolution of the population growth rate (Pelletier & Coulson, 2012).
Our models are deliberately very simplistic. In real settings, fitness functions are likely to include environmental variation and density dependence (Ellner et al., 2016;Simmonds et al., 2019Simmonds et al., , 2020. We also include only one phenotypic trait, but selection operates simultaneously on multiple traits (Lande & Arnold, 1983). Finally, mean fitness (the population growth rate) will fluctuate with time within a generation in iteroparous species (Coulson et al., 2005). EE-IPMs have been constructed for multivariate phenotypic traits, for iteroparous species, and in both variable and deteriorating environments Coulson et al., 2017;Simmonds et al., 2020). A wide range of more realistic settings on evolutionary dynamics can consequently be examined. In addition, IPMs can be used to simultaneously study not only evolutionary, phenotypic trait and population dynamics, but also the dynamics of life histories and interacting species (Bassar et al., 2017;Childs et al., 2004;Coulson et al., 2011;Ellner et al., 2016;Rees et al., 2014). These models are consequently quite flexible and can also be used to study ecoevolutionary feedbacks and may be particularly relevant for humaninduced deterioration in the environment. Moreover, they are easily parameterized from data routinely used to conduct statistical quantitative genetic analyses and explicit genotype-by-environment interactions (where the environmental component of the phenotype is impacted in different ways by environmental change within different genotypes) can be easily incorporated.
Despite the positives of IPMs, they are not a panacea. To date, no one has constructed EE-IPMs for environment-specific phenotypic traits. Traits that are only expressed at specific ages have been incorporated into models , and similar logic could be used for environment-specific traits (Wilson et al., 2006). Second, quantitative geneticists often treat fitness as a phenotypic trait and are interested in the additive genetics of fitness. The evolution of fitness that is not coupled directly via fitness function to a specific trait has not yet been incorporated into an IPM and doing so will not be entirely straightforward, but is theoretically feasible. But perhaps the biggest limitation of IPMs is they do become computationally cumbersome as the size of the multivariate distribution being modelled increases (Ellner et al., 2016). Once the number of dimensions exceeds 6-10, high-performance computing may be required to iterate models unless some way of avoiding multiplying large matrices together can be found (as we demonstrate in our models).
There are, of course, many other structured populations models that have been developed to examine evolutionary dynamics (Barfield et al., 2011;Charlesworth, 1994;Chevin et al., 2010;Lande, 1982) and nonstructured models assuming normality of the additive genetic variance (e.g., approach 4 in Table 1). Some consider reaction norm approaches to examine the effects of environmental change on dynamics (Lande, 2009). These valuable approaches have rarely been parameterized for real systems from statistical quantitative genetic analyses, and they do not link to explicit environmental drivers such as climatic variation as EE-IPMs have been (Simmonds et al., 2020). Instead, models have assumed that different breeding values express different phenotypes in contrasting environments, without the driver of the contrast necessarily being included in models (Lande, 2009). The major difference between our approach and these other theoretical models, is that we explicitly model the dynamics of the environmental component of the phenotype and how it is impacted by environmental variation. However, in doing this, our approach is more intuitive, as it is straightforward to decompose the effects of environmental change on population, phenotypic trait and F I G U R E 4 A hypothetic example of evolution in bivariate space helping summarize our results. The diagonal lines represent constant phenotypic trait value clines, with the darker colour representing larger trait values (and when fitness is directional and positive) higher fitness. Because the additive genetic variance equals the variance in the environmental component of the phenotype, the vector describing the direction of selection is at 45 degrees (red line). In both our models, selection is in this direction but the strength varies over time. The types of evolutionary dynamics model 1 produces are depicted by the purple arrow. In our second model, evolution is partly cryptic. When change in the dynamics of the mean breeding value is completely offset by nonadaptive change in the environmental component of the phenotype evolution is fast (not depicted), and the phenotypic trait does not change (remains on the same diagonal line) but its components change in opposing directions. The green line would represent a case where there is no additive genetic variance and all phenotypic change is attributable to the dynamics of the environmental component Constant E evolutionary dynamics, and on the feedbacks between these processes while still being consistent with the evolutionary assumptions incorporated into other modelling frameworks (e.g., Barfield et al., 2011;Charlesworth, 1994;Chevin et al., 2010;Lande, 1982).
It is our belief that structured models and statistical quantitative genetics are both powerful tools to study evolution. There are ways these approaches can be combined, and once they are they offer potential to shed light on evolutionary dynamics. Investigating both the statistical quantitative genetic and structured modelling literature is time consuming given that both are large, specialized and use different vocabularies. Nonetheless, collaboration rather than distrust between researchers in each discipline could pay dividends.

ACK N OWLED G EM ENTS
Tom Potter is joint funded by NERC DTP and Lamb and Flag studentships at Oxford University, and Anja Felmy is funded by an Early Postdoc Mobility Fellowship from the Swiss National Science Foundation (P2EZP3_181775). We thank Dylan Childs and Joe Travis for helpful comments on an earlier version of the manuscript, and Jarrod Hadfield and two anonymous reviewers who provided extremely useful reviewer comments.

CO N FLI C T O F I NTE R E S T
There are no conflicts of interest to declare.

DATA AVA I L A B I L I T Y S TAT E M E N T
There are no data used in this paper. Code to run models will be made available on Zenodo following acceptance.