# Toward a unified approach to modeling adaptation among demographers and evolutionary ecologists

Joanie van de Walle, Jimmy Garnier, Timothée Bonnet, Stéphanie Jenouvrier

<figure><img src="image_1.png" /></figure>

## To cite this version:

Joanie van de Walle, Jimmy Garnier, Timothée Bonnet, Stéphanie Jenouvrier. Toward a unified approach to modeling adaptation among demographers and evolutionary ecologists. 2024. hal-04851098

HAL Id: hal-04851098

https://hal.science/hal-04851098v1

Preprint submitted on 20 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
# Toward a unified approach to modeling adaptation among demographers and evolutionary ecologists

Joanie Van de Walle*<fnref n="1" />,<fnref n="2" />, Jimmy Garnier*<fnref n="3" />, Timothée Bonnet*<fnref n="4" />, and Stephanie Jenouvrier*†<fnref n="1" />

**Running headline:** Evo-Demo Hyperstate Matrix Population Model

**Corresponding author:** Stephanie Jenouvrier, sjenouvrier@whoi.edu

## Author's contributions

Timothée Bonnet (TB), Jimmy Garnier (JG), Stephanie Jenouvrier (SJ), and Joanie Van de Walle (JV) contributed equally to all aspects of this research. This work reflects a collaborative and interdisciplinary effort, in which all authors were equally involved in: conceptualization, formal analysis, investigation, methodology, visualization and writing. SJ and JV supervised and administrate the project.

## Fundings

This work was supported by the National Science Foundation under Grant No. ORCC NSF 2222057 No NSFGEO-NERC 1951500.

## Data accessibility

All of the MATLAB and R code for building the Evo-Demo Hyperstate Matrix Population Model and the Individual-Based Model are provided on https://gitfront.io/r/fledge-whoi/q4KN7RDbpwUq/Eco-EvoHyperstateModel/.

<footnote n="1">¹Biology Department, Woods Hole Oceanographic Institution, Woods Hole, MA, USA</footnote>

<footnote n="2">²Department of Fisheries and Oceans Canada, Maurice Lamontagne Institute, Mon-Joli, Canada</footnote>

<footnote n="3">³LAMA, CNRS-University Grenoble Alpes, University Savoie Mont-Blanc, UMR-5127, Chambéry, France</footnote>

<footnote n="4">⁴Centre d'Etudes Biologiques de Chizé, CNRS-La Rochelle University UMR-7372, Villiers en Bois, France</footnote>

<footnote>*All authors contributed equally to this work.</footnote>

<footnote>†Corresponding author: sjenouvrier@whoi.edu</footnote>
Abstract

1. Demographic and evolutionary modeling approaches are critical to understanding and projecting species responses to global environmental changes. Population matrix models have been a favored tool in demography, but until recently, they failed to account for short-term evolutionary changes. Evolutionary-explicit demographic models remain computationally intensive, difficult to use, and have yet to be widely adopted for empirical studies. Researchers focusing on short-term evolution often favor individual-based simulations, which are more flexible but less transferable and computationally efficient. Limited communication between fields has led to differing perspectives on key issues, such as how life-history traits affect adaptation to environmental change.

2. We develop a new Evo-Demo Hyperstate Matrix Population Model (EvoDemo-Hyper MPM) that incorporates the genetic inheritance of quantitative traits, enabling fast computation of evolutionary and demographic dynamics. We evaluate EvoDemo-Hyper MPM against individual-based simulations and provide analytical approximations for adaptation rates across six distinct scales in response to selection. We show that different methods yield equivalent results for the same biological scenario, although semantic differences between fields may obscure these similarities.

3. Our results demonstrate that EvoDemo-Hyper MPM provides accurate, computationally efficient solutions, closely matching outcomes from individual-based simulations and analytical approximations under similar biological conditions. Adaptation rates per generation remain constant across species when selection acts on fertility but vary with other vital rates. Adaptation per time decreases with generation time unless selection targets adult survival, where intermediate life histories adapt fastest. Rates per generation, defined as the relative change in individual fitness, remain constant across species and vital rates.

4. We discuss that no general prediction emerges about which species or life-history traits yield higher adaptation rates, as outcomes depend on life cycles, vital rates, and the definition used. We provide Matlab and R code to support the application of our EvoDemo-Hyper MPM.

**Key-words:** Adaptive evolution; Evolutionary demography; Hyperstate matrix model; Life-history strategies; Quantitative genetics; Slow-fast continuum.
# Introduction

Understanding how populations respond to environmental changes is essential for effective conservation. Both demographic and evolutionary processes influence the dynamics of the population, but the integration of these perspectives remains a challenge. Bridging these perspectives is crucial to making accurate predictions about species' responses to global change.

Integral Projection Models (IPMs) based on continuous-variable matrix modeling have been widely used to model population responses to selection acting on continuous phenotypic traits in terms of future phenotypic distribution and stage abundance [Easterling et al., 2000]. IPMs allow efficient computation of population properties, but their original formulation did not account for evolution or include explicit genetic dimensions. Recent advances introduced inheritance and evolutionary dynamics into IPMs [Childs et al., 2016; Coulson et al., 2017]. However, these improvements involve nonlinear iterative steps (e.g., multiple for-loops) that reduce computational efficiency and limit their applicability.

Incorporating adaptive evolution into matrix population models (MPMs) requires additional dimensions to capture both phenotypic traits and genetic inheritance. Multitrait population projection matrices (MPPMs) address genetic variance and multiple vital rates [Coste & Pavard, 2020]. However, they often lack explicit trait inheritance within a quantitative genetic framework, limiting their applicability for evolutionary research. Here, we propose using hyperstate matrix models [Roth & Caswell, 2016], which extend vec-permutation models [Hunter & Caswell, 2005] to efficiently incorporate evolutionary processes and develop a new Evo-Demo Hyperstate Matrix Population Model (EvoDemo-Hyper).

In parallel, evolutionary ecologists and quantitative geneticists have focused on predicting phenotypic changes through genetic processes, using Individual-Based Models (IBMs) [Rees & Ellner, 2019]. Although IBMs can accommodate complex demographic processes, they are time consuming to develop, less transferable, and computationally demanding compared to matrix models.

These two modeling approaches, Individual-Based Models (IBMs) versus matrix-based approaches (MPMs and IPMs) have evolved independently, often appearing in conflict [Coulson et al., 2017; Pelletier, 2019]. Demographers emphasize the structure of the life cycle and the position along the slow-fast continuum, as the vital rates under selection determine adaptive responses [Salguero-Gómez et al., 2016]. Faster-lived species with shorter generation times and rapid population growth rates are expected to adapt more quickly to environmental change [Vedder et al., 2013]. In contrast, evolutionary ecologists focus on predicting adaptation through the variance-covariance properties of traits and individual fitness, using tools such as the breeder's equation, the secondary theorem of selection, or Fisher's fundamental theorem of natural selection [Queller, 2017]. Although these approaches may seem different, quantitative evolutionary genetics inherently incorporates life cycle dynamics [Barfield et al., 2011] and aligns with demographic principles [Fisher, 1930]. However, semantic differences between fields can cause confusion and hinder collaboration.

A key distinction between these frameworks lies in the time units used. Quantitative geneticists typically predict responses to selection per generation [Lynch & Walsh, 1998], whereas demographers focus on population dynamics over intervals corresponding to life-stage transitions, often annually [Caswell, 2000]. This difference comes from the tools used: the breeder's equation predicts changes from generation to generation, while Leslie matrices track life-stage transitions and reproductive events over time. Although converting between these time units is straightforward, it is not always applied in the academic literature, leading to misinterpretations.

Another difference involves how selection is quantified. In quantitative genetics, selection is expressed on the scale of relative fitness, defined as individual fitness divided by mean fitness. In
contrast, demography uses regression parameters to relate traits to vital rates (e.g., survival, fertility), treating the regression slopes as indicators of selection. Additionally, non-Gaussian link functions (e.g., logit or Poisson) complicate the back-transformation of regression parameters to relative fitness scales, adding further discrepancies between the fields [de Villemereuil et al., 2016]. In this study, we integrate perspectives from demography, evolutionary genetics, and mathematics by developing a matrix population model that explicitly incorporates evolutionary processes: EvoDemo-Hyper (Evo-Demo Hyperstate Matrix Population Model). Our model predicts genetic adaptation rates for five species along the slow-fast life history continuum, highlighting how evolutionary responses vary depending on which vital rate is under selection, with life history influencing both generation time and adaptation rates. We compare the predictions of our model with those of a similarly parameterized Individual-Based Model (IBM) and analytical approximations. Additionally, we express adaptation rates using six definitions, capturing different time frames and selection measures across research communities. Our findings indicate that the matrix model, IBM, and analytical approaches give similar results on a common scale, but no general prediction can be made on which species or life-history trait adapts faster; it depends on the life cycle, vital rates under selection, and definition choice.

## Methods

In the following sections, we present methods for the construction of the evolutionary matrix population model: EvoDemo-Hyper (Evo-Demo Hyperstate Matrix Population Model). Appendix S1 describes a corresponding Individual-Based Model. Both models are based on the same basic life cycle from which we simulated five different species with contrasted life histories.

### Life cycle

The stage-structured life cycle (Figure 1) consists of juveniles (J) and adults (A). Annual transitions occur from time $t$ to $t + 1$. The survival probabilities of juveniles and adults are $S_J$ and $S_A$. Juveniles at time $t$ mature and become adults with probability $\gamma$ at time $t + 1$, or remain juveniles with probability $(1 - \gamma)$ at time $t + 1$, both given survival. For a post-breeding census, fertility ($F$) indicates the *per capita* rate of contribution from adults to juveniles annually (Appendix S2). An equal offspring sex ratio at birth is assumed.

<figure><img src="image_3.png" /><figcaption>**Figure 1:** Life cycle graph of the population models</figcaption></figure>
This life cycle simulates species with diverse life history strategies [Neubert & Caswell, 2000]. We model five species with increasing generation time (Table 1), mapping them along the slow-fast life history continuum [Salguero-Gómez et al., 2016]. They vary in age at first reproduction (early vs. late) and lifespan (short vs. long) [Neubert & Caswell, 2000]. We also consider differences in reproductive strategies (semelparous vs. iteroparous) [Neubert & Caswell, 2000; Salguero-Gómez et al., 2016]. These species illustrate a range of lifespan and reproductive variations (Table 1). Species 1 features rapid development, a single lifetime reproductive event, and high fecundity, akin to annual plants and insects. Species 2 is a short-lived iteroparous species with limited offspring per time unit, similar to small mammals and birds. Species 3 spans several years, producing one offspring per time unit, like some deer. Conversely, species 4 and 5 are long-lived with delayed reproduction and low breeding output, resembling small primates or larger species like whales and albatross.

## Phenotypic Variation, heritability, and selective pressure on vital rates

Under a simple quantitative genetic infinitesimal model, the phenotype of an individual, $z$, can be decomposed into a population mean, an additive genetic deviation, termed the breeding value, and an environmental deviation [Falconer & Mackay, 1996]. The change in mean breeding values over time represents the genetic evolution of the quantitative trait considered. The variance in breeding values is the additive genetic variance ($V_a$), whereas the variance in environmental deviations is the environmental variance ($V_E$). The sum of the two is the phenotypic variance of a trait ($V_P$), with $V_P = V_a + V_E$. Narrow-sense heritability, denoted $h^2$, is the proportion of the phenotypic variance that can be attributed to inherited genetic factors, and is calculated as $h^2 = V_a/V_P$.

We explore four selection pathways targeting specific vital rates: juvenile survival $S_J$, adult survival $S_A$, fertility $F$, and maturation $\gamma$. Selection on juvenile survival $S_J$, adult survival $S_A$, or maturation $\gamma$ is modulated by parameter $\beta$ following a logit distribution (Appendix S3, Figure S4). For fertility $F$, selection is modulated by parameter $\beta$ using a Poisson distribution, except for long-lived species with one offspring annually, for which fertility is modeled with logit functions. In both scenarios, parameter $\beta$ indicates the strength of selective pressure on the link scale (logit or log).

