Abstract

Recent phylogenetic analyses indicate that RNA virus populations carry a significant deleterious mutation load. This mutation load has the potential to shape patterns of adaptive evolution via genetic linkage to beneficial mutations. Here, we examine the effect of deleterious mutations on patterns of influenza A subtype H3N2's antigenic evolution in humans. By first analyzing simple models of influenza that incorporate a mutation load, we show that deleterious mutations, as expected, act to slow the virus's rate of antigenic evolution, while making it more punctuated in nature. These models further predict three distinct molecular pathways by which antigenic cluster transitions occur, and we find phylogenetic patterns consistent with each of these pathways in influenza virus sequences. Simulations of a more complex phylodynamic model further indicate that antigenic mutations act in concert with deleterious mutations to reproduce influenza's spindly hemagglutinin phylogeny, co-circulation of antigenic variants, and high annual attack rates.

eLife digest

Each year, up to 15% of the world's population experience symptoms of an influenza infection, also commonly known as flu. The most common culprit is a strain of the virus called influenza type A subtype H3N2. One reason that so many people become infected each year is that this virus evolves rapidly. Within a few years, proteins on the surface of the virus known as antigens become less recognizable to the immune system of a person who has been previously infected. This means that the person can become ill with the virus again because their immune system cannot mount an effective response to the evolved virus strain.

Influenza virus strains evolve rapidly because their genetic material accumulates mutations quickly. Although some of these mutations are beneficial to the virus, other mutations are harmful and reduce the ability of the virus to spread. Sometimes beneficial mutations may occur alongside harmful ones, but it is not known how the harmful mutations affect the evolution of the virus.

Here, Koelle and Rasmussen used computer models of H3N2 influenza to examine the effect of harmful mutations on the evolution of this virus population. The models show that harmful mutations limit how quickly the antigens can evolve. Also, the presence of these harmful mutations effectively acts as a sieve: they allow only large changes in the antigens to establish in the virus population.

The models suggest that there are three routes by which large changes in the antigens on H3N2 viruses may occur. The first is by a single mutation that has a big effect on the antigens in viruses that only carry a few harmful mutations, but these large mutations would not happen very often. Another route may be through more common mutations that have only a small or moderate benefit, which would allow the virus to become more common in the population before it acquires a beneficial mutation with a much greater effect. The third possibility is that a large beneficial mutation may arise in viruses that have many harmful mutations. These harmful mutations may initially limit the ability of the virus to spread, but over time, some of these harmful mutations may then be lost.

Koelle and Rasmussen found that the computer models could recreate the patterns of virus evolution that have been observed in real strains of H3N2. Researchers use predictions of influenza evolution to help them decide which virus strains should be included in flu vaccines each year. Koelle and Rasmussen findings indicate that harmful mutations should be considered when making these predictions.

Distinct from these efforts, several phylogenetic analyses have indicated that influenza A/H3N2 in humans carries a deleterious mutation load (Fitch et al., 1997; Pybus et al., 2007; Strelkowa and Lässig, 2012). Specifically, early work by Fitch et al. (1997) found the number of nonsynonymous changes on tip branches of the HA phylogeny to be higher than expected, indicative of either strain selection bias or the presence of transiently circulating deleterious mutations in the influenza viral population. In more recent work, Pybus et al. (2007) performed a comprehensive phylogenetic analysis of over 140 viruses, including influenza A/H3N2. For H3N2's M1 protein, as well as for the majority of the other viral proteins examined in the study, they found heightened ratios of non-synonymous-to-synonymous substitutions on external tree branches relative to those found internally. This finding again points towards transiently circulating deleterious mutations in influenza and, more generally, across RNA virus populations. Other recent work on predicting the short-term evolution of influenza has highlighted the necessity of accounting for fitness costs associated with sublethal deleterious mutations when projecting the frequencies of influenza clades into the next season (Łuksza and Lässig, 2014). Together, these results indicate that purifying selection is not sufficiently strong to immediately eliminate deleterious mutations from the influenza A/H3N2 virus population that circulates among humans. As a result of genetic linkage within genes and, to a lesser extent, across genes, these deleterious mutations have the potential to interact with beneficial mutations in determining which viral lineages will persist and which ones will ultimately be lost. Indeed, a recent statistical analysis of HA sequences from influenza A/H3N2 has suggested that interference effects largely determine the fates of viral mutants, rather than their inherent selective effects (Illingworth and Mustonen, 2012). These interference effects are possible because of an extensive genetic linkage across influenza's HA (Strelkowa and Lässig, 2012) and arise from variation in the background fitness of viral strains and from variation in the fitness effects of subsequent mutations.

Taken together, this body of work indicates that sublethal deleterious mutations commonly arise and circulate for sufficiently long periods of time to be able to impact the trajectories of influenza A/H3N2 strains. However, the impact that these deleterious mutations have on the population dynamics and long-term evolutionary patterns of this subtype has to date not been explored. Here, we address this question with a set of increasingly complex population genetic and population dynamic models, under the common assumption that influenza's adaptive evolution is driven by antigenic changes that allow for escape from herd immunity. We start by extending classic population genetic models into an epidemiological context. As expected from previous analyses of these types of models (Fisher, 1930; Birky and Walsh, 1988; Peck, 1994; Barton, 1995; Orr, 2000), we find that circulating sublethal deleterious mutations in influenza A/H3N2's viral population reduce the rate of adaptive evolution and increase the average size of the beneficial mutants that fix. Extending this analysis to models explicitly incorporating epidemiological dynamics, we further show that the accumulation of deleterious mutations can contribute to explaining the invasion dynamics of new antigenic clusters, defined as sets of viral strains that are antigenically similar to one another (Smith et al., 2004). This model further predicts three distinct molecular pathways by which antigenic cluster transitions can occur, and we find empirical patterns consistent with each of these pathways in sequence data from 1992 to 2004.

Gaining intuition from these simple models, we then present more extensive phylodynamic simulations that incorporate the occurrence of both antigenic and non-antigenic (largely deleterious) mutations. This extension is critical given that antigenic mutations acquire their selective advantage through immune escape, which reduces competition for susceptible hosts and therefore in principle could allow for long-term coexistence of virus strains through niche partitioning of the host population (Cobey, 2014). Indeed, in the absence of other contributing processes, it has been shown that reduced between-strain competition resulting from antigenic evolution leads to explosive genetic and antigenic diversity (Ferguson et al., 2003), a pattern that is inconsistent with the long-term evolutionary dynamics of influenza A/H3N2 in humans. Intriguingly, the phylodynamic model we present robustly reproduces the spindly phylogeny of influenza A/H3N2's HA protein (Fitch et al., 1997) under parameterizations relevant to A/H3N2 in humans. It further reproduces the recently described patterns of co-circulation of minor antigenic variants (Strelkowa and Lässig, 2012), as well as high annual attack rates (World Health Organization, 2014). In the discussion, we situate these findings in the context of previously published models used to explain the characteristic evolutionary dynamics of influenza A/H3N2 in humans, speculate on the applicability of these findings to other influenza-host systems, and comment on the consequences of these findings for the predictability of influenza evolution.

Results

A deleterious mutation load modifies the tempo and nature of influenza's antigenic evolution

To develop an understanding for how circulating deleterious mutations will impact patterns of influenza's antigenic evolution, we first extend existing population genetic models (Haigh, 1978; Peck, 1994) to acute infectious diseases undergoing immune escape. As is common in many population genetic models, we assume an infinite population and consider this population subject to frequent sublethal deleterious mutations that act independently from one another in reducing fitness. In an explicit susceptible-infected-recovered (SIR) epidemiological context that does not yet incorporate immune escape, these assumptions lead to a deleterious mutation-selection balance in the virus population (‘Materials and methods’, Figure 1A). This balance is given by a Poisson distribution with mean λ/sd, where λ is per-genome deleterious mutation rate and sd is the transmission fitness cost of sublethal deleterious mutations. The virus population at epidemiological equilibrium will have an overall net reproductive rate (i.e., mean absolute fitness) of R = 1, with more transmissible viruses that carry fewer deleterious mutations having reproductive rates above one and less transmissible viruses that carry more deleterious mutations having reproductive rates below one (Figure 1A, inset). This within-population variation in viral transmission rates is also reflected in the distribution of infected individuals' basic reproductive rates (R0 values) (Figure 1A inset).

The effect of a deleterious mutation load on the fate of an antigenic mutant.

(A) A resident antigenic strain at its deleterious mutation-selection balance. The histogram shows pk, the proportion of infected individuals carrying a virus with k deleterious mutations at its endemic equilibrium. The viral class carrying the fewest number of deleterious mutations is defined as mutation class k = 0. Inset: variation in the basic reproductive rate of infected individuals (gray histogram) and variation in the net reproductive rate R of infected individuals (brown histogram) resulting from variation in the number of deleterious mutations carried by circulating viruses. Model parameters: λ = 0.10, sd = 0.008, μ = 1/30 years−1, γ = 1/4 days−1, R0,k = 0 = 2.25. (B) Simulations showing the three dynamical fates that an antigenic mutant can experience: rapid loss (blue, with expanded inset), transient circulation (green), and successful establishment (red). The blue antigenic mutant (with σ = 0.008) arises in a genetic background with k = 16 deleterious mutations. The green antigenic mutant (with σ = 0.04) arises in a background of k = 10 deleterious mutations. The red antigenic mutant (with σ = 0.06) arises in a background of k = 5 deleterious mutations. All other parameters are as in (A). (C) The proportion of antigenic mutants that result in each of the three dynamical fates shown in (B) as a function of their antigenic size σ. All parameters as in (A). (D) The average antigenic size of successfully establishing antigenic mutants under different deleterious mutation rates λ (x-axis) and for three transmission fitness costs (see (E) for legend). The size of arising antigenic mutations σ are assumed to be gamma distributed with mean of 0.012 and a shape parameter of 2. All other parameters are as in (A). (E) The relative rate of antigenic evolution under different deleterious mutation rates λ (x-axis) and for the three transmission fitness costs shown in (D). Other model parameters are as in (D). The relative rate of antigenic evolution is given by the fraction of arising antigenic mutants that establish under the deleterious mutation load relative to the fraction of arising antigenic mutants that would establish under a no-load scenario. See ‘Materials and methods’ for choice of model parameters.

