Phylodynamic theory of persistence, extinction and speciation of rapidly adapting pathogens

A model of an antigenically evolving pathogen population needs to account for cross-immunity between strains and the evolution of antigenically novel strains. We use an extension of the standard multi-strain SIR model (Gog and Grenfell, 2002). The fraction of individuals $I_a$ infected with viral strain $a$ changes according to

where $\beta$ is the transmissibility, $S_a$ is the population averaged susceptible to strain $a$, $\nu$ is the recovery rate, and $\gamma$ is the population turnover rate. The fraction $R_a$ of the population recovered from infection with strain $a$ changes according to

Our focus here is on antigenically evolving pathogens that reinfect an individual multiple times during its life-time, we shall ignore population turnover and set $\gamma =0$ right away to simplify presentation.

The dynamics of $I_a$ depends on the average susceptibility of the host population $S_a={⟨S_a(i)⟩}_i$, while the susceptibility $S_a(i)$ of host $i$ depends on the host’s history of previous infections. A plausible representation of the history dependence of susceptibility at the level of individuals has a product form (Wikramaratna et al., 2015)

where ${\sigma }_b(i)$ is one or zero depending on whether host $i$ has or has not been previously infected with strain $b$. Matrix $K_{ab}\le 1$ quantifies the cross-immunity to strain $a$ due to prior infection with strain $b$. Thus, Equation 3 expresses the susceptibility $S_a$ in terms of a product of attenuation factors each arising from a prior infection by a different strain $b$. A simple, but adequate approximation for the population averaged susceptibility is provided by replacing ${\sigma }_b(i)$ in the product in Equation 3 by the fraction of the population $R_b$ that recovered from infection with strain $b$:

This corresponds to the ‘order one independence closure’ by Kryazhimskiy et al. (2007) and is known as Mean-Field approximation in physics (Weiss, 1907; Landau and Lifshitz, 2013). The Mean-Field approximation here corresponds to ignoring correlations between subsequent infection in the individual histories. Approximating the product by the exponential is justified because the total fraction of the host population infected by any single strain in the endemic regime is typically small (Yang et al., 2015). A detailed derivation of Equation 4 and more detailed discussion of approximations is given in Appendix 1. While the original formulation of immunity in Equation 3 is based on the infection history of individuals (Andreasen et al., 1997), the population average over the factorized distribution of histories relates the model to status based formulations (Gog and Grenfell, 2002). While some differences between status and history-based models have been reported (Ballesteros et al., 2009), others have shown that different model types have similar properties (Ferguson and Andreasen, 2002). The differences between these models and approximations are small compared to the crudeness with which these simple mathematical models capture the complex immunity profile of the human population. A model similar to ours has been successfully applied to influenza virus evolution (Luksza and Lässig, 2014).

We note that differentiating Equation 4 with respect to time defines the equation governing the dynamics of population average susceptibility

which is exactly the same as the dynamics of susceptibility in Gog and Grenfell (2002) and Luksza and Lässig (2014) in the limit of negligible population turnover $\gamma /\nu \ll 1$.

New strains are constantly produced by mutation with rate $m$. The novel strain will differ from its parent at one position in its genome. Following Luksza and Lässig (2014), we assume that cross-immunity decays exponentially with the number of mutations that separate two strains:

where $|a-b|$ denotes the mutational distance between the two strains, $d$ denotes the radius of cross-immunity measured in units of mutations. Antigenic space is thereby assumed to be high dimensional and antigenic distance is proportional to genetic distance in the phylogenetic tree (Neher et al., 2016). The parameter $\alpha \le 1$ quantifies the reduction of susceptibility to reinfection by the same strain and hence the overall strength of protective immunity. We shall set $\alpha =1$ corresponding to perfect protection here for simplicity of presentation. Our analysis below applies equally well to the more realistic case of $\alpha <1$, since in our approximation this parameter can be eliminated by rescaling $R_a$ and $I_a$ and ultimately merely renormalizes the host population size, which serves as one of the ‘control parameters’ in our analysis.

Cross-immunity and the mutation/diversification process are illustrated in Figure 2. An infection with a particular strain (center of the graph) generates a cross-immunity footprint (shaded circles). Mutation away from the focal strain reduces the effect of existing immunity in the host population, but complete escape requires many mutations. Hence closely related viruses compete against each other for susceptible individuals.

Figure 2

Download asset Open asset

The above model was formulated in terms of the deterministic Equations 1-4. The actual dynamics, however, is stochastic in two respects: (i) antigenic mutations are generated at random with rate $m$ and (ii) stochasticity of infection and transmission. The latter can be captured by interpreting the terms in Equation 1 as rates of discrete transitions in a total population of $N_h$ hosts. This stochasticity is particularly important for novel mutant strains that are rare. Most rare strains are quickly lost by chance even if they have a growth advantage due to antigenic novelty. To account for stochasticity in a computationally efficient way, we employ a clone-based hybrid scheme where mutation and the dynamics of rare mutants are modeled stochastically, while common strains follow deterministic dynamics, see Materials and methods (Clone-based simulations).

We will use the recovery rate $\nu$ to set the unit of time, fixing $\nu =1$ in rescaled units. The remaining parameters of the model are (1) the transmission rate $\beta$ - in our units the number of transmission events per infection and hence equal to the basic reproduction number $R_0$, (2) the mutation rate $m$, (3) the range of cross-immunity $d$ measured as the typical number of mutations needed for an $e$-fold drop of cross-inhibition, and (4) the host population size $N_h$.

