Complete the following paragraph to describe the two processes by which the influenza virus evolves.

Our within-host model is a simple target cell-limited model, but the decrease in infectible cells as the infection proceeds can be qualitatively interpreted as any and all antigenicity-agnostic limiting factors that come into play as the virion and infected cell populations grow. This could include the action of innate immunity, which acts in part by killing infected cells and by rendering healthy cells difficult to infect via inflammation (see immunological review in Appendix Section A2). The key mechanistic role played by the target cells in our model is to introduce a non-antigenic limiting factor on infections. This prevents an infection from repeatedly evolving out from under successive well-matched antibody responses, as occurs in HIV and in influenza patients with compromised immune systems (Xue et al., 2017). These factors limit the virus regardless of whether the infection remains confined to the upper respiratory tract or also infect the lower respiratory tract, thereby gaining access to additional target cells (Koel et al., 2019). The key is that non-antigenic factors prevent persistent large virus populations.

Similarly, the antibody response we introduce at 48 hours is qualitative—it is modeled simply as an increase in the virion decay rate for antibody-matching virions. This response could represent IgA targeted at HA, but other antigenicity-specific modes of virus control could also be subsumed under the increased virion decay rate, for instance IgG antibodies, antibody dependent cellular cytotoxicity (ADCC), antibodies against the neuraminidase protein, or others. The key point we establish is that none of these mechanisms efficiently replication select because they all emerge once non-antigenic limiting factors have come into play. They speed clearance, but should not substantially alter virus evolution. A corollary to this point is that there are many mechanisms—and interventions—that could reduce the severity of influenza infections without substantially speeding up antigenic evolution, including universal vaccines.

A key proposal of our study is that population-level antigenic selection and homotypic protection are mediated by antibody neutralization (likely sIgA) at the point of transmission. Currently, empirical evidence for antibody protection at the point of transmission is mostly indirect. Most of this evidence comes from human and animal challenge studies (Clements et al., 1986; Memoli et al., 2020; Le Sage et al., 2020). In these studies, individuals who are challenged with the same antigenic variant sometimes display apparent sterilizing immunity, but other times develop detectable infections. The study by Memoli et al., 2020 is notable for having used very large inocula—106 or 107 TCID50. Despite these high doses, two of the challenge subjects had neither detectable virus nor seroconversion. Similarly, ferret experiments (Le Sage et al., 2020) found that many experienced ferrets developed sufficiently sterilizing immunity to prevent the virus from ever being detected, while some experienced ferrets showed briefly detectable infections that were then cleared. While we cannot rule out a powerful immediate cellular response that was differentially evaded in the various subjects, we believe that our model, coupled with existing understanding of the timing of cellular responses and the speed of influenza virus replication, provides a more parsimonious explanation.

Another limitation of our study is that, while we put forward mucosal sIgA as a biologically-documented potential mechanism of immune protection at the point of inoculation that would not lead to strong selection during early viral replication, no modeling study can establish such a mechanism without empirical investigation. Our study reveals that such an empirical investigation would be of substantial scientific value.

We model neutralization at the point of transmission as a binomial process. Each virion is independently neutralized with a probability κi that depends on its antigenic phenotype and the host’s immune history. As we discuss in Appendix Section A5.2, this assumption of independence may be violated in practice. Careful experiments are required to develop a more realistic model of neutralization at the point of transmission.

Moreover, individual variation in immune system properties and complex effects of host immune history (Fonville et al., 2014; Lee et al., 2019) mean that even a pair of hosts who have both been previously exposed to the currently circulating variant may exert different selection pressures at the point of transmission. Modeling neutralization as a series of independent events that depend only on host history and virus phenotype is a baseline: it allows us to establish that antigenic selection at the point of transmission is possible and to show what its consequences might be. But a more realistic model will be required to predict the selective pressures imposed by real hosts with real immune histories on real virions.

Finally, while we believe neutralization at the point of transmission is a crucial mechanism of protection against detectable reinfection for influenza, this may not be true for all RNA viruses. Some, such as measles and varicella, have long incubation periods even in naive hosts and induce reliable, long-lasting immune protection against detectable reinfection. This may be because they replicate slowly enough that they cannot ‘outrun’ the adaptive response as influenza can. Neither shows influenza virus-like patterns of clocklike immune escape, suggesting that (1) escape mutants may be less available and (2) the adaptive response acts on a small population and is forceful.

We parameterized our models based on estimates from previous studies, but there are not good estimates for several important quantities. There are no high quality estimates of the rate of antibody-mediated neutralization in the presence of a homotypic antibody response or of how much this rate is reduced by particular antigenic substitutions, and there may be substantial inter-individual variation (Lee et al., 2019). Within-host timeseries data from antigenically heterogeneous infections are needed to estimate these quantities. Similarly, there are no good empirical data on the size of the bottlenecks that precede or follow the sIgA bottleneck (Figure 3A). Better estimates of these bottlenecks, and of the probabilities of neutralization for individual virions encountering mucosal sIgA antibodies, would give more certainty about the strength of inoculation selection relative to neutral founder effects.

We also do not have a clear sense of exactly how neutralization probability in the respiratory tract mucosa (parameter κ in our model) and neutralization rate during replication (parameter k in our model) relate. We expect them to be positively related, but the exact strength and shape of this relationship is unknown. Knowing whether major antigenic changes reduce both equally or reduce one more than the other could help us better quantify the potential strength of inoculation selection and replication selection. In short, better mechanistic understanding of mucosal antibody neutralization could be extremely valuable for understanding and potentially predicting influenza virus evolution.

How readily a particular individual host or host population helps new antigenic variants reach within-host consensus depends upon several unknown quantities: (1) how host susceptibility changes with extent of antigenic dissimilarity, (2) the ratio of virions that encounter sIgA to virions that found the infection (v/b), (3) the probability that a single new antigenic variant virion inoculated alongside old antigenic variant virions evades neutralization 1-κm (Figure 3E–H), (4) the duration τ from the onset of the antibody response to the time of transmission. Better empirical estimates of these quantities could shed light on how the distribution of host immunity shapes antigenic evolution. However, over a range of biologically plausible parameter values, our model contradicts the existing hypothesis that antigenic novelty appears when moderately immune hosts fail to block transmission and then select upon a growing virus population (Grenfell et al., 2004; Volkov et al., 2010). For influenza viruses, hosts whose mucosal immunity regularly blocks old antigenic variant transmission may be crucial. Mucosal immunity not only produces a population-level advantage for new variants but may also play a role in their within-host emergence (Figure 3E,H).

Experienced hosts are undoubtedly heterogeneous in their immunity to a given influenza variant (Lee et al., 2019), so the overall population average protection against homotypic reinfection with variant i, zi, is in fact an average over experienced hosts. Our model implies that the degree of neutralization difference between ancestral variant virions and new antigenic variant virions at the point of transmission strongly affects the probability of inoculation selection. Hosts with more focused immune responses—highly-specific antibodies that neutralize old antigenic variant virions well and new antigenic variant virions poorly—could be especially good inoculation selectors and important sources of population-level antigenic selection. Hosts who develop less specific memory responses, such as very young children (Neuzil et al., 2006), could be less important. Similarly, immune-compromised hosts are excellent replication selectors (Xue et al., 2017; Lumby et al., 2020), and so their role in the generation of antigenic novelty and their impact on overall population-level diversification rates deserve further study.

Prior modeling has suggested that despite repeated tight bottlenecks at the point of transmission, evolution of influenza viruses should resemble evolution in idealized large populations (Sigal et al., 2018). In large populations, advantageous variants with small selective advantages should gradually fix and weakly deleterious variants should be purged. In Sigal et al., 2018, diversity is rapidly generated and fit variants are selected to frequencies at which they are likely to pass through even a tight bottleneck. This is likely true of the phenotypes modeled, which include receptor-binding avidity and virus budding time (Sigal et al., 2018). These phenotypes affect virus fitness throughout the timecourse of infection, so they can be efficiently replication-selected (where selection is manifest in the direct competition to infect cells rather than the indirect competition to escape antibodies). Indeed, next-generation sequencing studies have found observable adaptative evolution of non-antigenic phenotypes in individual humans infected with avian H5N1 viruses (Welkers et al., 2019).

Seasonal influenza antigenic evolution does not resemble idealized large population evolution. Within an antigenic cluster, influenza viruses acquire substitutions that change the antigenic phenotype by small amounts. Given the large influenza virus populations within individual hosts, we might expect a quasi-continuous directional pattern of evolution away from prior population immunity. Between clusters, evolution is indeed strongly directional: only ‘forward’ cluster transitions are observed. These are jumps—large antigenic changes. But incremental within-cluster evolution is not directional: the virus often evolves ‘backward’ or ‘sideways’ in antigenic space toward previously circulated variants (see Figures 1 and 2 of Smith et al., 2004).

This noisy jump pattern is easy to explain in light of the weakness of replication selection and the importance of antigenic founder effects. Selection acts on the small sub-sample of donor-host diversity that passes through the excretion, inter-host, and mucus bottlenecks to encounter the sIgA bottleneck. Evolution via inoculation selection is therefore slower and more affected by stochasticity than evolution via replication selection (Figure 3E–F). It resembles evolution in small populations—weakly adaptive and weakly deleterious substitutions become nearly neutral (Kimura, 1968; Ohta, 1992). If influenza virus evolution were not nearly neutral for small-effect substitutions at the within-host scale, it would be surprising to observe ‘backward’ antigenic changes and noisy evolution at higher scales (Figure 5A–D). In fact, there may be analogous ‘neutralizing’ population dynamics at higher scales as well, and those may also be needed to explain population-level noisiness. But whatever happens at higher scales, within-host replication selection creates a strong directional bias in population-level antigenic diversification, introducing many small forward antigenic changes (Figure 5A–D). Inoculation selection does not necessarily do this.

Local influenza virus lineages rarely persist between epidemics (Russell et al., 2008; Bedford et al., 2015), and so new antigenic variants must establish chains of infections in other geographic locations in order to survive. New antigenic variant chains are most often founded when inoculations are common—that is, when existing variants are causing epidemics. Epidemics result in high levels of local competition between extant and new antigenic variant viruses for susceptible hosts (Hartfield and Alizon, 2015) as well as metapopulation-scale competition to found epidemics in other locations. These dynamics could create tight bottlenecks between epidemics similar to those that occur between hosts, resulting in dramatic epidemic-to-epidemic diversity losses. That said, if immune hosts are present at the start of an epidemic, there will not be asynchrony between diversity and selection pressure, so new variants may pass through between-epidemic bottlenecks more readily than through between-host bottlenecks. Further work is needed to elucidate mechanisms at the population and meta-population scales.

Our work suggests a simple mechanism by which accumulating immunity to an antigenic variant could produce punctuated population-level antigenic evolution. Population-level modeling has shown that influenza virus global epidemiological and phylogenetic patterns can be reproduced if new antigenic variants emerge at the population level with increasing frequency the longer an old antigenic variant circulates (Koelle et al., 2009).

If immunity to the old variant only gives new variants a population-level transmission advantage (population-level selection), we anticipate a constant rate of population-level antigenic diversification, with selective sweeps once a new variant has a sufficient population-level advantage over the old variant. If population immunity also helps new variants become the dominant within-host variant through inoculation selection, increasing population immunity to an old variant can produce increasing rates of new variant emergence (Figure 6A). Whether this occurs in practice depends on the ecology of hosts and the relative strength of inoculation selection versus drift in partially and fully immune hosts (Figure 4).

Potential synergy between brief antigenic replication selection late in infection and subsequent inoculation selection (Figure 4) could further promote new variant emergence as population immunity accumulates. However, it is difficult to estimate how much this synergy matters in practice without knowing more about the kinetics of homotypic and heterotypic re-inoculation and reinfection.

