# Eco-evolutionary dynamics of adapting pathogens and host immunity

Pierre Barrat-Charlaix, Richard A Neher

► To cite this version:

Pierre Barrat-Charlaix, Richard A Neher. Eco-evolutionary dynamics of adapting pathogens and host immunity. eLife, 2024, 10.7554/eLife.97350.2. hal-04876052

HAL Id: hal-04876052

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

Submitted on 9 Jan 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>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Eco-evolutionary dynamics of adapting pathogens and host immunity

Pierre Barrat-Charlaix, Richard A Neher

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

Biozentrum, Universität Basel, Switzerland • Swiss Institute of Bioinformatics, Switzerland • DI SAT, Politecnico di Torino, Italy

**Reviewed Preprint**

v2 • October 18, 2024

Revised by authors

**Reviewed Preprint**

v1 • July 3, 2024

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

https://en.wikipedia.org/wiki/Open_access

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

Copyright information

## Abstract

As pathogens spread in a population of hosts, immunity is built up and the pool of susceptible individuals is depleted. This generates selective pressure, to which many human RNA viruses, such as influenza virus or SARS-CoV-2, respond with rapid antigenic evolution and frequent emergence of immune evasive variants. However, the host's immune systems adapt and older immune responses wane, such that escape variants only enjoy a growth advantage for a limited time. If variant growth dynamics and reshaping of host-immunity operate on comparable time scales, viral adaptation is determined by eco-evolutionary interactions that are not captured by models of rapid evolution in a fixed environment. Here, we use a Susceptible/Infected model to describe the interaction between an evolving viral population in a dynamic but immunologically diverse host population. We show that depending on strain cross-immunity, heterogeneity of the host population, and durability of immune responses, escape variants initially grow exponentially, but lose their growth advantage before reaching high frequencies. Their subsequent dynamics follows an anomalous random walk determined by future escape variants and results in variant trajectories that are unpredictable. This model can explain the apparent contradiction between the clearly adaptive nature of antigenic evolution and the quasi-neutral dynamics of high frequency variants observed for influenza viruses.

## eLife Assessment

This **important** study provides a new perspective on how human immunity shapes the antigenic evolution of pathogens. By combining theory and simulation the authors make a **compelling** case for the importance of eco-evolutionary interactions in population-level virus-host dynamics, which arise due to coupling between the dynamics of immune memories and viral variants. Although the work does not propose improved data-driven viral forecasting methods, it makes a conceptual contribution that advances the field's understanding of this problem's intrinsic difficulty.

https://doi.org/10.7554/eLife.97350.2sa3
Introduction

Many human RNA viruses adapt rapidly to evade pre-existing immunity and re-infect humans
multiple times over their lifetime. The most prominent examples of this evolution are influenza
virus and SARS-CoV-2 [1, 2], for which the changing virus population is surveilled in great
detail and vaccines are updated regularly. To improve the match between the virus population
and the vaccine, several groups are working on predictive models to anticipate the variants that
dominate future viral populations [3, 4],

A common framework to model the rapid evolutionary dynamics of RNA viruses is to consider a
population located away from the fitness optimum and with many accessible beneficial mutations
[5]. In this setting, clones compete to accumulate beneficial mutations as quickly as possible. In
a process called selective sweep, successful variants expand and tend to be the ancestors of the
future population while less successful mutants eventually disappear. The resulting fitness
distribution is a wave traveling along the fitness axis, the so-called *traveling fitness wave* [6–
8]. As the pathogen circulates, hosts develop immunity which leads to a ‘deterioration of the
environment’ for the pathogen which approximately balances the increase in average fitness due
to adaptation.

The traveling wave framework has been extensively used in this context as it allows for a
straightforward way to approach the prediction problem: each variant is assumed to have a fixed
fitness relative to other variants, and inferring the fitness of all competing variants should allow
prediction of the future composition of the population. Indeed, current methods typically infer the
instantaneous growth advantage of a strain based on past and present circulation and then project
this growth advantage forward in time [9–11], While future mutations can reshuffle the
relative fitness of lineages and thereby limit predictability, in these models a lineage that is most
fit at any given time is most likely to dominate in the long run.

One short-coming of the traveling wave approach is the lack of explicit representation of the
epidemiological dynamics and of the host's immunity. Indeed, fitness is only an effective
parameter that summarizes the complex interplay between viral antigenic properties and the
hosts' immune systems. As such, it cannot explicitly describe important phenomena such as the
build-up of host immunity to new variants, variant specific immunity, or the interaction between
strains through antigenic cross-reactivity. Taking the hosts' immunity and viral cross-immunity
into account has the potential to strongly improve predictions [4] or explain why prediction is
difficult [12],

The interaction between epidemiological dynamics and hosts' immunity are often modeled using
generalizations of the Susceptible-Infected-Recovered model (SIR) to include multiple viral strains
[13, 14]. In this setting, the natural question is that of the ultimate fate of the pathogen: will it
go extinct, diversify to the point of speciation, or reach the so-called Red Queen State where new
strains continuously replace old ones [15–18]. To remain tractable, these studies typically
approximate population immunity as a low-dimensional landscape in which the viral population
evolves and ignore the complex heterogeneity in the immunity of different individuals.
Furthermore, immunity is often assumed to be long-lived and evolution of the pathogen in a stable
low dimensional landscape gives rise to traveling waves.

Here, we study how novel variants of a virus shape the host population's immunity, which in turn
changes their own growth dynamics. To do so, we use a multi-strain SIR model combining immune
waning and heterogeneous immunity of the hosts. Such heterogeneity has been demonstrated for
influenza virus in individuals of different ages [19, 20]. We show that this model generically
leads to a situation where novel immune evasive variants emerge. In a homogeneous population
of hosts, this leads to a succession of selective sweeps where novel variants compete against each other and replace previously circulating variants. However, in a heterogeneous population with more rapid waning of immunity, initially growing variants lose their selective advantage before reaching fixation due to immunological adjustment of the host population. The phenomenology of our epidemiological model is reminiscent of ecological systems such as consumer resource models, where adaptation by one species shifts the global equilibrium and distribution of other species but does not necessarily result in a selective sweep [21], In these systems, adaptation can usually not be modeled by a fixed fitness parameter for each strain but rather depends on the composition of the population [22].

Strain dynamics in our model differ qualitatively from what is expected in the traveling wave scenario. While adaptive mutations are highly overrepresented in genetic diversity, they cease having a growth advantage when reaching intermediate frequencies, a process we call “expiring fitness”. Once the fitness effect of a mutation has expired, its frequency randomly changes up or down as subsequent adaptive mutations occur on the same or on different genomic backgrounds.

This resemblance to neutral evolution could have important consequences for predictability of viral evolution. It is interesting to relate this to the recent observations that the evolution of influenza is not as predictable as one would expect from typical models [11, 12], In particular, we observed in [12] that the frequency trajectories of mutants of A/H3N2 influenza show features that are expected in neutral evolution but hard to explain in a traveling wave framework.

# Results

## 1. Multi-strain SIR model

We describe the interaction of several viral strains and host immunity using a Susceptible/Infected compartmental model, similar to those used in [14, 15]. In the most general form, the model describes N variants of the virus labeled a ∈ {1... N} circulating among M groups of hosts with distinct exposure histories labeled i ∈ {1... M} (immune groups). These groups could be different age cohorts or could be geographically separated. For each group i, we define compartments $I_i$ and $S_i$ as respectively the number of individuals of this group infected or susceptible to strain a. We assume that the total population of hosts is 1 so that we always have $0 \le I_i \le 1$, $S_i \le 1$, and values of $I_i$ and $S_i$ can be interpreted as fractions of the host population.

As with usual compartmental models, we assume that the dynamics are driven by the interaction of susceptible and infected hosts, leading to infections and gains of immunity. The rate at which hosts of group i susceptible to variant a are infected by a is $\alpha \sum_{j=1}^{M} c_{ij} I_j$. Here, $\alpha$ is an overall infection rate while $c_{ij}$ represents the probability of an encounter between individuals of group i and j. Thus, the above rate takes into account infections with strain a caused by hosts of all groups. Considering that infected hosts recover at rate $\delta$, we can thus write the dynamics for $I_i$:

$$ \dot{I}_i = \alpha S_i^a \sum_{j=1}^{M} C_{ij} I_j^a - \delta I_i^a. \qquad (1) $$

When a host of group i is infected by strain b, it gains immunity against the infecting strain b, but also to other strains a with probability $0 \le \kappa_i^a \le 1$. Thus, $S_i^a$ decreases at a rate proportional to $\kappa_i^a$ and to the number of hosts infected by b for every strain b. Since susceptibles to a are depleted by infections from other strains, the dynamics of all strains are coupled. This coupling is determined
by the matrices $K_i$ of dimension $N \times N$, which in general differ between groups $i$ with different prior exposure history. Additionally, waning of immunity at a rate $\gamma$ causes immune hosts to re-enter the susceptible compartment. We can now write the dynamics of $S_i^\gamma$ as

$$ \dot{S}_i^a = -\alpha \sum_{b=1}^{N} \sum_{j=1}^{M} S_i^a K_i^{ab} C_{ij} I_j^b + \gamma(1 - S_i^a), \quad (2) $$

where the first term accounts for immunity gains (or loss of susceptibility) due to infections or cross-immunity while the second represents immune waning. This model introduced by Gog and Grenfell [14] assumes that immunity builds up through exposure and not only through infection. This explains that the change in $S_i^\gamma$ is simply proportional to $S_i^\gamma \cdot I_i$, regardless of the susceptibility of hosts to strain $b$. Alternative models that require infection for acquisition of immunity have qualitatively similar dynamics, but are mathematically more complex (see supplementary material). We also represented loss of susceptibility to $a$ due to exposure to $a$ using a trivial cross-immunity term $\kappa_i^{ab} = 1$.

An important property of our model is that the probability to generate cross-immunity can differ between groups. The motivation is that strains $a$ and $b$ may be perceived as antigenically different by some immune systems, leading to a low $\kappa_i^{ab}$, but as highly similar by others, leading to $\kappa_i^{ab} \approx 1$. Such a heterogeneous response by different immune systems has been observed experimentally in the case of influenza: in [19, 20] for instance, it was found that a given mutation in an influenza strain may allow escape from the antibodies of some individuals, i.e. low $\kappa_i^{ab}$. while it had little effect for the serum of other individuals, i.e. high $\kappa_i^{ab}$. Heterogeneous immune response could be caused by varying histories of strain exposure for different individuals, for instance due to difference in age or geographical region. If immune groups correspond to age cohorts, mixing between groups is rapid, and we can simplify the connectivity between groups to $C_{ij} = 1/M$. If immune groups are shaped by geographic differences in exposure, the connectivity would be close to 1 on the diagonal ($1 - C_{ii} \ll 1$) while off-diagonal terms would be small ($C_{ij} \ll 1$ for $i \neq j$).

## 2. Invasion of an adaptive variant

Hosts' immune heterogeneity and strain cross-immunity play two different roles in the model. The latter allows the model to reach a non-trivial equilibrium where multiple strains co-exist, while the former affects the convergence to the equilibrium.

To illustrate this, we design a simple scenario with only two strains: a wild-type and a variant. Accordingly, indices $(a,b)$ describing strains will take values $\{wt,v\}$. We consider that the two strains share the same infectivity rate $\alpha$, which amounts to say that they would have the same reproductive rate in a fully naive population. The case where the two strains differ in intrinsic fitness is explored in detail in the SI. In brief, as long as the difference in intrinsic fitness is not too large compared to cross-immunity effects, the qualitative results given below hold, while larger intrinsic fitness differences lead to more classical dynamics like selective sweeps.

In the first version of this scenario, we will only consider one immune group, that is $M = 1$. We can thus skip the indices $i, j \in \{1... M\}$, and we only have one cross-immunity matrix $K$ that we parametrize as

$$ K = \begin{bmatrix} 1 & b \\ f & 1 \end{bmatrix}, \quad (3) $$
with $0 \le b, f \le 1$. b quantifies the amount of “backward”-immunity to the wild-type caused by the variant: a large b means that it is likely that an infection by the variant causes immunity to the wild-type. Conversely, f quantifies the “forward”-immunity: infections by the wild-type causing immunity to the variant. If $f = b = 1$, the two strains are antigenically indistinguishable, and thus essentially identical for the model. Conversely, if $f = b = 0$, the two strains are completely distinct and do not interact through cross-immunity.

The dynamical equations now take a simplified form:

$$ \dot{S}^a = -\alpha S^a \sum_{b \in \{wt,v\}} K^{ab} I^b + \gamma(1 - S^a), \quad (4) $$

$$ \dot{I}^a = (\alpha S^a - \delta) I^a. $$

We can immediately derive the equilibrium state for this simplified case. We first define the reproductive number of strain a as $R^a = \alpha S^a/\delta$: $I^a$ grows when $R^a > 1$ and declines when $R^a < 1$. The equilibrium susceptibility is therefore $S^a = \delta/\alpha$, such that $R^a = 1$. On the other hand, the equilibrium prevalence is determined by the inverse of the cross-immunity matrix K:

$$ I_{eq}^a = \frac{\gamma}{\delta} \left(1 - \frac{\delta}{\alpha}\right) \left(\mathbf{K}^{-1}\vec{1}\right)_a, \quad (5) $$

with $\vec{1}$ being the vector [1; 1], The order of magnitude of the prevalence is given by the ratio of the rate of waning $\gamma(1 - \delta/\alpha)$ and the recovery rate $\delta$. In the following, we frequently use values $\alpha = 3$ and $\gamma = 5 \cdot 10^{-3}$ in units of inverse generations $\delta$, i.e. we set $\delta = 1$. At equilibrium, this corresponds to a fraction ~0.003 of the host population being infected at any time. If generation time is a week, which is roughly the case for respiratory viruses such as influenza virus or SARS-CoV-2, the fraction of hosts infected in any year is ~ 0.15, which is of similar magnitude as empirical estimates for influenza [23].

It is also straightforward to compute the fraction of infections caused by the variant at equilibrium, thereafter referred to as the frequency of the variant. We find that this frequency is

$$ \beta = \frac{1 - f}{(1 - b) + (1 - f)}. \quad (6) $$

In the case where $b = f$, the variant will ultimately settle at frequency 1/2. This includes the case where $b = f = 0$, where the two strains are completely independent and do not interact. On the contrary if $b \ne f$, the final frequency of the variant can in principle be anywhere between 0 and 1. For example if $b > f$, the variant is more likely to cause immunity to the wild-type than the wild-type is to cause immunity to the variant. In this case, $\beta > 1/2$ and the variant will be the dominant variant.

We are primarily interested in an “invasion” scenario where only the wild-type is initially present in the population, that is $I^v = 0$ at $t < 0$. Cross-immunity with the resident strain reduces the fraction of hosts susceptible to the variant below one even though it has not circulated yet. But the number of susceptible hosts is always larger than the equilibrium value $\delta/\alpha$ unless $f = 1$, As a result, the growth rate of the variant is initially positive and Aven bv

$$ s(t = 0) = \frac{(1 - f)(\alpha - \delta)}{\delta + f(\alpha - \delta)} \quad (7) $$
The variant thus increases initially exponentially until it has become sufficiently frequent that it starts having a substantial effect on the immunity landscape, before eventually settling into an equilibrium with the wild-type. The details of the equilibrium reached by the system in the absence of additional mutant variants is given in the SI. **Figure 1** explores different scenarios numerically.

The top row of **Figure 1** shows the dynamics of the model after the introduction of the variant in a homogeneous population ($M = 1$). As expected, the number of infections by the variant initially rises while the number of susceptibles $S^v$ decreases. However, as $S^v$ goes below the critical value $\delta/\alpha$, $I^v$ starts to decline and then oscillates around the equilibrium value before finally converging to it. The mathematical properties of these oscillations are discussed in the SI.

However, these strong and slowly damped oscillations are not what is observed in circulating viruses. For instance, in the first oscillation in the specific example of **Figure 1** , the prevalence of the wild-type $I^{wt}$ goes down to microscopic levels and the frequency of the variant approaches one, see **Figure 1** . During stochastic circulation in a finite population of hosts the wild-type would likely be lost. The theoretical equilibrium that is reached at long times is thus not very relevant, and what would be actually observed in reality is a selective sweep by the variant.

Oscillations are the consequence of the rapid rise of the variant followed by an overshoot. This effect is mitigated by immunological heterogeneity, as shown in the following example with $M = 10$ groups. For group $i = 1$, the cross-immunity matrix $K_1$ takes the same form as in the previous scenario, given by **Equation 3** . However, for other groups, we assume that the two strains are virtually identical, with the cross-immunity having the form

$$K_{i>1} = \begin{bmatrix} 1 & 1-\varepsilon \\ 1-\varepsilon & 1 \end{bmatrix}, \qquad (8)$$

where $\varepsilon \ll 1$. Our reasoning is that we expect an adaptive variant to escape the existing immunity for part of the host population, here immune group 1, while having little effect on the rest of the hosts.

One consequence of many groups that are indifferent to the variant is that globally the excess susceptibility to the variant is lower. If mixing is rapid, the initial growth rate of the variant is smaller by a factor $M$ compared to the one-group case. If mixing is slow, the initial growth of the variant is as rapid as in the one-group case, but then spreads only slowly across groups. Globally, the frequency of the variant thus never reaches values close to one and population wide oscillations are reduced.

