1. Physics of Living Systems
Download icon

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

  1. Le Yan  Is a corresponding author
  2. Richard A Neher  Is a corresponding author
  3. Boris I Shraiman  Is a corresponding author
  1. University of California, Santa Barbara, United States
  2. University of Basel, Swiss Institute for Bioinformatics, Switzerland
Research Article
  • Cited 0
  • Views 720
  • Annotations
Cite this article as: eLife 2019;8:e44205 doi: 10.7554/eLife.44205

Abstract

Rapidly evolving pathogens like influenza viruses can persist by changing their antigenic properties fast enough to evade the adaptive immunity, yet they rarely split into diverging lineages. By mapping the multi-strain Susceptible-Infected-Recovered model onto the traveling wave model of adapting populations, we demonstrate that persistence of a rapidly evolving, Red-Queen-like state of the pathogen population requires long-ranged cross-immunity and sufficiently large population sizes. This state is unstable and the population goes extinct or ‘speciates’ into two pathogen strains with antigenic divergence beyond the range of cross-inhibition. However, in a certain range of evolutionary parameters, a single cross-inhibiting population can exist for times long compared to the time to the most recent common ancestor (TMRCA) and gives rise to phylogenetic patterns typical of influenza virus. We demonstrate that the rate of speciation is related to fluctuations of TMRCA and construct a ‘phase diagram’ identifying different phylodynamic regimes as a function of evolutionary parameters.

https://doi.org/10.7554/eLife.44205.001

Introduction

In a host population that develops long-lasting immunity, a pathogen can persist by infecting immunological naive individuals such as children, or through rapid antigenic evolution that enables the pathogen to evade immunity and re-infect individuals. Childhood diseases like measles or chicken pox fall into the former category, while influenza virus adapts rapidly and re-infects most humans multiple times during their lifespan. The continuous adaptation of influenza is facilitated by high mutation rates resulting in diverse populations of co-circulating viral strains. Nevertheless, almost always a single variant eventually outcompetes the others such that diversity within one subtype or lineage remains limited (Petrova and Russell, 2018).

The contrast of rapid evolution while maintaining limited genetic diversity is most pronounced for the influenza virus subtype A/H3N2. Figure 1 shows a phylogenetic tree of HA sequences of type A/H3N2 with the characteristic ‘spindly’ shape. The most recent common ancestor of the population is rarely more than 3–5 years in the past (Rambaut et al., 2008). Other pathogenic RNA viruses that typically do not reinfect the same individual, (measles, mumps, HCV, or HIV) diversify for decades or centuries (Grenfell et al., 2004). Interestingly, influenza B has split into two co-circulating lineages in the 1970s which by now are antigenically distinct (Rota et al., 1990) and maintain intermediate levels of diversity (see Figure 1).

Spindly phylogenies and speciation in different human seasonal influenza virus lineages.

The top left panel shows a phylogeny of the HA segment of influenza A virus of subtype H3N2 from its emergence in 1968 to 2018. The virus population never accumulates much diversity but is rapidly evolving. The lower left panel shows a phylogeny of the HA segment of influenza B viruses from 1940 to 2018. In the 70s, the population split into two lineages known as Victoria (B/Vic) and Yamagata (B/Yam) in the 1970s. The graphs on the right quantify diversity via the time to the most recent common ancestor TMRCA for different influenza virus lineages. Influenza B viruses harbor more genetic diversity than influenza A viruses. The subtype A/H3N2 in particular coalesces typically in 3y while deeps splits in excess of 5y are rare.

https://doi.org/10.7554/eLife.44205.002

Influenza virus infections elicit lasting immunity rendering most individuals non-susceptible to viruses that circulated during their lifetime (Fonville et al., 2014). The virus population escapes collective human immunity by accumulating amino acid substitutions in its surface glycoproteins (Koel et al., 2013; Wilson and Cox, 1990). Extensive genetic characterizations have shown that within each subtype many HA sequence variants co-circulate (Rambaut et al., 2008; Fitch et al., 1997). These variants differ from each other by ∼10 substitutions and compete for susceptible hosts (Strelkowa and Lässig, 2012). The rapid sequence evolution results in a decay of immune cross-reactivity over ∼10 years (Smith et al., 2004; Bedford et al., 2014; Fonville et al., 2014; Neher et al., 2016).

Epidemiological dynamics of influenza is often modeled using generalizations of the classic Susceptible-Infected-Recovered (SIR) model to multiple antigenically distinct viral strains (Kermack and McKendrick, 1927; Gog and Grenfell, 2002). Such models need to capture (i) how the infection with one strain affects susceptibility to other strains and (ii) how novel strains are generated from existing strains by mutations. A common approach has been to impose a discrete one-dimensional strain space in which new strains are generated by mutation of adjacent strains. Infection results in a reduction of susceptibility in a manner that depends on the distance in this one-dimensional strain space (Andreasen et al., 1996; Gog and Grenfell, 2002). Such models naturally result in ‘traveling waves’ in the sense that the pathogen population moves through strain space by recurrent emergence of antigenically advanced variants produced by mutation from neighboring strains (Lin et al., 2003).

These models of antigenically evolving populations are related to general models of rapid adaptation in which populations form a traveling wave moving towards higher fitness (Tsimring et al., 1996; Rouzine et al., 2003; Desai and Fisher, 2007; Neher et al., 2014), reviewed in Neher (2013a). Recently, Rouzine and Rozhnova (2018) described an explicit mapping between a SIR model in a one-dimensional antigenic space and traveling wave models in fitness.

Traveling wave (TW) models in a one-dimensional antigenic space naturally result in spindly phylogenies: There is only one possible direction for immune escape and the fastest growing most antigenically advanced strain grows drives all other strains extinct. Influenza viruses, however, can escape immunity by mutations at a large number of positions (Wilson and Cox, 1990), suggesting antigenic space is high dimensional (Perelson and Oster, 1979). In many dimensions, different viral strains can escape immunity via different paths and diverge sufficiently from each other until they no longer compete for hosts and thereafter propagate independently evolve. A satisfactory explanation of spindly phylogenies therefore has to describe how evolution in a high dimensional space reduces to an effectively one-dimensional path without persistent branching or rapid extinction. Several computational studies have addressed this question and identified cross-immunity (Bedford et al., 2012; Tria et al., 2005; Koelle et al., 2011; Ferguson et al., 2003; Sasaki and Haraguchi, 2000) as well as deleterious mutations (Koelle and Rasmussen, 2015; Gog and Grenfell, 2002) as critical parameters. We will discuss this earlier work at greater length below.

Our work aims to examine the conditions under which the evolving pathogen can maintain a spindly phylogeny with an approximately constant level of diversity – sufficient to avoid extinction, yet constrained from further branching by cross-inhibition between not too distant strains. We show that long range cross immunity in generic stochastic models of antigenic evolution generates such phylogenies. However, in the long term the viral population either ‘speciates’ into weakly interacting diverging lineages or goes extinct with rates that are controlled by three dimensionless combinations of model parameters. While the relation of these parameters to the known characteristics of influenza epidemiology and evolution is not direct, the general ‘phase diagram’ captured by the parameters of the simple model illustrates the key competing factors governing expected long-term dynamics.

Results

Model

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 Ia infected with viral strain a changes according to

(1) ddtIa=βSaIa-(ν+γ)Ia

where β is the transmissibility, Sa is the population averaged susceptible to strain a, ν is the recovery rate, and γ is the population turnover rate. The fraction Ra of the population recovered from infection with strain a changes according to

(2) ddtRa=νIaγRa

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 γ=0 right away to simplify presentation.

The dynamics of Ia depends on the average susceptibility of the host population Sa=Sa(i)i, while the susceptibility Sa(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)

(3) Sa(σ)=b(1-Kabσb(i))i

where σb(i) is one or zero depending on whether host i has or has not been previously infected with strain b. Matrix Kab1 quantifies the cross-immunity to strain a due to prior infection with strain b. Thus, Equation 3 expresses the susceptibility Sa 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 σb(i) in the product in Equation 3 by the fraction of the population Rb that recovered from infection with strain b:

(4) Sab(1-KabRb)e-bKabRb

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

(5) ddtSa=-νSabKabIb

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 γ/ν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:

(6) Kab=αe-|a-b|d

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 α1 quantifies the reduction of susceptibility to reinfection by the same strain and hence the overall strength of protective immunity. We shall set α=1 corresponding to perfect protection here for simplicity of presentation. Our analysis below applies equally well to the more realistic case of α<1, since in our approximation this parameter can be eliminated by rescaling Ra and Ia 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.

Viral populations escape adaptive immunity by accumulating antigenic mutations.

Via cross-reactivity, the immunity foot-print of ancestral variants (center of the graph) mediates competition between related emerging viral strains and can drive all but one of the competing lineages extinct. At high mutation rates and relatively short range of antigenic cross-reactivity, different viral lineages can escape inhibition and continue to evolve independently.

