We describe the interaction of several viral strains and host immunity using a Susceptible/Infected compartmental model, similar to those used in Gog and Grenfell, 2002; Yan et al., 2019. In the most general form, the model describes ${\displaystyle N}$ variants of the virus labeled ${\displaystyle a\in \{1\dots N\}}$ circulating among ${\displaystyle M}$ groups of hosts with distinct exposure histories labeled ${\displaystyle i\in \{1\dots M\}}$ (immune groups). These groups could be different age cohorts or could be geographically separated. For each group ${\displaystyle i}$, we define compartments ${\displaystyle I_i^a}$ and ${\displaystyle S_i^a}$ as, respectively, the number of individuals of this group infected or susceptible to strain ${\displaystyle a}$. We assume that the total population of hosts is 1 so that we always have ${\displaystyle 0\le I_i^a,S_i^a\le 1}$, and values of ${\displaystyle I_i^a}$ and ${\displaystyle S_i^a}$ 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 ${\displaystyle i}$ susceptible to variant ${\displaystyle a}$ are infected by ${\displaystyle a}$ is ${\displaystyle \alpha S_i^a\sum _{j=1}^M{C}_{ij}{I}_j^a}$. Here, ${\displaystyle \alpha }$ is an overall infection rate while ${\displaystyle C_{ij}}$ represents the probability of an encounter between individuals of groups ${\displaystyle i}$ and ${\displaystyle j}$. Thus, the above rate takes into account infections with strains ${\displaystyle a}$ caused by hosts of all groups. Considering that infected hosts recover at rate ${\displaystyle \delta }$, we can thus write the dynamics for ${\displaystyle I_i^a}$:

When a host of group ${\displaystyle i}$ is infected by strain ${\displaystyle b}$, it gains immunity against the infecting strain ${\displaystyle b}$, but also to other strains ${\displaystyle a}$ with probability ${\displaystyle 0\le K_i^{ab}\le 1}$. Thus, ${\displaystyle S_i^a}$ decreases at a rate proportional to ${\displaystyle K_i^{ab}}$ and to the number of hosts infected by ${\displaystyle b}$ for every strain ${\displaystyle b}$. Since susceptibles to ${\displaystyle a}$ are depleted by infections from other strains, the dynamics of all strains are coupled. This coupling is determined by the matrices ${\displaystyle {\mathbf{K}}_i}$ of dimension ${\displaystyle N\times N}$, which in general differ between groups ${\displaystyle i}$ with different prior exposure history. Additionally, the waning of immunity at a rate ${\displaystyle \gamma }$ causes immune hosts to re-enter the susceptible compartment. We can now write the dynamics of ${\displaystyle S_i^a}$ 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, 2002 assumes that immunity builds up through exposure and not only through infection. This explains that the change in ${\displaystyle S_i^a}$ is simply proportional to ${\displaystyle S_i^a\cdot I_j^b}$, regardless of the susceptibility of hosts to strain ${\displaystyle b}$. Alternative models that require infection for acquisition of immunity have qualitatively similar dynamics, but are mathematically more complex (Appendix 1.5). We also represented loss of susceptibility to ${\displaystyle a}$ due to exposure to ${\displaystyle a}$ using a trivial cross-immunity term ${\displaystyle K_i^{aa}=1}$.