One suggestive population-level pattern is that new antigenic variants frequently are observed quasi-simultaneously on multiple branches of the influenza virus phylogeny shortly prior to sweeping (see Appendix Section A6). This suggests that emergence rate and population-level selection pressure do increase together. That said, alternative explanations are possible, such as reduced rates of stochastic loss of new variants at higher scales (e.g. the epidemic scale) with increased population immunity.

Regardless of the strength of these effects, however, our prediction is that new antigenic variant infection foundation events should constitute a rare but non-negligible fraction of transmission events: on the order of 1 in 105 or 1 in 104. This parsimoniously explains why new antigenic variants lineages are hard to observe prior to undergoing positive population-level selection but are readily available to be selected upon (see Appendix Section A6) once that population-level selection pressure has become sufficiently strong. Before that point, they could be rendered unobservable by population-level neutralizing dynamics or by clonal interference from non-antigenic weakly adaptive mutants (Strelkowa and Lässig, 2012).

The slow rate of antigenic evolution of A/H3N2 viruses in swine lends further support to this argument. A/H3N2 viruses accumulate genetic mutations at a similar pace in swine and humans, but antigenic evolution is much slower in swine (ESNIP3 consortium et al., 2016; de Jong et al., 2007). Slaughter for meat means that pig population turnover is high. It follows that the frequency of experienced hosts rarely becomes sufficient to facilitate the appearance of observable new antigenic variants.

The population-level emergence of new antigenic variants, in other words, tracks the accumulation of immunity in the population, not the accumulation of genetic diversity. This suggests that A/H3N2 evolution is indeed selection-limited, not diversity-limited. But much generated antigenic diversity is invisible to surveillance: in the absence of positive selection, it is likely to be lost at bottlenecks.

The inoculation and replication selection paradigm has implications for the understanding and management of other pathogens. For example, HIV does not readily evolve resistance to contemporary pre-exposure prophylaxis (PrEP) antiviral drugs, but it can do so when these antivirals are taken by an individual who is already infected with HIV (Partners PrEP Study Team et al., 2016). Developing resistance at the moment of exposure is a difficult problem of inoculation selection for the virus, but developing resistance during an ongoing infection is an easier problem of replication selection. Selection on small un-diverse introduced populations may also be of interest in invasion biology and island biogeography.

The within-host frequency of the new variant, fm, obeys a replicator equation of the form:

(12) d⁢fmd⁢t=fm⁢(1-fm)⁢δ⁢(t)

where δ⁢(t) is the fitness advantage of the new antigenic variant over the old antigenic variant at time t (see Appendix Section A3 for a derivation).

If the variant is neutral in the absence of antibodies, then δ⁢(t)=k⁢(cw-cm) if t>tM and δ⁢(t)=0 otherwise. Let te denote the first time during the infection that a de novo mutation produces a surviving new antigenic variant lineage. If te≤tM, then at a time t≥tM:

(13) fm⁢(t)=eδ⁢(t-tM)eδ⁢(t-tM)+fM-1-1

where δ=k⁢(cw-cm) and fM=fm⁢(tM).

When additional mutations after the first cannot be neglected, we add a correction term to d⁢fmd⁢t for te<t<min⁡{tM,tpeak} (Figure 7), where tpeak is the time of peak virus population:

(14) d⁢fmd⁢t=fm⁢(1-fm)⁢δ⁢(t)+μ⁢ℛ0⁢dv⁢(1-2⁢fm+fm2)

which for δ≠0 yields:

(15) fm(t)≈Aexp⁡(δt)+μg0Aexp⁡(δt)+μg0−δA=f01−f0(μg0−δ)−μg0f0

And when δ=0:

(16) fm(t)≈11−f0+μg0t−111−f0+μg0t

where g0=ℛ0dv and f0=fm(te). See Appendix Section A3.2 for derivations and discussion.

In our stochastic model, new variant lineages that survive stochastic extinction are produced by de novo mutation according to a continuous-time, variable-rate Poisson process. The cumulative distribution function for the time of the first successful mutation, te, depends on the mutation rate µ and the per-capita rate at which old antigenic variant virions are produced, gw⁢(t)=rw⁢β⁢C⁢(t). It also depends on psse, the probability that the generated new antigenic variant survives stochastic extinction. Denoting the new variant per-capita virion production rate gm⁢(t)=rm⁢β⁢C⁢(t), we calculate psse using a branching process approximation (Ball et al., 2016).

(17) psse=gm⁢(t)gm⁢(t)+dv+k⁢M⁢(t)⁢max⁡{Em,c⁢Ew}

Surviving mutants therefore occur at a rate λm⁢(t):

(18) λm⁢(t)=μ⁢gw⁢(t)⁢Vw⁢(t)⁢psse

We define the cumulative rate Λm⁢(x):

(19) Λm⁢(x)=∫0xλm⁢(x)⁢𝑑x

It follows that the CDF of the first mutation time te is:

This expression is exact for any given realization of the stochastic model if the realized values of the random variables Vw⁢(t), C⁢(t), and gw⁢(t) are used. In practice, we mainly use it to get a closed form for the CDF of te by making the approximation that C⁢(t)≈Cmax early in infection. This yields approximations for gw⁢(t), gm⁢(t), and Vw⁢(t):

(21) gw⁢(t),gm⁢(t)≈g0≡ℛ0⁢dv

