# Eco-evolutionary dynamics of adapting pathogens and host immunity

Pierre Barrat-Charlaix  Biozentrum, Universität Basel, Switzerland  Swiss
Institute of Bioinformatics, Switzerland  DISAT, Politecnico di Torino, Italy
Richard A. Neher  Correspondence to: Richard Neher, Biozentrum, Spitalstrasse
41, 4056, Basel, Switzerland<EMAIL_ADDRESS>Biozentrum, Universität
Basel, Switzerland  Swiss Institute of Bioinformatics, Switzerland

###### 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.

first keyword, second keyword, third keyword

###  .1 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 (Roemer et al., 2023;
Petrova and Russell, 2018), 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 (Morris et al., 2018; Meijers et al., 2023).

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 (Tsimring et al., 1996). 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_ Rouzine et al. (2003);
Desai and Fisher (2007); Neher (2013). 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 Luksza and
Lässig (2014); Neher et al. (2014); Huddleston et al. (2020).

In terms of modeling, a drawback 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 Meijers et al. (2023) or explain why prediction
is difficult (Barrat-Charlaix et al., 2021).

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 Gupta et al. (1998); Gog and Grenfell
(2002). 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 Yan et al. (2019); Marchi et al. (2021); Chardès et al. (2023); Rouzine
and Rozhnova (2018). 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 Lee et al. (2019); Welsh et al. (2023). 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 model is
reminiscent of ecological systems, where adaptation by one species shifts the
global equilibrium and distribution of other species but does not necessarily
result in a selective sweep. In these systems, adaptation can usually not be
modeled by a scalar fitness parameter for each strain Tikhonov and Monasson
(2018).

Strain dynamics in our model differ qualitatively from what is expected in the
traveling wave scenario. While adaptive mutations are highly overrepresented,
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 Huddleston et al. (2020); Barrat-Charlaix
et al. (2021). In particular, we observed in Barrat-Charlaix et al. (2021)
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.

###  .2 Results

####  .2.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 Gog and
Grenfell (2002); Yan et al. (2019). In the most general form, the model
describes $N$ variants of the virus labeled $a\in\\{1\ldots N\\}$ circulating
among $M$ groups of hosts with distinct exposure histories labeled
$i\in\\{1\ldots M\\}$ (immune groups). These groups could be different age
cohorts or be geographically separated. For each group $i$, we define
compartments $I_{i}^{a}$ and $S_{i}^{a}$ 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\leq
I^{a}_{i},S^{a}_{i}\leq 1$, and values of $I^{a}_{i}$ and $S^{a}_{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 S_{i}^{a}\sum_{j=1}^{M}C_{ij}I_{j}^{a}$.
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}^{a}$:

$\dot{I}_{i}^{a}=\alpha S_{i}^{a}\sum_{j=1}^{M}C_{ij}I_{j}^{a}-\delta
I_{i}^{a}.$ (1)

In the rest of this work, we will often consider the case $C_{ij}=1/M$,
corresponding to the case where hosts of different groups are equally likely
to encounter members of any other group. Note that in this case, if we define
the total number of susceptible or infected by $a$ as
$S^{a}=\sum_{i}S_{i}^{a}$ and $I^{a}=\sum_{i}I_{i}^{a}$, we can sum over the
immune groups in the above equation to obtain $\dot{I}^{a}=(\alpha
S^{a}-\delta)I^{a}$. We consequently 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 $S^{a}=\delta/\alpha$, such
that $R^{a}=1$.

When a host of group $i$ is infected by strain $b$, it also has a chance
$0\leq K_{i}^{ab}\leq 1$ to gain immunity to strain $a$. Thus, $S_{i}^{a}$
decreases at a rate proportional to $K_{i}^{ab}$ 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 $\mathbf{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}^{a}$ 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}),$
(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 (2002) assumes that immunity builds
up through exposure and not only through infection. This explains that the
change in $S_{i}^{a}$ is simply proportional to $S_{i}^{a}\cdot I_{j}^{b}$,
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 $K_{i}^{aa}=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 $K_{i}^{ab}$, but as highly similar by others, leading to
$K_{i}^{ab}\simeq 1$. Such a heterogeneous response by different immune
systems has been observed experimentally in the case of influenza: in Lee et
al. (2019); Welsh et al. (2023) 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 $K_{i}^{ab}$, while it had little effect for the serum
of other individuals, _i.e._ high $K_{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.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$. In this case, we can skip the indices $i,j\in\\{1\ldots M\\}$,
and we only have one cross-immunity matrix $\mathbf{K}$ that we parametrize as

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

with $0\leq b,f\leq 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:

$\begin{split}\dot{S}^{a}&=-\alpha
S^{a}\sum_{b\in\\{wt,v\\}}K^{ab}I^{b}+\gamma(1-S^{a}),\\\ \dot{I}^{a}&=\alpha
S^{a}I^{a}-\delta I^{a}.\\\ \end{split}$ (4)

We can immediately derive the equilibrium state for this simplified case. For
the susceptibles we find $S^{a}=\delta/\alpha$ for all strains $a$, while the
equilibrium prevalence is determined by the inverse of the cross-immunity
matrix $\mathbf{K}$:

$I^{a}_{eq}=\frac{\gamma}{\delta}\left(1-\frac{\delta}{\alpha}\right)\left(\mathbf{K}^{-1}\vec{1}\right)_{a},$
(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 $\sim 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 $\sim 0.15$, which is of similar magnitude as
empirical estimates for influenza (Kucharski et al., 2018).

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)}.$ (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\neq 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.

In our invasion scenario, only the wild-type is initially present in the
population, that is $I^{v}=0$ at $t<0$. The details of the equilibrium reached
by the system in the absence of the mutant is given in the SI. Importantly, in
this initial state, the number of hosts susceptible to the variant is always
larger than the equilibrium value $\delta/\alpha$ unless $f=1$, but smaller
than $1$ unless $f=0$. This comes from the fact that due to cross-immunity,
part of the hosts are immune to the variant even if it is not yet circulating.
As a result, when the variant is introduced at $t=0$, for instance because it
emerged by mutation, its growth rate is initially positive.

Figure 1: 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.99$. 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.

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},$
(7)

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. As a result, a
larger fraction of the hosts will initially be immune to the variant. If
mixing is rapid, the initial growth rate of the variant is smaller by a factor
$M$ compared to the one-group case. For $\varepsilon=0$, the equilibrium
frequency does not depend on $M$ and Equation 6 remains valid, while for small
$\varepsilon$ the equilibrium shifts slightly. If mixing is slow, the initial
growth of the variant is rapid as in the one-group case, but then spreads only
slowly across groups.

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$, and accordingly the initial growth rate of $I^{v}$
is smaller. 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. While this simple two-strain
model predicts that the two strains come to an equilibrium at frequency
$\beta$, 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. What allows the
initial growth is better described as 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 scalar fitness Tikhonov and Monasson (2018). An
important implication is that predicting the equilibrium frequency reached by
the variant and its ultimate fate is hard from the observation of a partial
frequency trajectory.

####  .2.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 $\beta_{2},\beta_{3},\ldots$ 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.

Figure 2: 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 8, all
starting at $x\simeq 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
scalar fitness parametrization of Equation 10. 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.

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 $\rho$, meaning
that they reach equilibrium before the next variant emerges, and if new
variants have a similar cross-immunity with all existing variants (see SI),
the dynamics of $x(t)$ are described by a particular random walk: in each time
interval $\text{d}t$, a partial sweep of amplitude $\beta$ occurs with
probability $\rho\text{d}t\cdot P_{\beta}(\beta)$, changing $x$ in the
following way:

$x(t+\text{d}t)=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}$ (8)

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
$\text{d}t$, that is with probability $1-\rho\text{d}t$, $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 (Gillespie, 2000).

