# Phylodynamics of rapidly adapting pathogens: extinction and speciation of a
Red Queen.

Le Yan<EMAIL_ADDRESS>Kavli Institute for Theoretical Physics, UC Santa
Barbara  Richard A. Neher  Biozentrum, University of Basel, Switzerland  Boris
I. Shraiman  Kavli Institute for Theoretical Physics, UC Santa Barbara

###### Abstract

Rapidly evolving pathogens like influenza viruses can persist by accumulating
antigenic novelty fast enough to evade the adaptive immunity of the host
population, yet without continuous accumulation of genetic diversity. This
dynamical state is often compared to the Red Queen evolving as fast as it can
just to maintain its foothold in the host population: Accumulation of
antigenic novelty is balanced by the build-up of host immunity. Such Red Queen
States (RQS) of continuous adaptation in large rapidly mutating populations
are well understood in terms of Traveling Wave (TW) theories of population
genetics. Here we shall make explicit the mapping of the established Multi-
strain Susceptible-Infected-Recovered (SIR) model onto the TW theory and
demonstrate that a pathogen can persist in RQS if cross-immunity is long-
ranged and its population size is large populations allowing for rapid
adaptation. We then investigate the stability of this state focusing on the
rate of extinction and the rate of “speciation” defined as antigenic
divergence of viral strains beyond the range of cross-inhibition. RQS states
are transient, but in a certain range of evolutionary parameters can exist for
the time long compared to the typical time to the most recent common ancestor
($T_{MRCA}$). In this range the steady TW is unstable and the antigenic
advance of the lead strains relative to the typical co-circulating viruses
tends to oscillate. This results in large fluctuations in prevalence that
facilitate extinction. We shall demonstrate that the rate of TW fission into
antigenically uncoupled viral populations is related to fluctuations of
$T_{MRCA}$ and construct a “phase diagram” identifying different regimes of
viral phylodynamics as a function of evolutionary parameters.

In a host population that develops long lasting immunity against a pathogen,
the pathogen can persist either 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
populations adapt rapidly and reinfect most humans multiple times during their
lifespan. Continuous adaptation of influenza viruses is facilitated by high
mutation rates and it is common that many different variants of the same
subtype co-circulate. Nevertheless, almost always a single variant eventually
outcompetes the others such that diversity within one subtype 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. Fig. 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. Other pathogenic RNA viruses that typically
don’t 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 1970ies which by now are
antigenically distinct (Rota  _et al._ , 1990) and maintain intermediate
levels of diversity. Thus, one may wish to understand under what conditions a
virus can continuously evolve in competition with host immunity, neither going
extinct nor spawning diverging lineages.

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 70ies, the population split into
two lineages known as Victoria (B/Vic) and Yamagata (B/Yam). The graphs on the
right quantify diversity via the time to the most recent common ancestor
$T_{MRCA}$ 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.

A variety of mechanisms have been proposed that can lead to spindly
phylogenies in rapidly evolving viral populations. A low dimensional antigenic
space might limit the number of distinct directions in which the population
can evolve thereby preventing diversification (Andreasen  _et al._ , 1997;
Andreasen and Sasaki, 2006; Gog and Grenfell, 2002). This scenario seems
inconsistent with estimates of the dimensionality of “antigenic shape space”
($d=5$ (Perelson and Oster, 1979)) and the number of distinct positions in
surface proteins that evolve under immune selection (Koel  _et al._ , 2013;
Neher  _et al._ , 2016; Bhatt  _et al._ , 2011). Others have shown that
competition between lineages and long range cross-immunity between strains can
prevent diversification, effectively canalizing the population into a single
lineage (Bedford  _et al._ , 2012; Tria  _et al._ , 2005; Ferguson  _et al._ ,
2003). Deleterious mutational load also facilitates the maintenance of a
single lineage (Koelle and Rasmussen, 2015).

While these previous studies have identified competition mediated by cross-
immunity as important determinants of lineage structure, these prior works
mostly relied on simulations with the aim to recapitulate patterns observed
for human seasonal influenza viruses. Here, we show that generic stochastic
models of antigenic evolution with finite cross-immunity are compatible with
the spindly phylogenies of influenza viruses if the range of cross-immunity is
large compared to population diversity. If cross-immunity decays more rapidly
with mutational distance, viral population becomes prone to speciation.
Finally, if antigenic evolution is too slow, the virus will go extinct after a
brief pandemic. This rich behavior is controlled by three dimensionless
parameters that have a direct relationship to the parameters of the multi-
strain SIR model. We show how multi-strain SIR models relate to traveling wave
models of adaptive evolution (Tsimring  _et al._ , 1996; Rouzine  _et al._ ,
2003; Desai and Fisher, 2007; Neher, 2013) and how extinction and speciation
relate to oscillations of prevalence and genetic diversity.

### 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 $I_{a}$ infected with viral
strain $a$ changes according to

$\frac{d}{dt}I_{a}=\beta S_{a}I_{a}-(\nu+\gamma)I_{a}$ (1)

where $\beta$ is the transmissibilty, $S_{a}$ is the fraction of the
population that is susceptible to strain $a$, $\nu$ is the recovery rate, and
$\gamma$ is the population turnover rate. Susceptibility of strain $S_{a}$
depends on the fraction $R_{b}$ of the population recovered from infections
with strain $b$

$\displaystyle S_{a}$ $\displaystyle=$ $\displaystyle
e^{-\sum_{b}K_{ab}R_{b}}$ (2) $\displaystyle\frac{d}{dt}R_{a}$
$\displaystyle=$ $\displaystyle\nu I_{a}-\gamma R_{a}$ (3)

The matrix $K_{ab}$ quantifies the cross-immunity to strain $a$ elicited by
infection with strain $b$, while Eq. 3 describes recovery. This model
approximates the population susceptibility by the average number of infections
with each strain $b$ and ignores the explicit infection history of
individuals.

Figure 2: 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.

New strains are constantly produced by mutation with rate $m$. The novel
strain will differ from its parent at one position in its genome. We shall
consider only mutations that contribute to the loss of immune recognition and
assume that cross-immunity decays exponentially with the number of mutations
that separate two strains:

$K_{ab}=e^{-\frac{|a-b|}{d}}$ (4)

where $|a-b|$ denotes the mutational distance between the two strains and $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 distance between two contemporaneous strains is on average twice
the distance to their common ancestor.

Cross-immunity and the mutation/diversification process are illustrated in
Fig. 2. An infection with a virus (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.

The above model was formulated in terms of the deterministic Eqs. 1 and 2. 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. This stochasticity can be captured by interpreting
the terms in Eq. 1 as rates of discrete transitions in a total population of
$N_{h}$ hosts. The latter manifestation of stochasticity is particularly
important for novel mutant strains that are rare. The great majority of novel
strains are quickly lost to stochastic extinction 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 is modeled stochastically, while
common strains follow the deterministic dynamics, see Clone-based simulation
in Methods. To simplify the analysis, we will assume that the population turn-
over rate $\gamma$ is small compared to other time scales of the system and we
will set $\gamma=0$. We will use the recovery rate $\nu$ to set the unit of
time, fixing $\nu=1$ in rescaled units. The remaining parameters of the model
are 1) the transmission rate $\beta$ \- in our units the number of
transmission events per infection and hence equal to the basic reproduction
number $R_{0}$, 2) the mutation rate $m$, 3) the range of cross-immunity $d$
measured as a typical number of antigenic mutations needed for an $e$-fold
drop of cross-inhibition, and 4) the host population size $N_{h}$, which
controls the number of opportunities for adaptive evolution and the time it
takes for a newly mutated strain to reach population frequency (i.e.
prevalence) of order one.

Before proceeding with a quantitative analysis we discuss different behaviors
qualitatively. Fig. 3A shows several trajectories of prevalence
$I_{tot}=\sum_{a}I_{a}$ (i.e. total actively infected fraction) for several
different parameters. Depending on the range of cross-immunity, the pathogen
either goes extinct after a single pandemic (red lines) or settles into a
persistently evolving state, the Red Queen State (RQS) traveling wave Van
Valen (1973). In large populations the RQS exhibits oscillation in prevalence.
The RQS is only transient, but its lifetime increases with the host population
size $N_{h}$ and the mutation rate $m$. To quantitatively understand the
dependence on parameters (shown schematically in Fig. 3BC), we will further
simplify the model and establish a connection to models of rapid adaptation in
population genetics.

Figure 3: Phases in the multi-strain SIR simulation and the transition from
the pandemic to endemic. (A) Typical trajectories of infection prevalence in
the regime of extinction (red), traveling wave RQS (yellow) and oscillatory
RQS (green) which are realized for different values of model parameters
($N_{h}$ and $m$). (B,C) The extended multi-strain SIR model supports a long-
lived but transient RQS regime in an intermediate range of parameters, flanked
be the regime of deterministic extinction (red) and the regime of continuous
branching and diversification – the “speciation” regime (blue). Panel (B)
shows dependence on population size and mutation rate: small populations at
low mutation rates go extinct, while large and rapidly mutating populations
continuously produce divergent lineages (speciation). Simulation results
supporting this phase diagram are shown in S2. Panel (C) presents dependence
on the range of antigenic cross-inhibition ($d$) and the population size. The
range of cross inhibition controls the crossover between RQS and speciation
regimes. The RQS regime itself undergoes a transition from a steady travelling
wave (yellow) to a limit cycle oscillation (green) with increasing population
size. Panel (D) is a schematic illustrating the dynamics of the RQS state. A
novel pandemic strain (red) initially expands in fully susceptible population.
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 (grey). The pathogen may then settle into an
endemic, Red Queen-type state where new antigenic variants are continuously
produced. 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 $x_{n}$. Different lineages are
related via their phylogenetic tree embedded in the fitness distribution in
the population.

### Large effect antigenic mutations allow transition from pandemic to
seasonal dynamics

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

Persistence and transition to an endemic state is only possible if the
pathogen can evade the rapid build-up of immunity via a small number of large
effect antigenic mutations. This process is indicated in Fig. 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 and Holt, 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 Fig. 3), the
average rate of new infections $\beta\sum_{a}I_{a}S_{a}/I_{tot}$ fluctuates
around the rate of recovery $\nu=1$ (in our time units). This balance is
maintained by the steady decrease in susceptibility due to rising immunity
against resident strains and the emergence of antigenically novel strains, see
Fig. 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 differentiate Eq. 2 with respect to
$t$ and divide by $S_{a}$:

$S_{a}^{-1}\frac{d}{dt}S_{a}(t)=-\sum_{b}e^{-\frac{|a-b|}{d}}I_{b}\approx-
I_{tot}+\sum_{b}\frac{|a-b|}{d}I_{b}$ (5)

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

