# Mutant emergence timing and population immunisation status impact epidemiological dynamics

Bastien Reyné, Ramsès Djidjou-Demasse, Mircea T Sofonea, Samuel Alizon

► To cite this version:

Bastien Reyné, Ramsès Djidjou-Demasse, Mircea T Sofonea, Samuel Alizon. Mutant emergence timing and population immunisation status impact epidemiological dynamics. Journal of Theoretical Biology, 2025, 608, pp.112140. 10.1016/j.jtbi.2025.112140. hal-05069967

HAL Id: hal-05069967

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

Submitted on 16 May 2025

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_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

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

# Mutant emergence timing and population immunisation status impact epidemiological dynamics

Bastien Reyné ID a,b,*, Ramsès Djidjou-Demasse ID a,c,1, Mircea T. Sofonea ID d,e,1,
Samuel Alizon ID f,1

## ARTICLE INFO

Keywords:

Evolutionary epidemiology

Evolutionary invasion analysis

SARS-CoV-2

Immune escape

Non-equilibrium dynamics

Immunity waning

## ABSTRACT

A key question in evolutionary epidemiology is to determine differences in the conditions that may allow some mutant strains to spread in a population where a resident strain is already circulating. Evolutionary invasion analyses assume that the immunity is long-lasting for previously infected individuals making it difficult to study traits such as immune escape. We relax this last assumption and allow the environment faced by the mutant to fluctuate outside of any epidemiological equilibrium. We introduce an original two-strains non-Markovian model that accounts for realistic immunity waning and cross-immunity, inspired by the case of SARS-CoV-2 variants. We show that mutants with increased contagiousness or with some immune escape abilities are more likely to invade the population. We also show that the timing of the introduction of mutant strain in the population is key because it is associated with the population's immunisation status. Our results underline the importance of immune waning and non-equilibrium dynamics on infectious disease evolution.

## 1. Introduction

Evolutionary and epidemiological processes frequently overlapped during the SARS-CoV-2 pandemic. It started as soon as mid-2020 with the D614G strain, which was associated with increases in viral loads (Volz et al., 2021) and continued with the evolution of the variants of concern (VOCs) (WHO, 2023). These VOCs had increased contagiousness (e.g. Davies et al. (2021)), higher virulence (e.g. Fisman and Tuite (2021)), and vaccination or infection-derived immune escape (e.g. An-deweg et al. (2022); Faria et al. (2021)). Years later, many variants with different phenotypic traits continue to impact the overall epidemiolog-ical dynamics (Markov et al., 2023; Lythgoe et al., 2023), as seen for other infectious diseases such as HIV (Wymant et al., 2022).

As illustrated by the case of SARS-CoV-2, a variant's evolutionary success can stem from different life-history traits, such as an increase in contagiousness (the number of secondary cases and the generation time distribution (Blanquart et al., 2022)) or some immune escape properties (Markov et al., 2023). Interactions between these traits or with the envi-ronment in which variants evolved can further contribute to increasing

invasion fitness. Prospective SARS-CoV-2 evolutionary epidemiology works mainly focused on competition between variants on a short-term scale (e.g. Park et al. (2022); Blanquart et al. (2022); Otto et al. (2021)) as the adaptive landscape was remodelled by each new variant, vaccine campaigns, or non-pharmaceutical intervention. It is possible to retro-spectively explain the appearance of new VOCs if there is data available (e.g. Benhamou et al. (2023)), although the analyses are rapidly quantitatively limited when immune escape is at play (Pearson et al., 2021). However, determining in advance what will be the phenotypic characteristics of a future circulating variant is less obvious and remains highly prospective. Yet, this endeavour is necessary to anticipate the emergence of future variants to optimise public health decisions.

Evolutionary invasion analyses are often used to determine under which conditions a variant can spread in a population as a function of its traits. These typically consider a population where an already es-tablished strain, the 'resident', circulates and assesses the ability of a newly introduced strain, generally called the 'mutant' (or *variants* in the SARS-CoV-2 case), to replace it (Geritz et al., 1998; Brannstrom et al., 2013). One of the strong assumptions made by many of these

<footnote>ªMIVEGEC, Univ. Montpellier, IRD, CNRS, Montpellier, France</footnote>

<footnote>ᵇUniv. Bordeaux, INSERM, INRIA, BPH, U1219, Bordeaux, F-33000, France</footnote>

<footnote>ᶜÉcole Polytechnique de Thiès, Thiès, Sénégal</footnote>

<footnote>ᵈPCCEI, Univ. Montpellier, INSERM, Montpellier, France</footnote>

<footnote>ᵉDepartment of Anesthesiology, Critical Care, Intensive Care, Pain and Emergency Medicine, CHU Nîmes, Nîmes, France</footnote>

<footnote>fCIRB, Collège de France, CNRS, INSERM, Université PSL, Paris, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: bastien.reyne@ird.fr (B. Reyné).</footnote>

<footnote n="1">¹ Equal contribution.</footnote>
frameworks that we need to alleviate is that of the separation between epidemiological and evolutionary timescales, meaning that the resident is assumed to reach an epidemiological equilibrium before the mutant appears (Brannstrom et al., 2013; Buckingham and Ashby, 2023). In the SARS-CoV-2 pandemic, for example, such an equilibrium has not yet been reached since variants with different traits still emerge frequently (Fig. A.5). Furthermore, each epidemic wave, whether or not it was caused by a new variant, reshaped the population immunisation level (Hozé et al., 2021).

Modelling population immunity is challenging for two main reasons. First, quantifying immunity is not feasible *per se*. In particular, we may only observe correlates of protection at the individual level (e.g. antibody titres), which are highly variable (Zachreson et al., 2023) or differences between different groups at a population level (through hazard or risk ratios) where observational studies may fail to identify causality (Monge et al., 2023). Some SARS-CoV-2-focused studies eventually provided efficacy estimations (at a population level) of vaccine and infection-derived immunity two or three years after the beginning of the pandemic (Powell et al., 2022; Stein et al., 2023). Second, incorporating empirical decline of immunity in mathematical models requires having a formalism that accounts for the time since clearance to make the immunity efficacy function of that time. We have previously shown that non-Markovian formalisms are convenient to capture dynamics of within-host processes showing non-stationary rates (Sofonea et al., 2021; Reyné et al., 2022), and these are well-adapted to model immunity efficacy.

In addition to immunity waning, it is important to also take into account immune escape to anticipate mutant invasion. The evolution of this phenotypic trait has been extensively studied at the molecular level (Vossen et al., 2002; Schmid-Hempel, 2008), including for SARS-CoV-2 (Starr et al., 2022; Volz et al., 2021), but it can have major epidemiological consequences. In the case of SARS-CoV-2, this was again illustrated by the high prevalence provoked by Omicron first sub-variants, which was mainly attributable to immune escape abilities that were later documented (Stein et al., 2023).

Here, we study how the interplay between population immunity and the timing of a mutant appearance shapes its evolutionary success using an evolutionary invasion analysis model that accounts for non-equilibrium phenomena. For simplicity, we consider two broad categories of infection life-history traits, namely intrinsic transmissibility and immunity, each encompassing a variety of sub-traits. An advantage in contagiousness can be due to a shorter generation time, a higher growth rate, or a higher reproduction number. Note that the latter two are commonly used as fitness metrics (Keeling and Rohani, 2008; Blanquart et al., 2022). Regarding immunity, a strain advantage may originate from different levels of immune waning. The strains are also distinguished by their immune escape properties from one another.

We define mutant evolutionary success as the ability to produce the next epidemic wave following its introduction. This departs from the canonical choice for the usual fitness metrics, which assumes that resident epidemiological dynamics are at equilibrium, thereby guaranteeing an equivalence between having the mutant reproduction number higher than one and its growth rate being non-negative (Hurford et al., 2010). We chose to depart from this core assumption because in the context we study, the population immune landscape is perpetually changing. These changes over time can be seen as environmental feedbacks, which are a key component in evolutionary epidemiology (Lion and Metz, 2018).

Our model allows us to show that contagiousness and immune escape both impact the risk of mutant emergence, providing an adequate timing in the mutant arrival regarding the population immunisation status.

## 2. The model

We assume a compartmental model with two strains, the resident and the mutant, respectively referred to by subscripts 1 and 2, using partial differential equations (PDEs, Fig. 1). This allows us to implement a

<figure><img src="image_4.png" /><figcaption>**Fig. 1.** Model flowchart. Susceptible individuals *S* are infected by an individual infected with resident strain *I*₁ or the mutant strain *I*₂. We record the time since infection *τ* for infected individuals. They eventually recover and end up in the recovered compartments *R*₁ or *R*₂. We record the time since clearance *s* for recovered individuals, which might be reinfected depending on the age of their immunity. Grey arrows show the demography-related transitions (births or deaths).</figcaption></figure>