Examples of trajectories from the random walk are shown on panel C of Figure
2, all initially starting around $x_{0}\simeq 1/2$. Two trajectories are seen
converging monotonously to $0$ and $1$. This is a consequence of one
interesting property of Equation 8: 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 (Hallatschek, 2018). 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{split}\langle\Delta x\rangle&=0\\\ \langle\Delta
x^{2}\rangle&=\rho\langle\beta^{2}\rangle_{P_{\beta}}x(1-x).\end{split}$ (9)

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 $x_{0}$ to reach $1$ and fix, and a
chance $1-x_{0}$ to reach $0$ and vanish.

On the other hand, the second moment closely reminds of the theory of neutral
drift Kimura (1964): 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_{e}^{-1}=\rho\langle\beta^{2}\rangle_{P_{\beta}}$ completely independent of
the size of the viral population. However, there are important difference 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 $\langle\beta^{k}\rangle_{P_{\beta}}$.
Depending on higher moments of $P_{\beta}$, allele dynamics will deviate
qualitatively from neutral behavior.

####  .2.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\times 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. For this reason, we propose a simpler
empirical model for partial sweeps where 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\mathrm{and}\quad\dot{s}=-\nu xs.$ (10)

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. 10).

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}.$ (11)

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. Due to its simplicity, this model does not show the
oscillatory behavior of the SIR model, but otherwise recapitulates the salient
feature of invading immune evasive variants: (i) initial exponential growth,
and (ii) eventual saturation at an intermediate frequency. This simplified
model of partial sweeps leads to the same random walk dynamics as the SIR
model above and allows us to study the implications of this dynamics
systematically.

####  .2.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 (Barrat-
Charlaix et al., 2021).

In Figure 3, we reproduce some of our results of Barrat-Charlaix et al. (2021)
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,\ldots,t$, what can we say about variant frequencies at times
$t+1,\ldots$? 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^{\star}$. Adaptive beneficial mutations
should in principle be overrepresented in this group. Looking at the
corresponding frequency trajectories, we can see on average whether they keep
on growing beyond $x^{\star}$. 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^{\star}=0.4$, with their average
displayed in black. The panels on the right show the average trajectory for
different threshold values $x^{\star}$ 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^{\star}$ instead of growing. This surprising
result is in line with the study in Barrat-Charlaix et al. (2021) 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^{\star}$.
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.

Figure 3: Retrospective analysis of predictability of viral evolution:
frequency trajectories of all amino acid substitutions that are observed to
rise from frequency $0$ to $x^{\star}$ for Top: influenza virus A/H3N2 from
2000 to 2023, and Bottom: SARS-CoV-2 from 2020-2023. Left: all trajectories
for $x^{\star}=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^{\star}$ (grey lines).

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. 8, 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. 10, 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 so as to favor mutations. In practice, we use an
exponential distribution for $P_{s}$, meaning that the typical magnitude of
initial fitness effects is described by only one 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_{\beta}$ of
partial sweeps size depending on $\nu/s_{0}$, and (ii) the ratio of the
variant emergence rate and their growth rate $\rho/s_{0}$, which determines
how often sweeps overlap and interact.

