Our model reflects the following biological realities: (1) Seasonal influenza virus infections of otherwise healthy individuals typically last 5–7 days (Suess et al., 2012); (2) In influenza virus-naive individuals, it can take up to 7 days for anti-influenza virus antibodies to start being produced (Wrammert et al., 2008), effectively resulting in no selection (Figure 1A); (3) In previously infected (‘experienced’) individuals, sIgA antibodies can neutralize inoculated virions before they can infect host cells; (Wang et al., 2017) (4) However, if an inoculated virion manages to cause an infection in an experienced individual, it takes 2–5 days for the infected host to mount a recall adaptive immune response, including producing new antibodies (Coro et al., 2006; Zuccarino-Catania et al., 2014) (see Appendix Section A2 for further discussion of motivating immunology). Importantly, this contrasts with previous within-host models of virus evolution, which have assumed that antibody-mediated neutralization of virions during virus replication is strong from the point of inoculation onward and is the mechanism of protection against reinfection (Luo et al., 2012; Volkov et al., 2010). It also reflects new animal model evidence of sterilizing antibody immunity (Le Sage et al., 2020). We discuss existing models and hypotheses for the rarity of population-level influenza antigenic variation in Appendix Section A7.

Our model can be parameterized to reflect different hypothesized immune mechanisms, different host immune statuses, and different durations of infection. In the model, virions $V_i$ of antigenic type i infect target cells ${\displaystyle C}$, replicate, mutate to a new antigenic type j at a rate ${\mu }_{ij}$, and decay at a rate d_v. We model the innate immune response implicitly as depletion of infectible cells. We model the antibody-mediated immune response as an increase k in the virion decay rate in the presence of well-matched antibodies. To model partial antibody cross reactivity, we scale k by a parameter $c_i\in [0,1]$; c_i reflects the binding strength of the host’s best-matched antibodies to antigenic type i. So in the presence of an antibody response, virions $V_i$ of type i decay at a rate $d_v+c_{i}k$.

The model can accommodate $N_v$ antigenic variants $i=1,2,\mathrm{\dots }{N}_v$ linked by an arbitrary network of possible substitutions and corresponding mutation probabilities ${\mu }_{ij}$, but in practice we typically consider two, the new variant m and the old variant w, and neglect back-mutation from new variant to old variant (${\mu }_{wm}>0$, but ${\mu }_{mw}=0$).

To assess the importance of transmission bottlenecks, initial virus diversity, and sIgA antibody neutralization in virus evolution, we model the point of transmission as a series of stochastic events which may ultimately lead to one of more virions invading cells and initiating an infection. The recipient host is inoculated with a random sample of within-host virus diversity from the transmitting host. In experienced hosts, this inoculum is probabilistically thinned by host antibodies. The founding population that initiates the infection is then randomly sampled from among any remaining virions.

Mathematically, we model the number of inoculated virions as Poisson-distributed with a mean v, so if variant i has frequency f_i within the transmitting host, the number of variant i virions inoculated is Poisson-distributed with mean $v{f}_i$. The virions then encounter antibodies, which we interpret as sIgA but can be understood to be any antigen-specific antibody-mediated protection that precedes cell infection; each virion of variant i is independently neutralized with a probability ${\kappa }_i$. This probability depends upon the strength of protection against homotypic reinfection $\kappa$ and the sIgA cross immunity between variants ${\sigma }_{ij}$, ($0\le {\sigma }_{ij}\le 1$). So if a host with antibodies to variant j is challenged with variant i, those virions will be neutralized with a probability ${\kappa }_i=\kappa {\sigma }_{ij}$. For simplicity, we assume the same homotypic protection level across all variants and hosts, though in practice there may be variation in the immunogenicity of individual variants and in the strength of responses generated by individual hosts. We typically fix host immune histories to test the effect of host immune history on selective dynamics. When necessary, we can model a novel (non-recall) antibody response to a strain i by designating the host as experienced to i at some time ${\displaystyle t_N^i}$ post-infection (see Materials and methods).

The model is continuous-time and stochastic: cell infection, virion production, virus mutation, and virion decay are stochastic events that occur at rates governed by the current state of the system, with exponentially distributed (memoryless) waiting times. The system is approximately well-mixed: we track counts of virions and cells without an explicit account of space within the upper respiratory tract. We treat infected and dead or removed cells implicitly.

Parameterized as in Table 1, the model captures key features of influenza infections: rapid peaks approximately 2 days post-infection and slower declines with clearance approximately a week post-infection, with faster clearance in experienced hosts.

We give a full mathematical description of the model in the Materials and methods.

In addition to analyzing this within-host model, we explored the between-host and population-level implications of these within-host dynamics using a simple transmission chain model and an SIR-like population-level model, which we describe in the Materials and methods.

In our model, sufficiently strong antibody neutralization during the virus exponential growth period can potentially stop replication and block onward transmission, but this mechanism of protection results in detectable new antigenic variants in each observable homotypic reinfection, since the infection is terminated rapidly unless it generates a new antigenic variant that substantially escapes antibody neutralization (Figure 2).

Figure 2

Download asset Open asset

If there is antibody neutralization throughout virus exponential growth and it is not sufficiently strong to control the infection, this facilitates the establishment of new antigenic variants: variants can be generated de novo and then selected to detectable and easily transmissible frequencies (Figure 1C,D,F, sensitivity analysis in Appendix 1—figure 3). We term selection on a replicating within-host virus population ‘replication selection’. Virus phenotypes that directly affect fitness independent of immune system interactions are likely to be subject to replication selection.