https://doi.org/10.7554/eLife.44205.003

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 Nh 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 ν to set the unit of time, fixing ν=1 in rescaled units. The remaining parameters of the model are (1) the transmission rate β - in our units the number of transmission events per infection and hence equal to the basic reproduction number R0, (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 Nh.

Phenomenology

Before proceeding with a quantitative analysis we discuss different behaviors qualitatively. Figure 3A shows several trajectories of prevalence Itot=aIa (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 Nh, 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
Extinction, speciation, and oscillations in multi-strain SIR models.

(A) Typical trajectories of infection prevalence in the parameter regimes corresponding to extinction (red), traveling wave RQS (yellow) and oscillatory RQS (green). Panels B and C schematically show parameter regimes corresponding to these qualitatively different behaviors. As explained in the text the boundaries qex,qsp of the RQS regime depend explicitly on the time scale considered (see also Figure 7). Simulation results supporting this diagram are shown in Figure 3—figure supplements 1 and 2. Panel (D) is a schematic illustrating how a novel pandemic strain (red) can settle into an endemic RQS state. As the cumulative number of infected individuals increases, the susceptible fraction decreases, and survival of the strain depends on the emergence of antigenic escape mutations (gray). The top part of the panel illustrates the population composition at a particular time point. Rare pioneering variants are q mutations ahead of the dominant variant and grow with rate xn. (Note, that boundaries of the ‘extinction’ regime in (B,C) correspond to q value close to one.) Different lineages are related via their phylogenetic tree embedded in the fitness distribution in the population.

https://doi.org/10.7554/eLife.44205.004

Large effect antigenic mutations allow transition from pandemic to seasonal dynamics

A novel virus in a completely susceptible population will initially spread with rate β-1 and the pandemic peaks when susceptible fraction falls to β-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 β-1 – the propagation threshold for the virus – and without rapid antigenic evolution the pathogen will go extinct after a time tβ-1logNh. 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.

Long-range cross-immunity results in evolving but low diversity pathogen populations

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 βaIaSa/Itot fluctuates around the rate of recovery ν=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 Kab in Equation 5

(7) Sa-1ddtSa(t)=-be-|a-b|dIb-Itot+b|a-b|dIb

where we have used that |a-b|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 βSa1. We can hence define the instantaneous growth rate of strain xa=(βSa-1)1 as its effective fitness. In this limit, the model can be simplified to

(8) ddtIa=xaIaddtxaItot

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 Kcb=Kabe-d-1Kab(1-d-1). The population susceptibility of strain c is therefore increased to

(9) Sce-(1-d-1)bKabRbSa(1-logSad)

Since the typical susceptibility is of order β-1, the growth rate of the mutant strain c is s=d-1logβ 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 β, 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, xn, is much higher than antigenic mutation rate, xnm, it takes typically

(10) τa=log(xn/m)xn

generations before a novel antigenic mutation arises in a newly arisen pioneer strain that grows exponentially with rate xn. 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 (Nm1), the fitness distribution has an approximately Gaussian shape with a variance σ22s2log(Ns)/log2(xn/m). The wave is σ/s mutations wide, while the most advanced strains are approximately q=2log(Ns)/log(xn/m) ahead of the mean (Desai and Fisher, 2007). Two contemporaneous lineages coalesce on a time scale τsw=sq/σ2=s-1log(xn/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 q1 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 Np which is related to the (fixed) host population size Nh by Np=NhItot. The average Itot depends on other parameters of the model, scaling in particular with I¯s2. Hence, it will be convenient for us to use Nhs2 as one of the relevant ‘control parameters’, replacing N of the standard TW model.

Stability and fluctuations of the RQS

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 xn 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 τsw to the ecological oscillation.

Figure 4 with 2 supplements see all
Oscillations in antigenically evolving populations.

(A) An example of the stochastic limit cycle trajectory from the fitness-class simulation. Note the rapid rise and fall of infection prevalence (blue), which causes a drop in nose fitness (yellow) which subsequently recovers (approximately linearly) during the remainder of the cycle. Fluctuations in Itot(t) and xn(t) from cycle to cycle are caused by the stochasticity of xn, that is antigenic evolution in pioneer strains. A particularly large fluctuation about τsw prior to the end, caused a large spike in prevalence, followed by the collapse of xn below zero and complete extinction. Inset (red) shows the cross-correlation between xn and σ2 which peaks with the delay τ=τsw (additional peaks reflect the oscillatory nature of the state and are displaced by integer multiples of mean period). (B) A family of limit cycles in the infection prevalence/mean fitness plane as described by Equation (11) with fixed variance. The variation of σ governed by the Equations 11 and 12 (in the deterministic limit) reduces the family to a single limit cycle (red); (C) Trajectories in the infection prevalence/nose fitness generated by the stochastic DD system in the regime above (right panel) and below (left panel) the oscillatory instability of the deterministic DD system.

https://doi.org/10.7554/eLife.44205.007

To recognize the ecological aspect of the oscillatory tendency, consider the total prevalence Itot and the mean fitness of the pathogen X=axaIa/Itot 

(11) ddtItot=XItot;ddtX=σ2-Itot

which follows directly from Equation (8). Selection on fitness variance σ2 increases X, while prevalence Itot reduces susceptibility and hence X. At fixed variance σ=σ¯ this system is equivalent to a non-linear oscillator describing a family of limit cycles oscillating about Itot=σ¯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 xn that in a continuum limit (suitable for the limit of high mutation rate) reads:

(12) ddtxn=τsw-1xn-Itot+sξ(t)

The first term on the right hand side represents the rate at which antigenic pioneer strains enter the population, τa-1, advancing the nose fitness by an increment s (with τa-1s=τsw-1xn ). The second term on the right hand side of Equation (12) represents gradual reduction of susceptibility of the host population, and ξ(t) is a random noise variable representing the stochasticity of the establishment of new strains. The Gaussian white noise ξ(t) is defined statistically by its correlation function ξ(t)ξ(0)=τa-1δ(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 xn accelerates its rate of advancement. The stabilizing factor is the subsequent increase in Itot, but to see how that comes about we must connect Equation (12) to Equation (11). The connection is provided by σ2 since it is controlled by the emergence of novel strains, that is the dynamics of the ‘nose’ xn, which impacts the bulk of the distribution after a delay τsw. Based on the analysis detailed in the Appendix 2, we approximate

(13) σ2(t)τsw-1xn(t-τsw)

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 xn(t) and σ2(t) 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 τsw-1x¯n=σ¯2=2τsw-2log(NhI¯). Small deviations decay in underdamped oscillations with frequency ω=σ¯=τsw-12log(NhI¯) if ωτsw<2π. For ωτsw>2π, 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 τsw. In Appendix 3, we derive the threshold of oscillatory instability to lie at log(NhI¯oscs)8.3 (leading to limit cycle period T1.5τsw, see Figure 4—figure supplement 1). We also find that the amplitude of the oscillation log(Imax/I¯) scales as log(NhI¯) 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 ω 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 xn increases until rising Itot catches up with it (when Itot=τsw-1xn) driving it down rapidly. Once this ‘mini-pandemic’ burns out, the population returns to the low prevalence part of the cycle Itot<τsw-1xn, when xn begins to increase again.

The rate of extinction

While in the deterministic limit the differential-delay system predicts a stable steady TW for q>qex,I¯<I¯osc and a limit cycle above 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(NhI¯) as shown in Figure 5A. The time to extinction increases dramatically in the range of q1-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 xn advance leads, after a time τsw to a rise in prevalence Itot, 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 δx of the fitness of pioneer strains becomes of the order of the mean x¯n. New mutations at the nose accumulate with rate 1/τa such that at short times t we expect δxst/τa. Hence δx becomes of the order of the mean x¯n at times τextqτ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(NhI¯). Numerical simulations, Figure 5B, confirm the dependence of τext on q and log(NhI¯). We note that the rate increase in τext with increasing q slows down in the oscillatory regime and appears to approach a power law dependence τext/τswq2.5 (albeit over a limited accessible range): presently we do not have an analytic understanding of this specific functional form.

Extinction and speciation dynamics.

(A) Simulation results for the average extinction time τext (open symbols) and the average speciation time τsp (filled symbols) as a function of pathogen diversity q. The life time of the endemic RQS state is limited by the smaller of τext and τsp. If τsp<τext, the population tends to speciate and persist, while the population is more likely to go extinct if τsp>τext. The graph shows τsp and τext rescaled with the sweep time τsw as a function of genetic diversity measured by the number of mutations q. The speciation time τsp increases with the range of cross-inhibition (d) and decreases with q. Note the agreement between the results of the fitness class-based simulation (black line in (A)) and the clone-based simulation (colored squares in (A)). (B) Extinction time over a broad range of parameters, obtained via fitness class-based simulation of population dynamics, confirms its primary dependence on q for large population sizes.

https://doi.org/10.7554/eLife.44205.010

The rate of speciation

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/q1. 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 TMRCA growing without bounds, see Figure 1 and Figure 6A. This situation contrasts the phylogeny of the single competing population, where TMRCA 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 δ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.

Speciation into antigenically distinct lineages.

(A) To speciate, two lineage have to diverge enough to substantially reduce cross-reactivity, that is T needs to be comparable to d. Inset: Illustration of the definition of time to most recent common ancestor T and the time interval δT by which T advances. (B) If such speciation happens, the host capacity - the average number of infected individuals increases two-fold. (C) The probability of such deep divergences decreases exponentially with the ratio d/q*, where effective antigenic diversity is q*=2log(Nhs2)/log(s/m). In the presence of deleterious mutations, the relevant q is not necessarily the total advance of the pioneer strains, but only the antigenic contribution. This antigenic advance q* can be computed as q*=2log(Nhs2)σag2 with antigenic variance σag2=σ2-σβ2, where σβ2 is fitness variance due to deleterious mutations. With this correction, speciation times agree with the predicted dependence (colored lines).

https://doi.org/10.7554/eLife.44205.011

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 TMRCA 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 τsw*ef(CI/q*) with ‘effective’ τsw*=τsw/(1+logq/log(s/m)) and q*=q(1+logq/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, τsw*0, when mutation rate become comparable to the selection strength m/s1.

Red Queen State is transient

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 τext and τ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 τ as a function of genetic diversity q for the case of d=50 and log(NhI¯)=6.5.

Lifetime of the RQS state.

This schematic diagram based on Figure 5A defines the ‘boundaries’ of the transient RQS qext(τ) (red line) and qsp(τ) (blue line). qext(τ) gives the value of q for which the average time to extinction is equal to τ is defined similarly but for the speciation process. Average times to speciation and extinction become equal at the critical value qc at which RQS persists the longest. For d=50, qc2.4 and this value, along with the RQS lifetime, increases with d.

https://doi.org/10.7554/eLife.44205.012

The regime of a single persistent lineage shrinks with increasing τ, for example after τ=10τsw the RQS state likely prevails between q=1.5 and 4, while (for d=50 and log(NhI¯)=6.5) it is unlikely to persist beyond τ100τsw for any q. Both the maximal RQS lifetime and corresponding critical qc, increase with increasing d.

Discussion

The epidemiological and evolutionary dynamics of human RNA viruses show a number of qualitatively distinct patterns (Grenfell et al., 2004; Koelle et al., 2011). Agents of classical childhood diseases like measles or mumps virus show little antigenic evolution, other viruses like dengue- or norovirus exist in distinct serotypes, while seasonal influenza viruses undergo continuous antigenic evolution enabling viruses of the same lineage to reinfect the same individual.

Here, we have integrated classical multi-strain SIR models with stochastic models of adaptation to understand the interplay between the epidemiological dynamics and the accumulation of antigenic novelty. The former is dominated by the most prevalent strains, while the latter depends critically on rare pioneer strains that become dominant at later times. Our model differs from that of Rouzine and Rozhnova (2018) in two aspects that are crucial to questions addressed here: To meaningfully study speciation and diversification, the model needs to allow for an high dimensional antigenic space. Similarly, fluctuations in pathogen population size determine the dynamics of extinction and this aspect can not be studied in models with constant population size. Including these aspects of the epi-evolutionary dynamics allowed to define a ‘phase’ diagram that summarizes qualitatively different behavior as a function of the relevant parameter combinations, see Figure 3B and C.

The phase diagram shows which combinations of key parameters lead to three distinct outcomes: (1) extinction (red), (2) an evolving but low diversity pathogen population (yellow and green), (3) a deeply branching and continuously diversifying pathogen population (blue). The key parameters are the size of the population log(Nhs2), the ratio of mutational effects to mutation rate log(s/m), and the cross-immunity range d. In particular, large d prevents speciation, while rapid mutation and large population sizes facilitate speciation.

In regime (2) of a low diversity but rapidly evolving pathogen population, incidence is determined by the range of cross-immunity d and by the speed of antigenic evolution which itself depends on the pathogen population size, mutation rates, and the fitness effect of novel mutations. A consistent solution of these dependencies shows that average incidence Itot decreases as d-2, while weakly depending on population size and mutation rates (see Equation A2.11), consistent with results by Rouzine and Rozhnova (2018). Typical values of the coalescent time of influenza A (2-4y), an infectious period of 5d, and a human population size 1010 result in an average annual incidence of 3–10%. This number is consistent with previous estimates of the annual attack rate of influenza (Yang et al., 2015) (which typically do not differentiate the different influenza lineages).

Of the different regimes, only extinction (1) and speciation (3) are truly asymptotic. The intermediate regimes of continuously evolving low diversity pathogen population - the Red Queen State (RQS) - are strictly speaking metastable states which eventually either go extinct or undergo branching, but in a certain regime of parameters are very long lived. In our simple model, stability against speciation on the time scale >10τsw required d10q (while stability against extinction requires q>2). These results are consistent with earlier studies that have shown that competition between lineages mediated by long-range cross-immunity can prevent diversification, effectively canalizing the population into a single lineage (Tria et al., 2005; Ferguson et al., 2003).

In practice, the range of cross-immunity required to prevent speciation might be smaller than the idealized model. Our model assumes that the pathogen population can escape immunity via many equivalent mutational path. But in reality, the number of path to escape will be limited and some path more accessible than others, which will reduce the tendency to speciate and the necessity for large d. Similarly, other factors such as population turn over and geographic heterogeneity can delay extinction.

Previous studies have shown that the rate of branching in the speciation regime increases with population size and mutation rate consistent with the phase diagram (Sasaki and Haraguchi, 2000; Koelle et al., 2011). Bedford et al. (2012) have used large-scale individual-based simulations to explore structure of influenza viruses phylogenies. Consistent with our results, they found that the speciation rate increases with the mutation rate (lowering logs/m and thereby facilitating speciation) and increasing standard deviation of mutational effects. The latter increases the typical antigenic effect of successful mutations, which decreases the radius of cross-immunity when measured in units of mutations making the population more prone to speciate.

Koelle and Rasmussen (2015) have implicated deleterious mutation load as a cause of spindly phylogenies. Deleterious mutations increase fitness variation, which results in more rapid coalescence and less antigenic diversity, which in turn reduces speciation rates. Our model can readily incorporate deleterious effects of antigenic mutations on transmission β. Such deleterious mutations reduce the selection coefficient of antigenic mutations, which in turn reduces the fitness variance σ2, see Appendix 6. After subtracting the contribution of deleterious mutations from the the fitness variance, the times to speciation follow the predicted dependence on q and d, see Figure 6C.

Outbreaks of emerging viruses that quickly infect a large fraction of the population, as for example the recent Zika virus outbreak in the Americas, fall into regime (1): In 2–3 years, large fractions of the population were infected and have developed long-lasting immunity. As far as we know, the viral population did not evolve antigenically to escape this build up of herd immunity and the virus population is not expected to continue to circulate in the Americas (O'Reilly et al., 2018).

Different influenza virus lineages, in contrast, persist in the human population, suggesting that they correspond to parameters that fall into the RQS region of the phase diagram. Furthermore, the different subtypes display quantitatively different circulation and diversity patterns that allow for a direct, albeit limited, comparison to theoretical models: subtype A/H1N1 circulated with interruption from 1918 to 2009, A(H2N2) circulated for about 10 years until 1968, A/H3N2 emerged in 1968 and is still circulating today, and the triple reassortant 2009 H1N1 lineage, called A(H1N1pdm), settled into a seasonal pattern following the pandemic in 2009. Influenza B viruses have split into two separate lineages (B/Victoria and B/Yamagata) in the 1970s (Rota et al., 1990). Phylogenetic trees of A/H3N2 and the influenza B lineages are shown in Figure 1.

The influenza B lineages tend to be more genetically diverse than the influenza A lineages with a typical time to the most recent common ancestor of around 6 compared to 3 years, see Figure 1. Influenza A/H3N2 tends to have the lowest diversity and most rapid population turnover. This difference in diversity is consistent with influenza B lineages being more prone to speciation.

The typical diversity of these viruses needs to be compared to their rate of antigenic evolution. Hemagglutination inhibition titers drop by about 0.7–1 log2 per year in A/H3N2 compared to 0.1–0.4 log2 per year for influenza B lineages (Smith et al., 2004; Bedford et al., 2014; Neher et al., 2016). Hence the ratio of the time required to lose immunity and TMRCA is similar for the different lineages, suggesting that the distinct rates of genetic and antigenic evolution can not be used as a straight forward rationalization of the speciation event of Influenza B and the lack of speciation of influenza A lineages. Nor should such an explanation be expected as there is only a single observation of speciation. We note that currently circulating A/H3N2 viruses are exceptionally diverse with a common ancestor that existed about 8 years in the past. Furthermore, the cocirculating 3c.3a and 3c.2a are antigenically distinct and it is conceivable that further antigenic evolution will result in speciation of A/H3N2 viruses.

While we have shown that the natural tendency of SIR models to oscillate couples to the instability of the nose of the pathogen fitness distribution, making a quantitative link to the observed epidemiological dynamics of the flu is difficult on account of seasonal oscillation in transmissivity. The latter confounding factor is widely believed to be the cause behind observed seasonality of the flu. Including explicit temporal variation (in β) in our model would lock the frequency of the prevalence oscillation to the seasonal cycle, possibly resulting in subharmonic modulation, yet distinguishing such a modulation on top of an already stochastic process is hard. Much remains to be done: finite birth rates, distinct age distributions (as for example is the case for the two influenza B lineages), realistic distribution of antigenic effect sizes, or very long range T-cell-mediated immunity would all be interesting avenues for future work.

Materials and methods

Clone-based simulations

Request a detailed protocol

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 Δ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 NhIa1. If the random number is smaller than mNhIaΔ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 Sb=(Sa)e-1/d. In most of the simulations, the transmissibility of different strains is held constant β. Otherwise we allow for a strain specific transmissibility that is to its parent: βb=βa-δβ with δβ>0 for the deleterious effect of antigenic mutations and βb=βmax if the mutation is compensatory. The new strain grows deterministically only if βbSb>1.

This simplified model contains six relevant parameters: transmissibility β, recovery rate ν, mutation rate of the virus m, birth/death rate of the hosts γ, the effective cross-immunity range d, and the effective size of the hosts Nh, whose empirical ranges are summarized in the Table 1. For flu and other asexual systems in RQS, βνm,γ, d1, and Nh1.

Table 1
Relevant quantities of influenza virus and parameters in multi-strain SIR model.
https://doi.org/10.7554/eLife.44205.013
SymbolDefinitionTypical values for influenzaRange in simulations
IaFraction of population infected with strain a
SaPopulation average susceptibility to strain a∼0.5
Itot=aIiTotal prevalence
I¯Average total prevalence0.005
dCross-immunity range[50, 200]
β=νR0Transmission rate∼0.5/day2
νRecovery rate∼0.2/day1 (sets unit of time)
γHost birth/death rate∼0.01/year0
Kab=e-|b-a|/dCross-immunity of strains a, b
τswCoalescent time scale/sweep time∼2-6 years
TMRCATime to most recent common ancestor∼2-10 years
δTFluctuations of TMRCA∼2-6 years
sSelection coefficient∼0.03/week[0.003, 0.05]
mBeneficial mutation rate per genome10-3/week[10-7,10-3]
NhHost population1010[106, 1012]

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

Fitness-class-based simulations

Request a detailed protocol

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 xi. At each time step, the population in each bin Ii updates to a number sampled from the Poisson distribution with parameter λi=NhIi(1+(xi-x¯)Δt) determined by mean fitness xi and a dynamic mean fitness x¯, which increases by ΔtItot, where Itot is the total infected fraction summed over all bins. When x¯ becomes larger than one bin size s, we shift the all populations to left by one bin and reset 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 mIiΔt. The typical ranges of the three parameters s, m, and Nh follow the parameters in the genealogical simulation, as documented also in Table 1.

Stochastic differential-delay simulation

Request a detailed protocol

To simulate the differential delay equations Equations (11-13), we discretize time in increments of Δt=τsw/k and update the dynamical variables χi=xn(ti) and ηi=Itot(ti) via the simple Euler scheme:

(14) χi+1=χi+Δt(χiηi)+χiqsΔtξi;
(15) ηi+1=I¯exp(τswχikτsw2kj=0kjηij),

where ξi is a Gaussian random variable with zero mean and unit variance. Mean prevalence, I¯, enters as the control parameter (which defines the time average of ηi).

Influenza phylogenies

Request a detailed protocol

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.

Appendix 1

Approximation of susceptibility

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 

(A1.1) Sa=b(1Kabσb(i))i=1bKabσb(i)+12c,bKabKacσb(i)σc(i)+...

where angular brackets denote the average over all individuals i in the population and σb(i)[0,1] denotes whether individuals i was infected with strain b in the past. This expansion assumes |Kab|1 which would hold uniformly for weak inhibition α1 but also holds for perfect inhibition for sufficiently distant strains a,b. For α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 σb(i)=Rb, that is the fraction of the population recovered from b, and σb(i)σc(i)=RbRc+ρbc where ρ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 ρbc compared to RbRc (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 Ra 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 Kab 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

(A1.2) S(x)=m=0(1K(x+mϵd))=(1K(x))em=1log1K(x+mϵd)(1K(x))eαex/dm=1emϵ=(1K(x))eαexd(eϵ1)

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 logSWI=-αe-xd(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 S1-K(x), thus keeping only the smallest term of the product representation of Sa(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 α and ϵ=1. For α=1, the ‘most recent’ approximation is better than the ‘weak inhibition’ approximation for x<d/2 but worse otherwise. For α<1, the weak inhibition approximation improves further.

Appendix 1—figure 1
Approximations to the susceptibility.

Panel A shows the effect of approximating the multiplicative effect of all past infection by logS-αe-xd(1-e-ϵ) (weak inhibition) and by ignoring all by the most recent infection. Panel B shows the population average susceptibility assuming every individual gets reinfected once the virus has evolved by dϵ for the multiplicative model, the mean field approximation, and the most recent infection approximation for different values of α as a function of the interval between infections ϵ.

https://doi.org/10.7554/eLife.44205.016

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

https://doi.org/10.7554/eLife.44205.015

Appendix 2

Differential-delay approximation of RQS dynamics

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:

(A2.1) G(λ,t)=iIi(t)eλxi(t)

where xi(t)=xn(ti)-tit𝑑tItot(t) is the current fitness of the pioneer strain that first appeared at time ti and Ii(t) is the fraction of the hosts infected by it:

(A2.2) Ii(t)=Nh-1etit𝑑txt(t)=Nh-1exn(ti)(t-ti)-0t-ti𝑑ttItot(t-t)

We next take a coarse grained view of pioneer strain establishment replacing the sum in Equation (13) by an integral over initial times tit-τ 

(A2.3) G(λ,t)=1Nh0dττa(t-τ)e(τ+λ)xn(t-τ)-0τ𝑑t(t+λ)Itot(t-t),

where 1/τa(t-τ) is the rate at which new clones are seeded at time t-τ. Let us evaluate the integral in the saddle approximation which is dominated by τ=τ* corresponding to the maximum in the exponential

(A2.4) τ*+λ=xn(t-τ*)xn(t-τ*)+Itot(t-τ*)τsw

where we have used the deterministic limit of Equation 12. To simplify presentation we shall ignore the time dependence of τsw=s-1log(xn/m) replacing xn(t-τ*) in the logarithm by the time average x¯n.

Within the saddle approximation we then have

(A2.5) logNhG(λ,t)xn(t-τsw+λ)τsw-0τsw-λ𝑑t(t+λ)Itot(t-t),

where we have omitted the logarithmic corrections for simplicity. Note that by definition G(0,t)=Itot(t).

We can now estimate fitness mean

(A2.6) x¯(t)=ddλlogG(λ,t)|λ=0=τsw[xn(tτsw)+Itot(tτsw)]0τswdtItot(tt)=xn(tτsw)0τswdtIt(tt)

and variance

(A2.7) σ2(t)=d2dλ2logG(λ,t)|λ=0=τsw[xn(tτsw)+Itot(tτsw)]+Itot(tτsw)

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

(A2.8) σ2(t)=τsw-1xn(t-τsw)

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

(A2.9) ddtx¯=xn(tτsw)+It(tτsw)It(t)=σ2(t)It(t)

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.

(A2.10) d2dt2logI(t)=τsw1xn(tτsw)Itot(t)ddtxn(t)=τsw1xn(t)Itot(t)

This system admits a family of fixed points of the form τswItot=xn=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 xn and Itot for the steady traveling wave is readily generalized to limit cycle states. Integrating the differential-delay system over one cycle yields xn=τswI. An additional relation is provided by integrating logNhG(0,t) over the cycle:

(A2.11) logNhItot=τsw22Itot

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.

https://doi.org/10.7554/eLife.44205.017

Appendix 3

Stability analysis of the differential-delay approximation

In the traveling wave case, it is natural to measure time in the units of the delay time scale τsw. The therefore define a time variable ζ via t=τswζ, the fitness variable χ via xn=τsw-1χ and the rescaled log-prevalence u via u=logτsw2I to obtain

(A3.1) d2dζ2u(ζ)=χ(ζ1)eu(ζ)ddζχ(ζ)=χ(ζ)eu(ζ)

As before, this system has a one parameter family of fixed points χ=χ¯,u=logχ¯. Note that from the traveling wave model (Desai and Fisher, 2007), we have χ¯=x¯nτsw=qlog(xn/m)=2log(Nhs2). To analyze fixed point stability we linearize and Laplace transform, yielding

(A3.2) z2δu^(z)=ezδχ^(z)χ¯δu^(z)+zδu(0)+δu(0)zδχ^(z)=δχ^(za)χ¯δu^(z)+zδχ(0)

Stability is governed by the poles of the Laplace transformed response to the initial perturbation δu(0),δu(0),δχ(0) and these poles are at the complex z that solve:

(A3.3) z=1+χ¯(1-z-e-z)/z2

Fixed point - and hence steady RQS - stability requires (z)<0 which is found for 2<χ¯<χ¯c. For χ¯>2.845 one finds (z)0 corresponding to the onset of oscillatory relaxation which turns into a limit cycle for χ¯>χ¯c16.6. The period of the limit cycle is well approximated by (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 q1. 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 χ(ti) and discrete approximation of χ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 τa. Including both corrections the Langevin equation becomes:

(A3.4) ddζχ(ζ)=11-1qχ(ζ)-11-12qd2dζ2u(ζ)+1lξ(ζ).

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

(A3.5) z2+χ¯[1+1z-1e-z+1qz2(z-1)2e-z]=0.

Solving for the onset of stability (z)=0, we find the extinction boundary qex(logNhs2) from the relation

(A3.6) χ¯=21-2qex=2log(Nhs2)

We observe that qex2 asymptotically for large logNhs2.

https://doi.org/10.7554/eLife.44205.018

Appendix 4

Stochastic form of the differential-delay approximation

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

(A4.1) ddtxn=τsw-1xn-It(t)+sξ(t),

combined with Equation (A2.5) at λ=0

(A4.2) logI(t)=τswxn(t-τsw)-0τsw𝑑ttIt(t-t).

Note that in this derivation we have avoided the need for explicitly approximating σ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.

https://doi.org/10.7554/eLife.44205.019

Appendix 5

Speciation rate as a stochastic ‘First Passage’ problem

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 x1, and the nose of branch 2 is at x2. Before the sweep time t<τsw, the cross-immunity grows mainly from the prevalent strains in the common ancestors of the two branches,

(A5.1) ddtxi=τsw-1xi-Itot+sξi. i=1,2
Appendix 5—figure 1
Process of speciation.

Left: Sketch of a branching event at t=0 with two branches 1 and 2. The fitnesses of the most fittest strains (noses) in branch 1 and 2 are x1 and x2. Branch 1 is the fitter one x1>x2. The antigenic distances from the cross-immune bulk to the noses of the two branches are d1 and d2. The Gaussian profile in fitness is illustrated in blue. Right: The fitness difference between the two branches x1-x2 is doing a biased random walk in time t of step size s with a reflecting boundary at x=0 and an absorbing boundary at x=xn.

https://doi.org/10.7554/eLife.44205.021

Later when t>τ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:

(A5.2) ddtx1=τsw1x1I1ed11/dI2ed21/d+Sξ1;ddtx2=τsw1x2I1ed12/dI2ed22/d+Sξ2;

where d11 and d22 scale roughly as q, the typical antigenic distance to the nose. In the limit d21d12d, 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 d21d12d, where I1+I2Itot 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 ζ(x1-x2)/xn. As illustrated in Appendix 5—figure 1, an establishment of nose one is a positive step of δζ=s/xn, while an establishment of nose two results in a backward step of the same size. As the mutations arrive in characteristic times τ1 and τ2 depending on the nose fitnesses, in the continuum limit, we have

(A5.3) ddtζ=τsw-1ζ+sxnξ,

where ξ is a random noise. There are two relevant boundaries: a reflecting boundary at ζ=0 where two branches switch roles in leading the fitness, and an absorbing boundary at ζ=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 ρ(ζ,t) 

(A5.4) tρ(ζ,t)=-ζ[v(ζ)ρ(ζ,t)]+ζ2[D(ζ)ρ(ζ,t)],

where the drift v and diffusivity D depend on ζ,

(A5.5) v(ζ)=1τswζ;D(ζ)1qτsw.

Solving with boundary and initial conditions,

(A5.6) ζρ(ζ,t)|ζ=0=0;ρ(ζ=1,t)=0;ρ(ζ,t=0)=δ(ζ),

we have

(A5.7) ρ(ζ,t)=n=1e-λnt/τswcnF11(1-λn2,12,q2ζ2),

where F11 is the generalized hypergeometric function, λn is the nth smallest values solving F11(1-λ2,12,q2)=0, and coefficient cn is determined by the initial condition. In long time t, the slowest mode dominates the dynamics. In the large q limit, we have λ1=1. Since F11const for ζ(0,1), the persistence probability is

(A5.8) P(T>t)ce-t/τsw.

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

(A5.9) P(d2>d/2)e-dτa/2τsw=e-d/2q.

we find the probability of a successful branching p1 to be proportional to e-d/2q.

In the phylogenetic tree, t/τa trial branchings from the backbone arrive in time t. The probability that none of them successfully speciate is thus

(A5.10) Pnsp(t)=(1-p1)t/τa=e-t/τsp,

where the waiting time for speciation event is

(A5.11) τspτswqed/2q,

as numerically verified in Figure 6.

https://doi.org/10.7554/eLife.44205.020

Appendix 6

Effect of mutations on infectivity

Suppose an antigenic mutation has a deleterious effect on infectivity reducing the latter by δβd on average. This would effectively reduce the fitness gain of antigenic innovation from s to sd(β)=s(β)-Δd, with Δd=β-1δβd. In addition let us assume that there also are compensatory mutations which restore maximal infectivity βmax. These compensatory mutation thus have a beneficial effect on fitness Δb(β)=β-1βmax-1. We assume that these mutations occur with rate mβb. In a dynamic balance state the rate of fixation of compensatory mutations would exactly balance the deleterious mutation effect on β so that τb-1Δd=τa-1Δd with the fixation rate controlled by the fitness of the leading strain via τb-1=xn/log(xnmβb). This dynamic balance is achieved at a certain value of β*<βmax, specifically βmax-β*=δβdτbτa-1 or β*=βmax-δβdr where r=log(xnm)/log(xnmβb).

The fitness of the nose of the distribution obeys

(A6.1) dxndt=sd(β)τa-1+Δbτb-1-Itot

where the 1 st term on the RHS is rate of nose advancement due to antigenetic mutations τa-1=xn/log(xnm) as before, but with reduced fitness gain sd(β). The 2nd term describes the contribution of compensatory mutations. However in the dynamic equilibrium (at β*) compensatory mutations exactly cancel the contribution the deleterious mutation contribution to s so that for the steady state we recover

(A6.2) Itot=s(β*)τag-1=s(β*)xnlogxnm

as we had for the TW driven by antigenic advancement only. The only effect is the reduction of s from s(βmax) to s(β*)=d-1logβ*.

The sweep time, τsw, upon which the fitness of the former pioneer strain comes down to the mean fitness and the nose fitness, xn, retain the TW form

(A6.3) τsw=xnItot=s(β*)log(xnm)

Following TW approximation to estimate infection prevalence ItotNh-1exp(xnτsw/2) as before one finds

(A6.4) xn=2τsw-1logCNh=2s(β*)log[Nhs2/log(xnm)]log(xnm)

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

(A6.5) σ2=τsw(sd2τa-1+Δb2τb-1)=xn[s-2Δd+Δd2s+ΔdΔbs]

but under conditions of Δd,Δbs(β*) 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:

(A6.6) qag=τswτag-1=xns(β*)=2logNhs2(β*)clogxnm

which is basically unchanged in the presence of infectivity mutations except for the expected reduction in the magnitude of s2 factor inside the logarithm. Therefore, speciation rate would be reduced, but rather weakly, via a contribution subleading in o(logNh)

https://doi.org/10.7554/eLife.44205.022

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
    Evaluation of Bovine, Cold-Adapted human, and Wild-Type human parainfluenza type 3 viruses in adult volunteers and in chimpanzees
    1. ML Clements
    2. RB Belshe
    3. J King
    4. F Newman
    5. TU Westblom
    6. EL Tierney
    7. WT London
    8. BR Murphy
    (1991)
    Journal of Clinical Microbiology 29:1175–1182.
  7. 7
  8. 8
  9. 9
  10. 10
    The Influence of Different Forms of Cross-Protective Immunity on the Population Dynamics of Antigenically Diverse Pathogens
    1. N Ferguson
    2. V Andreasen
    (2002)
    In: C Castillo-Chavez, S Blower, P van den Driessche, D Kirschner, A. -A Yakubu, editors. Mathematical Approaches for Emerging and Reemerging Infectious Diseases: Models, Methods, and Theory. New York: Springer. pp. 157–169.
    https://doi.org/10.1007/978-1-4613-0065-6_9
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
    A contribution to the mathematical theory of epidemics
    1. WO Kermack
    2. AG McKendrick
    (1927)
    Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 115:700–721.
    https://doi.org/10.1098/rspa.1927.0118
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
    Statistical Physics
    1. LD Landau
    2. EM Lifshitz
    (2013)
    Elsevier.
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
    A new evolutionary law
    1. L Van Valen
    (1973)
    Evol Theory 1:1–30.
  51. 51
  52. 52
  53. 53
  54. 54

Decision letter

  1. Katia Koelle
    Reviewing Editor; Emory University, United States
  2. Patricia J Wittkopp
    Senior Editor; University of Michigan, United States

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Phylodynamic theory of persistence, extinction and speciation of rapidly adapting pathogens" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom served as a guest Reviewing Editor, and the evaluation has been overseen by Patricia Wittkopp as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

This manuscript analyzes a model of virus evolution in a host population in response to accumulating immune memory in previously infected individuals. The main result is a phase diagram that delineates qualitatively different modes of evolution (rapid extinction, strain proliferation, and metastable traveling fitness wave dynamics) as a function of evolutionary and epidemiological parameters.

Essential revisions:

The reviewers all agreed that the presented analyses were thorough and that the results were interesting. However, they also felt that several essential revisions were required:

1) The manuscript makes use of an existing status-based SIR model when mapping epidemiological dynamics on to the traveling wave evolutionary model. This model, first, should be explained in greater detail and with more clarity in the manuscript text. Further, this multi-strain model is one of two general types of multi-strain models (the other being a history-based model formulation). Are the phase diagram results robust to other SIR model formulations, including history-based formulations and the model formulation by Lin et al., 2003? Previous relevant work (Ballesteros et al. PLOS One) indicates that this might not be the case.

2) Text should be added to refer to several previous analyses that focus on very similar questions. Of particular importance is incorporating (to a much greater extent, both in the Introduction and starting around Equation 6) text that relates to the results presented recently in Rouzine and Rozhnova, 2018, ensuring that the overlap and the differences between these two analyses is accurately described. Koelle et al., 2011,) and Andreasen and Sasaki, 2006, also address similar questions about how certain epidemiological factors (population size, breadth of cross-immunity, etc.) affect whether and how quickly antigenic diversification will occur. The results presented here should be compared to those found earlier, even if these former approaches do not consider explicitly a traveling wave model. Finally, the way in which the epidemiological dynamics are mapped to fitness and how cross-immunity is quantified is identical to the approach outlined in Luksza et al., 2014, and this paper therefore needs to be cited, most notably in the context of Equations 1-3 and 4.

3) The mathematical conditions for the mapping onto fitness waves should be made more precise. This mapping is used throughout to describe the endemic regime. However, the traveling fitness wave formalism is derived under more specific assumptions, namely, the coexistence of many small-effect mutations (which corresponds to large values of q). Some parts of the phase diagram clearly fall outside this regime. While it may still be permissible to extend the asymptotic formulae, this should be discussed. We suggest to mark the boundary of the many-mutations regime, say, given by the condition U_b \gtrsim s, as a dotted line in the phase diagram.

4) In the Discussion, it would be important to emphasize that the metastability of the endemic regime is a result of the specific assumptions of this model and to discuss potential biological effects that alter the phase diagram. In particular, much work has been devoted to discuss mechanisms that stabilize the TW regime, e.g. the ideas of short-term broad cross-immunity (Ferguson et al., 2003, or random fitness components (Tria et al., 2005).

[Editors’ note: further revisions were requested before acceptance.]

Thank you for sending your article entitled "Phylodynamic theory of persistence, extinction and speciation of rapidly adapting pathogens" for peer review at eLife. Your article is being evaluated by two peer reviewers, and the evaluation is being overseen by a guest Reviewing Editor and Patricia Wittkopp as the Senior Editor.

The primary concern stems from details that are now provided in the revised manuscript that you submitted. Given these new details, reviewer #2 is particularly concerned that the epidemiological model does not incorporate infection histories appropriately, such that the results of this analysis do not advance the literature. The reviewer's request is that you re-do the entire analysis, using an epidemiological model structure that is appropriate. We would like to give you an opportunity to respond to this request, given the set of highly divergent approaches for modeling multi-strain dynamics.

Finally, the second reviewer requested that his/her entire initial review is transmitted in full. We will follow up shortly on the transmission of this entire review, as well as those from the two other reviewers.

We appreciate the value of simplified models; however, our previous comments we designed primarily towards making the simplifying assumptions explicit and to embed the epidemiological model better into the context of previous work on the subject, which we also discussed in the review consultation. Both issues are not yet adequately addressed in the current revision, and we regard the following points as essential for publication of this manuscript:

1) The assumptions underlying the epidemiological model, in particular with respect to the applicability to influenza, should be made explicit. It is not clear to us whether the model used here is indeed a generalization of previous models, as claimed. Specific points:

a) Some justification should be given for the steps leading from a general multi-strain immunity model to their Equations 1 – 3, e.g. along the lines of Rozhnova and Rouzine. For example, the authors could try to estimate, by the order of magnitude, the error of this approximation, at least, in a simple population configuration.

b) Equation 3 should be linked to the underlying dynamical model.

c) We also note again that the first application of this model to influenza data analysis (Lukza et al., 2014), which contains a very similar form of the equations and discusses their application to influenza, should be acknowledged in the context of Equations 1 – 3.

2) A quantitative comparison of the results of this paper to Bedford et al., 2012, and to Rozhnova and Rouzine should be given, for example in a supplementary figure. Specifically:

a) The fraction N_inf/N and average selection coefficient \σ, which allow the mapping to traveling wave theory, should be compared with previous work.

