Integrating eco‐evolutionary dynamics into matrix population models for structured populations: Discrete and continuous frameworks

State‐structured populations are ubiquitous in biology, from the age‐structure of animal societies to the life cycles of parasitic species. Understanding how this structure contributes to eco‐evolutionary dynamics is critical not only for fundamental understanding but also for conservation and treatment purposes. Although some methods have been developed in the literature for modelling eco‐evolutionary dynamics in structured population, such methods are wholly lacking in the G function evolutionary game theoretic framework. In this paper, we integrate standard matrix population modelling into the G function framework to create a theoretical framework to probe eco‐evolutionary dynamics in structured populations. This framework encompasses age‐ and stage‐structured matrix models with basic density‐ and frequency‐dependent transition rates and probabilities. For both discrete and continuous time models, we define and characterize asymptotic properties of the system such as eco‐evolutionary equilibria (including ESSs) and the convergence stability of these equilibria. For multistate structured populations, we introduce an ergodic flow preserving folding method for analysing such models. The methods developed in this paper for state‐structured populations and their extensions to multistate‐structured populations provide a simple way to create, analyse and simulate eco‐evolutionary dynamics in structured populations. Furthermore, their generality allows these techniques to be applied to a variety of problems in ecology and evolution.


| INTRODUC TI ON
Many populations in biology are structured in ways such as age, physiology or habitat. This state-structure may have large impacts on the eco-evolutionary dynamics of populations (Tuljapurkar & Caswell, 1997). For example, it is well known that the age-structure of human societies plays a large role on long-term population trends: an aging population with the majority of its members at postreproductive ages will experience different population dynamics than one in which most individuals are in a pre-reproductive stage (Micó et al., 2006). These trends in turn inform policy making decisions such as allocation of investments in schooling or social security (Herrero et al., 2019;Lutz & Samir Kumar, 2013). As such, it is critical to develop modelling tools to examine the eco-evolutionary dynamics of structured populations.
Matrix population modelling (Caswell, 2001) is well-suited to study these sorts of problems. The goal of matrix population models is to understand the ecological dynamics in populations that are structured in a discrete way, such as by age or habitat classes, by capturing the transitions among key stages in the life cycle of organisms and projecting the long-term dynamics of such structured populations. Although not discussed further here, continuous state variables such as size can be incorporated by using an integral projection model approach (Merow et al., 2014;Rees et al., 2014). Evolutionary dynamics do not play a central role in matrix population and integral projection models as they do in evolutionary game theory. Yet, methods have been developed to understand the evolutionary dynamics in these models via analysis of selection gradients of asymptotic growth rates to perturbations of matrix entries via sensitivity and elasticity measures (Caswell, 2012;Caswell et al., 2018;Caswell & Salguero-Gómez, 2013;Caswell & Shyu, 2012) or adapted invasion analyses for structured populations (Barfield et al., 2011;Knight et al., 2015;Metcalf et al., 2015;Shefferson et al., 2014). By drawing on matrix population modelling techniques, we hope to expand the scope of the G function framework, previously confined to unstructured populations, to structured ones. Simultaneously, although the evolutionary tools in matrix population modelling are powerful and broadly applicable across a wide range of problems, they do not lend themselves easily to direct ODE simulation sensu the G function framework. Thus, there is a critical need to merge techniques from evolutionary game theory with those of matrix population theory to create a simple, easy-to-use framework within which to investigate a wide range of problems pertaining to eco-evolutionary dynamics in structured populations.
In this paper, we create such a modelling approach under the purview of G functions Vincent & Brown, 2005) and matrix population theory (Caswell, 2001;Tuljapurkar & Caswell, 1997). We provide an exposition of the theoretical grounding behind this approach, focusing on equilibrium concepts and stability. We then extend this eco-evolutionary framework to multistate-structured populations, that is populations that are structured in more than one way, and provide a tool for analysing such models. This leads to multistate population matrix models, also called hyperstate matrix models or age x stage models, a sparsely investigated area altogether.

| BACKG ROU N D: G FUN C TI ON S AND MATRIX P OPUL ATION MODEL S
To set the stage for our structured eco-evolutionary modelling approach, we first introduce the key tenets of G functions and matrix population modelling theory. For simplicity of exposition, we focus on the case of a single evolving continuous strategy in a monomorphic, discretely structured population.