## Parameter values

Vital rate values (Table 1) ensured stable equilibrium growth rates for each species ($\lambda = 1$) before selective pressures. We projected the population assuming a phenotypic variance of 1, additive genetic variance 0.2, and selective pressure 0.15. These values for $V_P$ and $V_a$ suggest moderate heritability per studies [Postma, 2014]. The value of $\beta$ isn't directly comparable to empirical standardized selection gradients due to scale differences across species and vital rates. With heritability 0.2, observed adaptation rates per generation align with empirical gradient estimates 0.25, fitting within empirical ranges [Kingsolver et al., 2012]. Results for other parameters are described in Appendix S3. Adaptation rate magnitudes depend on $V_a$ and $\beta$, yet species and vital rate differences are consistent (Figure S11).

## Principle of adaptation rate: theoretical derivations

Adaptation rate is determined by the change in mean phenotype $\bar{z}$ over time $\Delta\bar{z}_t$. We analyze and compare adaptation rates under selective pressures on the four vital rates for each of five species. This involved eight calculation methods, varying by the time unit (year vs or generation) and selection unit.
<table><caption>**Table 1:** Demographic parameters for the five species used in the simulations along with demographic characteristics. $\lambda$ is population growth rate and $T$ is generation time</caption><thead><tr><th>Parameter</th><th>Species 1</th><th>Species 2</th><th>Species 3</th><th>Species 4</th><th>Species 5</th></tr></thead><tbody><tr><td>Fertility (F)</td><td>10.9868</td><td>4.9958</td><td>0.9987</td><td>0.3004</td><td>0.2286</td></tr><tr><td>Adult survival (S<sub>A</sub>)</td><td>0.0365</td><td>0.2000</td><td>0.8000</td><td>0.9300</td><td>0.9500</td></tr><tr><td>Juvenile survival (S<sub>J</sub>)</td><td>0.0965</td><td>0.2500</td><td>0.3850</td><td>0.5050</td><td>0.8000</td></tr><tr><td>Maturation rate (γ)</td><td>0.9000</td><td>0.5720</td><td>0.4000</td><td>0.3000</td><td>0.0700</td></tr><tr><td>λ</td><td>1.0000</td><td>1.0000</td><td>1.0000</td><td>1.0000</td><td>1.0000</td></tr><tr><td>T</td><td>2.0476</td><td>2.3698</td><td>6.3004</td><td>15.7521</td><td>23.9042</td></tr></tbody></table>

**Rate of adaptation per unit of time, RA<sub>T</sub>.** In the Fisher infinitesimal model Fisher [1918], the additive genetic variance stays nearly constant over time, close to $V_a$. As changes in $V_e$ were not modeled, the phenotypic variance also remains stable, close to $V_P$. Thus, we calculate the adaptation rate per time unit RA<sub>T</sub> using Lande's equation. [Bulmer, 1980; Lande, 1979]:

$$RA_T = \Delta \bar{z}_t \approx V_a \delta, \qquad (1)$$

with

$$\delta = \frac{\partial \ln(\bar{\lambda})}{\partial \bar{z}_t} \qquad (2)$$

which can also be written as:

$$RA_T = h^2 V_P \frac{1}{\bar{\lambda}} \frac{\partial \bar{\lambda}}{\partial \bar{z}_t} \qquad (3)$$

In addition, when the phenotypic variance in the population is not too large, the mean growth rate $\bar{\lambda}$ is well approximated by the growth rate $\lambda$ at equilibrium of a monomorphic population with trait $\bar{z}$. The selection gradient in the population is thus approximated using the sensitivity matrix S associated to the growth rate $\lambda$ and the variation of the vital rate with respect to the phenotype, so that $(S = \frac{\partial \bar{\lambda}}{\partial \bar{z}_t})$. The sensitivity matrix S is described by

$$S = \frac{vw^T}{v^T w}, \quad \text{with} \quad v = \begin{pmatrix} S_J \gamma \\ \lambda - S_J(1 - \gamma) \end{pmatrix}, \quad w = \begin{pmatrix} F \\ \lambda - S_J(1 - \gamma) \end{pmatrix}$$

where $v$ and $w$ are the left and right eigenvectors associated with the eigenvalue $\lambda$ of the projection matrix $\tilde{A} = \tilde{U} + \tilde{F}$ at the mean trait $\bar{z}$. The projection matrices $\tilde{U}$ and $\tilde{F}$ are defined by the life cycle in Figure 1 (see equations (13),(16)).

**Rate of adaptation per generation unit, RA<sub>G</sub>.** The evolution of the mean phenotype can also be measured per generation time $T$. For a monomorphic population, generation time is defined by Bienvenu & Legendre [2015]

$$T = \frac{\lambda v^T w}{v^T R w} = \frac{\lambda v^T w}{v_1 w_2 F} \qquad (4)$$

The rate of adaptation per generation time is linked to the rate of adaptation per unit of time as follows

$$RA_G = RA_T T \qquad (5)$$

Here, the sensitivity matrix S is replaced by the following sensitivity matrix per generation $S_T$ defined by

$$S_T = S_T = \bar{\lambda} \frac{vw^T}{v^T R w} = \bar{\lambda} \frac{vw^T}{v_1 w_2 F}$$
**Rates of adaptation per unit of vital rate, RA<sub>TSθ</sub> and RA<sub>GSθ</sub>.** The rates of adaptation can also be measured in terms of units of vital rate, i.e., based on the response in the vital rate on which selection acts, denoted by RA<sub>TSθ</sub>, and RA<sub>GSθ</sub> when expressed in terms of generation time. They are calculated as follows

$$RA_{TS\theta} = \frac{RA_T}{V_{\theta}} \quad \text{and} \quad RA_{GS\theta} = \frac{RA_G}{V_{\theta}}, \quad \text{where} \quad V_{\theta} = \frac{\partial \ln(\theta)}{\partial \bar{z}_t} = \frac{1}{\theta} \frac{\partial \theta}{\partial \bar{z}_t}. \qquad (6)$$

where $V_θ$ is the relative rate of change in the vital rate $θ$ of interest.

**Rates of adaptation per unit of fitness.** The selection units can be characterized as the proportional rate of variation in an individual's fitness. Individual fitness can be measured as the average lifetime reproductive success, or the net reproductive rate $R_0$, which is defined for our simple life cycle as

$$R_0 = \frac{FS_J\gamma}{(1 - S_A)(1 - S_J(1 - \gamma))}. \qquad (7)$$

The rate of adaptation over time, denoted by RA<sub>TSR<sub>0</sub></sub> and the rate of adaptation per generation, denoted by RA<sub>GSR<sub>0</sub></sub>, measured per unit of individual fitness are given by

$$RA_{TSR_0} = \frac{RA_T}{V_{R_0}} \quad \text{and} \quad RA_{GSR_0} = \frac{RA_G}{V_{R_0}}, \quad \text{where} \quad V_{R_0} = \frac{\partial \ln(R_0)}{\partial \bar{z}_t} = \frac{1}{R_0} \frac{\partial R_0}{\partial \bar{z}_t} \qquad (8)$$

and $V_{R_0}$ is the relative change in individual fitness.

## Evo-Demo Hyperstate Matrix Population Model (EvoDemo-Hyper)

Our population model, based on the hyperstate model formulation Roth & Caswell [2016], includes three dimensions ($m = 3$) (Figure 2). Individuals are classified by stage ($i$), breeding value ($j$), and phenotype ($k$). The first dimension ($i$) ranges from 1 to $s$, the second ($j$) from 1 to $b$, and the third ($k$) from 1 to $p$. This allows tracking of breeding values and phenotypes. The life cycle has two stages ($s = 2$), with 40 classes for both breeding value and phenotype ($b = p = 40$). Breeding values and phenotypes are assumed to be centered at zero, ranging from -4 to 4.

The population vector $n_t$ consists of the number of individuals $n_{i,j,k}$ in each stage $i$, categorized by breeding value class $j$ and phenotype class $k$. Phenotypes are grouped within breeding values, and breeding values are grouped within stages. Our model projects this population vector from time $t$ to $t + 1$ using the projection matrix $\tilde{A}[n_t]$. Although the model is centered on females, it extends to a two-sex formulation to include the distribution of male breeding values, thus accounting for genetic transmission, assuming a similar structure for both male and female populations. The dynamic of the population is given by

$$n_{t+1} = \tilde{A}[n_t]n_t \qquad (9)$$

The matrix $\tilde{A}$ can be broken down into the matrices $\tilde{U}$ and $\tilde{F}$ (Figure 2), which represent the transition of living individuals and the production of offspring, respectively.

$$\tilde{A}[n_t] = \tilde{U} + \tilde{F}[n_t] \qquad (10)$$

The matrices $\tilde{U}$ and $\tilde{F}$ can also be broken down into sub-processes (Figure 2), which represent different transition processes happening within each dimension:

$$\tilde{U} = (K_{b,p}K_{s,b})^T P K_{b,p} B K_{s,b} U \qquad (11)$$

$$\tilde{F}[n_t] = (K_{b,p}K_{s,b})^T M K_{b,p} H[n_t] K_{s,b} R \qquad (12)$$
<figure><img src="image_4.png" /><figcaption>**Figure 2:** Schematic representation of the $\tilde{A}$ matrix and its constituents. Note that $\emptyset$ means that those matrices contain no specific parameters; they are identity matrices</figcaption></figure>

The matrices $\mathbf{P}$, $\mathbf{B}$, $\mathbf{U}$, $\mathbf{M}$, $\mathbf{H}$ and $\mathbf{R}$ consist of block diagonal matrices with matrices $\mathbf{P}_{ij}$, $\mathbf{B}_{ik}$, $\mathbf{U}_{jk}$, $\mathbf{H}_{ij}$, $\mathbf{M}_{ik}$ and $\mathbf{R}_{jk}$ placed on the diagonal, respectively. The $\mathbf{K}$ matrices are vec-permutation matrices that rearrange the structure of the $\mathbf{n}$ vector after each sub-process is applied, in order to align with the organization in the subsequent processes.

Equation (11) involves a series of operations on the population vector $\mathbf{n}$. First, it is multiplied by the block matrix $\mathbf{U}$, which represents the stage within the breeding value within the phenotype. The resulting vector is then rearranged using the matrix $\mathbf{K}_{s,b}$. Next, it is multiplied by the matrix $\mathbf{B}$, which represents the breeding value within stage within the phenotype. Again, the resulting vector is rearranged, this time using the matrix $\mathbf{K}_{b,p}$. Finally, it is multiplied by the matrix $\mathbf{P}$, which represents the phenotype within the stage within the breeding value. The vector $\mathbf{n}$ is then rearranged back to its original structure, i.e., stage within breeding value within phenotype, using the transpose of the product of matrices $\mathbf{K}_{b,p}$ and $\mathbf{K}_{s,b}$.

In Equation (12), the population vector $\mathbf{n}$ is subject to a comparable sequence of operations, although it involves distinct processes $\mathbf{R}$, $\mathbf{H}[\mathbf{n}_t]$, and $\mathbf{M}$ that serve, respectively, as substitutes for $\mathbf{U}$,
B, and P. Each of these subprocesses is detailed in the following sections.

## Transition of living individuals

The transitions of living individuals from time t to time t + 1 are captured by three sets of matrices, each of which represents a distinct process along a single dimension.

1. The stage transition matrices $U_{jk}$ (dimension $s \times s$) capture the transition of individuals between juvenile and adult stages for each combination of breeding value j and phenotype k. Transitions between stages include survival S and maturation probabilities $\gamma$;

$$ U_{jk} = \begin{pmatrix} S_{J_{jk}}(1 - \gamma_{jk}) & 0 \\ S_{J_{jk}}\gamma_{jk} & S_{A_{jk}} \end{pmatrix} \quad (13) $$

where $S_J$ and $S_A$ are respectively the juvenile and adult survival rate corresponding to the breeding value j and phenotype class k.