Adding a realistic delay to antibody production of two days post-infection (Lam and Baumgarth, 2019) curtails antigenic replication selection. There is no focused antibody-mediated response during the virus exponential growth phase, and so the infection is dominated by old antigenic variant virus (Figure 1C,D,G, sensitivity analysis in Appendix 1—figure 3). Antigenic variant viruses begin to be selected to higher frequencies late in infection, once a memory response produces high concentrations of cross-reactive antibodies. But by the time this happens in typical infections, both new antigenic variant and old antigenic variant populations have peaked and begun to decline due to innate immunity and depletion of infectible cells, so new antigenic variants remain too rare to be detectable with NGS (Figure 1G).

We find that replication selection of antigenic novelty to detectable levels becomes likely only if infections are prolonged, and virus antigenic diversity and antibody selection pressure therefore coincide (Figure 1G, see also Figure 7). This can explain existing observations: within-host adaptive antigenic evolution can be seen in prolonged infections of immune-compromised hosts (Xue et al., 2017), and prolonged influenza infections show large within-host effective population sizes (Lumby et al., 2020).

Adding antibody neutralization of virions at the point of inoculation (e.g. by mucosal sIgA) to our model produces realistic levels of protection against reinfection, and when reinfections do occur, they are overwhelmingly old antigenic variant reinfections. New antigenic variants that arise during these reinfections remain undetectably rare, reproducing observations from natural human infections (Figures 1A–D,G–H, Figure 3B–C; Debbink et al., 2017; Dinis et al., 2016; McCrone et al., 2018; Sobel Leonard et al., 2016). The combination of mucosal sIgA protection and a realistically-timed antibody recall response explains how there can exist immune protection against reinfection—and thus a population-level selective advantage for new antigenic variants—without observable within-host antigenic selection in typical infections of experienced hosts.

Figure 3

Download asset Open asset

Regardless of host immune status, an antigenic variant that has been generated de novo within a host must survive a series of population bottlenecks if it is to infect other individuals. To found a new infection, virions must be expelled from a currently infected host (excretion bottleneck), must enter another host (inter-host bottleneck), must escape mucus on the surface of the airway epithelium (mucus bottleneck), must avoid neutralization by sIgA antibodies on mucosal surfaces (sIgA bottleneck), and must infect a cell early enough to form a detectable fraction of the resultant infection (cell infection bottleneck) (Figure 3A). The sum of all of these effects is the net bottleneck and typically results in infections being initiated by a single genomic variant (McCrone et al., 2018; Xue and Bloom, 2019; Ghafari et al., 2020). That said, bottlenecks resulting from direct contact transmission may be substantially wider than those associated with respiratory droplet or aerosol transmission (Varble et al., 2014) and more human studies are required to quantify these differences.

We find that because antigenic variants appear at very low within-host frequencies when generated de novo and undergo minimal or no replication selection, new antigenic variants most commonly reach detectable levels within hosts through founder effects at the point of inter-host transmission: a low-frequency antigenic variant generated in one host survives the net bottleneck to found the infection of a second host (Figure 3D).

Given that influenza bottlenecks are thought to be on the order of a single virion (McCrone et al., 2018), any replication-competent mutant that founds an infection should occur at NGS-detectable levels, and likely at consensus. But for the same reason, these founder effects are rare events. The sampling process that produces these founder effects could be a purely neutral. It is likely quite close to neutral in a truly naive recipient host who does not possess well-matched antibodies to the inoculated old variant: all inoculated virions, regardless of antigenic phenotype, have an equal chance of becoming part of the new infection’s founding population. If there are v virions that compete to found the infection and $b=1$, then each virion founds the infection with probability $1/v$.

In our model, new antigenic variants therefore survive the transmission bottleneck upon inoculation into a naive host with a probability approximately equal to the donor-host variant frequency $f_{\text{mt}}$ times the bottleneck size b (see Materials and methods). New antigenic variant infections of naive hosts should then occur on the order of 1 in 10⁵ or 1 in 10⁴ such infections given biologically plausible parameters (Figure 3E–H).

But if different virions have different chances of being neutralized at the point of transmission, the founding process may be selective. Among the virions that encounter antibodies, those that are less likely to be neutralized have a higher than average chance of undergoing stochastic promotion to consensus while those that are more likely to be neutralized have a lower than average chance. A new antigenic variant may then be disproportionately likely to survive the net bottleneck (Figure 3, Appendix Section A4.2). We term this potential selection on inoculated diversity ‘inoculation selection’. Neutralization at the point of transmission thus not only gives new antigenic variant infections their transmission advantage (population-level selection) but may also increase the rate at which these new antigenic variant infections arise (inoculation selection).

There is some suggestive evidence of differential survival of particular (not necessarily antigenic) influenza genetic variants at the point of transmission from experiments in ferrets (Wilker et al., 2013; Moncla et al., 2016). But as Lumby and colleagues (Lumby et al., 2018) point out in a reanalysis of those experiments, it is difficult empirically to distinguish selection that occurs at the point of transmission from selection that occurs during early replication in the recipient host because of the challenges associated with sampling the small virus populations present at the earliest stages of infection. Here we define inoculation selection as selection on the bottlenecked virus population that establishes infection in the recipient host before any virus replication has taken place in that host.

The strength of inoculation selection depends on the ratio of the number of virions that compete to found an infection in the absence of well-matched sIgA antibodies (the mucus bottleneck size v) to the number of virions that actually found an infection (the final cell infection bottleneck size b). The larger this $v/b$ ratio is, the more inoculation selection in experienced hosts facilitates the survival of new antigenic variants (Figure 3F, Appendix 1—figure 1).