We use this simulation to address the question of predictability: given the
state of the population at generations $0,1,\ldots,t$, can we predict its
state at future times $t+1,\ldots$? Specifically, we ask whether we can
predict the frequency $x(t+\Delta t)$ of a variant $A$, given it is at
frequency $x$ at time $t$, as we did previously for influenza virus Barrat-
Charlaix et al. (2021), see Fig. 3. The dynamics of isolated selective sweeps
($\rho/s_{0}\ll 1$, $\nu/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
$\rho/s_{0}$ due to clonal interference (Strelkowa and Lässig, 2012), for
example when an adaptive variant is outcompeted by an even more adaptive one.
We also expect predictability to decrease with increasing $\nu/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:
$\nu\in[0,s_{0}/3,s_{0},3s_{0}]$ and low clonal interference
$\rho/s_{0}=0.05$. The case $\nu=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 $\nu/s_{0}$, corresponding to a quicker decay of fitness, predictability
gradually declines and becomes negligible for $\nu/s\gg 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.

Figure 4:  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 $\nu/s_{0}$), the
trajectories crossing $x^{*}$ have little inertia, while stable growth
advantage (small $\nu/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 $\nu/s_{0}$ both result in
gradual loss of predictability. 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\langle\beta^{2}\rangle$. Points
correspond to different choices of parameters $s_{0}$

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 $\rho/s_{0}$ and four values of $\nu/s_{0}$. Clonal
interference increases when going from left to right among these panels
(increasing $\rho/s_{0}$), while the intensity of fitness decay increases when
going from blue to red curves (increasing $\nu$). Increasing either
$\rho/s_{0}$ or $\nu/s_{0}$ reduces $p_{fix}$ towards the dashed diagonal
corresponding to $p_{fix}=x^{*}$. However, as observed previously (Barrat-
Charlaix et al., 2021), in the classic scenario with stable fitness effects
$\nu/s_{0}=0$ considerable predictability remains even in cases of strong
interference (blue curve in panel B).

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/\rho\langle\beta^{2}\rangle_{P_{\beta}}$, 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 $\Lambda$-coalescent Schweinsberg
(2000); Berestycki (2009) (see 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_{\beta}$ and
a given $\rho$ setting $T_{2}$. 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 $\beta$. 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 $\beta$ which slows the dynamics.

###  .3 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 (Bush et al.,
1999; Luksza and Lässig, 2014; Neher et al., 2014; Huddleston et al., 2020).
The A/H3N2 subtype in particular undergoes rapid antigenic change through
frequent substitutions in prominent epitopes on its surface proteins (Smith et
al., 2004; Bhatt et al., 2011; Neher et al., 2016; Kistler and Bedford, 2023).
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 (Barrat-Charlaix
et al., 2021; Huddleston et al., 2020). 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 Meijers et al.
(2023).

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 _(iii)_ waning or turn-over 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
(Bhatt et al., 2013). On the other hand, most substitutions do not sweep to
fixation but tend to meander in a quasi-neutral fashion (Barrat-Charlaix et
al., 2021). 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 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” (Fonville et
al., 2014)), 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/H1N1pdm, which emerged in humans in 2009,
exhibited more consistent trajectory dynamics than A/H3N2 (Barrat-Charlaix et
al., 2021).

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 (Pelletier et al., 2009). 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.jl

  * •

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

  * •

Code to generate empirical frequency trajectories and their averages is
available as scripts ‘flu_fixation.py’ and ‘sc2_fixation.py’ in
https://github.com/nextstrain/flu_frequencies on branch ‘fixation’.

## Data availability

Sequence data of influenza viruses was obtained from GISAID (Shu and McCauley,
2017). 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.