b) Also, it remains important to quantify the behaviour of the number of competing strains in the phase diagram (Figure 3B, C) in some fashion (see previous comment 3). The authors map the line q=1. Their reply otherwise refers to Figure 5, but it is not clear to us to which numbers q this refers to in Figure 3 (e.g., where is the locus q=10). Figure 3B, C shows two quantities as formulas in white font which we are not sure to which lines they refer to; please clarify and give units and numbers for these quantities.

Reviewer #2:

Unfortunately, not tracking the history of individual patients, i.e., not classifying patients according to previous infecting strains, is not a biologically correct approach, even though it has been done by two groups. This is not how the immune system works. One must track the memory cells from, at least, last infection. Virus infecting an individual reacts to memory in that individual, and not in other individuals. The oversimplification changes the results substantially and cannot be relied upon. For example, the turnover rate of population should not be an important parameter of the model. The dependence of the speed on parameters changes as well. We cannot be sure about the rest.

To avoid huge phase space, the simplest meaningful approximation is to track the last infection of an individual, i.e., to introduce the recovered uninfected individuals density and classify then according to memory cells left from their last infecting strain. One can show that older infection are a small correction. Then, consider multiple dimensions (analytically or numerically does not matter) and demonstrate that one-dimensional path arises automatically. After 1D path is assured, solve the 1D model analytically. Rouzine and Rozhnova did exactly that, in the case of the long-range immunity. Their multi-dimensional simulation is located in the end of Results and Supplementary Information. I also recommend to consult the previous numeric work of Bedford.

