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
9 figures, 1 table and 1 additional file


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.
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.
Figure 3 with 2 supplements
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.
Figure 3—figure supplement 1
‘Phase diagram’ in the clone-based simulation.

The data points are assigned to phases using the following criteria. The extinction phase (red): τext<10τsw; The speciation phase (blue): τsp<10τsw; The transient RQS regime (yellow and green): otherwise. Red and blue lines correspond to the extinction boundary qext and speciation boundary qsp respectively, defined for τ=10τsw. The black solid line corresponds to the onset of spontaneous oscillation log(NhI¯)=8.3. Color in the transient phase labels the amplitude of prevalence var(logItot), increasing from yellow to green. Simulation set d=50.
Figure 3—figure supplement 2
Value of q across the ‘Phase diagram’.

This figure shows the same data as one with each data point colored by the corresponding value of q, which is constant along the straight of different angle.
Figure 4 with 2 supplements
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.
Figure 4—figure supplement 1
Amplitude and time scale of oscillations.

Top: Amplitude of epidemic circulations in logarithm of total prevalence. The amplitude predicted in the fitness class-based (FC) simulation is consistent with the clone-based simulation, bounded from below by the amplitude in the deterministic differential-delay approximation (DDE), which sets on at the spontaneous oscillation threshold log(NhI¯s)osc=8.3, indicated by the dashed line. When the noise is properly considered as in Equation (15), the amplitude can also be predicted from the stochastic differential-delay equations (DDE). Bottom: Period of epidemic oscillation, also bounded below by the limit cycle period of the deterministic DDE.
Figure 4—figure supplement 2
Speed of the traveling wave.

(A) In our model, the speed of traveling wave varies in time. The instantaneous speeds measured from immune adaptation, Fisher’s theorem, and nose advance could be different at any given moment. (B) But in the phase with steady traveling waves, the time-averaged speeds are consistent for different parameters. (C) The average nose position and (D) the average speed of adaptation (in the clone-based simulation) are consistent with the prediction of the traveling wave theory.
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.
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).
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.
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 ϵ.
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.


Table 1
Relevant quantities of influenza virus and parameters in multi-strain SIR model.
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]

Additional files

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)

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

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

  1. Le Yan
  2. Richard A Neher
  3. Boris I Shraiman
Phylodynamic theory of persistence, extinction and speciation of rapidly adapting pathogens
eLife 8:e44205.