When new antigenic variant immune escape is incomplete due to partial cross-reactivity with previous antigenic variants, increased antibody neutralization is a double-edged sword for new antigenic variant virions. Competition to found the infection from old antigenic variant virions is reduced, but the new antigenic variant is itself at greater risk of being neutralized. The impact of inoculation selection therefore depends on the degree of similarity between previously encountered viruses and the new antigenic variant. An experienced recipient host could facilitate the survival of large-effect antigenic variants (like those seen at antigenic cluster transitions [Smith et al., 2004]) while impeding the survival of variants that provide less substantial immune escape (Figure 3E, Appendix 1—figure 1).

Inoculation selection is limited by the low frequency $f_{\text{mt}}$ of new antigenic variants in transmitting donor hosts (due to weak replication selection), the potentially small mucus bottleneck size v, and the fact that some hosts previously infected with the old variant or similar antigenic variants might 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 of individual-specific variation in antibody production (Lee et al., 2019). These factors combined make selection and onward transmission of new variants rare.

Onward transmission of new variants can be facilitated by natural selection—replication selection, inoculation selection, or both. The degree of facilitation depends principally on four quantities: (1) $\delta \tau$, the product of the replication selection fitness difference $\delta =k(c_w-c_m)$ and the time under replication selection $\tau$. This determines the degree to which the new variant is promoted by replication selection prior to transmission (increasing $f_{\text{mt}}$). (2) ${\kappa }_w$, the sIgA neutralization probability for the old variant. This must be large enough to reduce competition for the final bottleneck. (3) $v/b$, the ratio of the number of virions that encounter sIgA v to the cell infection (final) bottleneck size b. This determines the degree of competition to found the infection, and thus sets the maximum potential strength of inoculation selection when ${\kappa }_w$ is large: a $\frac{v}{b}$-fold improvement over drift for small b. (4) $1-{\kappa }_m$, how likely the new variant is to avoid neutralization at the sIgA bottleneck; this scales down the maximum inoculation selective strength set by $v/b$. Inoculation selection can impede new variant survival relative to drift if ${\displaystyle 1-{\kappa }_m}$ is small enough.

When $\delta \tau$ is small and ${\kappa }_w$, $v/b$, and $1-{\kappa }_m$ are large, new variant survival is most facilitated by inoculation selection. When the opposite is true, replication selection is most important. And there are parameter regimes in which both replication and inoculation selection provide a substantial improvement over drift (Figure 4, see Appendix Section A4.6 for mathematical intuition for these results).

Figure 4

Download asset Open asset

At realistic parameter values and assuming all individuals develop well-matched antibodies to previously encountered antigenic variants, only ∼1 to 2 in 10⁴ inoculations of an experienced host results in a new antigenic variant surviving the bottleneck (Figure 4E,H, Figure 4). This rate is likely to be an overestimate due to the factors mentioned above. Moreover, it is only about 2- to 3-fold higher than the rate of bottleneck survival in naive hosts, where new antigenic variant infections should occasionally occur via neutral stochastic founder effects. In short, even in the presence of experienced hosts, antigenic selection is inefficient and most generated antigenic diversity is lost at the point of transmission. Because of these inefficiencies, new antigenic variants can be generated in every infected host without producing explosive antigenic diversification at the population level.

To investigate the between-host consequences of adaptation given weak replication selection, tight bottlenecks, and possible inoculation selection, we simulated transmission chains according to our within-host model, allowing the virus to evolve in a 1-dimensional antigenic space (Smith et al., 2004; Bedford et al., 2012) until a generated antigenic mutant became the majority within-host variant. When all hosts in a model transmission chain are naive, antigenic evolution is non-directional and recapitulates the distribution of within-host mutations (Figure 5A). When antigenic selection is constant throughout infection and even a moderate fraction of hosts are experienced, antigenic evolution is unrealistically adaptive: the virus evolves directly away from existing immunity and large-effect antigenic changes are observed frequently (Figure 5B,C). When the model incorporates both mucosal antibodies and realistically-timed recall responses, major antigenic variants appear only rarely and the overall distribution of emerged variants better mimics empirical observations (Figure 5D)—most notably, the phenomenon of quasi-neutral diversification within an antigenic cluster seen in Figures 1 and 2 of Smith et al., 2004. A simple analytical model (see Methods) in which generated antigenic mutants fix according to their replication and inoculation-selective advantages also displays this behavior (Figure 5B–H).

Figure 5

Download asset Open asset

In particular, we note that whereas an immediate recall response would predict strong near-constant directed evolution of virus antigenic phenotypes away from existing immunity (Figure 5B,C), a realistically-timed recall response predicts that small-effect, drift-like antigenic substitutions will be observed. Even substitutions that move a virus ‘backward’ in antigenic space—so that it is more readily neutralized by existing antibodies than the ancestral variant—can be observed thanks to the large role of stochasticity at the point of transmission. That said, there is a slight bias favoring forward substitutions, especially those of sufficiently large effect to create a substantial selective advantage over the ancestral variant (Figure 5D). Coupled with the plausible assumption that large-effect substitutions are rarer than small-effect substitutions (here captured qualitatively by the Gaussian mutation kernel), this predicts the observed pattern of quasi-neutral diversification within antigenic clusters followed by rarer directional ‘jumps’ in phenotype. Exact rates of antigenic evolution will depend upon how these emergence processes intersect with population-level epidemic dynamics and competition among variants.