Following Peck (1994), we first examine the fate of a single advantageous mutant arising in such a population. This advantageous mutant will necessarily arise in a genetic background with a certain number of deleterious mutations and, in our case, carry an immune-escape mutation that is beneficial to its spread. The genetic background in which the antigenic mutation arises and the size of the antigenic change jointly determine the mutant's initial reproductive rate Rm(t = 0) in the virus population (‘Materials and methods’). If this initial reproductive rate is less than one, the antigenic mutant is likely to be rapidly lost from the virus population. If the mutant's initial reproductive rate is instead greater than one, the mutant will invade the virus population if it is not initially stochastically lost. In the case of invasion, the mean reproductive rate of the antigenic mutant lineage will necessarily decline because it will accumulate its own set of deleterious mutations. (The mutant lineage's mean reproductive rate will also necessarily decline because the size of its susceptible host pool will decline over time, a factor we for now ignore but return to in later models.) For this invading antigenic mutant lineage, we can calculate a final reproductive rate Rm(t = ∞), which is the mean reproductive rate of this lineage once it has reached its own mutation-selection balance (‘Materials and methods’). If this final reproductive rate falls below one, an invading antigenic mutant will therefore only transiently circulate before deterministically declining as a result of deleterious mutation accumulation. If this final reproductive rate exceeds one, however, the invading antigenic mutant lineage, under the assumptions of this model, will successfully establish. Antigenic mutants therefore experience one of three possible fates (Figure 1B): rapid loss, transient circulation, or successful establishment.

Which of these three fates awaits an antigenic mutant depends in part on the number of deleterious mutations carried by the strain in which the antigenic mutation arises: the lower the number of background deleterious mutations, the higher the antigenic mutant's chances are of at least transient establishment, consistent with the background fitness interference effects found in Illingworth and Mustonen (2012). Which fate occurs also depends on the extent to which the antigenic mutant escapes immunity (Figure 1C). The vast majority of small-sized antigenic mutants are rapidly lost; the remaining ones only circulate transiently before the accumulation of deleterious mutations results in their ultimate loss. Antigenic mutants that significantly escape immunity are less subject to rapid loss and also less likely to circulate only transiently. Given a specified size distribution for antigenic mutations, the average size of antigenic mutations that will successfully establish can be calculated. Figure 1D shows that this average antigenic size increases with increases in the deleterious mutation rate λ. The magnitude of the transmission fitness cost sd also affects the average size of antigenic mutants that will successfully establish. For any given deleterious mutation rate λ, as sd decreases the average size of antigenic mutants that fix increases (Figure 1D). These results can be interpreted in the context of findings from the population genetics literature: increases in λ and decreases in sd similarly increase the virus population's fitness variance (as quantified by the variance in net reproductive rate R, Figure 1A inset, ‘Materials and methods’). Increases in the fitness variance of asexual populations makes fixation of a beneficial mutant increasingly dependent on genetic background; only beneficial mutants that exceed a characteristic large size will have a high probability of fixing in populations with substantial fitness variance (Peck, 1994; Barton, 1995; Schiffels et al., 2011; Good et al., 2012).

In addition to their effect on the sizes of successfully establishing antigenic mutants, circulating deleterious mutations will act to slow the tempo of antigenic evolution (Figure 1E); that is, they will reduce the number of antigenic mutants that go to fixation in a given amount of time. This particular effect has previously been remarked upon in the context of a population genetics model for influenza's HA protein (Strelkowa and Lässig, 2012). Again, increases in the fitness variance of the viral population is the culprit: increases in the deleterious mutation rate λ and decreases in the fitness cost of deleterious mutations sd similarly act to increase fitness variance; with increased fitness variance in the population, the genetic background in which an antigenic mutant needs to arise in to have a chance at fixation will be increasingly restrictive. This leads to a reduced tempo of antigenic change, consistent with a reduction in the rate of adaptation that is known from the population genetics literature (Peck, 1994; Barton, 1995).

Our model's findings can now be situated in the context of influenza A/H3N2's characterized evolutionary dynamics in humans. In particular, detailed antigenic analyses have demonstrated that this virus undergoes punctuated antigenic evolution, with predominantly single amino acid changes of large antigenic effect being responsible for the occurrence of antigenic cluster transitions (Smith et al., 2004; Koel et al., 2013). This dynamic is consistent with our results that antigenic evolution—as traced by lineages that ultimately fix—should occur via mutations of characteristically large size. These same analyses, as well as others (Plotkin et al., 2002), have further indicated that antigenic cluster transitions occur only every 2 to 6 years, an incredibly slow pace given the virus's high mutation rate and the need for only a single amino acid to substantially alter antigenicity. Again, this dynamic is consistent with our results that the tempo of antigenic evolution should be severely reduced by circulating deleterious mutations. Thus, our model posits that the punctuated and surprisingly slow nature of influenza A/H3N2's antigenic evolution are related features of this largely asexual, adapting population subject to a deleterious mutation load: adaptive evolution requires not only that large antigenic mutants occur, but also that they occur in good genetic backgrounds.

The above model contains a number of simplifying assumptions. Among these are that the susceptible host pool is negligibly affected over the time period of the antigenic mutant's establishment and that the successful establishment of an antigenic mutant will lead to the fixation of the mutant lineage in the virus population. Although these assumptions may be reasonable under some scenarios, they may not always be in an epidemiological context. For example, an antigenic mutant might invade sufficiently slowly to erode its frequency-dependent advantage prior to the exclusion of the existing antigenic lineage. Due to only partial cross-immunity between these lineages, this would lead to long-term coexistence of the variants. Another scenario is one of an antigenic mutant with a particularly large selective advantage: this mutant might burn through its susceptible host population so rapidly that its net reproductive rate drops significantly below one, leading to the possibility of its own extinction along with that of the previously circulating variant (Ballesteros et al., 2009). This scenario underscores the importance of considering the possibility of a variable virus population size; population genetic models that assume a constant population size may not be appropriate under certain epidemiological conditions.

To relax both the assumption of a time-invariant selective advantage and the assumption that successful establishment of a mutant leads to fixation of the mutant lineage (including replacement of the resident lineage), we now consider a more complex epidemiological model that explicitly incorporates the dynamics of susceptible hosts (‘Materials and methods’). Figure 2 shows the dynamics of the resident strain and the antigenic mutant under four distinct scenarios: a scenario in which a small antigenic mutant arises in a low-load (‘good’) genetic background (Figure 2A); a scenario in which a small antigenic mutant arises in an average-load genetic background (Figure 2B); a scenario in which a large antigenic mutant arises in a good genetic background (Figure 2C); and a scenario in which a large antigenic mutant arises in an average genetic background (Figure 2D). These simulations indicate that the three fates predicted in the simpler model still play out in the explicit context of epidemiological dynamics when parameterized for influenza A/H3N2 in humans. Rapid loss is expected to occur when a small antigenic mutant arises in an average genetic background (Figure 2B). As this combination occurs commonly, rapid loss is the most frequent fate experienced by antigenic mutants. Transient circulation occurs when a small antigenic mutant arises in a good genetic background (Figure 2A), provided that it survived genetic drift. Transient circulation also occurs when a large antigenic mutant arises in an average genetic background (Figure 2D), again provided that it survived genetic drift. Successful establishment can only occur when a large antigenic mutation arises in a good genetic background (Figure 2C), a ‘jackpot’ combination. In this case, the resident strain is competitively excluded as a result of strain cross-immunity. Intriguingly, the presence of deleterious mutations can also affect the invasion dynamics of antigenic mutants having this ‘jackpot’ combination: because offspring of antigenic mutants accumulate deleterious mutations, and these deleterious mutations reduce viral transmissibility, the invasion dynamics of particularly large antigenic mutants are considerably less explosive than would be expected in the absence of deleterious mutation accumulation (Figure 3). Consequently, large antigenic mutants do not readily burn themselves out during attempted establishment, as might in theory be expected (Ballesteros et al., 2009).

Epidemiological dynamics following the emergence of an antigenic mutant.

First row: antigenic mutants arising in low-load genetic backgrounds (k = 3 deleterious mutations). Second row: antigenic mutants arising in average-load genetic backgrounds (k = 13 deleterious mutations). Left column: antigenic mutants that are of a small size (σ = 0.004). Right column: antigenic mutants that are of a large size (σ = 0.045). (A) Small antigenic mutants arising in good genetic backgrounds transiently circulate. (B) Small antigenic mutants arising in average genetic backgrounds are rapidly lost. (C) Large antigenic mutants arising in good genetic backgrounds can successfully establish and exclude the resident antigenic strain, resulting in an antigenic cluster transition. (D) Large antigenic mutants arising in average genetic backgrounds transiently circulate. All simulations assume a host population size of N = 4 billion and start with a single antigenic strain at evolutionary and epidemiological equilibrium. The remaining parameters are as in Figure 1A.

Explosiveness in cluster invasion dynamics in the presence and absence of deleterious mutation accumulation.

(A) A representative example of the population dynamics of an antigenic mutant in the presence of deleterious mutation accumulation. The new antigenic strain invades and successfully establishes, while excluding the resident strain, characteristic of a successful antigenic cluster transition. Model parameters are N = 4 billion, μ = 1/30 years−1, γ = 1/4 days−1, λ = 0.10, sd = 0.008, R0,k = 0 = 2.25, and σ = 0.05. In this simulation, the antigenic mutant arises in a genetic background with k = 4 deleterious mutations. (B) A representative example of the population dynamics of an antigenic mutant in the absence of deleterious mutations. The new antigenic strain invades explosively, leading to its own burn-out along with exclusion of the resident strain. Model parameters are N = 4 billion, μ = 1/30 years−1, γ = 1/4 days−1, R0 = 2.04, and σ = 0.121. The value of R0 was chosen such that, prior to the invasion of the antigenic mutant, the fraction of the host population susceptible to infection and the number of infected hosts was the same across the two simulations. In (B), the value of σ was chosen such that Rm(t = 0) was the same across the two simulations (at a value of 1.13).

Figure 2C shows that an antigenic mutant can exclude a resident antigenic strain, provided that it arises in a good genetic background and carries an antigenic mutation of large effect. This is consistent with detailed molecular studies of influenza A/H3N2 that have shown that single amino acid changes of large antigenic effect can precipitate an antigenic cluster transition (Smith et al., 2004; Koel et al., 2013), although the importance of the genetic background in which the antigenic mutation arises has not been discussed in the context of this work. The BE92-to-WU95 cluster transition, precipitated by an amino acid change from N to K at site 145, is a good example of a cluster transition occurring via a single mutational step (Koel et al., 2013) (Figure 4A). Of note, several clades witnessed the 145NK amino acid substitution before it occurred in the clade that founded the WU95 viral lineage. Our models above indicate that the failure of these early 145K clades to establish could in principle be a consequence of this N-to-K amino acid change occurring in insufficiently good genetic backgrounds, as also suggested in Neher et al. (2014). Other explanations for the failure of early 145K clades to establish are of course possible. One such explanation is that these clades might have circulated in spatial locations that are not sufficiently well-connected globally. Recent work has further emphasized the important role that spatial ecology plays in the global establishment of antigenic variants (Russell et al., 2008; Lemey et al., 2014; Bedford et al., 2015), with findings suggesting that Asia plays a dominant role in sourcing these new variants. Thus, if the early 145K clades were not geographically well-situated, the spatial context (rather than the genetic context) of these clades may have led to their failure in establishing. We thus examined the spatial locations of the sequences from the three largest 145K clades that failed to establish: all three clades were geographically widespread, spanning at least two continents. Furthermore, two of the three clades contained sequences from Asia, despite Asia being undersampled during this time period. It is thus unlikely that these early 145K clades were geographically restricted to ‘sink’ populations. An alternative explanation for the failure of these early 145K clades to establish is that herd immunity levels against BE92 may not yet have been high enough to result in a sufficiently large selective advantage for these WU95-like lineages. Given these alternative explanations, an in vitro experimental study that quantifies relative viral fitness of these 145K lineages would thus be necessary to confirm our genetic background hypothesis.

While hitting the ‘jackpot’ combination is a viable molecular pathway for precipitating an antigenic cluster transition, the combination of a large antigenic mutation arising in a low-load genetic background is expected to occur infrequently. We can therefore consider whether there might be alternative molecular pathways open to antigenic mutants that would similarly yield a successful cluster transition. One possibility is for a more common small- to medium-sized antigenic mutation to first occur in a good genetic background. This would result in a transient rise of this mutant (Figure 2A), thereby increasing the number of individuals infected with the virus carrying few deleterious mutations. A less common, large-sized antigenic mutation could then occur, effectively piggy-backing on the good genetic background that the smaller antigenic mutant inflated. Because smaller antigenic mutations can only circulate transiently, and accumulate deleterious mutations during their circulation, the large antigenic mutation must not only follow, but also rapidly follow, the rise of the smaller antigenic mutation for this molecular pathway to yield an antigenic cluster transition. This scenario is depicted in Figure 5A, and phylogenetically would result in a sudden appearance of a viral lineage carrying two antigenic mutations. A phylogenetic analysis of the WU95-SY97 antigenic cluster transition provides an empirical example that is consistent with this molecular pathway of antigenic turnover, with a seemingly simultaneous accumulation of two antigenic amino acid changes (156KQ and 158EK) occurring on a short branch of the reconstructed phylogeny (Koel et al., 2013) (Figure 4B).

A third molecular pathway that could in principle precipitate an antigenic cluster transition is for a large antigenic mutation to first arise in an average genetic background and, during its transient circulation (Figure 2D), to purge itself of a large number of deleterious mutations. This purging could arise through within-subtype viral reassortment taking place in an individual coinfected with a strain belonging to the resident cluster and a strain belonging to the transiently circulating antigenic mutant. Even though both of the strains infecting this individual would likely carry an average deleterious mutation load, it is highly unlikely that they will carry the same set of deleterious mutations if phylogenetically sufficiently far apart. Reassortment of the eight gene segments within the coinfected host could therefore significantly lower the number of deleterious mutations carried by viral progeny characterized as belonging to the new antigenic cluster. Once the deleterious mutation load has largely been shed, the reassortant virus would quickly rise and cause an antigenic cluster transition (Figure 5B). Indeed, many historical instances of intrasubtypic reassortment contributing to antigenic turnover have been documented (Morens et al., 2009); these instances have been associated with high incidence levels as would be expected by purging of deleterious mutations. A phylogenetic analysis of the SY97-FU02 cluster transition provides an especially compelling example that is consistent with this molecular pathway of antigenic turnover. In this cluster transition, a virus antigenically characterized as FU02 reassorted with a genetically distant virus antigenically characterized as SY97 (Barr et al., 2005; Holmes et al., 2005) (Figure 4C). This reassortant viral lineage circulated extensively in Australia and New Zealand in 2003 and subsequently in the US in the 2003–2004 influenza season, causing substantial morbidity and mortality (Barr et al., 2005). Although this viral lineage may ultimately have led to the replacement of the non-reassortant FU02 viral lineage, both of these FU02 lineages were excluded by the subsequent CA04 antigenic cluster, which appears to have originated from the non-reassortant FU02 lineage (Figure 4C). Intriguingly, the CA04 lineage carried with it not only an HA mutation of large antigenic size, but also two amino acid changes in its polymerase gene segment that enhanced replicative fitness (Memoli et al., 2009).

Antigenic evolution in the context of fitness variation generated by deleterious mutations results in a spindly HA phylogeny, antigenic variant co-circulation, and high annual attack rates

Our above analyses have relied on simple epidemiological models to gain intuition for how circulating sublethal deleterious mutations would impact patterns of influenza A/H3N2's antigenic evolution. These analyses indicate that deleterious mutation loads should lower the rate of antigenic evolution (Figure 1E). Based on previous modeling work (Koelle et al., 2006, 2009, 2010; Zinder et al., 2013), a lower rate of antigenic evolution is known to constrain genetic and antigenic diversity and thus might lead to a spindly HA phylogeny. Our analyses also indicate that deleterious mutation loads increase the average size of antigenic variants that establish in the long run (Figure 1D); observed patterns of punctuated antigenic evolution (Smith et al., 2004) may thus be better reproduced with a model that integrates sublethal deleterious mutations than one that ignores these mutations. Furthermore, our above analyses indicate that antigenic variants can reach appreciable numbers even when their ultimate fate is one of only transient circulation. This pattern of transient circulation points towards the possibility of co-circulation of a substantial number of antigenic variants, as argued for in statistical analyses of influenza's HA sequences (Strelkowa and Lässig, 2012). In our model, these antigenic variants need not necessarily strongly compete with one another for susceptible hosts for only a single lineage to persist. In the absence of exceptionally strong competition for susceptible hosts, the co-circulation of these antigenically distinct variants may thus be capable of reproducing empirically observed high annual attack rates (World Health Organization, 2014).

To determine whether spindly HA phylogenies, co-circulation of antigenic variants, and high annual attack rates indeed come out of a model that incorporates deleterious mutations, we implemented a phylodynamic model that simulates the occurrence of both non-antigenic (largely deleterious) mutations and antigenic mutations (‘Materials and methods’). When simulated under parameters appropriate for influenza A/H3N2 in humans, this model yields a spindly HA phylogeny with low-load viruses populating the trunk and higher load viruses populating the tips of the phylogeny (Figure 6A). This distribution of deleterious mutation loads on the simulated phylogeny is consistent with the excess number of non-synonymous substitutions on external tree branches that was found for human influenza A/H3N2 (Fitch et al., 1997; Pybus et al., 2007) and arises because deleterious mutations contribute a substantial fraction of the fitness variance in the viral population (Figure 7A).

Viral phylogenies from a simulation of the phylodynamic model incorporating antigenic and non-antigenic mutations.

(A) Simulated phylogeny reproducing H3N2's spindly HA phylogeny, with low levels of genetic and antigenic diversity over the long run. Lineages are colored according to their deleterious mutation loads. (B) Simulated phylogeny shown in (A) with lineages colored according to their antigenic type. Similarly colored lineages that are genetically distinct are antigenically distinct (colors were re-used due to their limited number). (C) Simulated phylogeny shown in (A) with lineages colored according to the fraction of the host population susceptible to infection with that lineage (Seff/N). (D) Simulated phylogeny shown in (A) with lineages colored according to their net reproductive rate R.

Fitness variance dynamics, epidemiological dynamics, and times to most recent common ancestor for a simulation of the phylodynamic model incorporating antigenic and non-antigenic mutations.

(A) The fraction of (log) viral fitness variation explained by deleterious mutations over time for the simulation whose phylogenies are plotted in Figure 6. The fraction not explained by mutation load is due to antigenic variation in the population. (B) Simulated epidemiological dynamics showing co-circulation of multiple antigenic variants and sustained prevalence levels over time for this same simulation. (C) Times to the most recent common ancestor (tMRCAs), computed over time from the phylogenies shown in Figure 6.

Despite extensive co-circulation of antigenic variants, the overall phylogeny remains ladder-like due to the selective sweeps initiated by rare large-sized antigenic mutations that arise in good genetic backgrounds. This dynamic is evident by jointly considering Figure 6A and Figure 6C, with Figure 6C showing, for each lineage, the fraction of the host population that is susceptible to infection with that lineage. From these figures, it is clear that the trunk of the phylogeny carries low-load viruses (Figure 6A) that have a high number of susceptible hosts (Figure 6C). This combination together yields high-fitness viruses, as epidemiologically given by their net reproductive rates R (Figure 6D). It is these viruses that establish and thus form the trunk of the tree. Neither a low mutation load alone nor a high number of susceptible hosts alone suffice in generating a sufficiently high-fitness viral lineage that will ensure its long-term evolutionary success.

Inspection of Figure 6C also shows that trunk lineages abruptly gain susceptible hosts. These abrupt gains are a result of single large antigenic mutations that, when occurring in good genetic backgrounds, initiate the selective sweeps that ultimately limit influenza's genetic diversity. How these selective sweeps affect levels of standing genetic diversity in the viral population is shown in Figure 7C, where we use the time to the most recent common ancestor (tMRCA) of all circulating lineages as a proxy for total genetic diversity. This tMRCA plot reproduces quantitative features of influenza A/H3N2's tMRCA plot presented in Bedford et al. (2011), including its interannual variation and the observed major drops in tMRCA following the emergence of new and antigenically very distinct clusters.

The importance of deleterious mutations in constraining the genetic and antigenic diversity of influenza can be further illustrated by simulating the phylodynamic model under the assumption of their absence. These simulations very rapidly yield explosive genetic and antigenic diversity (Figure 8A) and, as a consequence, prevalence levels that generally increase over time (Figure 8B). Reassuringly, this finding is consistent with predictions from previous influenza modeling work incorporating only strain-specific immunity (Ferguson et al., 2003).

Given that purifying selection alone is known to reduce genetic diversity (Charlesworth et al., 1993; Walczak et al., 2012), we further simulated the phylodynamic model under the assumption of no antigenic mutations. In these simulations, we phenomenologically incorporated antigenic drift by simulating susceptible-infected-recovered-susceptible (SIRS) dynamics such that prevalence levels were similar to those shown in Figure 7B. Compared with the results shown in Figure 6, these simulations gave rise to significantly higher levels of genetic diversity (Figure 9A) and longer tMRCAs (Figure 9B). This indicates that purifying selection alone does not account for the spindly phylogenies shown in Figure 6. Rather, it is antigenic evolution in the context of fitness variation generated by deleterious mutations that constrains the viral phylogeny.

Simulations of the phylodynamic model in the absence of antigenic mutations.

(A) Simulated phylogeny under a parameterization with no antigenic mutations, showing genetic diversity generally increasing over time. (B) Times to the most recent common ancestor (tMRCAs) computed over time from the phylogeny shown in (A).

While the above simulations demonstrate that a model incorporating both antigenic and deleterious mutations can reproduce influenza A/H3N2's spindly phylogeny, its antigenic co-circulation patterns, and its high annual attack rates, a number of the parameters that required specification are empirically not well characterized. Most notably, these are the evolutionary parameters of the model: the deleterious mutation rate λ, the fitness cost of deleterious mutations sd, and the antigenic mutation rate λantigenic. In Figure 10, we show how the evolutionary and epidemiological dynamics of the model simulations depend on these three parameters. Specifically, we vary one parameter while keeping the remaining two parameters fixed. For each parameter set considered, we perform 20 simulations to determine the range of dynamics predicted under a single parameterization. For each simulation, we quantify the virus's evolutionary dynamics by plotting the minimum and maximum tMRCA calculated over a continuous 10-year period (years 15–25 of the simulation) as a measurement of the extent of genetic diversity in the viral population. To quantify the virus's epidemiological dynamics, we plot the minimum and maximum annual attack rates over this same time period. Figure 10A shows that in the absence of deleterious mutations (λ = 0) the model simulations yield explosive viral diversity, with the maximum tMRCA ranging between 5 and 25 years. The observed increases in viral genetic and antigenic diversity result in unrealistically high maximum annual attack rates, in the range of 15–45% (Figure 10B). Simulated under this parameterization, the results shown in Figure 8A,B are representative of these results. With increasing deleterious mutation rates, the maximum tMRCAs decline as do maximum annual attack rates (Figure 10A,B). For λ values of 0.10 and higher, empirically documented annual attack rates can be consistently reproduced. Maximum and minimum tMRCAs are best reproduced for λ values between 0.10 and 0.15.

Sensitivity of evolutionary and epidemiological dynamics to parameters of the phylodynamic model.

Subplots (A, B) show model sensitivity to the deleterious mutation rate λ. Subplots (C, D) show model sensitivity to the fitness cost of deleterious mutations sd. Subplots (E, F) show model sensitivity to the antigenic mutation rate λantigenic. The top row shows maximum (red dots) and minimum (black dots) times to the most recent common ancestor (tMRCAs) for 20 independent simulations. The red dashed line indicates the maximum tMRCA inferred from a phylogenetic analysis of influenza A/H3N2's HA (Bedford et al., 2011); the black dashed line indicates the minimum tMRCA inferred from this same analysis. The bottom row shows maximum (red dots) and minimum (black dots) annual attack rates for the same 20 simulations. The red dashed line indicates an estimate of the maximum annual attack rate for influenza A/H3N2; the black dashed line indicates an estimate of the minimum annual attack rate for influenza A/H3N2. These values are based on annual attack rate estimates in adults of 5–10%, such that the maximum annual attack rate is on the order of 10%, and the minimum annual attack rate is shown at 1% (which would correspond to years of negligible circulation of this influenza subtype). Each simulation was run for 28 years, and minimum and maximum tMRCAs and attack rates were computed from years 15–25 of the simulation. In subplots (A) and (B), λ is varied, sd = 0.008 and λantigenic = 0.00075. In subplots (C) and (D), λ = 0.10, sd is varied, and λantigenic = 0.00075. In subplots (E) and (F), λ = 0.10, sd = 0.008, and λantigenic is varied. All other parameter values are as listed in Figure 6.

Figure 10C,D shows the evolutionary and epidemiological effects of the fitness cost of deleterious mutations, sd. It is clear from these plots that neither the mean maximum nor the mean minimum tMRCA across the simulations depends strongly on sd (Figure 10C). However, the range of maximum tMRCA values is considerably higher at lower sd values (Figure 10C). This may be because selective sweeps are expected to occur more rarely at lower sd values (Figure 1E), such that the viral population is homogenized less frequently at these lower values, leading to higher tMRCA values. This explanation is consistent with the slightly lower annual attack rates at lower sd values (Figure 10D): when the rate of antigenic evolution is slower, individuals cannot be reinfected as rapidly and thus annual attack rates would be lower. Despite the dependency of evolutionary and epidemiological dynamics on sd, Figure 10C,D shows that a broad range of sd values yields dynamics that are consistent with influenza A/H3N2 dynamics in humans.

Finally, Figure 10E,F shows the sensitivity of the model to the antigenic mutation rate λantigenic. In the absence of antigenic evolution (λantigenic = 0), maximum tMRCAs are significantly higher than empirically documented (Figure 10E) and maximum annual attack rates do not exceed 3.8% (Figure 10F). These findings are consistent with, Figure 9, which indicates that a spindly viral phylogeny cannot be reproduced under purifying selection alone in a model that is parameterized for influenza. (Figure 9's model differs slightly from the model parameterized with λantigenic = 0, whose simulations are shown in Figure 10E,F: Figure 9 shows simulation results from an SIRS model that yields annual attack rates consistent with those empirically observed for flu; the results shown in Figure 10E,F under λantigenic = 0 use the same parameterization as in Figure 6, with the exception of λantigenic, and thus simulate a simple SIR model. In either case, purifying selection alone cannot reproduce a spindly phylogeny.) At increasing λantigenic values, tMRCAs decrease and annual attack rates increase to empirically documented values. These plots further indicate that influenza's evolutionary and epidemiological dynamics can be reproduced over a wide range of antigenic mutation rates as long as the rate exceeds a certain minimum value. Finally, taken together, Figure 10A,B,E,F again indicates that it is the interaction between deleterious and advantageous immune escape mutations that consistently yields a spindly phylogeny and high annual attack rates. Neither deleterious mutations nor immune escape mutations alone succeed in reproducing these dynamic features of influenza.

Discussion

Here, we have shown that population genetic and population dynamic models incorporating sublethal deleterious mutations can reproduce the characteristic features of influenza A/H3N2's evolutionary dynamics in humans. These include the virus's rare punctuated antigenic evolution and the low genetic diversity of its hemagglutinin protein. The low genetic diversity of the virus's HA, reflected in the spindliness of its phylogeny, has been a particular evolutionary characteristic that influenza modelers have sought to reproduce. To date, three other models exist that can explain this evolutionary characteristic of influenza (Ferguson et al., 2003; Koelle et al., 2006; Bedford et al., 2012). All three of these models, however, are subject to criticism. Influenza's epochal evolution model (Koelle et al., 2006) assumes that neutral or nearly neutral mutations accumulate at HA epitopes and that these changes enable a previously neutral mutation to exact a large antigenic effect and thereby to precipitate an antigenic cluster transition. Criticisms of this model are several. First, recent work indicates that amino acid changes that are responsible for cluster transitions have large antigenic effects in consensus sequences (Koel et al., 2013), such that genetic context is unlikely to be of utmost importance in determining the antigenic effect of mutations. Second, a molecular evolutionary analysis following the publication of the epochal evolution model has indicated that positive selection acts not only between antigenic clusters but also within antigenic clusters (Suzuki, 2008). In support of this finding, a more recent statistical analysis has shown that multiple antigenic variants co-circulate (Strelkowa and Lässig, 2012). Both of these empirical analyses are inconsistent with the assumptions of the epochal evolution model and indicate that the evolution of influenza is unlikely to be limited by the occurrence of antigenic mutations.

The two other existing models that can reproduce influenza's spindly HA phylogeny are the generalized cross-immunity model put forward by Ferguson et al. (2003) and the canalization model put forward by Bedford et al. (2012). When simulated, both of these models yield antigenic variant co-circulation and are thus consistent with the analyses detailed above that have shown that influenza evolution is not antigenic mutation limited. In Ferguson et al. (2003), generalized (strain-transcending) cross-immunity lasting on the order of 6 months was invoked to limit the genetic and antigenic diversity of influenza A/H3N2 in humans. The major criticism of this model is lack of empirical support for this duration and form of generalized cross-immunity: while there is evidence for long-lasting cross-immunity between more genetically distant strains (including heterologous influenza A subtypes), it appears to reduce pathology and possibly accelerate viral clearance rather than prevent infection (Grebe et al., 2008). A recent experimental study in ferrets indicates that generalized cross-immunity that prevents infection appears to exist, but that it lasts for less than a week (Laurie et al., 2015), which is insufficiently long to limit genetic and antigenic diversity in the model of Ferguson et al. The canalization model by Bedford et al. does not invoke generalized cross-immunity but instead assumes that antigenic mutations move viral strains in a two-dimensional (or higher) antigenic space. While these antigenic spaces or maps have been used to visualize the trajectories of flu's antigenic evolution (Smith et al., 2004; de Jong et al., 2007), models that start with the assumption of these maps considerably inflate the degree of competition between antigenic variants. This inflation of competition results from antigenic distances between daughter variants (variants that are produced from the same parent) necessarily being subadditive in this space. Inflated competition for susceptible hosts between antigenic variants is expected to lead to a more spindly phylogeny due to increased ‘niche overlap’ and therefore more frequent occurrences of between-strain competitive exclusion.

In light of our findings presented here and existing knowledge on the theory of asexual evolution, we can reflect on why the deleterious mutation and canalization models can reproduce limited genetic and antigenic diversity in simulations of influenza A/H3N2 in humans, whereas Ferguson et al. (2003) initially found that strain-specific immunity did not suffice in reproducing these patterns. From the population genetics literature on asexual populations in which interference effects are at play, it is well known that an increase in the fitness variance of a population acts to slow adaptive evolution and make the characteristic size of adaptive mutations that fix larger. While here we have invoked deleterious mutations in the generation of fitness variation, beneficial mutations can similarly contribute to inherited variation in fitness (Birky and Walsh, 1988; Barton, 1995). Thus, in population genetic models, interference between beneficial mutations (or a combination of beneficial and deleterious mutations) similarly slows down the rate of adaptive evolution and increases the size of adaptive mutants that fix (Birky and Walsh, 1988; Barton, 1995; Gerrish and Lenski, 1998; Rouzine et al., 2003, 2008; Park et al., 2010; Sniegowski and Gerrish, 2010; Schiffels et al., 2011; Good et al., 2012). However, all of these population genetic models assume a constant population size and therefore full resource competition between individuals. In the context of influenza, however, it is a change in antigenicity that is thought to provide the selective advantage. Because antigenic changes allow for escape from herd immunity, these changes by definition reduce competition for susceptible hosts (the virus resource). As such, antigenic mutants create a partially new niche and so do not necessarily lead to competitive exclusion. The establishment of antigenic mutants can therefore instead lead to long-term coexistence and, as a result of only partial cross-immunity, a larger infected population size. This is why explosive genetic and antigenic diversity is expected in the presence of only strain-specific immunity (Ferguson et al., 2003). In this context, generalized cross-immunity acts to considerably increase the competition between strains for susceptible hosts (as well as to reduce the overall infected population size). As such, the fitness variance generated by beneficial antigenic mutations in this model results in a decrease in the rate of antigenic evolution, an increase in the size of antigenic mutations that establish, and limited diversity in the long run, as in population genetic models that assume full competition between beneficial mutations by considering populations of constant size (Good et al., 2012). Similarly, the canalization model by Bedford et al. likely reproduces punctuated antigenic evolution and long-term limited genetic and antigenic diversity as a result of interference effects generated by inflated competition between antigenic mutations.

The generalized cross-immunity model, the canalization model, and the deleterious mutations model presented here therefore share fundamental similarities: they can reproduce the characteristic features of influenza evolution in humans by generating enough fitness variation among competing strains in the viral population. However, the first two of these models create fitness variation by beneficial antigenic mutations alone; because these mutations obtain their selective advantages by reducing competition for susceptible hosts, these models need other components (generalized cross-immunity or mutations that are necessarily subadditive in effect) to augment strain competition. In contrast, the deleterious mutation model we present here does not need to invoke processes to augment competition between antigenic strains because a large proportion of fitness variation in the virus population arises from differences in deleterious mutation loads (Figure 7A). Letting fitness variance be generated by deleterious mutations allows for limited diversity in the long run despite the co-circulation of antigenically very distinct variants that do not necessarily compete strongly for susceptible hosts. As such, our model can reproduce empirically observed high annual attack rates on the order of 10–15% (Figure 10B for λ = 0.10). Because deleterious mutations are known to circulate in the influenza A/H3N2 virus population (Fitch et al., 1997; Pybus et al., 2007; Strelkowa and Lässig, 2012), and in light of existing criticisms of the other two models, we therefore argue that the model presented here provides a more plausible mechanistic explanation for influenza's characteristic evolutionary features.

Our finding that specifically non-antigenic fitness variation is an important contributing driver in shaping the characteristic features of influenza's evolutionary dynamics in humans sheds light on other recent virological findings. One such finding is that cellular receptor binding avidity is an important phenotype that impacts viral fitness (Hensley et al., 2009). In a naive host, influenza viruses with low receptor binding avidities have a selective advantage, whereas, in a more immune host, influenza viruses with high receptor binding avidities have a selective advantage (Hensley et al., 2009). Which receptor binding avidity phenotype can be considered the ‘deleterious’ variant is therefore subject to which individual the virus finds itself in, as well as the overall degree of herd immunity in the host population (Yuan and Koelle, 2013). Because changes in receptor binding avidity frequently also alter antigenicity (Hensley et al., 2009), it is therefore also easily conceivable that certain of these changes increase virus fitness via two distinct mechanisms. As such, even though the genetic change is one and the same, the change in binding avidity that results can be thought of as increasing the background fitness of a virus strain, whereas the change in antigenicity that results can be considered as in this paper. Successful cluster transitions via the ‘jackpot’ strategy, as depicted in Figure 2C, may therefore preferentially involve mutations that simultaneously affect binding avidity and antigenicity. Indeed, the majority of cluster-transition mutations that have been recently characterized fall near the receptor binding site of influenza's HA (Koel et al., 2013), suggesting that changes in the receptor binding avidity phenotype that occur with these antigenic mutations may improve virus fitness.

In addition to cellular receptor binding avidity, glycosylation of influenza's HA is known to impact virus fitness. Whether glycosylation sites accumulate over evolutionary time (as in the case of H3N2 in humans [Blackburne et al., 2008]) or do not (as in the case of H1N1 in humans [Das et al., 2011]) therefore will depend on how these sites impact the background fitness of virus strains. In the case of H1N1, for example, glycosylation has been shown to significantly reduce receptor binding avidity and therewith to lower overall virus fitness, despite the beneficial effect of glycosylation on escape from antibody-mediated neutralization (Das et al., 2011). Compensatory mutations are therefore needed to restore virus fitness following the addition of a glycosylation site (Das et al., 2011). Other phenotypes that impact virus fitness include those that influence protein stability (Bloom and Glassman, 2009) and those that confer resistance to antivirals (Herlocher et al., 2004). In the case of permissive mutations that influence protein stability (Bloom et al., 2010; Gong et al., 2013), these mutations may not directly influence virus fitness; their occurrence, however, may impact the viability of other mutations that have fitness consequences and thus may alter the size distribution of viable beneficial mutations, including those that allow for immune escape. Epistatic interactions such as these can therefore be accommodated within this general framework of fitness variation generated by phenotypes other than antigenicity in driving patterns of influenza's antigenic evolution. While we here simply model this fitness variation as arising from circulating deleterious mutations, any of these non-antigenic phenotypes can similarly contribute to this fitness variation. Indeed, deleterious mutations are necessarily deleterious as a consequence of some phenotype, whether it is susceptibility to antivirals, a suboptimal receptor binding avidity, a protein with reduced stability, or simply another phenotype that reduces viral replication. The complementary phenotypes to these can conversely be considered as beneficial mutations that contribute to fitness variation. As such, there is not always a clear source of fitness variation in the influenza virus population. What is critical, however, is that a significant proportion of the virus's fitness variation has to arise from non-antigenic phenotypes that do not reduce competition between virus strains, as antigenic mutants alone result in explosive genetic and antigenic diversity as a result of effective niche partitioning. Our choice of modeling simply deleterious mutation accumulation, instead of specific non-antigenic phenotypes, stems from the finding that virus populations, and more specifically influenza virus populations, carry substantial deleterious mutation loads (Fitch et al., 1997; Pybus et al., 2007).

Returning to our model, given that an antigenic mutation arises in a strain carrying some deleterious mutation load, one might expect the virus population's deleterious mutation load to increase in the long run, causing a long-term decline in influenza A/H3N2 fitness. Two processes exist, however, that can keep this long-term accumulation of deleterious mutations from occurring. First, as the virus population becomes less fit the proportion of non-antigenic mutations that are beneficial is likely to increase, such that the mutation-selection balance becomes an evolutionary attractor (Goyal et al., 2012). Second, within-subtype viral reassortment could occur sufficiently frequently to keep any long-term decline in viral fitness at bay by combining segments with low mutational load onto the same genetic background. The possibility that influenza A/H3N2 has been declining in fitness over time may, however, also be entertained: a recent virological analysis of human H3N2 viruses points towards a long-term decrease (since 1968) in the propensity of the virus to bind human sialic acid receptors (Lin et al., 2012). This decrease has been invoked to explain the reduction in this virus's disease impact over the last 10 years (Lin et al., 2012). Whether this finding can be interpreted as evidence for a ‘weakening’ of the virus over time is unclear, especially because ‘weakening’ by any epidemiological measure has not been empirically demonstrated. Clearly, more virological studies are needed to determine the fitness trajectory of influenza A/H3N2 in humans over the past decades.

While we focused on the role of deleterious mutations in shaping influenza A/H3N2's antigenic evolution, the presence of circulating deleterious mutations should also impact the adaptive evolutionary dynamics of other flu types/subtypes in humans. For example, influenza B is known to have a lower mutation rate than influenza A/H3N2 (Nobusawa and Sato, 2006; Sanjuán et al., 2010). Given that the deleterious mutation rate λ is likely to decrease with the overall mutation rate, we would expect influenza B to carry a lower mean mutation load and further to have lower fitness variance. As such, we would expect smaller antigenic mutants to be able to successfully establish (Figure 1D). Our model would therefore predict that influenza B evolves antigenically in a less punctuated manner than influenza A/H3N2—a pattern that has been recently documented (Bedford et al., 2014). All else equal, we would also expect influenza B to have a faster rate of antigenic evolution (Figure 1E). However, a lower mutation rate would also surely reduce the rate at which antigenic mutations occur. The slower rate of antigenic evolution observed for influenza B (Bedford et al., 2014) is therefore not inconsistent with our model. Relative to these two influenza viruses, human influenza A/H1N1 shows a similar pattern to H3N2 in terms of punctuated antigenic evolution (Bedford et al., 2014). Despite this similarity, H1N1's rate of antigenic evolution is considerably slower than H3N2's, although it is faster than influenza B's rate of antigenic evolution (Bedford et al., 2014, 2015). The observed difference in rate of antigenic evolution between H3N2 and H1N1 is unlikely to reflect differences in these virus's mutation rates, since both are influenza A subtypes. Instead, these differences may stem from differences in their basic reproductive rates and, therefore, differences in selection pressures. This explanation is consistent with the finding that H1N1 appears to experience weaker antigenic selection than H3N2 in humans (Bhatt et al., 2011). Alternatively, the lower rate of antigenic evolution in H1N1 relative to H3N2 may be a consequence of differences in these viruses' global circulation patterns, which have only recently been described (Bedford et al., 2015).

Differences in the evolutionary dynamics of influenza viruses also exist across host species. For example, the same H3N2 virus that is circulating in humans emerged in pigs in the early 1970s. This swine influenza A/H3N2 evolves genetically at a rate that is similar to that of human influenza A/H3N2 (de Jong et al., 2007). However, its rate of antigenic evolution is approximately six times slower than the same virus's in humans (de Jong et al., 2007). This difference in the rate of antigenic evolution likely stems more from the stark ecological differences between the two hosts, rather than from differences in their deleterious mutation loads, with escape from humoral immunity being a less important evolutionary driver of HA in short-lived hosts than in humans.

Beyond the influenza viruses, circulation of deleterious mutations has been established in a wide range of RNA viruses (Pybus et al., 2007). The evolutionary dynamics of many of these viruses are characterized by spindly phylogenies and punctuated phenotypic changes. Whether deleterious mutations can account for these apparently similar patterns remains an open question. One especially intriguing case is that of norovirus, for which punctuated antigenic evolution has been documented (Lindesmith et al., 2008, 2011) and for which deleterious mutations along with differential binding to histoblood group antigens may contribute to fitness variation (Donaldson et al., 2008; Lindesmith et al., 2008). Similar to the case of influenza, an interplay between antigenic and non-antigenic/deleterious mutations may alone be sufficient to explain this viral evolutionary pattern, obviating the need to invoke the mutation-limited process of epochal evolution (Siebenga et al., 2007; Lindesmith et al., 2008).

In addition to helping us understand patterns of viral adaptive evolution, acknowledging the ‘rubbish around the ruby’—to paraphrase (Peck, 1994)—may also help us predict the course of adaptive evolution. Indeed, two recent publications have made great strides in predicting the genetic evolution of influenza A/H3N2 in humans over the short term (Neher et al., 2014; Łuksza and Lässig, 2014). Łuksza and Lässig's approach incorporated knowledge of antigenic and non-antigenic sites in developing a fitness model for this virus. With the assumption that mutations at non-antigenic sites were weakly deleterious, the authors were able to predict which influenza clades would grow and which ones would decline with an accuracy of 93% and 76%, respectively (Koelle and Rasmussen, 2014; Łuksza and Lässig, 2014). Instead of using knowledge specific to influenza A/H3N2, Neher et al.’s approach relied on branching patterns present in the HA phylogeny to predict strain evolution, with a success rate comparable to that of Łuksza and Lässig. The ability of both of these models to predict influenza evolution rests on the presence of substantial fitness variation in the influenza A/H3N2 virus population. The work here contributes to this understanding by reconciling the presence of fitness variation in the virus population with the observation of limited genetic and antigenic diversity of influenza A/H3N2's HA in the long run: cluster transitions will only succeed when large antigenic mutations find themselves in good genetic backgrounds, which must be largely determined by non-antigenic phenotypes. Successfully predicting cluster transitions will therefore require not only better characterizing the antigenic effects of mutations but also characterizing relative virus fitness, as mediated by differential deleterious mutation loads and contributions of non-antigenic phenotypes, in influenza's HA and other gene segments. An integrative understanding of these non-antigenic components of viral fitness will require the continued work of virologists and modelers alike, and ideally their interaction, for predicting viral evolution.

Materials and methods

The deleterious mutation-selection balance in an epidemiological context

The deleterious mutation-selection balance was classically derived in the context of a discrete generation, constant-size Wright–Fisher population (Haigh, 1978). We here briefly re-derive the deleterious mutation-selection balance for a virus population in an explicit epidemiological context. To do so we consider a virus population subject exclusively to deleterious mutations. Because we are not considering antigenic variation at this point, the epidemiological dynamics are governed by a basic SIR model. In this model, the number of susceptible hosts increases only through births into the host population and decreases through background mortality and infection of susceptible hosts. The number of infected hosts increases through infection of susceptible hosts and decreases through recovery from infection and through natural mortality of infected hosts. The number of recovered hosts increases through recovery of infected individuals and decreases through natural mortality of recovered hosts. To extend this basic SIR model to allow for a virus population undergoing deleterious mutations, we classify infected individuals according to the number of deleterious mutations harbored by the virus they carry. This classification makes the implicit assumption that an infected host harbors a genetically homogeneous virus population; that is, in this case, that there is no variation in the number of deleterious mutations carried by distinct virions that make up the intrahost viral population. Although the assumption of a genetically homogeneous within-host virus population is clearly a simplifying assumption, it is an assumption that is commonly made in population-level models of influenza evolution (Andreasen et al., 1997; Gog and Grenfell, 2002; Ferguson et al., 2003; Koelle et al., 2006; Bedford et al., 2012). As in Haigh (1978), we assume that mutations occur at birth, which, in the context of a virus population, are disease transmission events. We choose this approach to introduce new deleterious mutations over an approach that assumes that infected individuals can ‘mutate’ from being infected with a virus having a specified number of deleterious mutations to being infected with a virus having a higher number of deleterious mutations. This is because it is unlikely that a deleterious mutation that arises within a host could rapidly sweep to fixation in that host once the intrahost viral population is large. Introducing deleterious mutations at disease transmission events is therefore a more biologically reasonable assumption. Again as in Haigh (1978), we further let the number of deleterious mutations that arise at ‘birth’ be Poisson distributed with mean λ, where λ is the per-genome per-transmission deleterious mutation rate. With these assumptions, the rate of change in the number of infected individuals carrying a virus with k deleterious mutations is given by:

(1)dIkdt=SN∑j=0k(βk−je−λλjj!Ik−j)−(μ+γ)Ik,

where the first term in this equation captures the increase in the number of individuals infected with virus in mutation class k arising from the transmission process and the second term captures the decrease in the number of infected individuals resulting from background mortality (at per capita rate μ) and recovery from infection (at per capita rate γ). In the first term, S is the number of susceptible hosts, N is the (constant) host population size, and βi is the transmission rate of a virus carrying i deleterious mutations. The Poisson term e−λλjj! provides the probability that j deleterious mutations occur at transmission. As in Haigh (1978), we assume multiplicative fitness effects of deleterious mutations, with each deleterious mutation exacting a fitness cost of size sd. We can thus write βi in terms of the transmission rate of a virus in the highest fitness class (i.e., the lowest deleterious mutation class):

(2)βi=β0(1−sd)i.

The dynamics of susceptible hosts are given by:

(3)dSdt=μN−μS−SN∑k=0∞βkIk,

where the first term captures births into the host population (at a per capita rate μ, equal to the background mortality rate), the second term captures background mortality of susceptible hosts, and the third term captures depletion of susceptible hosts through the transmission process. The dynamics of recovered individuals are simply given by dRdt=γ∑k=0∞Ik−μR, where the first term captures the increase in the number of recovered individuals through recovery of infected individuals and the second term captures background mortality. For this set of equations, we can define the basic reproductive rate R0,k as the expected number of secondary infections (belonging to any mutation class) produced by a single infected individual carrying a virus with k deleterious mutations in a completely susceptible population. This mutation class-specific basic reproductive rate is given by R0,k=βk(μ+γ).

To solve for epidemiological equilibrium, we can now consider the dynamics of the first infected mutation class k = 0: dI0dt=β0e−λSNI0−(μ+γ)I0. Setting this equation to zero, the equilibrium fraction of the population susceptible to infection can be solved for:

(4)S^N=(μ+γ)β0e−λ,

S^N is greater than 1R0,k=0 by a factor of 1e−λ. That the fraction of susceptible hosts exceeds the inverse of the basic reproductive number of the highest-fitness virus class makes sense as the highest-fitness virus class occupies only a fraction of the total virus population.

The equilibrium number of infected individuals carrying virus with k deleterious mutations, I^k, can now be solved. Substituting Equations 2, 4 into Equation 1 and simplifying yields:

dIkdt=(μ+γ)Itot∑j=0k((1−sd)k−jλjj!pk−j)−(μ+γ)Ik,

where the total number of infected individuals is given by Itot=∑i=0∞Ii and pi is the proportion of infected individuals carrying virus with i deleterious mutations, pi=IiItot. Setting this equation to 0 and solving for pk yields:

where θ = λ/sd. The total number of infected individuals at equilibrium, I^tot, can now be solved by substituting Equations 2, 4 into Equation 3, and replacing Ik with pkItot. Setting equal to 0 and solving for Itot yields:

I^tot=e−λμN(1−μ+γβ0e−λ)(μ+γ)∑k=0∞(1−sd)kpk.

Substituting Equation 5 into the above expression and simplifying yields:

(6)I^tot=μN(μ+γ)(1−μ+γβ0e−λ).

The equilibrium number of mutation class-specific infected individuals can then be calculated using I^k=pkI^tot for any deleterious mutation class k.

Absolute fitness in epidemiological models is provided by the net reproductive rate R. With a population experiencing deleterious mutations, each mutation class k will have its own net reproductive rate Rk at equilibrium. Given the above definition of the basic reproductive rate for viruses carrying k deleterious mutations, the fitness cost associated with deleterious mutations (Equation 2), and the equilibrium fraction of susceptible hosts (Equation 4), the mutation class k net reproductive rate is given by:

(7)Rk=R0,k(S^N)=eλ(1−sd)k.

As expected, the mean net reproductive rate of the virus population (∑k=0∞(Rkpk)) is 1 when the population is at epidemiological equilibrium, with some viruses having a net reproductive rate above 1 and other viruses having a net reproductive rate below 1. The variance in the net reproductive rate is given by ∑k=0∞(Rk2pk)−1, which increases with increases in λ and increases with decreases in sd.

Initial and final reproductive rates of an antigenic mutant evolving de novo from a resident viral population

To compute the antigenic mutant's initial and final net reproductive rates, we first use an epidemiological model to mathematically describe the interaction between the resident antigenic strain and the antigenic mutant. We denote the resident strain with super- and subscripts r and the antigenic mutant with super- and subscripts m, and use a history-based model (Andreasen et al., 1997) to specify the immunological interaction between these two strains. Note here that we again make the implicit assumption that the intrahost viral population is genetically homogeneous with respect to both antigenicity and the number of deleterious mutations carried. With deleterious mutations accumulating at transmission, we have:

dS0dt=μN−μS0−S0N∑k=0∞(βk(I0,kr+IB,kr+I0,km+IA,km)),

dI0,krdt=S0N∑j=0k(βk−je−λλjj!(I0,k−jr+Im,k−jr))−(μ+γ)I0,kr,

dI0,kmdt=S0N∑j=0k(βk−je−λλjj!(I0,k−jm+Ir,k−jm))−(μ+γ)I0,km,

(8)dSrdt=γ∑k=0∞I0,kr−μSr−σSrN∑k=0∞(βk(I0,km+Ir,km)),

dSmdt=γ∑k=0∞I0,km−μSm−σSmN∑k=0∞(βk(I0,kr+Im,kr)),

dIm,krdt=σSmN∑j=0k(βk−je−λλjj!(I0,k−jr+Im,k−jr))−(μ+γ)Im,kr,

dIr,kmdt=σSrN∑j=0k(βk−je−λλjj!(I0,k−jm+Ir,k−jm))−(μ+γ)Ir,km,

dS{r,m}dt=γ∑k=0∞(Ir,km+Im,kr)−μSr,m,

where Sx is the number of uninfected individuals who have previously been infected with strain(s) x, and Ix,ky is the number of individuals previously infected with strain x who are currently infected with a virus of strain y carrying k deleterious mutations. The parameter σ quantifies the extent to which susceptibility to infection with one strain is affected if the individual has previously been infected with the other strain. With σ = 1 a previous infection does not reduce susceptibility to infection with the second strain, whereas with σ = 0 a previous infection results in complete protection from reinfection with a second infection. A higher value of σ therefore corresponds to a greater degree of immune escape (i.e., a larger antigenic change). As expected, at the time immediately prior to the antigenic mutant's emergence, Equation 8 reduce to Equations 1, 3.

Given Equation 8, an antigenic mutant that arises in a background with i deleterious mutations initially has a net reproductive rate of:

(9)Rm(t=0)=(βiμ+γ)(S0N+σSrN),

where we have defined time t in terms of the time since the antigenic mutant's emergence.

With S^0N=(μ+γ)β0e−λ (Equation 4) and with S^rN≈1−S^0N, the mutant's initial net reproductive rate is given by:

(10)Rm(t=0)=eλ(1−sd)i(1+σ(β0e−λμ+γ−1)).

This expression would be exact if we allowed for coinfection, treating individuals currently infected with the resident strain similarly to individuals previously infected with the resident strain. Equation 10 consists of three components: (i) eλ is the net reproductive rate of resident strain viruses that are in mutational class k = 0; (ii) (1 − sd)i quantifies the extent to which the antigenic mutant's initial reproductive rate is reduced as a result of the deleterious mutations it carries; and (iii) the term (1+σ(β0e−λμ+γ−1)) quantifies the extent to which the antigenic mutant's initial reproductive rate is increased as a result of immune escape. To make the link stronger between this model and traditional population genetic models, we can define the selective advantage of an antigenic mutant at the time of its emergence as sb=σ(β0e−λμ+γ−1).

The reproductive rate of a strain carrying an antigenic mutation necessarily decreases as it establishes through its own accumulation of deleterious mutations. Neglecting any changes in the host immune landscape, the final mean reproductive rate of the antigenic mutant lineage is:

(11)Rm(t=∞)=eλ(∑j=i∞pj−i(1−sd)j)(1+σ(β0e−λμ+γ−1)),

where the second term of the product quantifies the extent to which the antigenic mutant lineage's reproductive rate is reduced as a result of the deleterious mutations it carries once it has reached its own deleterious mutation-selection balance. Because this term ignores the possibility of stochastic loss of mutation class i during the strain's establishment, Equation 11 is an upper estimate for Rm(t = ∞). Changes in the host immune landscape would change the third term of the product.

Under this simple model, successful establishment can only occur when Rm(t = ∞) > 1. From Equation 11, it is therefore clear that the probability of an antigenic mutant's establishment is higher the lower its original mutation load i. This is because, in the absence of any compensatory or back-mutations, i provides the lower limit to the deleterious mutation load carried by the new antigenic strain. It is also clear from this equation that the probability of an antigenic mutant's establishment is higher the greater the ability of the mutant strain to reinfect previously infected individuals, reflected in a higher σ.

Model parameters

Two evolutionary parameters determine influenza's deleterious mutation load: the per-genome per-transmission deleterious mutation rate λ and the transmission fitness cost of a deleterious mutation sd. Neither of these parameters have been estimated specifically for influenza. We therefore do our best with estimating them using existing data from other viruses or via rough estimates that incorporate existing knowledge about influenza. The parameter λ was roughly computed by first multiplying the genome size of influenza's major coding regions (12,741 nucleotides [Holmes et al., 2005]) with the empirical estimate of 2.3 × 10−5 for the number of substitutions that occur per nucleotide per cell infection in influenza A viruses (Sanjuán et al., 2010). This yielded a per-genome, per-transmission substitution rate of 0.293. With roughly 2/3 of mutations being non-synonymous and roughly 50% of non-lethal, non-synonymous mutations being deleterious (Sanjuán et al., 2004), we calculated the per-genome, per-transmission deleterious mutation rate λ as the product: 0.293×23×0.5≈0.10. Given the uncertainty in λ, we look across a range of λ values from λ = 0 to λ = 0.20 in Figure 1D,E.

To get a rough estimate of the transmission fitness cost of deleterious mutations for influenza, we start with empirical estimates of mean fitness cost values for the five viruses studied in Sanjuán (2010): 0.103, 0.107, 0.112, 0.126, and 0.132. These fitness costs are in vitro fitness costs associated with viral growth in cells. These costs therefore differ from transmission fitness costs at the population level, which is how sd is defined in the models we present. To obtain a rough estimate for transmission fitness cost from these within-host viral growth fitness costs, we used a published within-host model of influenza dynamics (Baccam et al., 2006) to simulate viral load dynamics under two scenarios: a ‘wild-type’ parameterization and a ‘deleterious mutant’ parameterization, where the mutant carries a specified fitness cost. The ‘wild-type’ scenario used parameter values that were estimated in Baccam et al. (2006) using viral load data from six individuals who were experimentally infected with influenza. To parameterize the ‘deleterious mutant’ scenario, we assumed that the in vitro fitness cost manifested itself through a reduction in the within-host viral production rate. We therefore set all of the within-host parameters of the ‘deleterious mutant’ scenario to be equal to the ones used in the ‘wild-type’ scenario, with the exception of the viral production rate. This parameter we set to the product of the ‘wild-type’ scenario's viral production rate and the quantity 1 minus the in vitro fitness cost. For the six patients in Baccam et al. (2006), we simulated the ‘wild-type’ and the ‘deleterious mutant’ scenarios for each of the five in vitro fitness costs listed above. We then mapped these simulated viral load dynamics to between-host transmission rates by assuming that viral transmission rates are proportional to the log of viral load: β(τ)∝log10(V(τ)). This assumption is commonly used (see Handel et al., 2013) and has some empirical support, particularly from studies of HIV (Quinn et al., 2000; Fraser et al., 2007). The population-level transmission fitness costs sd for the five empirical within-host fitness cost estimates from Sanjuán (2010) could thus be inferred for each of the patients studied in Baccam et al. (2006). Our simulations indicate, first and foremost, that relatively large fitness costs within a host translate into much smaller fitness effects at the population level. This is in part unsurprising given that we assume that transmission rates scale with the log of the viral load. We further found that in vitro fitness costs in the last three patients studied in Baccam et al. (2006) would translate into fitness benefits in terms of transmission. Rather than trusting this to be the case biologically, this paradoxical result is likely to be a result of the rough mapping that we (and others) use for translating between viral load and transmissibility. We therefore restricted ourselves to the first three patients studied in Baccam et al. (2006). For them, we found that within-host fitness costs of 0.103–0.132 yielded transmission fitness costs between 0.002 and 0.018, with a mean transmission fitness cost sd of 0.008, which we use in Figure 1A–C. Given the uncertainties present in our estimate of sd, we instead show three different values for sd in Figure 1D,E. These indicate that the qualitative effects of circulating deleterious mutations (larger average antigenic sizes of mutants that establish and slower rates of antigenic evolution) are robust to specific values for sd.

In addition to these two evolutionary parameters, the models we present contain three epidemiological parameters: the birth/death rate μ, the recovery rate γ, and the basic reproductive rate of a virus in the k = 0 mutation class (R0,k = 0). We use a birth/death rate of μ = 1/30 years−1, based on crude birth rate estimates of approximately 33 per 1000 individuals per year (United Nations, 2011). This value for μ is also consistent with a recent flu modeling study (Bedford et al., 2012). The recovery rate γ was chosen based on an in-depth meta-analysis of 56 human challenge studies (Carrat and Flahault, 2007). In this study, the authors found an average duration of viral shedding of 4.8 days for influenza viruses (regardless of subtype). They also calculated an average generation time of 3.1 days for influenza A/H3N2, where generation time is defined as the time between an individual becoming infected and transmitting the virus. Because, in an epidemiological model with a constant transmission rate and a constant recovery rate, the duration of infectiousness (taken to be the duration of viral shedding) and the generation time take the same value when a virus is endemically circulating, we chose an intermediate value between these estimates, letting the recovery rate γ = 1/4 days−1. We chose an R0,k = 0 of 2.25. This value falls within the range of recent R0 estimates for influenza A/H3N2 (range = 1.21–3.58 for influenza A/H3N2's second pandemic wave [Jackson et al., 2010]).

Figure 1D,E require specification of a distribution for the degree of immune escape σ. We assume that σ is gamma distributed with mean of 0.012 and a shape parameter of 2. The choice of a gamma distribution is based on a virological study showing that the distribution of beneficial mutation effects appear to be gamma distributed (Sanjuán et al., 2004). Our choice of mean σ value is not based on an independent empirical estimate as relevant data to parameterize this value to our knowledge do not exist. Assuming an exponential distribution with a mean of 0.012 for σ instead of a gamma distribution did not change the qualitative findings that circulating deleterious mutations act to slow antigenic evolution and make it more punctuated in nature.

Figures 2, 3 in the main text show stochastic simulations of the epidemiological model given by Equation 8. To simulate the model stochastically, we used Gillespie's τ-leap method with a time step of 1 hr. These stochastic simulations required further specification of a host population size N. We used a population size of N = 4 billion hosts, corresponding to the human population size around 1980, which can also be roughly considered today's tropical population size.

Description and implementation of the phylodynamic model

To explore the full evolutionary dynamics of our model with non-antigenic and antigenic mutations, we implemented an individual-based model that allowed for an arbitrary number of new strains to enter the population and co-circulate. Individuals in this model were categorized as either currently infected or currently uninfected. Using the same notation as in history-based multi-strain models (Andreasen et al., 1997), each currently uninfected individual carries a list of antigenic types with which he has previously been infected. Each infected individual similarly carries this list of previously experienced antigenic types. In addition to this list, each infected individual carries a current viral infection that is characterized by the infecting virus's deleterious mutation load and its antigenic type. Again, we assume for the sake of model simplicity that the intrahost viral population is genetically homogeneous such that a single deleterious mutation load and a single antigenic type suffice in characterizing a viral infection. Upon recovery, infected individuals add to their strain history list the antigenic type of the viral infection from which they are recovering.

Viral antigenic types are defined by their relationship to one another in terms of their pairwise cross-immunity values σ. The minimum antigenic distance between the challenging strain and the host's repertoire of previously encountered strains (as provided by the host's antigenic type list) determines the probability of a currently uninfected host becoming infected by the challenging strain. Specifically, the probability of infection given contact with a host infected with strain i was given by pij = min(1.0, σij), where strain j is the strain in the host's repertoire that is antigenically closest to strain i. This formulation leads to complete immunity to reinfection with antigenic types that a host has previously experienced and complete susceptibility to infection for naive hosts. To compute any σij value, the model uses tracked parent–offspring relationships of distinct antigenic variants, with mother–daughter variants differing by one antigenic mutation having a degree of cross-immunity that is drawn from the specified antigenic size distribution (see below). The degree of cross-immunity σij between two strains i and j is assumed to be additive; for example, if strain i gives rise to strain l and strain l gives rise to strain j, then σij = σil + σlj.

When they occurred, both antigenic and non-antigenic mutations occurred during transmission events. A virus in an infected host therefore never changed antigenically or changed in mutation load during the course of infection. At a transmission event, the number of non-antigenic mutations was drawn from a Poisson distribution with mean λ. Each of the non-antigenic mutations was characterized as beneficial (with probability εc) or deleterious (with probability 1 − εc), where εc << 1. The net change in the number of deleterious mutations was then calculated and the virus's mutation load was appropriately increased or decreased. Our phylodynamic simulations required that a proportion of non-antigenic mutations were beneficial rather than deleterious because of the finite number of hosts we were computationally able to simulate. Specifically, in finite populations it is well known that the lowest-load mutation classes will at some point be stochastically lost, initiating Muller's ratchet (Muller, 1964). Recent work has indicated, however, that if a fraction of mutations are beneficial rather than deleterious, the average mutation load in a population can be maintained indefinitely over time (Goyal et al., 2012). In our simulations, we set εc to a value of 0.16, which resulted in a stable deleterious mutation load over time.

The number of antigenic mutations occurring at a transmission event was drawn from a Poisson distribution with mean λantigenic, independently of the number of non-antigenic mutations. We used a λantigenic of 0.00075. The size of each of these antigenic mutations was drawn from a gamma distribution with a mean of 0.012 and a shape parameter of 2. The antigenic difference σij between the virus i in the infecting host and the virus j in the host becoming infected was calculated as the sum of the sizes of the antigenic mutations occurring at transmission.

The simulations were run using an individual-based stochastic simulation algorithm based on Gillespie's tau-leap algorithm, using a time step τ of 1 day. Experimentation with smaller time steps yielded similar results, such that we used a τ of 1 day for computational efficiency. Over each time interval τ, the number of events occurring during that interval was drawn from a Poisson distribution. The modeled events were births, deaths, recoveries, and infectious contacts. Individuals born into the population were considered naive to infection (such that their strain history lists were empty). Deaths occurred from both uninfected and currently infected hosts, indiscriminately. For each infected host death, an infected individual was randomly chosen to be discarded from the population, independent of the virus currently infecting the individual. Similarly, for each recovery, an infected individual was chosen at random to recover. For each infectious contact from individuals infected with a class-k virus, an infected individual was chosen at random from the current class-k infected population and a currently uninfected individual was also chosen at random. The uninfected host was then infected with probability pij, where, as described above, pij is given by min(1.0, σij), where σij is the antigenic distance between infecting strain i and the strain j in the host's repertoire that is antigenically closest to strain i. We did not allow for coinfection, so our phylodynamic simulations did not include the possibility of viral reassortment leading to antigenic cluster transitions. For these individual-based simulations we used a population size of N = 40 million individuals, 1/100th of the human population size in 1980. Our choice of this limited population size was due to the large amount of memory required to keep track of the immune histories of each individual host. The relationships of who-infected-whom and at what time were tracked such that the true phylogenetic tree could be reconstructed and compared against influenza A/H3N2's HA phylogeny reconstructed from empirical sequences.

Preliminary simulations revealed highly volatile population dynamics, with large peaks in prevalence followed by extinction or long periods of low prevalence. To stabilize the dynamics we allowed a small number of infectious contacts to occur from outside of the focal population by having uninfected individuals experience an additional force of infection from an external pool of infected individuals (M = 200 for all simulations). So as not to alter the evolutionary dynamics, this external pool of infected individuals was assumed to have the same mutational load distribution and antigenic type frequencies as the focal population.

All remaining evolutionary and epidemiological parameters used in these individual-based simulations were the same as the ones listed in the model parameterization section above. The simulations were implemented in Java using a modified version of the program Antigen (http://bedford.io/projects/antigen/). The original code was modified to track the mutation load and the antigenic phenotypes of the viral population, which in our case requires us to track the pairwise distances between all antigenic types instead of viral locations on a two-dimensional antigenic surface. Source code for our full phylodynamic model is available on the GitHub repository at https://github.com/davidrasm/MutAntiGen.git.

In our simulations, we tracked the fitness of individual viruses in terms of their net reproductive rate R. For a particular virus i carrying k deleterious mutations,

R=βk(μ+ν)SeffN,

where Seff is the effective number of susceptible hosts available to a particular virus. Because the susceptibility of a host to a particular virus depends on the host's detailed immune history, to compute Seff we first compute the probability ρij of strain i infecting each host in the population upon contact, and then sum ρij over all uninfected hosts in the population to arrive at Seff. In practice, for computationally tractability, we only sum over a large number (n = 10,000) of randomly sampled hosts to approximate Seff. Given the distribution of R in the viral population, we can then compute the fraction of viral fitness variation attributable to deleterious mutations (i.e., variation in βk) vs antigenic differences among strains (i.e., variation in Seff). To decompose the total variance in fitness into these components, we first log-transform R, so that

ln(R)=ln(βk)+ln(SeffN)+ln(1(μ+ν)).

Dropping the constant 1(μ+ν) term, we can then decompose the total variance in log R into its component parts using the well-known relation that the variance of the sum of two random variables is equal to the sum of their variances plus their covariance,

var(R)=var(βk)+var(SeffN)+2cov(βk,SeffN).

Note that the covariance term cannot be ignored because βk and SeffN are not independent random variables and can be correlated due to the shared common ancestry of viral strains (i.e., phylogenetic correlations). To compute the fraction of total viral fitness variation attributable to βk and thus deleterious mutations, we subtract the variance attributable to the covariance:

var explained byβk=var(βk)var(R)−2cov(βk,SeffN).

We also conducted additional simulations that did not include any antigenic evolution but did include deleterious mutations to determine if the spindly phylogeny shown in Figure 6 arose simply due to the presence of purifying selection rather than being driven by antigenic change. Simulations without antigenic evolution were conducted with the same evolutionary and epidemiological parameters; average viral R0 was as in the simulations with antigenic evolution except λantigenic was set to zero. However, because there was no antigenic evolution, the average number of infected humans over time and the annual attack rates were significantly lower than in the simulations with antigenic evolution. To compensate for this, in Figure 9, we allowed immunity to wane over time, as in standard SIRS epidemiological models. The rate of immune waning was set such that immunity from a previous infection lasted 12 years on average before the host became completely susceptible again. This rate of immune waning was chosen to produce an average annual attack rate of ∼5–10%, consistent with what was observed in the simulations with antigenic and non-antigenic mutations.

Decision letter

Arne Traulsen

Reviewing Editor; MPI Ploen, Germany

eLife posts the editorial decision letter and author response on a selection of the published articles (subject to the approval of the authors). An edited version of the letter sent to the authors after peer review is shown, indicating the substantive concerns or comments; minor concerns are not usually shown. Reviewers have the opportunity to discuss the decision before the letter is sent (see review process). Similarly, the author response typically shows only responses to the major concerns raised by the reviewers.

Thank you for sending your work entitled “The effects of a deleterious mutation load on patterns of influenza A/H3N2's antigenic evolution in humans” for consideration at eLife. Your article has been favorably evaluated by Ian Baldwin (Senior Editor) and three reviewers, one of whom served as Guest Reviewing Editor.

The paper by Koelle and Rasmussen treats the effects of deleterious mutations in the evolution of the influenza virus. This is an important aspect of influenza evolution, because traditional theory and data analysis of influenza evolution has focused on antigenic effects alone and neglected other evolutionary forces such as the genetic load caused by deleterious mutations. The dynamics analysed in this paper are complex, because they are shaped by the interplay of beneficial antigenic mutations, deleterious background mutations, and by large fluctuations in population size and absolute fitness caused by the epidemiology of the virus. The main findings are that deleterious mutations (a) change the speed and character of antigenic evolution, and (b) affect the epidemiology of influenza; in particular, they can explain the linear shape of the influenza tree.

Explaining the elegant ladder-like (‘spindly’) phylogenetic pattern of H3N2 evolution in humans has been a prize aim of infectious disease modelers for at least a decade. The authors present a thought-provoking model that emphasizes the importance of considering key antigenic substitutions within the context of a high background deleterious mutation load. The need for a very fit mutation (that causes a large antigenic change) to overcome a deleterious background also explains the punctuated nature of antigenic changes.

Some parts of the article are quite difficult to follow, and it is not clear if it is possible to reproduce the results without further explanation. Please provide a more thorough description of the model.

Essential revisions:

1) Connection to the theory of asexual evolution:

The results in the first part of the paper should be put in context of the general theory of asexual evolution. In light of this theory, the result that deleterious mutations reduce the probability of fixation for beneficial mutations is well known. This effect is contained in the theory of background selection (e.g. in the work of B. Charlesworth), as well as in travelling wave models: deleterious mutations generate fitness variance, which makes the fixation probability of a beneficial mutation strongly background-dependent (Good et al., Distribution of fixed beneficial mutations and the rate of adaptation in asexual populations, PNAS 2012). Beneficial mutations with effect size below a certain threshold fix with a reduced rate close to neutral mutations, which results in an increased average effect of fixed mutations (Good et al., PNAS 2012; Schiffels et al., Emergent neutrality in adaptive asexual evolution, Genetics 2011); specifically for influenza, the reduction of the adaptive rate due to linked deleterious mutations in non-adaptive sequence has also been observed in Strelkowa and Lässig, 2012.

It would also be good to discuss the possible consequences of epistasis in this context (if applicable).

2) Stressing the interplay between evolution and epidemiology:

In view of 1), it would be good to focus the paper more on the interplay between the evolutionary and the epidemiological dynamics, which is the novel and most interesting part of the paper. In particular, the authors find that the reduced rate and increased average effect of antigenic mutations may contribute to the linear shape of the influenza tree and to increased attack rates. However, it should become clearer that a strong pruning of antigenic mutations by the joint dynamics with deleterious mutations requires that the latter contribute a substantial fraction of the average fitness variance in the population. This may be the case in model simulations (and the results on tree shape should be compared with the relative contributions of mutation classes to fitness variance). But it is not clear for the actual system. Furthermore, we suggest the authors juxtapose their finding of the influence of deleterious mutations on tree shape with previous explanations of this characteristic.