non-Markovian structure of the infected (*I*₁ and *I*₂) and recovered (*R*₁ and *R*₂) compartments to account for time since infection in the former and since clearance in the latter. Practically, we define parametric generation time distributions in the infected compartments and time-dependent immunity efficacy in the recovered compartments. Based on longitudinal studies of patients infected by SARS-CoV-2 (Lind et al., 2023), we further assume immunity to be ‘leaky’, meaning recovered individuals may be reinfected by the two strains, although at a lower rate than susceptible individuals.

Mathematically, if we denote the time since infection by $\tau$, infected individuals enter the $I_1$ compartment according to the following boundary condition

$$I_1(t, \tau = 0) = \lambda_1(t) \left[ S(t) + \int_0^\infty R_1(t, s)[1 - \xi_{1,1}(s)]ds + \int_0^\infty R_2(t, s)[1 - \xi_{2,1}(s)]ds \right], \quad (1)$$

where $\lambda_1$ is the force of infection of the resident strain, *S* the susceptible individuals, $R_1$ (resp. $R_2$) individuals that were previously infected by the resident (resp. the mutant). The function $s \mapsto 1 - \xi_{x,y}(s)$ represents the decrease in susceptibility against strain *y* provided by a previous infection by strain *x*, *s* days after the clearance. Put in another way, the integral

$$\int_0^\infty R_2(t, s)[1 - \xi_{2,1}(s)]ds$$

represents individuals previously infected by strain 2 (the mutant), weighted by the efficacy of cross-immunity against the resident strain.

Symmetrically, individuals enter the $I_2$ compartment under the following boundary condition

$$I_2(t, \tau = 0) = \lambda_2(t) \left[ S(t) + \int_0^\infty R_1(t, s)[1 - \xi_{1,2}(s)]ds + \int_0^\infty R_2(t, s)[1 - \xi_{2,2}(s)]ds \right]. \quad (2)$$

Infected individuals clear the infection with a recovery rate, $\gamma_1$ or $\gamma_2$ depending on the strain, and enter recovered compartments at time $s = 0$ (i.e., newly recovered individuals), as formalised by the boundaries conditions

$$R_1(t, s = 0) = \int_0^\infty I_1(t, \tau) \gamma_1(\tau) d\tau, \quad (3)$$

$$R_2(t, s = 0) = \int_0^\infty I_2(t, \tau) \gamma_2(\tau) d\tau. \quad (4)$$
These boundary conditions Eqs. (1)–(4) are coupled with a dynamical system for the departures of the different compartments:

$$ 
\begin{aligned}
\frac{\partial S(t)}{\partial t} &= \mu N(t) - \lambda(t) S(t) - \mu S(t), \\
\left( \frac{\partial I_1(t, \tau)}{\partial t} + \frac{\partial I_1(t, \tau)}{\partial \tau} \right) &= -\gamma_1(\tau) I_1(t, \tau) - \mu I_1(t, \tau), \\
\left( \frac{\partial I_2(t, \tau)}{\partial t} + \frac{\partial I_2(t, \tau)}{\partial \tau} \right) &= -\gamma_2(\tau) I_2(t, \tau) - \mu I_2(t, \tau), \\
\left( \frac{\partial R_1(t, s)}{\partial t} + \frac{\partial R_1(t, s)}{\partial s} \right) &= -\lambda_1(t) [1 - \xi_{1,1}(s)] R_1(t, s) \\
&\quad - \lambda_2(t) [1 - \xi_{1,2}(s)] R_1(t, s) \\
&\quad - \mu R_1(t, s), \\
\left( \frac{\partial R_2(t, s)}{\partial t} + \frac{\partial R_2(t, s)}{\partial s} \right) &= -\lambda_1(t) [1 - \xi_{2,1}(s)] R_1(t, s) \\
&\quad - \lambda_2(t) [1 - \xi_{2,2}(s)] R_2(t, s) \\
&\quad - \mu R_2(t, s),
\end{aligned}
$$

where $\mu$ is the birth/death rate, $N(t)$ the population size, and $\lambda$ the overall force of infection. In particular,

$$ N(t) = S(t) + \int_{0}^{\infty} [I_{1}(t, \tau) + I_{2}(t, \tau)] d\tau + \int_{0}^{\infty} [R_{1}(t, s) + R_{2}(t, s)] ds, $$

and

$$ \lambda_x(t) = \beta_x \int_0^\infty \omega_x(\tau) I_x(t, \tau) d\tau, \quad \lambda(t) = \lambda_1(t) + \lambda_2(t), $$

where $\beta_x$ is the transmission rate for strain x and $\omega_x(\cdot)$ is the generation time for strain x.

## 2.1. Model parameterization

### 2.1.1. The resident strain

The resident strain's key parameter values are chosen to reflect those of Omicron BA.2. This is justified by the fact that this sub-variant was the first to circulate in populations with high immunity (induced by past infections, vaccination or both) and limited non-pharmaceutical interventions. Based on estimates from empirical distributions ((UKHSA, 2022), see Appendix B and Fig. B.8), we assume a Gamma distribution for the generation time,

$$ \omega_1 \sim \text{Gamma}(\text{shape} = 1.67, \text{scale} = 1.94). \quad (5) $$

We also assume a basic reproduction number ($R_{0,1}$) of 4 for the resident (cf. Appendix C). Importantly, we introduce an upper limit to the force of infection to account for the fact that, in the majority of countries, health policies were implemented to mitigate the propagation and avoid the collapse of the healthcare system (more details in Appendix C). As shown in Fig. 2(a), the resident strain infection-derived immunity against itself (i.e. the immunity waning) is modelled with linear decreasing functions:

$$ \xi_{1,1}(s) \mapsto \max \{1 - s \cdot \alpha_{1,1} \cdot X_{1,1}, \zeta\}, \\ \text{with } \alpha_{1,1} = 0.001, X_{1,1} \sim \mathcal{U}[0.85, 1.15], \text{ and } \zeta \in [0; 0.37], $$

where $\alpha_{1,1}$ is a coefficient chosen to qualitatively reproduce the immunity waning reported in Stein et al. (2023), $X_{1,1}$ a random variable to allow up to 15% variations around the baseline value to conduct later the sensitivity analysis and $\zeta$ represents the long-term residual immunity efficacy, which is observed in longitudinal follow-ups (Stein et al., 2023) (see Appendix B for details).

The resulting dynamic for the resident strain alone is shown in Fig. 2(c).

### 2.1.2. The mutant strain

We explore three different basic reproduction numbers for the mutant, with $R_{0,2} \in \{2, 3, 4\}$. We assume the mutant generation time follows a generation time that may differ from that of the resident and is defined by

$$ \omega_2 \sim \text{Gamma}(\text{shape} = 1.67 \cdot G_1, \text{scale} = 1.94 \cdot G_2), \text{ with} $$

$$ G_1, G_2 \sim \mathcal{U}[0.75, 1.25]. $$

For parsimony reasons, we assume the same parametrization for the mutant's immunity waning than for the resident, namely:

$$ \xi_{2,2}(s) \mapsto \max \{1 - s \cdot \alpha_{2,2} \cdot X_{2,2}, \zeta\}, \\ \text{with } \alpha_{2,2} = 0.001, X_{2,2} \sim \mathcal{U}[0.85, 1.15], \text{ and } \zeta \in [0; 0.37]. $$

### 2.1.3. Cross-immunity and immune escape

One of the reasons why investigating the role of immune escape in variant emergence is so delicate is that, in addition to capturing waning immunity, it requires defining how this immunity reacts to the different strains.

For simplicity, we assume cross-immunity to be the protection against a strain conferred by the immunity developed after the infection by another strain (the blue curve in Fig. B.9). Immune escape can then be defined as the difference between the level of immunity provided by the strain that originated the last infection and the cross-immunity provided against the other strain.

Regarding the parameterization of the level of protection conferred by a mutant infection against the resident, we assume the same function as for the natural waning:

$$ \xi_{2,1}(s) \mapsto \max \{1 - s \cdot \alpha_{2,1} \cdot X_{2,1}, \zeta\}, \\ \text{with } \alpha_{2,1} = 0.001, X_{2,1} \sim \mathcal{U}[0.85, 1.15], \text{ and } \zeta \in [0; 0.37]. $$

Finally, for cross-immunity derived from a resident infection, we again assume the same general form of immunity efficacy decay:

$$ \xi_{1,2}(s) \mapsto \max \{1 - s \cdot \alpha_{1,2} \cdot X_{1,2}, \zeta\}, \\ \text{with } X_{1,2} \sim \mathcal{U}[0.85, 1.15], \text{ and } \zeta \in [0; 0.37]. $$

We explore three scenarios for the parameter $\alpha_{1,2}$ to account for immune escape possibilities suggested by empirical data (Stein et al., 2023). We name each scenario according to the level of immune escape subsequently implied: None ($\alpha_{1,2} = 0.001$), Moderate ($\alpha_{1,2} = 0.003$) and High ($\alpha_{1,2} = 0.006$). The different levels of protection conferred by immunity are summarised in Fig. 2(b) and our assumptions are further detailed in Appendix B.

### 2.1.4. Mutant introduction

The mutant is introduced after the resident's first wave but before the population reaches an endemic equilibrium. This assumption is important to understand the role of waning immunity and is consistent with the fact that no SARS-CoV-2 variant has yet reached an equilibrium before being replaced by another variant.

We further assume that the mutant is introduced once the immunity in the population is below the herd immunity threshold, given by $1 - 1/R_0^m$. This threshold is purely theoretical (see Appendix C) and should be considered as transient in a population with waning immunity. The level of immunity in the host population against a mutant strain m at time t is

$$ p_m(t) := 1 - S(t) - \int_0^\infty [1 - \xi_{1,2}(s)] R_1(t, s) ds, $$

and, by definition, the mutant cannot theoretically invade if it is introduced when $p_m(t) > 1 - 1/R_{0,2}$. Practically, we first perform a run with only the resident strain (Fig. 2(c)) to identify the time when $p_m$ passes below the herd immunity threshold (Fig. 2(d)). Mutants are then introduced at that time by default.

There are many ways to define the evolutionary success of a strain, often with respect to another. Here, we define mutant success as the ability to produce the next epidemic wave.

## 2.2. Sensitivity analyses

To further delineate the relative importance of the studied parameters, we perform global sensitivity analyses on the ones that drive the phenotypic traits, namely the mutant generation time and the immunity
for the two strains. To this end, we use Sobol indices that decompose the output's variance according to each input parameter (Saltelli et al., 2008). More specifically, the Sobol index

$$S_X = \frac{\mathrm{Var}(E[Y|X])}{\mathrm{Var}(Y)},$$

corresponds to the variance of the output Y which is due solely to the parameter X. In our case, Y is the mutant density and X can be any of the input parameters related to the strains' phenotypic traits. These indices are computed at each time point and for each scenario. Additional details are provided in Appendix B & Appendix D.

Given that the mutant introduction date is fixed in our default scenarios, we also test different timings for the mutant introduction with respect to the resident's dynamics.

# 3. Results

## 3.1. Temporal herd immunity

The nature of the environment experienced by a mutant is essential to predict its ability to invade. Before the mutant introduction in the population, the level of immunity it faces depends on the level of population immunisation status (against the resident) which vary over time. In our model, the immunity waning (Fig. 2(a)) induces various epidemic waves in the resident population dynamics (Fig. 2(c)) that generates temporal variation in this environment through the level of cross-immunity against a mutant (Fig. 2(d)). The intensity of the variations in the temporal immunity against the mutant (Fig. 2(d)) depends on the level of cross-immunity assumed (Fig. 2(b)). This implies the requirement for a mutant to invade are dependant of the introduction date.

Theoretically, to invade, a mutant must face a host's herd immunity below some threshold at the time of invasion. If it has high immune escape properties, the mutant may invade the population with an $R_{0,2}$ of 2 during extended time windows (purple line in Fig. 2(d)). Conversely, without immune escape, even with an $R_{0,2}$ of 4, the mutant can only invade the population during very small time window (yellow line in Fig. 2(d)).

Points in Fig. 2(d) indicate the timing of mutant introductions in the default scenario for varying levels of cross-immunity with the resident. As explained in the methods, we assume that the mutant is introduced as soon as the (maximum) immunity it faces crosses the herd immunity threshold.

## 3.2. Next epidemic wave

We study the ability of a mutant to produce the next epidemic wave, which is our metric of evolutionary success, under three immune escape scenarios (None, Moderate, High) and three mutant reproduction numbers $R_{0,2} \in \{2, 3, 4\}$. As expected, mutants with a lower reproduction number that do not present any immune escape properties do not produce any epidemic wave (Fig. 3(a) and (b)).

Conversely, if its reproduction number is high enough ($R_{0,2} = 4$), the mutant is always present, but it may not be dominant in the epidemic wave if its immune escape is limited (Fig. 3(c)).

Finally, mutants with immune escape produce the next epidemic wave, even with a low reproduction number (Fig. 3(g) and (h)).

In all the scenarios, the uncertainty is high among the replicates, i.e. depending on the parameter combinations used (see Appendix B for the origin of the variance and quantiles intervals construction). The sensitivity analyses further show the variance is explained by different parameters depending on the scenario (see Appendix D).

Finally, note that in some cases, even though the mutant generates the second epidemic wave, the resident may produce later waves (e.g. Fig. 3(d) and (g)).

## 3.3. Impact of the timing of mutant introduction

Until now, we assumed that the introduction of the mutant occurred as soon as the population immunity it would face waned below the herd immunity threshold. We now introduce the mutant each day between day 66 (the earliest introduction date previously considered) and day 350 (during the resident's second wave in the absence of the mutant), with the baseline parameters.

The timing of introduction strongly impacts the height of the mutant density peak, as well as its date of occurrence (Fig. 4). A sub-optimal date of introduction, whether it is too early or too late, leads to a lower epidemic peak in all cases (Fig. 4, top panel). In some cases, at the time of the mutant peak, the resident has a higher density, suggesting that the mutant did not succeed although it did in the baseline scenario (e.g when $R_{0,2} = 2$ with high immune escape).

The time for the mutant to reach its epidemic peak also depends on its introduction date (Fig. 4, bottom panel). In particular, if the mutant is introduced too late, there is an interference with the resident second wave that strengthens the competition between strains. If the resident second wave happens, the mutant may be delayed and produce its epidemic peak only after the resident's second wave, leading to discontinuities (high immune escape with $R_{0,2} = 2$ or moderate immune escape with $R_{0,2} = 3$).

Note that these results cannot capture the initial mutant density decrease that can be observed when the mutant is introduced before the immunity crosses the herd immunity threshold. This numerical phenomenon is also responsible for the increase in the date of epidemic wave when introductions occur before 100 days (it is especially visible for low $R_{0,2}$). Additional details about deterministic modelling can be found in Appendix C.

# 4. Discussion

The ability of a viral mutant to invade a population already exposed to a resident strain depends on several life-history traits. Among these, the immunological state of the population is particularly important but has rarely been investigated because it raises methodological challenges and requires specific biological data. Here, we use models and data generated during the SARS-CoV-2 pandemic to determine how mutant emergence timing and host immunity interact to shape invasion probabilities.

We used numerical simulations to assess the ability of a mutant to invade a population with realistic waning immunity. We considered mutants with a variety of traits and different timing of their introduction to determine how mutant traits and the timing of their introduction affect their evolutionary success.

One of our first findings is that variable immunity in the SARS-CoV-2 case could be explained by the sole effect of immunity waning, which induces recurrent epidemic waves (Fig. 2(c)) that would perpetually reshape the population immunisation level. This variable resident-derived immunity provides different levels of immunity against potential new mutants depending on the level of cross-immunity assumed (Fig. 2(d)). Overall, an increase in reproduction number or a decrease in cross-immunity level favours the mutant invasion fitness.

In our framework, the timing of the arrival of a mutant determines the environment it faces, thereby impacting the evolutionary outcome (Lion and Metz, 2018). The presented results could only be obtained thanks to a non-Markovian formalism with a realistic modelling of the immunity component. The use of PDEs to account for immunity waning dates back to the early days of infectious disease modeling (Kermack et al., 1932), and continued over the following decades (White and Medley, 1998) up to the present day (Reyné et al., 2022). A comparison between our model and a Markovian model is available in Appendix E). However, the importance of the timing must be minimised in a real-world setting where the immunisation status of the different populations has no reason to be synchronised, meaning the same mutant could
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Temporal resident-driven cross-immunity against the mutant before its introduction. (a) Waning of the immunity against the resident following its clearance, $\xi_{1,1}(s)$. (b) Cross-immunity against the mutant after a resident infection for the three immune escape scenarios, $\xi_{1,2}(s)$. (c) Resident infected dynamics, $I_1(t, \tau)$, assuming the immunity waning in Panel A. (d) Population immunisation against the mutant, $p_m(t)$, as a function of strains cross-immunity (Panel B), in a population with a resident-derived immunity. Grey dashed lines show the theoretical herd immunity thresholds corresponding to the required levels of population immunisation to prevent a mutant invasion depending on $R_{0,2}$. Points show the default mutant introduction dates, which correspond to the time the herd immunity passes below the threshold for each scenario (or at the lowest immunity for $R_{0,2} = 2$ or 3 when there is no immune escape). In all panels, thick lines correspond to median trajectories and shaded areas to 95% quantiles interval. See Fig. A.7 for interpretation.</figcaption></figure>