We used an epidemic-level model to study the consequences of individual-level inoculation selection for population-level antigenic selection. If inoculation selection is efficient, an intermediate initial fraction of immune hosts maximizes the probability that a new antigenic variant infection is created during an epidemic (Figure 6). This is due to a trade-off between the frequency of naive or weakly immune ‘generator’ hosts who can propagate the epidemic and produce new antigenic variants through de novo mutation, and the frequency of strongly immune ‘selector’ hosts who, if inoculated, are unlikely to be infected, but can facilitate the survival of these new antigenic variants at the point of transmission. As selector host frequency increases, epidemics become rarer and smaller, eventually decreasing opportunities for evolution, but moderate numbers of efficient selectors can substantially increase the rate at which new antigenic variants reach within-host consensus (Figure 6B,C).

Figure 6

Download asset Open asset

We consider several possible explanations for the observed phenomenon that new antigenic variants are rare within experienced infected hosts and at the population level prior to cluster transitions. But among these candidate hypotheses, only the mechanism of small inocula, transmission-blocking mucosal antibodies, and a slow-to-mount recall adaptive immune response can explain all the aforementioned empirical observations simultaneously.

Alternative possibilities include strong immune protection through antibody neutralization during early viral replication (an immediate recall response), heterogeneous neutralization rates during early viral replication, new antigenic variants that are deleterious in the absence of antibody selection, and the need for new variants to emerge against a favorable genetic background.

As shown above, protection through antibody neutralization early in replication can result in rare within-host observation of new antigenic variants, but it contradicts understanding of antibody kinetics and makes other empirical predictions that are unrealistic. It implies that homotypic challenge has a binary outcome: either it results in an undetectable infection that is rapidly cleared or it results in a visible infection dominated by an escape mutant (Figure 2). This is how influenza viruses have been hypothesized to behave in previous models of immune escape (Luo et al., 2012). But empirical work (Debbink et al., 2017; McCrone et al., 2018; Javaid et al., 2020) and human challenge studies (Clements et al., 1986; Memoli et al., 2020) have shown that detectable reinfection of experienced hosts can occur without observable immune escape. Another empirical prediction of such a model is that intermediately immune hosts should be efficient selectors, since they neutralize the old variant virus poorly enough to allow it to grow, but strongly enough to impose antigenic selection upon that growing population (Volkov et al., 2010). While such intermediately immune hosts should be present from the beginning of a new antigenic cluster’s circulation (Fonville et al., 2014), new variants are rarely observed until a new variant has circulated for multiple years (Smith et al., 2004).

Antibody neutralization during early replication could avoid binary outcomes if individual hosts are heterogeneous in the strength of their neutralizing response, so some individuals clear the infection rapidly while others barely exert antibody selection upon it. But while heterogeneity in immunity exists (Lee et al., 2019), this explanation requires extreme, bimodal hetegoneity to avoid the intermediately immunity regime in which replication selection is efficient (Volkov et al., 2010) and again requires an unrealistically early antibody response.

New antigenic variants could be replication-competent, but weakly deleterious within-host in the absence of immune selection and/or compensatory substitutions. Two studies (Gog, 2008; Kucharski and Gog, 2012) have invoked this hypothesis to explain population-level antigenic dynamics. However, a realistically strong antibody response during virus replication could still promote new variants during infections of experienced hosts, even in the absence of compensatory mutations (see Appendix Section A7.4 for a model analysis). So a hypothesis to explain the relative weakness of selection during replication is still required, especially for weakly deleterious mutants that offer substantial immune escape. That said, under our hypothesis of founder effects and possible inoculation selection, weak new variant deleteriousness could further limit the rate of antigenic evolution by reducing the probability that new variants are inoculated into hosts (see Appendix Section A7.4).

Finally, it has been hypothesized that antigenic mutants can only proliferate at the population level if they arise against a favorable genetic background (Koelle and Rasmussen, 2015)—in the absence of deleterious substitutions elsewhere in the genome. But antigenic cluster transitions are frequently polyphyletic (see Appendix Section A6): a new variant emerges quasi-simultaneously in multiple virus lineages. Since these lineages should have different genetic backgrounds, this suggests that favorable backgrounds are readily available, and that emergence is limited instead by the presence or absence of selection pressure.

We discuss these alternative explanations in more depth in the Appendix (Section A7.4).

Previous literature on influenza virus transmission mentions a ‘selective bottleneck’ (Wilker et al., 2013; Moncla et al., 2016; Sobel Leonard et al., 2016), but those studies do not refer to antigenic inoculation selection. Rather, a ‘selective bottleneck’ typically refers either to a tight neutral bottleneck that leads to stochastic loss of diversity or to non-antigenic factors that lead to preferential transmission of certain variants (Moncla et al., 2016). An important exception is (Lumby et al., 2018). Those authors studied ferret transmission experiments and partitioned selection for adaptive mutants (not necessarily antigenic) into selection for transmissibility (acting at a potentially tight bottleneck) and selection during exponential growth. Subsequently, several studies have hypothesized that influenza virus antigenic selection might be weak in short-lived infections of individual experienced hosts and might occur at the point of transmission (Petrova and Russell, 2018; Han et al., 2019; Lumby et al., 2020).

To our knowledge, however, ours is the first study to undertake a mechanistic, model-based comparison between the role of antigenic selection at the point of transmission and that of antigenic selection during replication, to show that immunologically plausible mechanisms could make the former more salient than the latter, and to connect that finding to the rarity of observable new antigenic variants in homotypically reinfected human hosts. We discuss the relationship between this paradigm and those put forward in previous modeling studies at further length in Appendix Section A7.