3) The effect of deleterious mutations:

In the Discussion, the authors argue that the presence of deleterious mutations reduce clonal competition compared to populations without mutational load. We think this argument is incorrect. According to all models cited above, deleterious mutations should add to the fitness variance in the population and, hence, increase the amount of interference selection. This does not contradict the suggested pruning of antigenic mutations, because in the joint model, antigenic competition is no longer equivalent to fitness competition. The antigenic variance in the population may decrease with increasing rate or effect of deleterious mutations, even if the total fitness variance increases.

4) Connection to empirical observations and previous models:

One concern is that the model is not particularly data-driven. Although empirical trees are presented in Figure 4 to illustrate 3 pathways to antigenic change that are consistent with the model, these trees are more useful as illustrative of concepts than as good evidence for the model. The authors should stress the illustrative character of that figure.

As for (virtually) all theoretical models, it is possible that the patterns in the trees could alternatively be explained by other mechanisms. For example, in the BE92 to WU95 cluster transition, the 145N to 145K substitutions that did not take off globally (that are dead-ends) could have occurred in locations that are thought to not be global source regions (i.e. could the failure of those viruses to persist globally despite beneficial mutations relate to their emergence in less inter-connected non-source regions (i.e., other than SE Asia))? A more thorough discussion of the limits of the model would be in place. For example, in the WU95-SY97 transition, would a scenario not predicted by the model be a longer branch length? How long?

