Introduction

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

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

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

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

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

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

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

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

Results

1. Multi-strain SIR model

We describe the interaction of several viral strains and host immunity using a Susceptible/Infected compartmental model, similar to those used in [14, 15]. In the most general form, the model describes N variants of the virus labeled a ∈ {1… N} circulating among M groups of hosts with distinct exposure histories labeled i ∈ {1… M} (immune groups). These groups could be different age cohorts or could be geographically separated. For each group i, we define compartments and 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 ≤ , ≤ 1, and values of and 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 . Here, α 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 δ, we can thus write the dynamics for :

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

where the first term accounts for immunity gains (or loss of susceptibility) due to infections or cross-immunity while the second represents immune waning. This model introduced by Gog and Grenfell [14] assumes that immunity builds up through exposure and not only through infection. This explains that the change in is simply proportional to · , 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 .

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 , but as highly similar by others, leading to ≃ 1. Such a heterogeneous response by different immune systems has been observed experimentally in the case of influenza: in [19, 20] for instance, it was found that a given mutation in an influenza strain may allow escape from the antibodies of some individuals, i.e. low . while it had little effect for the serum of other individuals, i.e. high . 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 ≪ 1) while off-diagonal terms would be small (C_ij ≪ 1 for i ≠ j).

2. Invasion of an adaptive variant

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

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

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

with 0 ≤ b,f ≤ 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:

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

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

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

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 ≠ 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, β > 1/2 and the variant will be the dominant variant.

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

The variant thus increases initially exponentially until it has become sufficiently frequent that it starts having a substantial effect on the immunity landscape, before eventually settling into an equilibrium with the wild-type. The details of the equilibrium reached by the system in the absence of additional mutant variants is given in the SI. Figure 1 explores different scenarios numerically.

FIG. 1.

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 δ/α, 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₁ 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

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

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

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

Notably, the equilibrium frequency in the above examples does not depend on M and Equation 6 remains valid for ε = 0. This invariance is a consequence of the fact that for ε = 0, the variant and wild-type strains are completely equivalent in immune groups i > 1 and equilibrium is only determined by cross-immunity in group i = 1 (see SI). For small ε the equilibrium shifts slightly, but Equation 6 remains a good approximation.

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

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

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

3. Ultimate fate of the invading variant

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

FIG. 2.

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

The scenario illustrated in Fig. 2 suggests that many aspects of the variant dynamics can be approximated by a simple abstraction: If x is the frequency of a mutation A, a new variant has a probability x to appear on the background of A and thus carry A, and a probability 1 − x to not carry A. If new mutants emerged well separated in time with rate p, meaning that they reach equilibrium before the next variant emerges, and if new variants have a similar cross-immunity with all existing variants (see SI), the dynamics of x(t) are described by a particular random walk: in each time interval dt, a partial sweep of amplitude β occurs with probability ρdt · P_β(β), changing x in the following way:

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

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

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

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₀ should on average stay at this frequency, which means in practice that in a finite population it has a chance to reach 1 and fix, and a chance 1 − x₀ to reach 0 and vanish.

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

4. Abstraction as ‘expiring’ fitness advantage

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

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

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