be reintroduced at many dates in a given population. Nevertheless, this work could provide a way of assessing the invasion of a mutant in a given population depending on the population's immunisation status and the potential need to strengthen hypothetical non-pharmaceutical interventions.

The choice of the fitness metric is always question-dependent but needs to be clearly stated (Roff, 2008). Two classical metrics in evolutionary epidemiology are the reproduction number and the epidemic growth rate. In fact, when considering the canonical framework of evolutionary invasion analysis, the resident reaches an endemic equilibrium before the mutant is introduced, which guarantees an equivalence between a positive growth rate and reproduction above unity (Hurford et al., 2010). However, we assumed a setting outside the equilibrium

regime, which imposed an arbitrary choice. We opted for the ability to produce the next epidemic wave, which was motivated by the dynamics of SARS-CoV-2 variants since Omicron BA.1 (Fig. A.6). It seemed more appropriate than the instantaneous growth rate at the time of introduction (there may not be a peak) or reproduction number (we work out of equilibrium). Our criterion to assess evolutionary success remains nonetheless debatable. It is not sufficient to determine a winner strain as there are some patterns of co-circulation. For example, in some settings, the mutant and the resident can produce an epidemic wave at the same time (Fig. 3(f)). It occurred during the SARS-CoV-2 pandemic in some countries, e.g. with the shared epidemic wave between Omicron BA.5 and Omicron BQ.1 (Fig. A.6). Another issue is that in some other cases, although the mutant produces the second epidemic wave, the resident
<figure><img src="image_6.png" /><figcaption>**Fig. 3.** Mutant dynamics in a population previously infected by the resident strain. Mutant dynamics are displayed in colours, and resident dynamics are always shown in black. Columns represent the different mutant reproduction numbers ($R_{0,2}$) considered. Rows represent the different immune escape considered. Shaded areas represent the 95% quantiles intervals.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. 4.** Effect of the mutant introduction date. The top panel shows the maximum mutant density reached in function of its introduction date. The solid lines represent the mutant strain, and the dotted lines represent the resident's density at the date of the mutant peak. The lower panel displays the date the mutant reaches its maximum, given its initial introduction date. The black dotted line displays the line $y = x$.</figcaption></figure>
causes the third one (Fig. 3(g)). Our success metric does not deal with these two cases.

Comparisons with other SARS-CoV-2 evolutionary modelling stud-
ies can help to anticipate how some assumptions we made (or did not
make) could impact our results. For instance, **Bushman et al.** (2021)
found that a variant with higher transmissibility was more likely to in-
fect more people than a partial immune escape. It could be explained
by a differentiated level of immunity in their model, with a distinction
between infection-derived immunity and vaccination, the latter being
assumed to be perfect. We did not make this assumption in our model
because, since then, both vaccine and infection-derived immunity effi-
cacies have been shown to decline over time (Powell et al., 2022; Stein
et al., 2023).

On another note, Blanquart et al. (2022); Park et al. (2022); Ben-
hamou et al. (2023) assumed a perfect immunity for infected individu-
als, but these works focused on the Alpha VOC earlier in the pandemic
when there was less than 15% that had been infected (Hozé et al., 2021).
Three years later, this assumption is not realistic any more. Our results
are different and less focused on the differences in growth rates as the
immune escape completely changed the nature of the implicit 'limited
resources assumption'.

Day et al. (2022) assumed a structured population between vacci-
nated and unvaccinated individuals. By doing so, they can differentiate
the ability of mutants to invade, given differences in uptake in the two
sub-populations and the level of vaccination in the population. They did
not consider any waning properties regarding immunity. In comparison,
our work neglects the differences between recovered and vaccinated in-
dividuals in favour of the age of immunity that presents empirically a
disparity between newly immunised individuals and the ones that de-
veloped immunity weeks ago (cf. Fig. B.10). More generally, our work
neglects any structure (besides age of immunity) in the population that
might allow a strain to subsist in sub-populations before establishing
itself in the overall population.

All the above-cited studies, ours included, make some additional
simplifications such as the lack of spatial structure, the homogeneity
in contact rates, or the lack of stochasticity that could happen in the
establishment of a strain (see also Appendix C). The consequence of al-
leviating these assumptions on mutant invasion has been investigated
in simpler contexts, i.e. without immune waning and at equilibrium,
which could be used as a basis to study their effect on these results in the
future.

All the works discussed so far considered only two strains, where
there were always many strains competing in reality. Some frameworks
allow for the whole diversity of strains to coexist. It is the case, for
instance, of work building on the Price equation that allows an arbitrary
number of strains to evolve in a population (Price, 1970). Some of these
models have focused in particular on variations in hosts (Gandon and
Day, 2009; Day et al., 2020). However, we do not believe it is adapted
in our case because it assumes traits to be defined by a parameter where
we considered functions for immunity.

Another future approach could be to focus on antigenic drift (Sasaki et al., 2022; Haraguchi and Sasaki, 1997) and to define the pathogens on an immunity continuum. A strain's immune evasion properties are determined by its antigenic distance from the particular strain that previously infected the hosts. This could explain the flu epidemics over the years (Gog and Grenfell, 2002) and could be of interest in the case of SARS-CoV-2.

In the line of Starr et al. (2022), improving evolutionary predictions
could benefit from estimating phenotypic changes caused by specific
mutations. It could allow us to assess the uptake of new mutants almost
instantaneously provided the current characteristics of the population,
such as the level of immunity in the population. Furthermore, assessing
the level of immunity of the population could be the first step through
robust surveillance system, because, as shown in this study, the adap-
tive landscape is evolving rapidly and yet important in the ability to
invade.

CRediT authorship contribution statement

**Bastien Reyné**: Writing – review & editing, Writing – original draft, Software, Methodology, Formal analysis, Data curation, Conceptualization; **Ramsès Djidjou-Demasse**: Writing – review & editing, Supervision, Software, Methodology; **Mircea T. Sofonea**: Writing – review & editing, Supervision, Methodology, Conceptualization; **Samuel Ali-zon**: Writing – review & editing, Writing – original draft, Supervision, Methodology, Conceptualization.

Code availability

All the scripts and data used in this study are available on the
following git repository: https://gitlab.com/reyb/sc2_invasion_analysis.
A permanent version of the code is also available on Zenodo: https:
//zenodo.org/doi/10.5281/zenodo.11519413.

Declaration of competing interest

The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.

Acknowledgments

We thank the ETE modelling team for fruitful discussions. The au-
thors acknowledge the ISO 9001 certified IRD i-Trop HPC (member of
the South Green Platform) at IRD Montpellier for providing HPC re-
sources that have contributed to the research results reported within
this paper. URL: https://bioinfo.ird.fr/ — http://www.southgreen.fr.
This work received funding from the Agence Nationale de la Recherche
(ANR) as part of the investment programme France 2030 under the ref-
erence ANR-21-EXES-0005, from the Occitanie Region; and from the
ExposUM Institute of the University of Montpellier (NEXUS EMIPSA),
as well as from the ANRS|Maladies Infectieuses Émergentes and France
2030 (ANRS-24-PEPRMIE-0003 PReViX).

**Appendix A. Supplementary figures**

**Appendix B. Technical details**

B.1. Parametrization

We detail here the parameters used to perform the results presented
in this study. We distinguished the parameters into two categories: the
ones with a direct interest in our topic that are allowed to vary to take
into account the variability in the results and some parameters that are
fixed in all simulations.

In the following, we try to list the different parameters in the same
order as they appear in the dynamical system.

B.1.1. *Fixed parameters*

The demography $\mu$ is fixed to $1/(365 \times 80)$ which corresponds to an average lifetime of 80 years for an individual. Notice that given the duration of the runs (700 days), it has very little impact on the dynamics. The population size $N(t)$ is fixed to one, for all $t \in [0; 700]$.

The recovery rates $\gamma_1(\cdot)$ and $\gamma_2(\cdot)$ are the same from a previous study of ours, which already applied to SARS-CoV-2 (Reyné et al., 2022). The median infection duration is ~7.4 days and all individuals recover in a maximum ~26 days.