The central and bottom rows of **Figure 1** show the result of the invasion for $M$ groups respectively for the rapid and slow mixing cases. In both scenarios, the initial number of hosts susceptible to the variant is now closer to $\delta/\alpha$. When mixing is fast, the frequency of the variant initially resembles a standard selective sweep (dashed line in **Figure 1** ) before saturating, while dynamics are more complicated for the slow mixing case. Either way, the main effect of the immune groups is that the overshoot past the equilibrium is much smaller and dampening of the oscillations stronger. As a result, the frequency of the variant approaches its equilibrium value without effectively sweeping to fixation before.

Notably, the equilibrium frequency in the above examples does not depend on $M$ and **Equation 6** remains valid for $\varepsilon = 0$. This invariance is a consequence of the fact that for $\varepsilon = 0$, the variant and wild-type strains are completely equivalent in immune groups $i > 1$ and equilibrium is only
<figure><img src="image_6.png" /><figcaption>FIG. 1.</figcaption></figure>

Invasion scenario for **Top row:** one immune group, **Middle row:** $M = 10$ immune groups and fast mixing $C_{ij} = 1/M$ and **Bottom row** $M = 10$ immune groups and slow mixing $C_{ij} = 1/10M$. Other parameters are the same for all rows: in units of $\delta$. We set $\alpha = 3$ and $\gamma = 5 \cdot 10^{-3}$, and $f = 0.65$, $b = 0.8$, $\varepsilon = 0.01$. For both rows, graphs represent: **Left:** number of hosts infectious with the wild-type and the variant; **Middle:** number of hosts susceptible to the wild-type and the variant, with the equilibrium value $\delta/\alpha$ as a gray line; **Right:** fraction of the infections due to the variant. The thick gray line shows the expected equilibrium frequency $\beta$ in the case with one immune group, given in **Eq. 6**. The dashed line shows the trajectory of a constant fitness logistic growth with the same initial growth rate.
determined by cross-immunity in group *i* = 1 (see SI). For small *ε* the equilibrium shifts slightly, but **Equation 6** ☐ remains a good approximation.

While this simple two-strain model predicts that the two strains come to an equilibrium at frequency *β*, their frequency will of course continue to change due to emergence of additional strains, which we will discuss below.

Even though the variant has a clear growth rate advantage when it appears, this does not result in it replacing the wild-type. This contrasts with the typical “selective sweep” that occurs when the growth rate advantage stays constant, which is illustrated as a dashed line in the figure. We refer to frequency trajectories of a variant that at first rise exponentially before settling at an intermediate frequency as *partial sweeps*. As we will discuss below, such partial sweeps can lead to qualitatively different evolutionary dynamics and its predictability.

If the initial growth is due to higher susceptibility, it is misleading to think of it as an intrinsic fitness advantage of the variant. Instead, the initial growth is the result of an imbalance in the immune state of the host population, which gradually disappears as the variant becomes more frequent, as shown in the central panels of **Figure 1** ☐. In this sense, our model is comparable to ecological systems where interaction between organisms cannot be fully explained using a fixed scalar fitness for each strain but rather depends on the composition of the viral and host population. In particular, the stalling of frequency increase giving rise to the partial sweep is reminiscent of consumer resource models [21 ☐, 22 ☐], highlighting the link between ecological and epidemiological models. An important consequence of these dynamics is that predicting the equilibrium frequency reached by the variant and its ultimate fate is hard from the observation of a partial frequency trajectory.

**3. Ultimate fate of the invading variant**

In the invasion scenario discussed above, dynamics stop after the initial variant has reached an equilibrium frequency. However, as the viral population evolves, new adaptive mutants can appear. In the framework of the SIR model, a new strain translates into extending the cross-immunity matrix by one row and one column. Each new variant will perform its own partial sweep, and saturate at frequencies $β_2, β_3, ...$ sampled from some distribution $P_β$. This process is shown in panel A of **Figure 2** ☐, using the SIR model to simulate up to $N = 7$ variants. For the sake of illustration, it shows a simple scenario where the initial variant appears at time 0 in a homogeneous wild-type population, and subsequent mutants appear at regular time intervals. Simulations are performed using $M = 10$ immune groups, resulting in a slight overshoot of the equilibrium frequency for each trajectory.

Here, we focus on the mutation or set of mutations *A* that defines the initial variant. The initial growth rate advantage given by *A* eventually disappears, meaning that after some time we can consider it as neutral. As subsequent mutants appear, they either do so on the background of the wild-type, in which case they do not carry *A*, or on the background of the initial variant in which case they do carry *A*. If we suppose that recombination is negligible, the frequency of *A* increases or decreases as each new variant undergoes its own partial sweep. This process is shown in panel **B** of **Figure 2** ☐, with shades of red (resp. blue) indicating a variant carrying *A* (resp. not carrying *A*). The thick line in between the red and blue surfaces indicate the frequency at which mutation *A* is found, and in practice moves up and down randomly.

The scenario illustrated in **Fig. 2** ☐ suggests that many aspects of the variant dynamics can be approximated by a simple abstraction: If *x* is the frequency of a mutation *A*, a new variant has a probability *x* to appear on the background of *A* and thus carry *A*, and a probability $1 - x$ to not carry *A*. If new mutants emerged well separated in time with rate *p*, meaning that they reach equilibrium before the next variant emerges, and if new variants have a similar cross-immunity
<figure><img src="image_7.png" /><figcaption>FIG. 2.</figcaption></figure>

**A:** Simulation of SIR Equations 1 & 2 with additional strains appearing at regular time intervals. The fraction of infections (frequency) caused by each strain is shown as a function of time. The first strain to appear at $t = 0$ is the variant of interest, and curves are shown in shades of red if they appear on the background of this variant, and of blue if they appear on the background of the wild-type. **B:** Same as **A** but with frequencies stacked vertically. The black line delimiting the red and blue areas represents the frequency at which the mutations defining the original variant are found. **C:** Three realizations of the random walk of Equation 9, all starting at $x = 0.5$. Two instances converge rapidly to frequency 0 and 1, corresponding to apparent selective sweeps, while the remaining one oscillates for a longer time. **D:** Representation of a partial sweep using the expiring fitness parametrization of Equation 11. The frequency $x$ of the variant is shown as a blue line saturating at value $\beta$ (gray line). The thin dashed line shows a selective sweep with constant fitness advantage $s_0$. The fitness $s$ is a red dashed line, using the right-axis.
with all existing variants (see SI), the dynamics of $x(t)$ are described by a particular random walk: in each time interval $dt$, a partial sweep of amplitude $\beta$ occurs with probability $\rho dt \cdot P_{\beta}(\beta)$, changing $x$ in the following way:

$$x(t + dt) = x(t) + \Delta x, \text{ where } \Delta x = \begin{cases} -\beta x & \text{with probability } (1-x), \\ \beta(1-x) & \text{with probability } x. \end{cases} \qquad (9)$$

For example, if a new mutant appearing on the background of A does a partial sweeps of amplitude $\beta$, the frequency of A among the fraction of strains $(1 - \beta)$ not concerned by the sweep will still be $x$, and its frequency among the fraction $\beta$ of strains concerned by the sweep will be 1. Overall, this gives a frequency change of $\Delta x = (1 - x)\beta$. A similar reasoning gives us the frequency change when the new mutant appears on the wild-type background. Finally, if no sweep occurs in the time interval $dt$, that is with probability $1 - \rho dt$, $x$ remains unchanged. The resulting frequency dynamics of mutations has many similarities to the effect of 'genetic draft', that is the frequency dynamics of neutral mutations due to linked selective sweeps [24],

Examples of trajectories from the random walk are shown on panel C of Figure 2, all initially starting around $x_0 \approx 1/2$. Two trajectories converge monotonically to 0 and 1. This is a consequence of one interesting property of Equation 9: the probability for $\Delta x$ to be positive increases with $x$, but the magnitude of the upwards steps decreases as $1 - x$, and symmetrically with downward steps. This leads to trajectories leading almost exponentially to 0 and 1: it can in fact be shown that trajectories that always go downwards or upwards represent a finite and relatively large fraction of all possible trajectories (see SI). On the other hand, steps away from the closest boundary are unlikely but much larger, resulting in 'jack-pot' events [25]. This can be seen for the blue trajectory in Figure 2, which oscillates for a longer time.

It is also interesting to look at the moments of the step size $\Delta x$. The first two are easily computed, and we find

$$\begin{align} \langle \Delta x \rangle &= 0 \tag{10} \\ \langle \Delta x^2 \rangle &= \rho \langle \beta^2 \rangle P_{\beta} x (1-x). \nonumber \end{align}$$

The first moment being 0 means that for the random walk, increasingly probable but small steps towards the closest boundary (0 or 1) are exactly compensated by rarer but larger steps away from the boundary. Importantly, this means that on average, the trajectory of mutation A is not biased towards either fixation or loss, regardless of the frequency that the initial partial sweep brought it to. For instance, a mutation seen at frequency $x_0$ should on average stay at this frequency, which means in practice that in a finite population it has a chance to reach 1 and fix, and a chance $1 - x_0$ to reach 0 and vanish.

On the other hand, the second moment resembles neutral drift [26]: in neutral evolution, allele frequency also undergoes a zero-average random walk with the second moment having the form $x(1-x)/N$ with $N$ being the population size. Therefore, this model would predict an “effective population size” as $N^{-1} = \rho(\beta^2)_{\rho}$ completely independent of the size of the viral population. However, there are important differences to neutral drift: in neutral evolution, higher moments of order $k > 2$ decay as $N^{1-k}$ and are thus negligible in large populations, whereas here they are independent of $N$ and scale as $(\rho\beta)^k$. Depending on higher moments of $P_{\beta}$, allele dynamics will deviate qualitatively from neutral behavior.
## 4. Abstraction as ‘expiring’ fitness advantage

In general, the dynamics of the SIR model proposed in Equations 1 & 2 depend on the interactions between N strains through an N × N cross-immunity matrix. While this model is useful to give a mechanistic explanation of partial sweeps, it is in general impractical to analyze and numerically simulate for a large number of variants. The random walk model introduced above is simple to analyze and simulate, but assumes that variants rise to their equilibrium frequency instantaneously.

To explore the consequences of partial sweeps over broader parameter ranges, we propose an empirical model that has the same qualitative properties as the overdamped SIR, namely a growth rate that decreases as a strain becomes more frequent and partial sweep trajectories, but is simpler to analyze and simulate on a large scale. In this effective model, the growth rate s of the variant is not explicitly set by the susceptibility dynamics in the host population, but instead decays at a rate proportional to the frequency x of the variant:

$$ \dot{x} = sx(1-x) \quad \text{and} \quad \dot{s} = -\nu xs. \qquad (11) $$

The dynamic of x in the first equation is simply given by the usual logistic growth with fitness s. To mimic increasing immunity against the invading variant, the growth advantage s decreases proportionally to the abundance of the variant (second part of Eq. 11). The initial value of $s_0$ is connected to the invasion rate of the SIR models given in Eq. A26 of the SI.

The dynamics of this new model are represented in panel D of Figure 2, with an initial frequency $x_0 \ll 1$ and an initial growth rate $s_0 = 0.05$. The initial growth of x is identical to a classical selective sweep of fitness $s_0$ (represented as a dashed line). However, its fitness advantage gradually “expires”, as shown by the red line in the figure. As the variant progressively “runs out of steam”, its frequency finally saturates at a value $\beta$ given by (see SI)

$$ \beta = 1 - e^{-s_0/\nu}. \qquad (12) $$

This final value $\beta$ depends only on the ratio between the initial fitness advantage $s_0$ and the rate of fitness decay $\nu$. For a large enough $s_0$, $\beta$ can be arbitrarily close to 1, meaning that this model still accommodates for full selective sweeps as a special case. In the general case, x reaches its final value $\beta < 1$ and remains there forever unless other variants appear.

It is important to state that the main aim of this effective model is to qualitatively reproduce the phenomenology of the SIR, and in particular the partial sweeps, while being easier to simulate. It recapitulates the salient feature of invading immune evasive variants: (i) initial exponential growth, and (ii) eventual saturation at an intermediate frequency. We can thus use it to analyze the long term consequences of the random walk dynamics of Figure 2. However, we do not expect the frequency of the variant x to have quantitatively equivalent dynamics in the two models. In particular, due to its simplicity, this model does not show the complex oscillatory behavior of the SIR model. Section (A 9) of the SI discusses in more detail the links between the parameters of the two models and the fundamental differences. While we can express the rate $\nu$ at which the growth rate declines in terms of the parameters of the simplest SIR models, for models with many groups or with oscillatory dynamics, the decay rate of the growth advantage should be interpreted as an effective parameter that captures a generic effect of reduced growth with increasing circulation.
# 5. Consequences for predictability and population dynamics

Accurate prediction of dominant viral variants of the future could improve the choice of antigens in vaccines against rapidly evolving viruses. Specifically, if a potentially adaptive mutation is observed in a viral population, one would want to know if the corresponding variants will grow in frequency, and if yes to what point? The typical traveling wave framework would predict that fast-growing variants should keep on growing until an even fitter one appears. This way of thinking about the prediction problem has shown mixed results. In the case of A/H3N2 influenza for instance, we showed that there are few signatures that suggest fit variants grow in frequency consistently [12](Figure 3),

In **Figure 3** (Figure 3), we reproduce some of our results of Barrat-Charlaix et al. [12](Figure 3) and extend them to SARS-CoV-2. To quantify predictability, we ask the following question: given the state of a viral population at times 0, 1,..., *t*, what can we say about variant frequencies at times *t* + 1,...? We performed a retrospective analysis of viral evolution and identified all amino-acid mutations that were observed to grow from frequency 0 to an arbitrary threshold $x^*$. Adaptive beneficial mutations should in principle be overrepresented in this group and if they provide a persistent fitness advantage, we would expect them on average to keep on growing beyond $x^*$. **Figure 3** shows these trajectories for the amino acid substitutions in the HA protein of A/H3N2 influenza, using data from 2000 to 2023, and the SARS-CoV-2 genome using data from 2020 to 2023. Panels on the left show all trajectories that reached $x^* = 0.4$, with their average displayed in black. The panels on the right show the average trajectory for different threshold values $x^*$ between 0.1 and 0.8.

While the dynamics of the variants of the two viruses can not be compared directly due to vastly different sampling intensities and different rates of adaptation, the qualitative patterns differ strikingly. In the case of influenza, trajectories of seemingly adaptive mutations show little inertia and on average hover around $x^*$ instead of growing. This surprising result is in line with the study in [12](Figure 3) which used data from the period 2000–2018. On the other hand, trajectories of SARS-CoV-2 mutations show a much smoother behavior with steady growth beyond $x^*$. On longer timescales however, we observe a systematic decrease in frequency: this is explained by the particular initial dynamics of SARS-CoV-2, where new variants arose at a rapid pace and replaced old ones. This process is often called clonal interference and reduces long term predictability.

In our setting of eco-evolutionary adaptation, the random walk model predicts that the probability of fixation of an immune evasive variant is given by the final frequency $\beta$ of its initial partial sweep. Subsequent allele dynamics and diversity are governed by an anomalous coalescent process driven by the random walk defined in Eqs. 9 and 11, leading to little predictability of evolution. This abstraction should hold when partial sweeps are instant and do not overlap, meaning that the rate $\rho$ at which new variants emerge is small compared to their initial growth rate $s_0$.

To explore the behavior of our partial sweep model in a more general setting, we simulate evolutionary dynamics of a population under a Wright-Fisher model with expiring fitness dynamics. Simulations involve a population of *N* viruses with a genotype where each position can be in one of two possible states $\sigma_i \in \{0,1\}$. Fitness effects $S_i$ are associated with mutations at each position, and the total fitness of a virus is given by $F = \sum_i \sigma_i S_i$. At each generation, viruses with a fitness *F* expand by a factor $e^F$, and the next generation is constructed by sampling *N* individuals from the previous one. Following Eq. 11, mutational effects $S_i$ decrease by an amount $\nu x_i \cdot S_i$, where $x_i$ is the frequency at which mutation *i* is found in the population.

We simulate the emergence of adaptive variants in the following way. At a constant rate $\rho$, we pick one sequence position *i* that has no polymorphism and set the fitness effect of the corresponding mutation to an initial value $S_i$, with an amplitude drawn from probability distribution $P_s$ and the sign chosen such that the mutation is adaptive. In practice, we use an exponential distribution $P_s \propto e^{-s/s_0}$, meaning that the typical magnitude of initial fitness effects is described by only one
<figure><img src="image_8.png" /><figcaption>FIG. 3.</figcaption></figure>