$\begin{split}\frac{d}{dt}I_{a}&=x_{a}I_{a}\\\ \frac{d}{dt}x_{a}&\approx-
I_{tot}\end{split}$ (6)

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

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

$S_{c}\approx e^{-(1-d^{-1})\sum_{b}K_{ab}R_{b}}\approx
S_{a}\left(1-\frac{\log S_{a}}{d}\right)$ (7)

Since the typical susceptibility is of order $\beta^{-1}$, the growth rate of
strain of the mutant $c$ is $s=d^{-1}\log\beta$ higher than that of its
parent. The growth rate increment plays the role of a selection coefficient in
typical population genetic models and corresponds to the step size of the
fitness distribution in Fig. 3D. In this model, 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 Methods.

The simplified model in Eq. 6 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, 2013) 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 ahead of the most common strains. This “nose” of the fitness
distributions contains the strains that dominate in the future and the only
adaptive mutations that fixate in the population arise in pioneer strains in
the nose. Consequently, the rate with which antigenic mutations establish in
the population is controlled by the rate at which they arise in the nose
(Desai and Fisher, 2007). If the growth rate at the nose of the distribution,
$x_{n}$, is much faster than antigenic mutation rate, $x_{n}\gg m$ it takes
typically

$\tau_{a}=\frac{\log(x_{n}/m)}{x_{n}}$ (8)

generations before a novel antigenic mutation arises in a newly arisen pioneer
strain that grows exponentially with rate $x_{n}$. The advancement of the nose
is balanced rapidly by the increasing population mean fitness. Previous
analyses have shown that in the limit of very large population $N\gg 1$, the
fitness distribution has an approximately Gaussian shape with a variance
$\sigma^{2}\approx 2s^{2}{\log(Ns)}/\log^{2}(x_{n}/m)$. The wave is $\sigma/s$
mutations wide, while the most advanced strains are approximately
$q=2\log(Ns)/\log(x_{n}/m)$ ahead of the mean (Desai and Fisher, 2007). Two
contemporaneous lineages coalesce on a time scale $\tau_{\rm
sw}=sq/\sigma^{2}=s^{-1}\log(x_{n}/m)$ and the branching patterns of the tree
resemble a Bolthausen-Sznitman coalescent rather than a Kingman coalescent
(Desai  _et al._ , 2013; Neher and Hallatschek, 2013). We note that, while
parameter $N$ in the TW analysis summarized above is the fixed population
size, the corresponding entity in our SIR model is the pathogen population
size $N_{p}$, which is related to the (fixed) host population size $N_{h}$ by
$N_{p}=N_{h}{\bar{I}}$ where $\bar{I}$ is the average viral prevalence, which
itself depends on other parameters of the model, scaling in particular with
$s^{2}$. Hence, it will be convenient for us to use $N_{h}s^{2}$ as one of the
relevant “control parameters”, replacing $N$ of the standard TW model. A
recent related work, that also explicitly maps a multi-strain SIR models to
the TW models, but does not consider the role of population size fluctuations
(Rouzine and Rozhnova, 2018).

Figure 4: Oscillatory RQS. (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 $I_{tot}(t)$ and $x_{n}(t)$ from cycle to cycle are
caused by the stochasticity of $x_{n}$, i.e. antigenic evolution in pioneer
strains. A particularly large fluctuation about $\tau_{\rm sw}$ prior to the
end, caused a large spike in prevalence, followed by the collapse of $x_{n}$
below zero and complete extinction. Inset (red) shows the cross-correlation
between $x_{n}$ and $\sigma^{2}$ which peaks with the delay $\tau=\tau_{\rm
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 infection prevalence/mean fitness plane as described by Eq. (9) with fixed
variance. The variation of $\sigma$ governed by the Eqs. 9-10 (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.

### 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: The first effect is ecological: a bloom of pathogen
depletes susceptible hosts leading to a crash in pathogen population and a
tendency of the population size to oscillate London and Yorke (1973). The
second effect is evolutionary: higher nose fitness $x_{n}$ begets faster
antigenic evolution and vice versa, resulting in an apparent instability in
the advancement of the antigenic pioneer strains. This instability has been
identified in the study of adaptive traveling waves (Fisher, 2013). In our
epidemiological model, fluctuations in the rate of antigenic advance of the
pioneer strains couple, with a delay of $\tau_{\rm sw}$, to the ecological
oscillation.

To recognize the origin of the oscillatory tendency, consider the total
prevalence $I_{tot}$ and the mean fitness of the pathogen
$X=\sum_{a}x_{a}I_{a}/I_{tot}$

$\frac{d}{dt}I_{tot}=XI_{tot};\ \ \ \ \ \frac{d}{dt}X=\sigma^{2}-I_{tot}$ (9)

At fixed variance $\sigma={\bar{\sigma}}$ this system is equivalent to a non-
linear oscillator, describing a family of limit cycles oscillating about
$I_{tot}={\bar{\sigma}}^{2}$ and $X=0$ as shown in Fig. 4B.

While Eq. (9) describes the behavior of common strains, the dynamics of the
antigenic pioneer strains is governed by the equation for $x_{n}$ that in a
continuum limit (suitable for the limit of high mutation rate) reads:

$\frac{d}{dt}x_{n}=\tau_{\rm sw}^{-1}x_{n}-I_{tot}+s\xi(t)$ (10)

The first term on the right hand side represents the rate with which antigenic
pioneer strains enter the population, $\tau_{a}^{-1}$, advancing the nose
fitness by an increment $s$ ($\tau_{a}^{-1}s=\tau_{\rm sw}^{-1}x_{n}$). The
second term on the right hand side of Eq. (10) represents gradual reduction of
susceptibility of the host population, and $\xi(t)$ is a random noise variable
representing the stochasticity of the establishment of new strains. (The
Gaussian white noise $\xi(t)$ is defined statistically by its correlation
function $\langle\xi(t)\xi(0)\rangle=\tau_{a}^{-1}\delta(t)$, see Methods. )

The first term of Eq. (10) captures the apparent instability of the nose: an
advance of the nose to higher $x_{n}$ accelerates its rate of advancement. The
stabilizing factor is the subsequent increase in $I_{tot}$, but to see how
that comes about we must connect Eq. (10) to Eq. (9). The connection is
provided by $\sigma^{2}$ since it is controlled by the emergence of novel
strains, i.e. the dynamics of the “nose” $x_{n}$, which impacts the bulk of
the distribution after a delay $\tau_{\rm sw}$. Based on the analysis detailed
in the Appendix A we approximate

$\sigma^{2}(t)\approx\tau_{\rm sw}^{-1}x_{n}(t-\tau_{\rm sw})$ (11)

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

Figure 5: (A) Simulation results for the average extinction time and the
average speciation time obtained in the clone-based simulation. Extinction
time $\tau_{\rm ext}$, scaled with the sweep time $\tau_{\rm sw}$, increases
with the depth of the genealogy measured by the number of mutations $q$
separating the most antigenically advanced strain from the most common strain.
Also shown is the characteristic time to speciation $\tau_{\rm sp}$, which
increases with the range of cross-inhibition and decreases with $q$. The
crossover of the two time scales defines the transition from transient RQS to
speciation. (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. 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)).

The deterministic limit of the DD system (obtained by omitting the noise term
in Eq. (10)) has two qualitatively different regimes that correspond to the TW
and oscillatory regimes. Small deviations from the steady state with
$\tau_{\rm sw}^{-1}\bar{x}_{n}=\bar{\sigma}^{2}=2\tau_{\rm
sw}^{-2}\log(N_{h}\bar{I})$, are underdamped and oscillate with frequency
$\omega={\bar{\sigma}}={\tau_{\rm sw}}^{-1}{\sqrt{2\log(N_{h}\bar{I})}}$ (as
determined by linearizing Eq. (9)). When $\omega\tau_{\rm sw}<2\pi$, the
delayed feedback via $I_{tot}$ stabilizes the steady state, while in the
opposite regime, 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. As the nonlinearity of Eq. (9) implies a longer period with
increasing amplitude, the system is stabilized at a limit cycle with the
period long enough compared to the feedback delay $\tau_{\rm sw}$. In Appendix
A we derive the threshold of oscillatory instability to lie at
$\log(N_{h}\bar{I}_{osc}s)\approx 8.3$ (leading to limit cycle period
$T\approx 1.5\tau_{\rm sw}$, see Fig. S1 in SI). We also find that the
amplitude of the oscillation $\log(I_{max}/{\bar{I}})$ scales as
$\log(N_{h}\bar{I})$ for large values of the later. This transition defines
quantitatively the boundary between the TW RQS and the Oscillatory RQS regimes
that appear on the phase diagrams in Fig. 3(BC).

The distinction between the TW and Oscillatory RQS is obscured by the
stochasticity of antigenic advance, Eq. (10), which continuously feeds the
underdamped relaxation mode, generating a noisy oscillation with the frequency
$\omega$ defined above. The difference between the two regimes is illustrated
by Fig. 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 Fig. 4A, can
be understood in terms of a non-linear relaxation oscillator. At relatively
low infection prevalence nose fitness $x_{n}$ increases until rising $I_{tot}$
catches up with it (when $I_{tot}=\tau_{\rm sw}^{-1}x_{n}$) driving it down
rapidly. Once this “mini-pandemic” burns out, the population returns to the
low prevalence part of the cycle $I_{tot}<\tau_{\rm sw}^{-1}x_{n}$, when
$x_{n}$ begins to increase again. This relaxation oscillator approximation is
discussed in more detail in the Supplementary Information.

### The rate of extinction

While in the deterministic limit the differential-delay system predicts a
stable steady TW (for $q>q_{ex},\ \bar{I}<\bar{I}_{osc}$) and a limit cycle
(above $\bar{I}_{osc}$), (see SI) fluctuations in the establishment of the
antigenic pioneer strains (Eq. (10)) can lead to stochastic extinction. In
fact, both the TW and Oscillatory RQS (see Fig. 3BC) are transient, subject to
extinction due to a sufficiently large stochastic fluctuation. (Note however
the contrast with the “Extinction” state in Fig. 3BC, where extinction is
deterministic and rapid.) The rate of extinction depends on $q$ and
$\log(N_{h}\bar{I})$ as shown in Figure 5A, which compares results of the
clone-based simulations of the multi-strain SIR model with the fitness class-
based simulation. Although extinction is fluctuation driven, the mechanism of
extinction in the oscillatory state is related closely to the deterministic
dynamics, according to which large amplitude excursion in infection prevalence
can lead to extinction. A large $x_{n}$ advance leads, after a time $\tau_{\rm
sw}$ to a rise in prevalence $I_{tot}$, followed by the rapid fall in the
number of susceptible hosts and hence loss of viral fitness. This turns out to
be the main mode of fluctuation driven extinction as illustrated by Fig. 4C.
One expects extinction to take place when the fluctuation induced deviation of
$x_{n}$ (from its mean) $\delta x\approx s\sqrt{\tau_{\rm ext}/\tau_{a}}$
becomes of the order of the mean $\delta
x\approx{\bar{x}}_{n}\theta(\log(N_{h}\bar{I}))$ with the extinction threshold
$\theta(\log(N_{h}\bar{I}))$ dependent on the shape of the oscillatory limit
cycle (as it depends on the minimum of infection prevalence during the cycle).
This argument suggests $\tau_{\rm ext}/\tau_{\rm sw}\sim
f(\sqrt{q}\theta(\log(N_{h}\bar{I})))$ – a functional relation borne out by
the results of numerical simulations in Fig. 5. We note that the rate increase
in $\tau_{\rm ext}$ with increasing $q$ slows down in the oscillatory regime
and appears to approach a power law dependence $\tau_{\rm ext}/\tau_{\rm
sw}\sim q^{2.5}$ (albeit over a limited accessible range): presently we do not
have analytic understanding of this specific functional form.

Figure 6: Speciation into antigenically distinct lineages. (A) To speciate,
two lineage have to diverge enough to substantially reduce cross-reactivity,
i.e., $T$ needs to be comparable to $d$. Inset: Illustration of the definition
of time to most recent common ancestor $T$ and the time interval $\delta 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^{*}=2\log(N_{h}s^{2})/\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^{*}=\sqrt{2\log(N_{h}s^{2})\sigma_{ag}^{2}}$ with antigenic variance
$\sigma_{ag}^{2}=\sigma^{2}-\sigma^{2}_{\beta}$, where $\sigma^{2}_{\beta}$ is
fitness variance due to deleterious mutations. With this correction,
speciation times agree with the predicted dependence (colored lines).

### 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 populations, i.e., $d/q\gg 1$. In this case, all members of the population
compete against each other for the same susceptible hosts. Conversely, if the
viral population were to split into two sub-populations separated by antigenic
distance greater than the range of cross-inhibition $d$, these sub-population
would no-longer compete for the hosts, becoming effectively distinct viral
“species” that propagate (or fail) independently of each other. Such a split
has for example occurred among influenza B viruses, see Fig. 1.

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

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 B we reformulate this intuition
mathematically as a “first passage”-type problem which shows that $T_{MRCA}$
distribution has an exponential tail which governs the probability of
speciation events. The propensity to speciate depends on the radius of cross-
immunity $d$ and the typically genetic diversity $q$. Fig. 6 shows that the
time to speciation increases approximately exponentially with the ratio $d/q$.
More precisely we found that average simulated speciation time behaves as
$\tau_{\rm sw}^{*}e^{f(d/q^{*})}$ with “effective” $\tau_{\rm
sw}^{*}=\tau_{\rm sw}/(1+\log q/\log(s/m))$ and $q^{*}=q(1+\log q/\log(s/m))$
picking up an additional logarithmic dependence on parameters, the exact
origin of which is beyond our current approximations. This correction
plausibly suggests rapid speciation, $\tau_{\rm sw}^{*}\rightarrow 0$, when
mutation rate become comparable to the selection strength $m/s\rightarrow 1$.

We emphasize that the line in the “phase diagram” (Fig. 3BC) separating the
RQS domain from the Speciation domain refers to the crossover between the
regime where RQS is more likely to go extinct before speciating and vice
versa. Fig. 5A shows average speciation and extinction times for different
parameters – the intersection of the two times defines the boundary in the
phase diagram.

## Discussion

The epidemiological and evolutionary dynamics of human RNA viruses show a
number of qualitatively distinct patterns. While agents of classical childhood
diseases like measles or mumps virus show little antigenic evolution, others
viruses like dengue- or norovirus exist in distinct serotypes, while seasonal
influenza viruses undergo continuous antigenic evolution enabling the 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. The integration of these
two different crucial 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 Fig. 3B&C.

The phase diagram shows different combinations of key parameters that 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(N_{h}s^{2})$, the ratio of mutational effects
and mutation rate $\log(s/m)$, and the cross-immunity range $d$. Large $d$
prevents speciation, while rapid mutation and large population sizes
facilitate speciation. 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 can be very long
lived.

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 didn’t 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. We know of four
seasonal influenza A lineages: 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) around 1983 (Rota  _et al._ , 1990).
Phylogenetic trees of A/H3N2 and the influenza B lineages are shown in Fig. 1.
In addition, the figure shows diversity of five lineages as measured by the
instantaneous $T_{MRCA}$ through time.

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 Fig. 1. 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 loose immunity and $T$ 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.

From the phase diagram, we found the most relevant parameters that determine
the fate of a pathogen are the antigenic diversity $q$ and the range of cross-
immunity $d$. A previous study by Koelle and Rasmussen (2015) has implicated
deleterious mutation load as a cause of spindly phylogenies. Our model can
readily incorporate the effect of deleterious mutations affecting transmission
coefficient $\beta$ by $\delta\beta$ and subsequent compensatory mutations
that increase $\beta$. Such a modification is expected to reduce the average
$\beta$ and reduce the selection coefficient of antigenic mutations, which in
turn reduces the fitness variance $\sigma^{2}$, as derived in Supplementary
Information. 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 Fig. 6C.

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 $\beta$) 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.