## References

  * Roemer et al. [2023] Cornelius Roemer, Daniel J. Sheward, Ryan Hisner, Federico Gueli, Hitoshi Sakaguchi, Nicholas Frohberg, Josette Schoenmakers, Kenta Sato, Áine O’Toole, Andrew Rambaut, Oliver G. Pybus, Christopher Ruis, Ben Murrell, and Thomas P. Peacock.  SARS-CoV-2 evolution in the Omicron era.  _Nature Microbiology_ , October 2023.  ISSN 2058-5276.  doi: 10.1038/s41564-023-01504-w.
  * Petrova and Russell [2018] Velislava N. Petrova and Colin A. Russell.  The evolution of seasonal influenza viruses.  _Nature Reviews Microbiology_ , 16(1):47–60, January 2018.  ISSN 1740-1534.  doi: 10.1038/nrmicro.2017.118.  URL https://www.nature.com/articles/nrmicro.2017.118.  Number: 1 Publisher: Nature Publishing Group.
  * Morris et al. [2018] Dylan H. Morris, Katelyn M. Gostic, Simone Pompei, Trevor Bedford, Marta Łuksza, Richard A. Neher, Bryan T. Grenfell, Michael Lässig, and John W. McCauley.  Predictive Modeling of Influenza Shows the Promise of Applied Evolutionary Biology.  _Trends in Microbiology_ , 26(2):102–118, February 2018.  ISSN 0966-842X.  doi: 10.1016/j.tim.2017.09.004.  URL https://www.sciencedirect.com/science/article/pii/S0966842X17302093.
  * Meijers et al. [2023] Matthijs Meijers, Denis Ruchnewitz, Jan Eberhardt, Marta Łuksza, and Michael Lässig.  Population immunity predicts evolutionary trajectories of SARS-CoV-2.  _Cell_ , 186(23):5151–5164.e13, November 2023\.  ISSN 0092-8674, 1097-4172.  doi: 10.1016/j.cell.2023.09.022.  URL https://www.cell.com/cell/abstract/S0092-8674(23)01076-0.  Publisher: Elsevier.
  * Tsimring et al. [1996] Lev S. Tsimring, Herbert Levine, and David A. Kessler.  RNA Virus Evolution via a Fitness-Space Model.  _Physical Review Letters_ , 76(23):4440–4443, June 1996.  ISSN 0031-9007, 1079-7114.  doi: 10.1103/PhysRevLett.76.4440.  URL https://link.aps.org/doi/10.1103/PhysRevLett.76.4440.
  * Rouzine et al. [2003] Igor M. Rouzine, John Wakeley, and John M. Coffin.  The solitary wave of asexual evolution.  _Proceedings of the National Academy of Sciences_ , 100(2):587–592, January 2003.  doi: 10.1073/pnas.242719299.  URL https://www.pnas.org/doi/abs/10.1073/pnas.242719299.  Publisher: Proceedings of the National Academy of Sciences.
  * Desai and Fisher [2007] Michael M Desai and Daniel S Fisher.  Beneficial Mutation–Selection Balance and the Effect of Linkage on Positive Selection.  _Genetics_ , 176(3):1759–1798, July 2007.  ISSN 1943-2631.  doi: 10.1534/genetics.106.067678.  URL https://doi.org/10.1534/genetics.106.067678.
  * Neher [2013] Richard A. Neher.  Genetic Draft, Selective Interference, and Population Genetics of Rapid Adaptation.  _Annual Review of Ecology, Evolution, and Systematics_ , 44(1):195–215, 2013.  doi: 10.1146/annurev-ecolsys-110512-135920.  URL https://doi.org/10.1146/annurev-ecolsys-110512-135920.  _eprint: https://doi.org/10.1146/annurev-ecolsys-110512-135920.
  * Luksza and Lässig [2014] Marta Luksza and Michael Lässig.  A predictive fitness model for influenza.  _Nature_ , 507(7490):57–61, March 2014.  ISSN 1476-4687.  doi: 10.1038/nature13087.
  * Neher et al. [2014] Richard A Neher, Colin A Russell, and Boris I Shraiman.  Predicting evolution from the shape of genealogical trees.  _eLife_ , 3, November 2014.  ISSN 2050-084X.  doi: 10.7554/eLife.03568.  URL https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4227306/.
  * Huddleston et al. [2020] John Huddleston, John R Barnes, Thomas Rowe, Xiyan Xu, Rebecca Kondor, David E Wentworth, Lynne Whittaker, Burcu Ermetal, Rodney Stuart Daniels, John W McCauley, Seiichiro Fujisaki, Kazuya Nakamura, Noriko Kishida, Shinji Watanabe, Hideki Hasegawa, Ian Barr, Kanta Subbarao, Pierre Barrat-Charlaix, Richard A Neher, and Trevor Bedford.  Integrating genotypes and phenotypes improves long-term forecasts of seasonal influenza A/H3N2 evolution.  _eLife_ , 9:e60067, September 2020.  ISSN 2050-084X.  doi: 10.7554/eLife.60067.  URL https://doi.org/10.7554/eLife.60067.  Publisher: eLife Sciences Publications, Ltd.
  * Barrat-Charlaix et al. [2021] Pierre Barrat-Charlaix, John Huddleston, Trevor Bedford, and Richard A Neher.  Limited Predictability of Amino Acid Substitutions in Seasonal Influenza Viruses.  _Molecular Biology and Evolution_ , (msab065), March 2021\.  ISSN 0737-4038.  doi: 10.1093/molbev/msab065.  URL https://doi.org/10.1093/molbev/msab065.
  * Gupta et al. [1998] Sunetra Gupta, Neil Ferguson, and Roy Anderson.  Chaos, Persistence, and Evolution of Strain Structure in Antigenically Diverse Infectious Agents.  _Science_ , 280(5365):912–915, May 1998.  ISSN 0036-8075, 1095-9203.  doi: 10.1126/science.280.5365.912.  URL https://www.science.org/doi/10.1126/science.280.5365.912.
  * Gog and Grenfell [2002] Julia R. Gog and Bryan T. Grenfell.  Dynamics and selection of many-strain pathogens.  _Proceedings of the National Academy of Sciences_ , 99(26):17209–17214, December 2002.  doi: 10.1073/pnas.252512799.  URL https://www.pnas.org/doi/abs/10.1073/pnas.252512799.  Publisher: Proceedings of the National Academy of Sciences.
  * Yan et al. [2019] Le Yan, Richard A Neher, and Boris I Shraiman.  Phylodynamic theory of persistence, extinction and speciation of rapidly adapting pathogens.  _eLife_ , 8:e44205, September 2019.  ISSN 2050-084X.  doi: 10.7554/eLife.44205.  URL https://doi.org/10.7554/eLife.44205.  Publisher: eLife Sciences Publications, Ltd.
  * Marchi et al. [2021] Jacopo Marchi, Michael Lässig, Aleksandra M. Walczak, and Thierry Mora.  Antigenic waves of virus–immune coevolution.  _Proceedings of the National Academy of Sciences_ , 118(27):e2103398118, July 2021.  doi: 10.1073/pnas.2103398118.  URL https://www.pnas.org/doi/full/10.1073/pnas.2103398118.  Publisher: Proceedings of the National Academy of Sciences.
  * Chardès et al. [2023] Victor Chardès, Andrea Mazzolini, Thierry Mora, and Aleksandra M. Walczak.  Evolutionary stability of antigenically escaping viruses.  _Proceedings of the National Academy of Sciences_ , 120(44):e2307712120, October 2023.  doi: 10.1073/pnas.2307712120.  URL https://www.pnas.org/doi/10.1073/pnas.2307712120.  Publisher: Proceedings of the National Academy of Sciences.
  * Rouzine and Rozhnova [2018] Igor M. Rouzine and Ganna Rozhnova.  Antigenic evolution of viruses in host populations.  _PLOS Pathogens_ , 14(9):e1007291, September 2018\.  ISSN 1553-7374.  doi: 10.1371/journal.ppat.1007291.  URL https://journals.plos.org/plospathogens/article?id=10.1371/journal.ppat.1007291.  Publisher: Public Library of Science.
  * Lee et al. [2019] Juhye M Lee, Rachel Eguia, Seth J Zost, Saket Choudhary, Patrick C Wilson, Trevor Bedford, Terry Stevens-Ayers, Michael Boeckh, Aeron C Hurt, Seema S Lakdawala, Scott E Hensley, and Jesse D Bloom.  Mapping person-to-person variation in viral mutations that escape polyclonal serum targeting influenza hemagglutinin.  _eLife_ , 8:e49324, August 2019.  ISSN 2050-084X.  doi: 10.7554/eLife.49324.  URL https://doi.org/10.7554/eLife.49324.  Publisher: eLife Sciences Publications, Ltd.
  * Welsh et al. [2023] Frances C. Welsh, Rachel T. Eguia, Juhye M. Lee, Hugh K. Haddox, Jared Galloway, Nguyen Van Vinh, Andrea N. Loes, John Huddleston, Timothy C. Yu, Mai Quynh Le, Nguyen TD Nhat, Nguyen Thi Le Thanh, Alexander L. Greninger, Helen Y. Chu, Janet A. Englund, Trevor Bedford, Frederick A. Matsen, Maciej F. Boni, and Jesse D. Bloom.  Age-dependent heterogeneity in the antigenic effects of mutations to influenza hemagglutinin, December 2023.  URL https://www.biorxiv.org/content/10.1101/2023.12.12.571235v1.  Pages: 2023.12.12.571235 Section: New Results.
  * Tikhonov and Monasson [2018] Mikhail Tikhonov and Remi Monasson.  Innovation Rather than Improvement: A Solvable High-Dimensional Model Highlights the Limitations of Scalar Fitness.  _Journal of Statistical Physics_ , 172(1):74–104, July 2018.  ISSN 1572-9613.  doi: 10.1007/s10955-018-1956-6.  URL https://doi.org/10.1007/s10955-018-1956-6.
  * Kucharski et al. [2018] Adam J. Kucharski, Justin Lessler, Derek A. T. Cummings, and Steven Riley.  Timescales of influenza A/H3N2 antibody dynamics.  _PLoS biology_ , 16(8):e2004974, August 2018.  ISSN 1545-7885.  doi: 10.1371/journal.pbio.2004974.
  * Gillespie [2000] John H Gillespie.  Genetic Drift in an Infinite Population: The Pseudohitchhiking Model.  _Genetics_ , 155(2):909–919, June 2000.  ISSN 1943-2631.  doi: 10.1093/genetics/155.2.909.  URL https://doi.org/10.1093/genetics/155.2.909.
  * Hallatschek [2018] Oskar Hallatschek.  Selection-Like Biases Emerge in Population Models with Recurrent Jackpot Events.  _Genetics_ , 210(3):1053–1073, November 2018\.  ISSN 1943-2631.  doi: 10.1534/genetics.118.301516.  URL https://doi.org/10.1534/genetics.118.301516.
  * Kimura [1964] Motoo Kimura.  Diffusion models in population genetics.  _Journal of Applied Probability_ , 1(2):177–232, December 1964.  ISSN 0021-9002, 1475-6072.  doi: 10.2307/3211856.  URL https://www.cambridge.org/core/product/identifier/S0021900200108368/type/journal_article.
  * Strelkowa and Lässig [2012] Natalja Strelkowa and Michael Lässig.  Clonal Interference in the Evolution of Influenza.  _Genetics_ , 192(2):671–682, October 2012.  ISSN 1943-2631.  doi: 10.1534/genetics.112.143396.  URL https://doi.org/10.1534/genetics.112.143396.
  * Schweinsberg [2000] Jason Schweinsberg.  Coalescents with Simultaneous Multiple Collisions.  _Electronic Journal of Probability_ , 5(none):1–50, January 2000.  ISSN 1083-6489, 1083-6489.  doi: 10.1214/EJP.v5-68.  URL https://projecteuclid.org/journals/electronic-journal-of-probability/volume-5/issue-none/Coalescents-with-Simultaneous-Multiple-Collisions/10.1214/EJP.v5-68.full.  Publisher: Institute of Mathematical Statistics and Bernoulli Society.
  * Berestycki [2009] Nathanaël Berestycki.  Recent progress in coalescent theory.  _Ensaios Matemáticos_ , 16(1), 2009.  ISSN 21750432.  doi: 10.21711/217504322009/em161.  URL http://ensaios.sbm.org.br/wp-content/uploads/sites/14/2019/04/EM_16.pdf.
  * Bush et al. [1999] Robin M. Bush, Catherine A. Bender, Kanta Subbarao, Nancy J. Cox, and Walter M. Fitch.  Predicting the Evolution of Human Influenza A.  _Science_ , 286(5446):1921–1925, December 1999\.  doi: 10.1126/science.286.5446.1921.  URL https://www.science.org/doi/10.1126/science.286.5446.1921.  Publisher: American Association for the Advancement of Science.
  * Smith et al. [2004] Derek J. Smith, Alan S. Lapedes, Jan C. De Jong, Theo M. Bestebroer, Guus F. Rimmelzwaan, Albert D. M. E. Osterhaus, and Ron A. M. Fouchier.  Mapping the Antigenic and Genetic Evolution of Influenza Virus.  _Science_ , 305(5682):371–376, July 2004.  ISSN 0036-8075, 1095-9203.  doi: 10.1126/science.1097211.  URL https://www.science.org/doi/10.1126/science.1097211.
  * Bhatt et al. [2011] Samir Bhatt, Edward C. Holmes, and Oliver G. Pybus.  The Genomic Rate of Molecular Adaptation of the Human Influenza A Virus.  _Molecular Biology and Evolution_ , 28(9):2443–2451, September 2011.  ISSN 0737-4038.  doi: 10.1093/molbev/msr044.  URL https://doi.org/10.1093/molbev/msr044.
  * Neher et al. [2016] Richard A. Neher, Trevor Bedford, Rodney S. Daniels, Colin A. Russell, and Boris I. Shraiman.  Prediction, dynamics, and visualization of antigenic phenotypes of seasonal influenza viruses.  _Proceedings of the National Academy of Sciences_ , 113(12):E1701–E1709, March 2016.  doi: 10.1073/pnas.1525578113.  URL https://www.pnas.org/doi/full/10.1073/pnas.1525578113.  Publisher: Proceedings of the National Academy of Sciences.
  * Kistler and Bedford [2023] Kathryn E. Kistler and Trevor Bedford.  An atlas of continuous adaptive evolution in endemic human viruses.  _Cell Host & Microbe_, October 2023.  ISSN 1931-3128.  doi: 10.1016/j.chom.2023.09.012.  URL https://www.sciencedirect.com/science/article/pii/S1931312823003803.
  * Bhatt et al. [2013] S. Bhatt, T. T. Lam, S. J. Lycett, A. J. Leigh Brown, T. A. Bowden, E. C. Holmes, Y. Guan, J. L. N. Wood, I. H. Brown, P. Kellam, Combating Swine Influenza Consortium, and O. G. Pybus.  The evolutionary dynamics of influenza A virus adaptation to mammalian hosts.  _Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences_ , 368(1614):20120382, March 2013\.  ISSN 1471-2970.  doi: 10.1098/rstb.2012.0382.
  * Fonville et al. [2014] J. M. Fonville, S. H. Wilks, S. L. James, A. Fox, M. Ventresca, M. Aban, L. Xue, T. C. Jones, Le N. M. H., Pham Q. T., Tran N. D., Y. Wong, A. Mosterin, L. C. Katzelnick, D. Labonte, Le T. T., G. van der Net, E. Skepner, C. A. Russell, T. D. Kaplan, G. F. Rimmelzwaan, N. Masurel, J. C. de Jong, A. Palache, W. E. P. Beyer, Le Q. M., Nguyen T. H., H. F. L. Wertheim, A. C. Hurt, A. D. M. E. Osterhaus, I. G. Barr, R. A. M. Fouchier, P. W. Horby, and D. J. Smith.  Antibody landscapes after influenza virus infection or vaccination.  _Science_ , 346(6212):996–1000, November 2014\.  doi: 10.1126/science.1256427.  URL https://www.science.org/doi/10.1126/science.1256427.  Publisher: American Association for the Advancement of Science.
  * Pelletier et al. [2009] F. Pelletier, D. Garant, and A.p. Hendry.  Eco-evolutionary dynamics.  _Philosophical Transactions of the Royal Society B: Biological Sciences_ , 364(1523):1483–1489, June 2009.  doi: 10.1098/rstb.2009.0027.  URL https://royalsocietypublishing.org/doi/full/10.1098/rstb.2009.0027.  Publisher: Royal Society.
  * Shu and McCauley [2017] Yuelong Shu and John McCauley.  GISAID: Global initiative on sharing all influenza data – from vision to reality.  _Eurosurveillance_ , 22(13):30494, March 2017\.  ISSN 1560-7917.  doi: 10.2807/1560-7917.ES.2017.22.13.30494.  URL https://www.eurosurveillance.org/content/10.2807/1560-7917.ES.2017.22.13.30494?crawler=true.  Publisher: European Centre for Disease Prevention and Control.
  * Bolthausen and Sznitman [1998] E. Bolthausen and A.-S. Sznitman.  On Ruelle’s Probability Cascades and an Abstract Cavity Method.  _Communications in Mathematical Physics_ , 197(2):247–276, October 1998.  ISSN 1432-0916.  doi: 10.1007/s002200050450.  URL https://doi.org/10.1007/s002200050450.
  * Brunet et al. [2007] É. Brunet, B. Derrida, A. H. Mueller, and S. Munier.  Effect of selection on ancestry: An exactly soluble case and its phenomenological generalization.  _Physical Review E_ , 76(4):041104, October 2007\.  doi: 10.1103/PhysRevE.76.041104.  URL https://link.aps.org/doi/10.1103/PhysRevE.76.041104.  Publisher: American Physical Society.
  * Neher and Hallatschek [2013] R. A. Neher and O. Hallatschek.  Genealogies of rapidly adapting populations.  _Proceedings of the National Academy of Sciences_ , 110(2):437–442, January 2013.  ISSN 0027-8424, 1091-6490.  doi: 10.1073/pnas.1213113110.  URL http://www.pnas.org/cgi/doi/10.1073/pnas.1213113110.