Therefore, I have to insist that the authors redo the work properly, with tracking the last memory of infection in individual. Otherwise, no numeric comparison is possible and cannot be in the future used for data comparison.

The original review follows for the authors' information:

The manuscript analyzes a model of the virus evolution in a host population due to accumulating immune memory in previously infected individuals. The authors use the SIR model by Gog et al., 2002, to map it to results of the traveling wave theory of evolution. If the cross-immunity between the virus and the memory is long-range, the authors demonstrate that the virus persists indefinitely. The state is a Red Queen process, a never-ending chase between virus and immune system in the antigenic space. If the cross-immunity is short-range, they find out that persistent infection is either unstable or splits into new states. An effective selection coefficient which makes the mapping to traveling wave possible is calculated.

The topic of the manuscript is important and the problem is challenging. The novel part of this work, compared to a recent paper on the same topic (below), is the comparison between long-range and short-range cross-immunity, and predicting the existence of a phase diagram of various behaviors including instability and oscillations.

I have some questions regarding the choice of the initial model, the sensitivity of results to its assumptions, and the connection to the previous work, as follows:

Major comments:

1) The SIR model is not explained in the manuscript, not the original paper by Gog et a. My questions are as follows.

a) According to my understanding, a typical infection is a stochastic event. An individual exposed to virus is either infected, at the systemic level, or not. If the individual is infected, the virus reaches high loads, causes a strong immune response, leaves high numbers of memory cells, and can be transmitted with appreciable probability to another individual. If the exposed individual is not infected at the systemic level, none of these events takes place. The probability of each of the two outcomes, given the exposure dose, depends on the presence of memory cells left from the previous infections, and their genetic distance from the infecting strain. Is this the scenario that the authors had in mind?