2. The transition matrices for stage i and phenotypic class k, $B_{i,k}$ (dimensions $b \times b$) contain transitions of the breeding value class for individuals in each stage i and phenotypic class k. Breeding values are assumed fixed from birth, hence $B_{i,k}$ are identity matrices:

$$ B_{ik} = \begin{pmatrix} 1 & 0 & \dots & 0 \\ 0 & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \dots & 0 & 1 \end{pmatrix} \quad (14) $$

3. Transition matrices for stage i and breeding value class j, $P_{i,j}$ (dimensions $p \times p$) contains transitions between phenotypic classes for individuals of stage i and in breeding value class j. Individuals are assumed to stay within their phenotypic class throughout their life, so that the matrices $P_{i,j}$ are also identity matrices:

$$ P_{ij} = \begin{pmatrix} 1 & 0 & \dots & 0 \\ 0 & \ddots & \ddots & \vdots \\ \vdots & \ddots & \ddots & 0 \\ 0 & \dots & 0 & 1 \end{pmatrix} \quad (15) $$

However, this assumption can be relaxed to account for phenotypic changes over an individual's lifetime, such as traits that improve with age.

## Production of new individuals

The production of new individuals with their own breeding values and phenotypes is represented by the following set of matrices:

1. Reproduction matrices ($R_{jk}$) contain the state-specific fertilities of individuals within the breeding values and phenotypic classes j and k, respectively:

$$ R_{jk} = \begin{pmatrix} 0 & F_{jk} \\ 0 & 0 \end{pmatrix} \quad (16) $$

where F is the fertility rate of adults carrying breeding value j and phenotype k.
2. The breeding values of newly produced offspring are determined by the matrices $H_{i,k}[\mathbf{n}]$ (dimensions $b \times b$). Those matrices contain the probabilities that newly produced offspring receive a given breeding value depending on the breeding values of their parents and are assumed fixed across all $i$ and $k$ values.

Individuals are assumed to reproduce sexually with random mating. The breeding values of the parents are denoted $x_{jf}$ and $x_{jm}$. The breeding value of the offspring is sampled from a normal distribution with a mean equal to the average of the parental breeding values $(x_{jf} + x_{jm})/2$ and a variance equal to $V_a = h^2V_P$. Here, $h^2$ represents the heritability of the phenotypic trait of interest in the population. This model is commonly referred to as the Fisher infinitesimal model [Bulmer, 1980; Fisher, 1918].

In other words, the probability that a mother with a breeding value in class $j_f$ finds a male in class $j_m$ and gives birth to an individual with a breeding value in class $j$ is given by

$$G_{V_a} \left( x_j - \frac{x_{j_f} + x_{j_m}}{2} \right) \delta_{j_m}(\mathbb{R}\mathbf{n}), \qquad (17)$$

where $G_{V_a}$ represents the Gaussian distribution with a mean of 0 and a variance of $V_a$, and $\delta_{j_m}$ denotes the frequency of reproductive male individuals with a breeding value in class $j_m$ within the population. We make the assumption that both males and females experience similar selective pressures and have an even sex ratio. Therefore, the frequency of alive reproductive males is equivalent to the frequency of alive reproductive females, which can be defined as follows:

$$\delta_{j_m}(\mathbb{R}\mathbf{n}) = \frac{\sum_{i,k=1}^{s,g} F_{j_m,k} n_{i,j_m,k}}{\sum_{i,j,k=1}^{s,b,g} F_{j,k} n_{i,j,k}}. \qquad (18)$$

Then the transition matrices for a breeding value class $j$, $H_{i,k}[\mathbf{n}]$ is defined for all breeding value classes $j_f$ by

$$(H_{i,k}[\mathbf{n}])_{j,j_f} = \sum_{j_m=1}^{b} G_{V_a} \left( x_j - \frac{x_{j_f} + x_{j_m}}{2} \right) \delta_{j_m}(\mathbb{R}\mathbf{n}) \qquad (19)$$

3. Matrices for stage $i$ and breeding value class $j$, $M_{i,j}$ (dimensions $p \times p$), assigns new offspring, with breeding value $j$, to their phenotypic class. The matrices $M_{i,j}$ give the probability that a newborn with breeding value $x_j$ will express the phenotype $z_k$ at birth. Importantly, this probability is solely dependent on the breeding value of the newborn and is not influenced by the phenotype of the parent. For a newborn with breeding value class $j$, the probability to fall within the phenotypic class $k$ is

$$G_{V_E}(z_k - x_j),$$

where $G_{V_E}$ is the Gaussian distribution with mean 0 and variance $V_E = V_P - V_a = (1 - h^2)V_P$. Thus the matrix $M_{i,j}$, is defined by

$$M_{ij} = (G_{V_E}(z_k - x_j))_{k,l \in \{1, \dots, p\}} \qquad (20)$$
## Construction of block diagonal matrices

From the matrices $U_{j,k}$, $R_{j,k}$, $B_{i,k}$, $H_{i,k}$, $P_{i,j}$ and $M_{i,j}$, we construct block diagonal matrices; e.g.,

$$ U = \begin{pmatrix} \mathbf{U}_{11} & \mathbf{0} & \dots & \mathbf{0} \\ \mathbf{0} & \mathbf{U}_{21} & \dots & \mathbf{0} \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{0} & \mathbf{0} & \dots & \mathbf{U}_{bg} \end{pmatrix}. \qquad (21) $$

with similar block-diagonal constructions for $\mathbb{R}$, $\mathbb{B}$, $\mathbb{H}$, $\mathbb{P}$, and $\mathbb{M}$.

## Results

### Theoretical derivations of rates of adaptation

When selection acts through fertility $F$, that is, $F(z) = \exp(\log(F) + \beta z)$, the adaptation rates can be computed with respect to the parameters of our model. The selection gradient described by (2) takes the form

$$ \frac{\partial \ln(\bar{\lambda})}{\partial \bar{z}_t} = \frac{1}{\bar{\lambda}} \frac{\partial \bar{\lambda}}{\partial \bar{z}_t} = \frac{1}{\bar{\lambda}} \frac{\partial \bar{\lambda}}{\partial F} \frac{\partial F}{\partial \bar{z}_t} \qquad (22) $$

Using the sensitivity matrix $S$ associated with the matrix $(\tilde{\mathbf{U}} + \tilde{\mathbf{F}})$ and the generation time $T$, defined by (4) then

$$ \frac{\partial \bar{\lambda}}{\partial F} = \sum_{i_1, i_2=1}^{2} S_{i_1, i_2} \frac{\partial (\mathbf{R} + \mathbf{U})_{i_1, i_2}}{\partial F} = S_{1, 2} = \frac{v_1 w_2}{v^\top w} \qquad (23) $$

$$ \frac{\partial F}{\partial \bar{z}_t}(0) = \beta F \qquad (24) $$

Equations (3) and (5) for the rate of adaptation become

$$ \mathbf{R} \mathbf{A}_{\mathbf{T}} = h^2 \mathbf{V}_P \beta \frac{1}{T} \quad \text{and} \quad \mathbf{R} \mathbf{A}_{\mathbf{G}} = h^2 \mathbf{V}_P \beta. \qquad (25) $$

Then, assuming similar selective pressure ($\beta$), heritability ($h^2$), and phenotypic variance ($\mathbf{V}_P$) among the five species, the adaptation rates per generation are also all equal. In contrast, the rate of adaptation per unit of time will decrease with the species' generation time. When selection acts on fertility, the relative rate of change in fertility $V_\theta$ and the relative change in individual fitness $V_{R_0}$ are equal

$$ V_{\theta} = V_{R_0} = \beta. $$

The rates of adaptation measured in terms of unit of selection (6) and in terms of unit of individual selection (8) are equal

$$ \mathbf{R} \mathbf{A}_{\text{TS}\theta} = \mathbf{R} \mathbf{A}_{\text{TSR}_0} = h^2 \mathbf{V}_P \frac{1}{T} \quad \text{and} \quad \mathbf{R} \mathbf{A}_{\text{GS}\theta} = \mathbf{R} \mathbf{A}_{\text{GSR}_0} = h^2 \mathbf{V}_P. $$

In addition, the rate of adaptation per generation and per unit of individual selection only depends on the heritability and the phenotypic variance, when selection acts on the fertility. The theoretical formulas for fertility, along with those for juvenile survival ($S_J$), maturation rate ($\gamma$) and adult survival ($S_A$) are provided in Table 2, and their derivations are detailed in Appendix S4.
<table><caption>**Table 2:** Rates of adaptation for different measures of time and different measures of selection when selection acts on different vital rate of the life cycle</caption><thead><tr><th>Rate</th><th>Fertility</th><th>Juvenile survival</th><th>Maturation</th><th>Adult survival</th></tr></thead><tbody><tr><td>RA<sub>T</sub></td><td>h<sup>2</sup> V<sub>P</sub> β <sup>1</sup>&frasl;<sub>T</sub></td><td>h<sup>2</sup> V<sub>P</sub> β <sup>1</sup>&frasl;<sub>T</sub> <sup>λ(1 - S<sub>J</sub>)</sup>&frasl;<sub>λ - S<sub>J</sub>(1 - γ)</sub></td><td>h<sup>2</sup> V<sub>P</sub> β <sup>1</sup>&frasl;<sub>T</sub> <sup>(λ - S<sub>J</sub>)(1 - γ)</sup>&frasl;<sub>λ - S<sub>J</sub>(1 - γ)</sub></td><td>h<sup>2</sup> V<sub>P</sub> β <sup>1</sup>&frasl;<sub>T</sub> <sup>S<sub>A</sub>(1 - S<sub>A</sub>)</sup>&frasl;<sub>λ - S<sub>A</sub></sub></td></tr><tr><td>RA<sub>G</sub></td><td>h<sup>2</sup> V<sub>P</sub> β</td><td>h<sup>2</sup> V<sub>P</sub> β <sup>λ(1 - S<sub>J</sub>)</sup>&frasl;<sub>λ - S<sub>J</sub>(1 - γ)</sub></td><td>h<sup>2</sup> V<sub>P</sub> β <sup>(λ - S<sub>J</sub>)(1 - γ)</sup>&frasl;<sub>λ - S<sub>J</sub>(1 - γ)</sub></td><td>h<sup>2</sup> V<sub>P</sub> β <sup>S<sub>A</sub>(1 - S<sub>A</sub>)</sup>&frasl;<sub>λ - S<sub>A</sub></sub></td></tr><tr><td>RA<sub>TSθ</sub></td><td>h<sup>2</sup> V<sub>P</sub> <sup>1</sup>&frasl;<sub>T</sub></td><td>h<sup>2</sup> V<sub>P</sub> <sup>1</sup>&frasl;<sub>T</sub> <sup>λ</sup>&frasl;<sub>λ - S<sub>J</sub>(1 - γ)</sub></td><td>h<sup>2</sup> V<sub>P</sub> <sup>1</sup>&frasl;<sub>T</sub> <sup>(λ - S<sub>J</sub>)</sup>&frasl;<sub>λ - S<sub>J</sub>(1 - γ)</sub></td><td>h<sup>2</sup> V<sub>P</sub> <sup>1</sup>&frasl;<sub>T</sub> <sup>S<sub>A</sub></sup>&frasl;<sub>λ - S<sub>A</sub></sub></td></tr><tr><td>RA<sub>GSθ</sub></td><td>h<sup>2</sup> V<sub>P</sub></td><td>h<sup>2</sup> V<sub>P</sub> <sup>λ</sup>&frasl;<sub>λ - S<sub>J</sub>(1 - γ)</sub></td><td>h<sup>2</sup> V<sub>P</sub> <sup>(λ - S<sub>J</sub>)</sup>&frasl;<sub>λ - S<sub>J</sub>(1 - γ)</sub></td><td>h<sup>2</sup> V<sub>P</sub> <sup>S<sub>A</sub></sup>&frasl;<sub>λ - S<sub>A</sub></sub></td></tr><tr><td>RA<sub>TSR<sub>0</sub></sub></td><td>h<sup>2</sup> V<sub>P</sub> <sup>1</sup>&frasl;<sub>T</sub></td><td>h<sup>2</sup> V<sub>P</sub> <sup>1</sup>&frasl;<sub>T</sub> <sup>λ(1 - S<sub>J</sub>(1 - γ))</sup>&frasl;<sub>λ - S<sub>J</sub>(1 - γ)</sub></td><td>h<sup>2</sup> V<sub>P</sub> <sup>1</sup>&frasl;<sub>T</sub> <sup>(λ - S<sub>J</sub>)(1 - S<sub>J</sub>(1 - γ))</sup>&frasl;<sub>(1 - S<sub>J</sub>)(λ - S<sub>J</sub>(1 - γ))</sub></td><td>h<sup>2</sup> V<sub>P</sub> <sup>1</sup>&frasl;<sub>T</sub> <sup>(1 - S<sub>A</sub>)</sup>&frasl;<sub>λ - S<sub>A</sub></sub></td></tr><tr><td>RA<sub>GSR<sub>0</sub></sub></td><td>h<sup>2</sup> V<sub>P</sub></td><td>h<sup>2</sup> V<sub>P</sub> <sup>λ(1 - S<sub>J</sub>(1 - γ))</sup>&frasl;<sub>λ - S<sub>J</sub>(1 - γ)</sub></td><td>h<sup>2</sup> V<sub>P</sub> <sup>(λ - S<sub>J</sub>)(1 - S<sub>J</sub>(1 - γ))</sup>&frasl;<sub>(1 - S<sub>J</sub>)(λ - S<sub>J</sub>(1 - γ))</sub></td><td>h<sup>2</sup> V<sub>P</sub> <sup>(1 - S<sub>A</sub>)</sup>&frasl;<sub>λ - S<sub>A</sub></sub></td></tr></tbody></table>