##  Appendix A SIR model

###  A.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},\ldots,S^{N}]$ and
$\vec{I}=[I^{1},\ldots,I^{N}]$ will respectively represent the compartments of
hosts susceptible and infectious to each strain; the matrix $\mathbf{K}$
describes the cross-immunity; the vector $\vec{1}$ 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{1}$. Introducing this in the equation for the
derivative of $S^{a}$, we obtain the following equilibrium:

$\displaystyle\vec{S}_{eq}$ $\displaystyle=\frac{\delta}{\alpha}\vec{1},$ (12)
$\displaystyle\vec{I}_{eq}$
$\displaystyle=\frac{\gamma}{\delta}\left(1-\frac{\delta}{\alpha}\right)\mathbf{K}^{-1}\vec{1}.$

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\leq K^{ab}<1$
for $a\neq b$.

###  A.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\times 2$ matrix

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

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

$\displaystyle S^{wt}=S^{m}=\frac{\delta}{\alpha}$ (14)
$\displaystyle\vec{I}=\frac{\gamma}{\delta}\left(1-\frac{\delta}{\alpha}\right)\mathbf{K}\vec{1},$

where $\vec{I}$ stands for $[I^{wt},I^{m}]$ and $\vec{1}$ for $[1,1]$. It is
straightforward to invert the cross-immunity matrix, and we obtain