b) The model in MS considers a population structured into recovered and infected individuals classified by genetic variants of memory cells and virus, respectively. Are these population groups mutually excluding. What is their sum, the total population?

c) Can the authors draw a multi-compartment flow diagram of the model in supplement to show the processes they have included in the model?

d) An average adult person is infected by influenza virus more than once during lifetime. Indeed, between 4% and 20% individuals are infected annually. Therefore, all infections occur in previously infected (recovered) individuals. Yet, I do not see any infection of recovered individuals in model's equation. Who is infected then? This is especially confusing given that memory cells left from previous infections are the force that drives virus evolution.

f) What is the meaning of the exponential term in susceptibility S?

2) After Gog et al., 2002, Lin et al., 2003 have proposed an alternative, more transparent SIR model. Rouzine and Rozhnova, 2018 (RR) mapped that model to the traveling wave theory.

a) How does the change to Lin et al.'s version of SIR would affect the results on the stability of infection and the oscillatory states?

b) What is the difference with RR's results in the long range immunity case?

Additional comments:

3) In contrast to authors' statement, neither them nor RR's included the fluctuations of population size. If the authors implied that RR substituted the total population size instead of infected population to the traveling wave theory, they are mistaken: RR did the same rescaling.

4) The main difference between two models is in the choice of the initial SIR model (see above). RR considered the case of long-range cross-immunity only.