It's useful to explore how robust a flu model is to other subtypes and hosts, and the authors make an admirable attempt to consider their model in the context of influenza B and swine. However, the most natural comparison would be with seasonal H1N1, and its omission is striking here. The authors also mistakenly assume that H3N2 has similar mutation rates in different hosts, and as a result miss the opportunity to examine how well their model works in a system such as swine (a particularly good example because versions of the human H3N2 virus also circulates in pigs) which have evolutionary rates that are higher than in humans for non-surface proteins (Worobey paper) and should carry a higher mutational load than humans (Discussion, third paragraph). The swine comparison raises ecological considerations that may also relate to human demography.

In the Discussion, it would be important to describe in more detail if the competing models have fundamental issues in capturing the current knowledge or if they would be valid for other parameter choices. It seems odd to write “2%” (Discussion, first paragraph) for a model that is not described in detail and that has presumably free parameters.

5) Technical issues:

In general, it is difficult to follow the mathematical description of the model. The description is not particularly complete. For example, in Equation 3: Why does the general background death rate lead to a positive term for the change in the number of susceptible hosts? Or is it only the background death rate of infected individuals? Why does recovery not lead to a positive term here? Another example is the definition of R0,k=0, which differs from the inverse of Equation 5 only by a factor in the exponential function, or Equation 13, which is not a mere modification of Equation 6, but an extension.