## Method comparison

The adaptation rates per unit of time are almost identical in the three methods: theoretical derivations, IBM, or EvoDemo-Hyper MPM (Figure 3). Minor deviations between the theoretical derivations and MPM occur (for example, species 1 in the juvenile survival selection scenario, Figure 3), validating our theoretical approximations. Although IBM results embody demographic stochasticity and genetic drift, they match the deterministic results on average. However, the confidence interval around the mean change in breeding value after 100 years tends to be larger for scenarios with lower rates of adaptation (Figure 3) because of lower population sizes and higher genetic drift occurring in those scenarios.

MPM results align with IBM as they cover the observed phenotypic and breeding value ranges over time. When observed values reach class boundaries, adaptation rates decrease (see the Appendix S3).

Computation times differ significantly between methods (see Appendix S5). Theoretical and MPM computations are fast (0.04 to 49 seconds), independent of species or vital rate. In contrast, IBM simulation times vary widely with scenarios and initial population size (1.45 to 18,500 seconds). Multiple IBM simulation runs (e.g., 100 replicates) are needed to estimate adaptation, affecting dynamics and extinction probability, amplifying differences from deterministic simulations.

## Definition comparison

Determining which species have the greatest adaptation capacity and identifying vital rates with the highest adaptation potential under similar selection pressures are essential questions. Our theoretical findings show that responses vary based on the adaptation rate definition. Adaptation rates differ among species and depend on the selected vital rate, all influenced by how adaptation rate is defined.

The greatest differences in adaptation rate between species are found when adaptation is expressed per unit of time (RA<sub>T</sub>), compared to per unit of generation (RA<sub>G</sub>). The exception is when selection acts on adult survival (Figure 4). However, the pattern of variation in adaptation rates between species and vital rates remains consistent (top row on Figure 4), unless adaptation rates are calculated per selection units identified as the relative rate of change in an individual's fitness (RA<sub>TSθ</sub> and RA<sub>TSR<sub>0</sub></sub>, middle and bottom rows on Figure 4). For RA<sub>TSR<sub>0</sub></sub> and RA<sub>GSR<sub>0</sub></sub>,
<figure><img src="image_5.png" /><figcaption>**Figure 3:** Comparison of the breeding value changes after 100 years for the three approaches: theoretical derivations (THEO, blue bars), Individual-Based Model (IBM, orange bars), and EvoDemo-Hyper Matrix Projection Model (MPM, brown bars). Theoretical derivations $RA_T$ (THEO) are multiplied by 100 to match projection time of the IBM and MPM. Each panel corresponds to a particular vital rate under selection. Color bars illustrate the outcomes for the five species, ranging from the species with the shortest generation time to the species with the longest generation time. Error bars for IBMs represent the 95% confidence interval for the mean of 100 simulation replicates. Note the difference in y-axis scales among panels. Table S1 details the values.</figcaption></figure>

there are no longer variations in the adaptation rates between the vital rates. Furthermore, when adjusted for generation time, adaptation rates are the same across all vital rates and species.

**Comparison of the rate at which organisms adapt across different life histories.**

The rate of adaptation per time $RA_T$ is slower for species 5 than other species for all vital rates except adult survival. It decreases when the generation time increases if the pressure of selection acts on fecundity and juvenile survival. However, when the selection pressure acts on maturation rate and adult survival this monotonicity vanishes and intermediate species (2 and 3) have a more rapid rate of adaptation. Highest adaptation rates per time occur when selection acts on juvenile survival $S_J$, and remain high under fecundity selection. This pattern slightly differs for species 5, whose rate of adaptation is higher under fecundity selection (0.0013) than juvenile survival selection (0.0012). Furthermore, the main influence of juvenile survival is based on our assumption that the newborn survival rate $S_0$ equals the juvenile survival rate ($S_0 = S_J$). Under the opposite assumption, adaptation rate per time may be greater when selection influences fecundity (Appendix S2).

The rate of adaptation per generation $RA_G$ decreases along the slow-fast continuum if the pressure of selection acts on juvenile survival, but remains constant across species if it acts on fecundity. Furthermore, intermediate species (2 and 3) no longer have a faster rate of adaptation when selection pressure acts on the maturation rate or adult survival. Instead, the rate of adaptation per generation increases with the generation time when selection pressure acts on maturation and
<figure><img src="image_6.png" /><figcaption>**Figure 4:** Comparison of adaptation rates based on the six definitions (Equations (3)-(5)-(6)-(8)). Each panel presents a different definition. The rate of adaptation can be indicated per unit of time $RA_T$ (first column) or per unit of generation $RA_G$ (second column). These rates can also be indicated per unit of selection (second and third row) or not (first row). The second row is represented by the unit of selection expressed by the relative rate of change in the vital rate $\theta$: $RA_{TS\theta}$ left column and $RA_{GS\theta}$ right column. The third row corresponds to $RA_{TSR_0}$ (left column) and $RA_{GSR_0}$ (right column), where the selection units are defined by the relative rate of change in an individual's fitness, quantified by the mean lifetime reproductive success $R_0$. The color bar represents a specific vital rate that undergoes selection for the five species (x-axis) arranged from the species with the shortest generation time to the species with the longest generation time.</figcaption></figure>

adult survival. Similarly to the adaptation rate per time, the adaptation rate per generation exhibits the highest values under juvenile survival selection and remains high when selection acts on fecundity, with a switch for species 5 (Figure 4).

**The rates of adaptation per unit of vital rate** over time $RA_{TS\theta}$ exhibits a comparable trend to the rates of adaptation per time $RA_T$ when selection influences fecundity and the survival of juveniles. However, it increases with the generation time when selection acts on adult survival and it decreases when selection pressure impacts the maturation rate. The highest $RA_{TS\theta}$ values are still evident under juvenile survival selection in species 1 and 2. However, in species 3 to 5 the highest $RA_{TS\theta}$ values occur under selection on adult survival. Differences in vital rates and species are comparable when the rates of adaptation per unit of vital rate are expressed in terms of generation time, $RA_{GS\theta}$. However, adaptation rates exhibit significantly higher values when selection influences adult survival.
**The rate of adaptation per unit of fitness** over time $RA_{TSR_0}$ shows no differences among the vital rates. Additionally, when the data is normalized for generation time, the rate of adaptation $RA_{GSR_0}$ remains constant across both vital rates and species.

## Discussion

Here, we present an Evo-Demo Hyperstate Matrix Population Model (EvoDemo-Hyper MPM) that incorporates quantitative genetic processes to predict the evolution of breeding values and phenotypes within a population. We calculate theoretical adaptation rates using the Lande equation [Bulmer, 1980; Lande, 1976] for the rate of adaptation and demonstrate that our theoretical predictions, along with simulation results from the MPM, align with those from a more conventional evolutionary Individual-Based Model (IBM). This collaboration between evolutionary biologists and demographers highlighted different possible definitions of the rate of adaptation and their implications for assessing variations in the rate of adaptation across life-history traits and species' generation times. We found no overarching prediction as to which species or life-history trait yields a higher rate of adaptation; this is contingent on the life cycle and the specific vital rates subject to selection pressure, as well as the choice of definition.

## Converging Dynamics Across Diverse Methodologies

The debates surrounding the use of evolutionary IPMs in the years 2010 [Chevin, 2015; Janeiro et al., 2017; van Benthem et al., 2017] cast population matrix models against quantitative genetic models, and may have suggested that one or the other approach was superior. Our work illustrates the rather self-evident fact that different approaches lead to similar results when they actually represent the same biological processes. In particular, population matrix models can incorporate quantitative trait inheritance correctly and thus provide consistent results with other methods.

While it is relatively trivial to create models of evolution and demography using IBMs or iterative equations, empirical applications remain relatively rare. Based on our own experience, we believe that it is in part because the range of possible models and assumptions is very large for any given system, requiring intense research to develop a useful empirical model. Furthermore, IBMs can be particularly long to develop, check, and analyze. Using matrix projection model with a more constrained structure and ease of creation, e.g. Integral Projection Models (IPM), may be a good option to tackle evolution-demography questions.

Regrettably, the majority of IPMs that include evolutionary processes, subsequently referred to as EvoDemo-IPMs, have predominantly employed insufficient models for genetic transmission. For example, EvoDemo-IPMs incorrectly dampen evolutionary responses or they introduce an unjustified environmental influence on inheritance [for instance see inheritance function in Clark-Wolf et al., 2024].

An explicit evolutionary matrix model is not novel, but our Evo-Demo Hyperstate Matrix Population Model is a rare option to include quantitative genetics accounting for two-sex inheritance while keeping the flexibility of EvoDemo-IPMs. Although Coulson et al. [2017] introduced a broad framework for EvoDemo-IPMs with quantitative genetic inheritance, potentially with two sexes, they did not provide a code. To our knowledge, Simmonds et al. [2020] is the only application of the proper EvoDemo-IPMs Coulson et al. [2017], but it only incorporates the breeding values of the females within the inheritance function. Indeed, this model assumes that offspring breeding values are independent of paternal contributions, which, in the absence of stabilizing selection, may drive an increase in phenotypic variance over time. Furthermore, if environmental factors influence annual variations in the sex ratio of recruits by favoring males due to competitive advantages in resource-limited conditions, as suggested by Oddie [2000], excluding male contributions in the
model may lead to biases in population estimates.

For simple applications, the analytical results we derived may prove particularly useful. The the-
oretical formulas allow us to compute adaptation rates for all species, vital rates, and definitions
within a short CPU time of 0.04 s. This paves the way for further theoretical exploration across a
wider range of life histories, allowing for an examination of a more gradual spectrum of genera-
tion times.

Which species adapt faster? And what is adaptation anyway?

Lots of research currently aims to assess the potential of species to adapt to rapid anthropogenic environmental change [Radchuk et al., 2019; Urban et al., 2023]. More specifically, scientists and wildlife managers would benefit from knowledge on the potential speed of adaptation of species to given selection pressure. Unfortunately, the present work shows that there cannot be a one-line answer to the question of which type of species will adapt slower or faster.

First, differences in the rates of adaptation among species do not always fall along a continuum. For instance, when using $RA_T$, in the cases of maturation and adult survival, the highest rates are observed for species with intermediate generation times (Figure 4). The ordering of species from slowest to highest adaptation rate also varies depending on which vital rate is under selective pressure.