The study presented here nonetheless has important limitations that suggest opportunities for future investigation. This is a modeling study, and a mainly theoretical one. That is out of necessity. Quality experiments of the kind that are necessary to observe selection at the point of transmission directly and measure its strength have not been published, and we were unable to find any experimental measurements of within-host competition between known antigenic variants. For additional discussion of unmodeled biological realities and mechanistic uncertainties, see Appendix Sections A2 and A5.

Our study has a number of implications for the study and control of influenza viruses.

The asynchrony between within-host virus diversity and antigenic selection pressure provides a simple mechanistic explanation for the phenomenon of weak within-host selection but strong population-level selection in seasonal influenza virus antigenic evolution. Measuring or even observing antibody selection in natural influenza virus infections is likely to be difficult because it is inefficient and consequently rare. Theoretical studies are therefore essential for understanding these phenomena and for determining which measurable quantities will facilitate influenza virus control. Our study highlights a critical need for new insights into sIgA neutralization and IgA responses to natural influenza virus infection and vaccination. Cross-scale dynamics can decouple selection and diversity, introducing randomness into otherwise strongly adaptive evolution.

In all model descriptions, $X+=y$ and $X-=y$ denote incrementing and decrementing the state variable ${\displaystyle X}$ by the quantity ${\displaystyle y}$, respectively, and $X=y$ denotes setting the variable ${\displaystyle X}$ to the value ${\displaystyle y}$. $\dot{X}$ denotes the rate of event ${\displaystyle X}$.

The within-host model is a target cell-limited model of within-host influenza virus infection with three classes of state variables:

• ${\displaystyle C}$: Target cells available for the virus to infect. Shared across all virus variants.

• $V_i$: Virions of virus antigenic variant i

• $E_i$: Binary variable indicating whether the host has previously experienced infection with antigenic variant i ($E_i=1$ for experienced individuals; $E_i=0$ for naive individuals).

New virions are produced through infections of target cells by existing virions, at a rate $\beta C{V}_i$. Infection eventually renders a cell unproductive, so target cells decline at a rate ${\displaystyle \ell \beta C\overline{V}}$, where $\overline{V}$ is the total number of virions of all variants. The model allows mutation: a virion of antigenic variant i has some probability ${\mu }_{ij}$ of producing a virion of antigenic variant j when it reproduces.

Virions have a natural per-capita decay rate d_v. A fully active specific antibody-mediated immune response to variant i increases the virion per-capita decay rate for variant i by a factor k (assumed equal for all variants). The degree of activation of the antibody response during an infection is given by a function $M(t)$, where t is the time since inoculation.

We use a parameter $c_{ij}$ to denote the protective strength of antibodies raised against a strain j against a different strain i. $c_{ii}=1$ by definition and $c_{ij}=0$ indicates complete absence of cross-protection. So if host has antibodies against a strain j but not against the infecting strain i, a fully active antibody response raises the virion decay rate for strain i by $c_{ij}k$. If there are multiple candidate forms of cross protection $c_{ij}$ and $c_{ik}$, we choose the strongest. We typically assume that $c_{ij}=c_{ji}$.

We assume that an antibody immune response is raised whenever the host has experienced a prior infection with a partially cross-reactive strain. For notational ease, we define the host’s strongest cross-reactivity against strain i, c_i, by:

So a recall antibody response is raised during an infection with strains $i,j,\mathrm{\dots }$ whenever one of $c_i,c_j,\mathrm{\dots }>0$.

For this study, we consider a two-antigenic variant model with an ancestral ‘old antigenic variant’ virions $V_w$ and novel ‘new antigenic variant’ virions $V_m$, though the model generalizes to more than two variants. The state variables are charged by the following stochastic events:

These events occur at the following rates:

where $M(t)$ is a minimal model of a time-varying antibody response given by:

For simplicity, the equations are symmetric between old antigenic variant and new antigenic variant viruses, except that we neglect back mutation, which is expected to be rare during a single infection, particularly before mutants achieve large populations. The parameters r_w and r_m allow the two variants optionally to have distinct within-host replication fitnesses; for all results shown, we assumed no replication fitness difference ($r_w=r_m=r$) unless otherwise stated. A de novo antibody response raised to a not-previously-encountered variant i can be modeled by setting $E_i=1$ at a time $t_N^i\ge t_M$ post-infection. By default, we model such a de novo response only for fully naive hosts, and assume that it is mounted only against the variant that was most common at the start of infection, which is typically w.

We characterize a virus variant i by its within-host basic reproduction number ${ℛ}_0^i$, the mean number of progeny virions produced by a single virion at the start of infection in a naive host:

When parametrizing our model, we fixed the within-host basic reproduction number ${ℛ}_0^i$, the initial target cell population $C_{\max }$, the virus reproduction rate r_i, and the shared virion decay rate d_v. We then calculated the implied $\beta$ according to Equation (5).

Another useful quantity is the within-host effective reproduction number ${ℛ}^i(t)$ of variant i at time t: the mean number of progeny virions produced by a single virion of variant i at a given time t post-infection.

Note that ${ℛ}_0^i$ is ${ℛ}^i(0)$ in a naive host, and that if ${ℛ}^i<1$, the variant i virus population will usually decline.

We denote the frequency of new variant virions at time t by ${\displaystyle f_m(t)=\frac{V_m(t)}{V_w(t)+V_m(t)}}$.