Symbol | meaning | comment/range | numerical range
---|---|---|---
$I_{a}$ | number of individuals infected with strain $a$ |  |
$S_{a}$ | weighted fraction of individuals susceptible to strain $a$ | $\sim 0.7$ |
$d$ | cross-immunity range |  $\sim 10$y | $[50,200]$
$\beta=\nu R_{0}$ | transmission rate | $\sim 3$ | 2
$\nu$ | recovery rate |  $\sim 1$w | 1
$\gamma$ | birth/death rate of people | $\sim 0.01y^{-1}$ | 0
$K_{ab}=e^{-|b-a|/d}$ | cross-immunity of strains $a$, $b$ |  $\sim 20$ mutations |
$\tau_{\rm sw}$ | coalescent time scale/sweep time | $2-6y$ |
$T$ | $T_{MRCA}$ |  $\sim 3-10$y |
$\delta T$ | $T_{MRCA}$ |  $\sim 2-6$y |
$s$ | Selection coefficient | $\sim 0.03w^{-1}$ | $[0.003,0.05]$
$\Delta$ | antigenic effect size |  $\sim 0.5$ log2 titer approx $0.1d$ |
$m$ | mutation rate | beneficial $10^{-3}$ per week and genome | $[10^{-7},10^{-3}]$
$I_{tot}=\sum_{a}I_{i}$ | total prevalence | $0.005$ |
$\bar{I}$ | average prevalence |  |
$N_{h}$ | host population | human population $10^{10}$ | $[10^{6},10^{12}]$
Table 1: Relevant quantities of influenza virus and parameters in multi-strain
SIR model.

## Methods

###  .1 Clone-based simulations

We simulate the original model on a genealogical tree in two phases: one on
the deterministic SIR-type epidemics and one on the stochastic mutation
introducing new strains. In each time step $\Delta t<1$, we apply the mid-
point method to advance the epidemic equations Eqs. (1,2,3). We then generate
a random number uniformly sampled between zero and one for each surviving
strain with $N_{h}I_{a}\geq 1$. If the random number is smaller than
$mN_{h}I_{a}\Delta t$ for strain $a$, we append a new strain $b$ as a
descendent to $a$. The susceptibility to strain $b$ is related to
susceptibility to strain $a$ via $S_{b}=(S_{a})^{e^{-1/d}}$. In most of the
simulations, the transmissibility of different strains is held constant
$\beta$. Otherwise we allow for a strain specified transmissibility that is
drawn from its parent $\beta_{b}=\beta_{a}-\delta\beta$ with $\delta\beta>0$
for the deleterious effect of antigenic mutations and $\beta_{b}=\beta_{\rm
max}$ if the mutation is compensatory. The new strain grows deterministically
only if $\beta_{b}S_{b}>1$.

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

###  .2 Fitness class-based simulations

The stability of the RQS and the extinction dynamics is fully captured by the
traveling wave equations (6). We simulate the traveling wave by gridding the
fitness space $x$ into bins of step size $s$ around zero. The infections of
different strains correspond to a natural number in each bin $x_{i}$. At each
time step, the population in each bin $I_{i}$ updates to a number sampled from
the Poisson distribution with parameter
$\lambda_{i}=I_{i}(1+(x_{i}-\bar{x})\Delta t)$ determined by mean fitness
$x_{i}$ and a dynamic mean fitness $\bar{x}$, which increases by $\Delta
tI_{tot}/N_{h}$, where $I_{tot}$ is the total population summed over all bins
and $N_{h}$ is the parameter giving the host population. When $\bar{x}$
becomes larger than one bin size $s$, we shift the all populations to left by
one bin and reset $\bar{x}$ to $0$, a trick to keep only a finite number of
bins in the simulation. At the same time, antigenic mutation is represented by
moving a the mutated fraction in each bin to the adjacent bin on the right.
The faction is determined by a random number drawn from the Poisson
distribution with the mean $mI_{i}\Delta t$. The typical ranges of the three
parameters $s$, $m$, and $N_{h}$ are according to the parameters in the
genealogical simulation, as documented also in Table. 1.

###  .3 Stochastic differential-delay simulation

To simulate the differential delay equations Eqs. (9,10,11), we discretize
time in increments of $\Delta t=\tau_{\rm sw}/k$ and update the dynamical
variables $\chi_{i}=x_{n}(t_{i})$ and $\eta_{i}=I_{tot}(t_{i})$ via the simple
Euler scheme:

$\displaystyle\chi_{i+1}$ $\displaystyle=$ $\displaystyle\chi_{i}+\Delta
t(\chi_{i}-\eta_{i})+\frac{\chi_{i}}{qs}\sqrt{\Delta t}\xi_{i};$ (12)
$\displaystyle\eta_{i+1}$ $\displaystyle=$
$\displaystyle\bar{I}\exp\left(\tau_{\rm sw}\chi_{i-k}-\frac{\tau_{\rm
sw}^{2}}{k}\sum_{j=0}^{k}j\eta_{i-j}\right),$ (13)

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

###  .4 Influenza phylogenies

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._ , 2017) and
TreeTime (Sagulenko  _et al._ , 2018).

## Appendix

###  .1 Differential-delay approximation of RQS dynamics

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

$G(\lambda,t)=\sum_{i}I_{i}(t)e^{\lambda x_{i}(t)}$ (14)

where $x_{i}(t)=x_{n}(t_{i})-\int_{t_{i}}^{t}dt^{\prime}I_{t}(t^{\prime})$ is
the current fitness of the pioneer strain that first appeared at time $t_{i}$
and $I_{i}(t)$ is the fraction of the hosts infected by it:

$I_{i}(t)=N^{-1}e^{\int_{t_{i}}^{t}dt^{\prime}x_{t}(t^{\prime})}=N^{-1}e^{x_{n}(t_{i})(t-t_{i})-\int_{0}^{t-t_{i}}dt^{\prime}t^{\prime}I_{tot}(t-t^{\prime})}$
(15)

We next take a coarse grained view of pioneer strain establishment replacing
the sum in Eq. (13) by an integral over initial times $t_{i}\rightarrow
t-\tau$

$G(\lambda,t)=\int_{0}^{\infty}\frac{d\tau
N^{-1}}{\tau_{a}(t-\tau)}e^{(\tau+\lambda)x_{n}(t-\tau)-\int_{0}^{\tau}dt^{\prime}(t^{\prime}+\lambda)I_{tot}(t-t^{\prime})}$
(16)

Let us evaluate the integral in the saddle approximation which is dominated by
$\tau=\tau^{*}$ corresponding to the maximum in the exponential

$\tau^{*}+\lambda=\frac{x_{n}(t-\tau^{*})}{x^{\prime}_{n}(t-\tau^{*})+I_{tot}(t-\tau^{*})}\approx\tau_{\rm
sw}$ (17)

where we have used the deterministic limit of Eq (9). To simplify presentation
we shall ignore the time dependence of $\tau_{\rm sw}=s^{-1}\log(x_{n}/m)$
replacing $x_{n}(t-\tau^{*})$ in the logarithm by the time average
$\bar{x_{n}}$.

In the saddle approximation we then have

$\log NG(\lambda,t)\approx x_{n}(t-\tau_{\rm sw}+\lambda)\tau_{\rm
sw}-\int_{0}^{\tau_{\rm
sw}-\lambda}dt^{\prime}(t^{\prime}+\lambda)I_{tot}(t-t^{\prime})$ (18)

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

We can now estimate fitness mean

$\bar{x}(t)=\frac{d}{d\lambda}\log G(\lambda,t)|_{\lambda=0}\\\ =\tau_{\rm
sw}[x^{\prime}_{n}(t-\tau_{\rm sw})+I_{tot}(t-\tau_{\rm
sw})]-\int_{0}^{\tau_{\rm sw}}dt^{\prime}I_{tot}(t-t^{\prime})\\\
=x_{n}(t-\tau_{\rm sw})-\int_{0}^{\tau_{\rm
sw}}dt^{\prime}I_{t}(t-t^{\prime})$ (19)

and variance

$\sigma^{2}(t)=\frac{d^{2}}{d\lambda^{2}}\log G(\lambda,t)|_{\lambda=0}\\\
=\tau_{\rm sw}[x^{\prime\prime}_{n}(t-\tau_{\rm
sw})+I_{tot}^{\prime}(t-\tau_{\rm sw})]+I_{tot}(t-\tau_{\rm sw})$ (20)

Eq. (20) involves the second derivative $x^{\prime\prime}_{n}$ and expect
fluctuations in the establishment of new lineages (which contribute to
$x^{\prime}_{n}$) to be quite important. Yet we can get useful insight by
continuing to use the deterministic approximation to $x_{n}$ dynamics, in
which case we arrive at simple delay relation between the variance and $x_{n}$

$\sigma^{2}(t)=\tau_{\rm sw}^{-1}x_{n}(t-\tau_{\rm sw})$ (21)

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

$\frac{d}{dt}{\bar{x}}=x^{\prime}_{n}(t-\tau_{\rm sw})+I_{t}(t-\tau_{\rm
sw})-I_{t}(t)\\\ =\sigma^{2}(t)-I_{t}(t)$ (22)

Combining Eqs. (9,10,21) we arrive at the deterministic dynamical system
approximating coupled “ecological” SIR dynamics with the evolutionary dynamics
of antigenic innovation due to the pioneer strains.

$\frac{d^{2}}{dt^{2}}\log I(t)=\tau_{\rm sw}^{-1}x_{n}(t-\tau_{\rm
sw})-I_{tot}(t)$ (23) $\frac{d}{dt}x_{n}(t)=\tau_{\rm
sw}^{-1}x_{n}(t)-I_{tot}(t)$ (24)

This system admits a family of fixed points of the form $\tau_{\rm
sw}I_{tot}=x_{n}={\bar{x}}_{n}$, but as we show in the SI, the corresponding
steady TW states are not always stable giving rise to limit cycle oscillations
or leading to rapid extinction. Self-consistency condition relating $x_{n}$
and $I_{tot}$ for the steady traveling wave is readily generalized to limit
cycle states. Integrating the differential-delay system over one cycle yields
$\langle x_{n}\rangle=\tau_{\rm sw}\langle I\rangle$. An additional relation
is provided by integrating $\log NG(0,t)$ over the cycle:

$\langle\log N_{h}I_{tot}\rangle=\frac{\tau_{\rm sw}^{2}}{2}\langle
I_{tot}\rangle$ (25)

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

###  .2 Speciation rate as a stochastic “First Passage” problem.

Speciation occurs when two most distant clades persist to the antigenic
independence. This persistence problem can be formulated as a first passage
problem by including the second ”nose” in the TW approximation.

Figure 7: 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
$x_{1}$ and $x_{2}$. Branch 1 is the fitter one $x_{1}>x_{2}$. The antigenic
distances from the cross-immune bulk to the noses of the two branches are
$d_{1}$ and $d_{2}$. The Gaussian profile in fitness is illustrated in blue.
Right: The fitness difference between the two branches $x_{1}-x_{2}$ 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=x_{n}$.

We consider the births of two pioneer strains at time $t=0$, as illustrated in
Fig. 7. The descendants of the two strains forming two branches 1 and 2
diverge in the antigenic space as they persist in time. Suppose that at time
$t$, the nose of branch 1 is at fitness $x_{1}$, and the nose of branch 2 is
at $x_{2}$. Before the sweep time $t<\tau_{\rm sw}$, the cross-immunity grows
mainly from the prevalent strains in the common ancestors of the two branches,

$\frac{d}{dt}x_{i}=\tau_{\rm sw}^{-1}x_{i}-I_{\rm tot}+s\xi_{i}.\qquad i=1,2$
(26)

Later when $t>\tau_{\rm sw}$, the infection bulk splits and moves on to the
two branches. As the antigenic distances from the noses to the infection bulks
on different branches are different, cross-immunity effects to different noses
grow in different rates,

$\displaystyle\frac{d}{dt}x_{1}$ $\displaystyle=$ $\displaystyle\tau_{\rm
sw}^{-1}x_{1}-I_{1}e^{-d_{11}/d}-I_{2}e^{-d_{21}/d}+s\xi_{1};$ (27)
$\displaystyle\frac{d}{dt}x_{2}$ $\displaystyle=$ $\displaystyle\tau_{\rm
sw}^{-1}x_{2}-I_{1}e^{-d_{12}/d}-I_{2}e^{-d_{22}/d}+s\xi_{2},$

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

In this limit, the survival probability of the less fit nose maps to a first
passage problem in the random walk of relative fitness
$\zeta\equiv(x_{1}-x_{2})/x_{n}$. As illustrated in Fig. 7, an establishment
of nose 1 is a positive step of $\delta\zeta=s/x_{n}$, while an establishment
of nose 2 ends up into a backward step of the same size. As the mutations
arrive in characteristic times $\tau_{1}$ and $\tau_{2}$ depending on the nose
fitnesses, in the continuum limit, we have

$\frac{d}{dt}\zeta=\tau_{\rm sw}^{-1}\zeta+\frac{s}{x_{n}}\xi,$ (28)

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

The system can be solved for the probability density distribution
$\rho(\zeta,t)$, which obeys

$\partial_{t}\rho(\zeta,t)=-\partial_{\zeta}[v(\zeta)\rho(\zeta,t)]+\partial_{\zeta}^{2}[D(\zeta)\rho(\zeta,t)],$
(29)

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

$v(\zeta)=\frac{1}{\tau_{\rm sw}}\zeta;\qquad
D(\zeta)\approx\frac{1}{q\tau_{\rm sw}}.$ (30)

Solving with boundary and initial conditions,

$\displaystyle\partial_{\zeta}\rho(\zeta,t)|_{\zeta=0}$
$\displaystyle=0;\quad\rho(\zeta=1,t)=0;$ (31) $\displaystyle\rho(\zeta,t=0)$
$\displaystyle=\delta(\zeta),$

we have

$\rho(\zeta,t)=\sum_{n=1}^{\infty}e^{-\lambda_{n}t/\tau_{\rm
sw}}c_{n1}F_{1}(\frac{1-\lambda_{n}}{2},\frac{1}{2},\frac{q}{2}\zeta^{2}),$
(32)

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

$P(T>t)\approx ce^{-t/\tau_{\rm sw}}.$ (33)

The typical time interval between successive strains scales as
$\tau_{a}=\tau_{\rm sw}/q$. So the probability of the branch depth being
deeper than $a$ is then

$P(D>a)\approx e^{-a\tau_{a}/\tau_{\rm sw}}=e^{-a/q}.$ (34)

Recalling that the speciation, or escape of the cross-immunity occurs when the
branch depth is larger than $d$, we find the probability of a successful
branching $p_{1}$ is proportional to $e^{-d/q}$.

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

$P_{\rm nsp}(t)=(1-p_{1})^{t/\tau_{a}}=e^{-t/\tau_{\rm sp}},$ (35)

where the waiting time for speciation event is

$\tau_{\rm sp}\propto\frac{\tau_{\rm sw}}{q}e^{d/q},$ (36)

as numerically verified in Fig. 6.

Figure S1: 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(N_{h}\bar{I}s)_{osc}=8.3$, indicated by the dashed line. When the noise
is properly considered as in Eq. (13), 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.

## Supplementary Information

Figure S2: ”Phase diagram” in the clone-based simulation. The data points in
corresponding phases are determined by following criteria. The extinction
phase (red): $\tau_{\rm ext}<10\tau_{\rm sw}$; The speciation phase (blue):
$\tau_{\rm sp}<\tau_{\rm ext}$; The transient RQS regime (yellow and green):
otherwise. Color in the transient phase labels the amplitude of prevalence
${\rm var}(\log I_{tot})$, increasing from yellow to green.

### Stability analysis of the differential-delay approximation

It is natural to measure time in the units of $\tau_{\rm sw}$: $t=\tau_{\rm
sw}\zeta$, $x_{n}=\tau_{\rm sw}^{-1}\chi$ and define $u=\tau_{\rm sw}^{2}I$
and

$\frac{d^{2}}{d\zeta^{2}}u(\zeta)=\chi(\zeta-1)-e^{u(\zeta)}$ (S1)
$\frac{d}{d\zeta}\chi(\zeta)=\chi(\zeta)-e^{u(\zeta)}$ (S2)