Furthermore, and perhaps more importantly, determining which species adapts faster depends on how we define the rate of adaptation. Throughout this work we highlighted two sensitive dimensions in the definition: the time unit and the measure of the strength of selection. Our six definitions provide dramatically different answers to which species, if any, adapts faster. Specifically, when defining adaptation using measures of selection on a linear scale ($RA_T$ and $RA_G$) or per unit of vital rate ($RA_{TS\theta}$ and $RA_{GS\theta}$), the variation in adaptation rates among species is contingent upon the particular vital rates targeted by selection. Thus determining which species adapt faster requires precise knowledge on selection pressure experienced by the species. The various perspectives on selection strengths and time units are not mutually exclusive; rather, they can be viewed as complementary strategies to take into account the effect of life history on adaptive dynamics. For example, an increase in mean survival rates will extend generation time, which consequently will slow down adaptation when measured per unit of vital rate and reduce selection pressure when measured per unit of fitness.

Our study does not explore another crucial aspect of variation in the definition of adaptation: whether it exclusively refers to genetic changes resulting from natural selection (as we assume in our work) or encompasses broader processes that enhance the expected population growth rate, such as phenotypic plasticity or shifts in demographic structure. In evolutionary biology, adaptation (referring to adaptive evolution rather than the state of a trait being adaptive) is typically defined as an evolutionary response to natural selection. This distinction is important because alternative definitions can create apparent paradoxes in fundamental theorems of micro-evolution [Kokko, 2021].

However, works on demography and wildlife management can adopt a broader definition of
adaptation, encompassing any deterministic increase in population growth rate, including im-
provements in phenotypic plasticity and changes in stage structure [e.g., Clark-Wolf et al., 2024;
Fox et al., 2019]. Phenotypic plasticity, defined as the ability of a genotype to produce different
phenotypes in response to environmental variation, can enable organisms to cope with environ-
mental changes in the short term. For example, plasticity in the timing of breeding of great tits
(*Parus major*) allows reducing mismatches between hatching and food availability under changing
conditions [Simmonds et al., 2020]. Our EvoDemo-Hyper model can easily be extended to incor-
porate phenotypic plasticity by defining non-trivial transition probabilities between phenotypic
classes. This can be achieved through the transition matrices $P_{i,j}$ (dimensions $p \times p$), which contain transitions between phenotypic classes for individuals of stage $i$ with breeding value in class $j$.

The intricacy of defining adaptation is further complicated by sentences like "adaptive change due to phenotypic plasticity," that raise ambiguity about whether phenotypic plasticity is a subset of adaptation or a distinct biological phenomenon. While we do not aim to prescribe one definition over another, we emphasize the importance of recognizing the diversity of definitions across disciplines and the potential for misunderstanding. Clarifying these distinctions is particularly relevant when designing studies or interpreting results across ecological and evolutionary research fields.

## Adapting to Adaptation: How Conservationists, Demographers, and Geneticists Define and Measure Change

In turn, different definitions may be better suited depending on the scientific context. For instance, a conservation manager focused on implementing actions to aid species' adaptation over time, might analyze adaptation rates per unit of time. They could observe that adaptation rates decline along the slow-fast continuum if selection pressure affects juvenile survival and fecundity. While the patterns become subtler when selection pressure impacts maturation and survival, there are no significant differences in adaptation rates overall. Consequently, conservation efforts might prioritize enhancing juvenile survival and fecundity in short-lived species and focus on any vital rates for longer-lived species to facilitate adaptation to global changes.

Demographers with an interest in life history strategies might interpret the rate of adaptation per time per unit of selection in terms of its effect on the relative rate of change in an individual's vital rate. From this perspective, adaptation rates appear faster when selection impacts juvenile survival and fecundity, as well as maturation in short-lived species, whereas in long-lived species, adaptation is quicker when selection acts on adult survival.

Quantitative geneticists, on the other hand, prioritize the study of the causes and consequences of genetic change. They express adaptation rate per generation time and the response of selective pressure in standardized units of relative fitness or relative fitness component. Their approach is in line with the theoretical basis and conceptual focus on genetic change and is encapsulated by the breeder's equation. In one formulation, the breeder's equation expresses the per generation predicted change in the trait under selection as the product of additive genetic variance and a selection gradient (i.e., the average slope of relative fitness on the trait). When standardizing the response by the strength of selection, the predicted response depends only on additive genetic variance in the trait, which explains why the rate of adaptation $RA_{GSR_0}$ is equal for all species and all vital rates on Figure 4. The actual standardized strength of selection does depend on the details of the life-history and which fitness component selection acts on, but by measuring selection in a standardized way, quantitative genetics can take those details out of the equation. Instead of the details of life-history and on which fitness component selection acts, quantitative geneticists therefore stress the importance of differences in generation time, strength of selection and additive genetic variance (or, equivalently, heritabilities) across species [e.g. Carlson et al., 2014; Kingsolver et al., 2012; Postma, 2014].

## Conclusion

Here, we demonstrate that both Evo-Demo Hyperstate Matrix Population Model and Individual-Based Models yield consistent results regarding the average dynamics of breeding values and phenotypic change. Thus, the choice of method is not crucial when analyzing the rate of evolution. Researchers may select a method based on computational efficiency, their background in
evolution or demography, or the degree of detail about the life cycle they wish to incorporate. We also highlight that patterns observed among species and vital rates are influenced by definitions, underscoring the importance of maintaining awareness of this context when drawing broad conclusions about species adaptation, particularly in response to global changes. Overall, we recommend researchers interested in which species adapt faster in a given scenario do not assume an answer based on some verbal expectation, but instead model the rate of adaptation with the definition that matters to them, and find out a quantitative answer.

## References

* Barfield, M., Holt, R.D. & Gomulkiewicz, R. (2011) Evolution in stage-structured populations. *The American Naturalist*, **177**, 397–409.

* Bienvenu, F. & Legendre, S. (2015) A new approach to the generation time in matrix population models. *The American Naturalist*, **185**, 834–843.

* Bulmer, M.G. (1980) *The Mathematical Theory of Quantitative Genetics*. Oxford, Clarendon Press.

* Carlson, S.M., Cunningham, C.J. & Westley, P.A.H. (2014) Evolutionary rescue in a changing world. *Trends in Ecology & Evolution*, **29**, 521–530.

* Caswell, H. (2000) *Matrix population models*, volume 1. Sinauer Sunderland, MA.

* Chevin, L.M. (2015) Evolution of adult size depends on genetic variance in growth trajectories: A comment on analyses of evolutionary dynamics using integral projection models. *Methods in Ecology and Evolution*, **6**, 981–986.

* Childs, D.Z., Sheldon, B.C. & Rees, M. (2016) The evolution of labile traits in sex- and age-structured populations. *Journal of Animal Ecology*, **85**, 329–342.

* Clark-Wolf, T.J., Boersma, P.D., Plard, F., Rebstock, G.A. & Abrahms, B. (2024) Increasing environmental variability inhibits evolutionary rescue in a long-lived vertebrate. *Proceedings of the National Academy of Sciences*, **121**, e2406314121.

* Coste, C.F. & Pavard, S. (2020) Analysis of a multitrait population projection matrix reveals the evolutionary and demographic effects of a life history trade-off. *Ecological Modelling*, **418**, 108915.

* Coulson, T., Kendall, B.E., Barthold, J., Plard, F., Schindler, S., Ozgul, A. & Gaillard, J.m. (2017) Modeling Adaptive and Nonadaptive Responses of Populations to Environmental Change. *The American Naturalist*, **190**.

* de Villemereuil, P., Schielzeth, H., Nakagawa, S. & Morrissey, M. (2016) General methods for evolutionary quantitative genetic inference from generalized mixed models. *Genetics*, **204**, 1281–1294.

* Easterling, M.R., Ellner, S.P. & Dixon, P.M. (2000) Size-specific sensitivity: applying a new structured population model. *Ecology*, **81**, 694–708.

* Falconer, D.S. & Mackay, T.F.C. (1996) *Introduction to quantitative genetics*, chapter Variance. Pearson, 4th edition.

* Fisher, R.A. (1918) The correlation between relatives on the supposition of Mendelian inheritance. *Trans R Soc Edinburgh*, **52**, 399–433.

* Fisher, R.A. (1930) *The genetical theory of natural selection*. Oxford University Press.

* Fox, R.J., Donelson, J.M., Schunter, C., Ravasi, T. & Gaitán-Espitia, J.D. (2019) Beyond buying time: the role of plasticity in phenotypic adaptation to rapid environmental change. *Philosophical Transactions of the Royal Society B: Biological Sciences*, **374**, 20180174.

* Hunter, C.M. & Caswell, H. (2005) The use of the vec-permutation matrix in spatial matrix population models. *Ecological Modelling*, **188**, 15–21.
Janeiro, M.J., Coltman, D.W., Festa-Bianchet, M., Pelletier, F. & Morrissey, M.B. (2017) Towards robust evolutionary inference with integral projection models. *Journal of Evolutionary Biology*, **30**, 270–288.

Jenouvrier, S., Long, M.C., Coste, C.F.D., Holland, M., Gamelon, M., Yoccoz, N.G. & Sæther, B.E. (2022) Detecting climate signals in populations across life histories. *Global Change Biology*, **28**, 2236–2258.

Kendall, B.E., Fujiwara, M., Diaz-Lopez, J., Schneider, S., Voigt, J. & Wiesner, S. (2019) Persistent problems in the construction of matrix population models. *Ecological modelling*, **406**, 33–43.

Kingsolver, J.G., Diamond, S.E., Siepielski, A.M. & Carlson, S.M. (2012) Synthetic analyses of phenotypic selection in natural populations: lessons, limitations and future directions. *Evolutionary Ecology*, **26**, 1101–1118.

Kokko, H. (2021) The stagnation paradox: the ever-improving but (more or less) stationary population fitness. *Proceedings of the Royal Society B: Biological Sciences*, **288**, 20212145.

Lande, R. (1979) Quantitative genetic analysis of multivariate evolution, applied to brain: body size allometry. *Evol*, **33**, 402–416.

Lande, R. (1976) Natural selection and random genetic drift in phenotypic evolution. *Evolution*, **30**, 314–334.

Lynch, M. & Walsh, B. (1998) *Genetics and Analysis of Quantitative Traits*. Sinauer Associates.

Neubert, M.G. & Caswell, H. (2000) Density-dependent vital rates and their population dynamic consequences. *Journal of Mathematical Biology*, **41**, 103–121.

Oddie, K.R. (2000) Size matters: competition between male and female great tit offspring. *Journal of Animal Ecology*, **69**, 903–912.

Pelletier, F. (2019) Testing evolutionary predictions in wild mice. *Science*, **363**, 452–453.

Postma, E. (2014) Four decades of estimating heritabilities in wild vertebrate populations: improved methods, more data, better estimates? A. Charmentier, D. Garant & L.E.B. Kruuk, eds., *Quantitative genetics in the wild*, pp. 16–33. Oxford University Press, Oxford, 1st edition.

Queller, D.C. (2017) Fundamental Theorems of Evolution. *The American Naturalist*, **189**, 345–353.

Radchuk, V., Reed, T., Teplitsky, C., Pol, M.v.d., Charmantier, A., Hassall, C., Adamík, P., Adri-aensen, F., Ahola, M.P., Arcese, P., Avilés, J.M., Balbontin, J., Berg, K.S., Borras, A., Burthe, S., Clobert, J., Dehnhard, N., Lope, F.d., Dhondt, A.A., Dingemanse, N.J., Doi, H., Eeva, T., Fickel, J., Filella, I., Fossøy, F., Goodenough, A.E., Hall, S.J.G., Hansson, B., Harris, M., Hasselquist, D., Hickler, T., Joshi, J., Kharouba, H., Martínez, J.G., Mihoub, J.B., Mills, J.A., Molina-Morales, M., Moksnes, A., Ozgul, A., Parejo, D., Pilard, P., Poisbleau, M., Rousset, F., Rödel, M.O., Scott, D., Senar, J.C., Stefanescu, C., Stokke, B.G., Kusano, T., Tarka, M., Tarwater, C.E., Thonicke, K., Thor-ley, J., Wilting, A., Tryjanowski, P., Merilä, J., Sheldon, B.C., Møller, A.P., Matthysen, E., Janzen, F., Dobson, F.S., Visser, M.E., Beissinger, S.R., Courtiol, A. & Kramer-Schadt, S. (2019) Adaptive responses of animals to climate change are most likely insufficient. *Nature Communications*, **10**, 1–14.