The distribution of virions that encounter sIgA antibody neutralization depends on the mean mucosal bottleneck size v (i.e. the mean number of virions that would pass through the respiratory tract mucosa in the absence of antibodies) and on the frequency of new antigenic variant $f_{\text{mt}}=f_m(t_{\text{inoc}})$ in the donor host at the time of inoculation $t_{\text{inoc}}$. n_w old antigenic variant virions and n_m new antigenic variant virions encounter sIgA. The total number of virions $n_{\text{tot}}=n_w+n_m$ is Poisson-distributed with mean v and each virion is independently a new antigenic variant with probability $f_{\text{mt}}$ and otherwise an old antigenic variant. The principle of Poisson thinning then implies:

Note that since $f_{\text{mt}}$ is typically small, the results should also hold for a binomial model of n_w and n_m with a fixed total number of virions encountering sIgA antibodies: $v_{\text{tot}}=v$.

We then model the sIgA bottleneck—neutralization of virions by mucosal sIgA antibodies. Each virion of variant i is independently neutralized with a probability ${\kappa }_i$. This probability depends upon the strength of protection against homotypic reinfection $\kappa$ and the sIgA cross immunity between variants $\sigma$ ($0\le \sigma \le 1$):

Since each virion of strain i in the inoculum is independently neutralized with probability ${\kappa }_i$, then given n_w and n_m, the populations that compete the pass through the cell infection bottleneck x_w and x_m are binomially distributed:

By Poisson thinning, this is equivalent to:

At this point, the remaining virions are sampled without replacement to determine what passes through the cell infection bottleneck, b, with all virions passing through if $x_w+x_m\le b$, so the final founding population is hypergeometrically distributed given x_w and x_m.

If ${\kappa }_w={\kappa }_m=0$, $f_{\text{mt}}$ is small, and v is large, this is approximated by a binomially distributed founding population of size b, in which each virion is independently a new antigenic variant with (low) probability $f_{\text{mt}}$ and is otherwise an old antigenic variant. Alternatively, it can be approximated by a Poisson distribution with a small mean: $f_{\text{mt}}b\ll 1$. So the probability that a new variant survives the bottleneck in the absence of mucosal neutralization is:

When there is mucosal antibody neutralization, the variant’s survival probability can be reduced below this (inoculation pruning) or promoted above it (inoculation promotion), depending upon parameters. There can be inoculation pruning even when the new antigenic variant is more fit (neutralized with lower probability, positive inoculation selection) than the old antigenic variant (see Figure 3E).

Default parameter values for the minimal model and sources for them are given in Table 1.

In this section, we derive analytical expressions for the within-host frequency of the new variant over time in an infected host (Figure 7), the probability distribution for the time of the first de novo mutation to produce a surviving new variant lineage, and the approximate probability of replication selection to frequency x by time t given our parameters.

Figure 7

Download asset Open asset

Distribution of first mutation times

Request a detailed protocol

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, t_e, depends on the mutation rate µ and the per-capita rate at which old antigenic variant virions are produced, $g_w(t)=r_w\beta C(t)$. It also depends on $p_{\text{sse}}$, the probability that the generated new antigenic variant survives stochastic extinction. Denoting the new variant per-capita virion production rate $g_m(t)=r_m\beta C(t)$, we calculate $p_{\text{sse}}$ using a branching process approximation (Ball et al., 2016).

(17) $p_{\text{sse}}=\frac{g_m(t)}{g_m(t)+d_v+kM(t)\max \{E_m,c{E}_w\}}$

Surviving mutants therefore occur at a rate ${\lambda }_m(t)$:

(18) ${\lambda }_m(t)=\mu g_w(t)V_w(t)p_{\text{sse}}$

We define the cumulative rate ${\mathrm{\Lambda }}_m(x)$:

(19) ${\mathrm{\Lambda }}_m(x)={\int }_0^x{\lambda }_m(x)𝑑x$

It follows that the CDF of the first mutation time ${\displaystyle t_e}$ is:

(20) $P(t_e<x)=1-e^{-{\mathrm{\Lambda }}_m(x)}$

This expression is exact for any given realization of the stochastic model if the realized values of the random variables $V_w(t)$, $C(t)$, and $g_w(t)$ are used. In practice, we mainly use it to get a closed form for the CDF of ${\displaystyle t_e}$ by making the approximation that $C(t)\approx C_{\text{max}}$ early in infection. This yields approximations for $g_w(t)$, $g_m(t)$, and $V_w(t)$:

(21) $g_w(t),g_m(t)\approx g_0\equiv {ℛ}_0{d}_v$