The statement in the last paragraph of the subsection “Model parameters” is problematic. Equation 14 is a set of (deterministic) ODEs and thus does not give a stochastic model. It is well known that many stochastic models can be the basis of a single ODE model, so there is no unique model that the authors simulate. Are the rates in Equation 14 directly used, i.e. for a Gillespie algorithm? Or do they reflect net rates that describe the net effect of several rate processes?

[Editors' note: further revisions were requested prior to acceptance, as described below.]

Thank you for resubmitting your work entitled “The effects of a deleterious mutation load on patterns of influenza A/H3N2's antigenic evolution in humans” for further consideration at eLife. Your revised article has been favorably evaluated by Ian Baldwin (Senior Editor), a Reviewing editor, and two reviewers. The manuscript has been improved but there are some remaining issues that need to be addressed before acceptance, as outlined below:

In the revised version, the authors have successfully addressed most of the comments. In particular, they have substantially improved the population-genetic side of their findings, which was one of the main concerns of the reviewers.

However, we ask the authors to be more explicit about some of the assumptions in their model. The success or failure of a beneficial antigenic mutation is a complex process, and while the authors do not need to incorporate all aspects into their model, they should be explicit about what is not being factored in:

1) The authors refer to selection as occurring on a ‘strain’, whereas in reality mutations are occurring within a swarm (quasi species) of intrahost diversity. The frequency of mutation as well as back-mutation within this swarm means that beneficial and deleterious mutations are not only arising but also being removed. The impact of a mutation on viral fitness is therefore also a product of the frequency of that mutation within the swarm. The potential bottleneck at transmission also means that not all variants, particularly low-frequency variants, will be transmitted between hosts.