Retrospective analysis of predictability of viral evolution: frequency trajectories of all amino acid substitutions that are observed to rise from frequency 0 to x* for **Top**: influenza virus A/H3N2 from 2000 to 2023, and **Bottom**: SARS-CoV-2 from 2020-2023. **Left**: all trajectories for x* = 0.4, with blue ones ultimately vanishing and red ones ultimately fixing. The average of all trajectories is shown as a thick black line. **Right**: showing only the average trajectories for different values of x* (grey lines).
parameter $s_0$. The corresponding distribution of partial sweep size is described in the SI. At the same time, we introduce the corresponding mutant in the population at a low frequency, picking its background genotype from a random existing strain. The behavior of the model is determined by (i) the distribution $P_β$ of partial sweep size depending on $v/s_0$, and (ii) the ratio of the variant emergence rate and their growth rate $ρ/s_0$, which determines how often sweeps overlap and interact. The probability of two sweeps overlapping is defined in Section (B 2) of the SI.

We use this simulation to address the question of predictability: given the state of the population at generations 0, 1,..., *t*, can we predict its state at future times *t* + 1, ...? Specifically, we ask whether we can predict the frequency $x(t + Δt)$ of a variant *A*, given it is at frequency $x$ at time *t*, as we did previously for influenza virus [12], see Fig. 3. The dynamics of isolated selective sweeps ($ρ/s_0 \ll 1$, $v/s_0 \ll 1$) should be perfectly predictable: after an initial stochastic phase when the variant is very rare, its frequency grows monotonically to fixation. This predictability decreases with increasing $ρ/s_0$ due to clonal interference [27, 28], for example when an adaptive variant is outcompeted by an even more adaptive one. We also expect predictability to decrease with increasing $v/s_0$ since sweeps are then partial and their ultimate fixation determined by subsequent variants with dynamics that resemble a random walk.

To quantify these effects, we select from a long simulation all rising frequency trajectories of adaptive mutations that cross an arbitrary threshold $x^*$. The results are shown in panel A of **Figure 4**, where we show the average $x(t)$ of rising frequency trajectories after crossing the threshold $x^* = 0.5$. We use three rates of fitness decay: $v \in [0, s_0/3, s_0, 3s_0]$ and low clonal interference $ρ/s_0 = 0.05$. The case $v = 0$ corresponds to a classical traveling wave scenario with constant fitness effects, and, as expected, is the most predictable: the average trajectory rises well above 0.5. For larger values of $v/s_0$, corresponding to a quicker decay of fitness, predictability gradually declines and becomes negligible for $v/s \ll 1$. Note that this matches quite well with the predictions from the random walk model where the average change in frequency $\langle \Delta x \rangle$ is null.

To explore parameter space more systematically, we quantify predictability as the probability of fixation $p_{fix}$ of rising variants that cross threshold $x^*$. In a perfectly predictable scenario with well separated selective sweeps, $p_{fix}$ should be close to 1 regardless of $x^*$, while it should be equal to $x^*$ in an unpredictable setting such as neutral evolution.

In panels **B**, **C** and **D** of **Figures 4** ⚙️, probability of fixations are shown for three values of $ρ/s_0$ and four values of $v/s_0$. Clonal interference increases when going from left to right among these panels (increasing $ρ/s_0$), while the intensity of fitness decay increases when going from blue to red curves (increasing $v$). Increasing either $ρ/s_0$ or $v/s_0$ reduces $p_{fix}$ towards the dashed diagonal corresponding to $p_{fix} = x^*$. However, as observed previously [12], in the classic scenario with stable fitness effects $v/s_0 = 0$ considerable predictability remains even in cases of strong interference (blue curve in panel D and **Figure S3** ⚙️). The strong interference setting is explored in more detail in the SI up to values $ρ/s_0 \approx 30$, using similar simulations but without expiring fitness $v = 0$. **Figure S3** ⚙️ shows that even in these cases of strong interference, $p_{fix}$ remains significantly above the diagonal.

Finally, we use our simulation to investigate typical levels of diversity in the population and the time to the most recent common ancestor. One quantity that can easily be estimated from the random walk model is the average pairwise coalescence time $T_2$, that is the typical time separating two random strains from their most recent common ancestor (MRCA). In the SI, we show that under the random walk approximation, $T_2 = 1/ρ(β^2)P_β$, which in neutral models of evolution would correspond to the effective population size $N_e$. A more detailed analysis of the coalescent process reveals that the random walk approximation corresponds to the so-called A-coalescent [29, 30] (see SI).
<figure><img src="image_9.png" /><figcaption>FIG. 4.</figcaption></figure>

Simulations under the Wright-Fisher model with expiring fitness. **A:** Average frequency dynamics of immune escape mutations that are found to cross the frequency threshold $x^* = 0.5$, for four different rates of fitness decay. If the growth advantage is lost rapidly (high $v/s_0$), the trajectories crossing $x^*$ have little inertia, while stable growth advantage (small $v/s_0$) leads to steadily increasing frequencies. **B,C,D:** Ultimate probability of $P_{fix}(x)$ of trajectories found crossing frequency threshold $x$. Each panel corresponds to a different rate of emergence of immune escape variants, with four rates of fitness decay per panel. Increased clonal interference $\rho/s_0$ and fitness decay $v/s_0$ both result in gradual loss of predictability. We use $s_0 = 0.03$. **E:** Time to most recent common ancestor $T_{MRCA}$ for the simulated population, as a function of the prediction obtained using the random walk $N_e = 1/\rho(\beta^2)$. Points correspond to different choices of parameters $\rho$ and $P_\beta$, and a darker color indicates a higher probability of overlap as computed in the SI.
In panel E of Figure 4 ⚙️, the average time to the common ancestor of pairs of strains in the population is plotted as a function of $T_2$ predicted by the random walk model. Each point in the figure corresponds to one simulation of long duration with a given distribution of partial sweep size $P_β$ and a given $ρ$ setting $T_2$, with darker color indicating a higher probability of overlap as computed in the SI. We find a good agreement between the empirical time to MRCA and the estimation from the random walk, at least as long the probability of overlap between successive partial sweeps is small (indicated by shading, see SI). With increasing overlaps, coalescence slows down and diversity increases: points in darker shades of red tend to have larger time to MRCA than what is expected from the distribution of $β$. This is expected intuitively: if another adaptive variant emerges before the previous one has reached its final frequency, it has a lower probability to land on the same background and thus tends to be in competition with the first variant. This leads to a smaller effective $β$ which slows the dynamics.

Discussion

Evolutionary adaptation is often pictured as an optimization problem in a static environment. In many cases, however, this environment is changed by the presence of the evolving species, for example because a host population develops immunity or a dynamic ecology. Here, we have explored consequences of such eco-evolutionary dynamics in a case of host-pathogen co-evolution where different variants of a pathogen shape each other's environment through generation of cross-immunity.

Influenza virus evolution has been the subject of intense research with efforts to predict the composition of future viral populations [9 ⚙️–11 ⚙️, 31 ⚙️]. The A/H3N2 subtype in particular undergoes rapid antigenic change through frequent substitutions in prominent epitopes on its surface proteins [32 ⚙️–35 ⚙️]. Given the clear signal of adaptive evolution, one might expect A/H3N2 to be predictable in the sense that variants that grow keep growing. Yet, it has been difficult to find convincing signals of fit, antigenically novel, variants that consistently grow and replace their competitors [11 ⚙️, 12 ⚙️]. In contrast, SARS-CoV-2 evolution has been consistently predictable in the sense that dynamics are well modeled by exponentially growing variants that compete for a common pool of susceptible hosts. However, even in this case, taking into account the immune adaptation of hosts leads to a better description of variant dynamics [4 ⚙️].

We have shown that depending on (i) the heterogeneity of immunity in the population, (ii) the asymmetry between backward and forward cross-immune recognition, and (Hi) waning or turnover of immunity, the immune escape can either lead to dynamics dominated by selective sweeps, or to one were escape mutations have an initial growth advantage that dissipates before the variant fixes. The former scenario is observed when initial growth is fast, backward immunity high, and waning slow compared to variant dynamics. In this case, new variants can rise to high frequency driven by their own advantage and fix. Immunological heterogeneity slows down the initial rise, allowing for population immunity to respond and adjust before the variant has fixed.

This process of partial sweeps reconciles two seemingly contradicting observations: HA evolution in human influenza A virus is clearly driven by adaptive immune escape and most substitutions are clustered in epitope regions [36 ⚙️]. On the other hand, most substitutions do not sweep to fixation but tend to meander in a quasi-neutral fashion [12 ⚙️]. In the partial sweep scenario proposed here, diversity is dominated by immune escape mutations that are rapidly brought to macroscopic frequency by their initial growth advantage, but their ultimate fate is determined mostly by subsequent mutations.