$\displaystyle I^{wt}$
$\displaystyle=\frac{\gamma(1-\delta\alpha^{-1})}{\delta(1-bf)}(1-b)$ (15)
$\displaystyle I^{m}$
$\displaystyle=\frac{\gamma(1-\delta\alpha^{-1})}{\delta(1-bf)}(1-f).$

Note that without cross-immunity, the number of infected by either virus would
be $\frac{\gamma(1-\delta\alpha^{-1})}{\delta}$. Positive values of $b$ and
$f$ thus have the effect of lowering the equilibrium values of $I^{m}$ and
$I^{wt}$ 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}.$ (16)

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^{wt}$ to its threshold value $R_{0}^{-1}$, 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$.

###  A.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:

$\displaystyle S^{wt}$ $\displaystyle=\delta/\alpha=R_{0}^{-1}$ (17)
$\displaystyle I^{wt}$ $\displaystyle=\frac{\gamma}{\delta}(1-S^{wt}).$

We then set the derivative of $S^{m}$ to $0$:

$\displaystyle-\alpha fS^{m}I^{wt}+\gamma(1-S^{m})=0$ (18)
$\displaystyle\rightarrow$ $\displaystyle S^{m}=\frac{1}{1+f(R_{0}-1)}$
$\displaystyle=\frac{\delta}{\delta+f(\alpha-\delta)}>\frac{\delta}{\alpha}=S^{wt}.$

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).$
(19)

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
$\mathbf{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)}.$ (20)