The dynamics of this new model are represented in panel D of Figure 2, with an initial frequency x₀ ≪ 1 and an initial growth rate s₀ = 0.05. The initial growth of x is identical to a classical selective sweep of fitness s₀ (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 β given by (see SI)

This final value β depends only on the ratio between the initial fitness advantage s₀ and the rate of fitness decay ν. For a large enough s₀, β 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 β < 1 and remains there forever unless other variants appear.

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

5. Consequences for predictability and population dynamics

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

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

FIG. 3.

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

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

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 σ_i ∈ {0,1}. Fitness effects S_i are associated with mutations at each position, and the total fitness of a virus is given by F = ∑_iσ_is_i. At each generation, viruses with a fitness F expand by a factor e^F, and the next generation is constructed by sampling N individuals from the previous one. Following Eq. 11, mutational effects S_i decrease by an amount νx_i · 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 ρ, we pick one sequence position i that has no polymorphism and set the fitness effect of the corresponding mutation to an initial value S_i, with an amplitude drawn from probability distribution P_s and the sign chosen such that the mutation is adaptive. In practice, we use an exponential distribution P_s ∝ e^−s/s₀, meaning that the typical magnitude of initial fitness effects is described by only one parameter s₀. The corresponding distribution of partial sweep size is described in the SI. At the same time, we introduce the corresponding mutant in the population at a low frequency, picking its background genotype from a random existing strain. The behavior of the model is determined by (i) the distribution P_β of partial sweeps size depending on ν/s₀, and (ii) the ratio of the variant emergence rate and their growth rate ρ/s₀, which determines how often sweeps overlap and interact. The probability of two sweeps overlapping is defined in Section (B 2) of the SI.

We use this simulation to address the question of predictability: given the state of the population at generations 0, 1,…, t, can we predict its state at future times t + 1, …? Specifically, we ask whether we can predict the frequency x(t + Δt) of a variant A, given it is at frequency x at time t, as we did previously for influenza virus [12], see Fig. 3. The dynamics of isolated selective sweeps (ρ/s₀ ≪ 1, ν/s₀ ≪ 1) should be perfectly predictable: after an initial stochastic phase when the variant is very rare, its frequency grows monotonically to fixation. This predictability decreases with increasing ρ/s₀ due to clonal interference [27, 28], for example when an adaptive variant is outcompeted by an even more adaptive one. We also expect predictability to decrease with increasing ν/s₀ 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: ν ∈ [0, s₀/3, s₀, 3s₀] and low clonal interference ρ/s₀ = 0.05. The case ν = 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 ν/s₀, corresponding to a quicker decay of fitness, predictability gradually declines and becomes negligible for ν/s ≪ 1. Note that this matches quite well with the predictions from the random walk model where the average change in frequency 〈Δx〉 is null.

FIG. 4.

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

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

Finally, we use our simulation to investigate typical levels of diversity in the population and the time to the most recent common ancestor. One quantity that can easily be estimated from the random walk model is the average pairwise coalescence time T₂, 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₂ = 1/ρ〈β²〉P_β, which in neutral models of evolution would correspond to the effective population size N_e. A more detailed analysis of the coalescent process reveals that the random walk approximation corresponds to the so-called A-coalescent [29, 30] (see SI).

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₂ predicted by the random walk model. Each point in the figure corresponds to one simulation of long duration with a given distribution of partial sweep size P_β and a given ρ setting T₂, with darker color indicating a higher probability of overlap as computed in the SI. We find a good agreement between the empirical time to MRCA and the estimation from the random walk, at least as long the probability of overlap between successive partial sweeps is small (indicated by shading, see SI). With increasing overlaps, coalescence slows down and diversity increases: points in darker shades of red tend to have larger time to MRCA than what is expected from the distribution of β. This is expected intuitively: if another adaptive variant emerges before the previous one has reached its final frequency, it has a lower probability to land on the same background and thus tends to be in competition with the first variant. This leads to a smaller effective β which slows the dynamics.

Discussion

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

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

We have shown that depending on (i) the heterogeneity of immunity in the population, (ii) the asymmetry between backward and forward cross-immune recognition, and (Hi) waning or 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 [36]. On the other hand, most substitutions do not sweep to fixation but tend to meander in a quasi-neutral fashion [12], In the partial sweep scenario proposed here, diversity is dominated by immune escape mutations that are rapidly brought to macroscopic frequency by their initial growth advantage, but their ultimate fate is determined mostly by subsequent mutations.

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

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

Code availability

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

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

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

• Code to generate empirical frequency trajectories and their averages is available as scripts ‘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 [39]. We thank the teams involved in sample collection, sequencing, and processing of these data for their contribution to global surveillance of influenza virus circulation. A table acknowledging all originating and submitting laboratories is provided as supplementary information.

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

Acknowledgements

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

Declaration of Interests

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

A. SIR model

1. Equilibrium of the SIR model with one immune group

To help us compute the equilibrium reached by the SIR model, we introduce additional notation: the vectors = [S¹, …, S^N] and = [I¹, …, I^N] will respectively represent the compartments of hosts susceptible and infectious to each strain; the matrix K describes the cross-immunity; the vector 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 = δ/α. We thus have . Introducing this in the equation for the derivative of S^a, we obtain the following equilibrium:

An interesting remark is that at equilibrium, I is of order γ/δ ≪ 1. Note that the structure of K makes it invertible in most cases. Indeed, we impose K^aa = 1 and 0 ≤ K^ab < 1 for a ≠ b.

2. Equilibrium for two viruses with one immune group

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

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

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

where stands for [I^wt, I^m] and for [1,1]. It is straightforward to invert the crossimmunity matrix, and we obtain

Note that without cross-immunity, the number of infected by either virus would be 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:

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 , making it impossible for the wild-type to grow.

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

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

3. Equilibrium without the mutant

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

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:

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

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

If f = 0, the growth rate is α − δ, i.e. the one expected in a fully naive population.

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

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

In a given immune group i, the mutant growth rate is proportional to . 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 ≃ 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 as an exponential growth with a time dependent rate. In this case, the overall growth rate is given by the derivative of :

In particular, using the invasion scenario from the main text with ε = 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:

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

In the case of a non well-mixed population, i.e. arbitrary C_ij, it is not possible to write a pseudo-exponential growth rate as in Eq A10. However, it is clear that the initial growth rate will also be smaller as in the single group case since the mutant initially only grows in group i = l.

4. Equilibrium with M immune groups

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

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

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

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

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

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

where we have used C_ij = 1/M and to remove the sum on immune groups. Multiplying this equation by ν^a = I^a/I, we obtain

We now eliminate by using the expression = I_iν^a:

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

where the last expression is obtained by substracting the two previous ones. First, we see that for i > 1 any frequency vector ν is a solution since = 1 for all a, b. For i = 1 and defining ν^m = β and ν^wt = 1 − β and using the expression for K₁, we obtain

as claimed.

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

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

and ε > 0, the expressions above do not hold. However, if ε ≪ 1, the perturbation is small and we expect an equilibrium frequency close to β, which the case in Figure 1.

5. Realistic modeling of host’s immune state

The SIR model proposed in the main text relies on the assumption of immunity acquisition through exposure. This explains terms like −α ∑_bS^aK^abI^b in the derivative of S^a: acquiring immunity to a through cross-immunity K^abrequires 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⁰ = 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:

and three for the immune:

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 β is different but relatively close to the one from the main text.

6. Change in frequency when adding subsequent strains

This section shows that under certain condition, adding a new variant to the SIR model does not change the relative frequencies of previous variants. This is an important condition for the random walk of the main text to be valid.

Here is a quick summary of the results proved below. Adding a variant to the SI model involves adding a column and a row to the cross-immunity matrix, which can be given by two vectors. If these vectors only depend on one parameter, i.e. = b · and = f · , 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.

Figure S 1.

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

where 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

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 is now written as

where and 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 :

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

Below, we prove this equality under a condition for cross-immunity of the new variant and :

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 . To do that, we make use of the following formula for inverting a block matrix:

where we defined Λ = (D − CA⁻¹B)⁻¹. The following identities map to our problem:

We immediatly see that reduces to a scalar that we note A for more clarity. We also define the other scalar value A few manipulations give us the following for :

Multiplying this by results in

where we used the equalities .

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γ(1 − fμ). This implies that the relative frequencies of the original variants are conserved when adding a new one.

7. Case with intrinsic fitness effects

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

The equations of motion are now

Computing the equilibrium, we immediately obtain

where we have defined the following quantities:

The quantities h^a can be interpreted as a scaled growth rate of each variant given a fully susceptible population, and the matrix G combines the cross-immunity and the ratio of fitness values ϕ. 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 .

Inverting G and simplifying the equations a bit, we obtain

where we defined

Note the interesting structure of equations A22: for each variant, they involve a first term that depends only of the intrinsic growth rate of the mutant, and a second that involves cross-immunity and relative growth rate through ξ.

These equilibrium equations give us two conditions for the co-existence of the two variants,

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

• If the mutant has an intrinsic fitness disadvantage ξ < 1, it will only be able to invade if f < ξ. 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 ξ.

• Invertly, if the mutant is fitter and ξ > 1, the mutant is always able to invade. The wild-type only survives if b < ξ⁻¹, 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 ξ > 1 meaning ϕ^m > ϕ^wt, and the result is a full selective sweep.

8. Oscillations of the SI model

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

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

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

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

This is only valid if the terms in the square roots above are positive, which requires yto be small enough. In our setting, we assume γ ≪ αδ, so this will always hold.

From the eigenvalues we can compute:

• the rate of convergence to equilibrium: αγ/2δ. This means that convergence is slower for a smaller waning rate γ.

• the two oscillation frequencies that appear:

  

Note that since (1 − b)(1 − f)/(1 − bf) ≤ 1, we always have ω_low ≤ ω_high With the values of the main text α = 3, δ = 1, γ = 5 · 10⁻³ (units are inverse of generations), we obtain ω_high ≃ 0.016. If one generation is one week, this gives us a period , that is approximately one year.

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

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

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

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

which allows us to define the fitness in the SIR case: s = α(S^m − S^wt). At the beginning of the invasion, the initial growth rate is readily computed:

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

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

where the approximation is valid if 1 − b ≪ 1. This would give an expiry rate of fitness ν = α(I^wt + I^m) in the case of the SIR model.

These results can also be obtained in the case of immune groups. We then have the following expressions for s and ν at t ≃ 0:

where the second expression is again valid if 1 − b_i ≪ 1.

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

The expression of the partial sweep size becomes beta = ϵ_f/(ϵ_b + ϵ_f). Note that both e should remain smaller than 1, which is not guaranteed with exponential distributions. However, this is not problematic if μ and γ take small enough values.

These assumptions allow us to compute the distribution of beta:

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

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

Figure S 2.

B. Expiring fitness model and random walk

1. Clonal interference

In non recombining genomes with a large mutation rate, the appearance of many adaptive mutations in close succession leads to clonal interference. In this regime, beneficial mutations present on different background compete for fixation, and the success of a mutation does not depend only on its fitness effect but also on the global state of the population. For this reason, clonal interference causes a decrease in predictability: dynamics are not deterministic but rather depend on the precise structure of which mutation appeared on which background. For instance, a beneficial mutation that increases in frequency can be outcompeted by a fitter one before it has the time to fix, making the extrapolation of frequency trajectories difficult.

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

There is no fitness decay in this simulation, and the parameters N, s₀ and δf have numerical values such that neutral drift is mostly irrelevant to the dynamics of mutations. For example, without accounting for interference, the probability of fixation of a mutation of effect s₀ when it is introduced at frequency δf is P_fix(ϴf) ≃ 1 − e⁻⁴⁰. As a consequence, the two parameters governing the dynamics are the scale of fitness effects s₀ and the rate of introduction of beneficial mutations ρ. The ratio of these two quantities represents the amount of clonal interference: at low ρ/s₀, mutations are mostly independent, while at high ρ/s₀ they strongly interfere.

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

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

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

Figure S 3.

First, we prove the expression for the amplitude of the partial sweeps. We divide the equation for ẋ by the one for ṡ to obtain

This immediately gives us

with γ constant. At t = 0, we have s = s₀ and x = x₀ ≪ 1, while for t → ∞ we have s = 0 and x = β to be determined. From the t = 0 case we obtain γ = (x₀ − 1)e^s₀/ν, and from t → ∞ we get β = 1 + γ. Assuming x₀ ≪ 1, we obtain the result of the main text:

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:

where β is a function of s and ν and x₀ = x(t = 0). This is simply the expression of a logistic growth starting at x₀, and saturating at β. From there, we compute the time T_r(s) it takes a partial sweep of initial fitness s to reach a frequency rβ with x₀β⁻¹ < r < 1. We quickly find

We now consider that two consecutive partial sweeps of initial fitnesses s₁ and s₂ overlap if the first one is not yet at frequency rβ₁ while the second one is already at (1 − r)β₂. In the figure of the main text, we use r = 3/4: an overlap occurs if the first sweep is not yet at 3/4 of its final value while the second one is already at 1/4 of its final value. Thus, for an overlap to occur, we need the time τ between the two partial sweeps to be smaller than T_r(s₁) − T_1-r(s₂). For sweeps happening at rate σ, this has probability 1 − exp (−σ(T_r(s₁ − T_1−r(s₂))). Since the two sweeps have random initial fitness effects, we find that the overall probability for two consecutive sweeps to overlap is

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

3. Distribution of partial sweep size β

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

A first interesting case is when fitness effects are exponentially distributed, with parameter s₀:

This is the distribution we use in most of the population simulations. We compute the corresponding distribution of β in a straightforward way:

Taking the derivative with respect to x, we obtain

This distribution can accomodate various shape: for ν/s₀ > 1 it peaks at 0, and for ν/s₀ < 1 it peaks at 1. We can also compute the following formular for the moments of β:

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 β, and we decide to sample β 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:

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

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 = ε · m² with ε = 10⁻⁵

• 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 β for each new adaptive mutation, and then computing s by using Equation 12 from the main text. For each distribution P_β, the simulation is performed for 6 values of ρ ∈ [0.003, 0.01, 0.018, 0.032, 0.056, 0.1].

4. Random walk: monotonous trajectories

Here we compute the probability that in the random walk defined in the main text, a trajectory starting at x₀ 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 β.

If the random walk always goes down, its position at time step t will be x_t = (1 − β)^tx₀. Since the probability of going down is 1 − x_t, the probability of always going down is

We simplify this expression by taking the logarithm and assuming that (1 − β)^t ≪ 1 for t ≥ 1:

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

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

Figure S 4.

5. Coalescent

Consider a partial sweep happening between generations t and t + 1, with probability ρ. One individual A in generation t will then have βN children in generation t + 1. Any individual in generation t + 1 has a probability β of having A as a direct ancestor, and a probability 1 − β 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 β^k. Averaging over P^β, we find the probability of k specificlineages to have a common ancestor in the previous generation:

Another useful quantity is the probability γ_n(k) that given n lineages, a particular set of exactly k lineages merge one generation back. If a partial sweep of known amplitude β took place, this requires the set of k lineages to merge at this generation, with probability β^k, and that the other n − k do not merge, with probability (1 − β)^n−k.

This turns out to be the definition of the Λ-coalescent with Λ(β) ∝ β²P(β) [29, 30]. The A-coalescent is a general model for genealogies of multiple mergers. We mention two interesting subcases:

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

• if Λ(β) is uniform in [0,1], meaning P(β) ∝ β⁻², we obtain the Bolthausen-Sznitman coalescent [30, 41], which is used to describe the genealogy of populations under strong selection [41–43].

Finally, we derive a few more properties of the partial sweep coalescent and show the explicit link to Kingman’s when β ≪ 1. Using the γ_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 ν/(n), where ν/(n) is the total rate of coalescence given n lineages:

Since we have , we finally obtain

Figure S 5.

With n〈β〉 ≪ 1, we now exactly recover the Kingman coalescent. For simplicity, we assume a constant β and expand T_n up to the second order in nβ, to obtain

These are the times expected for the Kingman coalescent with population size N_e = 1/ρβ².

In the high n limit, we also obtain T_n → ρ⁻¹, since quantities of the type (1 − β)ⁿ vanish. This is expected as coalescences only take place when a partial sweep happens, with rate ρ. It is another qualitative difference with the Kingman coalescent: since T_n ≥ ρ⁻¹ for all n, one must wait a time ~ρ⁻¹ to observe the first coalescence even in large trees. The shortest branches will thus always be of order ρ⁻¹. In contrast, in the Kingman process, the shortest branches vanish when the number of lineages n increases. This difference is clearly visible when looking at terminal branches of trees in Figure S6.

C. Supplementary figures

Figure S 6.

Figure S 7.