2) We believe the authors are also assuming in this model a consistent immune landscape, for simplicity, which does not reflect variation in immune responses and prior exposures among hosts that affect the fitness of antigenic mutations.

Host variation can have a profound effect on the fitness of a particular mutation (as evidenced nicely by the antigenic mutations in receptor binding sites, which have very different fitnesses in hosts with different immune profiles based on whether antigenic escape or binding avidity is more important.

3) Successful H3N2 antigenic variants most frequently are produced in Southeast Asia before spreading globally. The spatial ecology of the virus has pronounced effects on which mutations succeed on a global scale and which are not sustained. It would have been simple to examine the proportion of early 145K mutations that arose in SE Asia on the phylogeny (suggested in the original reviewer comments), but short of that the importance of spatial ecology in the global success of a particular variant should be discussed (no additional analysis is required).

4) The authors should clarify that/why beneficial non-antigenic mutations were not included in the model.

5) There is still a lack of clarity and the role of competition between antigenic variants in the model. We find it hard to agree with the argument that 'rather, we propose that many of these circulating antigenic variants ultimately decline from the accumulation of deleterious mutations in the context of an only slowly changing herd immunity landscape”.

Author response

Essential revisions:

1) Connection to the theory of asexual evolution:

The results in the first part of the paper should be put in context of the general theory of asexual evolution. In light of this theory, the result that deleterious mutations reduce the probability of fixation for beneficial mutations is well known. This effect is contained in the theory of background selection (e.g. in the work of B. Charlesworth), as well as in travelling wave models: deleterious mutations generate fitness variance, which makes the fixation probability of a beneficial mutation strongly background-dependent (Good et al., Distribution of fixed beneficial mutations and the rate of adaptation in asexual populations, PNAS 2012). Beneficial mutations with effect size below a certain threshold fix with a reduced rate close to neutral mutations, which results in an increased average effect of fixed mutations (Good et al., PNAS 2012; Schiffels et al., Emergent neutrality in adaptive asexual evolution, Genetics 2011); specifically for influenza, the reduction of the adaptive rate due to linked deleterious mutations in non-adaptive sequence has also been observed inStrelkowa and Lässig, 2012.