In a given immune group $i$, the mutant growth rate is proportional to
$S_{i}^{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}\simeq 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}I^{m}_{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}I^{m}_{i}$:

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

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}.$
(22)

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 21. 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$.

###  A.4 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
immunty 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:

$\displaystyle I^{a}$ $\displaystyle=\alpha(1-R^{a}-R^{ab})I^{a}-\delta
I^{a},$ (23) $\displaystyle I^{b}$
$\displaystyle=\alpha(1-R^{b}-R^{ab})I^{b}-\delta I^{b},$

and three for the immune:

$\displaystyle R^{a}$ $\displaystyle=\alpha(1-K^{ba})R^{0}I^{a}-\gamma R^{a},$
(24) $\displaystyle R^{b}$ $\displaystyle=\alpha(1-K^{ab})R^{0}I^{b}-\gamma
R^{b},$ $\displaystyle R^{ab}$ $\displaystyle=\alpha
R^{0}(K^{ba}I^{a}+K^{ab}I^{b})+\alpha(R^{a}I^{b}+R^{b}I^{a})-\gamma R^{ab}.$

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.

Figure S 1: Left: Dynamics of the frequency of the variant for the SIR model
from Equations 23&24 using the invasion scenario from the main text. Two
$2\times 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 $\beta$ for this new SIR model ($y$-axis) versus the
$\beta$ from the main text ($x$-axis). Each point corresponds to a given pair
$(f,b)$.

###  A.5 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{1}$ and
$\vec{f}=f\cdot\vec{1}$, 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\ldots 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},$
(25)

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}}.$ (26)

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 $\mathbf{\tilde{K}}$ is now written as

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

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 $\tilde{I}$:

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

The question we ask here is whether the relative frequency of two variants
$1\leq a,b\leq 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{{\scriptstyle?}}{{=}}\frac{\tilde{I}_{a}}{\tilde{I}_{a}+\tilde{I}_{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}\;\;\;\text{and}\;\;\;\vec{f}=f\cdot\vec{1}_{N},$
(28)

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
$\mathbf{\tilde{K}}^{-1}\vec{1}_{N+1}.$ 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
CA^{-1}&-A^{-1}B\Lambda\\\ -\Lambda CA^{-1}&\Lambda\end{bmatrix},$

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

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

We immediatly see that
$\Lambda=\left(1-\vec{f}^{\intercal}\mathbf{K}^{-1}\vec{b}\right)^{-1}$
reduces to a scalar that we note $\lambda$ for more clarity. We also define
the other scalar value
$\mu=\vec{1}^{\intercal}_{N}\mathbf{K}^{-1}\vec{1}_{N}$. A few manipulations
give us the following for $\mathbf{\tilde{K}}^{-1}$:

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

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

$\displaystyle\tilde{I}$ $\displaystyle=\mathbf{\tilde{K}}^{-1}\vec{1}_{N+1}$
$\displaystyle=\left[\underbrace{\left(\mathbf{K}^{-1}+\lambda\mathbf{K}^{-1}\vec{b}\cdot\vec{f}^{T}\mathbf{K}^{-1}\right)\vec{1}_{N}-\lambda\mathbf{K}^{-1}\vec{b}}_{\text{size}\;N}\,;\underbrace{-\lambda\vec{f}^{\intercal}\mathbf{K}^{-1}\vec{1}_{N}+\lambda}_{\text{scalar}}\right]$
$\displaystyle=\left[I+bf\mu\lambda I-b\lambda I\,;\lambda(1-\mu f)\right]$
$\displaystyle=\left[(1-b\lambda+bf\lambda\mu)I;\lambda(1-\mu f)\right],$

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-b\lambda(1-f\mu)$. This implies that the relative frequencies of the
original variants are conserved when adding a new one.

###  A.6 Case with intrinsic fitness effects

In the SI model of the main text, we assume that the transmission rate
$\alpha$ 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 $\mathbf{K}$ with off-diagonal terms
$f$ and $b$.

The equations of motion are now

$\displaystyle\dot{S}^{a}$ $\displaystyle=-\alpha
S^{a}\sum_{b\in\\{wt,m\\}}\phi^{b}K^{ab}I^{b}+\gamma(1-S^{a})$ (29)
$\displaystyle\dot{I}^{a}$ $\displaystyle=\alpha\phi^{a}S^{a}I^{a}-\delta
I^{a}.$

Computing the equilibrium, we immediately obtain

$\displaystyle S^{a}$ $\displaystyle=\frac{\delta}{\alpha\phi^{a}}$ (30)
$\displaystyle I$
$\displaystyle=\frac{\gamma}{\delta}\cdot\mathbf{G}^{-1}\vec{h},$

where we have defined the following quantities:

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

The quantities $h^{a}$ can be interpreted as a scaled growth rate of each
variant given a fully susceptible population, and the matrix $\mathbf{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^{ab}=\frac{\phi^{b}}{\phi^{a}}K^{ab}$.

Inverting $\mathbf{G}$ and simplifying the equations a bit, we obtain

$\displaystyle I^{wt}$
$\displaystyle=\frac{\gamma}{\delta}\frac{\alpha\phi^{wt}-\delta}{\alpha\phi^{wt}}\cdot\frac{1-b\xi}{1-bf},$
(32) $\displaystyle I^{m}$
$\displaystyle=\frac{\gamma}{\delta}\frac{\alpha\phi^{m}-\delta}{\alpha\phi^{m}}\cdot\frac{1-f\xi^{-1}}{1-bf},$

where we defined

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

Note the interesting structure of equations 32: for each variant, they involve
a first term $\frac{\alpha\phi^{a}-\delta}{\alpha\phi^{a}}$ that depends only
of the intrinsic growth rate of the mutant, and a second $\frac{1-b\xi}{1-bf}$
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,$ (34)

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.

###  A.7 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

$\displaystyle Q$
$\displaystyle=\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}$ $\displaystyle g_{1}$
$\displaystyle=\frac{\gamma}{\delta}(\alpha-\delta)\frac{1-b}{1-bf}$
$\displaystyle g_{2}$
$\displaystyle=\frac{\gamma}{\delta}(\alpha-\delta)\frac{1-f}{1-bf}.$

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

$\displaystyle\lambda_{1}$
$\displaystyle=-\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}$
$\displaystyle\lambda_{2}$
$\displaystyle=-\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 $\gamma$to 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:

$\displaystyle\omega_{high}$
$\displaystyle=\frac{\sqrt{\gamma(\alpha-\delta)}}{2\pi},$
$\displaystyle\omega_{low}$
$\displaystyle=\omega_{high}\sqrt{\frac{(1-b)(1-f)}{1-bf}}.$

Note that since $(1-b)(1-f)/(1-bf)\leq 1$, we always have
$\omega_{low}\leq\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}\simeq 0.016$. If one generation is one week, this
gives us a period $\omega_{high}^{-1}\simeq 60w$, that is approximately one
year.

##  Appendix B Other

###  B.1 Differential equation with expiring fitness effects

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),\;\;\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{\text{d}x}{\text{d}s}=\nu^{-1}(x-1).$

This immediately gives us

$x(s)=1+\lambda e^{s/\nu}$

with $\lambda$ constant. At $t=0$, we have $s=s_{0}$ and $x=x_{0}\ll 1$, while
for $t\rightarrow\infty$ we have $s=0$ and $x=\beta$ to be determined. From
the $t=0$ case we obtain $\lambda=(x_{0}-1)e^{-s_{0}/\nu}$, and from
$t\rightarrow\infty$ we get $\beta=1+\lambda$. 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)\simeq\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}$. For
instance, taking $r=3/4$, this 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 $\rho$, this has probability
$1-\exp\left(-\rho(T_{r}(s_{1})-T_{1-r}(s_{2}))\right)$. Since the two sweeps
have random initial fitness effects, we find that the overall probability for
two consecutive sweeps to overlap is

$P(\text{overlap})=\int_{0}^{\infty}\text{d}s_{1}\text{d}s_{2}P_{s}(s_{1})P_{s}(s_{2})\left\\{1-e^{-\rho\left(T_{r}(s_{1})-T_{1-r}(s_{2})\right)}\right\\}.$

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.

###  B.2 Distribution of partial sweep size $\beta$

This section discusses the distribution of the size of partial sweeps $\beta$
in the context of Equation 10 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:

$\displaystyle P(\beta<x)$ $\displaystyle=P(s<-\alpha\log(1-x))$
$\displaystyle=s_{0}^{-1}\int_{0}^{-\alpha\log(1-x)}e^{-s/s_{0}}\text{d}s$
$\displaystyle=1-(1-x)^{\alpha/s_{0}}.$

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

$P(\beta)\propto(1-\beta)^{\alpha/s_{0}-1}.$

This distribution can accomodate various shape: for $\alpha/s_{0}>1$ it peaks
at $0$, and for $\alpha/s_{0}<1$ it peaks at $1$. We can also compute the
following formular for the moments of $\beta$:

$\langle\beta\rangle=\frac{s_{0}}{s_{0}+\alpha}\;\text{and}\;\langle\beta^{2}\rangle=\frac{s_{0}}{s_{0}+\alpha}\frac{2s_{0}}{2s_{0}+\alpha}.$

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
accomodate 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 11 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]$.

###  B.3 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}\left(1-(1-\beta)^{t}x_{0}\right).$ (35)

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

$P_{down}\simeq(1-x_{0})e^{x_{0}(1-\beta^{-1})}.$

Since the random walk is invariant by the change $x\rightarrow 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$.

Figure S 2: 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 35 up to $t=100$.

The quality of the approximation is quite good, as can be seen from Figure S2.
The same figure also shows that this probability is relatively high, even for
$x_{0}$ far from either boundary.

###  B.4 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 $\lambda_{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}$.

$\displaystyle\lambda_{n}(k)$
$\displaystyle=\rho\langle\beta^{k}(1-\beta)^{n-k}\rangle,$ (36)
$\displaystyle=\int_{0}^{1}\beta^{k}(1-\beta)^{n-k}P(\beta)\text{d}\beta,$
$\displaystyle=\int_{0}^{1}\beta^{k-2}(1-\beta)^{n-k}\frac{\beta^{2}P(\beta)}{\langle\beta^{2}\rangle}\text{d}\beta.$

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

  * •

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

  * •

if $\Lambda(\beta)$ is uniform in $[0,1]$, meaning
$P(\beta)\propto\beta^{-2}$, we obtain the Bolthausen-Sznitman coalescent
Bolthausen and Sznitman [1998], Berestycki [2009], which is used to describe
the genealogy of populations under strong selection Bolthausen and Sznitman
[1998], Brunet et al. [2007], Neher and Hallatschek [2013].

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
$\lambda_{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 coalescence will lower the number of lineages below
$n$. The time $T_{n}$ is thus exponentially distributed with rate $\nu(n)$,
where $\nu(n)$ is the total rate of coalescence given $n$ lineages:

$\nu(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\left(1-\langle(1-\beta)^{n}\rangle-n\langle\beta(1-\beta)^{n-1}\rangle\right).$
(37)

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}\rightarrow\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}\geq\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 S4.

Figure S 3:  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 S 4:  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}=1\,000$. Right: Kingman coalescent with population
size $N=N_{e}=1\,000$.