Before proceeding with a quantitative analysis we discuss different behaviors qualitatively. Figure 3A shows several trajectories of prevalence $I_{tot}={\sum }_a{I}_a$ (i.e. total actively infected fraction) for several different parameters. Depending on the range of cross-immunity, the pathogen either goes extinct after a single pandemic (red line) or settles into a persistently evolving state, the Red Queen State (RQS) traveling wave (Van Valen, 1973 In large populations the RQS exhibits oscillations in prevalence. As we will discuss further below, the RQS state is transient and either goes extinct after some time or splits into multiple antigenically diverging lineages that propagate independently. To quantitatively understand the dependence on parameters, we will further simplify the model and establish a connection to models of rapid adaptation in population genetics. Figure 3BC shows parameter regimes corresponding to distinct qualitative behaviors. The relevant parameters are three combinations of the population size $N_h$, the selection coefficient of novel mutations $s$, the mutation rate $m$, and the radius of cross-immunity $d$. A long-lived but transient RQS regime is flanked be the regime of deterministic extinction (red) and the regime of continuous branching and diversification – the ‘speciation’ regime (blue). The RQS regime itself undergoes a transition from a steady traveling wave (yellow) to a limit cycle oscillation (green) with increasing population size. The location of the boundaries depend on the time scale of observation as the cumulative probability of extinction and speciation increases with time.

Figure 3 with 2 supplements see all

Download asset Open asset

A novel virus in a completely susceptible population will initially spread with rate $\beta -1$ and the pandemic peaks when susceptible fraction falls to ${\beta }^{-1}$. The trajectory of such a pandemic strain in the time-susceptibility plane is indicated in red in Figure 3D. Further infections in the contracting epidemic will then push susceptibility below ${\beta }^{-1}$ – the propagation threshold for the virus – and without rapid antigenic evolution the pathogen will go extinct after a time $t\sim {\beta }^{-1}\log N_h$. Such boom-bust epidemics are reminiscent of the recent Zika virus outbreak in French Polynesia and the Americas where in a short time a large fraction of the population was infected and developed protective immunity (O'Reilly et al., 2018).

Persistence and transition to an endemic state is only possible if the pathogen can evade the rapid build-up of immunity via a small number of large effect antigenic mutations. This process is indicated in Figure 3D by horizontal arrows leading to antigenically evolved strains of higher susceptibility and bears similarity to the concept of ‘evolutionary rescue’ in population genetics (Gomulkiewicz, 1995). The parameter range of the idealized SIR model that avoid extinction after a pandemic resulting in persistent endemic disease is relatively small. Yet, various factors like geographic structure, heterogeneity of host adaptation and population turn-over slow down the pandemic and extinction, thereby increasing the chances of sufficient antigenic evolution to enter the endemic, RQS-type, regime. The 2009 pandemic influenza A/H1N1 has undergone such a transition from a pandemic to a seasonal/endemic state. We shall not investigate the transition process in detail here, but will assume that endemic regime has been reached.

Once the pathogen population has established an endemic circulation through continuous antigenic evolution (green and yellow regimes in Figure 3BC), the average rate of new infections $\beta {\sum }_a{I}_a{S}_a/I_{tot}$ fluctuates around the rate of recovery $\nu =1$ (in our time units). This balance is maintained by the steady decrease in susceptibility due to rising immunity against resident strains and the emergence of antigenically novel strains, see Figure 3D. If the typical mutational distance between strains is small compared to the cross-immunity range $d$, the rate at which susceptibility decreases is similar for all strains. To see this we expand $K_{ab}$ in Equation 5

where we have used that $|a-b|\ll d$ for all pairs of strains with substantial prevalence. In fact it will suffice to keep only the first, leading, term on the right hand side. Close to a steady state, prevalent strains obey $\beta S_a\approx 1$. We can hence define the instantaneous growth rate of strain $x_a=(\beta S_a-1)\ll 1$ as its effective fitness. In this limit, the model can be simplified to

The second equation means that effective fitness of all strains $a$ decreases approximately at the same rate since the pathogen population is dominated by antigenically similar strains.

If a new strain $c$ emerged from strain $a$ by a single antigenic mutation, its mutational distance from a strain $b$ is $|c-b|=|a-b|+1$ and $K_{cb}=K_{ab}{e}^{-d^{-1}}\approx K_{ab}(1-d^{-1})$. The population susceptibility of strain $c$ is therefore increased to

Since the typical susceptibility is of order ${\beta }^{-1}$, the growth rate of the mutant strain $c$ is $s=d^{-1}\log \beta$ higher than that of its parent. The growth rate increment, s, plays the role of a selection coefficient in typical population genetic models and corresponds to the step size of the fitness distribution in Figure 3D. In such models, individuals within a fitness class (bin of the histogram) are equivalent and different classes can be modeled as homogeneous populations which greatly accelerates numerical analysis of the model, see Materials and methods.

Rouzine and Rozhnova (2018) have recently formulated a similar model of antigenic evolution of rapidly adapting pathogens. Analogously to our model, Rouzine and Rozhnova couple strain dynamics to antigenic adaptation through mutations, albeit assuming a one-dimensional antigenic space. In agreement with Rouzine and Rozhnova, we find that selection coefficients of novel mutations are inversely proportional to the cross-immunity rate $d$ and increase with infectivity $\beta$, see Equation 9. Rouzine and Rozhnova, however, do not consider oscillations, extinction, and speciation (see below).

The simplified model in Equation 8, along with the model developed by Rouzine and Rozhnova (2018), is analogous to the traveling wave (TW) models of rapidly adapting asexual populations that have been studied extensively over the past two decades (Tsimring et al., 1996; Desai and Fisher, 2007; Rouzine et al., 2003; Hallatschek, 2011), see Neher (2013a) for a review. These models describe large populations that generate beneficial mutations rapidly enough that many strains co-circulate and compete against each other. The fittest (most antigenically advanced) strains are often multiple mutational steps, $q$, ahead of the most common strains, see Figure 3D. This ‘nose’ of the fitness distributions contains the strains that dominate in the future and the only adaptive mutations that fixate in the population arise in pioneer strains in the nose. Consequently, the rate with which antigenic mutations establish in the population is controlled by the rate at which they arise in the nose (Desai and Fisher, 2007). If the growth rate at the nose of the distribution, $x_n$, is much higher than antigenic mutation rate, $x_n\gg m$, it takes typically

generations before a novel antigenic mutation arises in a newly arisen pioneer strain that grows exponentially with rate $x_n$. The advancement of the nose is balanced rapidly by the increasing population mean fitness.

If beneficial mutations have comparable effects on fitness and population sizes are sufficiently large ($Nm\gg 1$), the fitness distribution has an approximately Gaussian shape with a variance ${\sigma }^2\approx 2s^2\log (Ns)/{\log }^2(x_n/m)$. The wave is $\sigma /s$ mutations wide, while the most advanced strains are approximately $q=2\log (Ns)/\log (x_n/m)$ ahead of the mean (Desai and Fisher, 2007). Two contemporaneous lineages coalesce on a time scale ${\tau }_{\mathrm{sw}}=sq/{\sigma }^2=s^{-1}\log (x_n/m)$ and the branching patterns of the tree resemble a Bolthausen-Sznitman coalescent rather than a Kingman coalescent (Desai et al., 2013; Neher and Hallatschek, 2013b).

In circulating influenza viruses, typically around 3–10 adaptive mutations separate pioneer strains from the most common variants (Strelkowa and Lässig, 2012; Neher and Bedford, 2015). While this clearly corresponds to a regime where multiple stains compete, it does not necessarily mean that asymptotic formulae assuming $q\gg 1$ are accurate. Nevertheless, many qualitative features of TW models have been shown to qualitatively extend into regimes where $q$ takes intermediate values (Neher and Hallatschek, 2013b).

While parameter $N$ in the TW models summarized above is a fixed population size, the corresponding entity in our SIR model is the fluctuating pathogen population size $N_p$ which is related to the (fixed) host population size $N_h$ by $N_p=N_h{I}_{tot}$. The average $I_{tot}$ depends on other parameters of the model, scaling in particular with $\overline{I}\sim s^2$. Hence, it will be convenient for us to use $N_h{s}^2$ as one of the relevant ‘control parameters’, replacing $N$ of the standard TW model.

In contrast to most population genetic models of rapid adaptation, our epidemiological model does not control the total population size directly. Instead, the pathogen population size (or prevalence) depends on the host susceptibility, which in itself is determined by recent antigenic evolution of the pathogen. The coupling of these two different effects results in a rich and complicated dynamics (see Figure 4A for an example trajectory): The first effect is ecological: a bloom of the pathogen depletes susceptible hosts leading to a crash in pathogen population and a tendency of the population size to oscillate London and Yorke, 1973 (blue line in Figure 4A). The second effect is evolutionary: higher nose fitness $x_n$ begets faster antigenic evolution and vice versa, resulting in an apparent instability in the advancement of the antigenic pioneer strains (Fisher, 2013) (yellow line and inset in Figure 4A). In our epidemiological model, as we shall show below, fluctuations in the rate of antigenic advance of the pioneer strains couple with a delay of ${\tau }_{\mathrm{sw}}$ to the ecological oscillation.

Figure 4 with 2 supplements see all

Download asset Open asset

To recognize the ecological aspect of the oscillatory tendency, consider the total prevalence $I_{tot}$ and the mean fitness of the pathogen $X={\sum }_a{x}_a{I}_a/I_{tot}$ 

which follows directly from Equation (8). Selection on fitness variance ${\sigma }^2$ increases $X$, while prevalence $I_{tot}$ reduces susceptibility and hence $X$. At fixed variance $\sigma =\overline{\sigma }$ this system is equivalent to a non-linear oscillator describing a family of limit cycles oscillating about $I_{tot}={\overline{\sigma }}^2$ and $X=0$ as shown in Figure 4B.

While Equation (11) describes the behavior of common strains, the mutation driven dynamics of the antigenic pioneer strains is governed by the equation for $x_n$ that in a continuum limit (suitable for the limit of high mutation rate) reads:

The first term on the right hand side represents the rate at which antigenic pioneer strains enter the population, ${\tau }_a^{-1}$, advancing the nose fitness by an increment $s$ (with ${\tau }_a^{-1}s={\tau }_{\mathrm{sw}}^{-1}{x}_n$ ). The second term on the right hand side of Equation (12) represents gradual reduction of susceptibility of the host population, and $\xi (t)$ is a random noise variable representing the stochasticity of the establishment of new strains. The Gaussian white noise $\xi (t)$ is defined statistically by its correlation function $⟨\xi (t)\xi (0)⟩={\tau }_a^{-1}\delta (t)$, see Materials and methods (Stochastic differential-delay simulation).

The first term of Equation (12) captures the apparent instability of the nose: an advance of the nose to higher $x_n$ accelerates its rate of advancement. The stabilizing factor is the subsequent increase in $I_{tot}$, but to see how that comes about we must connect Equation (12) to Equation (11). The connection is provided by ${\sigma }^2$ since it is controlled by the emergence of novel strains, that is the dynamics of the ‘nose’ $x_n$, which impacts the bulk of the distribution after a delay ${\tau }_{\mathrm{sw}}$. Based on the analysis detailed in the Appendix 2, we approximate

relating population dynamics, Equation (11), to antigenic evolution of pioneer strains described by Equation (12). Taken together Equations (11-13) define a Differential Delay (DD) system of equations. Sample simulations of this stochastic DD system are shown in Figure 4 BC. The delay approximation Equation (13) is supported by the cross-correlation of $x_n(t)$ and ${\sigma }^2(t^{\prime })$ measured using fitness-class simulations (see Figure 4A Inset).

The deterministic limit of the DD system (obtained by omitting the noise term in Equation (12)) has a fixed point at ${\tau }_{\mathrm{sw}}^{-1}{\overline{x}}_n={\overline{\sigma }}^2=2{\tau }_{\mathrm{sw}}^{-2}\log (N_h\overline{I})$. Small deviations decay in underdamped oscillations with frequency $\omega =\overline{\sigma }=\tau _{\mathrm{sw}}{}^{-1}\sqrt{2\log (N_h\overline{I})}$ if $\omega {\tau }_{\mathrm{sw}}<2\pi$. For $\omega {\tau }_{\mathrm{sw}}>2\pi$, the system fails to recover from a deviation of the nose in a single period and the steady state becomes unstable to a limit cycle oscillation. The nonlinearity of Equation (11) implies a longer period with increasing amplitude and the system is stabilized at a limit cycle with the period long enough compared to the feedback delay ${\tau }_{\mathrm{sw}}$. In Appendix 3, we derive the threshold of oscillatory instability to lie at $\log (N_h{\overline{I}}_{osc}s)\approx 8.3$ (leading to limit cycle period $T\approx 1.5{\tau }_{\mathrm{sw}}$, see Figure 4—figure supplement 1). We also find that the amplitude of the oscillation $\log (I_{max}/\overline{I})$ scales as $\log (N_h\overline{I})$ for large values of the later. This transition defines quantitatively the boundary between the TW RQS and the Oscillatory RQS regimes that appear on the phase diagrams in Figure 3 (BC). The validity of the predictions of standard TW theory for our adapting SIR system are explored in Figure 4—figure supplement 2.

The distinction between the TW and Oscillatory RQS is obscured by the stochasticity of antigenic advance, Equation (12), which continuously feeds the underdamped relaxation mode, generating a noisy oscillation with the frequency $\omega$ defined above. The difference between the two regimes is illustrated by Figure 4C: in the TW RQS noisy oscillation is about the fixed point, whereas in the Oscillatory RQS it is about deterministic limit cycle.

Interestingly, the dynamics of the Oscillatory RQS, as shown in Figure 4A, can be understood in terms of a non-linear relaxation oscillator. At relatively low infection prevalence nose fitness $x_n$ increases until rising $I_{tot}$ catches up with it (when $I_{tot}={\tau }_{\mathrm{sw}}^{-1}{x}_n$) driving it down rapidly. Once this ‘mini-pandemic’ burns out, the population returns to the low prevalence part of the cycle $I_{tot}<{\tau }_{\mathrm{sw}}^{-1}{x}_n$, when $x_n$ begins to increase again.

While in the deterministic limit the differential-delay system predicts a stable steady TW for $q>q_{ex},\overline{I}<{\overline{I}}_{osc}$ and a limit cycle above ${\overline{I}}_{osc}$, fluctuations in the establishment of the antigenic pioneer strains (Equation (12)) can lead to stochastic extinction. In fact, both the TW and Oscillatory RQS (see Figure 3BC) are transient, subject to extinction due to a sufficiently large stochastic fluctuation. (Note however the contrast with the ‘extinction’ state in Figure 3BC, where extinction is deterministic and rapid.) The rate of extinction depends on $q$ and $\log (N_h\overline{I})$ as shown in Figure 5A. The time to extinction increases dramatically in the range of $q\sim 1-2$ and more slowly thereafter. Although extinction is fluctuation driven, the mechanism of extinction in the oscillatory state is related closely to the deterministic dynamics, according to which large amplitude excursion in infection prevalence can lead to extinction. A large $x_n$ advance leads, after a time ${\tau }_{\mathrm{sw}}$ to a rise in prevalence $I_{tot}$, followed by the rapid fall in the number of susceptible hosts and hence loss of viral fitness. This turns out to be the main mode of fluctuation driven extinction as illustrated by Figure 4C. One expects extinction to take place when a fluctuation induced deviation $\delta x$ of the fitness of pioneer strains becomes of the order of the mean ${\overline{x}}_n$. New mutations at the nose accumulate with rate $1/{\tau }_a$ such that at short times $t$ we expect $\delta x\approx s\sqrt{t/{\tau }_a}$. Hence $\delta x$ becomes of the order of the mean ${\overline{x}}_n$ at times ${\tau }_{\mathrm{ext}}\sim q{\tau }_{\mathrm{sw}}$. However the probability of extinction will also depend on the shape of the oscillatory limit cycle (as it depends on the minimum of infection prevalence during the cycle), which in turn depends on $\log (N_h\overline{I})$. Numerical simulations, Figure 5B, confirm the dependence of ${\tau }_{\mathrm{ext}}$ on $q$ and $\log (N_h\overline{I})$. We note that the rate increase in ${\tau }_{\mathrm{ext}}$ with increasing $q$ slows down in the oscillatory regime and appears to approach a power law dependence ${\tau }_{\mathrm{ext}}/{\tau }_{\mathrm{sw}}\sim q^{2.5}$ (albeit over a limited accessible range): presently we do not have an analytic understanding of this specific functional form.

Figure 5

Download asset Open asset

The correspondence of the multi-strain SIR and the TW models discussed above assumes that cross-immunity decays slowly compared to the coalescent time of the population, that is $d/q\gg 1$. In this case, all members of the population compete against each other for the same susceptible hosts. Conversely, if the viral population were to split into two sub-populations separated by antigenic distance greater than the range of cross-inhibition $d$, these sub-population would no-longer compete for the hosts, becoming effectively distinct viral ‘species’ that propagate (or fail) independently of each other. Such a split has for example occurred among influenza B viruses, see Figure 1.

A ‘speciation’ event corresponds to a deep split in the viral phylogeny, with the $T_{MRCA}$ growing without bounds, see Figure 1 and Figure 6A. This situation contrasts the phylogeny of the single competing population, where $T_{MRCA}$ fluctuates with a characteristic ramp-like structure generated by stochastic extinction of one of the two oldest clades. In each such extinction event the MRCA jumps forward by $\delta T$. Hence the probability of speciation depends on the probability of the two oldest clades to persist without extinction for a time long enough to accumulate antigenic divergence in excess of $d$. The combined carrying capacity of the resulting independent lineages is then twice their original carrying capacity as observed in simulations, see Figure 6B.

Figure 6

Download asset Open asset

To gain better intuition into this process let’s follow two most antigenically advanced ‘pioneer strains’. In the TW approximation one of these will with high probability belong to the backbone giving the rise to the persisting clade, while the other clade will become extinct, unless it persist long enough to diverge antigenically beyond $d$, becoming a speciation event. As their antigenic distance gradually increases, the two clades are evolving to evade immunity built up against the common ancestor. The less advanced of the two clades is growing less rapidly and takes longer to generate antigenic advance mutations, resulting in still slower growth and slower antigenic advance. Deep splits are hence unstable and it is rare for a split to persist long enough for speciation. In Appendix 5, we reformulate this intuition mathematically as a ‘first passage’-type problem which shows that $T_{MRCA}$ distribution has an exponential tail which governs the probability of speciation events. Figure 6C shows that the time to speciation increases approximately exponentially with the ratio $d/q$. More precisely we found that average simulated speciation time behaves as ${\tau }_{\mathrm{sw}}^{*}{e}^{f(CI/q^{*})}$ with ‘effective’ ${\tau }_{\mathrm{sw}}^{*}={\tau }_{\mathrm{sw}}/(1+\log q/\log (s/m))$ and $q^{*}=q(1+\log q/\log (s/m))$ picking up an additional logarithmic dependence on parameters, the exact origin of which is beyond our current approximations. This correction plausibly suggests rapid speciation, ${\tau }_{\mathrm{sw}}^{*}\to 0$, when mutation rate become comparable to the selection strength $m/s\to 1$.

We emphasize that the RQS regime in Figure 3BC is only transient. For any given $q$ and $d$, the RQS is likely to persist for a time given by the smaller of ${\tau }_{\mathrm{ext}}$ and ${\tau }_{\mathrm{sp}}$, before undergoing either extinction or speciation. These two processes limit the range of $q$ corresponding to the RQS from both sides in a time-dependent manner. Figure 7 shows the likely state of an RQS system after time $\tau$ as a function of genetic diversity $q$ for the case of $d=50$ and $\log (N_h\overline{I})=6.5$.

Figure 7

Download asset Open asset

The regime of a single persistent lineage shrinks with increasing $\tau$, for example after $\tau =10{\tau }_{\mathrm{sw}}$ the RQS state likely prevails between $q=1.5$ and $\approx 4$, while (for $d=50$ and $\log (N_h\overline{I})=6.5$) it is unlikely to persist beyond $\tau \approx 100{\tau }_{\mathrm{sw}}$ for any $q$. Both the maximal RQS lifetime and corresponding critical ${\displaystyle q_c}$, increase with increasing $d$.

We simulate the original model on a genealogical tree by combining the deterministic update of SIR-type equations and the stochastic step introducing mutated strains. In each time step $\mathrm{\Delta }t<1$, we apply the mid-point method to advance SIR equations Equations (1,2,4). We then generate a random number uniformly sampled between zero and one for each surviving strain with $N_h{I}_a\ge 1$. If the random number is smaller than $m{N}_h{I}_a\mathrm{\Delta }t$ for strain $a$, we append a new strain $b$ as a descendent to $a$. The susceptibility to strain $b$ is related to susceptibility to strain $a$ via $S_b={(S_a)}^{e^{-1/d}}$. In most of the simulations, the transmissibility of different strains is held constant $\beta$. Otherwise we allow for a strain specific transmissibility that is to its parent: ${\beta }_b={\beta }_a-\delta \beta$ with $\delta \beta >0$ for the deleterious effect of antigenic mutations and ${\beta }_b={\beta }_{\max }$ if the mutation is compensatory. The new strain grows deterministically only if ${\beta }_b{S}_b>1$.

This simplified model contains six relevant parameters: transmissibility $\beta$, recovery rate $\nu$, mutation rate of the virus $m$, birth/death rate of the hosts $\gamma$, the effective cross-immunity range $d$, and the effective size of the hosts $N_h$, whose empirical ranges are summarized in the Table 1. For flu and other asexual systems in RQS, $\beta \gtrsim \nu \gg m,\gamma$, $d\gg 1$, and $N_h\gg 1$.

Simulation code and output are available on github in repository FluSpeciation of the neherlab organization (Neher and Yan, 2019; copy archived at https://github.com/elifesciences-publications/FluSpeciation).

The stability of the RQS and the extinction dynamics is fully captured by the traveling wave Equation (8). We simulate the traveling wave by discretizing the fitness space $x$ into bins of step size $s$ around zero. The number of individuals infected by different strains correspond to integers in each bin $x_i$. At each time step, the population in each bin $I_i$ updates to a number sampled from the Poisson distribution with parameter ${\lambda }_i=N_h{I}_i(1+(x_i-\overline{x})\mathrm{\Delta }t)$ determined by mean fitness $x_i$ and a dynamic mean fitness $\overline{x}$, which increases by $\mathrm{\Delta }t{I}_{tot}$, where $I_{tot}$ is the total infected fraction summed over all bins. When $\overline{x}$ becomes larger than one bin size $s$, we shift the all populations to left by one bin and reset $\overline{x}$ to , a trick to keep only a finite number of bins in the simulation. At the same time, antigenic mutation is represented by moving the mutated fraction in each bin to the adjacent bin on the right. The fraction is determined by a random number drawn from the Poisson distribution with the mean $m{I}_i\mathrm{\Delta }t$. The typical ranges of the three parameters $s$, $m$, and $N_h$ follow the parameters in the genealogical simulation, as documented also in Table 1.

To simulate the differential delay equations Equations (11-13), we discretize time in increments of $\mathrm{\Delta }t={\tau }_{\mathrm{sw}}/k$ and update the dynamical variables ${\chi }_i=x_n(t_i)$ and ${\eta }_i=I_{tot}(t_i)$ via the simple Euler scheme:

where ${\xi }_i$ is a Gaussian random variable with zero mean and unit variance. Mean prevalence, $\overline{I}$, enters as the control parameter (which defines the time average of ${\eta }_i$).

Influenza virus HA sequences for the subtypes A/H3N2, A/H1N1, A/H1N1pdm, as well as influenza B lineages Victoria and Yamagata were downloaded from fludb.org.

We aligned HA sequences using mafft (Katoh et al., 2002) and reconstructed phylogenies with IQ-Tree (Nguyen et al., 2015). Phylogenies were further processed and time-scaled with the augur (Hadfield et al., 2018) and TreeTime (Sagulenko et al., 2018). The analysis pipeline and scripts are available on github in repository 2019_Yan_flu_analysis of the neherlab organization.

A microscopic model that tracks the infection history of every individual in population is computationally costly and impossible to analyze analytically. To gain insight, we and other authors before us have used approximations that reduce the exploding combinatorial complexity of the state space (Kryazhimskiy et al., 2007). Here, we explore and justify the two separate approximations we have made to arrive at Equation 2: We ignore correlations between subsequent infections of the same individual and approximate the multiplicative effect of all subsequent infections by an exponential term.

To derive Equation 4 we start with Equation 3 and expand it in powers of $K$ 

where angular brackets denote the average over all individuals $i$ in the population and ${\sigma }_b(i)\in [0,1]$ denotes whether individuals $i$ was infected with strain $b$ in the past. This expansion assumes $|K_{ab}|\ll 1$ which would hold uniformly for weak inhibition $\alpha \ll 1$ but also holds for perfect inhibition for sufficiently distant strains $a,b$. For $\alpha \approx 1$ the greatest cause of concern is the contribution of the most proximal strain, to which we shall return later. To evaluate the terms on the right-hand-side we note that $⟨{\sigma }_b(i)⟩=R_b$, that is the fraction of the population recovered from $b$, and $⟨{\sigma }_b(i){\sigma }_c(i)⟩=R_b{R}_c+{\rho }_{bc}$ where ${\rho }_{bc}$, by definition, is the correlation between infection with $b$ and $c$ at the level of individuals. Our approximation – following the well established logic of ‘Mean Field’ theories – neglects ${\rho }_{bc}$ compared to $R_b{R}_c$ (Landau and Lifshitz, 2013; Weiss, 1907). In this case, correct to order ${|K|}^2$, we can re-exponentiate the right-hand-side obtaining Equation 4. This simple derivation effectively captures the content of the ‘order-1 independence closure’ in Kryazhimskiy et al. (2007).

Several facts about influenza in human populations suggest that the weak-correlation approximation is a reasonable starting point for modeling population scale behavior. (i) Seasonal flu epidemics involve a large number of strains, a particular strain infects only a small fraction of the population. Hence the $R_a$ are small and correlation effects are of minor importance. (ii) Challenge studies have shown that protection through vaccination or infection with antigenically similar strains is moderate and a large fraction of challenged individuals still shed virus (Clements et al., 1991). This possibility of homotypic re-infection shows that all $K_{ab}$ are substantially smaller than 1, supporting our approximation of population wide susceptibility, as discussed above. (iii) Antibody responses are polyclonal and differ between individuals such that the cross-immunity matrix is stochastic at the level of individuals. This variation in the cross-immunity matrix further reduces correlations in infection history at the population level and justifies the mean field approach taken here (Lee et al., 2019). (iv) Correlation in infection history induced by immunity are further reduced by the variation in exposure history through geography and variation in contact networks.

To quantify the error made by these approximations in the worst case scenario, we explore the case of a one-dimensional strain space with strictly periodic re-infection as soon as the virus population as evolved by $ϵd$ – the case of maximal correlation. The susceptibility of an individual last infected with a strain $x<ϵd$ mutations away from the current strain has susceptibility

where we separated the most recent infection from previous infections to explicitly compare our approximations to that of Rouzine and Rozhnova (2018). The only approximation so far happened between step 2 and 3. In the following, we will include the most recent infection in sum in the exponential to obtain the weak-inhibition approximation $\log S_{\mathrm{WI}}=-\alpha e^{-\frac{x}{d(1-e^{-ϵ})}}$.

Rouzine and Rozhnova (2018) follow Lin et al. (2003) in approximating an individual’s susceptibility by ignoring all but the most recent infection $S\approx 1-K(x)$, thus keeping only the smallest term of the product representation of $S_a(i)$ in Equation 3. This approximation is referred to as ‘minimum’ cross-immunity in Wikramaratna et al. (2015). Appendix 1—figure 1 compares the full expression and different approximations of $S(x)$ for different values of $\alpha$ and $ϵ=1$. For $\alpha =1$, the ‘most recent’ approximation is better than the ‘weak inhibition’ approximation for $x<d/2$ but worse otherwise. For $\alpha <1$, the weak inhibition approximation improves further.

Appendix 1—figure 1

Download asset Open asset

To investigate the effect of ignoring correlations, we now compare the most correlated case of strictly periodic re-infection as soon as the pathogen has evolved by $ϵd$. For simplicity, we assume a time invariant density of recovered $1/dϵ$ (as in the analysis by Rouzine and Rozhnova, 2018). To calculate the population susceptibility, we integrate the expression for $S(x)$ for the full model and the ‘most recent’ approximation over the interval $[0,dϵ]$ and compare it to the mean field approximation in Equation 4 with constant $R_b=1/dϵ$. Appendix 1—figure 1B shows that the ‘mean-field’ approximation is closer to the full model across the entire range of relevant $ϵ<1$. Note that $ϵ$ has to be determined self-consistently and will typically be of the order of the susceptibility.

Real-world influenza population are much less correlated then the extreme ‘periodic infection’ assumption used here for reasons listed above. The linearized mean-field approximation in Equation 4 is therefore justified and can be expected to give a qualitatively correct approximation to a full model that tracks all infection histories.

Here we derive the differential delay system of equations that relate the behavior of the pioneer strains to the bulk of the population. Let us consider the generating function associated with the virus fitness distribution at time $t$:

where $x_i(t)=x_n(t_i)-{\int }_{t_i}^t𝑑t^{\prime }{I}_{tot}(t^{\prime })$ is the current fitness of the pioneer strain that first appeared at time $t_i$ and $I_i(t)$ is the fraction of the hosts infected by it:

We next take a coarse grained view of pioneer strain establishment replacing the sum in Equation (13) by an integral over initial times $t_i\to t-\tau$ 

where $1/{\tau }_a(t-\tau )$ is the rate at which new clones are seeded at time $t-\tau$. Let us evaluate the integral in the saddle approximation which is dominated by $\tau ={\tau }^{*}$ corresponding to the maximum in the exponential

where we have used the deterministic limit of Equation 12. To simplify presentation we shall ignore the time dependence of ${\tau }_{\mathrm{sw}}=s^{-1}\log (x_n/m)$ replacing $x_n(t-{\tau }^{*})$ in the logarithm by the time average ${\overline{x}}_n$.

Within the saddle approximation we then have

where we have omitted the logarithmic corrections for simplicity. Note that by definition $G(0,t)=I_{tot}(t)$.

We can now estimate fitness mean

and variance

Equation (A2.7) involves the second derivative $x_n^{\prime \prime }$ and we therefore expect fluctuations in the establishment of new lineages (which contribute to $x_n^{\prime }$) to be quite important. Yet we can get useful insight by using the deterministic approximation to $x_n$ dynamics in Equation 12, in which case we arrive at simple delay relation between the variance and $x_n$ 

which is consistent with the variance calculated for the case of the steady TW and also satisfies the generalized Fisher theorem

Combining Equations 11, 12 and A2.8 we arrive at the deterministic dynamical system approximating coupled ‘ecological’ SIR dynamics with the evolutionary dynamics of antigenic innovation due to the pioneer strains.

This system admits a family of fixed points of the form ${\tau }_{\mathrm{sw}}{I}_{tot}=x_n={\overline{x}}_n$, but as we show in C, the corresponding steady TW states are not always stable giving rise to limit cycle oscillations or leading to rapid extinction. The self-consistency condition relating $x_n$ and $I_{tot}$ for the steady traveling wave is readily generalized to limit cycle states. Integrating the differential-delay system over one cycle yields $⟨x_n⟩={\tau }_{\mathrm{sw}}⟨I⟩$. An additional relation is provided by integrating $\log N_{h}G(0,t)$ over the cycle:

A great deal of insight into the behavior of the (deterministic) differential delay system defined above is provided by its deterministic limit (see Appendix 3) which defines the stability ‘phase diagram’ shown in Figure 3 (BC) that correctly captures key aspects of the behavior observed in fully stochastic simulations.

In the traveling wave case, it is natural to measure time in the units of the delay time scale ${\tau }_{\mathrm{sw}}$. The therefore define a time variable $\zeta$ via $t={\tau }_{\mathrm{sw}}\zeta$, the fitness variable $\chi$ via $x_n={\tau }_{\mathrm{sw}}^{-1}\chi$ and the rescaled log-prevalence $u$ via $u=\log {\tau }_{\mathrm{sw}}^{2}I$ to obtain

As before, this system has a one parameter family of fixed points $\chi =\overline{\chi },u=\log \overline{\chi }$. Note that from the traveling wave model (Desai and Fisher, 2007), we have $\overline{\chi }={\overline{x}}_n{\tau }_{\mathrm{sw}}=q\log (x_n/m)=2\log (N_h{s}^2)$. To analyze fixed point stability we linearize and Laplace transform, yielding

Stability is governed by the poles of the Laplace transformed response to the initial perturbation $\delta u(0),\delta u^{\prime }(0),\delta \chi (0)$ and these poles are at the complex $z$ that solve:

Fixed point - and hence steady RQS - stability requires $\mathrm{\Re }(z)<0$ which is found for $2<\overline{\chi }<{\overline{\chi }}_c$. For $\overline{\chi }>2.845$ one finds $\mathrm{\Im }(z)\ne 0$ corresponding to the onset of oscillatory relaxation which turns into a limit cycle for $\overline{\chi }>{\overline{\chi }}_c\approx 16.6$. The period of the limit cycle is well approximated by $\mathrm{\Im }(z)$, as the dashed line shown in the bottom panel of Figure 4—figure supplement 1.

The above stability analysis is done for the continuum limit $q\gg 1$. However the finiteness of $q$ does matter, especially close to extinction where only a small number of mutations separate most advanced strains from the bulk of the distribution. We shall now include the corrections to the first order in $1/q$. One such correction arises from the difference between the continuum $\chi (t_i)$ and discrete approximation of ${\chi }_i$, the position of the nose fitness bin relative to mean fitness at the time of its establishment. The other correction term comes via the establishment time ${\tau }_a$. Including both corrections the Langevin equation becomes:

To first order of $1/q$, the poles of the Laplace transform are determined by

Solving for the onset of stability $\mathrm{\Re }(z)=0$, we find the extinction boundary $q_{ex}(\log N_h{s}^2)$ from the relation

We observe that $q_{ex}\to 2$ asymptotically for large $\log N_h{s}^2$.

A sensible stochastic generalization is obtained by the stochastic approximation for the ‘nose’ dynamics in Equation (12)

combined with Equation (A2.5) at $\lambda =0$

Note that in this derivation we have avoided the need for explicitly approximating ${\sigma }^2$! (We have also neglected the effect of fluctuations arising from the logarithmic correction term effectively replacing it by its average value.) This stochastic differential delay (DD) system was used in simulations presented in Figure 4C.

Speciation occurs when two most distant clades persist until they are antigenically independent. This persistence problem can be formulated as a first passage problem by including the second ‘nose’ in the TW approximation.

We consider the birth of two pioneer strains at time $t=0$, as illustrated in Appendix 5—figure 1. The descendants of the two strains forming two branches 1 and 2 diverge in the antigenic space as they persist in time. Suppose that at time $t$, the nose of branch 1 is at fitness $x_1$, and the nose of branch 2 is at $x_2$. Before the sweep time $t<{\tau }_{\mathrm{sw}}$, the cross-immunity grows mainly from the prevalent strains in the common ancestors of the two branches,

Appendix 5—figure 1

Download asset Open asset

Later when $t>{\tau }_{\mathrm{sw}}$, the pathogen population splits and the different lineages evolve away from each other on two branches in the phylogeny. As the antigenic distances of each nose from the dominant strains on the own and the other branch differ, fitness of the two sets of pioneer strains changes at different rates:

where $d_{11}$ and $d_{22}$ scale roughly as $q$, the typical antigenic distance to the nose. In the limit $d_{21}\approx d_{12}\gtrsim d$, Equation (A5.2) reduce to two independent replicas of Equation (12) and the two branches are thus antigenically independent. What is the probability of reaching this limit? The approach to this question rather relies on the persistence probability of two branches in the other limit when $d_{21}\approx d_{12}\lesssim d$, where $I_1+I_2\approx I_{\mathrm{tot}}$ cross-immunity growth rate is approximately the same at both noses.

In this limit, the survival probability of the less fit nose maps to a first passage problem in the random walk of relative fitness $\zeta \equiv (x_1-x_2)/x_n$. As illustrated in Appendix 5—figure 1, an establishment of nose one is a positive step of $\delta \zeta =s/x_n$, while an establishment of nose two results in a backward step of the same size. As the mutations arrive in characteristic times ${\tau }_1$ and ${\tau }_2$ depending on the nose fitnesses, in the continuum limit, we have

where $\xi$ is a random noise. There are two relevant boundaries: a reflecting boundary at $\zeta =0$ where two branches switch roles in leading the fitness, and an absorbing boundary at $\zeta =1$ where the fitness of less fit nose drops below the mean fitness and becomes destined for extinction.

The Langevin Equation in Equation A5.3 corresponds to a diffusion equation for the probability density distribution $\rho (\zeta ,t)$ 

where the drift $v$ and diffusivity $D$ depend on $\zeta$,

Solving with boundary and initial conditions,

we have

where ${}_{1}F_1$ is the generalized hypergeometric function, ${\lambda }_n$ is the $n$th smallest values solving ${}_{1}F_1(\frac{1-\lambda }{2},\frac{1}{2},\frac{q}{2})=0$, and coefficient $c_n$ is determined by the initial condition. In long time $t$, the slowest mode dominates the dynamics. In the large $q$ limit, we have ${\lambda }_1=1$. Since ${}_{1}F_1\approx \mathrm{const}$ for $\zeta \in (0,1)$, the persistence probability is

The typical time interval between the establishment of successive pioneer strains at the nose scales as ${\tau }_a={\tau }_{\mathrm{sw}}/q$. We recall that speciation, or escape from cross-immunity, occurs when antigenic distance between the two branches in Appendix 5—figure 1 $d_1+d_2$ is larger than $d$. For that it suffices that the shorter branch $d_2>d/2$ which occurs with probability

we find the probability of a successful branching $p_1$ to be proportional to $e^{-d/2q}$.

In the phylogenetic tree, $t/{\tau }_a$ trial branchings from the backbone arrive in time $t$. The probability that none of them successfully speciate is thus

where the waiting time for speciation event is

as numerically verified in Figure 6.

Suppose an antigenic mutation has a deleterious effect on infectivity reducing the latter by $\delta {\beta }_d$ on average. This would effectively reduce the fitness gain of antigenic innovation from $s$ to $s_d(\beta )=s(\beta )-{\mathrm{\Delta }}_d$, with ${\mathrm{\Delta }}_d={\beta }^{-1}\delta {\beta }_d$. In addition let us assume that there also are compensatory mutations which restore maximal infectivity ${\beta }_{\max }$. These compensatory mutation thus have a beneficial effect on fitness ${\mathrm{\Delta }}_b(\beta )={\beta }^{-1}{\beta }_{\max }-1$. We assume that these mutations occur with rate $m_{\beta b}$. In a dynamic balance state the rate of fixation of compensatory mutations would exactly balance the deleterious mutation effect on $\beta$ so that ${\tau }_b^{-1}{\mathrm{\Delta }}_d={\tau }_a^{-1}{\mathrm{\Delta }}_d$ with the fixation rate controlled by the fitness of the leading strain via ${\tau }_b^{-1}=x_n/\log (\frac{x_n}{m_{\beta b}})$. This dynamic balance is achieved at a certain value of ${\beta }_{*}<{\beta }_{\max }$, specifically ${\beta }_{\max }-{\beta }_{*}=\delta {\beta }_d{\tau }_b{\tau }_a^{-1}$ or ${\beta }_{*}={\beta }_{\max }-\delta {\beta }_{d}r$ where $r=\log (\frac{x_n}{m})/\log (\frac{x_n}{m_{\beta b}})$.

The fitness of the nose of the distribution obeys

where the 1 st term on the RHS is rate of nose advancement due to antigenetic mutations ${\tau }_a^{-1}=x_n/\log (\frac{x_n}{m})$ as before, but with reduced fitness gain $s_d(\beta )$. The 2nd term describes the contribution of compensatory mutations. However in the dynamic equilibrium (at ${\beta }_{*}$) compensatory mutations exactly cancel the contribution the deleterious mutation contribution to $s$ so that for the steady state we recover

as we had for the TW driven by antigenic advancement only. The only effect is the reduction of $s$ from $s({\beta }_{\max })$ to $s({\beta }_{*})=d^{-1}\log {\beta }_{*}$.

The sweep time, ${\tau }_{\mathrm{sw}}$, upon which the fitness of the former pioneer strain comes down to the mean fitness and the nose fitness, $x_n$, retain the TW form

Following TW approximation to estimate infection prevalence $\sqrt{I_{tot}}\sim N_h^{-1}\exp (x_n{\tau }_{\mathrm{sw}}/2)$ as before one finds

The total fitness variance of the population contains a contribution, from antigenic mutations and the mutations in infectivity:

but under conditions of ${\mathrm{\Delta }}_d,{\mathrm{\Delta }}_b\ll s({\beta }_{*})$ total variance would also be decreasing.

Most relevant for our analysis however is not the typical, but the maximal antigenic distance within the viral population: