Abstract

Since its emergence in 1968, influenza A (H3N2) has evolved extensively in genotype and antigenic phenotype. However, despite strong pressure to evolve away from human immunity and to diversify in antigenic phenotype, H3N2 influenza shows paradoxically limited genetic and antigenic diversity present at any one time. Here, we propose a simple model of influenza that displays rapid evolution but low standing diversity and simultaneously accounts for the epidemiological, genetic, antigenic and geographical patterns displayed by the virus. In this model, antigenic phenotype is represented by a N-dimensional vector and virus mutations perturb phenotype within this continuous Euclidean space. We implement this model in a large-scale individual-based simulation, and in doing so, find a remarkable correspondence between model behavior and observed influenza dynamics. We find that evolution away from existing human immunity results in rapid population turnover in the influenza virus and that this population turnover occurs primarily along a single antigenic axis. Thus, selective dynamics induce a canalized evolutionary trajectory, in which the evolutionary fate of the influenza population is surprisingly repeatable and hence, in theory, predictable.

Introduction

Epidemic influenza is responsible for between 250,000 and 500,000 global deaths annually, with influenza A (and in particular, A/H3N2) having caused the bulk of human mortality and morbidity [1]. Influenza A (H3N2) has continually circulated within the human population since its introduction in 1968, exhibiting recurrent seasonal epidemics in temperate regions and less periodic transmission in the tropics. During this time, H3N2 influenza has continually evolved both genetically and antigenically. Most antigenic drift is thought to be driven by changes to epitopes in the hemagglutinin (HA) protein [2]. Phylogenetic analysis of the genetic relationships among HA sequences has revealed a distinctive genealogical tree showing a single predominant trunk lineage and side branches that persist for only 1–5 years before going extinct [3]. This tree shape is indicative of serial replacement of strains over time; H3N2 influenza shows rapid evolution, but low standing genetic diversity.

This observation has remained puzzling from an epidemiological standpoint. Antigenic evolution occurs rapidly and strong diversifying selection exists to escape from human immunity; why then do we see serial replacement of strains rather than continual genetic and antigenic diversification? Indeed, simple epidemiological models show explosive diversity of genotype and phenotype over time [4, 5]. Previous work has sought model-based explanations of the limited diversity of influenza, relying on short-lived strain-transcending immunity [4, 5], complex genotype-to-phenotype maps [6] or a limited repertoire of antigenic phenotypes [7].

Experimental characterization of antigenic phenotype is possible through the hemagglutination inhibition (HI) assay, which measures the cross-reactivity of HA from one virus strain to serum raised against another strain [8]. The results of many HI assays can be combined to yield a two-dimensional map, representing antigenic similarity and distance between strains as an easily visualized and quantified measure [9]. The path traced across this map by influenza A (H3N2) from 1968 until present is largely linear, showing serial replacement of one strain by another; there are no major bifurcations of antigenic phenotype [9].

Herein, we seek to simultaneously model the genetic and antigenic evolution of the influenza virus. We represent antigenic phenotypes as points in a N-dimensional Euclidean space. Based on the finding that a two-dimensional map adequately explains observed antigenic distance between strains [9], we begin with antigenic phenotypes as points on a plane, but relax this assumption later on in the analysis. After exposure to a virus, a host’s risk of infection is proportional to the Euclidean distance between the infecting phenotype and the closest phenotype in the host’s immune history. Mutations perturb antigenic phenotype, moving phenotype in a random radial direction and for a randomly distributed distance. We implemented this geometrical model in a large-scale individual-based simulation intended to directly model the antigenic map and genealogical tree of the global influenza population. The simulation includes multiple host populations with different seasonal forcing, hosts with complete immune histories of infection, and viruses with antigenic phenotypes. As the simulation proceeds, infections are tracked and a complete genealogy connecting virus samples is constructed. Results shown here are for simulations of 40 years of virus evolution in a population of 90 million hosts.

Results

At the end of most simulations, there remain only a few closely related viral lineages, indicating that genealogical diversity is restricted by evolution in the two-dimensional antigenic landscape. However, in 20 out of the 100 replicate simulations, we observed a major bifurcation of antigenic phenotype and a consequent increase in incidence and genealogical diversity. These simulations were removed from the analysis. Similar to Koelle et al. [10], we assume that although the historical evolution of H3N2 influenza followed the path of a single lineage, it could have included a major bifurcation. Here, we focus on one representative simulation out of the 80 low-diversity outcomes.

The virus persists over the course of the 40-year simulation, infecting a significant fraction of the host population through annual winter epidemics in temperate regions and through less periodic epidemics in the tropics (Figure 1A). Across replicate simulations, we observe average yearly attack rates of 6.8% in temperate regions and rates of 7.1% in the tropics, comparable with estimated attack rates of influenza A (H3N2) of 3–8% per year [11, 12]. Over the course of the simulation, the virus population evolves in antigenic phenotype exhibiting, at any point, a handful of highly abundant phenotypes sampled repeatedly and a large number of phenotypes appearing at low abundance (Figure 1B). By including measurement noise on antigenic locations, we approximate an experimental antigenic map of H3N2 influenza (Figure 1D). Over the 40-year simulation, antigenic drift moves the virus population at an average rate across replicates of 1.05 antigenic units per year, corresponding closely to the empirical rate of 1.2 units per year [9]. The appearance of clusters in the antigenic map comes from the regular spacing of high abundance phenotypes combined with measurement noise. Over time, clusters of antigenically similar strains are replaced by novel clusters of more advanced strains (Figure S1A). Across replicate simulations, clusters persist for an average of 5.0 years measured as the time it takes for a new cluster to reach 10% frequency, peak and decline to 10% frequency. The transition between clusters occurs quickly, taking an average of 1.8 years.

Remarkably, although antigenic phenotype is free to mutate in any direction in the two-dimensional space, selection pressures force the virus population to move in nearly a straight line in antigenic space (Figure 1B). Across replicate simulations, 94% of the variance of antigenic phenotype can be explained by a single dimension of variation. This mirrors the empirical results showing a largely linear antigenic map for H3N2 influenza isolates from 1968 to 2003 [9]. Because of the primarily one-dimensional movement, antigenic distance from the original phenotype increases nearly linearly with time (Figure S1B). Antigenic evolution occurs in a punctuated fashion; periods of relative stasis are interspersed with more rapid antigenic change (Figure S1B). Antigenic and epidemiological dynamics show a fundamental linkage, so that large jumps of antigenic phenotype result in increased rates of infection (Figure 1, Figure S2). In general, evolution via many smaller mutations results in a smoother antigenic map and less variation in yearly epidemics (Figure S3), while evolution via rare mutations of large effect exhibits a more clustered antigenic map and wider variation in seasonal incidence (Figure S4).

The genealogical tree connecting the evolving virus population appears characteristically sparse with pronounced trunk and short-lived side branches (Figure 1C). This tree shape is reflected in low levels of standing diversity; across replicates, an average of 5.68 years of evolution separate two randomly sampled viruses in the population. This result matches well with the average diversity observed in influenza A (H3N2) of 5.65 years separating randomly sampled contemporaneous viruses [13]. A spindly genealogical tree is indicative of population turnover, wherein novel antigenic phenotypes continually replace more primitive ‘spent’ phenotypes, purging their genealogical diversity. In general, natural selection reduces effective population size and genealogical diversity [13].

Selective pressures can be examined by comparing which mutations fix, i.e. are incorporated into the progenitor trunk lineage, and which mutations are lost, i.e. incorporated into side branches bound for extinction. This approach has shown that, in influenza A (H3N2), natural selection promotes mutations to epitope sites in the HA1 region [14, 15]. By examining antigenic mutations, we find a corresponding effect in simulated evolutionary trajectories (Table 1). Trunk mutations tend to push antigenic phenotype forward along the line of primary antigenic variation (Figure S5). We find a roughly linear relationship between the antigenic effect of a mutation and the likelihood that this mutation becomes incorporated into the trunk (Figure S6). Additionally, we find that trunk mutations occur at strikingly regular intervals, with less variation of waiting times than expected under a simple random process (Figure S7). There is a relative scarcity of mutation events occurring in intervals under 1 year and a relative excess of a mutation events occurring in 2–3 year intervals (Figure S7). This pattern arises from clonal interference between competing mutations which reduces variability in the fixation process of adaptive substitutions [16].

The genealogical tree also contains detailed information on the history of migration between regions. We find that, consistent with empirical estimates [17, 18], the trunk resides primarily within the tropics, where seasonal dynamics are less prevalent (Figure 2A). Across replicate simulations, we observe 72% of the trunk’s history within the tropics and 28% within temperate regions. With symmetric host contact rates and equal host population sizes, and without seasonal forcing, we would expect trunk proportions of one third for each region. We calculated rates of migration based on observed event counts across replicate simulations, separating region-specific rates on side branches from region-specific rates on trunk branches. We find that migration patterns on side branches are close to symmetric, with similar rates between all regions, while migration patterns on trunk branches are highly asymmetric, with high rates of movement between temperate regions and from temperate regions into the tropics (Figure 2B). Extrapolating from these rates, we arrive at an expected stationary distribution of trunk location of 76% tropics and 24% temperate regions, in line with the observed residency patterns of the trunk. It may at first seem counter-intuitive to see higher rates of movement from the temperate regions into the tropics along trunk branches, but it makes sense when thought of in terms of conditional probability. Only those lineages that remain in the tropics, migrate into the tropics or those lineages which rapidly migrate between the north and south have a chance at becoming the trunk lineage, while lineages that remain within the temperate regions are doomed to extinction. Along similar lines, Adams and McHardy [19] use a modeling approach to show the importance of non-seasonal transmission to the evolution of the virus.

These findings suggest that persistence and migration are fundamentally connected and have important implications for future phylogeographic analyses. Russell et al. [17] emphasize a source-sink model of movement of the HA protein of influenza A (H3N2) based on their finding of a trunk lineage residing within China and the Southeast Asian tropics. Whereas, Bedford et al. [18] emphasize a global metapopulation model based on phylogenetic inference of migration rates across the entire tree. Our results suggest that both scenarios are simultaneously possible; side branches may be highly volatile moving rapidly and symmetrically between regions, while the trunk lineage may be more stable remaining within a region (or within a highly connected network of regions) that has more persistent transmission. In light of these results, we suggest that future work on the phylogeography of influenza take into account trunk vs. side branch differences in migration patterns.

Discussion

Correspondence between model and data

Although multiple epidemiological/evolutionary mechanisms have been proposed to explain the restricted genetic diversity and rapid population turnover of influenza A (H3N2) [4, 5, 6, 7], our results show that a simple model coupling antigenic and genealogical evolution exhibits broad explanatory power. We find a strong correspondence between the antigenic and genealogical patterns generated by our model (Figure 1) and patterns of genetic and antigenic evolution exhibited by influenza A (H3N2) [3, 9]. Our model simultaneously captures seasonal attack rates, the rate and pattern of antigenic drift, genealogical diversity and geographic migration patterns.

Our model suggests that punctuated antigenic evolution need only be explained by a lack of mutational opportunity and predicts that more detailed classification of influenza strains will support a relatively small number of predominant phenotypes (Figure 1B). We suggest that a large proportion of intra-cluster variation in the observed antigenic map is due to experimental noise, rather than each strain possessing a unique antigenic location. Additionally, our model accurately predicts the contrasting dynamics of other types/subtypes of influenza. We find that lowering mutation size/effect or lowering intrinsic R₀ results in decreased incidence, slower antigenic movement and greater genealogical diversity, all distinguishing characteristics of H1N1 influenza and influenza B (Figure S8, Figure S9).