Rees, M. & Ellner, S.P. (2019) Why So Variable: Can Genetic Variance in Flowering Thresholds Be Maintained by Fluctuating Selection? *The American Naturalist*, **194**, E13–E29.

Roth, G. & Caswell, H. (2016) Hyperstate matrix models: extending demographic state spaces to higher dimensions. *Methods in Ecology and Evolution*, **7**, 1438–1450.

Salguero-Gómez, R., Jones, O.R., Jongejans, E., Blomberg, S.P., Hodgson, D.J., Mbeau-Ache, C., Zuidema, P.A., de Kroon, H. & Buckley, Y.M. (2016) Fast-slow continuum and reproductive strategies structure plant life-history variation worldwide. *Proceedings of the National Academy of
*Sciences*, **113**, 230–235.

Simmonds, E.G., Cole, E.F., Sheldon, B.C. & Coulson, T. (2020) Testing the effect of quantitative genetic inheritance in structured models on projections of population dynamics. *Oikos*, **129**, 559–571.

Urban, M.C., Swaegers, J., Stoks, R., Snook, R.R., Otto, S.P., Noble, D.W.A., Moiron, M., Hällfors, M.H., Gómez-Llano, M., Fior, S., Cote, J., Charmantier, A., Bestion, E., Berger, D., Baur, J., Alexander, J.M., Saastamoinen, M., Edelsparre, A.H. & Teplitsky, C. (2023) When and how can we predict adaptive responses to climate change? *Evolution Letters*, p. qrad038.

van Benthem, K.J.V., Bruijning, M., Bonnet, T., Jongejans, E., Postma, E. & Ozgul, A. (2017) Disentangling evolutionary, plastic and demographic processes underlying trait dynamics: A review of four frameworks. *Methods in Ecology and Evolution*, **8**, 75–85.

Vedder, O., Bouwhuis, S. & Sheldon, B.C. (2013) Quantitative Assessment of the Importance of Phenotypic Plasticity in Adaptation to Climate Change in Wild Bird Populations. *PLOS Biology*, **11**, e1001605.
# Supplementary information

## S1 Individual-based Model.

### Individual-based simulations

Individual-based simulations are based on the same life cycle (Figure 1), incorporating the same assumptions as the infinitesimal quantitative genetic model described for sexual reproduction in the EvoDemo-Hyper model and using vital rates specified in Table 1. Unlike the EvoDemo-Hyper model which tracks density of classes of phenotype-genotype-stage, the individual-based model (IBM) explicitly represents individuals at each time step. Individual-based simulations, more commonly used in evolutionary research, validate the results obtained from the EvoDemo-Hyper model and the analytical approximations. In addition, in the simulations survival, maturation, and reproduction are probabilistic events and thus subject to chance, introducing demographic stochasticity into the system (Figure S1). The EvoDemo-Hyper model does not take into account demographic stochasticity, although it has minimal impact on large populations.

The initial population consists of a combination of adult and juvenile individuals, whose breeding values are randomly drawn from a normal distribution with a mean of 0 and a variance of $V_a$, while the environmental values are drawn from a normal distribution with a mean of 0 and a variance of $V_P - V_a$.

At the beginning of each time step, new young individuals are born. The reproductive success of mature females at a given time step follows a Poisson distribution and is dependent on their phenotypes (see section "Phenotypic variation, heritability, and selective pressure on vital rates"). In contrast, the fathers of the offspring are randomly assigned from the pool of mature males at that time. Subsequently, adults and juveniles, including those born during the current time step, face a risk of mortality depending on their phenotype. Lastly, juveniles born before the current time step have the opportunity to mature into adults with a given probability dependent on their phenotype. If there are no females or no males present before the start of the next time step, the population is considered extinct, and the simulation ends. Otherwise, the life cycle is repeated for 100 time steps, i.e., 100 years.

### Comparison between methods

Using each of the three methods (theoretical derivations, the EvoDemo-Hyper model and individual-based simulations), we calculated the rate of adaptation according to the six definitions, and according to the four different pathways (i.e., vital rate) through which selection can act for each species separately (Table S1). For the EvoDemo-Hyper model and IBM, we estimate the rate of adaptation per time ($RA_T$) as the difference in mean breeding value in the population between time 101 and time 1, divided by 100. As we did not model plastic changes, the change in mean breeding values is very close to the change in mean phenotype.
<table><caption>**Table S1:** Difference in breeding value after 100 years of simulation according to the demographic pathway (i.e., vital rate) through which selective pressure was applied for each of the five fictive species considered. Results from the matrix population model (MPM), the theoretical analytical equations (THEO), and the individual-based model (IBM) are compared. IBM values are averages (standard errors) based on 100 simulation replicates, excluding replicates in which populations went extinct.</caption><thead><tr><th rowspan="2">Demographic pathway</th><th rowspan="2">Model</th><th colspan="5">Species</th></tr><tr><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td rowspan="3">Fertility f</td><td>MPM</td><td>1.460</td><td>1.280</td><td>0.507</td><td>0.196</td><td>0.144</td></tr><tr><td>THEO</td><td>1.465</td><td>1.266</td><td>0.476</td><td>0.190</td><td>0.127</td></tr><tr><td>IBM</td><td>1.550<br>(0.040)</td><td>1.280<br>(0.028)</td><td>0.514<br>(0.019)</td><td>0.267<br>(0.022)</td><td>0.162<br>(0.025)</td></tr><tr><td rowspan="3">Adult survival S<sub>A</sub></td><td>MPM</td><td>0.053</td><td>0.255</td><td>0.362</td><td>0.168</td><td>0.113</td></tr><tr><td>THEO</td><td>0.053</td><td>0.253</td><td>0.381</td><td>0.176</td><td>0.118</td></tr><tr><td>IBM</td><td>0.039<br>(0.081)</td><td>0.235<br>(0.057)</td><td>0.317<br>(0.030)</td><td>0.150<br>(0.023)</td><td>0.068<br>(0.023)</td></tr><tr><td rowspan="3">Juvenile survival S<sub>J</sub></td><td>MPM</td><td>2.422</td><td>1.874</td><td>0.681</td><td>0.246</td><td>0.128</td></tr><tr><td>THEO</td><td>2.660</td><td>2.013</td><td>0.674</td><td>0.240</td><td>0.124</td></tr><tr><td>IBM</td><td>2.740<br>(0.029)</td><td>2.070<br>(0.010)</td><td>0.770<br>(0.006)</td><td>0.415<br>(0.004)</td><td>0.139<br>(0.007)</td></tr><tr><td rowspan="3">Maturation rate γ</td><td>MPM</td><td>0.134</td><td>0.455</td><td>0.229</td><td>0.102</td><td>0.092</td></tr><tr><td>THEO</td><td>0.131</td><td>0.436</td><td>0.221</td><td>0.100</td><td>0.088</td></tr><tr><td>IBM</td><td>0.163<br>(0.099)</td><td>0.504<br>(0.050)</td><td>0.240<br>(0.033)</td><td>0.117<br>(0.023)</td><td>0.149<br>(0.025)</td></tr></tbody></table>
<figure><img src="image_7.png" /><figcaption>**Figure S1:** Results from the individual-based simulations for the changes in breeding values over a 100 years period after directional positive selection has been applied to a) fertility, b) adult survival, c) juvenile survival and d) maturation rate for five species with contrasting life cycles. Thick plain lines correspond to the average change over 100 replicates, and the dashed lines represent 80% intervals across those replicates. For the fertility selection case, species 1, 2 and 3 follow the fixed value model of selection, whereas species 4 and 5 follow the Bernoulli model of selection. Populations are initialised with 200 individuals at demographic equilibrium.</figcaption></figure>
## S2 Comparison of life cycle

Matrix Population Models (MPMs) are based on the assumption that the population census takes place either immediately before reproduction (referred to as a "prebreeding census") or right after reproduction (referred to as a "postbreeding census") [Caswell, 2000; Kendall et al., 2019]. During a prebreeding census, the youngest age group consists of individuals who are all approximately one year old, categorized as age one. Conversely, in a postbreeding census, the initial age group includes newborn individuals, all aged zero. In our model, age is not directly included, and we assume that individuals are classified as juveniles or adults at the start of the time step in the prebreeding census and at the end of the time step in a postbreeding census. As a result, juveniles would be considered one year old and older in both scenarios.

In the model, the production of new individuals in year $t + 1$ by a mature individual (adult) at time $t$ is denoted as the fertility rate $F$. The birth rate or fecundity $f$ represents the mean number of offspring produced at time $t$ by a mature individual alive at that time. This process is viewed as nearly instantaneous. The fertility rate $F$ is derived by multiplying $f$ by a survival factor.

In a prebreeding census model, the fertility coefficient is represented as $F = fS_0$: where an adult produces $f$ offspring right after the census. The offspring then survive until the end of the time-step with a survival rate of $S_0$. In a post-breeding census, adults must survive in order to reproduce, resulting in a birth rate at the end of the time interval. Therefore, the fertility coefficient can be calculated as $F = S_A f$.

In the paper, we demonstrate our findings using a prebreeding life cycle, where $S_0 = S_J$, assuming that the selective pressures affecting survival during the first year of life ($S_0$) are equivalent to those affecting the survival of juvenile individuals ($S_J$). We adopted this simplified hypothesis to delineate the selective pressures acting on newborns via the birth rate $f$ of those who influence their juvenile survival $S_0 = S_J$. In this context, the rate of adaptation when selective pressure influences $S_J$ is higher than if $S_0 \neq S_J$ because when selection acts on $S_J$, it affects both survival and fertility $F$ (Figure S2).

Distinguishing between selective pressures on birth rate and survival becomes even more critical in a postbreeding life cycle, where it is essential to differentiate between the birth rate and adult survival since the selective pressures on newborns and their parents diverge. Demographers frequently concentrate on fertility in theoretical research (e.g., Jenouvrier et al. [2022]; Neubert & Caswell [2000]), and in this appendix, we investigate a range of scenarios based on varying assumptions to assess the impact of choosing between prebreeding versus postbreeding life cycles, as well as the choice of selection acting on fertility versus the birth rate.

Specifically, we show the rate of adaptation over time when selection influences the four vital rates in five scenarios:

1. In the prebreeding census model, selection acts on fertility (referred as $F = fS_0$ on x-axis of Figure S2);

2. In the postbreeding census model, selection acts on fertility (referred as $F = fS_A$ on x-axis of Figure S2);

3. In the prebreeding census model, selection acts on birth rate $f$ but not $S_0$ (referred as $f$ on x-axis of Figure S2);

4. In the postbreeding census model, selection acts either on $f$ or $S_A$ referred as $f \& S_A$ on x-axis of Figure S2;

5. In the prebreeding census model, selection acts either on $f$ or $S_0$ with $S_0 = S_J$ referred as $f \& S_J$ on x-axis of Figure S2.
This fifth scenario is the one illustrated in the main manuscript.

Our findings remain consistent across different scenarios when selection influences reproduction, either through birth rate or fertility (Figure S2 upper left panel). In our model, selection affects a single vital rate, either *F* or *f*, at a time. Therefore, multiplying the birth rate by a constant rate does not alter the outcomes.

As expected, the rate of adaptation over time does not vary between different scenarios when selection influences maturation, as this parameter is not influenced by our various assumptions (Figure S2 lower left panel).

The rate of adaptation over time is higher under our assumption that $S_0 = S_J$ when selection acts on both parameters in the same way (Figure S2 upper right panel). In fact, in this case, the sensitivity of the population growth rate to juvenile survival is greater (Figure S3 upper right panel). However, the patterns of variations of the rate of adaptation over time remain the same between species.

