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

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 [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 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 [21].

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

B. 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 α ∈ {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 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 , 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, a 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 :

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 and , we can sum over the immune groups in the above equation to obtain . We consequently define the reproductive number of strain a as . I^a grows when R^a > 1 and declines when R^a < 1. The equilibrium susceptibility is S^a = δ/α, such that R^a = 1.

When a host of group i is infected by strain b, it also has a chance to gain immunity to strain a. 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 . 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 a, 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 ∈ {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. For the susceptibles we find S^a = δ/α for all strains a, while 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 [22].

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.

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 δ/α 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.

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.

FIG. 1.

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. 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 ε = 0, the equilibrium frequency does not depend on M and Equation 6 remains valid, while for small ε 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 δ/α, 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 β, 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 [21]. 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.

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 β₂, β₃,… This process is shown in panel A of Figure 2, using the SIR model to simulate up to N = 7 variants. For the sake of illustration, it shows a simple scenario where the initial variant appears at time 0 in a homogeneous wild-type population, and subsequent mutants appear at regular time intervals. Simulations are performed using M = 10 immune groups, resulting in a slight overshoot of the equilibrium frequency for each trajectory.

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

FIG. 2.

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 ρ, 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 [23].

Examples of trajectories from the random walk are shown on panel C of Figure 2, all initially starting around x₀ ≃ 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 Δ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 [24]. 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 x₀ to reach 1 and fix, and a chance 1 – x₀ to reach 0 and vanish.

On the other hand, the second moment closely reminds of the theory of neutral drift [25]: 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 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 . 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. 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:.

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₀ ≪ 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 v. 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. 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.

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. Looking at the corresponding frequency trajectories, we can see on average whether they keep on growing beyond x^*. Figure 3 shows these trajectories for the amino acid substitutions in the HA protein of A/H3N2 influenza, using data from 2000 to 2023, and the SARS-CoV-2 genome using data from 2020 to 2023. Panels on the left show all trajectories that reached x^* = 0.4, with their average displayed in black. The panels on the right show the average trajectory for different threshold values x^* between 0.1 and 0.8.

While the dynamics of the variants of the two viruses can not be compared directly due to vastly different sampling intensities and different rates of adaptation, the qualitative patterns differ strikingly. In the case of influenza, trajectories of seemingly adaptive mutations show little inertia and on average hover around x^* instead of growing. This surprising result is in line with the study in [12] 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. 8, 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 si are associated with mutations at each position, and the total fitness of a virus is given by . 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 vx_i · s_i, where x_i is the frequency at which mutation i is found in the population.

FIG. 3.

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 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₀. 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 v/s₀, and (ii) the ratio of the variant emergence rate and their growth rate ρ/s₀, 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, …, 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, v/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 [26], for example when an adaptive variant is outcompeted by an even more adaptive one. We also expect predictability to decrease with increasing v/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: v ∈ [0, s₀/3, s₀, 3s₀] and low clonal interference ρ/s₀ = 0.05. The case v = 0 corresponds to a classical traveling wave scenario with constant fitness effects, and, as expected, is the most predictable: the average trajectory rises well above 0.5. For larger values of v/s₀ , corresponding to a quicker decay of fitness, predictability gradually declines and becomes negligible for v/s ≫ 1. Note that this matches quite well with the predictions from the random walk model where the average change in frequency is null.

FIG. 4.

To explore parameter space more systematically, we quantify predictability as the probability of fixation pf_ix of rising variants that cross threshold x^*. In a perfectly predictable scenario with well separated selective sweeps, pf_ix 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 v/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 v). Increasing either ρ/s₀ or v/s₀ reduces pf_ix towards the dashed diagonal corresponding to pf_ix = x^*. However, as observed previously [12], in the classic scenario with stable fitness effects v/s₀ = 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₂, 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, , 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 [27, 28] (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₂. 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.

C. 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, 29]. The A/H3N2 subtype in particular undergoes rapid antigenic change through frequent substitutions in prominent epitopes on its surface proteins [30–33]. 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 (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 [34]. 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 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” [35]), 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 [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 [36]. 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 ‘flufixation.py’ and ‘sc2fixation.py’ in https://github.com/nextstrain/flu_frequencies on branch ‘fixation’.

Data availability

Sequence data of influenza viruses was obtained from GISAID [37]. 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.

Appendix

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 and 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 = 1.

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

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 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. and , 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 γ 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.

6. 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 A21: 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.

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

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 γ to 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.

B. Other

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:

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₀/v, 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 v 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 rft1 while the second one is already at (1 − r)β₂. 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 τ 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.

2. Distribution of partial sweep size β

This section discusses the distribution of the size of partial sweeps β 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₀:

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(1 − 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 11 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].

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₀ 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 S2. The same figure also shows that this probability is relatively high, even for x₀ far from either boundary.

4. 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:

Figure S 2.

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(β) [27, 28].

The Λ-coalescent is a general model for genealogies of multiple mergers. We mention two interesting subcases:

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

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

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 v(n), where v(n) is the total rate of coalescence given n lineages:

Since we have , we finally obtain

With , 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 S4.

Figure S 3.

Figure S 4.