The historical record of influenza evolution suggests that bifurcation of viral lineages is rare, but possible. We have observed no bifurcations in H1N1 influenza from 1918 to 1957 and again from 1977 to 2010, no bifurcations in H2N2 influenza from 1957 to 1968 and no bifurcations of H3N2 influenza from 1968 to 2010. We have observed one bifurcation in influenza B from 1980 to 2010. Thus, ignoring differences between influenza types and subtypes, we have very roughly observed a rate of one major bifurcation in 155 years of evolution. In our model, in 20 out of 100 replicate simulations, we observe a deep bifurcation in the viral genealogy, which translates to observing one deep bifurcation in 200 years of evolution. Thus, we suggest that the 20 of 100 simulations where deep branching occurs are not necessarily evidence of poor model fit.

In our model, when antigenic phenotype remains static, there may be multiple consecutive seasons without appreciable incidence (Figure 1A), a pattern apparently absent from H3N2 influenza [20]. Additionally, we observe antigenic trajectories that are more linear and deterministic than the highly clustered trajectory observed by Smith et al. [9]. We suggest that any model exhibiting punctuated evolution broadly consistent with the punctuated change seen in the experimental antigenic map will show similar patterns of incidence. We can ‘fix’ the incidence patterns, but at the cost of too smooth an antigenic map (Figure S3). Evolutionary patterns of the neuraminidase (NA) protein may provide an explanation. Epitopes in the HA and NA proteins are jointly responsible for determining antigenicity [2], and it is now clear that levels of adaptive evolution are similar between HA and NA [21]. Thus, changes in NA may be driving incidence patterns as well, resulting in an observed timeseries of incidence partially divorced from the antigenic map of HA. Incorporating antigenic evolution of NA could thus yield a rougher antigenic map for HA, more closely matching experimental results, while simultaneously yielding smoother year-to-year incidence patterns.

It remains a central question as to the extent that short-lived strain-transcending immunity is responsible for influenza’s limited diversity and spindly genealogical tree [4, 5]. Our findings suggest a possible resolution. Although lacking short-lived immunity, our model shows a detailed correspondence to both the antigenic map and genealogical tree of H3N2 influenza. If an antigenic map were to show a deep bifurcation, where two viral lineages move in different antigenic directions, then we would expect the same bifurcation to be evident in the genealogical tree. Short-lived strain-transcending immunity provides a mechanism by which lineages may diverge in antigenic phenotype, but still show epidemiological interference. This mechanism would explain a situation where bifurcations emerge in the antigenic map, but competition results in the extinction of divergent antigenic lineages. The empirical antigenic map [9] suggests that this is not the case; one cluster leads to another cluster in orderly succession and there is never competition between antigenically distant clusters. This supports the hypothesis that antigenic evolution is primarily limited by a lack of mutational availability. This is not to say that short-lived strain-transcending immunity is not present; observed interference between subtypes [4, 22] and evolution at CTL epitopes [23] provides substantial evidence for its existence. Instead, we suggest that short-lived strain-transcending immunity does not automatically generate antigenic maps and genealogical trees consistent with empirical evidence.

Linear antigenic movement

It would seem possible for one viral lineage to move in one antigenic direction, while another lineage moves tangentially, eventually resulting in two non-interacting viral lineages. Instead, we find that only movement in a single antigenic direction is favored. The origins of this pattern can be seen in the interaction between virus evolution and host immunity (Figure 3). As the virus population evolves forward it leaves a wake of immunity in the host population, and evolution away from this immunity results in the canalization of the antigenic phenotype; mutations that continue along the line of primary antigenic variation will show a transmission advantage compared to more tangential mutations.

Following the work of Smith et al. [9], it remained an open question of why a two-dimensional map should explain the antigenic variation of H3N2 influenza. Although the authors astutely speculated that “there is a selective advantage for clusters that move away linearly from previous clusters as they most effectively escape existing population-level immunity, and this is a plausible explanation for the somewhat linear antigenic evolution in regions of the antigenic map.” This hypothesis remained to be tested. Here, we show from a simple model of epidemiology and evolution that a linear trajectory of antigenic evolution dynamically emerges due to basic selective pressures. This result simultaneously explains the linear pattern of antigenic drift [9] and the characteristically spindly genealogical tree [3] exhibited by influenza A (H3N2).