5) I would write the equations for the effective selection coefficient and for the rescaling of population size in separate lines, since they are important mapping formulas.

Reviewer #3:

The revised version has adequately addressed most of the comments, except the following:

I still think it would be relevant to quantify the behaviour of the number of competing strains in the phase diagram (Figure 3B, C) in some fashion (see previous comment 3). The authors map the line q=1. Their reply otherwise refers to Figure 5, but it is not clear to me to which numbers q this refers to in Figure 3 (e.g., where is the locus q=10). Figure 3B, C shows two quantities as formulas in white font which I am not sure to which lines they refer to; please clarify and give units and numbers for these quantities.

With this amendment, I think the paper is ready for publication.

https://doi.org/10.7554/eLife.44205.024

Author response

Essential revisions:

The reviewers all agreed that the presented analyses were thorough and that the results were interesting. However, they also felt that several essential revisions were required:

1) The manuscript makes use of an existing status-based SIR model when mapping epidemiological dynamics on to the traveling wave evolutionary model. This model, first, should be explained in greater detail and with more clarity in the manuscript text. Further, this multi-strain model is one of two general types of multi-strain models (the other being a history-based model formulation). Are the phase diagram results robust to other SIR model formulations, including history-based formulations and the model formulation by Lin et al., 2003? Previous relevant work (Ballesteros et al. PLOS One) indicates that this might not be the case.