The key distinctions emerge in the prebreeding and postbreeding life cycles, with varying patterns of adaptation rates over time among species when selection influences adult survival (Figure S2 lower right panel). The influence of adult survival on the population growth rate varies significantly between the prebreeding and postbreeding stages, which affects the rate of adaptation over time (Figure S3 lower right panel). These differences in the rate of adaptation are particularly notable for short-lived species, with values of 0.05 for prebreeding compared to 1.44 for postbreeding for species # 1, and 0.25 compared to 1.23 for species 2, while they remain relatively consistent for other species. Consequently, the adaptation rate decreases with generation time when selection acts on adult survival in a postbreeding life cycle, whereas it peaks for species with intermediate generation times (species # 3) in a prebreeding life cycle.
<figure><img src="image_8.png" /><figcaption>**Figure S2:** Comparison of the adaptation rate per unit of time $RA_T$ among the five life cycle scenario. Each panel displays a comparison for a particular vital rate that is subject to selection. The color bar illustrates the outcomes for the five species, ranging from the species with the shortest generation time (blue representing fast species) to the species with the longest generation time (green representing slow species). Specifically, our five scenarios on the x-axis are as follows:</figcaption></figure>

1. In the prebreeding census model, selection acts on fertility (referred as $F = fS_0$);

2. In the postbreeding census model, selection acts on fertility (referred as $F = fS_A$);

3. In the prebreeding census model, selection acts on birth rate $f$ but not $S_0$ (referred as $f$);

4. In the postbreeding census model, selection acts either on $f$ or $S_A$ referred as $f\&S_0$);

5. In the prebreeding census model, selection acts either on $f$ or $S_0$ with $S_0 = S_J$ referred as $fS_J$ (scenario of the main manuscript).
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

**Figure S3:** Comparison of the sensitivity of the population growth rate with respect to vital rates among the five life cycle scenarios. Each panel displays a comparison for a particular vital rate that is subject to selection. The color bar illustrates the outcomes for the five species, ranging from the species with the shortest generation time (blue represents fast species) to the species with the longest generation time (green represents slow species).
## S3 Comparison of additive variance and functional relationships of the phenotypes.

In the manuscript, we present our findings using specific values for additive genetic variance ($V_a = 0.2$) and the slope of the functional relationship between phenotype and vital rates ($\beta = 0.15$). This appendix extends our analysis by exploring a variety of values for these two parameters. Additionally, we offer recommendations on choosing the suitable range of phenotypic categories to be incorporated into matrix models and the optimal time period for calculating the adaptation rate.

### Functional relationship

When selection acts on juvenile survival $S_J$, adult survival $S_A$ or maturation $\gamma$, the functional relationships are:

$$S_J(z) = \text{logit}^{-1}(\text{logit}(S_J) + \beta z)$$

$$S_A(z) = \text{logit}^{-1}(\text{logit}(S_A) + \beta z)$$

$$\gamma(z) = \text{logit}^{-1}(\text{logit}(\gamma) + \beta z)$$

When selection acts on fecundity $f$ the functional relationships are:

$$f(z) = \exp(\log(f) + \beta z).$$

The resulting functional relationships are shown in Figure S4 for $z = [-4, 4]$, $V_a = 0.2$, $\beta = 0.15$.

<figure><img src="image_13.png" /><figcaption>**Figure S4:** Functional relationships between the phenotype and vital rates for $\beta = 0.15$. Each panel displays a comparison for a particular vital rate. The colors illustrate the relationships for the five species, ranging from the species with the shortest generation time (blue represents fast species) to the species with the longest generation time (green represents slow species)</figcaption></figure>
**Phenotypic range and the optimal time period for calculating the adaptation rate.**

It is critical to note that because MPMs are structured by phenotypic class, adaptation is inher-
ently limited by the last phenotypic class. In our analysis, since selection is directional, both the
breeding values and phenotype inevitably reach a plateau at a maximum value (panels located
in the upper left and lower right of the Figures S5-S10). This is demonstrated by the shift in the
distribution of breeding values over time, as shown in the middle top panel of Figures S5-S10.
The initial step involves choosing a suitable range of phenotypic variations (z) to prevent unre-
alistic values of vital rates (panel located at the bottom left of the Figures S5-S10). This task can
be accomplished by visualizing the functional relationships (lower left panel of Figures S5-S10),
and the average vital rates over time (the lower middle panel of Figures S5-S10). For instance, for
the third species when selection affects adult survival, the average adult survival rate at time zero
is approximately 0.8 (Figure S7). The adult survival rate varies between 0.36 and 0.96 for differ-
ent phenotypes in the population. With intense and rapid selection (V_a = 0.8 and β = 0.5), the
average vital rate rises to 0.92 within a century.

The next stage involves accurately estimating the rate of adaptation $RA_T$ during the adapt-
tion period. Indeed, $RA_T$ decreases as the optimal phenotype approaches (panels located on the
upper left and lower right of the Figures S5-S7). Furthermore, additive genetic variance and phe-
notypic variance decrease as the optimal phenotype is closer (panel located on the top left of the
Figures S5-S7). When adaptation reaches a plateau, it indicates that the population has reached
an adaptive peak, where further gains in fitness are constrained by the selective pressures defined
by the functional relationships and the range of phenotypic classes in our model. For instance, in
three cases illustrated in Figures S5 to S7, the values $RA_T$ calculated over a 100-year period are
0.0328, 0.0271, and 0.0238, respectively. When calculated over the initial ten years, the values of
$RA_T$ are 0.0548, 0.0794, and 0.0427, respectively. In our main manuscript, we report adaptation
rates calculated over the primary adaptation period, excluding the period of slowdown near the
plateau, to capture how quickly a population's average fitness increases over time in response to
these selective pressures.

In this appendix, we have documented the rate of adaptation calculated over the initial ten years, because adaptation can proceed rapidly in the presence of high additive genetic variance $V_a$ and strong selective pressure $\beta$. The rate of adaptation calculated over a century is reported in the main body of the paper. We present identical comprehensive outcomes for $\beta$ = 0.15 and $V_a$ = 0.2 as discussed in the main text, focusing on the two highest $RA_T$ values seen in fast species one and the minimum $RA_T$ recorded for slow species four. Although the highest value $RA_T$ estimated for species one when selection acts on juvenile survival may be underestimated due to the slowing of selection at the end of the century, the entire range of vital rates and species studied indicates that achieving the optimal phenotype was far from being realized.

**Comparison of the rate of adaptation for different levels of additive genetic variance and the slope of the relationship between phenotype and vital rates.**

As expected, the rate of adaptation per unit of time, $RA_T$, increases when both $\beta$ and $V_a$ increase (Figure S11). The pattern of variations of $RA_T$ between species and vital rates is the same regardless of the values of additive genetic variance ($V_a$) and the slope of the functional relationship between phenotype and vital rates ($\beta$). The pattern of variations of $RA_T$ is described in detail in the main body of the article for $\beta$ = 0.15 and $V_a$ = 0.2.
<figure><img src="image_14.png" /><figcaption>**Figure S5:** Detailed eco-evolutionary outcomes of MPMs output for species one under strong selection on fecundity ($V_a = 0.2$ and $\beta = 0.5$). The top left panel illustrates the average breeding value ± one standard deviation across time. The top center panel displays the distribution of breeding values at different time points. The top right panel exhibits the additive genetic variance and phenotypic variance over time. The bottom right panel illustrates the functional relationship between the phenotype and vital rate. The middle bottom panel depicts the average vital rate over time. The bottom right panel shows the average phenotype ± one standard deviation over time.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure S6:** Detailed eco-evolutionary outcomes of MPMs output for species two under rapid selection on juvenile survival ($V_a = 0.8$ and $\beta = 0.15$). Same legends as previous figure.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure S7:** Detailed eco-evolutionary outcomes of MPMs output for species three under strong and rapid selection on adult survival ($V_a = 0.8$ and $\beta = 0.5$). Same legends as previous figure.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure S8:** Detailed eco-evolutionary outcomes of MPMs output for fast species one when subjected to selection based on juvenile survival ($V_a = 0.2$ and $\beta = 0.15$). Same legends as previous figure.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure S9:** Detailed eco-evolutionary outcomes of MPMs output for fast species one when subjected to selection based on fecundity ($V_a = 0.2$ and $\beta = 0.15$). Same legends as previous figure.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure S10:** Detailed eco-evolutionary outcomes of MPMs output for slow species four when subjected to selection based on adult survival ($V_a = 0.2$ and $\beta = 0.15$). Same legends as previous figure.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure S11:** Comparison of the rate of adaptation for various values for the additive genetic variance ($V_a$) and the slope of the functional relationship between phenotype and vital rates ($\beta$). Each panel presents a different combination of these two parameters. The rate of adaptation are indicated per unit of time $RA_T$. The color bar represents a specific vital rate that undergoes selection for the five species (x-axis) arranged from the species with the shortest generation time to the species with the longest generation time.</figcaption></figure>
# S4 Theoretical adaptation rates of the different vital rates

In the subsequent sections, we establish the rates of adaptation for each vital rate of the life cycle described in Figure 1.

## Juvenile survival $S_J$ (fertility case)

When selection acts on juvenile survival $S_J$, the selection gradient becomes

$$ \frac{\partial \ln(\lambda)}{\partial \bar{z}_t} = \frac{1}{\lambda} \frac{\partial \lambda}{\partial \bar{z}_t} = \frac{1}{\lambda} \frac{\partial \lambda}{\partial S_J} \frac{\partial S_J}{\partial \bar{z}_t} \quad (26) $$

Since the phenotypic trait acts on the logistic scale on the juvenile survival, i.e. $S_J(z) = \text{invlogit}(\text{logit}(S_J) + \beta z)$, the selection gradient can be computed using the sensitivity matrix associated to the matrix $(\mathbf{R} + \mathbf{U})$

$$ \frac{\partial \lambda}{\partial S_J} = \sum_{i_1, i_2=1}^{2} S_{i_1, i_2} \frac{\partial(\mathbf{R} + \mathbf{U})_{i_1, i_2}}{\partial S_J} = S_{1,1}(1-\gamma) + S_{2,1}\gamma = \frac{\lambda}{S_J} \frac{v_1 w_1}{v^\top w} \quad (27) $$

$$ \frac{\partial S_J}{\partial \bar{z}_t}(0) = \beta(1 - S_J)S_J \quad (28) $$

Thus, the rates of adaptation per unit of time and per generation satisfy

$$ \mathbf{RA}_{\mathbf{T}} = h^2 \mathbf{V}_P \beta \frac{1}{T} \frac{\lambda(1 - S_J)}{\lambda - S_J(1 - \gamma)} \quad \text{and} \quad \mathbf{RA}_{\mathbf{G}} = h^2 \mathbf{V}_P \beta \frac{\lambda(1 - S_J)}{\lambda - S_J(1 - \gamma)} \quad (29) $$

The relative rate of change in juvenile survival is given by $V_\theta = \beta(1 - S_J)$ and the rates of adaptation per unit of selection satisfy

$$ \mathbf{RA}_{\mathbf{TS}\theta} = h^2 \mathbf{V}_P \frac{1}{T} \frac{\lambda}{\lambda - S_J(1 - \gamma)} \quad \text{and} \quad \mathbf{RA}_{\mathbf{GS}\theta} = h^2 \mathbf{V}_P \frac{\lambda}{\lambda - S_J(1 - \gamma)} \quad (30) $$

The strength of individual selection $V_{R_0}$ is

$$ V_{R_0} = \beta(1 - S_J) \frac{1}{1 - S_J(1 - \gamma)} \quad (31) $$

and the rates of adaptation per units of individual selection are

$$ \mathbf{RA}_{\mathbf{TSR}_0} = h^2 \mathbf{V}_P \frac{1}{T} \frac{\lambda(1 - S_J(1 - \gamma))}{\lambda - S_J(1 - \gamma)} \quad \text{and} \quad \mathbf{RA}_{\mathbf{GSR}_0} = h^2 \mathbf{V}_P \frac{\lambda(1 - S_J(1 - \gamma))}{\lambda - S_J(1 - \gamma)} \quad (32) $$