It would also be good to discuss the possible consequences of epistasis in this context (if applicable).

We agree that situating the results of the first part of the paper in the context of the theory of asexual evolution would greatly improve the paper. To this end, we have first added the following text to the Introduction: “We start by extending classic population genetic models into an epidemiological context; as expected from previous analyses of these types of models (Barton, 1995; Birky and Walsh, 1988; Fisher, 1930; Orr, 2000; Peck, 1994), we find that circulating sublethal deleterious mutations in influenza A/H3N2’s viral population reduce the rate of adaptive evolution and increase the average size of the beneficial mutants that fix.”

Under the Results section, we have now fully situated the results shown in Figure 1 in the context of the theory of asexual evolution, and also cited Strelkowa and Lässig (2012): “Given a specified size distribution for antigenic mutations […] This leads to a reduced tempo of antigenic change, consistent with a reduction in the rate of adaptation that is known from the population genetics literature (Barton, 1995; Peck, 1994).”

Finally, in the Discussion section, we now very briefly mention the consistency of epistasis with the model we present.

2) Stressing the interplay between evolution and epidemiology:

In view of 1), it would be good to focus the paper more on the interplay between the evolutionary and the epidemiological dynamics, which is the novel and most interesting part of the paper. In particular, the authors find that the reduced rate and increased average effect of antigenic mutations may contribute to the linear shape of the influenza tree and to increased attack rates. However, it should become clearer that a strong pruning of antigenic mutations by the joint dynamics with deleterious mutations requires that the latter contribute a substantial fraction of the average fitness variance in the population. This may be the case in model simulations (and the results on tree shape should be compared with the relative contributions of mutation classes to fitness variance). But it is not clear for the actual system. Furthermore, we suggest the authors juxtapose their finding of the influence of deleterious mutations on tree shape with previous explanations of this characteristic.

This is again an excellent comment. We now first highlight the importance of epidemiological dynamics and the competition-reducing nature of antigenic change in the Introduction: “Gaining intuition from these simple models […] a pattern that is inconsistent with the long-term evolutionary dynamics of influenza A/H3N2 in humans.”

In Figure 2, which is the first model we present that is explicitly a multi-strain epidemiological model, we further emphasize the importance of considering the interplay between evolutionary and epidemiological dynamics, and also added the following excerpt to the Results: “The above model contains a number of simplifying assumptions […] population genetic models that assume a constant population size may not be appropriate under certain epidemiological conditions.”

In Figure 6, we have also added several subplots that show more clearly that in our model deleterious mutations contribute a substantial fraction of fitness variance in the viral population. In particular, Figure 6B shows that between 20-80% of the fitness variance is explained by deleterious mutations. How this proportion is calculated is detailed in added text under Materials and methods. Figure 6A (which shows the deleterious mutation load distributed on viral lineages present in the phylogeny) further illustrates that deleterious mutations contribute substantially to fitness variance. This figure can be compared to the added Figure 6E, which shows the distribution of the net reproductive rate (i.e., absolute fitness) on the same viral lineages.

In the Discussion, we now also in detail review the three published epidemiological models that can reproduce the spindly HA tree (Koelle et al. 2006, Bedford et al. 2012, and Ferguson et al. 2003), describing in detail their shortcomings and how their findings can be interpreted in the context of the findings presented here.

For the first model (Koelle et al. 2006), we have the following text: “Influenza’s epochal evolution model (Koelle et al., 2006) assumes that neutral or nearly neutral mutations accumulate at HA epitopes […] and indicate that the evolution of influenza is unlikely to be limited by the occurrence of antigenic mutations.”

For the latter two models, which both reproduce antigenic co-circulation, we have added the following passage: “The generalized cross-immunity model, the canalization model, and the deleterious mutations model presented here therefore share fundamental similarities […] and in light of existing criticisms of the other two models, we therefore argue that the model presented here provides a more plausible mechanistic explanation for influenza’s characteristic evolutionary features.”

3) The effect of deleterious mutations:

In the Discussion, the authors argue that the presence of deleterious mutations reduce clonal competition compared to populations without mutational load. We think this argument is incorrect. According to all models cited above, deleterious mutations should add to the fitness variance in the population and, hence, increase the amount of interference selection. This does not contradict the suggested pruning of antigenic mutations, because in the joint model, antigenic competition is no longer equivalent to fitness competition. The antigenic variance in the population may decrease with increasing rate or effect of deleterious mutations, even if the total fitness variance increases.

This comment reflects different interpretations of the term ‘clonal interference’ in the literature. If ‘clonal interference’ is used to capture all interference effects (both positive and negative), we completely agree with the reviewers. In the original submission, we were using the term ‘clonal interference’ to refer exclusively to interference caused by co-circulating beneficial mutants (as in the work by Gerrish and Lenski). We have clarified this throughout the text, using the term ‘interference effects’ over ‘clonal interference’. We also now emphasize in the Discussion that what appears to be needed to constrain antigenic diversification is that a substantial fraction of fitness variance is generated by deleterious mutations or phenotypes unrelated to antigenicity.

4) Connection to empirical observations and previous models:

One concern is that the model is not particularly data-driven. Although empirical trees are presented inFigure 4to illustrate 3 pathways to antigenic change that are consistent with the model, these trees are more useful as illustrative of concepts than as good evidence for the model. The authors should stress the illustrative character of that figure.

We have stressed the illustrative character of Figure 4 even more than before, indicating that the observed phylogenetic patterns shown in Figure 4are consistent with the proposed molecular pathways of antigenic cluster transitions.

As for (virtually) all theoretical models, it is possible that the patterns in the trees could alternatively be explained by other mechanisms. For example, in the BE92 to WU95 cluster transition, the 145N to 145K substitutions that did not take off globally (that are dead-ends) could have occurred in locations that are thought to not be global source regions (i.e. could the failure of those viruses to persist globally despite beneficial mutations relate to their emergence in less inter-connected non-source regions (i.e., other than SE Asia))? A more thorough discussion of the limits of the model would be in place. For example, in the WU95-SY97 transition, would a scenario not predicted by the model be a longer branch length? How long?

We have added text mentioning that the observed patterns in Figure 4’s trees could also be explained by other mechanisms. For example, for Figure 4A, we now have this paragraph: “Our models above indicate that the failure of these early 145K […] would be necessary to confirm our genetic background hypothesis.”

It's useful to explore how robust a flu model is to other subtypes and hosts, and the authors make an admirable attempt to consider their model in the context of influenza B and swine. However, the most natural comparison would be with seasonal H1N1, and its omission is striking here. The authors also mistakenly assume that H3N2 has similar mutation rates in different hosts, and as a result miss the opportunity to examine how well their model works in a system such as swine (a particularly good example because versions of the human H3N2 virus also circulates in pigs) which have evolutionary rates that are higher than in humans for non-surface proteins (Worobey paper) and should carry a higher mutational load than humans (Discussion, third paragraph). The swine comparison raises ecological considerations that may also relate to human demography.

We appreciate that the referees are interested in how deleterious mutations would further impact the evolutionary dynamics of H1N1 in humans and H3N2 in pigs. We kept this section brief in our first submission, but now extend it further. For seasonal H1N1 in humans, we have added the following paragraph: “Relative to these two influenza viruses, human influenza A/H1N1 shows a similar pattern to H3N2 […] a consequence of differences in these viruses’ global circulation patterns, which have only recently been described (Bedford et al., 2015).”

We have also added text on swine influenza H3N2, referencing de Jong et al.’s finding that the rate of genetic evolution of this virus’s HA is similar to that of human influenza A/H3N2, while its rate of antigenic evolution is 6x slower. Worobey et al.’s 2014 Nature paper infers a significantly higher substitution rate for swine than for humans. Neither gets at differences in l well, though, because both quantify substitution rates, not mutation rates. (The difference between these rates is described nicely in Belshaw, Sanjuan, & Pybus (2011) Current Opinion in Virology.) We now also mention the importance of ecological factors (host lifespan) in generating different selective pressures on the HA in humans vs. pigs (please see the following passage: “Differences in the evolutionary dynamics of influenza viruses also exist across host species. […] being a less important evolutionary driver of HA in short-lived hosts than in humans.”

In the Discussion, it would be important to describe in more detail if the competing models have fundamental issues in capturing the current knowledge or if they would be valid for other parameter choices. It seems odd to write “2%” (Discussion, first paragraph) for a model that is not described in detail and that has presumably free parameters.

We agree with this comment and have altered the text appropriately. We now discuss both the Bedford paper and the Ferguson paper in the context of these models creating fitness variance solely by the co-circulation of beneficial antigenic mutants. Generalized cross-immunity (in the Ferguson paper) and subadditive antigenic effects (in the Bedford paper) both augment competition between circulating beneficial antigenic mutants, and thus result in the generation of spindly HA phylogenies. We removed mention of any specific attack rates.

5) Technical issues:

In general, it is difficult to follow the mathematical description of the model. The description is not particularly complete. For example, inEquation 3: Why does the general background death rate lead to a positive term for the change in the number of susceptible hosts? Or is it only the background death rate of infected individuals? Why does recovery not lead to a positive term here? Another example is the definition of R0,k =0, which differs from the inverse ofEquation 5only by a factor in the exponential function, orEquation 13, which is not a mere modification ofEquation 6, but an extension.

We have reworked the entire Materials and methods section entitled ‘The deleterious mutation-selection balance in an epidemiological context’ to address this concern. This section now details the specific terms in the differential equations, and reorganizes the analytical calculations such that they are more straightforward to follow.

To address the specific questions here: The positive term mN in the equation for the rate of change of susceptible hosts reflects births into the host population (the birth rate is equal to the death rate, and both are quantified by m, such that the host population remains constant in size. Recovery from infection does not lead to a positive term in the equation for the rate of change of susceptible hosts because individuals recover into the recovered (and immune) class. This first model is without antigenic evolution, so is given by an SIR (not SIRS) model. The definition of R0,k=0 is now more straightforward because we now ignore lethal deleterious mutations (such that rho_L is 0).

The statement in the last paragraph of the subsection “Model parameters” is problematic.Equation 14is a set of (deterministic) ODEs and thus does not give a stochastic model. It is well known that many stochastic models can be the basis of a single ODE model, so there is no unique model that the authors simulate. Are the rates inEquation 14directly used, i.e. for a Gillespie algorithm? Or do they reflect net rates that describe the net effect of several rate processes?

We simulated the model stochastically using the Gillespie’s t-leap algorithm, and have added text to that effect.

[Editors' note: further revisions were requested prior to acceptance, as described below.]

In the revised version, the authors have successfully addressed most of the comments. In particular, they have substantially improved the population-genetic side of their findings, which was one of the main concerns of the reviewers.

We again thank the referees for suggesting that we place our findings in the context of population genetic theory on asexual evolution. We believe that this has considerably improved our manuscript.

However, we ask the authors to be more explicit about some of the assumptions in their model. The success or failure of a beneficial antigenic mutation is a complex process, and while the authors do not need to incorporate all aspects into their model, they should be explicit about what is not being factored in:

1) The authors refer to selection as occurring on a ‘strain’, whereas in reality mutations are occurring within a swarm (quasi species) of intrahost diversity. The frequency of mutation as well as back-mutation within this swarm means that beneficial and deleterious mutations are not only arising but also being removed. The impact of a mutation on viral fitness is therefore also a product of the frequency of that mutation within the swarm. The potential bottleneck at transmission also means that not all variants, particularly low-frequency variants, will be transmitted between hosts.

We fully agree with this comment. In recent years, there have been an increasing number of studies that examine influenza’s intrahost diversity and attempt to quantify the size of influenza’s transmission bottleneck (Varble et al. (2014) Cell Host Microbe, Murcia et al. (2010) J. Virology, Murcia et al. (2012) PLoS Pathogens, Hoelzer et al. (2010) J. Virology, Illingowrth et al. (2014) PLoS Computational Biology, Ghedin et al. (2011) J. Infectious Diseases, Hughes et al. (2012) PLoS Pathogens, among others). These include both experimental and observational studies, and consider influenza infections of both human and non-human animal hosts. While these studies are exceptionally valuable, no mathematical models that focus on understanding influenza evolution at the population-level over the long-term have to date integrated this extent of genetic complexity into their model simulations. Instead, the most well known, as well as the most recent, multi-strain models of influenza at the population level have all – without exception that we know of – made the model-simplifying assumption that hosts, when infected, are infected with a single genotype of the virus (e.g. Gog & Grenfell (2002) PNAS, Andreasen et al. (1997) J. Math Bio., Ferguson et al. (2003) Nature, Koelle et al. (2006) Science, Bedford et al. (2012) BMC Biology and papers based on this model, most recently Bedford et al. (2015) Nature). We similarly make this modeling assumption. Although not ‘true’, this assumption greatly simplifies the structure of these types of dynamical systems models. Whether explicitly including intrahost diversity in any of these population-level models would make any difference in terms of predicted dynamics should be theoretically considered at some point in a separate piece of work. The degree of realism that is incorporated into epidemiological models is always to some extent a subjective decision, although it should be based in large part on the question being addressed in the work.

In our manuscript, the fundamental question we address is the role that circulating deleterious mutations play in shaping influenza’s antigenic evolution. What is critical in terms of the model is therefore not how faithfully we model the exact manner in which deleterious mutations arise or are removed from the influenza viral population (this will depend on the size of the transmission bottleneck, the degree of intrahost diversity, and the relationship between within- host and population-level fitness). Instead, what is critical is that the formulation of our model allows for the transient circulation of deleterious mutations in influenza’s viral population. This transient circulation of deleterious mutations has been empirically established through phylogenetic analyses of influenza virus sequence data (Fitch et al. (1997) PNAS, Pybus et al. (2007) Molecular Biology and Evolution). While a highly complex model that incorporates intrahost genetic diversity and transmission bottlenecks of a specified size may more faithfully capture the mechanistic origins of deleterious mutations, this model’s complexity would very likely limit our ability to identify the specific roles that circulating deleterious mutations play in shaping patterns of influenza’s antigenic evolution. Furthermore, there are still too many unknowns related to intrahost influenza diversity and the size of the transmission bottleneck to even accurately include these complexities in model simulations.

While our model does to some extent simplify the manner in which deleterious mutations mechanistically arise and are removed from the influenza viral population, we make the most biologically reasonable assumptions possible given this simplification. Specifically, we assume that deleterious mutations arise at the time of transmission and that they are removed from the population with recovery events of infected individuals. We make the assumption that deleterious mutations arise at the time of transmission to be consistent with population genetic models which assume that deleterious mutations arise at birth. Presumably, deleterious mutations do arise during the period of infection. Indeed, all of the multi-strain influenza models referenced above assume that individuals ‘mutate’ from being infected with one genotype/antigenic phenotype of the virus to being infected with a different genotype/antigenic phenotype of the virus (e.g., I1 → I2). While this formulation may make sense for antigenic mutants, which are likely to have a selective advantage within infected hosts, this formulation is not reasonable for deleterious mutations. This is because if transmission-reducing deleterious mutations emerged de novo within hosts during their infectious period, it is unlikely that they would experience a strong intrahost selective advantage. As such, they would not rise to appreciable frequencies within hosts. Simulating the origination of deleterious mutations as coming from infected individuals ‘mutating’ from being infected with one viral strain carrying k deleterious mutations to being infected with a lower-fitness strain (in terms of transmission) carrying k+1 or more deleterious mutations is therefore not a great modeling choice. We therefore have deleterious mutations be introduced at the time of transmission (‘birth’). For consistency, we similarly model the introduction of antigenic mutations at transmission in the phylodynamic model. As the referees have noted, the assumption of a single genotype within an infected host is not biologically accurate (‘true’). However, short of modeling intrahost genetic diversity and transmission bottlenecks explicitly, we feel that the formulations of our models adopt the most parsimonious assumptions possible in terms of when mutations accumulate.

Although we strongly feel that intrahost genetic diversity or transmission bottlenecks do not need to be explicitly modeled in our submitted work, we now explicitly state our assumption of a genetically homogeneous intrahost viral population (i.e., we explicitly state that we do not incorporate intrahost influenza genetic diversity in our models). We have made this clarification under the Materials and methods section entitled ‘The deleterious mutation-selection balance in an epidemiological context’ (“To extend this basic SIR model to allow for a virus population undergoing deleterious mutations […] is the per-genome per-transmission deleterious mutation rate”).

Under the Materials and methods section entitled ‘Initial and final reproductive rates of an antigenic mutant evolving de novo from a resident viral population’, we have also added: “We denote the resident strain with super- and subscripts r and the antigenic mutant with super- and subscripts m, and use a history-based model (Andreasen et al., 1997) to specify the immunological interaction between these two strains. Note here that we again make the implicit assumption that the intrahost viral population is genetically homogeneous with respect to both antigenicity and the number of deleterious mutations carried.”

Finally, under the Materials and methods section entitled ‘Description and implementation of the phylodynamic model’, we have added the sentence: “Again, we assume for the sake of model simplicity that the intrahost viral population is genetically homogeneous such that a single deleterious mutation load and a single antigenic type suffice in characterizing a viral infection.”

We have further re-read the entirety of the manuscript to look for instances in which we explicitly refer to selection occurring on a ‘strain’, but have not been able to identify any such instances. We do refer to competitive exclusion of strains, as well as between-strain cross- immunity. This is terminology that is commonly used in epidemiological models of influenza; removing this terminology would significantly reduce the clarity of the presented work.

2) We believe the authors are also assuming in this model a consistent immune landscape, for simplicity, which does not reflect variation in immune responses and prior exposures among hosts that affect the fitness of antigenic mutations.