The resident generation time $\omega_1$ follows a Gamma(shape = 1.67, scale = 1.94). It reproduces roughly the empirical serial interval provided by UKHSA (2022) shown on Fig. B.8. Note our formalism could have benefited from the empirical generation time without using a regular distribution function. We choose a parametrical distribution to allow the mutant generation time to be close to the resident's, with on top the possibility to vary.
<figure><img src="image_8.png" /><figcaption>Fig. A.5. Relative frequence of the different VOCs over time. Data from GISAID and CoVariants.org.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. A.6. Number of daily hospital admissions in France. The area under the curve is colored with the relative frequency of the different circulating VOCs in France. Frequency data from GISAID and CoVariants.org. Incidence data from Santé Publique France and data.gouv.fr.</figcaption></figure>

### B.1.2. Immunity

The immunity in our model is considered leaky, which seems appropriate for SARS-CoV-2 (Lind et al., 2023). That implies that recovered individuals will benefit from a decrease in susceptibility compared to susceptible individuals; they may however be infected again.

We consider two levels of protection for an individual: the immunity provided by the strain which infected the individual and a cross-immunity for the other strain. The infection-derived immunity is higher than the cross-immunity, as schematized on Fig. B.9.

We assumed a linear decay for the immunity efficacy to remain simple:

$$ \xi(s) \mapsto \max\{1 - s \cdot \alpha, \zeta\}, $$

where s is the time since clearance, α the decay slope and ζ a residual immunity. However, we used values for our decay that are close to values reported by Stein et al. (2023), as illustrated by Fig. B.10.

### B.1.3. Varying parameters and confidence intervals

We list here the different parameters that vary in our simulations. We have:

$$ \omega_2 \sim \text{Gamma}(\text{shape} = 1.67 \cdot x_1, \text{scale} = 1.94 \cdot x_2), \\ \text{with } x_1, x_2 \sim \mathcal{U}[0.75, 1.25], $$

$$ \xi_{1,1}(s) \mapsto \max\{1 - s \cdot \alpha_{1,1} \cdot x_3, \zeta\}, $$

$$ \xi_{1,2}(s) \mapsto \max\{1 - s \cdot \alpha_{1,2} \cdot x_4, \zeta\}, $$

$$ \xi_{2,1}(s) \mapsto \max\{1 - s \cdot \alpha_{2,1} \cdot x_5, \zeta\}, $$

$$ \xi_{2,2}(s) \mapsto \max\{1 - s \cdot \alpha_{2,2} \cdot x_6, \zeta\}, $$

$$ \text{with } x_3, x_4, x_5, x_6 \sim \mathcal{U}[0.85, 1.15] \text{ and } \zeta \sim \mathcal{U}[0, 0.37], $$

and $\alpha_{1,1} = \alpha_{2,1} = \alpha_{2,2} = 0.001$. The coefficient $\alpha_{1,2}$ depends on the scenario chosen: 0.001 (Immune escape: None), 0.003 (Moderate), 0.006 (High).

The parameters $x_1, x_2, x_3, x_4, x_5, x_6, \zeta$ are sampled with a Latin Hypercube Sample of size $n = 900$. (We draw n.sim = 100 unique values for each one of the 7 parameters, and the total number of runs is given by $n = n.sim \times (nb.params + 2)$.) Note the selected ranges are arbitrary, and chosen according to what we believe could be a good compromise between not having too much variance (which would be uninformative) and high enough variability to assess the impact of each parameter on the output through the sensitivity analyses.

We then perform 900 runs for each pair of scenarios (None, Moderate, High) and mutant reproduction number (2, 3, 4), hence we performed 9 batches of 900 runs. Confidence intervals are the 2.5% and 97.5% quantiles from the 900 runs, and the median run is the 50% quantile. The same 900 runs will be reused to compute the global sensitivity analyses (which are Sobol indices performed at each time step, see below).

## B.2. Numerical implementation

The code is available on the git repository https://gitlab.com/reyb/sc2_invasion_analysis.

The model is implemented in R v4.3.2 (Team, 2023). We used the following R packages:

* here for file managing/reproducibility (Müller, 2020),

* data.table for managing data (Barrett et al., 2025),

* Rcpp for implementing the dynamical system (Eddelbuettel and Balamuta, 2018),

* ggplot2 and cowplot for graphics (Wickham, 2016; Wilke, 2024),
<figure><img src="image_10.png" /><figcaption>**Fig. A.7.** Theoretical invasion time windows illustration. This plot provides a notice for Fig. 2(d). The blue line represents the level of immunity the mutant faces before its introduction (it depends on the immune escape scenario in our study). Each time this immunity is above the herd immunity threshold (red zones), the mutant cannot theoretically invade. The herd immunity threshold depends on the mutant reproduction number (an $R_{0,2} = 2$ on the left panel and $R_{0,2} = 3$ on the right panel); the more the mutant reproduction number is high the more the mutant will have possibilities to invade. For instance, at day 400, the mutant can invade if it has an $R_{0,2} = 3$, but cannot if it has an $R_{0,2} = 2$. See also Appendix C for discussion about the theoretical invasion in a deterministic setting.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. B.8.** Omicron BA.2 serial interval. In red is shown the empirical serial interval provided by UKHSA (2022). In black corresponds to the resident generation time used in our model.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Fig. B.9.** Schematic representation of immunity in this study.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Fig. B.10. Immunity implemented compared to Stein et al. (2023). Note that we refer in Stein et al. (2023) the “Protection against ancestral, alpha, and delta re-infection” as the waning in our study and the “Protection against omicron BA.2 re-infection” as the immune escape.</figcaption></figure>

* Ihs for generating the latin hypercube sample (Carnell, 2022),

* multisensi and sensitivity for performing the sensitivity analyses (Lamboni et al., 2011; Iooss et al., 2024),

* fitdistrplus for fitting the generation time on UKHSA data (Delignette-Muller and Dutang, 2015).

## Appendix C. Model limitations and assumptions criticisms

Beyond the classical limitations/assumptions that come with compartmental models (homogeneity in contact, no spatial structure, no heterogeneity in immunity profile), we list hereafter the main caveats we think are related to this study.

### C.1. Reproduction number

We mention throughout this study the basic reproduction number, $R_0$, for both the resident ($R_{0,1}$) and the mutant ($R_{0,2}$). Note this does not make much sense since the $R_0$ is defined in a completely susceptible population. Many epidemic waves later in many countries, populations around the world developed at least a partial immunity against SARS-CoV-2. The notion of $R_0$ is no more than a convenient way to parametrize the model and define “simply” a phenotypic trait quantifying an “intrinsic transmissibility”, even if it interplays with others such as the immunity waning baseline. We assumed an $R_{0,1} := 4$ for the resident strain and $R_{0,2} \in \{2, 3, 4\}$ for the mutant. It simply means if the strains were hypothetically introduced in a completely susceptible population, they would have the defined basic reproduction numbers.

Note also that $R_{0,1} := 4$ is an assumption since there are no estimates for the Omicron BA.2 reproduction number in the literature for the reasons listed hereinabove. This value is then arbitrary. However, it remains informative relative to the mutant reproduction number (that could hence have an advantage or disadvantage in its “intrinsic transmissibility”).

### C.2. Force of infection

As mentioned in the main text, we limit the force of infection to 0.04 in our model. The true formula of the force of infection for the strain x implemented is

$$ \lambda_x^{\text{true}}(t) := \min\{\lambda_x(t), 0.04\}. $$

It is done to be realistic, reflecting the non-pharmaceutical implementations that would occur if too many individuals were infected altogether. It then prevents everyone from being infected at the same time and smoothing the age of clearance after the first epidemic wave. However, the selected value is arbitrary and reflects a force of infection that was never topped during the first epidemic wave according to one of our previous SARS-CoV-2 transmission models (in Reyné et al. (2022) to be precise, even if the force of infections were never published). In this study, it results that there are always less than 25% individuals infected at the same time (cf. Fig. 2(c)). We show on Fig. C.11 what would have been the baseline resident's dynamic (without the mutant) if our upper limit (0.04) was decreased or increased up to 20%. It does not appear to be determinant in the overall dynamic shape. Ideally, this parameter would have been included in the sensitivity analysis, but it would have increased the number of simulations. We believe it would not change qualitatively the results nor their interpretation.

### C.3. Herd immunity threshold and mutant introduction