(22) ${\displaystyle V_w(t)\{\begin{array}{ll}b\exp ((g_0-d_v)t)& t\le t_M\\ \\ b\exp ((g_0-d_v)t_M)\exp ((g_0-d_v-c_{w}k)(t-t_M))& t>t_M\end{array}}$

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

Figure 8

Download asset Open asset

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.

Probability of replication selection

Request a detailed protocol

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 ${ℛ}_0^w>1$ early in infection:

(27) $p_{\text{repl}}(a,t)=P(t_e<t^{*}(a,t))=1-e^{-\mathrm{\Lambda }(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)\approx C_{\text{max}}$ early infection approximations that give us a closed form for ${\mathrm{\Lambda }}_m(x)$.

Figure 9

Download asset Open asset

When $t^{*}>t_M$, the integral ${\int }_0^{t^{*}}{\lambda }_m(x)𝑑x$ can be evaluated piecewise, first from 0 to $t_M$, and then from $t_M$ 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 $t_e<t_M$, the expression for $p_{\text{repl}}$ in practice depends only on $\delta =(c_w-c_m)k$, not on c_w, c_m and k separately. In the absence of an antibody response, the mutant is generated with near certainty by 1 day post-infection ($P(t_e<1)\approx 1$, Figure 8B). So when $t_M\ge 1$, values for $p_{\text{repl}}$ calculated with ${\displaystyle c_w=1,c_m=0}$, and $\delta =k$—as in (Figure 1C,D)—will in fact hold for any c_w, c_m, and k that produce that fitness difference $\delta$.

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

(28) $p_{\text{repl}}\approx \frac{{ℛ}^w(0)}{1-{ℛ}^w(0)}b\mu p_{\text{sse}}$

See Appendix Section A3.6 for a derivation.

In this section, we describe our model of the point of transmission, including how sIgA neutralization may impose selection pressure.

With these analyses in hand, it is possible to combine the within-host and the point of transmission processes to calculate an overall probability that a new variant survives the transmission bottleneck. Equation 36 gives the probability of bottleneck survival given $f_{\text{mt}}$ and the properties of the recipient host. And given a time of transmission t_t, we can calculate the probability distribution of $f_{\text{mt}}$ using our expressions for the CDF of successful mutation times t_e (Equation 20) and for $f_m(t)$ (Equations 13, 15). To calculate the overall probability $p_{\text{nv}}$, we average over the possible values of $f_{\text{mt}}$, weighted by their probability:

To evaluate the relative probabilities of replication selection and inoculation selection for antigenic novelty and to check the validity of the analytical results, we simulated 10⁶ transmissions from a naive host to a previously immune host. The transmitting host was simulated for 10 days, which given the selected model parameters is sufficient time for almost all infections to be cleared. Time of transmission was randomly drawn from that period, weighted by transmission probability, and an inoculum was drawn from the within-host virus population at that time. Variant counts in the inoculum were Poisson-distributed with probability equal to the variant frequency within the transmitting host at the time of transmission. We simulated the recipient host until the clearance of infection and found the maximum frequency of transmissible variant: that is, the variant frequencies when the transmission probability was greater than $5\times {10}^{-2}$ (probabilistic model) or when the virus population was above the transmission threshold $\theta$ (threshold model). For Figure 3D, we defined an infection with an emerged new antigenic variant as an infection with a maximum transmissible new variant frequency of greater than 50%.

To study evolution along transmission chains with mixed host immune statuses, we modeled the virus phenotype as existing in a 1-dimensional antigenic space. Host susceptibility s_x to a given phenotype x was $s_x=\min \{1,\min \{\left|x-y_i\right|\}\}$ for all phenotypes y_i in the host’s immune history.

Within-host cross immunity $c_{ij}$ between two phenotypes y_i and y_j was equal to $\max \{0,1-\left|y_i-y_j\right|\}$. When mucosal antibodies were present, protection against infection z_x was equal to the strength of homotypic protection $z_{\text{max}}$ scaled by susceptibility: ${\displaystyle z_x=(1-s_x)z_{\text{max}}}$. ${\displaystyle {\kappa }_x}$ was calculated from z_x by $\exp (-v(1-{\kappa }_x))=z_x$. We used ${\displaystyle z_{\text{max}}=0.95}$. Note that this puts us in the regime in which intermediately immune hosts are the best inoculation selectors (Figure 3H). We set ${\displaystyle k=25}$ so that there would be protection against reinfection in the condition with an immediate recall response but without mucosal antibodies.

We then simulated a chain of infections as follows. For each inoculated host, we tracked two virus variants: an initial majority antigenic variant and an initial minority antigenic variant. If there were no minority antigenic variants in the inoculum, a new focal minority variant (representing the first antigenic variant to emerge de novo) was chosen from a Gaussian distribution with mean equal to the majority variant and a given variance, which determined the width of the mutation kernel (for results shown in Figure 5, we used a standard deviation of 0.08). We simulated within-host dynamics in the host according to our within-host stochastic model.

We founded each chain with an individual infected with all virions of phenotype 0, representing the current old antigenic variant.

We simulated contacts at a fixed, memoryless contact rate ${\displaystyle \rho =1}$ contacts per day. Given contact, a transmission attempt occurred with a probability proportional to donor-host viral load, as described above. If a transmission attempt occurred, we chose a random immune history for our recipient host according to a pre-specified distribution of host immune histories. We then simulated an inoculation and, if applicable, subsequent infection, according to the within-host model described above. If the recipient host developed a transmissible infection, it became a new donor host. If not, we continued to simulate contacts and possible transmissions for the donor host until recovery. If a donor host recovered without transmitting successfully, the chain was declared extinct and a new chain was founded.

We iterated this process until the first phenotypic change event—a generated or transmitted minority phenotype becoming the new majority phenotype. We simulated 1000 such events for each model and examined the observed distribution of phenotypic changes compared to the mutation kernel.

For the results shown in Figure 5, we set the population distribution of immune histories as follows: 20% −0.8, 20% −0.5, 20% 0.0, and the remaining 40% of hosts naive. This qualitatively models the directional pressure that is thought to canalize virus evolution (Bedford et al., 2012) once a cluster has begun to circulate.

To assess the causes of the observed behavior in our transmission chain model, we also studied analytically how replication and inoculation selection determine the distribution of observed fixed antigenic changes given the mutation kernel when a host with one immune history inoculates another host with a different immune history.

We fixed a transmission time $t=2$ days, roughly corresponding to peak within-host virus titers. For each possible new variant phenotype, we calculated $p_{\text{nv}}$ according to Equation 43, with parameters given by the old variant antigenic phenotype, new variant phenotype, and host immune histories. Finally, we multiplied each phenotype’s survival probability by the same Gaussian mutation kernel used in the chain simulations (with mean 0 and s.d. 0.08), and normalized the result to determine the predicted distribution of surviving new variants given the mutation kernel and the differential survival probabilities for different phenotypes.

To evaluate the probability of variants being selected and proliferating during a local influenza virus epidemic, we first noted that the per-inoculation rate of new antigenic variant infections for a population with n_s susceptibility classes (which can range from full susceptibility to full immunity) is:

where ${\kappa }_w(i)$ and ${\kappa }_m(i)$ are the mucosal antibody neutralization probabilities for the old antigenic variant and the new antigenic variant associated with susceptibility class i, and $s_i(0)=\frac{S_i(0)}{N}$ is the initial fraction of individuals in susceptibility class i.

We then considered a well-mixed population with frequency-dependent transmission, where infected individuals from all susceptibility classes are equally infectious if infected. Using an existing result from epidemic theory (Magal et al., 2018), we calculated $R_{\mathrm{\infty }}$, the average fraction of individuals who are infected if an epidemic occurs in such a population. During such an epidemic, each individual will on average be inoculated (challenged) ${\Re }_0{R}_{\mathrm{\infty }}$ times, where ${\Re }_0$ is the population-level basic reproduction number (Miller, 2012). We can then calculate the probability that a new variant transmission chain is started in an arbitrary focal individual:

We assessed the sensitivity of our results to parameter choices by re-running our simulation models with randomly generated parameter sets chosen via Latin Hypercube Sampling from across a range of biologically plausible values. Appendix 1—figure 3 gives a summary of the results.

We simulated 50,000 infections of experienced hosts ($c_w=1$) according to each of 10 random parameter sets. We selected parameter sets using Latin Hypercube sampling to maximize coverage of the range of interest without needing to study all possible permutations. We did this for the following bottleneck sizes: 1, 3, 10, 50.

We analyzed two cases: one in which the immune response is unrealistically early and one in which it is realistically timed. In the unrealistically early antibody response model, $t_M$ varied between ${\displaystyle t_M=0}$ and ${\displaystyle t_M=1}$. In the realistically-timed antibody response model, $t_M$ varied between ${\displaystyle t_M=2}$ and ${\displaystyle t_M=4.5}$. Other parameter ranges were shared between the two models (Table 2).

Discussion of sensitivity analysis results can be found in Appendix Section A8.

We downloaded processed variant frequencies and subject metadata from the two NGS studies of immune-competent human subjects naturally infected with A/H3N2 with known vaccination status (Debbink et al., 2017; McCrone et al., 2018) from the study Github repositories and journal websites. We independently verified that reported antigenic site and antigenic ridge amino acid substitutions were correctly indicated, determined the number of subjects with no NGS-detectable antigenic amino acid substitutions, and produced figures from the aggregated data.

For within-host model stochastic simulations, we used the Poisson Random Corrections (PRC) tau-leaping algorithm (Hu and Li, 2009). We used an adaptive step size; we chose step sizes according to the algorithm of Hu and Li, 2009 to ensure that estimated next-step mean values were non-negative, with a maximum step size of 0.01 days. Variables were set to zero if events performed during a timestep would have reduced the variable to a negative value. For the sterilizing immunity simulations in Figure 2, we used a smaller maximum step size of 0.001 days in recipient hosts to better handle mutation dynamics involving very small numbers of replicating virions.

We obtained numerical solutions of equations, including systems of differential equations and final size equations, in Python using solvers provided with SciPy (Jones et al., 2001).

All code, data, and other materials needed to reproduce the analysis in this paper are provided online on the project Github repository: https://github.com/dylanhmorris/asynchrony-influenza (Morris, 2020; copy archived at swh:1:rev:5a9796fa3ab7b8a86aeccd7c9353542f9409e215).

They are also available on OSF: https://doi.org/10.17605/OSF.IO/jdsbp.

Output data generated by stochastic simulations, within-host NGS meta-analysis, and phylogenetic analyses are archived on OSF: https://doi.org/10.17605/OSF.IO/jdsbp.

This appendix contains additional information for ‘Asynchrony between virus diversity and antibody selection limits influenza virus evolution’.

We provide a detailed review of virological and immunological evidence supporting our model of the interaction between influenza viruses and the (human) host immune system (Section A2). We give a more detailed mathematical analysis of our within-host model, deriving the approximate and exact analytical results that are used in the main text, and provide intuition to explain model behavior and limitations (Section A3). We do the same for our model of the point of transmission, and assess which hosts are most likely to act as inoculation selectors and which mutants are most likely to be inoculation selected (Section A4).

We next discuss parameter uncertainties and report the results of a sensitivity analysis for our within-host models (Section A5). We report an additional empirical analysis showing that new variants often emerge polyphyletically: they appear simultaneously on multiple branches of the influenza virus phylogeny. We discuss how this is more easily explained in light of our model than in light of previous models (Section A6).

We then review these prior studies and models of within-host influenza virus antigenic evolution and general pathogen immune escape in depth. We argue that they are insufficient to explain within-host influenza virus antigenic evolution and that inoculation selection accompanied by a realistically-timed antibody response is the most likely competing hypothesis (Section A7). We elaborate on the conclusions that can be drawn from our study and its implications for influenza virus control and for future basic research (Section A8). We conclude by providing complete step-by-step mathematical derivations of lemmas and approximations from elsewhere in the text (Section A9).

In this section, we review relevant immunology to establish the biological realities that we built our model to capture. We also discuss some unmodeled biological complexities, and why they are unlikely to change our qualitative conclusions.

Here we analyze the within-host model in depth, and find expressions for within-host variant frequencies over time, the probability distribution for the time of first successful mutation to a new variant, and an approximate probability of replication selection to a given frequency by a given time.