(22) Vw(t){bexp⁡((g0−dv)t)t ≤ tMbexp⁡((g0−dv)tM)exp⁡((g0−dv−cwk)(t−tM))t > tM

The resultant approximate solution for the CDF of new variant mutation times agrees well with simulations (Figure 8).

Black line shows analytically calculated CDF. Blue cumulative histogram shows distribution of new variant mutation times for 250,000 simulations from the stochastic within-host model with (A) an immediate recall response (tM=0) and (B) a realistic recall response at 48 hr post-infection (tM=2). Other model parameters as in Table 1. Note that time of first successful mutation tends to be later with an immediate recall response than with a delayed recall response. This occurs because the cumulative number of viral replication events grows more slowly in time at the start of the infection because of the strong, immediate recall response.

The slightly earlier simulated mutation times in the immediate recall response case (Figure 8A) would only make replication selection more likely in that case than our analytical approximation suggests.

By inverting the within-host replicator equation, we can also calculate the time t*⁢(x,t) by which a new variant must emerge if it is to reach at least frequency x by time t. We show (see Appendix Section A9.4 for derivation) that there are two candidate values for t*, depending on whether the time that the new variant first emerges (te) is before or after the onset of the antibody response (tM):

If te≤tM:

(23) t−∗(x,t)=ln⁡[1−xxexp⁡(δ(t−tM))+1]−ln⁡bg0−dv

If te>tM:

(24) t+*⁢(x,t)≈ln⁡(1-xx)-ln⁡(b)+δ⁢t-cw⁢k⁢tMg0-dv-cw⁢k+δ

It may be that t+*⁢(x,t)>t and t-*⁢(x,t)>t. This indicates that the mutant will be at frequency x if it emerges at t itself. In that case, we therefore have t*⁢(x,t)=t. So combining:

(25) t∗(x,t)={t−∗(x,t)t−∗(x,t) < t and t−∗(x,t) ≤ tMt+∗(x,t)t+∗(x,t) < t and t−∗(x,t) > tMtotherwise

Finally, it is worth noting that in the case of a complete escape mutant (cm=0, δ=cw⁢k), the approximate expression for t+* is exactly equal to the equivalent approximate expression for t-*:

(26) t*⁢(x,t)≈ln⁡(1-xx)-ln⁡(b)+cw⁢k⁢(t-tM)g0-dv

This is a linear function of t.

Given this and the new variant first mutation time CDF calculated in Equation 20, it is straightforward to calculate the probability of replication selection to a given frequency a by time t, assuming that ℛ0w>1 early in infection:

(27) prepl(a,t)=P(te<t*(a,t))=1-e-Λ⁢(t*⁢(a,t))

This analytical model agrees well with simulations (Figure 9). We use it used to calculate the heatmaps shown in Figure 1, with the C⁢(t)≈Cmax early infection approximations that give us a closed form for Λm⁢(x).

Probability of replication selection to one percent (left column) or consensus (right column) by t=3 days post-infection as a function of k and tM for 250,000 simulations from the stochastic within-host model. cw=1,cm=0 (so the fitness difference δ=k). Other model parameters as parameters in Table 1. Dashed lines show analytical prediction and dots show simulation outcomes.

When t*>tM, the integral ∫0t*λm⁢(x)⁢𝑑x can be evaluated piecewise, first from 0 to tM, and then from tM to t*. A similar approach can be used to evaluate the probability density function (PDF) of first mutation times, as needed.

Finally, note that for te<tM, the expression for prepl in practice depends only on δ=(cw-cm)⁢k, not on cw, cm and k separately. In the absence of an antibody response, the mutant is generated with near certainty by 1 day post-infection (P(te<1)≈1, Figure 8B). So when tM≥1, values for prepl calculated with cw=1,cm=0, and δ=k—as in (Figure 1C,D)—will in fact hold for any cw, cm, and k that produce that fitness difference δ.

When ℛ0w<1 early in infection, the probability of replication selection depends on the probability of generating an escape mutant before the infection is extinguished:

(28) prepl≈ℛw⁢(0)1-ℛw⁢(0)⁢b⁢μ⁢psse

See Appendix Section A3.6 for a derivation.

Given contact between an infected host and an uninfected host, we assume that a transmission attempt occurs with probability proportional to the current virus population size Vtot in the infected host:

(29) P(transmit)=1−exp⁡(log⁡(2)VtotV50)

The parameter V50 sets the scaling, and reflects virus population size at which there is a 50% chance of successful transmission to a naive host. We used a default value of V50=1×108 virions. In addition to this probabilistic model, we also consider an alternative threshold model in which a transmission attempt occurs with certainty if Vtot is greater than a threshold θ and does not occur otherwise.

If there are xw old antigenic variant virions and xm new variant virions competing to pass through the cell infection bottleneck, at least one new variant virion passes through with probability:

(30) pcib(xw,xm,b)=1−(xwb)(xw+xmb)

Summing pcib⁢(xw,xm) over the possible values of xw and xm weighted by their joint probabilities yields the new variant’s overall probability of successful onward transmission psurv⁢(fmt,v,b,κ,c):

(31) psurv=∑xm,xwP⁢(xm,xw)⁢pcib⁢(xw,xm,b)

P⁢(xw,xm) is the product of the probability mass functions for xw and xm:

(32) P⁢(xm,xw)=w¯xwxw!⁢m¯xmxm!⁢e-(xw+xm)

where w¯=v⁢(1-fmt)⁢(1-κw) and m¯=v⁢fmt⁢(1-κm).

At low donor-host variant frequencies fmt≪1psurv, psurv can be approximated using the fact that almost all probability is given to xm=0 or xm=1 (0 or 1 new antigenic variant after the sIgA bottleneck).

At least one new antigenic variant survives the sIgA bottleneck with probability:

(33) pinoc⁢(κm,v,fmt)=1-e-v⁢fmt⁢(1-κm)

If fmt is small, there will almost always be at most one such virion (xm=1). That new antigenic variant virion’s probability of surviving the cell infection bottleneck depends upon how many old antigenic variant virions are present after neutralization, xw:

(34) pcib(xw,1,b)=1−(xwb)(xw+1b)=1−xw−b+1xw+1=bxw+1

Summing over the possible xw, we obtain a closed form for the unconditional probability pcib⁢(κw,v,b) (see Appendix Section A9.5 for a derivation):

(35) pcib(κw,v,b)=(1−e−w¯)bw¯+e−w¯∑j=0b−1w¯jj!(1−bj+1)

where w¯=v⁢(1-fmt)⁢(1-κw) is the mean number of old antigenic variant virions present after sIgA neutralization. If b=1, this reduces the just the first term, and it is approximately equal to just the first term when w¯ is large, since e-w¯ becomes small.

It follows that there is an approximate closed form for the probability that a new variant survives the final bottleneck:

(36) psurv⁢(κm,κw,v,b,fmt)≈pinoc⁢(κm,v,fmt)*pcib⁢(κw,v,b)

We use this expression to calculate the analytical new variant survival probabilities shown in the main text, and to gain conceptual insight into the strength of inoculation selection relative to neutral processes (see Appendix Section A4.5).

A given per-virion sIgA neutralization probability κi implies a probability zi that a transmission event involving only virions of variant i fails to result in an infected cell.

Some transmissions fail even without sIgA neutralization; this occurs with probability exp⁡(-v). Otherwise, with probability 1-exp⁡(-v), ni virions must be neutralized by sIgA to prevent an infection. We define zi as the probability of no infection given inoculated virions in need of neutralization (i.e. given ni>0).

The probability that there are no remaining virions of variant i after mucosal neutralization is exp⁡(-v⁢(1-κi)). So we have:

(37) zi=exp⁡(-v⁢(1-κi))-exp⁡(-v)1-exp⁡(-v)

This can be solved algebraically for κi in terms of zi, but it is more illuminating to express κi in terms of the overall probability of no infection given inoculation pno and then find pno in terms of zi:

(38) pno=exp⁡(−v(1−κi))κi=1+ln⁡(pno)v

An infection occurs with probability (1-exp⁡(-v))⁢(1-zi) (the probability of at least one virion needing to be neutralized times the conditional probability that it is not) so pno in terms of zi is:

(39) pno=1-(1-exp⁡(-v))⁢(1-zi)

This yields the same expression for κi in terms of zi as a direct algebraic solution of Equation 37.

For moderate to large v, 1-exp⁡(-v) approaches 1, so pno approaches zi and κi approaches 1+ln⁡(zi)v. This reflects the fact that for even moderately large v, it is almost always the case that ni>0: a least one virion must be neutralized to prevent infection. In those cases, zi can be interpreted as the (approximate) probability of no infection given a transmission event (i.e. as pno).

Note also that seeded infections can also go stochastically extinct; this occurs with approximate probability 1ℛ0. At the start of infection, if there is no antibody response and ℛ0 is large (5 to 10), stochastic extinction probabilities should be low (15 to 110), and equal in immune and naive hosts. We have therefore parametrized our model in terms of zi, the probability that no cell is ever infected, as that probability determines to leading order the frequency with which immunity protects against detectable reinfection given challenge.

Translating host immune histories into old antigenic variant and new antigenic variant neutralization probabilities for the analysis in Figure 3G,H requires a model of how susceptibility decays with antigenic distance, which we measure in terms of the typical distance between two adjacent ‘antigenic clusters’ (Smith et al., 2004). Figure 3 shows results for two candidate models: a multiplicative model used in a prior modeling and empirical studies of influenza evolution (Boni et al., 2006; Asaduzzaman et al., 2018), and a sigmoid model as parametrized from data on empirical HI titer and protection against infection (Coudeville et al., 2010).

In the multiplicative model, the probability z⁢(i,x) of no infection with variant i given that the nearest variant to i in the immune history is variant x is given by:

(40) z⁢(i,x)=z0⁢(z1z0)d⁢(i,x)

where z0 is the probability of no infection given homotypic reinfection, d⁢(i,x) is the antigenic distance in antigenic clusters between i and x, and z1 is the probability of no infection given d⁢(i,x)=1.

In the sigmoid model:

(41) z⁢(i,x)=1-11+eb⁢(ln⁡(T⁢(i,x))-a)

where a=2.844 and b=1.299 (α and β estimated in Coudeville et al., 2010) and T⁢(i,x) is the individual’s HI titer against variant i (Coudeville et al., 2010). To convert this into a model in terms of z0 and z1 we calculate the typical homotypic titer T0 implied by z0 and the n-fold-drop D in titer per unit of antigenic distance implied by z1, since units in antigenic space correspond to a n-fold reductions in HI titer for some value n (Smith et al., 2004). We calculate T0 by plugging z0 and T0 into Equation 41 and solving. We calculate D by plugging z1 and T1=T0⁢D-1 into Equation 41 and solving. We can then calculate T⁢(i,x) as:

(42) T⁢(i,x)=T0⁢D-d⁢(i,x)

The establishment of successful influenza virus infection requires access to respiratory epithelial cells covered by a mucosal layer, which acts as a protective barrier for virus entry. Sialic acid residues present in the mucus act as decoys for virions and the enzymatic activity of the influenza virus neuraminidase (NA) protein is essential for penetration of the mucosal layer (Matrosovich et al., 2004). The mucosal layer does not act only as a mechanical barrier. Secretory IgA (sIgA) mucosal antibodies specific to influenza virus proteins can apply virus-specific neutralization by immobilizing virions and slowing down the process of virus infection (Wang et al., 2017). Depending on their specificity, sIgA can either play a role in reduction of total virus population size or in providing specific competitive advantage to mutant viruses. Antibodies against the influenza virus NA protein will slow down the virus passage through the mucosal barrier irrespective of the HA antigenic phenotype of the newly infecting virus. By contrast, anti-HA antibodies are essential for inoculation selection because they recognise previously encountered antigenic variants. This provides a competitive advantage to new, distinct antigenic variants. The new variants with strongest competitive advantage are those with large antigenic effects; they have the lowest probability to be neutralized by any cross-reactive anti-HA sIgAs.

Virions that successfully cross the mucosal barrier can infect epithelial cells in the respiratory tract. Upon initiation of replication in host cells, newly generated virions can be subjected to replication selection via virus-specific antibodies. Replication selection can take place only if influenza virus-specific antibodies are present in mucosa-associated lymphoid tissue (MALT) at the time of virus replication.

Influenza virus-specific antibodies in the MALT are secreted by local plasmablast populations, generated following activation of B naive cells or B memory cells. The exposure history of the host determines the timing and the immunological trajectory for the generation of influenza virus-specific antibodies. In individuals without previous exposure to influenza virus, newly generated virions are recognised by naive B cells which undergo rounds of affinity maturation and somatic hypermutation to produce highly-specific antibodies to the virus. This process typically occurs in germinal centers, requires T-cell help, and takes around 7 days for production of highly-specific antibodies (Wrammert et al., 2008).

In individuals with previous exposure history to influenza viruses, where the newly generated influenza virions are antigenically matched or have some degree of cross-reactivity to previously encountered variants, the generation of influenza virus-specific antibodies mainly results from B memory recall response. The activation of pre-existing memory pools enables the immune system to respond quicker to previously encountered antigens without the need to recruit naive B cell populations. The generation of antigen-specific antibodies from recalled B memory pools does not typically start until 3–5 days post infection (Coro et al., 2006; Lam and Baumgarth, 2019).

The B naive and B memory cells generating serological responses to influenza viruses are referred to as B-2 cells. The activation of these cell types depends on the recognition of virus protein via the B-cell receptor and sufficient affinity of binding to trigger stimulatory signal. Depending on their specificity to the recognised virus antigen, the antibodies produced by these cell types have the potential to apply replication selection to any generated new antigenic variant viruses and to old variant viruses during the course of infection. Antibodies against influenza viruses can also be generated from the so-called B-1 cells, which do not require specific recognition of virus antigen and are activated by innate mechanisms (Lam and Baumgarth, 2019).

B-1 cells produce natural IgM antibodies with low affinity which constitute the earliest humoral response generated in the first 48–72 hr of infection (Lam and Baumgarth, 2019). As the activation of B-1 cells is independent of the specific recognition of influenza virus virions, the natural IgM antibodies act as a non-specific limiting factor for virus population sizes and cannot apply replication selection to particular virus antigenic variants.

The peak of influenza virus replication occurs in the first 24–48 hr post infection (Hadjichrysanthou et al., 2016). Despite the presence of several immunological routes for generation of influenza virus-specific antibodies, the timing of the serological immune response is always later than the time of peak of virus titer (at least in typical infections of otherwise healthy individuals), thus limiting the opportunity for replication selection during the course of infection.

The design of the inoculation selection model focuses on selection applied by the mucosal layer upon the entry of influenza virus in naïve or previously exposed hosts. We do not take into account the potential bottleneck applied during exit of viruses from the respiratory tract of an infected host. Selection of newly generated virions at the point of exit of an infected host is technically possible as the virions released from infected epithelial cells need to cross the mucosal barrier again to reach the lumen of the respiratory tract. Such selection upon virus exit would have the same directional effect as inoculation selection; it would provide a competitive advantage for mutant viruses with large antigenic effects. The potential strength of such selection depends heavily the how many sIgA antibodies in the transmitting host’s mucosa are free to neutralize the virions passing through; this quantity is unknown. Incorporating such mucosal selection upon virus exit thus requires better knowledge of the quantities of sIgA antibodies in given volume of mucus and what proportion of the available sIgA antibodies are engaged with antigen when an infected host transmits. Such data is not currently available and would require well-designed cellular culture models which take into account the role of mucosal layer in the process of virus infection and release of newly generated viruses.

According to our model, the strongest inoculation selection can take place in infection of previously exposed individuals and the highest probability of antigenic diversification occurs in populations with an intermediate proportion of immune hosts. In deriving this result, we have assumed homogeneous and static quantities of sIgA antibodies among previously exposed individuals. This is an oversimplification. There is substantial variation in level of sIgA antibodies across individuals depending on the nature of their previous exposures (through vaccination or natural infection). For example, intranasal inoculation with influenza virus via live attenuated influenza virus vaccines leads to generation of better mucosal immunity (Hoft et al., 2017) than intramuscular administration of trivalent inactivated vaccines and thus previously exposed individuals are likely to vary in the degree of inoculation selection applied during homotypic infections depending on the route of administration of previous vaccinations.

The model also does not account for waning of sIgA antibody titers in the months and years following each exposure. The waning of sIgA titers is likely to be substantial in absence of re-exposure but the waning process should only decrease the efficiency of inoculation selection and make the accumulation of population level selection pressure—and thus the appearance of observable new antigenic variants—rarer.

The dynamics of mucosal immunity with age are also likely to have impact on the selector phenotype of previously exposed hosts. In older individuals, multiple previous exposures will eventually lead to preferential activation of recall responses for antibody production to new variants as a result of original antigenic sin (Davenport and Hennessy, 1956), immune backboosting (Fonville et al., 2014), or antigenic seniority (Lessler et al., 2012).

Even if the recall responses result in backboosting and higher levels of sIgA, these antibodies are unlikely to be well-matched to the new variant or succeeding variants, as continuous re-exposure favors responses to conserved epitopes (Krammer and Palese, 2019; Krammer, 2019) rather than receptor-binding site epitopes suspected to be most important for immune escape.

As a result of the limited ability to generate novel immune responses and the phenomenon of antigenic seniority, an elderly individual who has been exposed multiple times in the past is likely to exert very weak inoculation selection to newly infecting viruses. By contrast, a young person recently exposed to influenza virus is more likely to act as a strong selector, due to the availability of mucosal sIgA antibodies and the more limited impact of original antigenic sin. However, due to the acute nature of influenza virus infection, very young children are likely to require more than one exposure before they begin to develop highly-specific antibody responses to influenza virus infection able to exert inoculation selection (Neuzil et al., 2006).

Most current efforts for characterizing humoral immunity to influenza virus infections are focused on antibody responses in the serum. However, there is limited understanding of the relationship between antibodies in the serum measured via HI and the levels of mucosal antibodies able to exert inoculation selection. Better characterization of the dynamics of mucosal immunity with age and across individuals is needed to understand the implications of inoculation selection to influenza virus evolution depending on the host population structure.

In our minimal model, competition between new variant and old variant viruses obeys a simple replicator equation.

The within-host frequency of new variant virus is fm⁢(t)=Vm⁢(t)Vm⁢(t)+Vw⁢(t)=Vm⁢(t)Vtot⁢(t). Its rate of change is d⁢fmd⁢t. We observe that d⁢Vid⁢t is of the form d⁢Vid⁢t=Vi⁢[gi⁢(t)-di⁢(t)] for both the new variant and the old variant, where, neglecting mutation, gi⁢(t)=ri⁢β⁢C⁢(t) and di⁢(t) is a (possibly time-varying) neutralization or decay rate. Differentiating fm⁢(t) with respect to time demonstrates that the frequency of the new variant over time is governed by the following replicator equation (see Section A9 for an explicit derivation):

(A1) d⁢fmd⁢t=fm⁢(1-fm)⁢([gm⁢(t)-gw⁢(t)]-[dm⁢(t)-dw⁢(t)])

If the new variant is neutral in the absence of antibody-mediated selection, rw=rm=r, and therefore gm⁢(t)=gw⁢(t). It follows that despite changing target cell populations and virion population sizes, the variants compete in a manner independent of Vtot, provided each di⁢(t) is independent of Vtot. We have competition only through the respective death rates of old variant and new variant virions. These may be affected by antibody-mediated neutralization of virions.

If the fitness difference δ⁢(t)=[gw⁢(t)-gm⁢(t)]-[dm⁢(t)-dw⁢(t)] is a constant δ, this differential equation has a closed-form solution:

(A2) fm⁢(t)=eδ⁢teδ⁢t+f0-1-1

where f0=fm⁢(0)

If δ=0, fm⁢(t) is constant and equal to f0, as expected. If δ⁢(t) is piecewise-constant, we can apply the constant-δ solution iteratively to find the new variant frequency.

It is therefore straightforward to apply this expression to the within-host model with binary immunity by evaluating different immunity regimes piece-wise.

We consider the case of a host who is experienced to the old variant (Ew=1) but naive to the new variant (Em=0).

If first appears at a time te≥0 post-infection at an initial frequency fm⁢(te), then if te<tM:

(A3) fm(t)={0t<tefm(te)te≤t<tMeδ(t−tM)eδ(t−tM)+fm(te)−1−1tM≤t

if te>tM:

(A4) fm(t)={0t<teeδ(t−te)eδ(t−te)+fm(te)−1−1te≤t

where δ=k⁢(cw-cm) is the fitness difference between the variants due to the recall adaptive response.

Thus far we have neglected mutation the first mutation of interest (i.e. the first that produces a new variant lineage that evades stochastic extinction). In fact, with symmetric mutation between the two types of interest, we have:

(A5) dVwdt=Vw[(1−μ)gw(t)−dw(t)]+μgm(t)VmdVmdt=Vm[(1−μ)gm(t)−dm(t)]+μgw(t)Vw

The rate of fitness-irrelevant mutation is assumed to be equal for the two types, and to preserve antigenic type (w or m). The rate of lethal or deleterious mutants is assumed equal for the two types and is captured in g.

It can be shown (see Section A9) that the equation for d⁢fmd⁢t in this case is:

(A6) d⁢fmd⁢t=fm⁢(1-fm)⁢(αm-αw)+μ⁢(gw⁢(t)-2⁢gw⁢(t)⁢fm+gw⁢(t)⁢fm2-gm⁢(t)⁢fm2)

When gw⁢(t)=gm⁢(t)=g⁢(t):

(A7) d⁢fmd⁢t=fm⁢(1-fm)⁢(dm⁢(t)-dw⁢(t))+μ⁢g⁢(t)⁢(1-2⁢fm)

When fm≪0.5, the net effect of symmetric mutation on fm is to increase it at rate ≈μ (note that the symmetric mutation term becomes zero at fm=1/2). In other words, it is roughly equivalent to a case with no back mutation at all:

(A8) dVwdt=Vw[(1−μ)gw(t)−dw(t)]dVmdt=Vm[gm(t)−dm(t)]+μgw(t)Vw

in that case, we have:

(A9) dfmdt=fm(1−fm)([gm(t)−(1−μ)gw(t)]−[dm(t)−dw(t)])+μgw(t)(1−2fm+fm2)

In either case, we mainly study a situation in which gw⁢(t)=gm⁢(t)=g⁢(t) at all times (though note that g⁢(t) varies in time), dw⁢(t)=dm⁢(t)=d⁢(t)=dv prior to the onset of the adaptive immune response, and dw⁢(t)=dv+k, dm(t)=dv+cmk after the onset of the adaptive immune response.

If fm is large relative to µ, the mutant growth term Vm⁢g⁢(t) is large relative to the de novo mutation term μ⁢g⁢(t)⁢Vw. We can see this from the fact that Vm=fm⁢Vtot≥fm⁢Vw, with equality if and only if fm=0. It follows that if fm>μ>0, then:

g⁢(t)⁢Vm>g⁢(t)⁢fm⁢Vw>g⁢(t)⁢μ⁢Vw

But if the first mutation to produce a new variant is sufficiently late and there is little or no positive selection, we do not necessarily have fm≫μ, so ongoing mutation cannot be ignored in the dynamics of fm. In that case, we can consider either the μ⁢g⁢(t)⁢(1-2⁢fm+fm2) term in Equation A9 or the μ⁢g⁢(t)⁢(1-2⁢fm) term in Equation A7. Denote the new variant frequency at the time of first mutation te by f0=fm(te). We can approximate g⁢(t)≈g0=g⁢(0)=ℛ0⁢dv early in infection, since for small t, C⁢(t)≈Cmax and therefore g⁢(t)=r⁢β⁢C⁢(t)≈g0.

With this approximation for g⁢(t), the one-way mutation case has an analytical solution for δ≠0:

(A10) fm(t)≈Aexp⁡(δt)+μg0Aexp⁡(δt)+μg0−δA=f01−f0(μg0−δ)−μg0f0

And when δ=0:

(A11) fm⁢(t)≈11-f0+μ⁢g0⁢t-111-f0+μ⁢g0⁢t

fm must be small for ongoing mutation to matter; it ceases to matter when fm≫μ⁢g⁢(t), and μ⁢g⁢(t) is small by assumption. It follows that in the absence of a fitness difference between the types:

and therefore:

We can also see this by inspecting the other two expressions for fm⁢(t) and seeing that they are quasi-linear for small µ and small t. This further confirms that it is not crucial to decide whether symmetric or one-way mutation is more realistic, as they will be functionally the same from the point of view of practical mutant frequencies in the absence of positive selection.

These approximations break down once C⁢(t)≪Cmax and therefore gw⁢(t)≪g0; fewer wild-type replications are occurring, and rate of ongoing mutation therefore falls. So a reasonable approximation for the maximum value of fm⁢(t) in the absence of positive selection is given by fm⁢(tpeak). For a mutation rate of 0.33×10-5, typical values of fm⁢(tpeak) without selection range from 1×10-5 to 2×10-4, depending on ℛ0 and tpeak.

We use the expression in Equation A10 to plot the analytical predicted mutant frequencies shown in Figure 1 and to calculate the transmission frequencies for the analytical model of observed phenotypes in Figure 1.

Ongoing mutation enhances possibilities for both replication and inoculation selection. It truncates the left tail of the distribution of fm⁢(tM) and the distribution of fmt (the new variant frequency at the time of transmission). Whenever tM is sufficiently late that we have reached Vw≪1μ before tM, fm⁢(tM) should be on the order of the inverse mutation rate, or larger.

The consequence is that probabilities of replication selection and inoculation selection should both be somewhat higher than those estimated based on the time to the first non-extinct de novo mutant lineage, as in Figure 1. This further strengthens the case for the importance of an adaptive response at 48 hr or later in explaining the absence of replication selection.

For a fixed ℛ0, replication selection becomes more likely if the virions produced per infected cell r becomes large—that is, if every infected cell produces more individual virions.

For a given dv, and f, a virus can achieve a large ℛ0 either by having a high cell productivity r or by having a high cell infection rate β.

The expected number of virus replications per lost target cell is given by V¯˙+C˙−=rβCV¯ℓβCV¯=rℓ. The infection peaks and begins to decline when d⁢Vd⁢t=0, which occurs when g⁢(t)=d⁢(t), or CmaxC=ℛ0, provided d⁢(t) is fixed. It follows that Cpeak=Cmaxℛ0. So the virus peaks once Cmax⁢(1-1ℛ0) target cells have been consumed (assuming negligible target cell replenishment during the timecourse of infection), and rℓCmax(1−1ℛ0) viral replication events have occurred (this is a slight overestimate due to target cell depletion).

Fixing dv and ℛ0, as r→∞, β→0 and the number of viral replication events before peak viral load goes to ∞. Since each generated mutant survives stochastic extinction roughly proportional to 1/ℛ0, these earlier mutants also are lost stochastically less frequently. In other words, a mutant that survives stochastic extinction is generated with certainty well before the infection peaks and can spend arbitrarily long under selection. Conversely, as r→0, β→∞, the number of replication events that occur before the infection peaks goes to 0, and replication selection becomes impossible.

A virus of a given ℛ0 that takes a shotgun approach of producing many not-especially-infectious virions is more likely to result in the proliferation of a mutant of interest under replication selection than one that produces a small number of highly infectious virions.

In other words, the relatively large infected cell productivity of influenza viruses—perhaps as large as hundreds or thousands of viable virions per cell (Frensing et al., 2016)—makes potential replication selection more efficient. This in turn makes the absence of observed replication selection harder to explain unless another mechanism can be invoked, such as inoculation selection with replication selection only occurring 2–3 days post-infection.

One special case bears mentioning, because it corresponds to multiple existing models (Luo et al., 2012; Volkov et al., 2010; Kennedy and Read, 2017) and may be relevant for other systems. If there is an immediate recall response (tM=0), a founding population of size b, and k is large enough such that ℛw(t=0)<1 (the virus population is initially declining, in expectation), then the infection will go extinct after generating q new virions for some finite q.

Each virion has a probability µ of being an escape mutant with ℛm⁢(0)>1 and each mutant independently survives stochastic extinction with probability psse=(1-1/ℛm⁢(0)) (per the theory of supercritical branching processes). The conditional probability of replication selection given q is therefore:

(A14) prepl(q)=1−(1−μpsse)q

if μ⁢psse is small, this is approximately equal to:

(A15) prepl⁢(q)≈1-e-q⁢μ⁢psse≈q⁢μ⁢psse

The approximations use the Poisson approximation to the binomial and the approximation ex≈1+x when |x|≪1.

The unconditional probability of replication selection prepl is the expectation in q of prepl(q): Eq(prepl). But since prepl is approximately linear in q, the linearity of expectation implies:

prepl≈Eq⁢(q⁢μ⁢psse)=q¯⁢μ⁢psse

where q¯ is the expected value of q.

Given ℛw⁢(0)<1, branching process theory predicts that the virus will produce on average q¯=b1-ℛw⁢(0)-b new copies before the infection is cleared (expected total length of b independent subcritical branching processes, each with expected length 11-ℛw⁢(0), minus the initial virions present [Ball et al., 2016]). This can be simplified to q¯=b⁢ℛw⁢(0)1-ℛw⁢(0).

So we have

(A16) prepl≈q¯μpsse=ℛw(0)1−ℛw(0)bμpsse

Note that Iwasa et al., 2004 have previously derived this approximate result for the probability that a subcritical replicator mutates to a supercritical one prior to going extinct in the context of analyzing tumor cell mutations. They use a more rigorous generating function approach.

In our model, mucosal antibodies acting at the point of transmission provide a mechanism for population level selection pressure: some protection of experienced hosts against infection with old variant virus, and worse protection against infection with new variant virus. In particular, it provides a mechanism that produces selection pressure for new antigenic variants while predicting that reinfections with old variant viruses—without observable new variant viruses—should still be observed. Prior models (see Section A7) predict that in fully immune hosts, any observable reinfections will have new antigenic variant viruses at consensus.

Antibodies at the point of transmission appear to play a key role in population-level influenza virus evolution: they produce the selection pressure that allows new variants to spread more reliably and thus rapidly from host to host than old variants. But they may play a second role as well: promoting of new variants in frequency from the low frequencies at which they are typically generated to high frequencies at which they can be observed and reliably transmitted onward. This inoculation selection takes on the role of previously held by replication selection in explaining why inoculations of experienced hosts might produce new variant infections at a higher rate (per inoculation) than inoculations of naive hosts.

As noted in the main text, new variants can reach high within-host frequencies through founder effects. These events are both possible and rare due to influenza’s tight transmission bottleneck (McCrone et al., 2018; Xue and Bloom, 2019).

Whether this process is purely neutral or stochastically selective depends on whether the recipient host is naive, and if not how well they neutralize old variant versus new variant virions.

As mentioned in the Methods, in a fully naive host with large v and small transmitting host new variant frequency, the sampling process approximates a low-probability binomial or a low-frequency Poisson:

When fmt≪1, this is approximately equal to b⁢fmt

As noted in the main text Materials and methods, mucosal antibody neutralization can reduce the new variant’s survival probability relative to this neutral case (inoculation pruning) or increase it (inoculation promotion), depending upon parameters.

There can be inoculation pruning even when the new variant is more fit than the old variant (i.e. neutralized with lower probability). But sufficient subsequent bottlenecking after sIgA neutralization can create an inoculation promotion effect.

Whether inoculation selection produces pruning or promotion depends on the ratio of v/b, and on the mutant’s probability of avoiding neutralization 1-κm (main text Figures 4F,5, Appendix 1—figure 1, and Section A4.5 below).

As shown in the Materials and methods, the probability that a new variant survives the cell infection bottleneck of size b is well approximated for small fmt by:

(A18) psurv(κm,κw,v,b,fmt)≈pinoc(κm,v,fmt)∗pcib(κw,v,b)

pinoc is the probability that at least one mutant survives the sIgA bottleneck:

(A19) pinoc⁢(κm,v,fmt)=1-e-v⁢fmt⁢(1-κm)

Note that when fmt≪1,pinoc≈v⁢fmt⁢(1-κm)

pcib⁢(κw,v,b) is the probability that a mutant that survives the sIgA bottleneck is not lost at the final bottleneck (as given in Equation 35 above):

pcib⁢(κw,v,b)=(1-e-w¯)⁢bw¯+e-w¯⁢∑j=0b-1w¯jj!⁢(1-bj+1)

where w¯=v⁢(1-fmt)⁢(1-κw). If b=1, this expression reduces to the first term.

For given values of the other parameters, larger fmt implies larger psurv. Since both pinoc and pcib are increasing in fmt and are both always positive, it follows that psurv=pinoc⁢pcib is also increasing in fmt. This makes intuitive sense: all else equal, higher frequency mutants always have a better chance of surviving at the point of transmission, regardless of the effects of the sIgA bottleneck.

Consider psurv/ pdrift. This ratio quantifies the degree to which sIgA neutralization at the point of transmission facilitates (or hinders) mutant survival of the final bottleneck relative to a purely neutral founder effect.

For realistically small b and fmt, we can use the approximations pdrift≈b⁢fmt and pinoc≈v⁢fmt⁢(1-κm)

Then we have

(A22) psurv/pdrift=pinoc⁢pcibpdrift≈vb⁢(1-κm)⁢pcib⁢(κw,v,b)

Several intuitive results follow:

• If we hold κw constant, increasing κm makes inoculation selection less effective relative to drift, because the inoculated new variant is at greater risk of neutralization.

• Conversely, if we hold κm constant, increasing κw makes inoculation selection more effective relative to drift (this follows from the fact that pcib is increasing in κw, see Section A4.3.1 above).

• Since pcib is increasing in fmt (Section A4.3.1), psurv/pdrift is increasing in fmt. That is, inoculation selection is more efficient relative to drift when promoting higher frequency mutants, all else equal.

• psurv/pdrift is decreasing in b (see Section A9.9 for a derivation). In other words, when the cell infection bottleneck is wider, inoculation selection is less important relative to drift. The large bottleneck gives the mutant a reasonably good chance of surviving even in the absence of any neutralization of competing old variant virions.

For a bottleneck of size b=1, we can also see that inoculation selection becomes more efficient relative to drift as v increases, since:

psurv/pdrift≈vb(1−κm)pcib=v(1−κm)1−exp⁡(−v(1−fmt)(1−κw))v(1−fmt)(1−κw)=11−fmt1−κm1−κw[1−exp⁡(−v(1−fmt)(1−κw))]

This is increasing in v.

In a host strongly immune to the old variant (large κw), we have -v⁢(1-fmt)⁢(1-κw)≪1. Applying the linear approximation to the exponential to equation A4.5:

psurv/pdrift≈11−fmt1−κm1−κw(v(1−fmt)(1−κw))=(1−κm)v

This provides intuition for the effects seen in main text Figure 4, in which inoculation selection improves survival relative to drift by a factor of vb for a full escape mutant, and by a factor scaled down by 1-κm for a partial escape mutant.

Moreover, if we let fmt=fd for the drift case and fmt=fr for the selective case, the approximate ratio becomes

Our replicator equation implies that this ratio frfd will increase in expectation when more replication selection occurs prior to transmission (larger δ⁢τ), explaining the increasing ratio of the selective variant survival relative to drift observed in Figure 4 when replication selection is permitted prior to transmission.

Note that we did not integrate over the emergence times to obtain this ratio; instead, we derived principles that take the initial frequency as a given.

We also wish to know which hosts best promote antigenic novelty when inoculated, given a level of cross-immunity σ, or a degree of escape 1-σ. We will show that whereas replication selection suggests that intermediately immune hosts should be key selectors for antigenic novelty (Grenfell et al., 2004), in an inoculation selection regime, new variants may most often reach observable frequencies in fully immune hosts. We do this by analyzing an expression for psurv, the probability that a mutant survives transmission to found an infection in a new host, and showing how it changes with old variant neutralization probability κw and sIgA cross immunity σ=κm/κw. We find that transmission survival probability can be maximized in strongly immune or fully naive hosts, depending on parameters.

We optimize psurv⁢(κw) given some value of σ by differentiating.

(A23) dpsurvdκw=dpcibdκwpinoc+pcibdpinocdκw=dpcibdκwpinoc+pcib(vfmtσ)(pinoc−1)

Consider the endpoint at κw=1. pcib=1, pinoc≈fmt⁢v⁢(1-κm)=fmt⁢v⁢(1-σ), and w¯=v⁢(1-fmt)⁢(1-κw)=0.

So dpcibdκw=dpcibdw¯dw¯dκw=dpcibdw¯(−v)(1−fmt) can be evaluated (using the expressions for dpcibdw¯ derived in section A9.7 below), and it evaluates to 0 except if b=1 (when it is equal to 12v(1−fmt)).

Case 1: b>1. For b>1, it follows from the above that

(A24) dpsurvdκw=0−pcibdpinocdκw

Since d⁢pinocd⁢κw≥0, d⁢psurvd⁢κw<0 unless d⁢pinocd⁢κw=0, which occurs in cases of interest when σ=0 (complete escape mutant).

In other words, if immune escape is incomplete for b>1, there is always some less strongly immune host (though possibly very slightly less, see Appendix 1—figure 1) that improves new variant survival relative to a host who neutralizes old variant virions with 100% certainty. Note, however, that hosts with κw=1 are very unlikely to exist in nature. Our study is motivated precisely by the fact that even an experienced host encountering homotypic reinfection neutralizes virions with probability κw<1, and therefore can be productively reinfected (Clements et al., 1986; McCrone et al., 2018). It follows that in practice the most strongly immune hosts (with κw large but less than 1) could still be the best selectors.

Case 2: b=1. When the bottleneck is 1 and w¯=0, dpcibdκw=12v(1−fmt). So in this case, the derivative at κw=1 may be positive or negative, depending on whether:

(A25) d⁢pcibd⁢κw⁢pinoc+pcib⁢(v⁢fmt⁢σ)⁢(pinoc-1)>0

substituting dpcibdκw=12v(1−fmt) and pcib=1, we have:

12(1−fmt)pinoc+fmtσpinoc>fmtσ

12fmt(1−fmt)pinoc+σpinoc>σ

12fmt(1−fmt)pinoc>σ(1−p inoc)

When pinoc is small, we have approximately pinoc≈fmt⁢v⁢(1-σ) and 1-pinoc≈1, so:

In other words, with extremely small cell infection bottlenecks like those observed for influenza viruses, fully immune hosts are the best selectors provided that the number of virions v encountering IgA antibodies is sufficiently large given the degree of escape achieved. Less immune escape (larger sIgA cross immunity σ) necessitates a larger v to make fully immune hosts the best selectors. Since fmt≪1, a rule of thumb is that v>2⁢σ1-σ.

The intuition is that larger v means more competition for the cell infection bottleneck among virions that reach IgA, and thus a greater opportunity for the mutant’s selective advantage to be realized, but that this only works provided that this advantage is large enough so that the mutant is not itself at too large a risk of being neutralized.

A crucial parameter in our model is δ: the magnitude of the fitness advantage of the new variant over the old variant during viral replication in an infected individual. One possible objection to our analysis here is that the antibody-mediated virion neutralization rate k is low enough or the antigenic similarity between the variants is great enough to make the fitness difference δ=k⁢(cw-cm) small. In that case, replication selection to consensus will be rare even if tM is very small, and the adaptive response is mounted immediately (main text Figure 1C,D). Despite uncertainty about both k and cw-cm, k is likely to be high, perhaps extremely so. And given sufficiently high k and a homotypically reinfected host (cw=1), even moderate immune escape (0≪cm<1) produces a substantial fitness difference. Moreover, assuming small k early in infection grants our basic hypothesis: antibodies mediated selection is weak early in infection, neutralization during viral replication is not the mechanism of protection against reinfection, and there is asynchrony between antigenic diversity and meaningful antigenic selection.

As noted in the Discussion, the relationship between the mucosal antibody neutralization rate κ and the antibody neutralization rate during replication k is unknown, though we have every reason to expect it to be positive.

κ values are particularly difficult to calculate. In our model, we have generally assumed that all virions of type i inoculated into an experienced host are independently neutralized with probability κi, and therefore zi=e-v⁢(1-κi) for an inoculum consisting only of type i. But this assumption of independence may be violated in practice.

One reason statistical independence may be violated is that antibody numbers are finite, and an antibody that binds to one virion cannot bind to another. Consider a focal virion. Given that another virion has been neutralized, there are fewer antibodies remaining to neutralize our focal virion. This is particularly important for mixed inocula. Old variant virions, which have higher affinity for the inoculated host’s sIgA antibodies, may indirectly protect the new variant by competing with it for antibody-binding. A new variant virion may have a higher individual chance of being neutralized by those same antibodies if it is part of a monomorphic inoculum composed of other, identical new variants. Such an interaction would strengthen inoculation selection relative to the independent neutralization we have modeled here.

Inoculation selection may improve mutant bottleneck survival relative to neutral drift (see Figure 4F, Appendix 1—figure 1, and Section A4). For this to be true, a non-antigenic bottlenecking must follow IgA antibody neutralization, with a ratio v/b≫1

Evidence suggests that a non-antigenic bottleneck does occur after the sIgA bottleneck because bottlenecks measured in vaccinated and non-vaccinated hosts are of comparable size (indistinguishable from one), suggesting that founding virion population sizes are cut down to small numbers even in the absence of IgA antibodies, though this would also be consistent with the case of v=1 (i.e. IgA antibodies must neutralize on average a single virion to prevent infection, and this virion, if not neutralized or otherwise lost, uniquely founds the infection).

An experimental evolution study of avian influenza virus adaptation in ferrets found that transmission bottlenecks in naive hosts became tighter that as the virus adapted (Moncla et al., 2016). A re-analysis of that data found that bottlenecks were tight throughout (Lumby et al., 2018). The fact that bottlenecks do not appear to loosen (and may tighten) with adaptation for better transmission and replication is further evidence, albeit circumstantial, that when an influenza virus is well-adapted to its host and replicates rapidly, the first virion or first few virions to infect a cell will be the ancestor of the vast majority of progeny viruses.

One modeling study of influenza viruses proposes that infections should have a second peak after initial innate responses are overcome and some target cells are once again susceptible to infection (Pawelek et al., 2012). A relaxation of non-antigenic limiting factors later in infection can and should provide additional opportunities for replication selection, as occurs in immune-comprised human patients.

That said, the empirical data suggesting double-peaked kinetics comes from experimental inoculation of naive horses with a large quantity of influenza virus: 106 50% egg infectious dose (Quinlivan et al., 2007). The novel antibody response curbed the second peak of infection, which was much lower than the first.

In human kinetics data (Hadjichrysanthou et al., 2016) and animal transmission experiments in which one animal infects another (Le Sage et al., 2020; Canini et al., 2020) it is common to see clearly single-peaked infections.

Furthermore, for realistic (incomplete) degrees of antibody-binding escape, we expect both old and new antigenic variant population sizes to decline even due to antigenic limiting factors, though of course we expect the new variant to decline more slowly.

The upshot is that virus clearance should be rapid following the mounting of the adaptive response in the experienced hosts where we expect selection to take place, and thus in immune-competent hosts we should expect a limited window for antibody-mediated selection on transmissible virus populations.

In the simulation results shown in the main text (Figures 1,2,4,6) we use a probabilistic model of transmission: the donor host’s probability of inoculating the recipient at a given point during the infection is proportional to donor viral load in virions, Vtot⁢(t).

Since influenza virus population sizes grow and shrink rapidly around the within-host peak, however, this model should be qualitatively similar to a threshold model in which transmission occurs with certainty if the total viral load is above a certain threshold θ and does not occur if it is below that threshold. To confirm this, we here show corresponding transmission chain simulation results based on a threshold model, with the threshold as in Table 1.

Probability shown as a function of probability of no old variant infection (zw), degree of cross immunity between mutant and new variant σ=κm/κw, mucus bottleneck size v, and final bottleneck size b. Gray dotted line indicates probability that a new variant survives the transmission bottleneck in a host who is naive both to the old variant and to the new variant (i.e. drift). fmt=9×10−5, a typical value for a naive transmitting host in stochastic simulations.

Threshold version of Figure 4. Distribution of antigenic changes along 1000 simulated transmission chains (A–D) and from an analytical model (E–H). In (A,E) all naive hosts, in other panels a mix of naive hosts and experienced hosts. Antigenic phenotypes are numbers in a 1-dimensional antigenic space and govern both sIgA cross immunity, σ, and replication cross immunity, c. A distance of ≥1 corresponds to no cross immunity between phenotypes and a distance of 0 to complete cross immunity. Gray line gives the shape of Gaussian within-host mutation kernel. Histograms show frequency distribution of observed antigenic change events and indicate whether the change took place in a naive (gray) or experienced (green) host. In (B–D) distribution of host immune histories is 20% of individuals previously exposed to phenotype −0.8, 20% to phenotype −0.5, 20% to phenotype 0 and the remaining 40% of hosts naive. In (E), naive hosts inoculate naive hosts. In (F–H) hosts with history −0.8 inoculate hosts with history −0.8. Initial variant has phenotype 0 in all sub-panels. Model parameters as in Table 1, except k=25. Spikes in densities occur at 0.2 as this is the point of full escape in a host previously exposed to phenotype −0.8.

(A, B) Probability of a detectable infection per inoculation of an experienced host. (C, D) Probability of a detectable new variant infections per inoculation of an experienced host. (E, F) Fraction of detectable infections of experienced hosts that are new variant infections. Points colored according to whether new variant infections were more frequently caused by de novo generated (purple) or inoculated (green) new variant viruses. Each point represents a random parameter set; 10 random parameter sets generated for each bottleneck value shown, and 50,000 inoculations of an experienced host simulated for each parameter set. Two regimes of simulated parameter set shown: (A, C, E) an immediate recall response regime, in which tM varied from 0 to 1, and (B, D, F) a realistically-timed response regime, in which tM varied from 2 to 4.5. All other parameters varied across the same ranges in both regimes (see Table 2 for ranges). Parameters Latin Hypercube sampled from within ranges for each regime and bottleneck size.

As described in the Materials and methods, we further assessed sensitivity of the within-host model to variation in key parameters by studying random parameter sets chosen from biologically plausible parameter ranges.

We analyzed two cases: one in which the immune response is unrealistically early (0≤tM≤1, Appendix 1—figure 4), and one in which it is realistically-timed (2≤tM≤4.5) Appendix 1—figure 5, see also Appendix 1—figure 3 and other model parameters in Table 2.

We found that with early immunity, regardless of particular parameter values, new variants are frequently seen when experienced hosts are detectably reinfected, and the overall probability of new variant infections is unrealistically high. These observable new variants are most often generated de novo and replication-selected (Appendix 1—figure 4, see also Appendix 1—figure 3).

When immunity is realistically-timed, the pattern of much rarer replication selection is robust to variation in parameters. New variant infections compose a realistically small fraction of all detectable reinfections (Appendix 1—figure 5, Appendix 1—figure 3). Moreover, inoculation selection is sometimes more common than replication selection as a source of new antigenic variants (Appendix 1—figure 5, see also Appendix 1—figure 3).

Parameters randomly varied across the ranges given in Table 2, with tM varied between 0 and 1. Each point represents a parameter set; the rate of new variant infections per hundred 100 detectable infections is estimated from 50,000 simulated inoculations of an experienced host. A new variant infection is defined as one in which the new variant reached a transmissible frequency of at least 1% at any point in the infection. Points are colored according to whether new variant infections were more frequently caused by de novo generated (purple) or inoculated (green) new variant viruses.

Parameters randomly varied across the ranges given in Table 2, with tM varied between 2 and 4.5. Each point represents a parameter set; the rate of new variant infections per hundred 100 detectable infections is estimated from 50,000 simulated inoculations of an experienced host. A new variant infection is defined as one in which the new variant reached a transmissible frequency of at least 1% at any point in the infection. Points are colored according to whether new variant infections were more frequently caused by de novo generated (purple) or inoculated (green) new variant viruses.

Amino acid substitutions associated with antigenic “cluster transitions” have reached fixation in parallel—and near-synchronously in time—in genetically distinct co-circulating lineages. This has occurred both in A/H1N1 and in A/H3N2 seasonal influenza viruses.

Existing 'mutation-limited' or 'diversity-limited' model explanations of why observable influenza virus antigenic novelty is rare despite continuous genetic evolution are less convincing in light of this synchronous polyphyly (see Section A7). In particular, a number of models hypothesize that large effect antigenicity-altering substitutions can only emerge in specific genetic contexts (Koelle et al., 2006; Gog, 2008; Koelle and Rasmussen, 2015; Kucharski and Gog, 2012), and that this constraint limits the rate of population-level antigenic change.

Synchronous, polyphyletic cluster transitions cast doubt on these claims. If major antigenic changes appeared infrequently due to the rarity of 'jackpot' scenarios of a large-effect mutant in a favorable genetic background (Koelle and Rasmussen, 2015), it would be surprising to see two or three distinct and synchronous jackpots after years of none. So while genetic background may play an important role in determining which lineages are fittest, synchronous polyphyly suggests that it is unlikely to set the clock of antigenic evolution. Rather, it suggests that accumulating population immunity may be crucial in setting the clock, and that there may even be a rising generation rate of observable antigenic novelty as a cluster ages, rather than a fixed rate (see main text Section "Population immunity sets the clock of antigenic evolution").

To show that large effect antigenic changes occur near-synchronously in multiple backgrounds, we assessed evolutionary relationships and built phylogenetic trees for known recent polyphyletic antigenicity-altering substitutions that have emerged in A/H1N1 and A/H3N2. Our analysis rules out reassortment as an explanation for the apparent polyphyly, thus verifying that the branches represent distinct, independent, but simultaneous de novo lineages.

We analyzed polyphyletic antigenic changes using hemagglutinin (HA) gene nucleotide sequences deposited in the GISAID EpiFlu database.

We downloaded all seasonal A/H1N1 HA sequences from the period 1999–2008, to center on the antigenic cluster transition from the New Caledonia/1999-like phenotype to the Solomon Islands/2006-like phenotype, excluding pandemic A/H1N1 viruses. We downloaded A/H3N2 virus HA sequences for two periods: 2008 to 2011, to center on cluster transition from Brisbane/07 to the Perth/2009-like and Victoria/2009-like phenotypic split, and 2012 to 2014, to capture the co-circulation of Clade 3C.2a and 3C.3a (Appendix 1—table 1). We discarded all sequences with an incomplete HA1 domain or with more than 1% ambiguous nucleotides.

We aligned sequences using MAFFT v7.397 (Katoh et al., 2002) and reconstructed phylogenetic trees using RAxML 8.2.12 under the GTRGAMMA model (Stamatakis, 2014). We performed global optimization of branch length and topology on the RAxML reconstructed tree using Garli 2.01 with model parameters matching the RAxML reconstruction (500,000 generations) (Bazinet et al., 2014). It was computationally intractable to run Garli on the full phylogeny of A/H3N2 2012–2014 (n=10107). We visualized phylogenetic trees using Figtree (http://tree.bio.ed.ac.uk/software/figtree) and ggtree (Yu et al., 2017)

We mapped amino acid substitutions onto branches using custom Python scripts. All numbering complies with H1/H3 numbering scheme (Burke and Smith, 2014).

The HA K140E antigenicity-altering substitution first occurred in 2000 and was detected sporadically prior to its independent fixation in 2006–2007 in three phylogenetically distinct, geographically segregated lineages of A/H1N1 that diverged in 2004 (Appendix 1—figure 6; Bedford et al., 2015). The K140E substitution resulted in a cluster transition from the New Caledonia/1999-like antigenic phenotype to the Solomon Islands/2006-like phenotype, with lineage A emerging as the dominant lineage globally before the 2009 H1N1 pandemic (Bedford et al., 2015). Position 140 is located in the Ca2 antigenic site of the HA1 domain immediately adjacent to the receptor binding site, where amino acid composition mediates receptor binding function (Koel et al., 2013).

The three lineages have distinct HA1 mutational trajectories away from their most recent common ancestor, as defined by the set of amino acid substitutions that accumulate along the predominant trunk lineages both prior to and following the fixation of the K140E substitution (Appendix 1—table 2, Appendix 1—figure 6). All three lineages acquired amino acid substitutions at previously characterized antigenic sites, as well as substitutions with suggested functional consequences, including glycosylation gain and loss (Caton et al., 1982). The three lineages share the Y94H substitution prior to K140E fixation, with limited overlap of other substitutions: lineage A and B share changes at residue 188, whereas A189T occurs in lineage B and A prior to and post K140E emergence respectively.

Branch tip color indicates the amino acid identity at position 140. Co-circulating lineages defined by the K140E fixation are highlighted. Scale bar indicates the number of nucleotide substitutions per site. Tree rooted to A/New Caledonia/20/1999.

In the period 2008-2011, the K158N substitution was only detected once in 2009 before fixing in combination with N189K in the same year in two distinct co-circulating A/H3N2 lineages. N189K was not detected before fixing as a K158N/N189K double substitution. The combination of the substitutions resulted in an antigenic phenotype switch from Wisconsin/2005-like to the Perth/2009-like and Victoria/2009-like phenotypes respectively. Residues 158 and 189 are both located in antigenic site B adjacent to the receptor binding site, with substitutions at both residues characterized as cluster-transition substitutions in multiple historic A/H3N2 cluster transitions (Koel et al., 2013). The HA-defined evolutionary trajectories of viruses from the Perth/2009-like lineage and Victoria/2009-like lineages from their most recent common ancestor are characterized by distinct sets of substitutions (Appendix 1—table 3, Appendix 1—figure 7). Both lineages acquired substitutions at previously characterized antigenic sites (Koel et al., 2013), but only share the T212A substitution.

Branch tip color indicates the amino acid identity at position 158 and 189. Co-circulating lineages defined by the K158N/N189K fixation are highlighted. Scale bar indicates the number of nucleotide substitutions per site. Tree rooted to A/Brisbane/10/2007.

In 2014 two independent substitutions at position 159, F159S and F159Y, fixed in distinct A/H3N2 lineages co-circulating globally (Appendix 1—figure 8). The S159Y substitution was detected sporadically in 2012/2013 before fixing in 2014 to define the A/H3N2 clade 3C.2a, which circulated as the dominant clade globally for the next three years. The F159S substitution was not detected in 2012–2013 before reaching fixation in 2014 to define clade 3C.3a, which continued to circulate globally at low frequencies.

Residue 159 is located in antigenic site B. The S159Y substitution in combination with Y155H and K189R substitutions previously resulted in the transition of the A/H3N2 Bangkok/79-like antigenic phenotype to Sichuan/87-like phenotype (Koel et al., 2013).

The two lineages have independent mutational trajectories for their HA gene away from their most recent common ancestor, excluding both acquiring the N225D substitution prior to F159X-fixation, with acquired changes at the major antigenic epitopes (Appendix 1—table 4, Appendix 1—figure 8; Koel et al., 2013).

Branch tip color indicates the amino acid identity at position 159. Co-circulating lineages defined by the F159X fixation are highlighted. Scale bar indicates the number of nucleotide substitutions per site. Tree is rooted to A/Perth/16/2009.

The two within-host studies both make two key assumptions: (1) immune selection pressure is strong from the moment of inoculation in experienced hosts and (2) influenza virus transmission bottlenecks are wide (1 × 105 virions [Luo et al., 2012]; the three most frequent within-host variants deterministically transmitted onward in Volkov et al., 2010). As discussed in Section A2, antibody-mediated selection pressure is in fact likely to vary substantially over the course of infection. Recent empirical studies have found that transmission bottlenecks are in fact very tight (on the order of a single virion [McCrone et al., 2018; Xue and Bloom, 2019]).

Both within-host studies find that new antigenic variants are routinely generated de novo during an infection and undergo substantial positive selection during that same infection. NGS studies of natural human infections suggest that this is in fact uncommon (Debbink et al., 2017; McCrone et al., 2018; Xue and Bloom, 2020).

The most crucial assumption in prior work on influenza virus immune escape within-hosts (Luo et al., 2012; Volkov et al., 2010) or within-host pathogen immune escape generally (Kennedy and Read, 2017) is that immunity protects against observable reinfection by neutralizing pathogens during the timecourse of replication. This is immunologically unrealistic for a virus that replicates as rapidly as influenza virus; it 'outruns' the memory B-cell response.

What is more, our models show that it protection via neutralization produces binary outcomes: either no reinfection, or detectable reinfection with exclusively mutant viruses. This binary outcome was previously considered a feature, not a bug, because the frequency of homotypic reinfection was not yet known.

A model of immune and therapeutic escape for a generic pathogen (Kennedy and Read, 2017) studied ideas related to replication selection to argue that prophylactic anti-pathogen interventions (such as pre-exposure vaccination) could be less vulnerable to pathogen evolutionary escape than therapeutic interventions (such as post-symptomatic courses of antivirals). Therapeutic interventions occur after a number of pathogen replication events and therefore select upon on a larger, more diverse array of potential pathogen variants. Prophylactic interventions may limit the number of replications before clearance and thereby limit opportunities for diversification. Like the influenza-specific models discussed, this general model does not address the question of how evolution can be rare in symptomatic or otherwise observable infections, where substantial antigenic diversity can be generated.

Similarly, the antimicrobial resistance literature makes a distinction between 'acquired' and 'transmitted' (or 'primary') drug resistance (Bonhoeffer et al., 1997), but this should not be confused with replication and inoculation selection. Transmitted resistance refers to the acquisition of a resistant infection from a host infected with resistant microbes, without regard for whether those microbes were a minority or majority variant in the transmitting host. Inoculation selection, in contrast, refers to natural selection acting on inoculated diversity that favors transmitted drug resistance variants or transmitted new antigenic variants. One variant present in a mixed infection has a higher probability of surviving the transmission bottleneck or becoming the majority variant in the new host than its frequency in the transmitting host alone would imply, simply because the recipient host is more susceptible to that variant than to any competing inoculated variants.

Two studies (Luo et al., 2012; Hartfield and Alizon, 2014) have previously noted that the evolutionary emergence of escape mutants or otherwise fitter pathogens can be made more difficult due to ongoing competition from the old variant for a shared resource, whether susceptible cells within a host, as in Luo et al., 2012 and in our model, or susceptible hosts within a population, as in Hartfield and Alizon, 2014.

In our paradigm, non-variant-specific limiting factors such as target cell depletion play three roles. (1) Denying new variants the opportunity to rise in frequency once an antibody response has been mounted and fitness differences have become substantial. (2) Explaining why immune-compromised hosts, in whom these non-antigenic limiting factors are weaker, can select for antigenic mutants at the scale of a single infection. (3) Clearing infections of naive hosts well before an adaptive immune response can select an escape variant. Point (3) has been posited previously by [60], but to our knowledge (1) and (2) are novel.

Hartfield and Alizon, 2014 focus on the probability that a fit mutant variant the emerges at the host population scale goes stochastically extinct before causing a large epidemic, and show that competition for hosts with the old variant raises this extinction probability. While this may be relevant to influenza viruses at the population level (see Section A8.6 below), an analogous within-host argument is unlikely to be a sufficient explanation for the rarity of influenza virus antigenic mutants within human hosts. New variants are likely to be generated when an influenza virus infection is still well within its exponential growth phase, competition for susceptible cells is relatively unimportant, and probability of stochastic extinction for all incipient lineages is therefore low (due also to influenza’s high ℛ0). The absence of selection, rather than a high rate of stochastic loss of mutants during within-host replication, is the most important component of our explanation for the rarity of observable influenza virus antigenic evolution within individual human hosts.

Note that Hartfield and Alizon, 2015 also have a model of within-host immune escape for a general pathogen, but that model is designed to estimate the contribution of an active adaptive immune response to increasing or reducing the probability that a small escape mutant population goes stochastically extinct.

For influenza viruses, the de novo mutation that creates a stable new variant lineage of interest typically precedes time of adaptive immune proliferation, so a distinct model is required.

Previous population-level studies have argued that non-antigenic fitness costs associated with antigenic substitutions (Gog, 2008; Kucharski and Gog, 2012) or deleterious mutational load on the influenza virus genome (Koelle and Rasmussen, 2015) are necessary to explain the pace of influenza virus antigenic evolution in the human population. These models introduce population-level antigenic novelty at rates 5–10 times higher than predicted by inoculation selection, and at equal rates early and late during the circulation of a particular antigenic variant.

There are empirical reasons to doubt that influenza viruses are in fact antigenic diversity-limited due to fitness costs. We find multiple instances of polyphyletic cluster transitions: a cluster- transition amino acid substitution arises and proliferates on two or more branches of the virus phylogeny at this same time (see Section A6). This is inconsistent with the hypothesis that antigenic variants are either intrinsically unfit (deleterious substitutions) or incidentally unfit (poor background) prior to the moment that they proliferate at the population level (Gog, 2008; Kucharski and Gog, 2012; Koelle and Rasmussen, 2015), since the cluster-transition mutant can reach surveillance-detectable levels even before substantial population immunity has accumulated (suggesting that strong antigenic selection is not required to overcome intrinsic deleteriousness), and it can fully emerge when the time is ripe against multiple independent genetic backgrounds, suggesting background may not suffice to constrain diversification earlier in a cluster’s circulation.

In the absence of meaningful replication selection, the population level rate of antigenic diversification is the average rate at which new variants survive all transmission bottlenecks to found or co-found infections. In an inoculation selection regime, that rate is unlikely to surpass 2 in 104 inoculations (Figure 4A). We obtain that estimate assuming no within-host replication cost for antigenic mutants, weak cross immunity between mutant and old variant viruses, and many experienced hosts who reliably neutralize old variant virions at the point of transmission. Actual rates of new variant survival are almost certainly lower. First, we have assumed that a new antigenic variant can be produced by a single amino acid substitution (accessible via a single nucleotide substitution), however a recent detailed study of the antigenic evolution of A/H3N2 viruses shows that some new variants require more than one substitution (Koel et al., 2013). Second, some (possibly many) previously infected hosts will not possess well-matched antibodies due to original antigenic sin (Davenport and Hennessy, 1956), antigenic seniority (Lessler et al., 2012), immune backboosting (Fonville et al., 2014), or other sources individual-specific variation in antibody production (Lee et al., 2019). These factors individually, and particularly in combination, which have not been modeled in this study mean that the above 2 in 104 transmission estimate is likely to be unrealistically high. By comparison, existing population level models that invoke fitness costs introduce population-level antigenic diversity at much higher rates. One study (Koelle and Rasmussen, 2015) introduces new antigenic variants at a rate of 7.5 per 104 transmissions; another (Kucharski and Gog, 2012) introduces new population-level mutations at a continuous rate of 6.8 × 10−4 mutations per infected individual per day, which should produce new variant infections by 2 days post infection in at least 1 of every 103 infected hosts:

1−exp⁡(−6.8×10−4∗2)≈0.0013

Prior to the accumulation of population immunity, infections dominated by new variants should be rare, and new variants should be nearly neutral relative to the old variant at the population level. This may be sufficient to explain the absence of diversification early in the circulation of a cluster without invoking non-antigenic fitness costs. That said, additional constraints on the virus, particularly those that reduce the within-host fitness of new antigenic variants, should further slow virus antigenic diversification by reducing variant frequencies in donor hosts and thereby reducing the probability that the variants survive transmission bottlenecks. Furthermore, in the absence of substantial population immunity, adaptive substitutions may be lost at the population level due to competition from lineages with non-antigenic adaptive substitutions. Our work is thus consistent with a role for clonal interference in population level influenza virus antigenic evolution (Strelkowa and Lässig, 2012).

In this section, we expand on the discussion of alternative hypotheses from the main text (see Alternative explanations for rare new antigenic variants).

Kennedy and Read, 2017 and Luo et al., 2012 study immune escape with ℛw⁢(0)<1 but consider only diversity generated after inoculation. Our study complicates their results: we find that when ℛw⁢(0)<1 the dominant mode of antigenic evolution will be selection on inoculated diversity, not selection on generated diversity (Figure 2). The pace of evolution will depend on how likely an escape mutant is to be transmitted to an experienced host. If bottlenecks in experienced hosts are wide, evolution may be rapid. If these bottlenecks are tight or escape mutants are deleterious within untreated naive hosts (so fmt is low), then ℛw⁢(0)<1 at the start of infection can result in slow evolution, since both replication and inoculation selection will be rare.

In the absence of mucosal neutralization with small fmt, b≪1fmt, and small µ, the threshold for where inoculation selection becomes more common than replication selection with ℛw⁢(0)<1 is approximately

fmtbpsse>μbℛw(0)1−ℛw(0)psse

or simply:

(A29) fmt>μ⁢ℛw⁢(0)1-ℛw⁢(0)

the left-hand side comes from Equation A17, which gives an probability that a mutant with ℛm⁢(0)>1 is inoculated. The right-hand side comes from Equation 28, which gives the probability that an escape mutant is generated in a declining but replicating virus population before that population goes extinct. psse, which occurs on both sizes and cancels, is the probability that a mutant survives stochastic extinction. We have also ignored the fact the right-hand side should have a term representing the probability that inoculation selection does not occur, but as we are considering cases when both inoculation selection and replication selection are rare, that term is approximately one.

On average, fmt>μ due the asymmetry between forward and back-mutation rates when the mutant is rare. It follows that fmt>μ⁢ℛw⁢(0)1-ℛw⁢(0) should hold when (A) the mutant not is too deleterious in the absence of positive selection (since this reduces fmt) or (B) ℛw⁢(0) is not close to 1 (since then many copies are made on average before the virus goes extinct). Alternatively, if b is sufficiently large—on the order of 1fmt—selection on inoculated diversity will be the dominant mode of evolution simply because most inocula will include mutants.

Selection on inoculated diversity may be important in many host-pathogen systems, not just in influenza virus antigenic evolution. Some anti-microbial resistance in bacteria, for instance, is acquired through the uptake of preexisting plasmids, rather than through de novo mutation (Perron et al., 2015). Whether drug resistance emerges in an individual treated host will depend on whether any of the bacteria inoculated bear the needed plasmid.

Population size, population level antigenic diversification ('mutation') rate, and degree of population structure all substantially impact the proliferation of new influenza virus antigenic variants at the population level. It has been hard for population-level models to distinguish plausible hypotheses regarding evolutionary constraints on the virus, because achieving simultaneous realism across in all these effects while retaining model tractability has proven extremely difficult.

For epidemiological models of influenza virus evolution overall population size matters: epidemics in larger populations involve more total inoculations, and therefore more opportunities for new variant survival at the point of transmission to occur. Large host populations with high transmission rates, strong population connectivities, and recurrent epidemics should promote antigenic evolution. This helps explain why influenza virus antigenic evolution appears to occur disproportionately in east and southeast Asia (Russell et al., 2008).

It also reveals an important consideration for interpreting results from individual-based simulation models of global influenza virus evolution. Due to computational constraints, many such models must use host population sizes that are orders of magnitude smaller than the true global human population (e.g. N= 40 million [Koelle and Rasmussen, 2015], N= 90 million [Bedford et al., 2012]). Such models will have smaller numbers of inoculations than real-world influenza virus dynamics. They therefore run the risk of overestimating rates of strain extinction or, to avoid this, of overestimating per-infection virus diversification rates.

The observation that immune escape can be reliably be selected for in egg (Davis et al., 2018) and mouse (Hensley et al., 2009) passage experiments and yet appears rare in reinfected experienced humans is unsurprising in light of our view that influenza virus evolution is limited more by the failure of antibody-mediated selection pressure and antigenic diversity to coincide than by the absence of either.

The egg experiments (Davis et al., 2018) amount to strong replication selection: the viruses were passaged in the presence of sera with highly-specific neutralizing antibodies, with these sera present at all times, including during the exponential growth phase within each egg.

The mouse passage experiments (Hensley et al., 2009) showed that serial passage in immune mice led to the appearance and fixation of antigenic escape mutants, but serial passage in naive mice did not. Several details of the experiment suggest that both replication selection and inoculation selection should have been more possible than in typical infected experienced humans.

First, the mice were intranasally inoculated with 50 μ⁢l of homogenized lung isolate from the previous mouse in the passage chain at two days post-infection. This is likely a substantially more concentrated dose of virions than is inhaled through aerosol or even contact transmission (an inter-host bottleneck much wider than occurs in humans, to use the terminology of main text Figure 4A). This could produce a very large value of v, the number of virions encountering IgA, and potentially a large value of b, the final (cell infection) bottleneck as well—in the presence of sufficiently many virions, early cell infections might be sufficiently simultaneous as to produce have observably diverse within-host populations. All of this should tend to facilitate new variant survival and promotion to observable levels, but also to improve chances, relative to humans, that at least some old variant virions could survive alongside the new variant.

Second, inoculations were also carried out until a successful inoculation could be achieved. This further improves the chances of observing an escape mutant selected at the point of transmission. Finally, vaccinated mice in the passage chain were inoculated 10–21 days post-vaccination, meaning antibody levels could be higher than the memory baseline. This would be expected to strengthen inoculation selection and perhaps permit replication selection.

One detail in supplementary table 1 of the mouse passage study (Hensley et al., 2009) is particularly relevant. In the passaging of virus stock #3, two mutants were observed for the first time in the second vaccinated mouse in the chain, but neither at fixation. Both remained present for 5 further passages in experienced mice without either being lost or fixing. This suggests a much wider final bottleneck than occurs in naturally-infected humans (McCrone et al., 2018). It also suggests that replication selection is unlikely to have been at all strong if present at all: over repeated exponential growth phases of two days with wide bottlenecks, even very small fitness differences can be easily amplified, so the more fit escape mutants should have fixed if they were subject to replication selection. It also suggests, as expected, that inoculation selection is weak if final bottlenecks are wide and neither antigenic variant is reliably neutralized. In that case, a long intermittent period of coexistence is possible.

An example calculation illustrates this: if v=1000, b=10, fmt=0.5, κm/κw>0.9, and κw<0.8 the mutant has an inoculation-selective advantage, but the expected mutant frequency at after the next transmission event approximately 0.58 or less. Letting w¯=E⁢(xw) and m¯=E⁢(xm):

(A30) w¯>v(1−fmt)(1−κw)=(1000)(0.5)(0.20)=100m¯>v(fmt)(1−κm)=(1000)(0.5)(0.28)=140

And since w¯,m¯≫b, the hypergeometric sampling is approximately binomial, and therefore the expected frequency of mutant is just m¯m¯+w¯≈0.58

While it is difficult to compare fitness differences during replication to fitness differences in mucosal neutralization, a very small fitness difference of δ=0.5 should raise the fitter type from frequency 0.5 to frequency 0.73 over the course of two days.

When antigenic effect size of mutations approaches zero, the effects of stochastic loss and mistimed selection pressure dominate, and inoculation and replication selection become sufficiently weak that true drift dominates as a force for introducing new antigenic variants. In such a regime, we expect influenza viruses to evolve gradually, fail to cause large epidemics, and possibly diversify, not unlike influenza B viruses (Bedford et al., 2015). In particular, slow spread enables a kind of immunological niche partitioning: if population immune histories are not spatially uniform because epidemiological spread is slow relative to antigenic diversification, multiple variants that are each locally favored in distinct areas could co-circulate and form separate lineages, as has happened with B/Victoria and B/Yamagata.

When antigenic jump size approaches the maximum size possible, such that all population immunity disappears each time there is an antigenic cluster transition (Smith et al., 2004), we expect influenza viruses to follow a strongly clock-like pattern with high-amplitude oscillations in case numbers. They would cause massive punctuated epidemics during jump years and then in the next year either select for a new variant or go extinct. There would be huge booms followed by one or more years of bust. In between, at intermediate degrees of escape, a punctuated but somewhat less clocklike pattern—such as the pattern observed in nature for A/H3N2 viruses (Smith et al., 2004)—becomes possible.

'Preview' substitutions and polyphyletic cluster transitions also imply that the influenza virus is not typically mutation-limited at the population level—a substantial-effect escape mutant is usually accessible in sequence-space—but that the virus may nonetheless be highly constrained in its evolutionary trajectory. The virus repeatedly finds the same substitution as a solution to its antigenic evolutionary problem. This could occur either because only one escape mutant is available or because one of the available escape mutants provides substantially more immune escape (on average) than the others.

Without an explanation for why within-host dynamics do not more frequently promote antigenic mutants, it is difficult to explain the slow pace and noisy trajectory of influenza virus evolution. But such a within-host explanation, though likely necessary, may not be sufficient. It is conceivable that in a sufficiently large and well-connected global host population, population-level antigenic selection favoring mutants with higher ℜe would lead even small-effect antigenic mutants rapidly to emerge and fix.

One important mechanism by which evolution might be also slowed at higher scales is analogous to the within-host dynamic of founder effects and inoculation selection: mutant lineages may be frequently lost (though perhaps selectively favored) at the bottleneck that occurs between distinct influenza virus epidemics.

Prior modeling work has found that population-level host competition between pathogen variants reduces the probability that a fit mutant variant causes a large epidemic in the population in which it first emerges (Hartfield and Alizon, 2014). This is likely to be relevant in influenza virus epidemics, since the virus is thought to provoke short-term strain-transcending (i.e. variant-transcending) immunity (Ferguson et al., 2003).

We propose a previously unexplored consequence of this argument: successful establishment of a mutant lineage at the population level requires that the mutant lineage be exported to a new host subpopulation where susceptible hosts are common. This is rare but potentially selective sampling event.

When a new variant lineage emerges at the population level, it is likely to be surrounded by many propagating old variant lineages due to the generator-selector dynamic discussed in the main text (see main text Figure 4B,C). Moreover, most mutant lineage generation events will occur as an epidemic is peaking, since that is when the most inoculations occur. The consequence is that most generated population-level mutant lineages will encounter severe competition for susceptible hosts (especially if we consider realistic, spatial host contact networks rather than well-mixed epidemic models). So a generated mutant lineage is unlikely to account for many cases—or even necessarily more than one case—in the epidemic in which it emerges.

Many local influenza virus epidemics are likely to be evolutionary dead ends with no cases exported that establish chains of transmission in other locations. Because only a minority of the human population is likely to travel to or from the site of any particular epidemic, the majority of virus diversity generated within each epidemic is likely to be lost in between-epidemic bottlenecks. Conditional on being exported, new variant lineages could have competitive advantages over old variant lineages because they spread should spread more rapidly and go stochastically extinct less easily due to their higher ℜe. It follows that there is analogous dynamic to inoculation selection at the population level: mutant lineages are rarely exported from one sub-population (host, epidemic) to another, but conditional on being exported, they have an advantage in any potential competition with exported old variant lineages.

The rate of proliferation of mutants at the population level, then, may be limited by the rarity of early generation of mutant lineages during an epidemic (analogous to the rarity of very early within-host de novo generation of mutants) and by the rarity of successful mutant lineage exportation when generation is not early (analogous to the rarity of mutant virions surviving the transmission bottleneck between hosts).

We aim to explore this argument with a formal mathematical model in future work.

The derivation in this section establishes Equation 12:

dfmdt=fm(1−fm)([gm(t)−gw(t)]−[dm(t)−dw(t)])

With symmetric mutation at a rate µ, the replicator equation remains calculable.

This establishes the expression for the time t* by which a mutant must emerge to reach at least frequency x by time t given in Equation 25 of the Materials and methods.

t∗(x,t)={t−∗(x,t)t−∗(x,t)<t and t−∗(x,t)≤tMt+∗(x,t)t+∗(x,t)<t and t−∗(x,t)>tMtotherwise

where:

t-*⁢(x,t)=ln⁡[1-xx⁢exp⁡(δ⁢(t-tM))+1]-ln⁡bg0-dv

t+*⁢(x,t)≈ln⁡(1-xx)-ln⁡(b)+δ⁢t-cw⁢k⁢tMg0-dv-cw⁢k+δ

The derivation in this section establishes Equation 35:

pcib⁢(κw,fmt,w⁢b)=E⁢(pcib⁢(xw,1,b))=(1-e-w¯)⁢bw¯+e-w¯⁢∑j=0b-1w¯jj!⁢(1-bj+1)

where w¯=v⁢(1-fmt)⁢(1-κw) is the mean number of old variant virions present after IgA neutralization.

pcib<1 if κw<1 and fmt<1, and pcib=1 if κw=1 or fmt=1.

pcib is decreasing in w¯: the more competing old variant virions, the less likely the new variant is to pass through the cell infection bottleneck.