We have provided more details in the derivation of susceptibility model (Equation 3) adding a discussion which connects it to “status” and “history” based models. Our approximation is based on a factorization of the probability of different immune histories of individual hosts, which directly relates to the approach of Kryazhimskiy et al., 2007, in enabling a drastic reduction of the state space relative to the original “status” and “history” models. This approximation retains the dependence of the susceptibility of the host population to the prior history of infections, without tracking infection histories of individuals. Our resulting model of susceptibility is analogous to that derived in Kryazhimskiy et al., 2007, and the one used, prior to that, by Gog and Grenfell, 2002. These connections are fully acknowledged in references. We have also added references to and comments on Lin et al. and Ballesteros et al. Ballesteros et al. present four different two strain models (history vs status, reduced infectivity vs reduced susceptibility) and compare simulation results for these models. Only one (the status based reduced infectivity model) shows recurrent waves of infection by the mutant strain at high levels of cross-immunity. However, the parameters for cross-immunity of these models cannot be compared quantitatively. The differences between the simulations are thus due to choices of values of parameters that are phenomenological in nature and don’t have a one-to-one correspondence to reality. In the multi-strain case, we expect all models to exhibit qualitatively similar dynamics with appropriately chosen parameters – especially after reducing the high-dimensional history of status space to a linear number of variables, see Kryazhimsky et al.

2) Text should be added to refer to several previous analyses that focus on very similar questions. Of particular importance is incorporating (to a much greater extent, both in the Introduction and starting around Equation 6) text that relates to the results presented recently in Rouzine and Rozhnova, 2018, ensuring that the overlap and the differences between these two analyses is accurately described.

We now discuss the work by Rouzine and Rozhnova (R&R) at greater length. R&R map a multi-strain model in a one-dimensional antigenic landscape to a TW models of population genetics and the formulation of the model as well as the mapping to TW models is analogous to our approach. In contrast to R&R, our point of departure is a model in a high dimensional antigenic space and we use this model to show how the effectively one-dimensional TW emerges rather than introducing this as model assumption. R&R use their model to infer parameters through explicit comparison to influenza diversity data, which we don’t attempt. Instead, we use this model to explore the processes of extinction and speciation which can’t be studied in one-antigenic dimension with constant population size examined in R&R.

Koelle et al., 2011, and Andreasen and Sasaki, 2006, also address similar questions about how certain epidemiological factors (population size, breadth of cross-immunity, etc.) affect whether and how quickly antigenic diversification will occur. The results presented here should be compared to those found earlier, even if these former approaches do not consider explicitly a traveling wave model.

We have added the two suggested references along with a brief description in the Discussion section.

Finally, the way in which the epidemiological dynamics are mapped to fitness and how cross-immunity is quantified is identical to the approach outlined in Luksza et al., 2014, and this paper therefore needs to be cited, most notably in the context of Equations 1-3 and 4.

We now explicitly refer to Luksza et al. in this context.

3) The mathematical conditions for the mapping onto fitness waves should be made more precise. This mapping is used throughout to describe the endemic regime. However, the traveling fitness wave formalism is derived under more specific assumptions, namely, the coexistence of many small-effect mutations (which corresponds to large values of q). Some parts of the phase diagram clearly fall outside this regime. While it may still be permissible to extend the asymptotic formulae, this should be discussed. We suggest to mark the boundary of the many-mutations regime, say, given by the condition U_b \gtrsim s, as a dotted line in the phase diagram.

We have added a paragraph (bottom of column 1 on p 6) discussing typical values of q in influenza populations and making explicit comment that flu evolution is not likely to be in the asymptotic regime of large q. Nevertheless, we and others before us (references in the text), find that qualitative features of the TW models extend to modest values of q (as seen through comparison with direct simulation). Instead of adding a dotted line to indicate crossover in the phase diagrams of Figure 3B, C we added a note in the caption explaining that boundary of the”extinction” regime corresponds to qo(1). More detailed view of the range of q corresponding to the RQS state is provided by Figure 5. We have added a remark about the range of q to the discussion of Figure 5 in the text, noting that (asymptotic) region of large q can only be reached in the limit of long-range cross-inhibition.

4) In the Discussion, it would be important to emphasize that the metastabilty of the endemic regime is a result of the specific assumptions of this model and to discuss potential biological effects that alter the phase diagram. In particular, much work has been devoted to discuss mechanisms that stabilize the TW regime, e.g. the ideas of short-term broad cross-immunity (Ferguson et al., 2003) or random fitness components (Tria et al., 2005).

In a strict sense, the endemic regime is always metastable as explicitly shown in Figure 5 which shows that RQS always goes away eventually either through extinction or through speciation. We have added text to further explain the meaning of Figure 5 to the reader. However the dependence of the extinction rate on model parameters might change through the addition of different model components. Instead of the logarithmic dependence of extinction rate on population size, this dependence might become polynomial or exponential in models that dampen population size fluctuations more strongly. We have added a short paragraph emphasizing that long range immunity and other features that dampen oscillations (population turn over, geographic structure, etc.) will tend to stabilize the RQS state.

[Editors' note: further revisions were requested prior to acceptance.]

Response to query regarding additional revisions:

Thank you for giving us an opportunity to respond to the points that came up after re-review of our manuscript. We are glad to read the reviewer #3 considers our revisions satisfactory and we can readily address the remaining request for clarification.

However, in response to our revised manuscript reviewer #2 now requests that we redo the analysis using a model that explicitly keeps track of the infection histories of all individuals. Based on previous work by several authors and our current understanding of influenza epidemiology and immunology, we have argued in the manuscript that this level of detail is not necessary and that the essential aspects of the dynamics are captured by simpler models that instead of individual infection histories keep track of susceptibilities to different strains in the population. We stand by this argument.

While infection history of an individual is important for predicting her/his susceptibility to infection by a given strain, the effective rate of spreading of that strain depends on the average susceptibility of individuals. This averaging makes our “mean field”-type approximation for population-wide susceptibility a natural first step. Let us restate here that there is a direct relation between the form of susceptibility that we use in the our work and the infection history description. It is easy to see that our expression for S given in Equation 3 is essentially exact in the limit of weak inhibition. The probability pa,i of individuals to become infected by strain a can be expressed as pa,i = Πb(1 − αKabσbi) where binary σbi ∈ [0,1] denotes the infection history of individual i and α ≤ 1 sets the strength of inhibition (so that weak inhibition corresponds to α ≪ 1 is the cross-immunity kernel defined in the manuscript. Population wide susceptibility is the population average of pa,i of all individuals i:

⟨" close="⟩" separators="|">pa,i=⟨" close="⟩" separators="|">∏b1-αKabσbi=1-α∑bKab⟨" close="⟩" separators="|">σbi

+α22∑b∑c≠bKabKac⟨" close="⟩" separators="|">σbiσci

(1)

The term hσbii = Rb is the fraction of people recovered from strain b. Correlations ρbc between infections with strain b and c show up in the second order term hσbiσcii = RbRc + ρbc. This simple derivation effectively captures the content of Kryazhimskiy et al., 2007, order-1 independence closure which assumes ρbc = 0. We cite this work to make connection with prior work on the subject. In this case, we can simply exponentiate the expression to obtain our Equation 3, correct to order α2 (which is small for weak inhibition – we discuss the case of strong inhibition below):

Sa=⟨" close="⟩" separators="|">pai=e-∑bαKabRb (2)

We note that given this form of susceptibility and the homogeneity property of Equations 2-3 in the manuscript, parameter α can be eliminated by rescaling of R and I fractions, i.e. can be absorbed into effective host population size and does not explicitly appear in model analysis and simulation.

Several facts about influenza in human populations suggest that the weak inhibition approximation is a reasonable starting point for modeling population scale behavior.

- Seasonal flu epidemics involve a large number of strains, a particular strain infects only a small fraction of the population. Hence the Ra are small and correlation effects are of minor importance.

- 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., 1986]. This possibility of homotypic re-infection shows that αKab are substantially smaller than 1, supporting our approximation of population wide susceptibility, as discussed above.

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

- Correlation in infection history induced by immunity are further reduced by the variation in exposure history through geography and variation in contact networks.

Having provided the reasons why we think weak-inhibition approximation is appropriate, we note that the utility of (3) as a model of susceptibility does not end there! To wit, this approximation correctly captures the cross-inhibiting contribution of distant strains (on account of Kab for those strains being much less than 1, so that quadratic terms in K can be neglected compared to linear ones). Hence, even in the case of strong immunity α ≈ 1, only correlation terms involving close strains could contribute. In the traveling wave description of continuous adaptation, most relevant effects involve the spreading of the newly emerging antigenic variants in the “nose” of the fitness distribution. While these strains are antigenically close they occur at low frequencies, so that Ra ≪ 1 and the correlations can again be neglected, the population-wide susceptibility to these strains is dominated by the effect of more distant strains (from further in the past).

Last but not least, expression (3) is an example of a “Mean Field Theory” type of approximations that replace the average of an exponential by an exponential of the average, neglecting correlation effects. This type of an approximation has an illustrious record of providing valuable insight into complex phenomena and are universally accepted in Physics: they are well recognized as the starting point for mathematical modeling.

In contrast, the proposal by reviewer #2 to simplify the problem by only tracking the last infection (that is dropping all terms in Πb(1 − Kabσbi) but the most recent one) is completely arbitrary and counterfactual. Fonville et al., 2014, have shown that immunity is maintained over decades and new infection results in a back-boost rather than a reset of the immunity landscape. It is also problematic as it would artificially facilitate speciation: It would increase susceptibility to viruses from sister clades since infection with one virus “wipes out” immune memory induced by a common ancestor.

Other points raised by reviewer #2 include:

- population turnover rate:the population turnover rate γ is NOT an important parameter of the model (as the reviewer pointed out) – and we never claimed it is. In fact, we explicitly set it to zero and it doesn’t feature in any of our conclusions.

- one- vs multidimensional trajectories:The reviewer suggests to “consider a multi-dimensional model […] and assure that a one dimensional path arises automatically. […] solve the 1D model analytically”. However, we show that a 1D path does not arise automatically and delineate conditions in which it does. Within these parameter ranges, we then investigate the model analytically as suggested by the reviewer. We are well aware of the work by Bedford et al., 2012, and Rouzine and Rozhnova, 2018. In Rouzine and Rozhnova, 2018, the 1D traveling wave is inherent in the model by either allowing for only one antigenic direction or imposing a preferred direction of antigenic escape (section 2.1 of the appendix of Rouzine and Rozhnova, 2018, – there seems to be inconsistent notation and a confusion of recovered and susceptible classes in the appendix).

- Fluctuations in population size:Our model explicitly accounts for fluctuations in the total number of infected individuals Itot. We don’t understand how the reviewer got the impression that our model assumes a constant (viral) population size.

Independent of the details of the epidemiological model, our work makes novel and important contributions to our understanding of pathogen evolution and dynamics:

- We show how qualitative features of epidemiological and evolutionary dynamics of rapidly adapting pathogens depend on parameters in a generic model of evolution in a high dimensional antigenic and genetic space. We delineate parameter regimes corresponding to speciation/diversification, single strain persistence, and extinction after a pandemic in a unified frame work. All previous analysis of this problem were either restricted to low dimensional spaces, a few strains, or purely numerical in nature.

- We connect the population genetics of rapid adaptation and with multistrain models of pathogens in a population that builds up immunity and show how epidemiological oscillation couple to the evolutionary dynamics of the pathogens. Rouzine and Rozhnova, 2018, don’t account for this crucial interaction – they consider a time invariant traveling wave solution of constant size.

The editor and reviewers have again carefully studied the previously revised manuscript and the authors' response to previous criticism. We appreciate the value of simplified models; however, our previous comments we designed primarily towards making the simplifying assumptions explicit and to embed the epidemiological model better into the context of previous work on the subject, which we also discussed in the review consultation. Both issues are not yet adequately addressed in the current revision, and we regard the following points as essential for publication of this manuscript:

1) The assumptions underlying the epidemiological model, in particular with respect to the applicability to influenza, should be made explicit. It is not clear to us whether the model used here is indeed a generalization of previous models, as claimed.