In our simulations, we introduced the mutant strains as soon as the herd immunity against the mutant passes below (if it happens) the corresponding threshold. We state it wouldn't be possible for a strain to invade if it were introduced when above the said threshold. However, it is not what happen in a deterministic model. When introduced above the threshold, the mutant strain density decreases until the herd immunity against the mutant diminishes below the threshold. In particular, the mutant instantaneous growth rate remains below 0 and the mutant temporal reproduction number remains below 1 as long as the immunity against the mutant is above the herd immunity threshold (see an
<figure><img src="image_14.png" /><figcaption>Fig. C.11. Impact of the force of infection upper limit on the resident's dynamic. The bold line represent the baseline parameter (0.04). Thin lines represent decrease or increase of respectively -20%, -10%, +10% & +20%.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. C.12. Illustration of different metrics taken from a run where the mutant is introduced when the immunity is above the herd immunity threshold. The red line shows the immunity against the mutant, and the shaded areas show when this immunity is above the herd immunity threshold. The mutant is introduced at time 66 (black line). The instantaneous growth rate (in blue) remains below 0 and the temporal reproduction number (in green) remains below 1 while the immunity is above the threshold, meaning the strain is decreasing as expected. However, due to the deterministic nature of the model the strain never goes extinct.</figcaption></figure>

example on Fig. C.12). There are then some cases where it happens, and the mutant still provokes an epidemic wave as shown on Fig. 4. It can be explained by the fact any strain goes extinct in our model and there is an epidemic reprisal as soon as the immunity against the mutant cross the herd immunity threshold.

## Appendix D. Sensitivity analyses

We performed global sensitivity analyses for each scenario. We compute Sobol indices for each time point of our dynamics, thus decompos-

ing the variance output according to the input parameters (see the Appendix of Reyné et al. (2022) or the vignette of the multisensi package for more details). The Sobol index for the input parameter X is given by

$$S_X = \frac{\mathrm{Var}(E[Y|X])}{\mathrm{Var}(Y)},$$

where E[Y|X] is the conditional expectation of Y with respect to X. We recommend interested readers to read Saltelli et al. (2008) for further explanations.
<figure><img src="image_16.png" /><figcaption>Fig. D.13. Sensitivity analyses for the different scenarios shown on Fig. 3.</figcaption></figure>

There was not enough variance to compute Sobol indices in the first two scenarios. The total curve represents the total variance at the given time and the colors represent the fraction of variance explained by each variable. Note the x-axes differ from one panel to another since the Sobol indices were computed only where there was variance.

We present the results on Fig. D.13. There is a diversity of patterns in the variables creating more variance given the different scenarios. It does not allow a global interpretation. We may still interpret some scenarios if taken on their own. For instance, Scenario A with $R_{0,2} = 4$, clearly indicates that the immunity of both strains against the resident drives the variance (and probably the evolutionary success). It can be explained by the fact that the two strains have the same mean parametrization in this scenario, and deviating from the mean immunity baseline might confer an advantage or a disadvantage against the other strain.

## Appendix E. Comparison between Markovian and non-Markovian models

We aim in this section to highlight the usefulness of non-Markovian formalism (PDEs in our case) compared to a simple ODEs-based Markovian model.

We consider the following ODE model

$$ \frac{dS}{dt} = -\lambda(t) S(t) - \mu S(t) + \mu N(t), $$

<figure><img src="image_17.png" /><figcaption>Fig. E.14. Illustration of non-Markovian (PDE) and Markovian (ODE) models.</figcaption></figure>
$$ \frac{dI}{dt} = \lambda(t) S(t) + \lambda(t) [1 - \xi] R(t) - \gamma I(t) - \mu I(t), $$

$$ \frac{dR}{dt} = \gamma I(t) - \lambda(t) [1 - \xi] R(t) - \mu R(t), $$

where $\lambda$ is the force of infection, $N(t)$ the total population, $\mu$ the birth/death rate, $\xi$ the reduction of susceptibility conferred by an immunity, $\gamma$ the recovery rate.

Regarding the parametrisation, we tried to choose parameters close to the ones selected in the PDE model:

* $\lambda(t) = \min\{0.04, \beta I(t)\}$,

* $R_0 = 4$,

* $\mu = 1/(365 \times 80)$,

* $\xi = 0.8$,

* $\gamma = 1/7$.

The Fig. E.14 shows the resident dynamics (without mutant) for both the non-Markovian (with the default parameters) and Markovian models. Particularly, the Markovian model shows only one epidemic wave whereas the non-Markovian shows three epidemic waves. This difference is mainly due to the way immunity is modelled in the non-Markovian model (leaky immunity with a temporal decay) compared with the Markovian model (leaky, but a constant efficacy). Other differences play a role in these two simulations: a constant recovery rate in the Markovian model or a deSolve implementation against Euler scheme handcrafted-implementation.

Please note that an ODE-based model could be non-Markovian provided chaining compartments (e.g. various infected compartments before getting recovered) and a good choice for the parameters associated with each compartment.

## References

Andeweg, S.P., Vennema, H., Veldhuijzen, I., Smorenburg, N., Schmitz, D., Zwagemaker, F., van Gageldonk-Lafeber, A.B., Hahné, S. J.M., Reusken, C., Knol, M.J., Eggink, D., on behalf of the SeqNeth Molecular surveillance group and, RIVM COVID-19 MOLECULAR EPIDEMIOLOGY GROUP, 2022. Elevated risk of infection with SARS-CoV-2 Beta, Gamma, and Delta variants compared with Alpha variant in vaccinated individuals. Sci. Trans. Med. 15 (684), eabn4338. Publisher: American Association for the Advancement of Science. https://doi.org/10.1126/scitranslmed.abn4338

Barrett, T., Dowle, M., Srinivasan, A., Gorecki, J., Chirico, M., Hocking, T., Schwendinger, B., Krylov, I., (2025). data.table: Extension of 'data.frame'. R package version 1.17.99. https://Rdatatable.github.io/data.table, https://github.com/Rdatatable/data.table, https://r-datable.com

Benhamou, W., Lion, S., Choquet, R., Gandon, S., 2023. Phenotypic evolution of SARS-CoV-2: a statistical inference approach. Evolution 77 (10), 2213-2223. https://doi.org/10.1093/evolut/qpad133

Blanquet, F., Hozé, N., Cowling, B.J., Débarre, F., Cauchemez, S., 2022. Selection for infectivity profiles in slow and fast epidemics, and the rise of SARS-CoV-2 variants. eLife 11, e75791. https://doi.org/10.7554/eLife.75791

Brannstrom, A., Johansson, J., Von Festenberg, N., 2013. The Hitchhiker's guide to adaptive dynamics. Games 4 (3), 304–328. https://doi.org/10.3390/g4030304

Buckingham, L.J., Ashby, B., 2023. The evolution of the age of onset of resistance to infectious disease. Bull. Math. Biol. 85 (5), 42. https://doi.org/10.1007/s11538-023-01144-5

Bushman, M., Kahn, R., Taylor, B.P., Lipsitch, M., Hanage, W.P., 2021. Population impact of SARS-CoV-2 variants with enhanced transmissibility and/or partial immune escape. Cell 184 (26), 6229–6242.e18. Publisher: Elsevier. https://doi.org/10.1016/j.cell.2021.11.026

Carnell, R., 2022. lhs: Latin hypercube samples. R package version 1.1.6. https://CRAN.R-project.org/package=lhs

Davies, N.G., Abbott, S., Barnard, R.C., Jarvis, C.I., Kucharski, A.J., Munday, J.D., Pearson, C. A.B., Russell, T.W., Tully, D.C., Washburne, A.D., Wenseleers, T., Gimma, A., Waites, W., Wong, K. L.M., Zandvoort, K.V., Silverman, J.D., Group, C. C.-W., Consortium, C.-G. U.K. C.-U., Diaz-Ordaz, K., Keogh, R., Eggo, R.M., Funk, S., Jit, M., Atkins, K.E., Edmunds, W.J., 2021. Estimated transmissibility and impact of SARS-CoV-2 lineage B.1.1.7 in England. Science 372 (6538). https://doi.org/10.1126/science.abg3055

Day, T., Kennedy, D.A., Read, A.F., Gandon, S., 2022. Pathogen evolution during vaccination campaigns. PLOS Biol. 20 (9), e3001804. Publisher: Public Library of Science. https://doi.org/10.1371/journal.pbio.3001804

Day, T., Parsons, T., Lambert, A., Gandon, S., 2020. The price equation and evolutionary epidemiology. Philos. Trans. Royal Soc. B 375 (1797), 20190357. Publisher: Royal Society. https://doi.org/10.1098/rstb.2019.0357