Host variation can have a profound effect on the fitness of a particular mutation (as evidenced nicely by the antigenic mutations in receptor binding sites, which have very different fitnesses in hosts with different immune profiles based on whether antigenic escape or binding avidity is more important.

We are unsure of what the reviewers exactly mean by ‘a consistent immune landscape’. In the phylodynamic model (which we assume is the specific model being referenced in this comment), individuals can and do differ in their infection histories. We model these infection histories in the same way as standard, frequently-used history-based epidemiological models infection histories (see Andreasen et al. (1997) J. Math Bio.). One individual in the population at time t may have been infected with antigenic strain types 1, 2, and 3, for example (each strain also carrying a certain number of deleterious mutations). Another individual at time t may have been infected with only antigenic strain type 2. A third individual at time t may not have been previously infected with any influenza strain (i.e. is considered completely naive). When challenged with a strain through contact with an infected individual, individuals will become infected with a probability that depends on their infection history, as detailed under the ‘Description and implementation of the phylodynamic model’ section of the manuscript. Specifically, the following (already existing) text describes this process: “The model tracked the complete immune history of each host by recording the antigenic type of all previous infections […] The degree of cross-immunity sij between two strains i and j is assumed to be additive, e.g., if strain i gives rise to strain l and strain l gives rise to strain j, then sij = sil + slj.”

We do not include individual variation in host immune responses beyond variability in susceptibility due to infection history. Given contact with an infected individual, previous infection history alone determines susceptibility to infection in our model. The number of deleterious mutations harbored by the virus an individual is infected with alone affects the infectivity of the infected individual (i.e., the number of productive contacts made per unit time that may lead to infection). Both of these model specifications are already detailed under the Materials and methods section ‘Description and implementation of the phylodynamic model’. Neither the infection history nor the number of deleterious mutations harbored by the virus an individual is carrying affects the infected individual’s recovery rate. These are all modeling choices we made in an effort to adopt biologically reasonable assumptions. The modeling choice of previous infection history affecting a host’s susceptibility to infection upon challenge with a new strain is furthermore based on existing epidemiological influenza models present in the literature (Andreasen et al. (1997) J. Math Bio., Ferguson et al. (2003) Nature, Koelle et al. (2006) Science, Bedford et al. (2012) BMC Biology) which make use of this same assumption. Because our submitted work specifically considers the impact of deleterious mutations in shaping patterns of influenza’s antigenic evolution, we do not explicitly incorporate receptor binding avidity. Our first revision, however, does go into considerably detail on how this phenotype can be considered in the context of the model we present. Specifically, in the Discussion section, we have the following text: “Our finding that specifically non-antigenic fitness variation is an important contributing driver in shaping the characteristic features of influenza’s evolutionary dynamics in humans sheds light on other recent virological findings […] suggesting that changes in the receptor binding avidity phenotype that occur with these antigenic mutations may improve virus fitness.”

For increased emphasis, we have revised the Materials and methods section ‘Description and implementation of the phylodynamic model’ to ensure that the readers know that host differ in their immune histories in this full model (please see the following passage: “To explore the full evolutionary dynamics of our model with non-antigenic and antigenic mutations, we implemented an individual-based model […] infected individuals add to their strain history list the antigenic type of the viral infection from which they are recovering.”

3) Successful H3N2 antigenic variants most frequently are produced in Southeast Asia before spreading globally. The spatial ecology of the virus has pronounced effects on which mutations succeed on a global scale and which are not sustained. It would have been simple to examine the proportion of early 145K mutations that arose in SE Asia on the phylogeny (suggested in the original reviewer comments), but short of that the importance of spatial ecology in the global success of a particular variant should be discussed (no additional analysis is required).

First, we would like to apologize to the referees for not addressing this point to the extent that was desired. We are of course familiar with the work that demonstrates that the spatial ecology of the virus has pronounced effects on which mutations succeed on a global scale and the importance of E/SE Asia (and now also India) on the sourcing of antigenic variants (Russell et al. (2008) Science, Lemey et al. (2014) PLoS Pathogens, Bedford et al. (2015) Nature, among several others). The comments that we were provided in the first round of reviews on this point stated:

Connection to empirical observations and previous models:

One concern is that the model is not particularly data-driven. Although empirical trees are presented inFigure 4to illustrate 3 pathways to antigenic change that are consistent with the model, these trees are more useful as illustrative of concepts than as good evidence for the model. The authors should stress the illustrative character of that figure.

As for (virtually) all theoretical models, it is possible that the patterns in the trees could alternatively be explained by other mechanisms. For example, in the BE92 to WU95 cluster transition, the 145N to 145K substitutions that did not take off globally (that are dead-ends) could have occurred in locations that are thought to not be global source regions (i.e. could the failure of those viruses to persist globally despite beneficial mutations relate to their emergence in less inter-connected non-source regions (i.e., other than SE Asia))? A more through discussion of the limts of the model would be in place. For example, in the WU95-SY97 transition, would a scenario not predicted by the model be a longer branch length? How long?

We addressed this point in our first revision by stressing the illustrative character of Figure 4 and adding text that mentions that the observed patterns in Figure 4’s trees could also be explained by other mechanisms. With specific reference to the BE92-to-WU95 cluster transition, we added the following text: “Our models above indicate that the failure of these early 145K […] would be necessary to confirm our genetic background hypothesis.”

Given the reviewers’ comments on our original submission (reproduced above), we did not recognize that they were requesting that we examine the proportion of early 145K mutations that arose in SE Asia on the phylogeny. We have now performed a straightforward additional analysis to examine the geographic locations of the early HA sequences that carried the 145K mutation to determine whether or to what extent spatial ecology may have instead played in the failure of these early 145K lineages.

Of the 395 full-length HA nucleotide sequences in NCBI’s Influenza Virus Resource Database that were isolated between 1993 and 1997 (the date range used for Figure 4A, which examined the BE92-to-WU95cluster transition), 60% were from the US. Only 15% were from Asia (China, Hong Kong, Malaysia, Japan, Singapore). Another 15% were from Europe. The remaining 10% were largely from Oceania. Given such incredibly low sampling from Asian countries, it becomes somewhat meaningless to attempt to do any phylogeographic analyses on the early WU95-like 145K clades that did not succeed. That said, we nevertheless looked at the three largest unsuccessful WU95-like 145K clades present in Figure 4A. The largest of these three clades had 10 sequences. 7 of the sequences were from the US (ranging in geographic extent from Louisiana to Massachusetts, with no state represented more than once). 1 sequence came from Canada. Two sequences came from Europe (Spain and England). Although none of these sequences were from Asia, this does not mean that this clade did not circulate in Asia. (Remember, only 15% of the available sequences came from Asia, so Asia was heavily undersampled.) What this analysis does show is that this clade was geographically widespread, with circulation in both the US and Europe.

The second largest clade had 4 sequences, with Germany, W. Virginia, Finland, and Singapore as the geographic locations of these 4 sequences. The third largest clade had 3 sequences, with Nebraska, Alaska, and Kwangju (South Korea) as the geographic locations of these 3 sequences. Both of these clades were therefore also geographically widespread.

Given these findings, we have added text to the main manuscript in the Results section that alludes to Figure 4A.

4) The authors should clarify that/why beneficial non-antigenic mutations were not included in the model.

In the simpler models (Figures 1, 2, 3, 5), we did not include beneficial non-antigenic mutations. This is because the primary question our manuscript seeks to address is how deleterious mutations can shape patterns of influenza’s antigenic evolution. By including deleterious non- antigenic mutations, we here show that they can greatly impact influenza’s antigenic dynamics. Including beneficial non-antigenic mutations would not alter this result. This is because, similarly to deleterious mutations, beneficial mutations increase fitness variance. Increasing fitness variance in the virus population will make the fate of antigenic mutants even more context-dependent. Please see the following passage in the Discussion: “From the population genetics literature on asexual populations in which interference effects are at play […] The establishment of antigenic mutants can therefore instead lead to long-term coexistence and, as a result of only partial cross-immunity, a larger infected population size.”

Later in the Discussion, we also have written: “While we here simply model this fitness variation as arising from circulating deleterious mutations, any of these non-antigenic phenotypes can similarly contribute to this fitness variation […] influenza virus populations carry substantial deleterious mutation loads (Fitch et al., 1997; Pybus et al., 2007).”

We believe these paragraphs make clear that including beneficial non-antigenic mutations would not alter our qualitative results, and also clarifies why we focus on deleterious mutations in this work.

Note also that in the more complex phylodynamic model (Figures 6-10) we did include beneficial non-antigenic mutations, to avoid Muller’s ratchet. We introduced these mutations through a parameter e, which specified the proportion of non-antigenic mutations that were beneficial rather than deleterious. We added this complexity because it is well known that in a finite population, the most fit genotype can be (repeatedly) lost through a process called Muller’s ratchet, thereby initiating a long-term decrease in viral fitness. We have rewritten the Materials and methods section accordingly (“At a transmission event, the number of non-antigenic mutations was drawn from a Poisson distribution with mean l […] which resulted in a stable deleterious mutation load over time.”

5) There is still a lack of clarity and the role of competition between antigenic variants in the model. We find it hard to agree with the argument that “rather, we propose that many of these circulating antigenic variants ultimately decline from the accumulation of deleterious mutations in the context of an only slowly changing herd immunity landscape”.

In response to this concern, we have examined in greater detail the reasons for loss of the antigenically distinct clades shown in (new) Figure 6B. All of these clades (unsurprisingly) have reproductive values that transition from having a net reproductive rate above 1 to having a net reproductive rate below 1 (new Figure 6D), such that they are excluded. For these clades, we examined their loss of fitness as brought about by deleterious mutation accumulation (new Figure 6A), and we further examined their loss of fitness as brought about by a decline in their number of susceptible hosts (new Figure 6C). Some of these clades did not see a precipitous decline in their number of susceptible hosts, such that their decline in fitness was almost entirely attributed to deleterious mutation accumulation (e.g., the large clade originating in year 20 and going extinct in year 28). For other clades, however, the picture was less clear, and both a decline in the number of susceptible hosts and an increase in the deleterious mutation load were evident. We rephrased the text accordingly (please see the passage: “In our model, these antigenic variants […] suffice in generating a sufficiently high-fitness viral lineage that will ensure its long-term evolutionary success”).

Thus, our inclusion of the new Figure 6C, alongside Figure 6A, makes it more evident that antigenic evolution and purifying selection act in concert with one another in shaping influenza’s phylodynamics.

eLife is a non-profit organisation inspired by research funders and led by scientists. Our mission is to help scientists accelerate discovery by operating a platform for research communication that encourages and recognises the most responsible behaviours in science.eLife Sciences Publications, Ltd is a limited liability non-profit non-stock corporation incorporated in the State of Delaware, USA, with company number 5030732, and is registered in the UK with company number FC030576 and branch number BR015634 at the address:
eLife Sciences Publications, Ltd
1st Floor, 24 Hills Road
Cambridge CB2 1JP
UK