An important property of our model is that the probability of generating cross-immunity can differ between groups. The motivation is that strains ${\displaystyle a}$ and ${\displaystyle b}$ may be perceived as antigenically different by some immune systems, leading to a low ${\displaystyle K_i^{ab}}$, but as highly similar by others, leading to ${\displaystyle K_i^{ab}\simeq 1}$. Such a heterogeneous response by different immune systems has been observed experimentally in the case of influenza: in Lee et al., 2019; Welsh et al., 2023 for instance, it was found that a given mutation in an influenza strain may allow escape from the antibodies of some individuals, i.e., low ${\displaystyle K_i^{ab}}$, while it had little effect on the serum of other individuals, i.e., high ${\displaystyle K_i^{ab}}$. Heterogeneous immune response could be caused by varying histories of strain exposure for different individuals, for instance, due to differences 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 ${\displaystyle C_{ij}=1/M}$. If immune groups are shaped by geographic differences in exposure, the connectivity would be close to 1 on the diagonal (${\displaystyle 1-C_{ii}\ll 1}$) while off-diagonal terms would be small (${\displaystyle C_{ij}\ll 1}$ for ${\displaystyle i\ne j}$).

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 ${\displaystyle (a,b)}$ describing strains will take values ${\displaystyle \{wt,v\}}$. We consider that the two strains share the same infectivity rate ${\displaystyle \alpha }$, which amounts to say that they would have the same reproductive rate in a fully naive population. The case where the two strains differ in intrinsic fitness is explored in detail in Appendix 1.7. 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 ${\displaystyle M=1}$. We can thus skip the indices ${\displaystyle i,j\in \{1\dots M\}}$, and we only have one cross-immunity matrix ${\displaystyle \mathbf{K}}$ that we parametrize as