This system has a one parameter family of fixed points $\chi=\bar{\chi},\
u=\log\bar{\chi}$. To analyze fixed point stability we linearize and Laplace
transform, yielding

$z^{2}\delta{\hat{u}}(z)=e^{-z}\delta{\hat{\chi}}(z)-\bar{\chi}\delta{\hat{u}}(z)+z\delta
u(0)+\delta u^{\prime}(0)$ (S3)
$z\delta{\hat{\chi}}(z)=\delta{\hat{\chi}}(za)-\bar{\chi}\delta{\hat{u}}(z)+z\delta\chi(0)$
(S4)

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

$z=1+\bar{\chi}(1-z-e^{-z})/z^{2}$ (S5)

Fixed point - and hence steady RQS - stability requires $\Re(z)<0$ which is
found for $2<\bar{\chi}<\bar{\chi}_{c}$. For $\bar{\chi}>2.845$ one finds
$\Im(z)\neq 0$ corresponding to the onset of oscillatory relaxation which
turns into a limit cycle for $\bar{\chi}>\bar{\chi}_{c}\approx 16.6$. The
period of the limit cycle is well approximated by $\Im(z)$, as the dashed line
shown in the bottom panel of Fig. S1.

### Stochastic form of the differential-delay approximation

A sensible stochastic generalization is obtained by the stochastic
approximation for the “nose” dynamics (Eq. (10))

$\frac{d}{dt}x_{n}=\tau_{\rm sw}^{-1}x_{n}-I_{t}(t)+s\xi(t),$ (S6)

combined with Eq. (18) at $\lambda=0$

$\log I(t)=\tau_{\rm sw}x_{n}(t-\tau_{\rm sw})-\int_{0}^{\tau_{\rm
sw}}dt^{\prime}t^{\prime}I_{t}(t-t^{\prime}).$ (S7)

Note that in this derivation we have avoided the need for explicitly
approximating $\sigma^{2}$! (We have also neglected the effect of fluctuations
arising from the logarithmic correction term effectively replacing it by its
average value.)

### Effect of mutations in infectivity

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

The fitness of the nose of the distribution obeys

$\frac{dx_{n}}{dt}=s_{d}(\beta)\tau_{a}^{-1}+\Delta_{b}\tau_{b}^{-1}-I$ (S8)

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

$I_{tot}=s(\beta_{*})\tau_{ag}^{-1}=\frac{s(\beta_{*})x_{n}}{\log\frac{x_{n}}{m}}$
(S9)

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

The sweep time, $\tau_{\rm sw}$, upon which the fitness of the former pioneer
strain comes down to the mean fitness) and the nose fitness, $x_{n}$, retain
the TW form

$\tau_{\rm
sw}=\frac{x_{n}}{I_{tot}}=\frac{s(\beta_{*})}{\log(\frac{x_{n}}{m})}$ (S10)

Following TW approximation to estimate infection prevalence
$\sqrt{I_{tot}}\sim N^{-1}\exp(x_{n}\tau_{\rm sw}/2)$ as before one finds

$x_{n}=2\tau_{\rm sw}^{-1}\log
CN_{h}=2s(\beta_{*})\frac{\log[N_{h}s^{2}/\log(\frac{x_{n}}{m})]}{\log(\frac{x_{n}}{m})}$
(S11)

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

$\sigma^{2}=\tau_{\rm
sw}(s_{d}^{2}\tau_{a}^{-1}+\Delta_{b}^{2}\tau_{b}^{-1})=x_{n}[s-2\Delta_{d}+\frac{\Delta_{d}^{2}}{s}+\frac{\Delta_{d}\Delta_{b}}{s}]$
(S12)

but under conditions of $\Delta_{d},\ \Delta_{b}\ll s(\beta_{*})$ total
variance would also be decreasing.

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

$q_{ag}=\tau_{\rm sw}\tau_{ag}^{-1}=\frac{x_{n}}{s(\beta_{*})}=2\frac{\log
N_{h}s^{2}(\beta_{*})c}{\log\frac{x_{n}}{m}}$ (S13)

which is basically unchanged in the presence of infectivity mutations except
for the expected reduction in the magnitude of $s^{2}$ factor inside the
logarithm. Therefore, speciation rate would be reduced, but rather weakly, via
a contribution subleading in $o(\log N_{h})$

## References

  * Andreasen  _et al._ (1997) Andreasen, V., J. Lin, and S. A. Levin (1997), Journal of mathematical biology 35 (7), 825.
  * Andreasen and Sasaki (2006) Andreasen, V., and A. Sasaki (2006), Theoretical Population Biology 70 (2), 164.
  * Bedford  _et al._ (2012) Bedford, T., A. Rambaut, and M. Pascual (2012), BMC biology 10 (1), 38.
  * Bedford  _et al._ (2014) Bedford, T., M. A. Suchard, P. Lemey, G. Dudas, V. Gregory, A. J. Hay, J. W. McCauley, C. A. Russell, D. J. Smith, and A. Rambaut (2014), Elife 3, e01914.
  * Bhatt  _et al._ (2011) Bhatt, S., E. C. Holmes, and O. G. Pybus (2011), Molecular Biology and Evolution 28 (9), 2443.
  * Desai and Fisher (2007) Desai, M. M., and D. S. Fisher (2007), Genetics 176 (3), 1759.
  * Desai  _et al._ (2013) Desai, M. M., A. M. Walczak, and D. S. Fisher (2013), Genetics 193 (2), 565.
  * Ferguson  _et al._ (2003) Ferguson, N. M., A. P. Galvani, and R. M. Bush (2003), Nature 422 (6930), 428.
  * Fisher (2013) Fisher, D. S. (2013), Journal of Statistical Mechanics: Theory and Experiment 2013 (01), P01011.
  * Gog and Grenfell (2002) Gog, J. R., and B. T. Grenfell (2002), Proceedings of the National Academy of Sciences 99 (26), 17209.
  * Gomulkiewicz and Holt (1995) Gomulkiewicz, R., and R. D. Holt (1995), Evolution 49 (1), 201.
  * Grenfell  _et al._ (2004) Grenfell, B. T., O. G. Pybus, J. R. Gog, J. L. Wood, J. M. Daly, J. A. Mumford, and E. C. Holmes (2004), science 303 (5656), 327.
  * Hadfield  _et al._ (2017) Hadfield, J., C. Megill, S. M. Bell, J. Huddleston, B. Potter, C. Callender, P. Sagulenko, T. Bedford, and R. A. Neher (2017), bioRxiv , 224048.
  * Hallatschek (2011) Hallatschek, O. (2011), Proceedings of the National Academy of Sciences 108 (5), 1783.
  * Katoh  _et al._ (2002) Katoh, K., K. Misawa, K.-i. Kuma, and T. Miyata (2002), Nucleic Acids Research 30 (14), 3059.
  * Koel  _et al._ (2013) Koel, B. F., D. F. Burke, T. M. Bestebroer, S. v. d. Vliet, G. C. M. Zondag, G. Vervaet, E. Skepner, N. S. Lewis, M. I. J. Spronken, C. A. Russell, M. Y. Eropkin, A. C. Hurt, I. G. Barr, J. C. d. Jong, G. F. Rimmelzwaan, A. D. M. E. Osterhaus, R. A. M. Fouchier, and D. J. Smith (2013), Science 342 (6161), 976.
  * Koelle and Rasmussen (2015) Koelle, K., and D. A. Rasmussen (2015), eLife 4, e07361.
  * London and Yorke (1973) London, W. P., and J. A. Yorke (1973), American journal of epidemiology 98 (6), 453.
  * Neher (2013) Neher, R. A. (2013), Annual Review of Ecology, Evolution, and Systematics 44 (1), 195.
  * Neher  _et al._ (2016) Neher, R. A., T. Bedford, R. S. Daniels, C. A. Russell, and B. I. Shraiman (2016), Proceedings of the National Academy of Sciences 113 (12), E1701, http://www.pnas.org/content/113/12/E1701.full.pdf .
  * Neher and Hallatschek (2013) Neher, R. A., and O. Hallatschek (2013), Proceedings of the National Academy of Sciences of the United States of America 110 (2), 437.
  * Nguyen  _et al._ (2015) Nguyen, L.-T., H. A. Schmidt, A. von Haeseler, and B. Q. Minh (2015), Molecular Biology and Evolution 32 (1), 268.
  * O?Reilly  _et al._ (2018) O?Reilly, K. M., R. Lowe, W. J. Edmunds, P. Mayaud, A. Kucharski, R. M. Eggo, S. Funk, D. Bhatia, K. Khan, M. U. G. Kraemer, A. Wilder-Smith, L. C. Rodrigues, P. Brasil, E. Massad, T. Jaenisch, S. Cauchemez, O. J. Brady, and L. Yakob (2018), BMC Medicine 16 (1), 180.
  * Perelson and Oster (1979) Perelson, A. S., and G. F. Oster (1979), Journal of Theoretical Biology 81 (4), 645.
  * Petrova and Russell (2018) Petrova, V. N., and C. A. Russell (2018), Nature Reviews Microbiology 16 (1), 47.
  * Rota  _et al._ (1990) Rota, P. A., T. R. Wallis, M. W. Harmon, J. S. Rota, A. P. Kendal, and K. Nerome (1990), Virology 175 (1), 59.
  * Rouzine and Rozhnova (2018) Rouzine, I. M., and G. Rozhnova (2018), PLOS Pathogens 14 (9), e1007291.
  * Rouzine  _et al._ (2003) Rouzine, I. M., J. Wakeley, and J. M. Coffin (2003), Proceedings of the National Academy of Sciences 100 (2), 587.
  * Sagulenko  _et al._ (2018) Sagulenko, P., V. Puller, and R. A. Neher (2018), Virus Evolution 4 (1), 10.1093/ve/vex042.
  * Smith  _et al._ (2004) Smith, D. J., A. S. Lapedes, J. C. de Jong, T. M. Bestebroer, G. F. Rimmelzwaan, A. D. Osterhaus, and R. A. Fouchier (2004), Science .
  * Tria  _et al._ (2005) Tria, F., M. L ssig, L. Peliti, and S. Franz (2005), Journal of Statistical Mechanics: Theory and Experiment 2005 (07), P07008.
  * Tsimring  _et al._ (1996) Tsimring, L. S., H. Levine, and D. A. Kessler (1996), Physical review letters 76 (23), 4440.
  * Van Valen (1973) Van Valen, L. (1973), Evol Theory 1, 1.