Delignette-Muller, M.L., Dutang, C., 2015. Fitdistrplus: an R package for fitting distributions. J. Stat. Softw. 64 (4), 1–34. https://doi.org/10.18637/jss.v064.i04

Eddelbuettel, D., Balamuta, J.J., 2018. Extending R with C + : a brief introduction to Rcpp. Am. Stat. 72 (1), 28–36. https://doi.org/10.1080/00031305.2017.1375990

Faria, N.R., Mellan, T.A., Whittaker, C., Claro, I.M., Candido, D. d.S., Mishra, S., Crispim, M. A.E., Sales, F. C.S., Hawryluk, I., McCrone, J.T., Hulswit, R. J.G., Franco, L. A.M., Raimundo, M.S., de Jesus, J.G., Andrade, P.S., Coletti, T.M., Ferreira, G.M., Silva, C. A.M., Manuli, E.R., Pereira, R. H.M., Peixoto, P.S., Kraemer, M. U.G., Gaburo, N., Camilo, C. d.C., Hoeltgebaum, H., Souza, W.M., Rocha, E.C., de Souza, L.M., de Pinho, M.C., Araujo, L. J.T., Malta, F. S.V., de Lima, A.B., Silva, J. d.P., Zauli, D. A.G., de S., F. A.C., Schnekenberg, R.P., Laydon, D.J., Walker, P. G.T., Schlüter, H.M., dos Santos, A. L.P., Vidal, M.S., Del Caro, V.S., Filho, R. M.F., dos Santos, H.M., Aguiar, R.S., Proença-Modena, J.L., Nelson, B., Hay, J.A., Monod, M., Miscouridou, X., Coupland, H., Sonabend, R., Vollmer, M., Gandy, A., Prete, C.A., Nascimento, V.H., Suchard, M.A., Bowden, T.A., Pond, S. L.K., Wu, C.-H., Ratmann, O., Ferguson, N.M., Dye, C., Loman, N.J., Lemey, P., Rambaut, A., Fraiji, N.A., Carvalho, M. d. P.S.S., Pybus, O.G., Flaxman, S., Bhatt, S., Sabino, E.C., 2021. Genomics and epidemiology of the P.1 SARS-CoV-2 lineage in Manaus, Brazil. Science 372 (6544), 815–821. Publisher: American Association for the Advancement of Science. https://doi.org/10.1126/science.abh2644

Fisman, D.N., Tuite, A.R., 2021. Evaluation of the relative virulence of novel SARS-CoV-2 variants: a retrospective cohort study in Ontario, Canada. CMAJ 193 (42), E1619–E1625. Publisher: CMAJ Section: Research. https://doi.org/10.1503/cmaj.211248

Gandon, S., Day, T., 2009. Evolutionary epidemiology and the dynamics of adaptation. Evolution 63 (4), 826–838. https://doi.org/10.1111/j.1558-5646.2009.00609.x

Geritz, S.A.H., Kisdie, E., MezéNA, G., Metz, J.A.J., 1998. Evolutionarily singular strategies and the adaptive growth and branching of the evolutionary tree. Evol. Ecol. 12 (1), 35–57. https://doi.org/10.1023/A:1006554906681

Gog, J.R., Grenfell, B.T., 2002. Dynamics and selection of many-strain pathogens. Proc. Nat. Acad. Sci. 99 (26), 17209–17214. Publisher: Proceedings of the National Academy of Sciences. https://doi.org/10.1073/pnas.252512799

Haraguchi, Y., Sasaki, A., 1997. Evolutionary pattern of intra-host pathogen antigenic drift: effect of cross-reactivity in immune response. Phil. Trans. Royal Soc. London. Ser. B 352 (1349), 11–20. Publisher: Royal Society. https://doi.org/10.1098/rstb.1997.0002

Hozé, N., Paireau, J., Lapidus, N., Tran Kiem, M., Salje, H., Severi, G., Touvier, M., Zins, M., de Lamballerie, X., Lévy-Bruhl, D., Carrat, F., Cauchemez, S., 2021. Monitoring the proportion of the population infected by SARS-CoV-2 using age-stratified hospitalisation and serological data: a modelling study. Lancet Public Health 6 (6), e408–e415. https://doi.org/10.1016/S2468-2667(21)00064-5

Hurford, A., Cowden, D., Day, T., 2010. Next-generation tools for evolutionary invasion analyses. J. Royal Soc. Interface 7 (45), 561–571. Publisher: Royal Society. https://doi.org/10.1098/rsif.2009.0448

Iooss, B., Veiga, S.D., Janon, A., Pujol, G., Broto, w. c. f.B., Boumhaout, K., Clouvel, L., Delage, T., Amri, R.E., Fruth, J., Gilquin, L., Guillaume, J., Herin, M., Idrissi, M.I., Le Gratiet, L., Lemaître, P., Marrel, A., Meynaoui, A., Nelson, B.L., Monarí, F., Oomen, R., Rakovec, O., Ramos, B., Rochet, P., Roustant, O., Sarazin, G., Song, E., Staum, J., Sueur, R., Touati, T., Verges, V., Weber, F., 2024. Sensitivity: global sensitivity analysis of model outputs and importance measures. R package version 1.30.0. https://CRAN.R-project.org/package=sensitivity.

Keeling, M.J., Rohani, P., 2008. Modeling Infectious Diseases in Humans and Animals. Princeton University Press. https://doi.org/10.2307/j.ctvmc4gk0

Kermack, W.O., McKendrick, A.G., Walker, G.T., et al., 1932. Contributions to the mathematical theory of epidemics. II. – The problem of endemicity. Proc. Royal Soc. London. Ser. A 138 (834), 55–83. Publisher: Royal Society. https://doi.org/10.1098/rspa.1932.0171

Lamboni, M., Monod, H., Makowski, D., 2011. Multivariate sensitivity analysis to measure global contribution of input factors in dynamic models. Reliab. Eng. Syst. Saf. 96.

Lind, M.L., Dorion, M., Houde, A.J., Lansing, M., Lapidus, S., Thomas, R., Yildirim, I., Omer, S.B., Schulz, W.L., Andrews, J.R., Hitchings, M. D.T., Kennedy, B.S., Richeson, R.P., Cummings, D. A.T., Ko, A.I., 2023. Evidence of leaky protection following COVID-19 vaccination and SARS-CoV-2 infection in an incarcerated population. Nat. Commun. 14 (1), 5055. Number: 1 Publisher: Nature Publishing Group. https://doi.org/10.1038/s41467-023-40750-8

Lion, S., Metz, J. A.J., 2018. Beyond R0 maximisation: on pathogen evolution and environmental dimensions. Trends Ecol. Evol. 33 (6), 458–473. Publisher: Elsevier. https://doi.org/10.1016/j.tree.2018.02.004

Lythgoe, K.A., Golubchik, T., Hall, M., House, T., Cahuantzi, R., MacIntyre-Cockett, G., Fryer, H., Thomson, L., Nurtay, A., Ghafani, M., Buck, D., Green, A., Trebes, A., Piazza, P., Lonie, L.J., Studley, R., Rourke, E., Smith, D., Bashton, M., Nelson, A., Crown, M., McCann, C., Young, G.R., Andre Nunes dos Santos, R., Richards, Z., Tariq, A., null, n., null, n., Fraser, C., Diamond, I., Barrett, J., Walker, A.S., Bonsall, D., 2023. Lineage replacement and evolution captured by 3 years of the United Kingdom Coronavirus (COVID-19) Infection Survey. Proc. Royal Soc. B 290 (2009), 20231284. Publisher: Royal Society. https://doi.org/10.1098/rspb.2023.1284

Markov, P.V., Ghafari, M., Beer, M., Lythgoe, K., Simmonds, P., Stilianakis, N.I., Katzourakis, A., 2023. The evolution of SARS-CoV-2. Nat. Rev. Microbiol. 21 (6), 361–379. Number: 6 Publisher: Nature Publishing Group. https://doi.org/10.1038/s41579-023-00878-2

Monge, S., Pastor-Barriuso, R., Hernán, M.A., 2023. The imprinting effect of covid-19 vaccines: an expected selection bias in observational studies. BMJ 381, e074404. Publisher: British Medical Journal Publishing Group Section: Research Methods & Reporting. https://doi.org/10.1136/bmj-2022-074404

Müller, K., 2020. Here: a simpler way to find Your Files. R package version 1.0.1. https://CRAN.R-project.org/package=here.