with ${\displaystyle {\displaystyle 0\le b,f\le 1}}$ . ${\displaystyle b}$ quantifies the amount of ‘backward’-immunity to the wild-type caused by the variant: a large ${\displaystyle b}$ means that it is likely that an infection by the variant causes immunity to the wild-type. Conversely, ${\displaystyle f}$ quantifies the ‘forward’-immunity: infections by the wild-type causing immunity to the variant. If ${\displaystyle f=b=1}$, the two strains are antigenically indistinguishable, and thus essentially identical for the model. Conversely, if ${\displaystyle 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 ${\displaystyle a}$ as ${\displaystyle R^a=\alpha S^a/\delta }$. ${\displaystyle I^a}$ grows when ${\displaystyle R^a>1}$ and declines when ${\displaystyle R^a<1}$. The equilibrium susceptibility is, therefore, ${\displaystyle S^a=\delta /\alpha }$, such that ${\displaystyle {\displaystyle R^a=1}}$. On the other hand, the equilibrium prevalence is determined by the inverse of the cross-immunity matrix ${\displaystyle \mathbf{K}}$:

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

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

In the case where ${\displaystyle b=f}$, the variant will ultimately settle at frequency 1/2. This includes the case where ${\displaystyle b=f=0}$, where the two strains are completely independent and do not interact. On the contrary if ${\displaystyle b\ne f}$, the final frequency of the variant can in principle be anywhere between 0 and 1. For example if ${\displaystyle 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, ${\displaystyle \beta >1/2}$ and the variant will be the dominant variant.

We are primarily interested in an ‘invasion’ scenario where only the wild-type is initially present in the population, that is ${\displaystyle I^v=0}$ at ${\displaystyle 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 ${\displaystyle \delta /\alpha }$ unless ${\displaystyle f=1}$, As a result, the growth rate of the variant is initially positive and given by

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 Appendix 1.1. Figure 1 explores different scenarios numerically.

Figure 1

Download asset Open asset

The top row of Figure 1 shows the dynamics of the model after the introduction of the variant in a homogeneous population (${\displaystyle M=1}$). As expected, the number of infections by the variant initially rises while the number of susceptibles ${\displaystyle S^v}$ decreases. However, as ${\displaystyle S^v}$ goes below the critical value ${\displaystyle \delta /\alpha }$, ${\displaystyle 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 Appendix 1.8.

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 ${\displaystyle 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 ${\displaystyle M=10}$ groups. For group ${\displaystyle i=1}$, the cross-immunity matrix ${\displaystyle K_1}$ takes the same form as in the previous scenario, given by Equation 3. However, for other groups, we assume that the two strains are virtually identical, with the cross-immunity having the form

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

One consequence of many groups that are indifferent to the variant is that globally the excess susceptibility to the variant is lower. If mixing is rapid, the initial growth rate of the variant is smaller by a factor of ${\displaystyle 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 ${\displaystyle M}$ groups, respectively, for the rapid and slow mixing cases. In both scenarios, the initial number of hosts susceptible to the variant are now closer to ${\displaystyle \delta /\alpha }$. When mixing is fast, the frequency of the variant initially resembles a standard selective sweep (dashed line in Figure 1) before saturating, while dynamics are more complicated for the slow mixing case. Either way, the main effect of the immune groups is that the overshoot past the equilibrium is much smaller and the 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 ${\displaystyle M}$ and Equation 6 remains valid for ${\displaystyle \epsilon =0}$. This invariance is a consequence of the fact that for ${\displaystyle \epsilon =0}$, the variant and wild-type strains are completely equivalent in immune groups ${\displaystyle i>1}$ and equilibrium is only determined by cross-immunity in group ${\displaystyle i=1}$ (Appendix 1.4). For small ${\displaystyle \epsilon }$ 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 ${\displaystyle \beta }$, their frequency will of course continue to change due to the 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 gives rise to the partial sweep is reminiscent of consumer resource models Tikhonov and Monasson, 2018; Good et al., 2018, 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.

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 ${\displaystyle {\beta }_2,{\beta }_3,\dots }$ sampled from some distribution ${\displaystyle P_{\beta }}$. This process is shown in panel A of Figure 2, using the SIR model to simulate up to ${\displaystyle 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 ${\displaystyle M=10}$ immune groups, resulting in a slight overshoot of the equilibrium frequency for each trajectory.

Figure 2

Download asset Open asset

Here, we focus on the mutation or set of mutations ${\displaystyle A}$ that defines the initial variant. The initial growth rate advantage given by ${\displaystyle A}$ eventually disappears, meaning that after some times 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 ${\displaystyle A}$, or on the background of the initial variant in which case they do carry ${\displaystyle A}$. If we suppose that recombination is negligible, the frequency of ${\displaystyle 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 ${\displaystyle A}$ (resp. not carrying ${\displaystyle A}$). The thick line in between the red and blue surfaces indicate the frequency at which mutation ${\displaystyle A}$ is found, and in practice moves up and down randomly.

The scenario illustrated in Figure 2 suggests that many aspects of the variant dynamics can be approximated by a simple abstraction: if ${\displaystyle x}$ is the frequency of a mutation ${\displaystyle A}$, a new variant has a probability ${\displaystyle x}$ to appear on the background of ${\displaystyle A}$ and thus carry ${\displaystyle A}$, and a probability ${\displaystyle 1-x}$ to not carry ${\displaystyle A}$. If new mutants emerge well separated in time with rate ${\displaystyle \rho }$, meaning that they reach equilibrium before the next variant emerges, and if new variants have a similar cross-immunity with all existing variants (see Appendix 1.6), the dynamics of ${\displaystyle x(t)}$ are described by a particular random walk: in each time interval ${\displaystyle \text{d}t}$, a partial sweep of amplitude ${\displaystyle \beta }$ occurs with probability ${\displaystyle \rho \text{d}t\cdot P_{\beta }(\beta )}$, changing ${\displaystyle x}$ in the following way:

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

Examples of trajectories from the random walk are shown on panel C of Figure 2, all initially starting around ${\displaystyle x_0\simeq 1/2}$. Two trajectories converge monotonically to 0 and 1. This is a consequence of one interesting property of Equation 9: the probability for ${\displaystyle \mathrm{\Delta }x}$ to be positive increases with ${\displaystyle x}$, but the magnitude of the upward steps decreases as ${\displaystyle 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 Appendix 2.4). On the other hand, steps away from the closest boundary are unlikely but much larger, resulting in ‘jack-pot’ events (Hallatschek, 2018). This can be seen in 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 ${\displaystyle \mathrm{\Delta }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 ${\displaystyle A}$ is not biased toward either fixation or loss, regardless of the frequency that the initial partial sweep brought it to. For instance, a mutation seen at frequency ${\displaystyle x_0}$ should on average stay at this frequency, which means in practice that in a finite population, it has a chance ${\displaystyle x_0}$ to reach 1 and fix, and a chance ${\displaystyle 1-x_0}$ to reach 0 and vanish.

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

In general, the dynamics of the SIR model proposed in Equations 1&2 depend on the interactions between ${\displaystyle N}$ strains through an ${\displaystyle N\times N}$ cross-immunity matrix. While this model is useful to give a mechanistic explanation of partial sweeps, it is in general impractical to analyze and numerically simulate for many 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 over-damped 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 ${\displaystyle 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 ${\displaystyle x}$ of the variant:

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

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

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

It is important to state that the main aim of this effective model is to qualitatively reproduce the phenomenology of the SIR, and in particular the partial sweeps, while being is 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 ${\displaystyle 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. Appendix 1.9 discusses in more detail the links between the parameters of the two models and the fundamental differences. While we can express the rate ${\displaystyle \nu }$ at which the growth rate declines in terms of the parameters of the simplest SIR models, for models with many groups or with oscillatory dynamics, the decay rate of the growth advantage should be interpreted as an effective parameter that captures a generic effect of reduced growth with increasing circulation.

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

In Figure 3, we reproduce some of our results of Barrat-Charlaix et al., 2021 and extend them to SARS-CoV-2. To quantify predictability, we ask the following question: given the state of a viral population at times ${\displaystyle 0,1,\dots ,t}$, what can we say about variant frequencies at times ${\displaystyle t+1,\dots }$? 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 ${\displaystyle x^{\star }}$. 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 ${\displaystyle x^{\star }}$. Figure 3 shows these trajectories for the amino acid substitutions in the HA protein of A/H3N2 influenza, using data from 2000–2023,, and the SARS-CoV-2 genome using data from 2020–2023. Panels on the left show all trajectories that reached ${\displaystyle x^{\star }=0.4}$, with their average displayed in black. The panels on the right show the average trajectory for different threshold values ${\displaystyle x^{\star }}$ between 0.1 and 0.8.

Figure 3 with 1 supplement see all

Download asset Open asset

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 ${\displaystyle x^{\star }}$ instead of growing. This surprising result is in line with the study in Barrat-Charlaix et al., 2021 which used data from the period 2000–2018. On the other hand, trajectories of SARS-CoV-2 mutations show a much smoother behavior with steady growth beyond ${\displaystyle x^{\star }}$. On longer timescales, however, we observe a systematic decrease in frequency: this is explained by the particular initial dynamics of SARS-CoV-2, where new variants arose at a rapid pace and replaced old ones. This process is often called clonal interference and reduces long-term predictability.

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 ${\displaystyle \beta }$ of its initial partial sweep. Subsequent allele dynamics and diversity are governed by an anomalous coalescent process driven by the random walk defined in Equation 9, leading to little predictability of evolution. This abstraction should hold when partial sweeps are instant and do not overlap, meaning that the rate ${\displaystyle \rho }$ at which new variants emerge is small compared to their initial growth rate ${\displaystyle s_0}$.

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

We simulate the emergence of adaptive variants in the following way. At a constant rate ${\displaystyle \rho }$, we pick one sequence position ${\displaystyle i}$ that has no polymorphism and set the fitness effect of the corresponding mutation to an initial value ${\displaystyle s_i}$, with an amplitude drawn from probability distribution ${\displaystyle P_s}$ and the sign chosen such that the mutation is adaptive. In practice, we use an exponential distribution ${\displaystyle P_s\propto e^{-s/s_0}}$, meaning that the typical magnitude of initial fitness effects are described by only one parameter ${\displaystyle s_0}$. The corresponding distribution of partial sweep size is described Appendix 2.3. 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 ${\displaystyle P_{\beta }}$ of partial sweeps size depending on ${\displaystyle \nu /s_0}$, and (ii) the ratio of the variant emergence rate and their growth rate ${\displaystyle \rho /s_0}$, which determines how often sweeps overlap and interact. The probability of two sweeps overlapping is defined in Appendix 2, Equation 42.

We use this simulation to address the question of predictability: given the state of the population at generations ${\displaystyle 0,1,\dots ,t}$, can we predict its state at future times ${\displaystyle t+1,\dots }$? Specifically, we ask whether we can predict the frequency ${\displaystyle x(t+\mathrm{\Delta }t)}$ of a variant ${\displaystyle A}$, given it is at frequency ${\displaystyle x}$ at time ${\displaystyle t}$, as we did previously for the influenza virus Barrat-Charlaix et al., 2021, see Figure 3. The dynamics of isolated selective sweeps (${\displaystyle \rho /s_0\ll 1}$, ${\displaystyle \nu /s_0\ll 1}$) should be perfectly predictable: after an initial stochastic phase when the variant is very rare, its frequency grows monotonically to fixation. This predictability decreases with increasing ${\displaystyle \rho /s_0}$ due to clonal interference (Schiffels et al., 2011; Strelkowa and Lässig, 2012), for example when an adaptive variant is outcompeted by an even more adaptive one. We also expects predictability to decrease with increasing ${\displaystyle \nu /s_0}$ since sweeps are then partial and their ultimate fixation is 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 ${\displaystyle x^{\ast }}$. The results are shown in panel A of Figure 4, where we show the average ${\displaystyle x(t)}$ of rising frequency trajectories after crossing the threshold ${\displaystyle x^{\ast }=0.5}$. We use three rates of fitness decay: ${\displaystyle \nu \in [0,s_0/3,s_0,3s_0]}$ and low clonal interference ${\displaystyle \rho /s_0=0.05}$. The case ${\displaystyle \nu =0}$ corresponds to a classical traveling wave scenario with constant fitness effects, and, as expected, is the most predictable: the average trajectory rises well above 0.5. For larger values of ${\displaystyle \nu /s_0}$, corresponding to a quicker decay of fitness, predictability gradually declines and becomes negligible for ${\displaystyle \nu /s\gg 1}$. Note that this matches quite well with the predictions from the random walk model where the average change in frequency ${\displaystyle ⟨\mathrm{\Delta }x⟩}$ is null.

Figure 4 with 1 supplement see all

Download asset Open asset

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

In panels B, C, and D of Figure 4, the probability of fixations are shown for three values of ${\displaystyle \rho /s_0}$ and four values of ${\displaystyle \nu /s_0}$. Clonal interference increases when going from left to right among these panels (increasing ${\displaystyle \rho /s_0}$), while the intensity of fitness decay increases when going from blue to red curves (increasing ${\displaystyle \nu }$). Increasing either ${\displaystyle \rho /s_0}$ or ${\displaystyle \nu /s_0}$ reduces ${\displaystyle p_{fix}}$ towards the dashed diagonal corresponding to ${\displaystyle p_{fix}=x^{\ast }}$. However, as observed previously (Barrat-Charlaix et al., 2021), in the classic scenario with stable fitness effects ${\displaystyle \nu /s_0=0}$ considerable predictability remains even in cases of strong interference (blue curve in panel D and Figure 4—figure supplement 1). The strong interference setting is explored in more detail in Appendix 2.1 up to values ${\displaystyle \rho /s_0\simeq 30}$, using similar simulations but without expiring fitness ${\displaystyle \nu =0}$. Figure 4—figure supplement 1 shows that even in these cases of strong interference, ${\displaystyle 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 ${\displaystyle T_2}$, that is the typical times separating two random strains from their most recent common ancestor (MRCA). In Appendix 2.5, we show that under the random walk approximation ${\displaystyle T_2=1/\rho ⟨{\beta }^2{⟩}_{P_{\beta }}}$, which in neutral models of evolution would correspond to the effective population size ${\displaystyle N_e}$. A more detailed analysis of the coalescent process reveals that the random walk approximation corresponds to the so-called ${\displaystyle \mathrm{\Lambda }}$-coalescent Schweinsberg, 2000; Berestycki, 2009.

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

• Figures in the main text can be reproduced using a set of notebooks at https://github.com/PierreBarrat/ExpiringFitnessFigures (copy archived at Barrat, 2024a).

• Code for the simulations of the SIR model is available at https://github.com/PierreBarrat/PartialSweepSIR.jl (copy archived at Barrat, 2024b).

• Code for the population dynamics simulation is available at https://github.com/PierreBarrat/WrightFisher.jl (copy archived at Barrat, 2024c).

• 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’ (Neher, 2024).

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

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