In any real-world scenario, there will be a variety of mutations, including some mutations that perform complete selective sweeps, either because they escape immunity of a large fraction of the population (M small), because they generate robust immunity against previous strains ("back-
boost" [37 🤷 ], or because the increase the intrinsic transmissibility of the virus (for example reverting a previous escape mutation that had a deleterious effect on transmissibility). The degree to which partial sweeps matter will vary from virus to virus and will change over time. Recently emerged viruses circulate in a homogeneous immune landscape and adapt to the new host for some time, consistent with rapid and complete sweeps of variants in SARS-CoV-2. Similarly, influenza virus A/HINlpdm, which emerged in humans in 2009, exhibited more consistent trajectory dynamics than A/H3N2 [12 🤷 ],

More generally, qualitative features of the partial sweep dynamics investigated here are expected
to exist in any system where the environment responds to evolutionary changes on time scales
comparable to the time it takes for the adaptive variants to take over, leading to eco-evolutionary
dynamics [38 🤷 ]. In ecological systems involving eukaryotes, it is the evolutionary part of this
interaction that is thought of as slow, while ecology is fast. In the case of rapidly adapting RNA
viruses in human populations with long-lived immunological memory, models often assume that
viral adaptation is fast while hosts have long-lasting memory. The most complex and least
predictable dynamics is expected when the evolutionary and ecological time scales are similar and
different host pathogen systems will fall on different points along this axis.

Code availability

* Figures in the main text can be reproduced using a set of notebooks at https://github.com/PierreBarrat/ExpiringFitnessFigures 🤷 .

* Code for the simulations of the SI model is available at https://github.com/PierreBarrat/PartialSweepSIR.j1 🤷

* Code for the population dynamics simulation is available at https://github.com/PierreBarrat/WrightFisher.j1 🤷

* Code to generate empirical frequency trajectories and their averages is available as scripts 'fluoprogram' and 'sc2oprogram' in https://github.com/nextstrain/flu_frequencies on branch 'fixation'.

Data availability

Sequence data of influenza viruses was obtained from GISAID [39 🤷 ]. We thank the teams involved in sample collection, sequencing, and processing of these data for their contribution to global surveillance of influenza virus circulation. A table acknowledging all originating and submitting laboratories is provided as supplementary information.

Sequence data of SARS-CoV-2 viruses was obtained from NCBI and restricted to data from North America to ensure more homogeneous sampling. We are grateful to all teams involved in the collection and generation of these data for generously sharing these data openly.

Acknowledgements

We gratefully acknowledge research support from the University of Basel (core funding) and the
Swiss National Science Foundation (grant 310030-188547).
# Declaration of Interests

RAN consults for Moderna TX on matters of virus evolution. Other authors declare no competing interests.

## A. SIR model

### 1. Equilibrium of the SIR model with one immune group

To help us compute the equilibrium reached by the SIR model, we introduce additional notation: the vectors $\vec{s} = [S^1, ..., S^N]$ and $\vec{r} = [I^1, ..., I^N]$ will respectively represent the compartments of hosts susceptible and infectious to each strain; the matrix K describes the cross-immunity; the vector $\vec{r}$ is a vector of dimension N whose elements are all equal to 1.

We derive the equilibrium state for Equations 1 & 2. For each strain a, equilibrium for $I^a$ is reached when $S^a = \delta/\alpha$. We thus have $\vec{s}_eq = \delta/\alpha \vec{r}$. Introducing this in the equation for the derivative of $S^a$, we obtain the following equilibrium:

$$ 
\begin{aligned}
\vec{S}_{eq} &= \frac{\delta}{\alpha} \vec{1}, \\
\vec{I}_{eq} &= \frac{\gamma}{\delta} \left( 1 - \frac{\delta}{\alpha} \right) \mathbf{K}^{-1} \vec{1}.
\end{aligned} 
\quad (A1) $$

An interesting remark is that at equilibrium, $I$ is of order $\gamma/\delta \ll 1$. Note that the structure of $\mathbf{K}$ makes it invertible in most cases. Indeed, we impose $K^{aa} = 1$ and $0 \le K^{ab} < 1$ for $a \neq b$.

### 2. Equilibrium for two viruses with one immune group

We consider the case where two viruses are present, called wild-type (wt) and mutant (m).

The cross-immunity is represented by a 2 × 2 matrix

$$ \mathbf{K} = \begin{bmatrix} 1 & b \\ f & 1 \end{bmatrix}, \quad (A2) $$

As shown in the previous section, the equilibrium is given by

$$ 
\begin{aligned}
S^{wt} &= S^m = \frac{\delta}{\alpha} \\
\vec{I} &= \frac{\gamma}{\delta} \left( 1 - \frac{\delta}{\alpha} \right) \mathbf{K} \vec{1},
\end{aligned} 
\quad (A3) $$
where ⌜ stands for [I<sup>wt</sup>, I<sup>m</sup>] and ⌜ for [1, 1]. It is straightforward to invert the crossimmunity matrix, and we obtain

$$
\begin{align}
I^{wt} &= \frac{\gamma(1 - \delta\alpha^{-1})}{\delta(1 - b\alpha)}(1 - b) \nonumber \\
I^{m} &= \frac{\gamma(1 - \delta\alpha^{-1})}{\delta(1 - b\alpha)}(1 - f). \tag{A4}
\end{align}
$$

Note that without cross-immunity, the number of infected by either virus would be <u>1</u>. Positive values of b and f thus have the effect of lowering the equilibrium values of I<sup>m</sup> and I<sup>wt</sup> with respect to the absence of cross-immunity.

It is interesting to compute the fraction of infections due to the mutant at equilibrium. This is
easily derived from the relations above:

$$
\frac{I^m}{I^{wt} + I^m} = \frac{1 - f}{2 - b - f}. \tag{A5}
$$

A few observations can be made:

* if b = 1 and f < 1, then the wild-type vanishes at equilibrium and the mutant reaches frequency 1. In this case, the presence of the mutant alone is enough to keep S<sup>wt</sup> to its threshold value R<sub>0</sub><sup>-1</sup>, making it impossible for the wild-type to grow.

• inversely, if f = 1, then the mutant stays at frequency 0.

* if b, f < 1, the mutant will reach a finite frequency x, with x > 0.5 if b > f and x < 0.5 if b < f.

**3. Equilibrium without the mutant**

We first derive the equilibrium situation before the mutant virus is introduced in the case with
only one immune group. We remind that in this case there is only one cross-immunity matrix
which has the form

$$
K = \begin{pmatrix} 1 & b \\ f & 1 \end{pmatrix},
$$

where *b* is the immunity to the wild-type caused by an infection with the mutant, and *f* the reverse.

Since the mutant is absent from the host population, we assume I^m = 0. The equilibrium values for S^wt and I^wt are easily obtained from the dynamical equations:

$$
\begin{align}
S^{wt} &= \delta/\alpha = R_0^{-1} \nonumber \\
I^{wt} &= \frac{\gamma}{\delta}(1 - S^{wt}). \tag{A6}
\end{align}
$$
We then set the derivative of $S^m$ to 0:

$$ 
\begin{aligned}
& -\alpha f S^m I^{wt} + \gamma(1 - S^m) = 0 \\
\rightarrow & S^m = \frac{1}{1 + f(R_0 - 1)} \qquad (A7) \\
& = \frac{\delta}{\delta + f(\alpha - \delta)} > \frac{\delta}{\alpha} = S^{wt}.
\end{aligned}
$$

Since we assume $f < 1$, the initial number of susceptibles to the mutant will be larger than $\delta/\alpha$, allowing initial growth of the mutant. Using the dynamical equation for $I^m$, the initial growth rate of the mutant can be written as

$$ \dot{I}^m(t = 0) = \alpha S^m - \delta = \delta \left( \frac{\alpha}{\delta + f(\alpha - \delta)} - 1 \right). \qquad (A8) $$

If $f=0$, the growth rate is $\alpha - \delta$, i.e. the one expected in a fully naive population.

If $f=1$ however, the growth rate is 0 as the wild-type confers perfect immunity to the mutant.

The equations above generalize to more immune groups. Cross-immunity matrices $K_i$ now depend on parameters $f_i$ and $b_i$, and the initial number of susceptibles in immune group $i$ is given by

$$ S_i^m = \frac{\delta}{\delta + f_i(\alpha - \delta)}. \qquad (A9) $$

In a given immune group $i$, the mutant growth rate is proportional to $S^m - \delta/\alpha$. The growth rate of the mutant will thus be initially faster in immune groups for which it is antigenically different, i.e. $f_i < 1$, than in groups where it is similar to the wild-type, i.e. $f_i \approx 1$.

In the case of a well-mixed population, that is $C_{ij} = 1/M$, we can write the growth of the infections by the mutant $I^m = \sum_{i=1}^{M} I^i$ as an exponential growth with a time dependent rate. In this case, the overall growth rate is given by the derivative of $I^m = \sum_{i=1}^{M} I^i$:

$$ \dot{I}^m = \left( \frac{\alpha}{M} \sum_{i=1}^{M} S_i^a - \delta \right) I^m \qquad (A10) $$

In particular, using the invasion scenario from the main text with $\varepsilon = 0$ (i.e. $f_i = 1$) in $M-1$ groups and an arbitrary value $f$ in group 1, we obtain the following growth rate at $t = 0$:

$$ \dot{I}^m = \frac{\delta}{M} \left( \frac{\alpha}{\delta + f(\alpha - \delta)} - 1 \right) I^m. \qquad (A11) $$

That is, the initial growth rate for $M$ groups is $M$ times smaller than the one in the single group case.

In the case of a non well-mixed population, i.e. arbitrary $C_{ij}$, it is not possible to write a pseudo-exponential growth rate as in Eq A10. However, it is clear that the initial growth rate will also be smaller as in the single group case since the mutant initially only grows in group $i = 1$.
## 4. Equilibrium with M immune groups

For M immune groups and arbitrary cross-immunity matrices $K_i$, the equilibrium frequency of the two strains is not easy to compute. However, it is possible to give an analytical expression in the regime of fast mixing $C_{ij} = M^{-1}$ and when the two strains differ immunologically for only one group, i.e. for matrices

$$ K_1 = \begin{pmatrix} 1 & b \\ f & 1 \end{pmatrix} \quad \text{and} \quad K_j = \begin{pmatrix} 1 & 1 \\ 1 & 1 \end{pmatrix} \quad \text{for } j > 1. $$

Note that this corresponds to the situation studied in the main text with fast mixing and $\epsilon \to 0$. We show here that in this case, the equilibrium frequency for all immune groups is the same as the one obtained for only one immune group with matrix $K_1$. In other words, the expression for $\beta$ in Eq. 6 is still valid.

To prove this, we assume the following form for the solution of the dynamical equations:

$$ \forall i \in \{1...M\}, \forall a, S_i^a = M \frac{\delta I_i^a}{\alpha I^a} \quad \text{and} \quad \frac{I_i^a}{I_i} = \frac{I^a}{I} = \nu^a, \qquad (A12) $$

where the index $a$ runs over all strains (here wild-type and mutant), and where we have defined the infectious levels for group $i$, for strain $a$ and globally:

$$ I_i = \sum_a I_i^a, \quad I^a = \sum_i I_i^a, \quad I = \sum_{i,a} I_i^a. $$

Note that the second equation in Eq. SA12 means that the frequency $\nu^q$ of strain $a$ is the same across all immune groups and consequently also globally.

We now show that injecting these expressions of $S$ and $I$ in the dynamical system and solving for $\nu^a$ gives the expected result. First, note that with this choice of $S^a$, the derivative of $\nu^a$ given by Eq. 1 immediately vanishes. We thus concentrate on $\dot{S}_i^a$ given by Eq. 2. For any immune group $i$ and strain $a$, we have

$$ \dot{S}_i^a = 0 = -\delta \frac{I_i^a}{I^a} \sum_b K_i^{ab} I^b + \gamma \left( 1 - M \frac{\delta I_i^a}{\alpha I^a} \right). $$

where we have used $C_{ij} = 1/M$ and $\nu^a = \sum_i \nu_i^a$ to remove the sum on immune groups. Multiplying this equation by $\nu^a = I^a/I$, we obtain

$$ -\delta I_i^a \sum_b K_i^{ab} \nu^b - \gamma \left( \nu^a - M \frac{\delta I_i^a}{\alpha I^a} \right) = 0. $$
We now eliminate $\varepsilon$ by using the expression $\varepsilon = I_i\nu^a$:

$$\delta I_i \nu^a \sum_b K_i^{ab} \nu^b - \gamma \left(1 - M \frac{\delta I_i}{\alpha I}\right) \nu^a = 0,$$

$$\delta I_i \sum_b K_i^{ab} \nu^b - \gamma \left(1 - M \frac{\delta I_i}{\alpha I}\right) = 0.$$

Note that this last expression is true for all strains a. Considering any two strains a and b, we can thus write

$$\delta I_i \sum_c K_i^{ac} \nu^c - \gamma \left(1 - M \frac{\delta I_i}{\alpha I}\right) = 0,$$

$$\delta I_i \sum_c K_i^{bc} \nu^c - \gamma \left(1 - M \frac{\delta I_i}{\alpha I}\right) = 0,$$

$$\Rightarrow \forall a, b \sum_c (K_i^{ac} - K_i^{bc}) \nu^c = 0,$$

where the last expression is obtained by subtracting the two previous ones. First, we see that for $i > 1$ any frequency vector $\nu$ is a solution since $\kappa_i^a = 1$ for all a, b. For $i = 1$ and defining $\nu^m = \beta$ and $\nu^{wt} = 1 - \beta$ and using the expression for $K_1$, we obtain

$$(1 - f)(1 - \beta) + (b - 1)\beta = 0 \\ \Rightarrow \beta = \frac{1 - f}{(1 - b) + (1 - f)}$$

as claimed.

This result is not completely trivial and should be commented. In this setting, the mutant escapes immunity built by the wild-type for a fraction $1/M$ of the population, and yet it reaches the same frequency as in the case with one immune group. This can be rationalized as follows: for immune groups $i > 1$, the cross-immunity matrix is such that the wild-type and mutant strains are completely equivalent. If immune group 1 was not here, the mutant could thus equilibrate at any frequency between 0 and 1. Since it is initially introduced at very low frequency, it would remain marginal in immune groups $i > 1$. However, since its “natural” equilibrium frequency in group $i = 1$ is $\beta$ and since the groups are connected, equilibrium is reached when the mutant reaches frequency $\beta$ in all groups.

Note that if we take the situation of the main text with

$$K_i = \begin{pmatrix} 1 & \varepsilon \\ \varepsilon & 1 \end{pmatrix}$$

and $\varepsilon > 0$, the expressions above do not hold. However, if $\varepsilon \ll 1$, the perturbation is small and we expect an equilibrium frequency close to $\beta$, which the case in **Figure 1**.
## 5. Realistic modeling of host's immune state

The SIR model proposed in the main text relies on the assumption of immunity acquisition through exposure. This explains terms like $-\alpha \sum_b S^a K^{ab} I^b$ in the derivative of $S^a$: acquiring immunity to a through cross-immunity $K^{ab}$ requires a combination of prior susceptibility and exposure to strain b. Importantly, this does not depend on the immune state of the hosts with respect to strain b.

A more realistic representation would be one where acquiring immunity to strain a from exposure to b requires being infected by b. However, this would require keeping track more precisely of the immune state of hosts, as we would need to separate hosts into two groups, namely

* hosts who are susceptible to both a and b, and can thus acquire immunity to a through infection by b;

* hosts who are susceptible to a but immune to b, and can no longer acquire immunity to a through infection by b.

To test whether our results are robust to such changes of hypothesis, we write a simple SIR model with two strains a and b where cross-immunity is only activated through infection rather than exposure. To properly track the immune states of the hosts, we introduce the groups $R^a$ and $R^b$ respectively representing hosts immune to only a or only b, and $R^{ab}$ representing hosts immune to both a and b. The compartment $R^0 = 1 - R^a - R^b - R^{ab}$ groups hosts susceptible to both strains. It is simpler in this case to write the dynamics in terms of compartments I and R, rather than I and S as in the main text. For simplicity, we do not use immune groups here. The dynamics involve two equations for the infected:

$$ 
\begin{aligned}
I^a &= \alpha(1 - R^a - R^{ab})I^a - \delta I^a, \\
I^b &= \alpha(1 - R^b - R^{ab})I^b - \delta I^b,
\end{aligned} 
\qquad (A13) 
$$

and three for the immune:

$$ 
\begin{aligned}
R^a &= \alpha(1 - K^{ba})R^0 I^a - \gamma R^a, \\
R^b &= \alpha(1 - K^{ab})R^0 I^b - \gamma R^b, \\
R^{ab} &= \alpha R^0 (K^{ba} I^a + K^{ab} I^b) + \alpha(R^a I^b + R^b I^a) - \gamma R^{ab}.
\end{aligned} 
\qquad (A14) 
$$

In Figure S1, we show both that the dynamical and equilibrium properties of this model are qualitatively the same as the one from the main text. On the left panel, we show that the dynamics of this new model does not differ qualitatively from the model of the main text. In particular, in the invasion scenario, the frequency of the variant converges to some equilibrium value after some oscillations. On the right, we show that this equilibrium value $\beta$ is different but relatively close to the one from the main text.

## 6. Change in frequency when adding subsequent strains

This section shows that under certain condition, adding a new variant to the SIR model does not change the relative frequencies of previous variants. This is an important condition for the random walk of the main text to be valid.
Here is a quick summary of the results proved below. Adding a variant to the SI model involves adding a column $\vec{b}$ and a row $\vec{f}$ to the cross-immunity matrix, which can be given by two vectors. If these vectors only depend on one parameter, i.e. $\vec{b} = b \cdot \vec{f}$ and $\vec{f} = f \cdot \vec{b}$, then the relative frequency of previous strains is unchanged in the new equilibrium. What this means in practice is that the new strain must be at equal “antigenic distance” from all previous strains. A possible interpretation is an antigenic space of infinite dimensions: all mutations explore an antigenic region which is new.

We start from an initial situation where there are $N$ variants with an $N \times N$ cross-immunity matrix $\mathbf{K}$. At equilibrium, the number of hosts infected by each virus $a \in \{1... N\}$ is given by the elements of the vector $I$ that can be computed from the cross-immunity matrix and parameters of the model:

$$I = \frac{\gamma(1 - \delta/\alpha)}{\gamma + \delta} \mathbf{K}^{-1} \vec{1}_N, \qquad (A15)$$

where $\vec{1}_N$ is a vector containing only 1's and of length $N$. The relative frequency of variant $a$ with respect to variant $b$ is simply defined as

$$f_{ab} = \frac{I_a}{I_a + I_b}. \qquad (A16)$$

We assume that the initial population has reached equilibrium.

We now add a new virus to this population, with index $N + 1$. The new cross-immunity matrix $\tilde{\mathbf{K}}$ is now written as

$$\tilde{\mathbf{K}} = \begin{bmatrix} \mathbf{K} & \vec{b} \\ \vec{f}^\top & 1 \end{bmatrix}, \qquad (A17)$$

where $\vec{b}$ and $\vec{f}$ are two vectors of length $N$. This is a general way to write that the backward cross-immunity to variant $a$ caused by an infection with the new variant $N + 1$ is $b_a$. Inversely, the forward cross-immunity to variant $N + 1$ caused by an infection with an old variant $a$ is $f_a$.

This new cross-immunity matrix will of course result in a new equilibrium for the number of infected hosts, given by the vector $\vec{I}$:

$$\tilde{I} = \frac{\gamma(1 - \delta/\alpha)}{\gamma + \delta} \tilde{\mathbf{K}}^{-1} \vec{1}_{N+1}.$$

The question we ask here is whether the relative frequency of two variants $1 \le a,b \le N$ is changed by the addition of the new variant. In other words, we want to know whether the equality below holds:

$$\frac{I_a}{I_a + I_b} \stackrel{?}{=} \frac{\tilde{I}_a}{\tilde{I}_a + \tilde{I}_b}.$$
<figure><img src="image_10.png" /><figcaption>Figure S 1.</figcaption></figure>

**Left:** Dynamics of the frequency of the variant for the SIR model from **Equations A13** & **A14** using the invasion scenario from the main text. Two 2 × 2 cross-immunity matrices are used, with off-diagonal parameters *f* and *b* chosen to give the same equilibrium. The gray line represents the equilibrium that would be obtained using the model of the main text. **Right:** Equilibrium frequency *β* for this new SIR model (y-axis) versus the *β* from the main text (x-axis). Each point corresponds to a given pair (*f*, *b*).
Below, we prove this equality under a condition for cross-immunity of the new variant $\vec{b}$ and $\vec{f}$:

$$ \vec{b} = b \cdot \vec{1}_N \quad \text{and} \quad \vec{f} = f \cdot \vec{1}_N, \qquad (A18) $$

where $0 < b, f < 1$ are scalars. This amounts to say that cross-immunity is the same between the new variant $N+1$ and any old variant $a$, i.e. that the new variant is at equal antigenic distance from all previous variants.

To prove the equality, we perform the computation $\vec{K}^{-1} \cdot \vec{f} \cdot \vec{1}_N$. To do that, we make use of the following formula for inverting a block matrix:

$$ \begin{bmatrix} A & B \\ C & D \end{bmatrix}^{-1} = \begin{bmatrix} A^{-1} + A^{-1}B\Lambda C A^{-1} & -A^{-1}B\Lambda \\ -\Lambda CA^{-1} & \Lambda \end{bmatrix}, $$

where we defined $\Lambda = (D - CA^{-1}B)^{-1}$. The following identities map to our problem:

$$ A = \mathbf{K}, \quad B = \vec{b}, \quad C = \vec{F}^{\top}, \quad D = 1. $$

We immediately see that $\Lambda = (1 - \vec{f} \cdot \vec{K}^{-1} \vec{b})^{-1}$ reduces to a scalar that we note $A$ for more clarity. We also define the other scalar value $\mu = \vec{f} \cdot \vec{K}^{-1} \vec{1}_N$. A few manipulations give us the following for $\vec{K}^{-1}$:

$$ \tilde{\mathbf{K}}^{-1} = \begin{bmatrix} \mathbf{K}^{-1} + \lambda \mathbf{K}^{-1} \vec{b} \cdot \vec{f}^{\top} \mathbf{K}^{-1} & -\lambda \mathbf{K}^{-1} \vec{b} \\ -\lambda \vec{f}^{\top} \mathbf{K}^{-1} & \lambda \end{bmatrix}. $$

Multiplying this by $\vec{1}_{N+1}$ results in

$$ \begin{align*} \tilde{I} &= \tilde{\mathbf{K}}^{-1} \vec{1}_{N+1} \\ &= \underbrace{\left( \mathbf{K}^{-1} + \lambda \mathbf{K}^{-1} \vec{b} \cdot \vec{f}^{\top} \mathbf{K}^{-1} \right) \vec{1}_{N} - \lambda \mathbf{K}^{-1} \vec{b} ; \underbrace{-\lambda \vec{f}^{\top} \mathbf{K}^{-1} \vec{1}_{N} + \lambda}_{\text{scalar}}}_{\text{size } N} \\ &= [I + b f \mu \lambda I - b \lambda I ; \lambda(1 - \mu f)] \\ &= [(1 - b \lambda + b f \lambda \mu) I ; \lambda(1 - \mu f)], \end{align*} $$

where we used the equalities $\mathbf{K}^{-1} \vec{b} = b \mathbf{K}^{-1} \vec{1}_N = bI$.

This result essentially shows that after adding the new variant, the fraction of hosts infected by the previous variants if simply multiplied by a scalar value $1 - by(1 - f\mu)$. This implies that the relative frequencies of the original variants are conserved when adding a new one.
## 7. Case with intrinsic fitness effects

In the SI model of the main text, we assume that the transmission rate $a$ is the same for the different strains. It is also interesting to investigate the case where this transmission rate varies. Here, we study a simple extension of the SI model without immune groups where there are two variants - mutant and wild-type - with respective transmission rates $\alpha^{wt} = \alpha\phi^{wt}$ and $\alpha^m = \alpha\phi^m$. The quantities $\phi^{wt}$, $\phi^m \in [\delta/\alpha, \infty)$ [can be interpreted as intrinsic fitness values for the two strains. Note that if $\phi^a < \delta/\alpha$, the strain $a$ cannot grow even in a fully susceptible population. The cross-immunity is as usual defined by matrix $K$ with off-diagonal terms $f$ and $b$.

The equations of motion are now

$$ 
\begin{aligned}
\dot{S}^a &= -\alpha S^a \sum_{b \in \{wt,m\}} \phi^b K^{ab} I^b + \gamma(1 - S^a) \\
\dot{I}^a &= \alpha \phi^a S^a I^a - \delta I^a.
\end{aligned}
\qquad (A19) $$

Computing the equilibrium, we immediately obtain

$$ 
\begin{aligned}
S^a &= \frac{\delta}{\alpha\phi^a} \\
I &= \frac{\gamma}{\delta} \cdot \mathbf{G}^{-1}\vec{h},
\end{aligned}
\qquad (A20) $$

where we have defined the following quantities:

$$ h^a = \frac{\alpha\phi^a - \delta}{\alpha\phi^a}, \quad \mathbf{G} = \begin{pmatrix} 1 & bs \\ fs^{-1} & 1 \end{pmatrix}, \quad s = \frac{\phi^m}{\phi^{wt}}. \qquad (A21) $$

The quantities $h^a$ can be interpreted as a scaled growth rate of each variant given a fully susceptible population, and the matrix $G$ combines the cross-immunity and the ratio of fitness values $\phi$. Note that it is straightforward to generalize these equations to an arbitrary number of strains: the relevant quantity will be the scaled cross-immunity matrix defined by $G^a = \frac{\phi^a}{\phi^{wt}}$.

Inverting $G$ and simplifying the equations a bit, we obtain

$$ 
\begin{aligned}
I^{wt} &= \frac{\gamma}{\delta} \frac{\alpha\phi^{wt} - \delta}{\alpha\phi^{wt}} \cdot \frac{1 - b\xi}{1 - bf}, \\
I^{m} &= \frac{\gamma}{\delta} \frac{\alpha\phi^{m} - \delta}{\alpha\phi^{m}} \cdot \frac{1 - f\xi^{-1}}{1 - bf},
\end{aligned}
\qquad (A22) $$

where we defined

$$ \xi = \frac{\alpha\phi^m - \delta}{\alpha\phi^{wt} - \delta}. \qquad (A23) $$

Note the interesting structure of equations A22: for each variant, they involve a first term $\frac{\phi^a}{\phi^{wt}}$ that depends only of the intrinsic growth rate of the mutant, and a second $\frac{\phi^m}{\phi^{wt}}$ that involves cross-immunity and relative growth rate through $\xi$.
These equilibrium equations give us two conditions for the co-existence of the two variants,

$$b < \xi^{-1} \text{ and } f < \xi, \qquad (A24)$$

respectively corresponding to $I^{wt} > 0$ and $I^m > 0$. We mention three interesting cases below.

* If the mutant has an intrinsic fitness disadvantage $\xi < 1$, it will only be able to invade if $f < \xi$. Since $f$ represents the probability that a host becomes immune to the mutant if infected by the wild-type, this means that the immune “niche” of the mutant must be large enough when compared to $\xi$.

* Invertly, if the mutant is fitter and $\xi > 1$, the mutant is always able to invade. The wild-type only survives if $b < \xi^{-1}$, meaning that the immunity to the wild-type caused by the mutant must be small enough.

* If one considers a situation without total cross-immunity, *i.e.* $b = f = 1$, the only way a mutant invades is if $\xi > 1$ meaning $\phi^m > \phi^{wt}$, and the result is a full selective sweep.

## 8. Oscillations of the SI model

The SI model from the main text tends to oscillate while returning to equilibrium. Here, we study this behavior in the simple case of one immune group ($M = 1$) and two viruses (wild-type and variant).

The idea is to linearize the dynamical equations around the equilibrium. This gives us

$$\dot{X} = QX,$$

where $X = [S^{wt}, S^m, I^{wt}, I^m]$ and

$$Q = \begin{pmatrix} -\alpha\gamma/\delta & 0 & -\delta & -\delta b \\ 0 & -\alpha\gamma/\delta & -\delta f & -\delta \\ g_1 & 0 & 0 & 0 \\ 0 & g_2 & 0 & 0 \end{pmatrix}$$

$$g_1 = \frac{\gamma}{\delta}(\alpha - \delta)\frac{1 - b}{1 - b f}$$

$$g_2 = \frac{\gamma}{\delta}(\alpha - \delta)\frac{1 - f}{1 - b f}.$$

To quantify the convergence to equilibrium is the frequency of the oscillations, we need the eigenvalues of matrix $Q$. For low enough $\gamma$, we can prove that the four eigenvalues are

$$\lambda_1 = -\frac{1}{2}\frac{\alpha}{\delta}\gamma \pm i\left(\gamma(\alpha - \delta) - \frac{1}{4}\frac{\alpha^2\gamma^2}{\delta^2}\right)^{1/2}$$

$$\lambda_2 = -\frac{1}{2}\frac{\alpha}{\delta}\gamma \pm i\left(\gamma(\alpha - \delta)\frac{(1-b)(1-f)}{1-bf} - \frac{1}{4}\frac{\alpha^2\gamma^2}{\delta^2}\right)^{1/2}.$$

This is only valid if the terms in the square roots above are positive, which requires yto be small enough. In our setting, we assume $\gamma \ll \alpha\delta$, so this will always hold.
From the eigenvalues we can compute:

* the rate of convergence to equilibrium: $\alpha\gamma/2\delta$. This means that convergence is slower for a smaller waning rate $\gamma$.

* the two oscillation frequencies that appear:

$$\omega_{high} = \frac{\sqrt{\gamma(\alpha - \delta)}}{2\pi},$$

$$\omega_{low} = \omega_{high} \sqrt{\frac{(1 - b)(1 - f)}{1 - b f}}.$$

Note that since $(1 - b)(1 - f)/(1 - b f) \le 1$, we always have $\omega_{low} \le \omega_{high}$. With the values of the main text $\alpha = 3$, $\delta = 1$, $\gamma = 5 \cdot 10^{-3}$ (units are inverse of generations), we obtain $\omega_{high} \approx 0.016$. If one generation is one week, this gives us a period $\omega_{high}^{-1} \approx 60$ weeks, that is approximately one year.

## 9. Link between parameters of the SIR and expiring fitness models

The effective expiring fitness model used in the second part of this work is characterized by the system of differential equations

$$\dot{x} = sx(1 - x) \quad \text{and} \quad \dot{s} = -\nu xs,$$

where $x$ is the frequency of the mutant strain. Here, we try to express the dynamics of the SIR model in this form to find a link between its parameters and the quantities $s$ and $v$.

We first focus on the case with two strains and one immune group. The frequency of the mutant is $x = I^m/(I^m + I^{wt})$. Using the logit function $\psi(x) = \log(x/(1 - x))$ and the dynamical equations of the SIR, we find

$$\frac{d}{dt}\psi(x) = \alpha(S^m - S^{wt}), \qquad (A25)$$

which allows us to define the fitness in the SIR case: $s = \alpha(S^m - S^{wt})$. At the beginning of the invasion, the initial growth rate is readily computed:

$$s(t = 0) = \frac{(1 - f)(\alpha - \delta)}{\delta + f(\alpha - \delta)} \delta, \qquad (A26)$$

which is the same as the initial growth rate of $I^m$. Note that if $f = 1$, the initial growth rate is 0.

We then compute the time derivative of $s$ early in the invasion, when $I^{wt}$, $S^{wt}$ and $S^m$ are close to their equilibrium values. In this case, a straightforward calculation gives

$$\begin{aligned} \dot{s} &= -\alpha(I^{wt} + I^m) \cdot \alpha(S^m - bS^{wt}) \cdot x \\ &\approx -\alpha(I^{wt} + I^m)xs, \end{aligned} \qquad (A27)$$

where the approximation is valid if $1 - b \ll 1$. This would give an expiry rate of fitness $v = \alpha(I^{wt} + I^m)$ in the case of the SIR model.
These results can also be obtained in the case of immune groups. We then have the following expressions for s and v at t ≈ 0:

$$
\begin{align}
s &= \frac{\alpha}{M} \sum_{i=1}^{M} (S_{i}^{m} - S_{i}^{wt}) \nonumber \\
\dot{s} &= -\alpha(I^{wt} + I^{m})x \sum_{i=1}^{M} (S_{i}^{m} - b_{i}S_{i}^{wt}) \tag{A28} \\
&= -\alpha(I^{wt} + I^{m})xs, \nonumber
\end{align}
$$

where the second expression is again valid if $1 - b_i \ll 1$.

Another question is that of the link between cross-immunity parameters *b* and *f* and the distribution of partial sweep sizes *β*. The relation between cross-immunity and *β* given by **Eq. 6** shows that the distribution of partial sweep size depends on the distribution of both *b* and *f*. As we do not have a prior on how *b* and *f* should be distributed, we explore the case where 1 - *f* and 1 - *b* are exponentially distributed. In other words, we define *ε*<sub>*f*</sub> = 1 - *f* and *ε*<sub>*b*</sub> = 1 - *b*, with the following distributions

$$
P(\epsilon_f) \propto e^{-\epsilon_f / \lambda} \quad \text{and} \quad P(\epsilon_b) \propto e^{-\epsilon_b / \mu}.
$$

The expression of the partial sweep size becomes *beta* = *epsilon*<sub>*f*</sub>(*epsilon*<sub>*b*</sub> + *epsilon*<sub>*f*</sub>). Note that both *epsilon* should remain smaller than 1, which is not guaranteed with exponential distributions. However, this is not problematic if *mu* and *gamma* take small enough values.

These assumptions allow us to compute the distribution of *beta*:

$$
P(\beta) = \frac{\mu/\lambda}{\left(\frac{\mu}{\lambda}\beta + (1 - \beta)\right)^2}
$$

with support on the interval $0 \le \beta \le 1$. **Figure S2** shows the various shapes that $P(\beta)$ then takes for different values of the $\alpha/\gamma$ parameter. Note that if $\mu a > \gamma$, *b* tends to be higher than *f* and *β* is biased towards one. If $\mu = \gamma$, $P(\beta)$ becomes uniform on the [0,1] range.

The exponential distribution away from one of *f* and *b* is a reasonable assumption, and allows analytical derivation of *P(β)*. Of course, any other distribution of *f* and *b* could be considered. For this reason, we choose to use a Beta distribution for *P*<sub>β</sub> in the analysis of the main text, as it can accommodate various shapes.

B. Expiring fitness model and random walk

1. Clonal interference

In non recombining genomes with a large mutation rate, the appearance of many adaptive mutations in close succession leads to clonal interference. In this regime, beneficial mutations present on different background compete for fixation, and the success of a mutation does not depend only on its fitness effect but also on the global state of the population. For this reason, clonal interference causes a decrease in predictability: dynamics are not deterministic but rather
<figure><img src="image_11.png" /><figcaption>Figure S 2.</figcaption></figure>

Distribution of partial sweep size *β* if 1 – *f* and 1 – *b* are exponentially distributed. Left: Probability distribution function *P*(*β*) for various values of *μ*/*y*. Right: Mean and standard deviation of *β* as a function of *μ*/*y*.
depend on the precise structure of which mutation appeared on which background. For instance, a beneficial mutation that increases in frequency can be outcompeted by a fitter one before it has the time to fix, making the extrapolation of frequency trajectories difficult.

We conduct simulations to quantify how much predictability decreases because of clonal interference, based on the ones of the main text. We study a large population of $N = 10^5$ genomes of length $L$, where each genome position $i$ can be in either of two states $\sigma_i \in \{0,1\}$. Fitness effects $s_i \in \mathbb{R}$ are associated to each position, and the fitness of an individual is $F = \sum_i s_i \sigma_i$. To simulate the adaptation of the population, we proceed in the following way: at a constant rate $p$, we pick a position $i$ that is non polymorphic and set the fitness effect $S_i$ by sampling its magnitude from distribution $P_s$ and choosing its sign in a way that favors mutations (positive if $\sigma_i = 0$ is more frequent, negative otherwise). At the same time, the corresponding mutation ($\sigma_i = 0$ if $S_i < 0$ and inversely) is introduced in the population at a small frequency $\delta f = 0.01$. We choose $P_s$ to be an exponential distribution with scale $s_0 = 0.02$, in agreement with findings on the distribution of fitness effects in [27, 40].

There is no fitness decay in this simulation, and the parameters $N$, $s_0$ and $\delta f$ have numerical values such that neutral drift is mostly irrelevant to the dynamics of mutations. For example, without accounting for interference, the probability of fixation of a mutation of effect $s_0$ when it is introduced at frequency $\delta f$ is $P_{fix}(\theta f) \approx 1 - e^{-40}$. As a consequence, the two parameters governing the dynamics are the scale of fitness effects $s_0$ and the rate of introduction of beneficial mutations $\rho$. The ratio of these two quantities represents the amount of clonal interference: at low $\rho/s_0$, mutations are mostly independent, while at high $\rho/s_0$ they strongly interfere.

We measure the probability of fixation $p_{fix}(x)$ of mutations found in a frequency bin $[x - \delta x, x + \delta x]$ over a long simulation. We only consider mutations with increasing frequency, meaning that their frequency was below $x$ at all times and at some point was measured in the frequency bin. **Figure S3** shows $P_{fix}(x)$ as a function of $x$ for different values of $\rho/s_0$. According to intuition, a low clonal interference value leads to easily predictable fixations: whatever the frequency $x$ at which it is observed, a mutation that is increasing in frequency fixes with a very high probability. Increasing clonal interference clearly makes dynamics less predictable and closer to neutrality, with $p_{fix}$ approaching the diagonal line. However, even in a regime of strong interference, e.g. $\rho/s_0 = 32$, deviations from neutrality remain very clear.

## 2. Expiring fitness effects: sweep size and probability of overlap

This section gives a few results about the expiring fitness equations from the main text. We rewrite the equations here for reference:

$$ \dot{x} = sx(1-x), \quad \dot{s} = -\nu xs. $$

First, we prove the expression for the amplitude of the partial sweeps. We divide the equation for $\dot{x}$ by the one for $\dot{s}$ to obtain

$$ \frac{dx}{ds} = \nu^{-1}(x-1). $$

This immediately gives us

$$ x(s) = 1 + \lambda e^{s/\nu} $$
<figure><img src="image_12.png" /><figcaption>Figure S 3.</figcaption></figure>

Probability of fixation of mutations *P*<sub>fix</sub>(*x*) of mutation frequency trajectories found crossing the frequency threshold *x*. Fitness effects are exponentially distributed with fixed scale *s*<sub>0</sub> = 0.03. The blue to red gradient in colors corresponds to the increasing rate *ρ* at which mutations are introduced. Strong clonal interference regime is obtain when *ρ*/*s*<sub>0</sub> > 1, in which case good mutations are introduced in close succession and compete for fixation. At low *ρ*/*s*<sub>0</sub>, trajectories are very predictable and an increasing trajectory almost certainly fixes. Even for the highest *ρ*/*s*<sub>0</sub>, *P*<sub>fix</sub>(*x*) remains significantly larger than *x* and dynamics are visibly not neutral.
with $\gamma$ constant. At $t = 0$, we have $s = s_0$ and $x = x_0 \ll 1$, while for $t \to \infty$ we have $s = 0$ and $x = \beta$ to be determined. From the $t = 0$ case we obtain $\gamma = (x_0 - 1)e^{s_0/\nu}$, and from $t \to \infty$ we get $\beta = 1 + \gamma$. Assuming $x_0 \ll 1$, we obtain the result of the main text:

$$ \beta = 1 - e^{-s_0/\nu}. $$

We now try to find an expression for the probability that two partial sweeps overlap. First, we try to estimate the time it takes for one partial sweep to complete. While we could not solve the differential equations of the main text analytically, we can give an approximate expression for the time dependent frequency $x$ during the partial sweep:

$$ x(t) \approx \frac{x_0\beta}{x_0 + (\beta - x_0)e^{-st}}, $$

where $\beta$ is a function of $s$ and $\nu$ and $x_0 = x(t = 0)$. This is simply the expression of a logistic growth starting at $x_0$, and saturating at $\beta$. From there, we compute the time $T_r(s)$ it takes a partial sweep of initial fitness $s$ to reach a frequency $r\beta$ with $x_0\beta^{-1} < r < 1$. We quickly find

$$ T_r(s) = s^{-1} \log \left( \frac{\beta}{x_0} \frac{r}{1-r} \right). $$

We now consider that two consecutive partial sweeps of initial fitnesses $s_1$ and $s_2$ overlap if the first one is not yet at frequency $r\beta_1$ while the second one is already at $(1 - r)\beta_2$. In the figure of the main text, we use $r = 3/4$: an overlap occurs if the first sweep is not yet at 3/4 of its final value while the second one is already at 1/4 of its final value. Thus, for an overlap to occur, we need the time $\tau$ between the two partial sweeps to be smaller than $T_r(s_1) - T_{1-r}(s_2)$. For sweeps happening at rate $\sigma$, this has probability $1 - \exp(-\sigma(T_r(s_1) - T_{1-r}(s_2)))$. Since the two sweeps have random initial fitness effects, we find that the overall probability for two consecutive sweeps to overlap is

$$ P_r(\text{overlap}) = \int_0^\infty ds_1 ds_2 P_s(s_1) P_s(s_2) \{1 - e^{-\rho(T_r(s_1) - T_{1-r}(s_2))}\}. $$

This integrates over all possible pairs of sweep amplitudes (or initial fitnesses) and weighs them by the probability that the time between the two leads to an overlap. It is this quantity (computed numerically) that is used for the scale of the colorbar in panel E of Figure 4 ☐ of the main text.

### 3. Distribution of partial sweep size $\beta$

This section discusses the distribution of the size of partial sweeps $\beta$ in the context of Equation 11 ☐ of the main text as well as the choice of parameters for panel E of Figure 4 ☐.

A first interesting case is when fitness effects are exponentially distributed, with parameter $s_0$:

$$ P(s) = s_0^{-1} e^{-s/s_0}. $$
This is the distribution we use in most of the population simulations. We compute the corresponding distribution of $\beta$ in a straightforward way:

$$ 
\begin{aligned}
P(\beta < x) &= P(s < -\nu \log(1-x)) \\
&= s_0^{-1} \int_0^{-\nu \log(1-x)} e^{-s/s_0} ds \\
&= 1 - (1-x)^{\nu/s_0}.
\end{aligned}
$$

Taking the derivative with respect to $x$, we obtain

$$ P(\beta) \propto (1 - \beta)^{\nu/s_0 - 1}. $$

This distribution can accommodate various shape: for $v/s_0 > 1$ it peaks at 0, and for $v/s_0 < 1$ it peaks at 1. We can also compute the following formula for the moments of $\beta$:

$$ \langle \beta \rangle = \frac{s_0}{s_0 + \nu} \quad \text{and} \quad \langle \beta^2 \rangle = \frac{s_0}{s_0 + \nu} \frac{2s_0}{2s_0 + \nu}. $$

When investigating the coalescence time in the main text, we use a different distribution for fitness effects. In this case, we want a finer control over the second moment of $\beta$, and we decide to sample $\beta$ directly using a Beta distribution. The Beta distribution has a support over [0,1] and can accommodate many different shapes. It is defined by two parameters $a$ and $b$:

$$ P(\beta) \propto \beta^{a-1} (1 - \beta)^{b-1}. $$

In our case, it is more practical to parametrize it by its mean $m$ and variance $v$. For given $m$ and $v < m(1 - m)$ we have

$$ a = \left( \frac{m(1-m)}{v} - 1 \right) m, \quad b = \left( \frac{m(1-m)}{v} - 1 \right) (1-m) $$

In the case of panel E of Figure 4 and in order to explore a wide range of distributions, we used three values of $m$: {0.3, 0.6, 0.9}, and for each $m$

* a low variance $v = \varepsilon \cdot m^2$ with $\varepsilon = 10^{-5}$

* a high variance $v = m(1 - m)/3$

For a given set of parameters defining a Beta distribution, we decide on the fitness effects by first sampling a $\beta$ for each new adaptive mutation, and then computing $s$ by using Equation 12 from the main text. For each distribution $P_{\beta}$, the simulation is performed for 6 values of $\rho \in [0.003, 0.01, 0.018, 0.032, 0.056, 0.1]$.

## 4. Random walk: monotonous trajectories

Here we compute the probability that in the random walk defined in the main text, a trajectory starting at $x_0$ converges straight to 0 without ever taking a step up. While going to 0 requires an infinite amount of downward steps, the probability is still finite since the steps are increasingly likely to go down. For simplicity, we compute this for a fixed $\beta$.
If the random walk always goes down, its position at time step $t$ will be $x_t = (1 - \beta)^t x_0$. Since the probability of going down is $1 - x_t$, the probability of always going down is

$$P_{down} = \prod_{t=0}^{\infty} (1 - (1 - \beta)^t x_0). \qquad (B1)$$

We simplify this expression by taking the logarithm and assuming that $(1 - \beta)^t \ll 1$ for $t \ge 1$:

$$P_{down} \approx (1 - x_0)e^{x_0(1-\beta^{-1})}.$$

Since the random walk is invariant by the change $x - 1 - x$, we can easily compute the probability of a trajectory always going up, and thus of a monotonous trajectory going straight to either boundary 0 or 1.

The quality of the approximation is quite good, as can be seen from **Figure S4**. The same figure also shows that this probability is relatively high, even for $x_0$ far from either boundary.

## 5. Coalescent

Consider a partial sweep happening between generations $t$ and $t + 1$, with probability $\rho$. One individual $A$ in generation $t$ will then have $\beta N$ children in generation $t + 1$. Any individual in generation $t + 1$ has a probability $\beta$ of having $A$ as a direct ancestor, and a probability $1 - \beta$ of the opposite. If we consider $n$ lineages at generation $t + 1$ and look backward in time, the probability that at least $k$ out of $n$ have $A$ as an ancestor is $\beta^k$. Averaging over $P^\beta$, we find the probability of $k$ specificlineages to have a common ancestor in the previous generation:

$$q_k = \rho \langle \beta^k \rangle.$$

Another useful quantity is the probability $\gamma_n(k)$ that given $n$ lineages, a particular set of exactly $k$ lineages merge one generation back. If a partial sweep of known amplitude $\beta$ took place, this requires the set of $k$ lineages to merge at this generation, with probability $\beta^k$, and that the other $n-k$ do not merge, with probability $(1 - \beta)^{n-k}$.

$$\lambda_n(k) = \rho \langle \beta^k (1 - \beta)^{n-k} \rangle, \\ = \int_0^1 \beta^k (1 - \beta)^{n-k} P(\beta) d\beta, \quad = \int_0^1 \beta^{k-2} (1 - \beta)^{n-k} \frac{\beta^2 P(\beta)}{\langle \beta^2 \rangle} d\beta. \qquad (B2)$$

This turns out to be the definition of the $\Lambda$-coalescent with $\Lambda(\beta) \propto \beta^2 P(\beta)$ [29, 30]. The $\Lambda$-coalescent is a general model for genealogies of multiple mergers. We mention two interesting subcases:

* if $P(\beta) = \delta(\beta)$ where $5$ is the Dirac distribution with all of the mass at $0$, the only possible merge is $k = 2$ and we recover the Kingman coalescent [30]. For this reason, we expect our coalescent to approach Kingman's when $\beta \ll 1$, which will be shown explicitly below.

* if $\Lambda(\beta)$ is uniform in [0,1], meaning $P(\beta) \propto \beta^{-2}$, we obtain the Bolthausen-Sznitman coalescent [30, 41], which is used to describe the genealogy of populations under strong selection [41, 43].

Finally, we derive a few more properties of the partial sweep coalescent and show the explicit link to Kingman's when $\beta \ll 1$. Using the $\gamma_n(k)$'s, we can compute the times $T_n$: the time during which exactly $n$ lineages are present in parallel in the genealogy. If there are $n$ lineages present, any
<figure><img src="image_13.png" /><figcaption>Figure S 4.</figcaption></figure>

Probability of a strictly monotonous trajectory in the random walk of the main text, as a function of $\beta$ (fixed) and the initial value $x_0$. The "exact" solution is obtained by numerically computing the product in **Equation B1** ☐ up to $t = 100$.
coalescence will lower the number of lineages below $n$. The time $T_n$ is thus exponentially distributed with rate $v/(n)$, where $v/(n)$ is the total rate of coalescence given $n$ lineages:

$$v(n) = \sum_{k=2}^{n} \rho \binom{n}{k} \lambda_n(k).$$

Since we have $\sum_{k=0}^{n} \rho \binom{n}{k} \lambda_n(k) = \rho(1 - \beta + \beta)^n = \rho$, we finally obtain

$$T_n^{-1} = \rho (1 - \langle (1 - \beta)^n \rangle - n \langle \beta(1 - \beta)^{n-1} \rangle). \quad (B3)$$

With $n\langle\beta\rangle \ll 1$, we now exactly recover the Kingman coalescent. For simplicity, we assume a constant $\beta$ and expand $T_n$ up to the second order in $n\beta$, to obtain

$$T_n^{-1} = \rho\beta^2 \frac{n(n-1)}{2} = \frac{n(n-1)}{2N_e}.$$

These are the times expected for the Kingman coalescent with population size $N_e = 1/\rho\beta^2$.

In the high $n$ limit, we also obtain $T_n \to \rho^{-1}$, since quantities of the type $(1 - \beta)^n$ vanish. This is expected as coalescences only take place when a partial sweep happens, with rate $\rho$. It is another qualitative difference with the Kingman coalescent: since $T_n \ge \rho^{-1}$ for all $n$, one must wait a time $\sim \rho^{-1}$ to observe the first coalescence even in large trees. The shortest branches will thus always be of order $\rho^{-1}$. In contrast, in the Kingman process, the shortest branches vanish when the number of lineages $n$ increases. This difference is clearly visible when looking at terminal branches of trees in **Figure S6** ☐.

## C. Supplementary figures
Figure S 5.

Average coalescence times $\langle T_n \rangle$ for a partial sweep coalescent with effective population size $N_e$ and a Kingman coalescent with population size $N_e$. For simplicity, a constant $\beta$ is used: **Left:** a high value $\beta = 0.25$; **Right:** a low value $\beta = 0.05$. For low $\beta$, the two coalescent processes are very similar until a high $n$. They considerably differ if $\beta$ is larger. Note that for the partial sweep process, $T_n$ never goes below $\rho^{-1}$.

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

Figure S 6.

Realisations of different coalescence processes for 30 lineages (leaves). **Left:** Partial sweep coalescent, with constant $\beta = 0.4$ and $\rho = 0.00625$ such that $N_e = (\rho\beta^2)^{-1} = 1000$. **Right:** Kingman coalescent with population size $N = N_e = 1000$.

<figure><img src="image_15.png" /></figure>
<figure><img src="image_16.png" /><figcaption>Figure S 7.</figcaption></figure>

Example of mutation frequency trajectories that are increasing up to a frequency of 0.5 for H3N2/HA influenza and the expiring fitness model. For the latter, parameters used are $\alpha = s_0 = 0.03$ and three values of $\rho$ to illustrate different clonal interference regimes. In each case, 10 randomly selected trajectories are plotted, with blue color indicating final loss and red final fixation.
References

[1] Roemer Cornelius *et al*. (2023) **SARS-CoV-2 evolution in the Omicron era** *Nature Microbiology* https://doi.org/10.1038/s41564-023-01504-w

[2] Petrova Velislava N., Russell Cohn A. (2018) **The evolution of seasonal influenza viruses** *Nature Reviews Microbiology* Nature Publishing Group **16**:47–60 https://doi.org/10.1038/nrmicro.2017.118

[3] Morris Dylan H., Gostic Katelyn M., Pompei Simone, Bedford Trevor, Luksza Marta, Neher Richard A., Grenfell Bryan T., Lässig Michael, McCauley John W. (2018) **Predictive Modeling of Influenza Shows the Promise of Applied Evolutionary Biology** *Trends in Microbiology* **26**:102–118 https://doi.org/10.1016/j.tim.2017.09.004

[4] Meijers Matthijs, Ruchnewitz Denis, Eberhardt Jan, Luksza Marta, Lässig Michael (2023) **Population immunity predicts evolutionary trajectories of SARS-CoV-2** *Cell* Elsevier **186**:5151–5164 https://doi.org/10.1016/j.cell.2023.09.022

[5] Tsimring Lev S., Levine Herbert, Kessler David A. (1996) **RNA Virus Evolution via a Fitness-Space Model** *Physical Review Letters* **76**:4440–4443 https://doi.org/10.1103/PhysRevLett.76.4440

[6] Rouzine Igor M., Wakeley John, Coffin John M. (2003) **The solitary wave of asexual evolution** *Proceedings of the National Academy of Sciences* Proceedings of the National Academy of Sciences **100**:587–592 https://doi.org/10.1073/pnas.242719299

[7] Desai Michael M, Fisher Daniel S (2007) **Beneficial Mutation-Selection Balance and the Effect of Linkage on Positive Selection** *Genetics* **176**:1759–1798 https://doi.org/10.1534/genetics.106.067678

[8] Neher Richard A. (2013) **Genetic Draft, Selective Interference, and Population Genetics of Rapid Adaptation** *Annual Review of Ecology, Evolution, and Systematics* **44**:195–215 https://doi.org/10.1146/annurev-ecolsys-110512-135920

[9] Luksza Marta, Lässig Michael (2014) **A predictive fitness model for influenza** *Nature* **507**:57–61 https://doi.org/10.1038/nature13087

[10] Neher Richard A, Russell Cohn A, Shraiman Boris I (2014) **Predicting evolution from the shape of genealogical trees** *eLife* **3** https://doi.org/10.7554/eLife.03568

[11] Huddleston John *et al*. (2020) **Integrating genotypes and phenotypes improves long-term forecasts of seasonal influenza A/H3N2 evolution** *eLife* eLife Sciences Publications, Ltd **9** https://doi.org/10.7554/eLife.60067

[12] Barrat-Charlaix Pierre, Huddleston John, Bedford Trevor, Neher Richard A (2021) **Limited Predictability of Amino Acid Substitutions in Seasonal Influenza Viruses** *Molecular Biology and Evolution* https://doi.org/10.1093/molbev/msab065

[13] Gupta Sunetra, Ferguson Neil, Anderson Roy (1998) **Chaos, Persistence, and Evolution of Strain Structure in Antigenically Diverse Infectious Agents** *Science* **280**:912–915 https://doi.org/10.1126/science.280.5365.912
[14] Gog Julia R., Grenfell Bryan T. (2002) **Dynamics and selection of many-strain pathogens**
*Proceedings of the National Academy of Sciences* Proceedings of the National Academy of Sciences **99**:17209-17214 https://doi.org/10.1073/pnas.252512799

[15] Yan Le, Neher Richard A, Shraiman Boris I (2019) **Phylodynamic theory of persistence, extinction and speciation of rapidly adapting pathogens** *eLife* eLife Sciences Publications, Ltd 8 https://doi.org/10.7554/eLife.44205

[16] Marchi Jacopo, Lässig Michael, Walczak Aleksandra M., Mora Thierry (2021) **Antigenic waves of virus-immune coevolution** *Proceedings of the National Academy of Sciences* Proceedings of the National Academy of Sciences **118** https://doi.org/10.1073/pnas.2103398118

[17] Chardes Victor, Mazzolini Andrea, Mora Thierry, Walczak Aleksandra M. (2023) **Evolutionary stability of antigenically escaping viruses** *Proceedings of the National Academy of Sciences* Proceedings of the National Academy of Sciences **120** https://doi.org/10.1073/pnas.2307712120

[18] Rouzine Igor M., Rozhnova Ganna (2018) **Antigenic evolution of viruses in host populations**
*PLOS Pathogens* Public Library of Science **14** https://doi.org/10.1371/journal.ppat.1007291

[19] Lee Juhye M *et al*. (2019) **Mapping person-to-person variation in viral mutations that escape polyclonal serum targeting influenza hemagglutinin** *eLife* eLife Sciences Publications, Ltd 8 https://doi.org/10.7554/eLife.49324

[20] Welsh Frances C. *et al*. (2023) **Age-dependent heterogeneity in the antigenic effects of mutations to influenza hemagglutinin** *bioRxiv*

[21] Good Benjamin H., Martis Stephen, Hallatschek Oskar (2018) **Adaptation limits ecological diversification and promotes ecological tinkering during the competition for substitutable resources** *Proceedings of the National Academy of Sciences* **115** https://doi.org/10.1073/pnas.1807530115

[22] Tikhonov Mikhail, Monasson Remi (2018) **Innovation Rather than Improvement: A Solvable High-Dimensional Model Highlights the Limitations of Scalar Fitness** *Journal of Statistical Physics* **172**:74-104 https://doi.org/10.1007/s10955-018-1956-6

[23] Kucharski Adam J., Lessler Justin, Cummings Derek A. T., Riley Steven (2018) **Timescales of influenza A/H3N2 antibody dynamics** *PLoS biology* **16** https://doi.org/10.1371/journal.pbio.2004974

[24] Gillespie John H (2000) **Genetic Drift in an Infinite Population: The Pseudohitchhiking Model** *Genetics* **155**:909-919 https://doi.org/10.1093/genetics/155.2.909

[25] Hallatschek Oskar (2018) **Selection-Like Biases Emerge in Population Models with Recurrent Jackpot Events** *Genetics* **210**:1053-1073 https://doi.org/10.1534/genetics.118.301516

[26] Kimura Motoo (1964) **Diffusion models in population genetics** *Journal of Applied Probability* **1**:177-232 https://doi.org/10.2307/3211856

[27] Schiffels Stephan, Szöllosi Gergely J, Mustonen Ville, Lässig Michael (2011) **Emergent Neutrality in Adaptive Asexual Evolution** *Genetics* **189**:1361-1375 https://doi.org/10.1534/genetics.iii.132027
[28] Strelkowa Natalja, Lässig Michael (2012) **Clonal Interference in the Evolution of Influenza Genetics** 192:671-682 https://doi.org/10.1534/genetics.H2.143396

[29] Schweinsberg Jason (2000) **Coalescents with Simultaneous Multiple Collisions** *Electronic Journal of Probability* Institute of Mathematical Statistics and Bernoulli Society 5:1-50 https://doi.org/10.1214/EJP.v5-68

[30] Berestycki Nathanaël (2009) **Recent progress in coalescent theory** *Ensaios Matemáticos* 16:1-193 https://doi.org/10.21711/217504322009/em161

[31] Bush Robin M., Bender Catherine A., Subbarao Kanta, Cox Nancy J., Fitch Walter M. (1999) **Predicting the Evolution of Human Influenza A** *Science* American Association for the Advancement of Science 286:1921-1925 https://doi.org/10.1126/science.286.5446.1921

[32] Smith Derek J., Lapedes Alan S., Jong Jan C. De, Bestebroer Theo M., Rimmelzwaan Guus F., Osterhaus Albert D. M. E., Fouchier Ron A. M. (2004) **Mapping the Antigenic and Genetic Evolution of Influenza Virus** *Science* 305:371-376 https://doi.org/10.1126/science.1097211

[33] Bhatt Samir, Holmes Edward C., Pybus Oliver G. (2011) **The Genomic Rate of Molecular Adaptation of the Human Influenza A Virus** *Molecular Biology and Evolution* 28:2443-2451 https://doi.org/10.1093/molbev/msr044

[34] Neher Richard A., Bedford Trevor, Daniels Rodney S., Russell Cohn A., Shraiman Boris I. (2016) **Prediction, dynamics, and visualization of antigenic phenotypes of seasonal influenza viruses** *Proceedings of the National Academy of Sciences* Proceedings of the National Academy of Sciences 113 https://doi.org/10.1073/pnas.1525578113

[35] Kistler Kathryn E., Bedford Trevor (2023) **An atlas of continuous adaptive evolution in endemic human viruses** *Cell Host & Microbe* https://doi.org/10.1016/j.chom.2023.09.012

[36] Bhatt S. et al. (2013) **The evolutionary dynamics of influenza A virus adaptation to mammalian hosts** *Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences* 368 https://doi.org/10.1098/rstb.2012.0382

[37] Fonville J. M. et al. (2014) **Antibody landscapes after influenza virus infection or vaccination** *Science* American Association for the Advancement of Science 346:996-1000 https://doi.org/10.1126/science.1256427

[38] Pelletier F., Garant D., Hendry A.p. (2009) **Eco-evolutionary dynamics** *Philosophical Transactions of the Royal Society B: Biological Sciences* Royal Society 364:1483-1489 https://doi.org/10.1098/rstb.2009.0027

[39] Shu Yuelong, McCauley John (2017) **GISAID: Global initiative on sharing all influenza data - from vision to reality** *Eurosurveillance* European Centre for Disease Prevention and Control 22 https://doi.org/10.2807/1560-7917.ES.2017.22.13.30494

[40] Rice Daniel P, Good Benjamin H, Desai Michael M (2015) **The Evolutionarily Stable Distribution of Fitness Effects** *Genetics* 200:321-329 https://doi.org/10.1534/genetics.114.173815

[41] Bolthausen E., Sznitman A.-S. (1998) **On Ruelle's Probability Cascades and an Abstract Cavity Method** *Communications in Mathematical Physics* 197:247-276 https://doi.org/10.1007/s002200050450
[42] Brunet E., Derrida B., Mueller A. H., Munier S. (2007) **Effect of selection on ancestry: An exactly soluble case and its phenomenological generalization** *Physical Review E* **American Physical Society** **76** https://doi.org/10.1103/PhysRevE.76.041104

[43] Neher R. A., Hallatschek O. (2013) **Genealogies of rapidly adapting populations** *Proceedings of the National Academy of Sciences* **110**:437–442 https://doi.org/10.1073/pnas.1213113110

**Editors**

Reviewing Editor

**Armita Nourmohammad**

University of Washington, Seattle, United States of America

Senior Editor

**Aleksandra Walczak**

École Normale Supérieure - PSL, Paris, France

**Reviewer #1 (Public review):**

In this work, the authors study the dynamics of fast-adapting pathogens under immune pressure in a host population with prior immunity. In an immunologically diverse population, an antigenically escaping variant can perform a partial sweep, as opposed to a sweep in a homogeneous population. In a certain parameter regime, the frequency dynamics can be mapped onto a random walk with zero mean, which is reminiscent of neutral dynamics, albeit with differences in higher order moments. Next, they develop a simplified effective model of time dependent selection with expiring fitness advantage, and posit that the resulting partial sweep dynamics could explain the behaviour of influenza trajectories empirically found in earlier work (Barrat-Charlaix et al. Molecular Biology and Evolution, 2021). Finally, the authors put forward an interesting hypothesis: the mode of evolution is connected to the age of a lineage since ingression into the human population. A mode of meandering frequency trajectories and delayed fixation has indeed been observed in one of the long-established subtypes of human influenza, albeit so far only over a limited period from 2013 to 2020. The paper is overall interesting and well-written.

In the revised version, the authors have addressed questions on the role of clonal interference by new simulations in the SI, clarified the connection between the SIR model and vanishing-fitness models, and placed their analysis into the broader context of consumer resource dynamics.

However, the general conclusion, as stated in the abstract, that variant trajectories become unpredictable as a consequence of the SIR dynamics remains somewhat misleading. Two aspects contribute to this problem. (1) The empirical observation of quasi-neutrality, i.e. the absence of a net frequency increase inferred as an average of many trajectories at intermediate frequencies, does not imply that individual trajectories are neutral (i.e., fully stochastic and unpredictable) over the time span of observation. Rather, it just says that some have a positive and some have a negative selection coefficient over that time span. (2) As stated by the authors, the observation of average quasi-neutrality is indeed incompatible with the travelling wave model, where initially successful new variants are assumed to retain a fixed, positive selection coefficient from origination to fixation. This observation also limits predictions by extrapolation, where a positive selection coefficient inferred at small frequency is assumed to remain the same at later times and higher frequencies. However, predictions derived from Gog and Grenfell's multi-strain SIR model, as used by several authors, do not make the assumption of fixed selection coefficients and incorporate
trajectory-specific, time-dependent expiration effects into their model predictions. This
distinction remains blurred throughout the text of the paper.

https://doi.org/10.7554/eLife.97350.2.sa2

Reviewer #3 (Public review):

In this work the authors present a multi-strain SIR model in which viruses circulate in a heterogeneous population with different groups characterized by different cross-immunity structures. They reformulate the qualitative features of these SIR dynamics as a random walk characterized by new variants saturating at intermediate frequencies. Then they recast their microscopic description to an effective formalism in which viral strains lose fitness independently from one another. They study several features of this process numerically and analytically, such as the average variants frequency, the probability of fixation, and the coalescent time. They compare qualitatively the dynamics of this model to variants dynamics in RNA viruses such as flu and SARS-CoV-2

The idea that vanishing fitness mechanisms that produce partial sweeps may explain important features of flu evolution is very interesting. Its simplicity and potential generality make it a powerful framework. As noted by the authors, this may have important implications for predictability of virus evolution and such a framework may be beneficial when trying to build predictive models for vaccine design. The vanishing fitness model is well analyzed and produces interesting structures in the strains coalescent. Even though the comparison with data is largely qualitative, this formalism would be helpful when developing more accurate microscopic ingredients that could reproduce viral dynamics quantitatively.

This general framework has the potential to be more universal than human RNA viruses, in
situations where invading mutants would saturate at intermediate frequencies.

The qualitative connection between the coarse-grained features of these vanishing fitness
dynamics and structured SIR processes offers additional intuition relevant to host-pathogens
interactions, although as noted by the authors other ecological processes could drive similar
evolutionary patterns. The additions in the revised manuscript, substantiating more
thoroughly the connection between the SIR and the vanishing fitness description, are
important to better appreciate the scope of the work.

https://doi.org/10.7554/eLife.97350.2.sa1

Author response:

The following is the authors' response to the original reviews.

Response to reviewers

We thank the Editor and the Reviewers for their constructure review. In the light of this
feedback, we have made a number of changes and additions to the manuscript, that we think
improved the presentation and hopefully address the majority of the concerns by the
reviewers.

Main changes:

* We added a new SI section (B1) with a population dynamics simulation in the high clonal interference regime and without expiring fitness (see R1: (1)).
* We added a new SI section (A9) with the derivation of the equilibrium state of our SIR model in the case of *M* immune groups and in the limit ε → 0 (see R1: (5)).

* The text of the section Abstraction as “expiring” fitness advantage has been modified.

* We added a new SI section (A4) describing the links between parameters of the “expiring fitness” and SIR models.

All three reviewers had concerns about the relation between our SIR model and the “expiring fitness” model, that we hope will be addressed by the last two items listed above. In particular, we would like to underline the following points:

* The goal of our SIR model is to give a mechanistic explanation of partial sweeps using traditional epidemiological models. While ecological models (e.g. consumer resource) can give rise to the same phenomenology, we believe that in the context of host-pathogen interaction it is relevant to explicitly show that SIR models can result in partial sweeps.

* The expiring fitness model is mainly an effective model: it reproduces some qualitative features of the SIR but does not quantitatively match all aspects of the frequency dynamics in SIR models.

* It is possible to link the parameters of the SIR (*a*,*,b*,f) and expiring fitness (*s*,*x*,*v*) models at the beginning of the invasion of the variant (new SI section A4). However, the two models also differ in significant ways (the SIR model can for example oscillate, while the effective model can not). The correspondence of quantities like the initial invasion rate and the ‘expiration rate’ of fitness effects is thus only expected to hold for some time after the emergence of a novel variant.

**Public reviews:**

**Reviewer 1:**

Summary In this work, the authors study the dynamics of fast-adapting pathogens under immune pressure in a host population with prior immunity. In an immunologically diverse population, an antigenically escaping variant can perform a partial sweep, as opposed to a sweep in a homogeneous population. In a certain parameter regime, the frequency dynamics can be mapped onto a random walk with zero mean, which is reminiscent of neutral dynamics, albeit with differences in higher order moments. Next, they develop a simplified effective model of time dependent selection with expiring fitness advantage, and posit that the resulting partial sweep dynamics could explain the behaviour of influenza trajectories empirically found in earlier work (Barrat-Charlaix et al. Molecular Biology and Evolution, 2021). Finally, the authors put forward an interesting hypothesis: the mode of evolution is connected to the age of a lineage since ingression into the human population. A mode of meandering frequency trajectories and delayed fixation has indeed been observed in one of the long-established subtypes of human influenza, albeit so far only over a limited period from 2013 to 2020. The paper is overall interesting and well-written. Some aspects, detailed below, are not yet fully convincing and should be treated in a substantial revision.

We thank the reviewer for their constructive criticism. The deep split in the A/H3N2 HA segment from 2013 to 2020 is indeed the one of the more striking examples of such meandering frequency dynamics in otherwise rapidly adapting populations. But the up and down of H1N1pdm clade 5a.2a.1 in recent years might be a more recent example. We argue that such meandering dynamics might be a common contributor to seasonal influenza dynamics, even if it only spans 3-6 years.
(1) The quasi-neutral behaviour of amino acid changes above a certain frequency (reported in Fig. 3), which is the main overlap between influenza data and the authors' model, is not a specific property of that model. Rather, it is a generic property of travelling wave models and more broadly, of evolution under clonal interference (Rice et al. Genetics 2015, Schiffels et al. Genetics 2011). The authors should discuss in more detail the relation to this broader class of models with emergent neutrality. Moreover, the authors' simulations of the model dynamics are performed up to the onset of clonal interference $\rho/s_0 = 1$ (see Fig. 4). Additional simulations more deeply in the regime of clonal interference (e.g. $\rho/s_0 = 5$) show more clearly the behaviour in this regime.

We agree with the reviewer that we did not discuss in detail the effects of clonal interference on quasi-neutrality and predictability. As suggested, we conducted additional simulations of our population model in the regime of high clonal interference ($\rho/s_0 \gg 1$) and without expiring fitness effects. The results are shown in a new section of the supplementary information. These simulations show, as expected, that increasing clonal interference tends to decrease predictability: the fixation probability of an adaptive mutation found at frequency $x$ moves closer to $x$ as $\rho$ increases. However, even in a case of strong interference $\rho/s_0 = 32$, pfix remains significantly different from the neutral expectation. We conclude from this that while it is true that dynamics tend to quasi-neutrality in the case of strong interference, this effect alone is unlikely to explain observations of H3N2 influenza dynamics. In our previous publication (BarratCharlaix et al, MBE, 2021) we have also investigated the effect of epistatic interactions between mutations, along side strong clonal interference. We concluded that, while most of these processes make evolution less predictable and push pfix towards the diagonal, it is hard to reproduce the empirical observations with realistic parameters. The "expiring fitness" model, however, produces this quite readily.

But there are qualitative differences between quasi-neutrality in traveling wave models and the expiring fitness model. In the traveling wave, a genotype carrying an adaptive mutation is always fitter than if it didn't carry the mutation. Quasi-neutrality emerges from the accumulation of fitness variation at other loci and the fact that the coalescence time is not much bigger than the inverse selection coefficient of the mutation. In the expiring fitness model, the selective effect of the mutation itself goes away with time. We now discuss the literature on quasi-neutrality and cite Rice et al. 2015 and Schiffels et al. 2011.

> In this context, I also note that the modelling results of this paper, in particular the stalling of frequency increase and the decrease in the number of fixations, are very similar to established results obtained from similar dynamical assumptions in the broader context of consumer resource models; see, e.g., Good et al. PNAS 2018. The authors should place their model in this broader context.

We thank the reviewer for pointing out the link between consumer resource models and our work. We further strengthened our discussion of the similarity of the phenomenology to models typically used in ecology and made an effort to highlight the link between consumer-resource models and ours in the introduction and in the part on the SIR model.

(2) The main conceptual problem of this paper is the inference of generic non-predictability from the quasi-neutral behaviour of influenza changes. There is no question that new mutations limit the range of predictions, this problem being most important in lineages with diverse immune groups such as influenza A(H3N2). However, inferring generic non-predictability from quasi-neutrality is logically problematic because predictability refers to individual trajectories, while quasi-neutrality is a property obtained by averaging over many trajectories (Fig. 3). Given an SIR dynamical model for trajectories, as employed here and elsewhere in the literature, the up and down of
individual trajectories may be predictable for a while even though allele frequencies do not increase on average. The authors should discuss this point more carefully.

We agree with the reviewer that the deterministic SIR model is of course predictable. Similarly, a partial sweep is predictable. But we argue that expiring fitness makes evolution less predictable in two ways: (i) When a new adaptive mutation emerges and rises in frequency, we typically don't know how rapidly its fitness effect is 'expiring'. Thus even if we can measure its instantaneous growth rate accurately, we can't predict its fate far into the future. (ii) Compared to the situation where fitness effects are not expiring, time to fixation is longer and there are more opportunities for novel mutations to emergence and change the course of the trajectory. We have tried to make this point clearer in the manuscript.

(3) To analyze predictability and population dynamics (section 5), the authors use a Wright-Fisher model with expiring fitness dynamics. While here the two sources of the emerging neutrality are easily tuneable (expiring fitness and clonal interference), the connection of this model to the SIR model needs to be substantiated: what is the starting selection $s_0$ as a function of the SIR parameters ($f, b, M, \varepsilon$), the selection decay $v = v(f, b, M, \varepsilon, \gamma)$? This would enable the comparison of the partial sweep timing in both models and corroborate the mapping of the SIR onto the simplified W-F model. In addition, the authors' point would be strengthened if the SIR partial sweeps in Fig. 1 and Fig. 2 were obtained for a combination of parameters that results in a realistic timescale of partial sweeps.

We added a new section to the SI (A4) that relates the parameters of the SIR and expiring fitness models. In particular, we compute the initial growth rate $s_0$ and a proxy for the fitness expiry rate $v$ as a function of the SIR parameters $\alpha, \gamma, f, b, M$, at the instant where the variant is introduced. The initial growth rate depends primarily on the degree of immune escape $f$, while the expiration rate $v$ is related to incidence $I_{wt} + I_{m}$. However, as both models have fundamentally different dynamics, these relations are only valid on time scales shorter than potential oscillations of the SIR model. Beyond that, the connection between the models is mostly qualitative: both rely on the fact that growth rate of a strain diminishes when the strain becomes more frequent, and give rise to partial sweeps.

In Figure 1, the time it takes a partial sweep to finish is roughly 100–200 generations (bottom right panel). If we consider H3N2 influenza and take one generation to be one week, this corresponds to a sweep time of 2 to 4 years, which is slightly slower but roughly in line with observations for selective sweeps. This time is harder to define if oscillatory dynamics takes place (middle right panel), but the time from the introduction of the mutant to the peak frequency is again of about 4 years. The other parameters of the model correspond to a waning time of 200 weeks and immune escape on the order of 20-30% change in susceptibility.

**Reviewer 2:**

**Summary**

This work addresses a puzzling finding in the viral forecasting literature: high-frequency viral variants evince signatures of neutral dynamics, despite strong evidence for adaptive antigenic evolution. The authors explicitly model interactions between the dynamics of viral adaptations and of the environment of host immune memory, making a solid theoretical and simulation-based case for the essential role of host-pathogen eco-evolutionary dynamics. While the work does not directly address improved data-driven viral forecasting, it makes a valuable conceptual contribution to the key dynamical ingredients (and perhaps intrinsic limitations) of such efforts.
Strengths

This paper follows up on previous work from these authors and others concerning the problem of predicting future viral variant frequency from variant trajectory (or phylogenetic tree) data, and a model of evolving fitness. This is a problem of high impact: if such predictions are reliable, they empower vaccine design and immunization strategies. A key feature of this previous work is a "traveling fitness wave" picture, in which absolute fitnesses of genotypes degrade at a fixed rate due to an advancing external field, or "degradation of the environment". The authors have contributed to these modeling efforts, as well as to work that critically evaluates fitness prediction (references 11 and 12). A key point of that prior work was the finding that fitness metrics performed no better than a baseline neutral model estimate (Hamming distance to a consensus nucleotide sequence). Indeed, the apparent good performance of their well-adopted "local branching index" (LBI) was found to be an artifact of its tendency to function as a proxy for the neutral predictor. A commendable strength of this line of work is the scrutiny and critique the authors apply to their own previous projects. The current manuscript follows with a theory and simulation treatment of model elaborations that may explain previous difficulties, as well as point to the intrinsic hardness of the viral forecasting inference problem.

This work abandons the mathematical expedience of traveling fitness waves in favor of explicitly coupled eco-evolutionary dynamics. The authors develop a multi-compartment susceptible/infected model of the host population, with variant cross-immunity parameters, immune waning, and infectious contact among compartments, alongside the viral growth dynamics. Studying the invasion of adaptive variants in this setting, they discover dynamics that differ qualitatively from the fitness wave setting: instead of a succession of adaptive fixations, invading variants have a characteristic "expiring fitness": as the immune memories of the host population reconfigure in response to an adaptive variant, the fitness advantage transitions to quasi-neutral behavior. Although their minimal model is not designed for inference, the authors have shown how an elaboration of host immunity dynamics can reproduce a transition to neutral dynamics. This is a valuable contribution that clarifies previously puzzling findings and may facilitate future elaborations for fitness inference methods.

The authors provide open access to their modeling and simulation code, facilitating future applications of their ideas or critiques of their conclusions.

We thank the reviewer for their summary, assessement, and constructive critique.

(1) The current modeling work does not make direct contact with data. I was hoping to see a more direct application of the model to a data-driven prediction problem. In the end, although the results are compelling as is, this disconnect leaves me wondering if the proposed model captures the phenomena in detail, beyond the qualitative phenomenology of expiring fitness. I would imagine that some data is available about cross-immunity between strains of influenza and sarscov2, so hopefully some validation of these mechanisms would be possible.

We agree with the reviewer that quantitatively confronting our model with data would be very interesting. Unfortunately, most available serological data for influenza and SARS-CoV-2 is obtained using post-infection sera from previously naive animal models. To test our model, we would require human serology data, ideally demographically resolved, and a way to link serology to transmission dynamics. Furthermore, our model is mostly an explanation for qualitative features of variant dynamics and their apparent lack of predictability. We therefore considered that quantitative validation using data is out of scope of this work.
(2) After developing the SIR model, the authors introduce an effective "expiring fitness" model that avoids the oscillatory behavior of the SIR model. I hoped this could be motivated more directly, perhaps as a limit of the SIR model with many immune groups. As is, the expiring fitness model seems to lose the eco-evolutionary interpretability of the SIR model, retreating to a more phenomenological approach. In particular, it's not clear how the fitness decay parameter $v$ and the initial fitness advantage $\delta_0$ relate to the key ecological parameters: the strain cross-immunity and immune group interaction matrices.

The expiring fitness model emerges as a limiting case, at least qualitatively, of the SIR model when growth rate of the new variant is small compared to the waning rate and the SIR model does not oscillate. This can be readily achieved by many immune groups, which reconciles the large effect of many escape mutations and the lack of oscillation by confining the escape to some fraction of the population. Beyond that, the expiring fitness model is mainly an effective model that allows us to study the consequences of partial sweeps on predictability on long timescales. As stated in the "Main changes" section at the start of this reply, we added an SI section which links parameters of the two models. However, we underline the fact that beyond the phenomenon of partial sweeps, the dynamics of the two are different.

**Reviewer 3:**

**Summary**

In this work the authors start presenting a multi-strain SIR model in which viruses circulate in an heterogeneous population with different groups characterized by different cross-immunity structures. They argue that this model can be reformulated as a random walk characterized by new variants saturating at intermediate frequencies. Then they recast their microscopic description to an effective formalism in which viral strains lose fitness independently from one another. They study several features of this process numerically and analytically, such as the average variants frequency, the probability of fixation, and the coalescent time. They compare qualitatively the dynamics of this model to variants dynamics in RNA viruses such as flu and SARS-CoV-2.

**Strengths**

The idea that a vanishing fitness mechanism that produces partial sweeps may explain important features of flu evolution is very interesting. Its simplicity and potential generality make it a powerful framework. As noted by the authors, this may have important implications for predictability of virus evolution and such a framework may be beneficial when trying to build predictive models for vaccine design. The vanishing fitness model is well analyzed and produces interesting structures in the strains coalescent. Even though the comparison with data is largely qualitative, this formalism would be helpful when developing more accurate microscopic ingredients that could reproduce viral dynamics quantitatively. This general framework has a potential to be more universal than human RNA viruses, in situations where invading mutants would saturate at intermediate frequencies.

We thank the reviewer for their positive remarks and constructive criticism below.

**Weaknesses**

The authors build the narrative around a multi-strain SIR model in which viruses circulate in an heterogeneous population, but the connection of this model to the rest of the paper is not well supported by the analysis. When presenting the random walk coarse-grained description in section 3 of the Results, there is no quantitative relation
between the random walk ingredients importantly $P(\beta)$ - and the SIR model, just a qualitative reasoning that strains would initially grow exponentially and saturate at intermediate frequencies. So essentially any other microscopic description with these two features would give rise to the same random walk.

As also highlighted in the response to other reviewers, we now discuss how the parameter of the SIR model are related to the initial growth rate and the ‘expiration’ rate of the effective model. While the phenomenology of the SIR model is of course richer, this correspondence describes its overdamped limit qualitatively well.

> Currently it's unclear whether the specific choices for population heterogeneity and cross-immunity structure in the SIR model matter for the main results of the paper. In section 2, it seems that the main effect of these ingredients are reduced oscillations in variants frequencies and a rescaled initial growth rate. But ultimately a homogeneous population would also produce steady state coexistence between strains, and oscillation amplitude likely depends on parameters choices. Thus a homogeneous population may lead to a similar coarse-grained random walk.

The reviewer is correct that the primary effects of using many immune groups is to slow down the increase of novel variant, which in turn dampens the oscillations. Having multiple immune groups widens the parameter space in which partial sweeps without dramatic oscillations are observed. For slow sweeps, similar dynamics are observed in a homogeneous population.

> Similarly, it's unclear how the SIR model relates to the vanishing fitness framework, other than on a qualitative level given by the fact that both descriptions produce variants saturating at intermediate frequencies. Other microscopic ingredients may lead to a similar description, yet with quantitative differences.

Both of these points were also raised by other reviewers and we agree that it is worth discussing them at greater length. We now discuss how the parameters of the ‘expiring fitness’ model relate to those of the SIR. We also discuss how other models such as ecological models give rise to similar coarse grained models.

> At the same time, from the current analysis the reader cannot appreciate the impact of such a mean field approximation where strains lose fitness independently from one another, and under what conditions such assumption may be valid.

In the SIR model, the rate at which strains lose fitness does depend on the precise state of the host population through the quantities $S_m$ and $S_w$, which is apparent in equation (A27) of the new SI section. The fact that a new variant shifts the equilibrium frequencies of previous strains in a proportional way is valid if the “antigenic space” is of very high dimensions, as explained in section Change in frequency when adding subsequent strains of the SI. It would indeed be interesting to explore relaxations of this assumption by considering a larger class of cross immunity matrices $K$. However, in the expiring fitness model, the fact that strains lose fitness independently from each other is a necessary simplification.

> In summary, the central and most thoroughly supported results in this paper refer to a vanishing fitness model for human RNA viruses. The current narrative, built around the SIR model as a general work on host-pathogen eco-evolution in the abstract, introduction, discussion and even title, does not seem to match the key results and may mislead readers. The SIR description rather seems one of the several possible models, featuring a negative frequency dependent selection, that would produce coarse-grained dynamics qualitatively similar to the vanishing fitness description analyzed here.
We have revised the text throughout to make the connections between the different parts of the manuscript, in particular the SIR model and the expiring fitness model, clearer. We agree that the phenomenology of the expiring fitness model is more general than the case of human RNA viruses described by the SIR model, but we think this generality is an attractive feature of the coarse-graining, not a shortcoming. Indeed, other settings with negative frequency dependent selection or eco-systems that adapt on appropriate time scale generate similar dynamics.

**Recommendations for the authors:**

**Reviewer 1:**

(4) *Line 74: what does fitness mean?*

Many population dynamics models, including ones used for viral forecasting, attach a scalar fitness to each strain. The growth rate of each strain is then computed by subtracting the average population fitness to the strain's fitness. In this sentence, fitness is intended in this way.

(5) *Fig. 1: The equilibrium frequency in the middle and bottom rows is hardly smaller than the equilibrium frequency in the top row for one immune group. This is surprising since for M=10, the variant escapes in only 1/10th of the population, which naively should impact the equilibrium frequency more strongly. Could the authors comment on this?*

This is indeed non-trivial, and a hand-waving argument can be made by considering the extreme case $ε = 0$. The variant is then completely neutral for the immune groups $i > 1$, and would be at equilibrium at any frequency in these immune groups. Its equilibrium frequency is then only determined by group 1, which is the only one breaking degeneracy. For $ε > 0$ but small, we naturally expect a small deviation from the $ε = 0$ case and thus $β$ should only change slightly.

A more rigorous argument with a mathematical proof in the case $ε = 0$ is now given in section A4 of the supplementary information.

(6) *Fig. 1: In the caption, it is stated that the simulations are performed with $ε = 0.99$. Is this a typo? It seems that it should be $ε = 0.01$, as in and just below equation (7).*

This was indeed a typo. It is now fixed.

(7) *Fig. 3: The data analysis should be improved. In order to link the average frequency trajectories to standard population genetics of conditional fixation probabilities, the focal time should always be the time where the trajectory crosses the threshold frequency for the first time. Plotting some trajectories from a later time onwards, on their downward path destined to loss, introduces a systematic bias towards negative clonal interference (for these trajectories, the time between the first and the second crossing of the threshold frequency is simply omitted). The focal time of first crossing of the threshold frequency can easily be obtained, e.g., by linear interpolation of the trajectory between subsequent time points of frequency evalution. In light of the modified procedure, the statements on the on the inertia of the trajectories after crossing $x^*$ (line 356) should be re-examined.*

The way we process the data is already in line with the suggestions of the reviewer. In particular, we use as focal time the first time at which a trajectory is found in the threshold
frequency bin. Trajectories that are never seen in the bin because of limited time-resolution are simply ignored.

In Fig. 3, there are no trajectories that are on their downward path at the focal time and when crossing the threshold frequency. Our other work on predictability of flu Barrat-Charlaix et al. (2021) has a similar figure, which maybe created confusion.

(8) Fig. 4: authors write $\alpha/\sigma_0$ in the figure, but should be $\nu/\sigma_0$.

Fixed.

(9) Line 420: authors refer to the blue curve in panel B as the case with strong interference. However, strong interference is for higher $\rho/\sigma_0$, that is panel D (see point 1).

Fixed.

(10) Line 477: typo "there will a variety of mutations".

Fixed.

**Reviewer 2:**

Should $\alpha$ be $\nu$ in Figure 4 legends?

Thank you very much for spotting this error. We fixed it.

Equations 4-5 could be further simplified.

We factorised the $I$ term in equation 4. In equation 5, we prefered to keep the $1 - \delta/\alpha$ term as this quantity appears in different calculations concerning the model. For instance, $S = \delta/\alpha$ at equilibrium.

The sentence before equation 8 references $P_\delta(\beta)$, but this wasn't previously introduced.

We now introduce $Pb_\eta$ at the beginning of the section Ultimate fate of the variant.

In the last paragraph of page 12, "monotonously" maybe should be "monotonically".

Fixed.

For the supplement section B, you might want a more descriptive title than "other".

We renamed this section to Expiring fitness model and random walk.

**Reviewer 3:**

To expand on my previous comments, my main concerns regard the connection of section 2 and the SIR model with the rest of the paper.

In the first paragraph of page 9 the authors argue that a stochastic version of the SIR model would lead to different fixation dynamics in homogeneous vs heterogeneous populations due to the oscillations. This paragraph is quite speculative, some numerical simulations would be necessary to quantitatively address to what extent these two scenarios actually differ in a stochastic setting, and how that depends on parameters.
Likewise, the connection between the SIR model, the random walk coarse-grained description and the vanishing fitness model can be investigated through numerical simulations of a stochastic SIR given the chosen population and cross-immunity structures with i.e. 10-20 strains. This would allow for a direct comparison of individual strain dynamics rather than the frequency averages, as well as other scalar properties such as higher moments, coalescent, and fixation probability once reaching a given frequency. It would also be possible to characterize numerically the SIR P[beta] bridging the gap with the random walk description. It's not obvious to me that the SIR P[beta] would not depend on the population size in the presence of birth-death stochasticity, potentially changing the moments scalings. I appreciate that such simulations may be computationally expensive, but similar numerical studies have been performed in previous phylodynamics works so it shouldn't be out of reach.

An alternative, the authors should consider re-centering the narrative directly on the random walk of the vanishing fitness model, mentioning the SIR more briefly as a possible qualitative way to get there. Either way the authors should comment on other ways in which this coarse-grained dynamics could arise.

In the vanishing fitness model, where variants fitnesses are independent, is an infinite dimensional antigenic space implicitly assumed? If that's the case, it should be explained in the main text.

A long simulation of the SIR model would indeed be interesting, but is numerically demanding and our current simulation framework doesn't scale well for many strains and susceptibilities. We thus refrained from adding extensive simulations.

In Figure 2B of the main text, the simulation with 7 strains illustrates the qualitative match between the expiring fitness and the SIR model. However, it is clearly not long enough to discuss statistical properties of the corresponding random walk. Furthermore, we do not expect the individual strain dynamics of the SIR and expiring fitness models to match. The latter depends on few parameters (a, s0), while the former depends on the full state of the host population and of the previous variants.

In the section linking the parameters of the two models, we now discuss the distribution P(β) of the SIR model for two strains and a specific choice of distribution for the cross immunity b and f.

**Minor comments:**

There is some back and forth in the writing. For instance, when introducing the model, *Cij* is first defined as 1/ *M*, then a few paragraphs later the authors introduce that in another limit *Cij* is just much higher than any *Cij*, and finally they specify that the former is the fast mixing scenario.

Another example is in section 2, in the first paragraph they put forward that heterogeneity and crossimmunity have different impacts on the dynamics, but the meaning attributed to these different ingredients becomes clear only a while later after the homogeneous population analysis. Uniforming the writing would make it easier for the reader to follow the authors' train of thought.

We removed the paragraph below Equation (1) mentioning the $C_{ij} = 1/M$ case, which we hope will linearize the writing.

When mentioning geographical structure, why would geography affect how immunity sees pairs of viral strains (differences in K)?
Geographic structure could influence cross-immunity because of exposure histories of hosts. For instance in the case of influenza, different geographical regions do not have the same dominating strains in each season, and hosts from different regions may thus build up different immunity.

> In the current narrative there are some speculations about non-scalar fitness, especially in section 2. The heterogeneity in this section does not seem so strong to produce a disordered landscape that defies the notion of scalar fitness in the same way some complex ecological systems do. A more parsimonious explanation for the coexistence dynamics observed here may be a negative frequency dependent selection.

Our language here was not very precise and we agree that the phenomenology we describe is related to that of frequency dependent selection (mediated by via immunity of the host population that integrates past frequencies). Traveling wave models typically use fitness function that are independent of the population distribution and only account for the evolution via an increasing average fitness. We have made discussion more accurate by stating that we consider a case where fitness depends explicitly on present and past population composition, which includes the case of negative frequency dependent selection.

> I don't understand the comparison with genetic drift (typo here, draft) in the last paragraph of section 3 given that there is no stochasticity in growth death dynamics.

We compare the random walk to genetic drift because of the expression of the second moment of the step size. The genetic draft has the same functional form. If one defines the effective population size as in the text, the drift due to random sampling of alleles (neutral drift) and the changes in strain frequency in our model have the same first and second moments. The stochasticity here does not come from the dynamics, which are indeed deterministic, but from the appearance of new mutations (variants) on backgrounds that are randomly sampled in the population. This latter property is shared with genetic draft.

> In the vanishing fitness model, I think the reader would benefit from having *P(s)* in the main text, and it should be made more clear what simulations assume what different choice of *P(s)*.

We added the expression of *P(s)* in the main text. Simulations use the value *s*0 = 0.03, which we added in the caption of Figure 4.

> When comparing the model and data, is the point that COVID is not reproduced due to clonal interference? It seems from the plot that flu has clonal interference as well though. Why is that negligible?

A similar point has been raised by the first reviewer (see R1-(1)). Clonal interference is not negligible, but we find it to be insufficient to explain the observations made for H3N2 influenza, namely the lack of inertia of frequency trajectories or the probability of fixation. This is shown in the new section (B1) of the SI. Both SARS-CoV-2 and H3N2 influenza experience clonal interference, but the former is more predictable than the latter. Our point is that expiring fitness effects should be stronger in influenza because of the higher immune heterogeneity of the host population, making it less predictable than SARS-CoV-2.

> Does the fixation probability as a function of frequency threshold match the flu data for some parameter sets?
For H3N2 influenza, the fixation probability is found to be equal to the threshold frequency (see Barrat-Charlaix MBE 2021, also indirectly visible from Fig. 3). In Figure 4, we obtain that either a high expiry rate or intermediate expiry rates and clonal interference regimes match this observation.

> It would be instructive to see examples of the individual variant dynamics of the vanishing fitness model compared to the presented data.

We added an extra SI figure (S7) showing 10 randomly selected trajectories of individual variants in the case of H3N2/HA influenza and for the expiring fitness model with different parameter choices.

> Figure 4E has no colorbar label. The reader shouldn't have to look for what that means in the bottom of the SIs. In panels A and B the label should be $v$, not $\alpha$. Same thing in most equations of page 42.

We added the colorbar label to the figure and also updated the caption: a darker color corresponds to a higher probability of sweeps to overlap. We fixed the $v-\alpha$ confusion in the SI and in the caption of the figure.

https://doi.org/10.7554/eLife.97350.2sa0