To consider to what extent these results were contingent on the dimensionality of the underlying antigenic model, we further implemented our model in a 10-dimensional antigenic space. Here, mutations occur as 10-spheres, but the distance moved by a mutation is the same as in the previous two-dimensional formulation. We arrive at nearly the same results with this model; principal components analysis shows that the first and second dimensions of variation account for 87% and 7%, respectively, of the total variance (Figure S10). Thus, our model predicts that future work probing mutational effects will support an underlying high-dimensional antigenic space, even though a two-dimensional map is sufficient to explain observed antigenic relationships among strains.

Winding back the tape

It seems clear that, in our model, selection reduces the degrees of freedom of antigenic evolution. In light of this, we wanted to examine the degree of stochasticity in replicated evolutionary trajectories, and thereby test what happens when we “wind back the tape” [24] on the evolution of the virus. We ran 100 replicate simulations, each starting from the endpoint of the original 40-year simulation (Figure 4, Figure 5). Initially, we find a great detail of repeatability; during the first year of evolution, every replicate virus population undergoes a similar antigenic transition (Figure 4), resulting in a repeatable peak in northern hemisphere incidence (Figure 5). After three years, repeatability has mostly disappeared, with antigenic phenotype and incidence appearing highly variable across replicates (Figure 4, Figure 5). The 1–2 year timescale of repeatability can be explained by the presence of standing antigenic variation. In the initial virus population, there are several novel antigenic variants present at low frequency (Figure 3), one of which, without fail, comes to predominate the virus population.

We see that the initial evolutionary trajectory, during which time standing variation plays out, is highly repeatable, and thus predictable given enough information and the right methods of analysis. However, prediction of longer-term evolutionary scenarios will necessarily be difficult or impossible except in a vague sense. Through careful surveillance efforts and genetic and antigenic characterization of influenza strains, the World Health Organization makes twice-yearly vaccine strain recommendations [25]. It should be possible to combine these sorts of modeling approaches with surveillance data to gauge the likelihood that a sampled variant will spread through the population.

Recent work on empirical fitness landscapes has shown that natural selection follows few mutational paths [26]. The spindly genealogical tree and the almost linear serial replacement of influenza strains has remained a puzzling phenomenon. We suggest that the evolutionary and epidemiological dynamics displayed by the influenza virus may simply be explained as an outgrowth of selection to avoid host immunity. Natural selection can only ‘see’ one step ahead, and so sacrifices long-term gains for short-term advantages. The result is a canalized evolutionary trajectory lacking antigenic diversification.

Materials and Methods

Transmission model