Specific points:

a) Some justification should be given for the steps leading from a general multi-strain immunity model to their Equations 1 – 3, e.g. along the lines of Rozhnova and Rouzine. For example, the authors could try to estimate, by the order of magnitude, the error of this approximation, at least, in a simple population configuration.

We have added an Appendix 1 that explains in detail how a model based on individual infection histories reduces approximately to our “mean-field-theory” (MFT) type model of cross-immunity. We discuss the underlying assumptions in the light of known influenza immunology. Diversity of immune responses, non-perfect immune protection, and long-lasting immunity with known back-boost effects all suggest that the entire infection history is important at the level of individuals, but that population level dynamics is well described by a factorized distribution of histories. Following the suggestion of the decision letter, we evaluate the accuracy of our approximation (in Figure 8) for a specific “simple population configuration”, adopting for this purpose the scenario of periodic reinfection of individuals specifically considered by Rozhnova and Rouzine (R&R), which allows us to explicitly compare with the approximation used in their paper. In this R&R scenario, infection histories of individuals contain strong temporal correlations so that it may be expected to be a tough case for the MFT approximation (which neglects correlations in evaluating population averages). Nevertheless, we show (in Figure 8) that our approximation to the population averaged susceptibility is quite accurate in the relevant parameter regime. In particular, it is more accurate than the “most recent” approximation advocated for by one of the reviewers and used by Rouzine and Rozhnova.

b) Equation 3 should be linked to the underlying dynamical model.

Equation 3 is linked to the dynamical model by simple differentiation, which reproduces the corresponding equation in Gog and Grenfell, 2002 (in the limit of slow population turnover). This is now made explicit in Equation 4.

c) We also note again that the first application of this model to influenza data analysis (Lukza et al., 2014), which contains a very similar form of the equations and discusses their application to influenza, should be acknowledged in the context of Equations 1 – 3.

In response to the previous decision, we had included a reference to Lukzsa and Lassig, 2014, in the context of the cross-immunity function (prev Equation 4 as had been requested). The basic model defined in Equations 1-3 was introduced specifically in the context of influenza by Gog and Grenfell, 2002 and L&L refer to Gog and Grenfell for the definition of their model. We now explicitly point out that L&L used this model for influenza as well.

2) A quantitative comparison of the results of this paper to Bedford et al., 2012, and to Rozhnova and Rouzine should be given, for example in a supplementary figure. Specifically: a) The fraction N_inf/N and average selection coefficient \σ, which allow the mapping to traveling wave theory, should be compared with previous work.

We have added a discussion of the infected fraction predicted by our model. The values predicted by our model are compatible with observation and up to logarithmic factors agree with predictions by Rouzine and Rozhnova. Similar results hold for the average selection coefficient s (called σ in R&R), where we predict the same qualitative dependence on parameters. The main difference (we predict a weaker dependence on R0) can be traced to the “most-recent” approximation to the immunity history made by R&R. We show (see above) that our approximation is more accurate. We also discuss the results of Bedford et al., 2012m who studied dependence of genetic diversity and the tendency to speciate in large scale agent based simulations. Their observations are consistent with our analytic results.

We would like to stress, however, that the primary contribution of our work is notin an recapitulation of specific parameters of seasonal influenza virus, but in elucidating general properties of the coupling of evolutionary and epidemiological dynamics.

b) Also, it remains important to quantify the behaviour of the number of competing strains in the phase diagram (Figure 3B, C) in some fashion (see previous comment 3). The authors map the line q=1. Their reply otherwise refers to Figure 5, but it is not clear to us to which numbers q this refers to in Figure 3 (e.g., where is the locus q=10). Figure 3B, C shows two quantities as formulas in white font which we are not sure to which lines they refer to; please clarify and give units and numbers for these quantities.

We agree that our previous presentation on how q is related to the phase diagram was not optimal. We have now expanded our explanation of the regime boundaries in the phase diagram of Figure 3 specifically discussion the nature of the “critical” values of q (associated with these boundaries) and stressing how the transitions to extinction and speciation regimes depend on the time scale of observation. To this effect, we have added a subsection “Red Queen State is transient” (just before the Discussion section) and another figure (Figure 7) that explicitly shows the regimes of extinction and speciation in the plane of q and the observation time scale. We have also edited Figure 3B, C to clarify the labelling of the regime boundaries.

https://doi.org/10.7554/eLife.44205.025

Article and author information

Author details

  1. Le Yan

    Kavli Institute for Theoretical Physics, University of California, Santa Barbara, Santa Barbara, United States
    Contribution
    Conceptualization, Resources, Software, Formal analysis, Investigation, Methodology, Writing—original draft, Writing—review and editing
    For correspondence
    lyan@kitp.ucsb.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-1323-0545
  2. Richard A Neher

    Biozentrum, University of Basel, Swiss Institute for Bioinformatics, Basel, Switzerland
    Contribution
    Conceptualization, Resources, Software, Formal analysis, Investigation, Methodology, Writing—original draft, Writing—review and editing
    For correspondence
    richard.neher@unibas.ch
    Competing interests
    Reviewing editor, eLife
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-2525-1407
  3. Boris I Shraiman

    Kavli Institute for Theoretical Physics, University of California, Santa Barbara, Santa Barbara, United States
    Contribution
    Conceptualization, Formal analysis, Funding acquisition, Investigation, Methodology, Writing—original draft, Writing—review and editing
    For correspondence
    shraiman@kitp.ucsb.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0886-8990

Funding

Simons Foundation (326844)

  • Boris I Shraiman

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Senior Editor

  1. Patricia J Wittkopp, University of Michigan, United States

Reviewing Editor

  1. Katia Koelle, Emory University, United States

Publication history

  1. Received: December 7, 2018
  2. Accepted: September 14, 2019
  3. Accepted Manuscript published: September 18, 2019 (version 1)
  4. Version of Record published: October 23, 2019 (version 2)

Copyright

© 2019, Yan et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 720
    Page views
  • 135
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Neuroscience
    2. Physics of Living Systems
    Andrew R McKinstry-Wu et al.
    Research Article
    1. Physics of Living Systems
    2. Structural Biology and Molecular Biophysics
    Christina M Caragine et al.
    Research Article