Otto, S.P., Day, T., Arino, J., Colijn, C., Dushoff, J., Li, M., Mechai, S., Domselaar, G.V., Wu, J., Earn, D. J.D., Ogden, N.H., 2021. The origins and potential future of SARS-CoV-2 variants of concern in the evolving COVID-19 pandemic. Curr. Biol. 31 (14), R918–R929. Publisher: Elsevier. https://doi.org/10.1016/j.cub.2021.06.049
Park, S.W., Bolker, B.M., Funk, S., Metcalf, C. J.E., Weitz, J.S., Grenfell, B.T., Dushoff, J., 2022. The importance of the generation interval in investigating dynamics and control of new SARS-CoV-2 variants. J. Royal Soc. Interf. 19 (191), 20220173. Publisher: Royal Society. https://doi.org/10.1098/rsif.2022.0173

Pearson, C. A.B., Silal, S.P., Li, M. W.Z., Dushoff, J., Bolker, B.M., Abbott, S., Schalkwyk, C.V., Davies, N.G., Barnard, R.C., Edmunds, W.J., Bingham, J., Meyer-Rath, G., Jamieson, L., Glass, A., Wolter, N., Govender, N., Stevens, W.S., Scott, L., Misana, K., Moultrie, H., Pulliam, J. R.C., 2021. Bounding the levels of transmissibility & immune evasion of the Omicron variant in South Africa, https://doi.org/10.1101/2021.12.19. 21268038

Powell, A.A., Kirsebom, F., Stowe, J., Ramsay, M.E., Lopez-Bernal, J., Andrews, N., Ladhani, S.N., 2022. Protection against symptomatic infection with delta (B.1.617.2) and omicron (B.1.1.529) BA.1 and BA.2 SARS-CoV-2 variants after previous infection and vaccination in adolescents in England, August, 2021–March, 2022: a national, observational, test-negative, case-control study. Lancet Infect. Dis. 0 (0). Publisher: Elsevier. https://doi.org/10.1016/S1473-3099(22)00729-0

Price, G.R., 1970. Selection and covariance. Nature 227 (5257), 520–521. Number: 5257
Publisher: Nature Publishing Group. https://doi.org/10.1038/227520a0

Reyné, B., Richard, Q., Selinger, C., Sofonea, M.T., Djidjou-Demasse, R., Alizon, S., 2022. Non-Markovian modelling highlights the importance of age structure on Covid-19 epidemiological dynamics. Math. Model. Nat. Phenom. 17, 7. Publisher: EDP Sciences. https://doi.org/10.1051/mmnpp/2022008

Roff, D.A., 2008. Defining fitness in evolutionary models. J. Genet. 87 (4).

Saltelli, A., Ratto, M., Andres, T., Campolongo, F., Cariboni, J., Gatelli, D., Saisana, M., Tarantola, S., 2008. Global Sensitivity Analysis: The Primer. John Wiley & Sons.

Sasaki, A., Lion, S., Boots, M., 2022. Antigenic escape selects for the evolution of higher pathogen transmission and virulence. Nat. Ecol. Evol. 6 (1), 51–62. Number: 1 Publisher: Nature Publishing Group. https://doi.org/10.1038/s41559-021-01603-z

Schmid-Hempel, P., 2008. Parasite immune evasion: a momentous molecular war. Trends Ecol. Evol. 23 (6), 318–326. Publisher: Elsevier. https://doi.org/10.1016/j.tree.2008.02.011

Sofonea, M.T., Reyné, B., Elie, B., Djidjou-Demasse, R., Selinger, C., Michalakis, Y., Alizon, S., 2021. Memory is key in capturing COVID-19 epidemiological dynamics. Epidemics 35, 100459. https://doi.org/10.1016/j.epidem.2021.100459

Starr, T.N., Greaney, A.J., Stewart, C.M., Walls, A.C., Hannon, W.W., Veesler, D., Bloom, J.D., 2022. Deep mutational scans for ACE2 binding, RBD expression, and antibody escape in the SARS-CoV-2 Omicron BA.1 and BA.2 receptor-binding domains. PLOS Pathogens 18 (11), e1010951. Publisher: Public Library of Science. https://doi.org/10.1371/journal.ppat.1010951

Stein, C., Nassereldine, H., Sorensen, R. J.D., Amlag, J.O., Bisignano, C., Byrne, S., Castro, E., Coberly, K., Collins, J.K., Dalos, J., Daoud, F., Deen, A., Gakidou, E., Giles, J.R., Hulland, E.N., Huntley, B.M., Kinzel, K.E., Lozano, R., Mokdad, A.H., Pham, T., Pigott, D.M., Reiner, J. R.C., Vos, T., Hay, S.I., Murray, C. J.L., Lim, S.S., 2023. Past SARS-CoV-2 infection protection against re-infection: a systematic review and meta-analysis. Lancet 0 (0). Publisher: Elsevier. https://doi.org/10.1016/S0140-6736(22)02465-5

Team, R.C., 2023. R: a language and environment for statistical computing. R Foundation for Statistical Computing. Vienna, Austria. https://www.R-project.org/.

UKHSA, 2022. SARS-CoV-2 variants of concern and variants under investigation report 36. https://www.gov.uk/government/publications/investigation-of-sars-cov-2-variants-technical-briefings.

Volz, E., Hill, V., McCrone, J.T., Price, A., Jorgensen, D., OToole, A., Southgate, J., Johnson, R., Jackson, B., Nascimento, F.F., Rey, S.M., Nicholls, S.M., Colquhoun, R.M., da Silva Filipe, A., Shepherd, J., Pascall, D.J., Shah, R., Jesudason, N., Li, K., Jarrett, R., Pacchiarini, N., Bull, M., Geidelberg, L., Siveroni, I., Koshy, C., Wise, E., Cortes, N., Lynch, J., Kidd, S., Mori, M., Fairley, D.J., Curran, T., McKenna, J.P., Adams, H., Fraser, C., Golubchik, T., Bonsall, D., Moore, C., Caddy, S.L., Khokhar, F.A., Wantoch, M., Reynolds, N., Warne, B., Maksimovic, J., Spellman, K., Patel, G., Garcia-Casado, M.V., Dibley, T., McGuigan, S., Rogers, H.A., Hunter, A.D., Souster, E., Neaverson, A.S., Goodfellow, I., Loman, N.J., Pybus, O.G., Robertson, D.L., Thomson, E.C., Rambaut, A., Connor, T.R., 2021. Evaluating the effects of SARS-CoV-2 spike mutation D614G on transmissibility and pathogenicity. Cell 184 (1), 64–75. https://doi.org/10.1016/j.cell.2020.11.020

Vossen, M.T., Westerhout, E.M., Söderberg-Nauclér, C., Wiertz, E.J., 2002. Viral immune evasion: a masterpiece of evolution. Immunogenetics 54 (8), 527–542. https://doi.org/10.1007/s00251-002-0493-1

White, L.J., Medley, G.F., 1998. Microparasite population dynamics and continuous immunity. Proc. Royal Soc. London. Ser. B 265 (1409), 1977–1983. Publisher: Royal Society. https://doi.org/10.1098/rspb.1998.0528

WHO, 2023. Updated working definitions and primary actions for SARS-CoV2 variants. https://www.who.int/publications/m/item/updated-working-definitions-and-primary-actions-for-sars-cov-2-variants.

Wickham, H., 2016. ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org.

Wilke, C.O., 2024. cowplot: streamlined plot theme and plot annotations for 'ggplot2'. R package version 1.1.3. https://CRAN.R-project.org/package=cowplot.

Wymant, C., Bezemer, D., Blanquart, F., Ferretti, L., Gall, A., Hall, M., Golubchik, T., Bakker, M., Ong, S.H., Zhao, L., Bonsall, D., de Cesare, M., MacIntyre-Cockett, G., Abeler-Dörner, L., Albert, J., Bannert, N., Fellay, J., Grabowski, M.K., Gunsenheimer-Bartmeyer, B., Günthard, H.F., Kivelä, P., Kouyos, R.D., Laeyendecker, O., Meyer, L., Porter, K., Ristola, M., van Sighem, A., Berkhout, B., Kellam, P., Cornelissen, M., Reiss, P., Fraser, C., the Netherlands ATHENA HIV Observational Cohort, THE BEEHIVE COLLABORATION, 2022. A highly virulent variant of HIV-1 circulating in the Netherlands. Science 375 (6580), 540–545. Publisher: American Association for the Advancement of Science. https://doi.org/10.1126/science.abk1688

Zachreson, C., Tobin, R., Szanyi, J., Walker, C., Cromer, D., Shearer, F.M., Conway, E., Ryan, G., Cheng, A., McCaw, J.M., Geard, N., 2023. Individual variation in vaccine immune response can produce bimodal distributions of protection. Vaccine 41 (45), 6630–6636. https://doi.org/10.1016/j.vaccine.2023.09.025