To characterize the joint epidemiological, genealogical, antigenic and spatial patterns of influenza, we implemented a large-scale individual-based model (full simulation source code available at http://github.com/trvrb/antigen). This model consists of daily time steps, in which the states of hosts and viruses are updated. Hosts may be born, may die, may contact other hosts allowing viral transmission, or may recover from infection. Viruses may mutate in antigenic phenotype. Each simulation ran for 40 years of model time.

Hosts in this model are divided between three regions: North, South and Tropics. There are 30 million hosts within each of the three regions, giving N = 9 × 10⁷ hosts. Host population size remains fixed at this number, but vital dynamics cause births and deaths of hosts at a rate of 1∕30 years = 9.1 × 10^-5 per host per day. Within each region, transmission proceeds through mass-action with contacts between hosts occurring at a rate of β = 0.36 per host per day. Regional transmission rates in temperate regions vary according to sinusoidal seasonal forcing with amplitude ϵ = 0.15 and opposite phase in the North and in the South. Transmission rate does not vary over time in the Tropics. Transmission between region i and region j occurs at rate mβ_i, where m = 0.001 and is the same between each pair of regions and β_i is the within-region contact rate. Hosts recover from infection at rate ν = 0.2 per host per day, so that R₀ in a naive host population is 1.8. An individual cannot be simultaneously infected with multiple virus lineages; there is no super-infection in the model.

Each virus possesses an antigenic phenotype, represented as a location in Euclidean space. Here, we primarily use a two-dimensional antigenic location. After recovery, a host ‘remembers’ the phenotype of its infecting virus as part of its immune history. When a contact event occurs and a virus attempts to infect a host, the Euclidean distance from infecting phenotype ϕ_v is calculated to each of the phenotypes in the host immune history ϕ_h1,…,ϕ_hn. Here, one unit of antigenic distance is designed to correspond to a twofold dilution of antiserum in a hemagglutination inhibition (HI) assay [9]. The probability that infection occurs after exposure is proportional to the distance d to the closest phenotype in the host immune history. Risk of infection follows the form ρ = min{ds,1}, where s = 0.07. Cross-immunity σ equals 1 - ρ. The initial host population begins with enough immunity to slow down the initial virus upswing and place the dynamics closer to their equilibrium state; initial R was 1.28.

Our model follows Lin et al. [27] and Gog and Grenfell [28] in representing antigenic distance as distance between points in a geometric space. By forcing one-dimension to directly modulate β, Gog and Grenfell find an orderly replacement of strains. Kryazhimskiy et al. [29] use a two-dimensional strain-space, but enforce a cross-immunity kernel that directly favors moving along a diagonal line away from previous strains. Our model does not ‘build in’ the one-dimensional direction of antigenic drift, which instead emerges dynamically from competition among strains.

The initial virus population consisted of 10 infections each with the identical antigenic phenotype of {0,0}. Over time viruses evolve in antigenic phenotype. Each day there is a chance μ = 10^-4 that an infection mutates to a new phenotype. This mutation rate represents a phenotypic rate, rather than genetic mutation rate, and can be thought of as arising from multiple genetic sources. When a mutation occurs, the virus’s phenotype is moved in a completely random direction ~ Uniform(0,360) degrees. Mutation size is sampled from the distribution ~ Gamma(α,β), where α and β are chosen to give a mean mutation size of 0.6 units and a standard deviation of 0.4 units. This distribution is parameterized so that mutation usually has little effect on antigenic phenotype, but occasionally has a large effect. This is similar to the neutral networks implemented by Koelle et al. [6], wherein most amino acid changes result in little decrease to cross-immunity between strains, but some changes result in large jumps in cross-immunity. This mutation model implicitly assumes independence between mutations; each mutation’s effect is independent of genetic background.

Model output

Daily incidence and prevalence are recorded for each region. During the course of the simulation, samples of current infections are taken from the evolving virus population at a rate proportional to prevalence. Each viral infection is assigned a unique ID, and in addition, infections have their phenotypes, locations and dates of infection recorded. In this model, viruses lack sequences and so standard phylogenetic inference of the evolutionary relationships among strains is impossible. Instead, the viral genealogy is directly recorded. This is made possible by tracking transmission events connecting infections during the simulation; infections record the ID of their ‘parent’ infection. Proceeding from a sample of infections, their genealogical history can be reconstructed by following consecutive links to parental infections. During this procedure, lineages coalesce to the ancestral lineages shared by the sampled infections, eventually arriving at the initial infection introduced at the beginning of the simulation. Commonly, phylodynamic simulations generate sequences that are subsequently analyzed with phylogenetic software to produce an estimated genealogy [4, 6, 30]. This step of phylogenetic inference is imperfect and computationally intensive, and by side-stepping phylogenetic reconstruction we arrive at genealogies quickly and accurately. Other authors have implemented similar tracking of infection trees [31, 32]. This genealogy-centric approach makes many otherwise difficult calculations transparent, such as calculating lineage-specific region-specific migration rates (Figure 2) and lineage-specific mutation effects (Figure S5).

Infections are sampled at a rate designed to give approximately 6000 samples over the course of the simulation, with genealogies constructed from a subsample of approximately 300 samples. The results presented in Figure 1 represent a single representative model output; one hundred replicate simulations were conducted to arrive at statistical estimates.

Parameter selection and sensitivity analysis

Estimating what the basic reproductive number R₀ for seasonal influenza would be in a naive population is notoriously difficult. Season-to-season estimates of effective reproductive number R for the USA and France gathered from mortality timeseries display an interquartile range of 0.9–1.8 [33]. Geographic spread within the USA suggests an average seasonal R of 1.35 [34]. These estimates of R will be lower than the R₀ of influenza due to the effects of human immunity. We assumed R₀ of 1.8, consistent with the upper range of seasonal estimates. Duration of infection was chosen based on patterns of viral shedding shown during challenge studies [35]. The linear form of the risk of infection and its increase as a function of antigenic distance s was chosen as 0.07 based on experimental work on equine influenza [36] and from studies of vaccine effectiveness [37]. Between-region contact rate m was chosen to yield a trunk lineage that resides predominantly in the tropics. With much higher rates of mixing, the trunk lineage ceases to show a preference the tropics, and with much lower rates of mixing, particular seasons in the north and the south will often be skipped. The amplitude of seasonal forcing ϵ was chosen to be just large enough to get consistent fade-outs in the summer months and is consistent with empirical estimates [38].

Mutational parameters were based, in part, on model behavior. We assumed 10 amino acid sites involved in antigenicity, each mutating at a rate of 10^-5 [39] to give a phenotypic mutation rate μ = 10^-4 per infection per day. We chose mutational effect parameters (mean = 0.6, sd = 0.4) that would give suitably fast rates of antigenic evolution corresponding to approximately 1.2 units of antigenic change per year, while simultaneously giving clustered patterns of antigenic evolution [9]. Similar outcomes are possible under a variety of parameterizations. If mutations are more common (μ = 3 × 10^-4) and show less variation in effect size (mean = 0.6, sd = 0.2), then antigenic drift occurs in a more continuous fashion, resulting in less variation in seasonal incidence and a smoother distribution of antigenic phenotypes (Figure S3). If mutations are less common (μ = 5 × 10^-5) and show more variance in effect (mean = 0.7, sd = 0.5), then antigenic change occurs in a more punctuated fashion (Figure S4). Basic reproductive number R₀ can be traded off with mutational parameters to some extent. Less mutational input and higher R₀ will give similar patterns of antigenic drift and seasonal incidence. Similarly, Kucharski and Gog [40] find that increasing R₀ results in increased rates of emergence of antigenically novel strains.

Antigenic map

Antigenic phenotypes are modeled as discrete entities on the Euclidean plane; multiple samples have the same antigenic location. However, in the empirical antigenic map of influenza A (H3N2), each strain appears in a unique location [9]. We would argue that some of this pattern comes from experimental noise. Indeed, Smith et al. [9] find that observed measurements and measurements predicted from the map differ by an average of 0.83 antigenic units with a standard deviation of 0.67 antigenic units. We take this as a proxy for experimental noise and add jitter to each sampled antigenic phenotype by moving it in a random direction for an exponentially distributed distance with mean of 0.53 antigenic units. If two samples with the same underlying antigenic phenotype are jittered in this fashion, the distance between them averages 0.83 antigenic units with a standard deviation of 0.64 units.

We added noise to each of the 5943 sampled viruses in this fashion resulting in an approximated antigenic map (Figure 1D). Virus samples were then clustering following standard clustering algorithms. We tried clustering by the k-means algorithm and also agglomerative hierarchical clustering with a variety of linkage criterion. We found that clustering by Ward’s criterion consistently outperformed other methods, when measured in terms of within-cluster and between-cluster variances. However, the exact clustering algorithm had little effect on our overall results.

Acknowledgments

We would like to thank Sarah Cobey, Aaron King, Pejman Rohani and the attendees of the 2011 RAPIDD Workshop on Phylodynamics for helpful discussion. We would also like to thank Ed Baskerville and Daniel Zinder for programming advice. The term ‘canalization’ was originally suggested by Micaela Martinez-Bakker.

Funding

TB is supported by the Howard Hughes Medical Institute and by the European Molecular Biology Organization. AR works as a part of the Interdisciplinary Centre for Human and Avian Influenza Research (ICHAIR) and the University of Edinburgh’s Centre for Immunity, Infection and Evolution (CIIE). MP is an investigator of the Howard Hughes Medical Institute.

References

Supporting Information

---

---

---

---

---

---

---

---

---

---

---

---

---

---

---

---

---

---

---