| G functions
G functions use ordinary differential equations (ODEs) or difference equations (DEs) to simultaneously model the ecological (population) and evolutionary (strategy) dynamics of biological populations Vincent & Brown, 2005). In doing so, it lifts the restriction of discrete trait values characteristic of traditional matrix games, and allows for continuous, gradualistic evolution. The G function approach is closely related to the more well-known adaptive dynamics framework (Brännström et al., 2013;Diekmann, 2004;Kisdi & Geritz, 2010), individual fitness-function approach (Cohen et al., 1999), evolutionary distributions (Cohen, 2003a(Cohen, , 2005(Cohen, , 2011(Cohen, , 2003b, and fast-slow dynamical systems as applied to evolution (Cortez & Ellner, 2010).
The G function framework is built upon Darwin's tenets of natural selection: heritable variation, a struggle for existence, and heritable variation influencing the struggle for existence (Gause, 1935;Mallet, 2012;Turchin, 2001). Heritable variation is captured by allowing each individual to have a heritable strategy, v ∈ U, where U captures the set of strategies for each state allowed in the model.
Since we are dealing with a monomorphic structured population, the struggle for existence is incorporated through a fitness generating function, G(v, u, x), which defines the expected per capita growth rate of an individual as influenced by its own strategy across states, v, the mean strategy across each state in the life cycle, u, and the population densities of members in each of these states, x. Note that under this formulation, v, u and x are all vectors of dimension s, where s is the number of states in the life cycle.
This framework can be expanded to a polymorphic population, in which case u becomes a vector of vector-valued strategies, with u i ∈ u capturing the strategies across states for morph i in the population. Finally, the influence of heritable variation on the struggle of existence is given by the dependence of the G function on the individual's strategy, v.
With this understanding, we can now outline the key dynamics of the G function framework. We start with the population dynamics. Since the G function is defined as the per capita growth rate (or the finite rate of increase in the discrete case), to obtain the population dynamics within a life cycle state, we simply multiply it by the current population size as follows: Next, we build an equation for the strategy dynamics to determine how the population's trait value in each state evolves over time. Intuitively, this depends on two factors: (1) how the fitness generating function due to perturbations in trait values (i.e. the local gradient of the G function) and (2) how fast species can scale this fitness gradient (i.e. their evolvability). Mathematically, this is represented as follows: The natural logarithm in the discrete evolutionary equation is the result of scaling between discrete and continuous time (Van Tienderen, 2000;Vincent & Brown, 2005). Including v in the fitness generating function allows us to model individual selection. Its absence leads to frequency independent formulations among organisms with the same strategy (Roughgarden, 1976) and leads to group selection in a frequency dependent context (Abrams, 1987;Brown & Vincent, 1987;Taper & Case, 1992). For a more detailed dive into the G function methodology, we refer the reader to Vincent & Brown, 2005).

| Matrix population models
Matrix population modelling dates back to the 1940s and uses the power of matrix theory to model transient and asymptotic population dynamics. Although not restricted to modelling structured populations, matrix population models are particularly effective at integrating population dynamics and population structure. This is done by encoding transitions between states in a matrix and making projections for how the population size and state distribution change over time. The term "projection" is used to stress the fact that the matrix simply projects the current state of the system forward in time-it does not predict the future state of the system. Although initially developed for simple, linear, age-structured populations, the scope of matrix population models has expanded in recent decades to encompass various aspects of stage-structure, data-driven modelling, statistical inference, evolutionary demography, density dependence, periodic and stochastic environments, integration with population genetics, quantitative genetics, invasion analysis theory, and their application to solve problems in conservation and management (Caswell, 2001). Here, we present the basics of matrix population modelling in the context of structured population modelling.
The first step to create matrix population models is to identify the key states in the population and quantify transition rates or probabilities between these states. This is commonly done by constructing a life cycle graph and then translating this into a population projection matrix (Fujiwara & Diaz-Lopez, 2017;Hansen, 2009;Lefkovitch, 1965), similar to the extraction of adjacency matrices from finite graphs. The transitions between states can either be constants, often derived from demographic data, or equations, derived from mechanistic biological understanding. Equivalently, we can start from a set of ODEs/DEs that captures state transitions in a population and convert this into a population projection matrix. This process can be seen in Figure 1. we assume that A is diagonalizable, a biologically reasonable assumption since most structured populations in nature produce matrices that can be diagonalized. Thus, we can perform a spectral decomposition to derive: where i represents the i th eigenvalue of A, w i is the corresponding right eigenvector, and c i is a constant chosen such that x(t) = ∑ i c i w i . The long-term behaviour of A depends on the spectral radius (discrete (2) (3) is the spectrum of A.

Definition 2. (Spectral
is the spectrum of A, and ℜ( ) denotes the real part of .
As t → ∞, the asymptotic growth rate of the population is con- we will soon see, the spectral radius and bound provide convenient measures for long-term per capita growth rates that we can incorporate into the G function framework. For a more expansive exposition of matrix population models, we refer the reader to (Caswell, 2001;Tuljapurkar & Caswell, 1997).

| ECOLOG IC AL DYNAMIC S
To integrate matrix population modelling into the G function framework, we begin with the ecological (population) dynamics. We let A be the population projection matrix that captures transitions between states. Consequently, the population dynamics of the structured population follows Equation 7. This leads us to a natural definition of an ecological equilibrium.

Definition 3. (Ecological
Although Equation 7 tracks the number of individuals in each state over time, we are also interested in changes to the size of the entire population. This will become particularly important when we incorporate evolutionary dynamics in Section 4. Based on results in Section 2.2, the spectral radius and bound provide good measures of long-term per capita growth rates of the population, analogous to the fitness generating function in the G function framework. This leads us to the following definition:

state-structured populations, the G function is defined as
Note that our formulation will not be used to simulate population dynamics. However, it gives us crucial information on the dynamics of the entire population. Namely, it allows us to derive a necessary condition for the ecological equilibrium as follows. In line with (Metz et al., 2008), we choose to stress the limitations of this approach rather than "advertising positive predictions." The choice of fitness measure is a hotly contested one in the literature, with some suggesting that eco-evolutionary predictions must be based on individual lifetime reproductive success R 0 and others suggesting that the intrinsic rate of increase, r is a better proxy (Caswell, 2001;Charlesworth, 1994;Charnov, 1993;Roff, 1993;Stearns, 1992). However, work by Metz, Mylius, and Dieckmann shows that such optimization approaches (including the one we take here) are limited in scope. Namely, they show that r and R 0 are strictly appropriate proxies only when density (or frequency) dependence operates in a state-independent manner on survival and fecundity, respectively (Mylius & Diekmann, 1995).
In further work, they prove that any scalar-fitness measure is attainable only when frequency or density dependence operates in an effectively one-dimensional monotone fashion (Metz et al., 2008).
Ultimately, a careful invasion fitness approach must be employed when dealing with more complex life histories.

Lemma 1. (Ecological equilibrium characterization). At
Proof. It follows directly from the explanation outlined in the previous section. At equilibrium, we re- With these basic ecological equilibrium concepts, we can examine the local stability of the equilibria. First, we define two key stability concepts that we will consider: Lyupanov stability and asymptotic stability.

Definition 5. (Stability Concepts). An ecological
If can be made arbitrarily large, the equilibrium is globally stable; otherwise, it is locally stable.
Intuitively, Lyupanov stability requires that solutions that start near x stay near x. Asymptotic stability, which we will use in this paper, requires that solutions that start near an equilibrium eventually converge to it. These concepts naturally lead to the following definition of a stable ecological equilibrium.

Definition 6. (ESE).
An ecological equilibrium x that is locally (globally) asymptotically stable is called a locally (globally) ecological stable equilibrium (ESE).
The stability of an equilibrium point of a nonlinear dynamical system as is given by Equation 7 can be assessed by linearizing the system around the equilibrium point, given by the Jacobian matrix.
This gives the following.

Theorem 2. (ESE Theorem). If an ecological equilib-
Jacobian matrix associated with the linearization of the ecological system around x, defined as Proof. To assess the stability of our nonlinear system, we perform a perturbation analysis by Taylor expanding around x and examining the resulting dynamics of the perturbation vector. Here, we perform the analysis in the discrete-time case. The proof for the continuous-time case is given in the Supporting Information. We define a perturbation vector of dimension s: where A x denotes A evaluated at the ecological equilibrium. Expanding, simplifying, and ignoring higher order terms gives us: We can further simplify this by rewriting the second term on the right-hand side as ∑ is a square matrix with x in its i th column and zeros everywhere else. Thus, we have We note that that the Jacobian, J = controls the dynamics of the perturbation vector.
Since the Jacobian is (4)

| E VOLUTIONARY DYNAMIC S
Before describing how our framework can be expanded to include evolutionary dynamics, we briefly mention alternative methods to modelling evolution in structured populations. One such approach focuses on the frequency dynamics of a small number of genes that have large phenotypic effects. This work arguably originated from Fisher's "Malthusian parameter" (Fisher, 1930), which assumed a constant population size, endowing genotypes with fixed per-capita birth and death rates. This work was extended by many authors in the context of agestructured modelling (Charlesworth, 1970(Charlesworth, , 1993(Charlesworth, , 1994Charlesworth & Williamson, 1975;Moran, 1962;Pollak & Kempthorne, 1970, 1971Rousset, 2004), often using difference equation formulations rather than matrices to study selection processes, mostly for notational convenience. This work was further extended to stage-structured population genetics models to probe clonal reproduction and the evolution of senescence (Orive, 1995) as well as deriving invasion criteria for variable environments (Tuljapurkar, 1982). where B is a matrix with the evolutionary dynamics per state normalized by the current trait value of that state on the diagonal entries and zeros everywhere else, We first consider the evolutionarily stable strategy (ESS). The ESS is the central equilibrium concept in much of evolutionary game theory and is often defined as an equilibrium that is uninvadable by alternative strategies not of the ESS. In other words, when common, individuals with this strategy must have a higher fitness than any other rare mutants that exist in the population (Smith & Price, 1973). However, under the G function framework we additionally require convergence stability: natural selection must favour strategies robust to perturbations in v, u and x Eshel, 1983Eshel, , 1996Vincent et al., 1993;Vincent & Brown, 1984). In the words of Christiansen, "the concept of convergence stability may therefore be said to define an evolutionarily attainable stable trait" (Christiansen, 1991).
Mathematically, we can define and characterize an ESS as follows.

Definition 7. (ESS). The strategy associated with x is
an ESS if for any m ∈ U such that m ≠ u, x is an ESE. If x is a local (global) ESE, the ESS will be local (global).

Theorem 3. (ESS Characterization). If the strategy û is an ESS, then
Necessary conditions for an ESS are given by the ESS Maximum Principle in the G function framework Vincent et al., 1993). This principle ensures that the system is at equilibrium, both ecologically and evolutionarily. It also requires the strategies to be at maximum on the adaptive landscape (Cohen et al., 1999;Pigliucci, 2008;Svensson & Calsbeek, 2012). Here, we present the ESS maximum principle for structured populations.

Theorem 4. (ESS Maximum Principle). If û is an ESS,
then the following must be true, Proof. Part (i) follows directly from Lemma 1 and ensures that the system is at an ecological equilibr ium.
Part (ii) checks if the selection gradient of the fitness generating function vanishes at û,x , thereby ensuring that the system is at an evolutionary equilibrium. (7) Part (iii) ensures that the system is at an evolutionary peak of the adaptive landscape and not at an invadable convergent stable minimum that promotes evolutionary branching. ☐ With this ESS definition, we are now ready to define stability concepts for eco-evolutionary equilibria. We will need to expand the notion of asymptotic stability to one of convergence stability, in which the equilibrium is robust to perturbations in v, u and x. where K i is a square matrix with û in its i th column and zeros everywhere else. Putting the ecological and evolutionary perturbations together, we arrive at our desired Jacobian. ☐ Another important evolutionary equilibrium concept is the neighbourhood invader strategy (Apaloo, 1997;Apaloo et al., 2009).

Theorem 5. (Convergence Stability). If an equilibrium
This concept, closely related to the idea of invasion analysis in adaptive dynamics (Brännström et al., 2013;Diekmann, 2004;Kisdi & Geritz, 2010), assesses whether a strategy, û, can invade a resident population with a strategy close to û. We call this "neighbouring strategy" m and its respective population equilibrium ŷ. By definition, individuals with v =m will have a fitness of 0 within this resident population. Thus, for û to invade the resident population, it must have a positive fitness when introduced to the population at low numbers. We can provide a formal definition and characterization of the neighbourhood invader strategy for structured populations:

Definition 8. (Neighbourhood Invader Strategy). A strategy û is an neighbourhood invader strategy if
If can be made arbitrarily large, the neighbourhood invader strategy is global. Otherwise, it is local.

Theorem 6. (Neighbourhood Invader Strategy Characterization). The strategy û is an neighbourhood invader strategy iff it is the unique solution to
The proof of this theorem follows directly from the definition of a neighbourhood invader strategy. It is worth noting that the ESS Maximum Principle and convergence stability results hold similarly for the neighbourhood invader strategy, although the second-order condition differs. Namely, we must assess how the fitness of an individual using û changes with changes in m. If d 2 G dm 2 i > 0 for all states i , then û is a neighbourhood invader strategy candidate. It is worth noting that an equilibrium point can be an ESS without being a neighbourhood invader strategy or vice versa. When a strategy is both an ESS and a neighbourhood invader strategy (called an evolutionarily stable neighbourhood invader strategy Apaloo, 1997), it is convergence stable, can invade neighbouring strategies, and repel invading mutants.
Furthermore, this allows for mutational jumps towards the equilibrium that may not be possible solely with convergence stability. If a strategy is an ESS but not a neighbourhood invader strategy, then the equilibrium, once established, cannot be invaded. However, it may not . (8)

| MULTIS TATE S TRUC TURED P OPUL ATIONS
Thus far, we have discussed populations structured by a single characteristic such as age, stage, or habitat. However, many populations are structured in multiple ways simultaneously. For example, human populations are structured by age and location at the same time, two factors that can greatly impact demographic outcomes. For clarity, we call each set of states (e.g. age and location in this case) a compartment. In this section, we will outline ways to model and analyse the eco-evolutionary dynamics of multi-compartment populations using multistate population matrix models, also known as hyperstate matrix models or as age x stage models.
The theory that we have developed above applies carries over to multi-compartment structured populations. In order to construct multistate matrix population models, the first consideration concerns how compartments interact with each other. Formally, this requires defining the nested compartment structure. For example, is age nested within habitat or is habitat nested within age? Once this is decided, we can construct our multistate matrix population model as a block matrix encoding transitions between states and compartments. In our age × habitat example given above, we can nest age structure within habitat structure to derive the following block matrix for our multistate matrix population model: The age transition block matrix in the upper diagonal describes the transitions among age groups within the first habitat. The migration block matrix in the upper off-diagonal captures migratory transitions from the second habitat to the first habitat. The migration block matrix in the lower off-diagonal represents migratory transitions from the first habitat to the second habitat. The age transition block matrix in the lower diagonal describes transitions among age groups within the second habitat. This matrix can then be used to simulate eco-evolutionary dynamics and perform the same analyses as in the prior sections. One additional tool that we can use to analyse multistate matrix population models is ergodic flow preserving foldings, analogous to Caswell's notion of marginal distributions (Caswell et al., 2018). Ergodic flow preserving foldings serve as a compartment analogue for traditional sensitivity and elasticity analyses on matrix entries (Caswell, 2001) and allow us to understand the importance of each compartment on the ecological and evolutionary dynamics of the population. This method was originally proposed in Enright et al. (1995) and later formalized for age (Hooley, 2018) and stage (Salguero-Gómez & Plotkin, 2015) structured populations. Most recently, these methods were adapted to multistate matrix population models in a purely ecological context (Coste et al., 2017;Coste & Pavard, 2020). Here, we provide an intuitive exposition of ergodic flow preserving foldings for eco-evolutionary analysis, but we refer the readers to (Coste et al., 2017;Coste & Pavard, 2020;Hooley, 2018;Salguero-Gómez & Plotkin, 2015) for more formal presentations.
The idea behind an ergodic flow preserving folding analysis, also known as a trait-level analysis, is to "fold" across a compartment and assess the resulting changes to the eco-evolutionary dynamics. In ergodic flow preserving folding analysis, "folding" refers to merging states within a compartment in a way that preserves asymptotic flows of individuals between the states (Coste et al., 2017;Coste & Pavard, 2020). To do this, we sum the transitions associated with the

| E X AMPLE: LIFE-HISTORY TR ADE-OFFS
To illustrate the process of creating and simulating eco-evolutionary dynamics of structured populations using the theory developed in this paper, we construct and probe a simple toy model of life-history trade-offs. This idealized toy model is not meant to elucidate novel aspects of life history evolution of a particular organism, but to rather serve as an expository example of how to construct, simulate, and analyse structured G function models.
To begin, assume that an organism has two life history states: This life history scenario can be mathematically formalized as follows: Note that b(A) = 0 at the eco-evolutionary equilibrium, as expected.
This ODE system can be formally written as a matrix population model by appropriately placing state transitions into the following matrix: The ESS maximum principle requires that the derivative of b(A), computed below, is 0: Equations 11, 12 and 14 give us a system of three equations and three unknowns. Solving this with the parameter values in Table 1 gives our ESS as We now numerically simulate the model and show that the system does indeed converge to the analytical ESS candidate derived above ( Figure 4). This expository example shows how to construct a system of ODEs to represent a biological problem of interest, how to convert this into a MPM, and how to simulate and analyse the resulting ecological and evolutionary dynamics.

| CON CLUS ION
In this paper, we introduced a framework that draws on evolutionary game theory and matrix population theory to model, simulate, and analyse the eco-evolutionary dynamics of state-structured populations. This involves developing analogues for fitness in state-structured populations and provide definitions of key equilibrium concepts (ecological equilibrium, evolutionarily stable strategy, convergent stability, and neighbourhood invader strategy; Apaloo et al., 2009;Diekmann, 2004;Geritz et al., 1998;Kisdi & Geritz, 2010;McGill & Brown, 2007). We assess and prove properties related to the asymptotic and convergence stability of these systems and develop an extension of the ESS maximum principle  (2023). (Vincent et al., 1996) for structured populations. We expand the framework of evolutionary game theory to include multistate structured populations. A toy model of life history trade-offs involving an organism with juvenile and adult states provides an example for how to develop, analyse, and simulate these models.
Evolutionary dynamics in structured populations is mainly studied by two approaches: sensitivity analyses to measure the impact of key parameters or state transitions on fitness (Caswell, 2012;Caswell et al., 2018;Caswell & Salguero-Gómez, 2013;Caswell & Shyu, 2012) and modifications of invasion analyses from adaptive dynamics to structured populations (Barfield et al., 2011;Knight et al., 2015;Metcalf et al., 2015;Shefferson et al., 2014). The methods developed in this paper complement these techniques by putting eco-evolutionary dynamics in structured populations on firm theoretical ground in the context of the G function framework.
These methods allow for direct simulation of deterministic ecoevolutionary dynamics in the form of ODEs or difference equations and provide simple tools for analysing these models. Along with sensitivity and invasion analyses, our techniques allow us to understand the impact of environmental or manager-induced changes on a biological system, identify the parameters that have the biggest impact on eco-evolutionary dynamics (for policy-making purposes or allocation of sampling and measurement efforts), and understand on a more fundamental and theoretical level, the dynamics of the underlying system (Caswell, 2019). We have also extended our framework to allow for modelling multistate structured populations, a scarcely investigated area, and provide a way to extend ergodic flow preserving folding analysis, which has been restricted to analysis of ecological dynamics, to eco-evolutionary dynamics in a qualitative manner. A major caveat of our proposed optimizationlike framework is that it is strictly only applicable for cases of simple density and frequency dependence. Future work will extend this theoretical framework to include more complex density and frequency dependent interactions, continuously structured populations, and polymorphic populations. We also aim to apply this framework to a host of problems across ecology, evolution, and biomedicine.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no potential conflict of interest.

DATA AVA I L A B I L I T Y S TAT E M E N T
Codes used to produce Figure 2 can be found at the author's