If the growth rate $\lambda$ is equal to 1, then the rate of adaptation per selection units per generation is independent of the vital rates.

## Maturation $\gamma$

When selection acts on maturation at the logistic scale, $\gamma(z) = \text{invlogit}(\text{logit}(\gamma) + \beta z)$, the selection gradient is given by

$$ \frac{\partial \lambda}{\partial \gamma} = \sum_{i_1, i_2=1}^{2} S_{i_1, i_2} \frac{\partial(\mathbf{R} + \mathbf{U})_{i_1, i_2}}{\partial \gamma} = -S_{1,1}S_J + S_{2,1}S_J = \frac{\lambda}{T\gamma} \frac{\lambda - S_J}{\lambda - S_J(1 - \gamma)} \quad (33) $$

$$ \frac{\partial S_J}{\partial \bar{z}_t}(0) = \beta(1 - \gamma)\gamma \quad (34) $$

Thus the rates of adaptation per unit of time and per generation are:

$$ \mathbf{RA}_{\mathbf{T}} = h^2 \mathbf{V}_P \beta \frac{1}{T} \frac{(\lambda - S_J)(1 - \gamma)}{\lambda - S_J(1 - \gamma)} \quad \text{and} \quad \mathbf{RA}_{\mathbf{G}} = h^2 \mathbf{V}_P \beta \frac{(\lambda - S_J)(1 - \gamma)}{\lambda - S_J(1 - \gamma)} \quad (35) $$
The strength of selection is given by $V_{\theta} = \beta(1 - \gamma)$ and the rates of adaptation per selection unit are

$$ \mathbf{RA}_{\text{TS}\theta} = h^2 \mathbf{V}_P \frac{1}{T} \frac{(\lambda - S_J)}{\lambda - S_J(1 - \gamma)} \quad \text{and} \quad \mathbf{RA}_{\text{GS}\theta} = h^2 \mathbf{V}_P \frac{(\lambda - S_J)}{\lambda - S_J(1 - \gamma)} \qquad (36) $$

The strength of individual selection $V_{R_0}$ is

$$ V_{R_0} = \beta \frac{(1 - S_J)(1 - \gamma)}{1 - S_J(1 - \gamma)} \qquad (37) $$

and the rates of adaptation per units of individual selection are

$$ \mathbf{RA}_{\text{TSR}_0} = h^2 \mathbf{V}_P \frac{1}{T} \frac{(\lambda - S_J)(1 - S_J(1 - \gamma))}{(1 - S_J)(\lambda - S_J(1 - \gamma))} \quad \text{and} \quad \mathbf{RA}_{\text{GSR}_0} = h^2 \mathbf{V}_P \frac{(\lambda - S_J)(1 - S_J(1 - \gamma))}{(1 - S_J)(\lambda - S_J(1 - \gamma))} \qquad (38) $$

If the growth rate $\lambda$ equals 1, then the adaptation rate per selection units per generation does not depend on the vital rates.

## Adult survival $S_A$

When selection affects adult survival $S_A$ at the logistic scale, invlogit(logit($S_A$) + $\beta z$), the selection gradient can be expressed as

$$ \frac{\partial \lambda}{\partial S_A} = \sum_{i_1, i_2=1}^{2} S_{i_1, i_2} \frac{\partial(\mathbf{R} + \mathbf{U})_{i_1, i_2}}{\partial S_A} = S_{2,2} = \frac{\lambda}{T} \frac{1}{\lambda - S_A} \qquad (39) $$

$$ \frac{\partial S_A}{\partial \bar{z}_t}(0) = \beta(1 - S_A)S_A \qquad (40) $$

Thus the rates of adaptation per unit of time and per generation are

$$ \mathbf{RA}_{\text{T}} = h^2 \mathbf{V}_P \beta \frac{1}{T} \frac{S_A(1 - S_A)}{\lambda - S_A} \quad \text{and} \quad \mathbf{RA}_{\text{G}} = h^2 \mathbf{V}_P \beta \frac{S_A(1 - S_A)}{\lambda - S_A} \qquad (41) $$

The strength of selection is given by $V_{\theta} = \beta(1 - S_A)$ and the rates of adaptation per unit of selection are

$$ \mathbf{RA}_{\text{TS}\theta} = h^2 \mathbf{V}_P \frac{1}{T} \frac{S_A}{\lambda - S_A} \quad \text{and} \quad \mathbf{RA}_{\text{GS}\theta} = h^2 \mathbf{V}_P \frac{S_A}{\lambda - S_A} \qquad (42) $$

The strength of individual selection $V_{R_0}$ is

$$ V_{R_0} = \beta S_A \qquad (43) $$

and the adaptation rates per individual selection units are

$$ \mathbf{RA}_{\text{TSR}_0} = h^2 \mathbf{V}_P \frac{1}{T} \frac{(1 - S_A)}{\lambda - S_A} \quad \text{and} \quad \mathbf{RA}_{\text{GSR}_0} = h^2 \mathbf{V}_P \frac{(1 - S_A)}{\lambda - S_A} \qquad (44) $$

If the growth rate $\lambda$ equals 1, then the adaptation rate per selection units per generation does not depend on the vital rates.

## Adaptation rate under prebreeding census

Under a prebreeding census, the fertility ($F$) corresponds to the number of offspring produced per female per year (i.e., fecundity, denoted as $f$) that survive from time $t$ to $t+1$ (with probability $S_J$ in our case), and it is defined by

$$ F = f S_J. \qquad (45) $$

In this situation the rate of adaptation will change when selection acts on juvenile survival and we describe them when selection acts on fecundity $f$.
### S4.0.1 Juvenile survival $S_J$ (fecundity case)

When selection acts on juvenile survival $S_J$, the selection gradient will take the following form

$$ \frac{\partial \lambda}{\partial S_J} = \sum_{i_1, i_2=1}^{2} S_{i_1, i_2} \frac{\partial(\mathbf{R} + \mathbf{U})_{i_1, i_2}}{\partial S_J} = S_{1,1}(1-\gamma) + S_{2,1}\gamma + S_{1,2}f = \frac{\lambda}{S_J} \left( \frac{v_1 w_1}{v^\top w} + \frac{1}{T} \right) \quad (46) $$

$$ \frac{\partial S_J}{\partial \bar{z}_t}(0) = \beta(1 - S_J)S_J \quad (47) $$

Thus, the rates of adaptation per unit of time and per generation become

$$ \mathbf{RA}_{\mathbf{T}} = h^2 \mathbf{V}_P \beta \frac{1}{T} \left( 1 + \frac{\lambda}{\lambda - S_J(1-\gamma)} \right) (1-S_J) \quad \text{and} \quad \mathbf{RA}_{\mathbf{G}} = h^2 \mathbf{V}_P \beta \left( 1 + \frac{\lambda}{\lambda - S_J(1-\gamma)} \right) (1-S_J) \quad (48) $$

The strength of selection is unchanged, $V_\theta = \beta(1 - S_J)$, but the rates of adaptation per unit of selection become

$$ \mathbf{RA}_{\mathbf{TS}\theta} = h^2 \mathbf{V}_P \frac{1}{T} \left( 1 + \frac{\lambda}{\lambda - S_J(1-\gamma)} \right) \quad \text{and} \quad \mathbf{RA}_{\mathbf{GS}\theta} = h^2 \mathbf{V}_P \left( 1 + \frac{\lambda}{\lambda - S_J(1-\gamma)} \right) \quad (49) $$

The strength of individual selection $V_{R_0}$ is modified as follows

$$ V_{R_0} = \beta(1 - S_J) \left( 1 + \frac{1}{1 - S_J(1 - \gamma)} \right) \quad (50) $$

and the rates of adaptation per units of individual selection become

$$ \mathbf{RA}_{\mathbf{TSR}_0} = h^2 \mathbf{V}_P \frac{1}{T} \frac{\left( 1 + \frac{\lambda}{\lambda - S_J(1-\gamma)} \right)}{1 + \frac{1}{1 - S_J(1-\gamma)}} \quad \text{and} \quad \mathbf{RA}_{\mathbf{GSR}_0} = h^2 \mathbf{V}_P \frac{\left( 1 + \frac{\lambda}{\lambda - S_J(1-\gamma)} \right)}{1 + \frac{1}{1 - S_J(1-\gamma)}} \quad (51) $$

If the growth rate $\lambda$ is equal to 1, then the rate of adaptation per selection units per generation is independent of the vital rates.

### S4.0.2 Fecundity $f$

If the selection occurs through the fecundity of adults $f, f(z) = \exp(\log(f) + \beta z)$, then the selection gradient satisfies

$$ \frac{\partial \lambda}{\partial f} = \sum_{i_1, i_2=1}^{2} S_{i_1, i_2} \frac{\partial(\mathbf{R} + \mathbf{U})_{i_1, i_2}}{\partial f} = S_{1,2}S_J = \frac{S_J v_1 w_2}{v^\top w} \quad (52) $$

$$ \frac{\partial f}{\partial \bar{z}_t}(0) = \beta f \quad (53) $$

The generation time satisfies in this case the following formula

$$ T = \frac{\lambda v^\top w}{v^\top \mathbf{R} w} = \frac{\lambda v^\top w}{v_1 w_2 f S_J} \quad (54) $$

Thus the rates of adaptation per unit of time and per generation become

$$ \mathbf{RA}_{\mathbf{T}} = h^2 \mathbf{V}_P \beta \frac{1}{T} \quad \text{and} \quad \mathbf{RA}_{\mathbf{G}} = h^2 \mathbf{V}_P \beta \quad (55) $$

When we make the assumption that the selective pressure ($\beta$), heritability ($h^2$), and phenotypic variance ($V_P$) are the same between species, then the rates of adaptation per generation will be equal as well.

The relative rate of change in fecundity $f, V_\theta$ and the relative change in individual fitness $V_{R_0}$ are also equal in this case $V_\theta = V_{R_0} = \beta$, and the rates of adaptation per unit of selection and per individual selection unit are equal and they satisfy

$$ \mathbf{RA}_{\mathbf{TS}\theta} = \mathbf{RA}_{\mathbf{TSR}_0} = h^2 \mathbf{V}_P \frac{1}{T} \quad \text{and} \quad \mathbf{RA}_{\mathbf{GS}\theta} = \mathbf{RA}_{\mathbf{GSR}_0} = vT/V_\theta = h^2 \mathbf{V}_P, $$
## S5 Computation time

In this study, we conducted a comparative analysis of the computational efficiency between Individual-Based Models (IBM) and Matrix Population Models (MPM) by evaluating the time required to simulate population dynamics over a century. The results indicate that, on average, MPMs significantly outperform IBMs in terms of speed, except for species 5, which is characterized by the longest life expectancy and highest generation time among the species studied. In particular, MPMs were able to provide theoretical insights into the adaptation rate for four different species, in merely 0.04 seconds. This stark contrast in computational efficiency underscores the potential advantages of utilizing MPMs for extensive population studies.

<table><caption>**Table S2:** The time (seconds) taken to compute one population simulation over 100 years was calculated. IBM time was determined by averaging 20 replicates and initialization with 200 or 2000 individuals. Note that, due to their stochasticity, IBM simulations should typically be replicated 100 times or more to compute expected dynamics. MPM time refers to the duration for which the model runs for a single vital rate, calculated as the median across the four vital rates. It is worth mentioning that the theoretical findings provide the adaptation rate for the four species and significant demographic results in 0.04 seconds. Provided times are estimates, subject to variation with computers, coding languages, or optimization.</caption><thead><tr><td>Species</td><td>IBM (200 individuals)</td><td>IBM (2000 individuals)</td><td>MPM</td></tr></thead><tbody><tr><td>1</td><td>1840</td><td>18500</td><td>42.5</td></tr><tr><td>2</td><td>592</td><td>9260</td><td>42.5</td></tr><tr><td>3</td><td>23.0</td><td>1210</td><td>42.9</td></tr><tr><td>4</td><td>74.2</td><td>981</td><td>42.5</td></tr><tr><td>5</td><td>1.45</td><td>36.9</td><td>42.1</td></tr></tbody></table>