B}_{0,1}{\mathrm{L}_{3}^{+}}={\mathrm{Y}_{2}^{-}}+f_{\mathcal{W}_{4}^{-}}\,,$
(B.43a) where
${\mathrm{Y}_{2}^{-}}(x):=\begin{bmatrix}\dfrac{p_{1}^{[1]}\big{(}a_{2}^{[2]}{\mathtt{c}}_{\mathtt{h}}(3+{\mathtt{c}}_{\mathtt{h}}^{4})-2(3+5{\mathtt{c}}_{\mathtt{h}}^{4})p_{2}^{[2]}\big{)}}{16{\mathtt{c}}_{\mathtt{h}}^{9/2}}\sin(2x)\\\\[11.38109pt]
-\dfrac{3(1+3{\mathtt{c}}_{\mathtt{h}}^{4})p_{1}^{[1]}(a_{2}^{[2]}{\mathtt{c}}_{\mathtt{h}}-2p_{2}^{[2]})+a_{1}^{[1]}{\mathtt{c}}_{\mathtt{h}}\big{(}-a_{2}^{[2]}{\mathtt{c}}_{\mathtt{h}}(3+{\mathtt{c}}_{\mathtt{h}}^{4})+2(3+5{\mathtt{c}}_{\mathtt{h}}^{4})p_{2}^{[2]}\big{)}}{32{\mathtt{c}}_{\mathtt{h}}^{11/2}}\cos(2x)\end{bmatrix}\,.$
(B.43b) In view of (B.3) and Lemma B.2, one obtains by inspection
$\frac{\lambda
m_{1}^{+}}{\lambda^{2}+4{\mathtt{c}}_{\mathtt{h}}^{2}}(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{0,1}f_{-1}^{-}=\mathcal{O}(\lambda)\,,\quad(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{0,1}\mathcal{O}_{\mathcal{W}_{3}}(\lambda)=\mathcal{O}(\lambda)\,.$
(B.43c)

Then (B.42), Lemmata B.4, B.2 and (B.43) give

${\footnotesize\begin{matrix}(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{0,1}(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{0,2}f_{1}^{+}\end{matrix}}=\tau_{1}^{+}{\mathrm{B}_{2}^{+}}+\frac{\ell_{1}^{+}}{4\text{${\mathtt{c}}_{\mathtt{h}}^{2}$}}{\mathrm{A}_{-2}^{+}}+\mathscr{L}_{0,0}^{\,\rotatebox[origin={c}]{90.0}{$\perp$}}({\mathrm{Y}_{2}^{-}}+f_{\mathcal{W}_{4}^{-}})+\mathcal{O}(\lambda^{-1}:\lambda)\,,$

and, by (B.7), the term $\mathrm{V\\!I\\!I\\!I}$ in (B.34) is

$\mathrm{V\\!I\\!I\\!I}=\tau_{1}^{+}{\mathrm{B}_{2}^{+}}+\frac{\ell_{1}^{+}}{4{\mathtt{c}}_{\mathtt{h}}^{2}}{\mathrm{A}_{-2}^{+}}+\mathscr{L}_{0,0}^{\,\rotatebox[origin={c}]{90.0}{$\perp$}}{\mathrm{Y}_{2}^{-}}+\widetilde{f}_{\mathcal{W}_{4}^{+}}\,,$
(B.44)

where

$\mathscr{L}_{0,0}^{\,\rotatebox[origin={c}]{90.0}{$\perp$}}{\mathrm{Y}_{2}^{-}}(x)=\begin{bmatrix}\kappa_{3}\cos(2x)\\\\[11.38109pt]
\varsigma_{3}\sin(2x)\end{bmatrix}\,.$ (B.45)

with

$\displaystyle\kappa_{3}:$
$\displaystyle=\dfrac{a_{1}^{[1]}{\mathtt{c}}_{\mathtt{h}}\big{(}a_{2}^{[2]}{\mathtt{c}}_{\mathtt{h}}({\mathtt{c}}_{\mathtt{h}}^{4}+3)-2(5{\mathtt{c}}_{\mathtt{h}}^{4}+3)p_{2}^{[2]}\big{)}-a_{2}^{[2]}{\mathtt{c}}_{\mathtt{h}}({\mathtt{c}}_{\mathtt{h}}^{8}+13{\mathtt{c}}_{\mathtt{h}}^{4}+6)p_{1}^{[1]}+2({\mathtt{c}}_{\mathtt{h}}^{4}+3)(5{\mathtt{c}}_{\mathtt{h}}^{4}+2)p_{1}^{[1]}p_{2}^{[2]}}{32{\mathtt{c}}_{\mathtt{h}}^{19/2}}\,,$
$\displaystyle\varsigma_{3}:$
$\displaystyle=\dfrac{({\mathtt{c}}_{\mathtt{h}}^{4}+1)\big{[}a_{1}^{[1]}{\mathtt{c}}_{\mathtt{h}}\big{(}a_{2}^{[2]}{\mathtt{c}}_{\mathtt{h}}({\mathtt{c}}_{\mathtt{h}}^{4}+3)-2(5{\mathtt{c}}_{\mathtt{h}}^{4}+3)p_{2}^{[2]}\big{)}-2a_{2}^{[2]}{\mathtt{c}}_{\mathtt{h}}(5{\mathtt{c}}_{\mathtt{h}}^{4}+3)p_{1}^{[1]}+4(7{\mathtt{c}}_{\mathtt{h}}^{4}+3)p_{1}^{[1]}p_{2}^{[2]}\big{]}}{64{\mathtt{c}}_{\mathtt{h}}^{21/2}}\,.$

$\mathrm{I\\!X}\big{)}$ By (LABEL:2104_1749) we get

$\displaystyle[(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{0,1}]^{3}f_{1}^{+}=\frac{\zeta_{2}^{+}}{\lambda}(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{0,1}f_{1}^{-}$ (B.46a)
$\displaystyle+\frac{\alpha_{2}^{+}}{\lambda^{2}+4{\mathtt{c}}_{\mathtt{h}}^{2}}(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{0,1}f_{-1}^{+}+\zeta_{3}^{+}(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{0,1}f_{1}^{+}+(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{0,1}{\mathrm{A}_{3}^{+}}$ $\displaystyle\begin{aligned}
&+\lambda(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{0,1}\left(f_{\mathcal{W}_{1}^{-}}+f_{\mathcal{W}_{3}^{-}}\right)+\lambda^{2}(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{0,1}f_{\mathcal{W}_{1}^{+}}\\\
&+(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{0,1}\mathcal{O}_{\mathcal{W}_{3}}(\lambda^{2})+(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{0,1}\mathcal{O}(\lambda^{3})\,.\end{aligned}$ (B.46b)

In view of (B.3) and Lemma B.2 the terms in the two lines in (LABEL:2605_1113)
are

$\displaystyle\lambda(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{0,1}\left(f_{\mathcal{W}_{1}^{-}}+f_{\mathcal{W}_{3}^{-}}\right)=\mathcal{O}(\lambda)\,,\
\ \lambda^{2}(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{0,1}f_{\mathcal{W}_{1}^{+}}=\mathcal{O}(\lambda)\,,\ $ $\displaystyle\
(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{0,1}\mathcal{O}_{\mathcal{W}_{3}}(\lambda^{2})=\mathcal{O}(\lambda^{2})\,,\quad(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{0,1}\mathcal{O}(\lambda^{3})=\mathcal{O}(\lambda)\,.$

The remaining terms in (B.46), again by Lemma B.4, are

$\displaystyle[(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{0,1}]^{3}f_{1}^{+}=\zeta_{2}^{+}{\mathrm{B}_{2}^{+}}+\frac{\alpha_{2}^{+}}{4\text{${\mathtt{c}}_{\mathtt{h}}^{2}$}}{\mathrm{A}_{-2}^{+}}+\zeta_{3}^{+}\mathrm{A}_{2}^{+}+\mathscr{L}_{0,0}^{\,\rotatebox[origin={c}]{90.0}{$\perp$}}\mathcal{J}{\cal
B}_{0,1}{\mathrm{A}_{3}^{+}}+\mathcal{O}(\lambda^{-1}:\lambda)\,.$ (B.47)

By applying the operators ${\cal B}_{0,1}$ in (3.8b) and $\mathcal{J}$ in
(1.5) to the vector ${\mathrm{A}_{3}^{+}}$ in (B.11) one gets

$\mathcal{J}{\cal
B}_{0,1}{\mathrm{A}_{3}^{+}}={\mathrm{Z}_{2}^{-}}+f_{\mathcal{W}_{4}^{-}}$
(B.48)

where

${\mathrm{Z}_{2}^{-}}(x):=-\begin{bmatrix}\frac{p_{1}^{[1]}\left((a_{1}^{[1]})^{2}\left({\mathtt{c}}_{\mathtt{h}}^{4}+3\right){\mathtt{c}}_{\mathtt{h}}^{2}-2a_{1}^{[1]}p_{1}^{[1]}\left({\mathtt{c}}_{\mathtt{h}}^{8}+9{\mathtt{c}}_{\mathtt{h}}^{4}+6\right){\mathtt{c}}_{\mathtt{h}}+(p_{1}^{[1]})^{2}\left(11{\mathtt{c}}_{\mathtt{h}}^{8}+29{\mathtt{c}}_{\mathtt{h}}^{4}+12\right)\right)}{32{\mathtt{c}}_{\mathtt{h}}^{19/2}}\sin(2x)\\\\[11.38109pt]
\frac{(a_{1}^{[1]}{\mathtt{c}}_{\mathtt{h}}-p_{1}^{[1]})\left((a_{1}^{[1]})^{2}\left({\mathtt{c}}_{\mathtt{h}}^{4}+3\right){\mathtt{c}}_{\mathtt{h}}^{2}-2a_{1}^{[1]}p_{1}^{[1]}\left({\mathtt{c}}_{\mathtt{h}}^{8}+13{\mathtt{c}}_{\mathtt{h}}^{4}+6\right){\mathtt{c}}_{\mathtt{h}}+3(p_{1}^{[1]})^{2}\left(9{\mathtt{c}}_{\mathtt{h}}^{8}+15{\mathtt{c}}_{\mathtt{h}}^{4}+4\right)\right)}{64{\mathtt{c}}_{\mathtt{h}}^{21/2}}\cos(2x)\end{bmatrix}\,.$
(B.49)

Finally, by (B.47) and (B.48), we get that the term $\mathrm{I\\!X}$ in (B.34)
is

$\mathrm{I\\!X}=-\zeta_{2}^{+}{\mathrm{B}_{2}^{+}}-\frac{\alpha_{2}^{+}}{4\text{${\mathtt{c}}_{\mathtt{h}}^{2}$}}{\mathrm{A}_{-2}^{+}}-\zeta_{3}^{+}\mathrm{A}_{2}^{+}-\mathscr{L}_{0,0}^{\,\rotatebox[origin={c}]{90.0}{$\perp$}}\left({\mathrm{Z}_{2}^{-}}+f_{\mathcal{W}_{4}^{-}}\right)\,,$
(B.50)

where

$\mathscr{L}_{0,0}^{\,\rotatebox[origin={c}]{90.0}{$\perp$}}{\mathrm{Z}_{2}^{-}}=\begin{bmatrix}\kappa_{4}\cos(2x)\\\\[11.38109pt]
\varsigma_{4}\sin(2x)\end{bmatrix}$ (B.51)

where

$\displaystyle\kappa_{4}:=\frac{1}{64{\mathtt{c}}_{\mathtt{h}}^{29/2}}\big{(}-(a_{1}^{[1]})^{3}({\mathtt{c}}_{\mathtt{h}}^{4}+3){\mathtt{c}}_{\mathtt{h}}^{3}+(a_{1}^{[1]})^{2}p_{1}^{[1]}(3{\mathtt{c}}_{\mathtt{h}}^{8}+31{\mathtt{c}}_{\mathtt{h}}^{4}+18){\mathtt{c}}_{\mathtt{h}}^{2}$
$\displaystyle\qquad-
a_{1}^{[1]}(p_{1}^{[1]})^{2}(2{\mathtt{c}}_{\mathtt{h}}^{12}+49{\mathtt{c}}_{\mathtt{h}}^{8}+101{\mathtt{c}}_{\mathtt{h}}^{4}+36){\mathtt{c}}_{\mathtt{h}}+(p_{1}^{[1]})^{3}(11{\mathtt{c}}_{\mathtt{h}}^{12}+67{\mathtt{c}}_{\mathtt{h}}^{8}+86{\mathtt{c}}_{\mathtt{h}}^{4}+24)\big{)}$
$\displaystyle\varsigma_{4}:=\frac{({\mathtt{c}}_{\mathtt{h}}^{4}+1)}{{128{\mathtt{c}}_{\mathtt{h}}^{31/2}}}\big{(}-(a_{1}^{[1]})^{3}({\mathtt{c}}_{\mathtt{h}}^{4}+3){\mathtt{c}}_{\mathtt{h}}^{3}+2(a_{1}^{[1]})^{2}p_{1}^{[1]}({\mathtt{c}}_{\mathtt{h}}^{8}+14{\mathtt{c}}_{\mathtt{h}}^{4}+9){\mathtt{c}}_{\mathtt{h}}^{2}$
$\displaystyle\qquad-
a_{1}^{[1]}(p_{1}^{[1]})^{2}(31{\mathtt{c}}_{\mathtt{h}}^{8}+89{\mathtt{c}}_{\mathtt{h}}^{4}+36){\mathtt{c}}_{\mathtt{h}}+2(p_{1}^{[1]})^{3}(19{\mathtt{c}}_{\mathtt{h}}^{8}+37{\mathtt{c}}_{\mathtt{h}}^{4}+12)\big{)}$

In conclusion, by (B.34), (B.36), (B.40), (B.44) and (B.50) we deduce that

$P_{0,3}f_{1}^{+}=\mathscr{L}_{0,0}^{\,\rotatebox[origin={c}]{90.0}{$\perp$}}\left(-{\mathrm{W}_{2}^{-}}+{\mathrm{X}_{2}^{-}}+{\mathrm{Y}_{2}^{-}}-{\mathrm{Z}_{2}^{-}}\right)+(\tau_{1}^{+}-\zeta_{2}^{+}){\mathrm{B}_{2}^{+}}+\frac{\ell_{1}^{+}-\alpha_{2}^{+}}{4{\mathtt{c}}_{\mathtt{h}}^{2}}{\mathrm{A}_{-2}^{+}}-\zeta_{3}^{+}\mathrm{A}_{2}^{+}+\widetilde{f}_{\mathcal{W}_{4}^{+}}$

which proves the expansion of $P_{0,3}f_{1}^{+}$ in (3.21) with

$\displaystyle\begin{bmatrix}\mathfrak{a}_{0,3}\cos(2x)\\\
\mathfrak{b}_{0,3}\sin(2x)\end{bmatrix}:=\begin{bmatrix}(\kappa_{1}+\kappa_{2}+\kappa_{3}-\kappa_{4})\cos(2x)\\\
(\varsigma_{1}+\varsigma_{2}+\varsigma_{3}-\varsigma_{4})\sin(2x)\end{bmatrix}$
$\displaystyle\qquad+(\tau_{1}^{+}-\zeta_{2}^{+})\mathrm{B}_{2}^{+}(x)+\frac{\ell_{1}^{+}-\alpha_{2}^{+}}{4{\mathtt{c}}_{\mathtt{h}}^{2}}\mathrm{A}_{-2}^{+}(x)-\zeta_{3}^{+}\mathrm{A}_{2}^{+}(x)\,,$

with $\kappa_{i}$, $\varsigma_{i}$, $i=1,\dots,4$ in (B.37), (B.41), (B.45),
(B.51), $\mathrm{B}_{2}^{+}$, $\mathrm{A}_{-2}^{+}$ $\mathrm{A}_{2}^{+}$ in
(B.9), $\zeta_{2}^{+}$, $\alpha_{2}^{+}$, $\zeta_{3}^{+}$ in (B.11) and
$\ell_{1}^{+}$, $\tau_{1}^{+}$ in (B.15), resulting in the coefficients
$\mathfrak{a}_{0,3}$ and $\mathfrak{b}_{0,3}$ in (3.22).

Computation of $P_{1,2}f_{0}^{-}$. By (3.13), (3.14) and the fact that ${\cal
B}_{1,0}f_{0}^{-}={\cal B}_{0,1}f_{0}^{-}={\cal B}_{0,2}f_{0}^{-}=0$, the term
$P_{1,2}f_{0}^{-}$ reduces to

$\displaystyle P_{1,2}f_{0}^{-}$
$\displaystyle=-\frac{1}{2\pi\mathrm{i}\,}\oint_{\Gamma}\frac{(\mathscr{L}_{0,0}-\lambda)^{-1}}{\lambda}\mathcal{J}{\cal
B}_{1,2}f_{0}^{-}\mathrm{d}\lambda$ $\displaystyle\
+\frac{1}{2\pi\mathrm{i}\,}\oint_{\Gamma}\frac{(\mathscr{L}_{0,0}-\lambda)^{-1}}{\lambda}\mathcal{J}{\cal
B}_{0,1}(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{1,1}f_{0}^{-}\mathrm{d}\lambda=:\mathrm{X}+\mathrm{X\\!I}\,.$ (B.52)

We now compute the two terms.

$\mathrm{X}\big{)}$ By (3.46) we have $\mathcal{J}{\cal
B}_{1,2}f_{0}^{-}=\mathrm{i}\,a_{3}f_{0}^{-}+{\mathrm{i}\,\mathrm{W}_{-2}^{-}}$
with ${\mathrm{W}_{-2}^{-}}(x):=p_{2}^{[2]}{\footnotesize\begin{bmatrix}0\\\
\cos(2x)\end{bmatrix}}\,$ and, by (B.4),

$\mathrm{X}\stackrel{{\scriptstyle\eqref{residuethm}}}{{=}}-\mathrm{i}\,\mathscr{L}_{0,0}^{\,\rotatebox[origin={c}]{90.0}{$\perp$}}{\mathrm{W}_{-2}^{-}}\,,\quad\mathscr{L}_{0,0}^{\,\rotatebox[origin={c}]{90.0}{$\perp$}}{\mathrm{W}_{-2}^{-}}(x)\stackrel{{\scriptstyle\eqref{L0inv}}}{{=}}p_{2}^{[2]}\begin{bmatrix}{\mathtt{c}}_{\mathtt{h}}^{-4}\cos(2x)\\\
\dfrac{1+{\mathtt{c}}_{\mathtt{h}}^{4}}{2{\mathtt{c}}_{\mathtt{h}}^{5}}\sin(2x)\end{bmatrix}\,.$
(B.53)

$\mathrm{X\\!I}\big{)}$ By Lemmata B.7, B.4 one has

$(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{0,1}(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{1,1}f_{0}^{-}=\mathrm{i}\,{\mathtt{c}}_{\mathtt{h}}^{-\frac{1}{2}}{\mathrm{B}_{2}^{+}}+\mathrm{i}\,\frac{1}{2{\mathtt{c}}_{\mathtt{h}}^{3/2}}{\mathrm{A}_{-2}^{+}}+\mathcal{O}(\lambda^{-1}:\lambda)\,,$

and therefore

$\mathrm{X\\!I}=\mathrm{i}\,{\mathtt{c}}_{\mathtt{h}}^{-\frac{1}{2}}{\mathrm{B}_{2}^{+}}+\mathrm{i}\,\frac{1}{2{\mathtt{c}}_{\mathtt{h}}^{3/2}}{\mathrm{A}_{-2}^{+}}\,.$
(B.54)

In conclusion, by (B.52), (B.53) and (B.54)

$P_{1,2}f_{0}^{-}=-\mathrm{i}\,(\mathscr{L}_{0,0}^{\,\rotatebox[origin={c}]{90.0}{$\perp$}}){\mathrm{W}_{-2}^{-}}+\mathrm{i}\,{\mathtt{c}}_{\mathtt{h}}^{-\frac{1}{2}}{\mathrm{B}_{2}^{+}}+\mathrm{i}\,\frac{1}{2{\mathtt{c}}_{\mathtt{h}}^{3/2}}{\mathrm{A}_{-2}^{+}}\,,$
(B.55)

which, in view of (B.53), proves the expansion of $P_{1,2}f_{0}^{-}$ in (3.21)
with

$\displaystyle\begin{bmatrix}\mathfrak{a}_{1,2}\cos(2x)\\\
\mathfrak{b}_{1,2}\sin(2x)\end{bmatrix}:=\begin{bmatrix}-p_{2}^{[2]}{\mathtt{c}}_{\mathtt{h}}^{-4}\cos(2x)\\\
-p_{2}^{[2]}\frac{1+{\mathtt{c}}_{\mathtt{h}}^{4}}{2{\mathtt{c}}_{\mathtt{h}}^{5}}\sin(2x)\end{bmatrix}+{\mathtt{c}}_{\mathtt{h}}^{-\frac{1}{2}}\mathrm{B}_{2}^{+}(x)+\frac{1}{2{\mathtt{c}}_{\mathtt{h}}^{3/2}}\mathrm{A}_{-2}^{+}(x)\,,$

with $\mathrm{B}_{2}^{+}$ and $\mathrm{A}_{-2}^{+}$ in (B.9), resulting in the
coefficients $\mathfrak{a}_{1,2}$ and $\mathfrak{b}_{1,2}$ given in (3.22).

Computation of $P_{2,1}f_{0}^{-}$. By (3.13) and the fact that ${\cal
B}_{1,0}f_{0}^{-}={\cal B}_{0,1}f_{0}^{-}$ and ${\cal B}_{2,1}=0$ the term
$P_{2,1}f_{0}^{-}$ reduces to

$\displaystyle P_{2,1}f_{0}^{-}$
$\displaystyle=\frac{1}{2\pi\mathrm{i}\,}\oint_{\Gamma}\frac{(\mathscr{L}_{0,0}-\lambda)^{-1}}{\lambda}\mathcal{J}{\cal
B}_{0,1}(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{2,0}f_{0}^{-}\mathrm{d}\lambda$
$\displaystyle+\frac{1}{2\pi\mathrm{i}\,}\oint_{\Gamma}\frac{(\mathscr{L}_{0,0}-\lambda)^{-1}}{\lambda}\mathcal{J}{\cal
B}_{1,0}(\mathscr{L}_{0,0}-\lambda)^{-1}\mathcal{J}{\cal
B}_{1,1}f_{0}^{-}\mathrm{d}\lambda\,.$

By repeated use of (B.3) and Lemma B.2 one finds that
$P_{2,1}f_{0}^{-}\in\mathcal{W}_{1}^{-}$ as stated in (3.21).

## References

  * [1] C. J. Amick, L. E. Fraenkel, and J. F. Toland. On the Stokes conjecture for the wave of extreme form. Acta Math. 148, 193-214, 1982.
  * [2] P. Baldi, M. Berti, E. Haus and R. Montalto. _Time quasi-periodic gravity water waves in finite depth_. Inventiones Math. 214 (2): 739–911, 2018.
  * [3] T. Benjamin. _Instability of periodic wavetrains in nonlinear dispersive systems_ , Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. Volume 299 Issue 1456, 1967.
  * [4] T. Benjamin and J. Feir. _The disintegration of wave trains on deep water_ , Part 1. Theory. J. Fluid Mech. 27(3): 417-430, 1967.
  * [5] M. Berti, L. Franzoi and A. Maspero. _Traveling quasi-periodic water waves with constant vorticity_ , Archive for Rational Mechanics, 240: 99–202, 2021.
  * [6] M. Berti, L. Franzoi and A. Maspero. _Pure gravity traveling quasi-periodic water waves with constant vorticity_ , arXiv:2101.12006, 2021, to appear in Comm. Pure Appl. Math.
  * [7] M. Berti, A. Maspero and P. Ventura. _Full description of Benjamin-Feir instability of Stokes waves in deep water_ , Inventiones Math., 230, 651–711, 2022.
  * [8] M. Berti, A. Maspero and P. Ventura. _Benjamin-Feir instability of Stokes waves_ , Rend. Lincei Mat. Appl., 33, 399-412, 2022.
  * [9] M. Berti, A. Maspero and P. Ventura. _On the analyticity of the Dirichlet-Neumann operator and Stokes waves_ , Rend. Lincei Mat. Appl., 33, 611–650, 2022.
  * [10] M. Berti, A. Maspero and P. Ventura. Benjamin-Feir instability of Stokes waves in finite depth, https://arxiv.org/abs/2204.00809, 2022.
  * [11] T. Bridges and A. Mielke. _A proof of the Benjamin-Feir instability._ Arch. Rational Mech. Anal. 133(2): 145–198, 1995.
  * [12] T. Bridges and D. Ratliff. _On the elliptic-hyperbolic transition in Whitham modulation theory._ SIAM J. Appl. Math. 77 (6): 1989–2011, 2017.
  * [13] J. Bronski, V. Hur and M. Johnson. _Modulational Instability in Equations of KdV Type._ In: Tobisch E. (eds) New Approaches to Nonlinear Waves. Lecture Notes in Physics, vol 908. Springer, 2016.
  * [14] J. Bronski and M. Johnson. _The modulational instability for a generalized Korteweg-de Vries equation._ Arch. Ration. Mech. Anal. 197(2): 357–400, 2010.
  * [15] G. Chen and Q. Su. _Nonlinear modulational instabililty of the Stokes waves in 2d full water waves._ Commun. Math. Phys, doi.org/10.1007/s00220-023-04747-0, 2023.
  * [16] W. Craig and D. P. Nicholls, _Traveling gravity water waves in two and three dimensions._ European Journal of Mechanics - B/Fluids, Volume 21, Issue 6, Pages 615-641, 2002\.
  * [17] W. Craig and C. Sulem. _Numerical simulation of gravity waves._ J. Comput. Phys., 108(1): 73–83, 1993.
  * [18] R. Creedon, B. Deconinck. A High-Order Asymptotic Analysis of the Benjamin-Feir Instability Spectrum in Arbitrary Depth, J. Fluid Mech. 956, A29, 2023.
  * [19] R. Creedon, B. Deconinck., O. Trichtchenko. High-Frequency Instabilities of Stokes Waves, J. Fluid Mech. 937, A24, 2022.
  * [20] B. Deconinck and K. Oliveras. _The instability of periodic surface gravity waves._ J. Fluid Mech., 675: 141–167, 2011.
  * [21] B. Deconinck, S. A. Dyachenko, P. M. Lushnikov, A. Semenova. _The instability of near-extreme Stokes waves_ , arXiv:2211.05473, 2022.
  * [22] R. Feola and F. Giuliani. Quasi-periodic traveling waves on an infinitely deep fluid under gravity. arXiv:2005.08280, to appear on Memoires American Mathematical Society.
  * [23] J. D. Fenton. _A Fifth-Order Stokes Theory for Steady Waves._ Journal of Waterway, Port, Coastal, and Ocean Engineering, Volume 111, Issue 2: 216–234, 1985.
  * [24] T. Gallay and M. Haragus. _Stability of small periodic waves for the nonlinear Schrödinger equation._ J. Differential Equations, 234: 544–581, 2007.
  * [25] M. Haragus and T. Kapitula. _On the spectra of periodic waves for infinite-dimensional Hamiltonian systems_. Phys. D, 237: 2649–2671, 2008.
  * [26] V. Hur and M. Johnson. _Modulational instability in the Whitham equation for water waves._ Stud. Appl. Math. 134(1): 120–143, 2015.
  * [27] V. Hur and A. Pandey. _Modulational instability in nonlinear nonlocal equations of regularized long wave type._ Phys. D, 325: 98–112, 2016.
  * [28] V. Hur and Z. Yang. _Unstable Stokes waves_. arXiv:2010.10766.
  * [29] J. Jin, S. Liao and Z. Lin. _Nonlinear modulational instability of dispersive PDE models._ Arch. Ration. Mech. Anal. 231(3): 1487-–1530, 2019.
  * [30] M. Johnson. _Stability of small periodic waves in fractional KdV type equations._ SIAM J. Math. Anal. 45: 2529–3228, 2013.
  * [31] R. S. Johnson, On the Modulation of Water Waves in the Neighbourhood of $kh\approx 1.363$, Proceedings of the Royal Society of London. Series A, Math. and Physical Sciences, Vol. 357, No. 1689, pp. 131-141, 1977.
  * [32] T. Kakutani, K. Michihiro. Marginal State of Modulational Instability –Note on Benjamin-Feir Instability–. Journal of the Physical Society of Japan, 52(12), 4129–4137, 1983.
  * [33] G. Keady, J. Norbury, On the existence theory for irrotational water waves, Math. Proc. Cambridge Philos. Soc. 83, no. 1, 137-157, 1978.
  * [34] A. O. Korotkevich, A. I. Dyachenko and V. E. Zakharov, Numerical simulation of surface waves instability on a homogeneous grid, Physica D: Nonlinear Phenomena, Volumes 321-322, 51-66, 2016.
  * [35] K. Leisman, J. Bronski, M. Johnson, and R. Marangell. _Stability of Traveling Wave Solutions of Nonlinear Dispersive Equations of NLS Type._ Arch. Rational Mech. Anal., 240: 927-–969, 2021.
  * [36] T. Levi-Civita.  Détermination rigoureuse des ondes permanentes d’ ampleur finie, Math. Ann. 93, 264-314, 1925.
  * [37] M. J. Lighthill, Contribution to the theory of waves in nonlinear dispersive systems, IMA Journal of Applied Mathematics, 1, 3, 269-306, 1965.
  * [38] J. B. McLeod, The Stokes and Krasovskii conjectures for the wave of greatest height, Stud. Appl. Math. 98, no. 4, 311-333, 1997.
  * [39] A. Nekrasov. On steady waves. Izv. Ivanovo-Voznesenk. Politekhn. 3, 1921\.
  * [40] H. Nguyen and W. Strauss. _Proof of modulational instability of Stokes waves in deep water_. Comm. Pure Appl. Math., Volume 76, Issue 5, 899–1136, 2023.
  * [41] P. I. Plotnikov, Proof of the Stokes conjecture in the theory of surface waves. Dinamika Sploshn. Sredy, 57, 41-76, 1982.
  * [42] H. Segur, D. Henderson, J. Carter and J. Hammack. _Stabilizing the Benjamin-Feir instability_. J. Fluid Mech. 539: 229–271, 2005.
  * [43] G. Stokes. On the theory of oscillatory waves. Trans. Cambridge Phil. Soc. 8: 441–455, 1847.
  * [44] Y. Sedletsky. A Fifth-Order Nonlinear Schrödinger Equation for Waves on the Surface of Finite-Depth Fluid. Ukrainian Journal of Physics, 66(1), 41, 2021.
  * [45] A. Slunyaev. A high-order nonlinear envelope equation for gravity waves in finite-depth water. J. Exp. Theor. Phys. 101, 926–941, 2005
  * [46] D. Struik. Détermination rigoureuse des ondes irrotationelles périodiques dans un canal á profondeur finie. Math. Ann. 95: 595–634, 1926.
  * [47] J. F. Toland, On the existence of a wave of greatest height and Stokes conjecture, Proc. Roy. Soc. London Ser. A 363, 1715, 469-485, 1978.
  * [48] G.B. Whitham. Non-linear dispersion of water waves. J. Fluid Mech, volume 26 part 2 pp. 399-412, 1967\.
  * [49] G.B. Whitham. Linear and Nonlinear Waves. J. Wiley-Sons, New York, 1974.
  * [50] J. Wilkening, V. Vasan. Comparison of five methods of computing the Dirichlet-Neumann operator for the water wave problem. Contemporary Mathematics 635, 2015.
  * [51] V. Zakharov. The instability of waves in nonlinear dispersive media. J. Exp.Teor.Phys. 24 (4), 740-744, 1967.
  * [52] V. Zakharov. Stability of periodic waves of finite amplitude on the surface of a deep fluid. Zhurnal Prikladnoi Mekhaniki i Teckhnicheskoi Fiziki 9(2): 86–94, 1969.
  * [53] V. Zakharov and V. Kharitonov. Instability of monochromatic waves on the surface of a liquid of arbitrary depth. J Appl Mech Tech Phys 11, 747-751, 1970.
  * [54] V. Zakharov and L. Ostrovsky. _Modulation instability: the beginning._ Phys. D, 238(5): 540–548, 2009.