The variation of influenza by mutation, i.e., antigenic drift, was soon recognized as the basis for the difficulty of predicting epidemics [5], and better understanding the quantitative impact of drift on seasonal epidemics has become an important goal in influenza research. Given that varying levels of cross-protection by earlier viruses are likely to be an important factor in epidemic severity, one experimental strategy has been to immunize subjects with known strains and determine susceptibility to a challenge with another known strain. McLaren et al. [6] vaccinated ferrets with A/Hong Kong/1/68 or A/England/42/72, later challenged with a A/England/42/72-like virus, and indeed found greater protection by the matching strain. A similar, more extensive study has been performed on horses using equine influenza virus in order to estimate the parameters of an epidemiological model of outbreak risk as a function of immune escape [7]. In a study with human volunteers immunized with one of four influenza strains isolated in 1968, 1972, 1973, or 1974, and then challenged with a 1974 strain, Potter et al. [8] found that susceptibility increased linearly with the number of years separating the challenge and vaccine strains. Consistent results were obtained in a similar study by Larson et al. [9].

Field studies of vaccine efficacy have confirmed the results of these experimental approaches (c.f. [10]). As expected, efficacy is substantially higher when the vaccine strain matches the dominant seasonal strain than when there is a measurable mismatch between the vaccine strain and the dominant strain [11]. Gupta et al. [12] compiled an extensive set of vaccine efficacy data and found a strong correlation (R² = 0.81) between vaccine efficacy in a given season and the fraction of amino acid replacements in the antigenic regions of the hemagglutinin protein between the vaccine strain and the dominant circulating strain. Evidence of a more direct quantitative connection between the antigenic novelty of the viruses from a given season and susceptibility to influenza can be found in a prospective study of the 1976 seasonal epidemic by Gill and Murphy [13]. The researchers determined the most recent year of laboratory confirmed influenza for a set of volunteers and found a strong linear relationship between the probability of infection in the 1976 season and the number of years since the previous infection. The above studies suggest that the number of susceptible individuals for a given season, and hence the number of influenza cases, will be strongly correlated with measures of the antigenic distance between prevailing viruses for that season and viruses that circulated in previous seasons.

In light of the compelling evidence for a quantitative relationship between antigenic drift and seasonal severity, it is perhaps surprising that until recently no study demonstrated a positive correlation between drift and severity using data derived directly from influenza surveillance. A recent study by Wu et al. [14] focuses on  influenza H1N1 and considers major "antigenic strains" based on hemagglutination inhibition (HI) assays -- a method traditionally used by the Centers for Disease Control (CDC) to characterize seasonal viruses for vaccine selection and surveillance purposes. Wu et al. then estimate the excess mortality attributable to each antigenic strain using the fraction of isolates positive for that strain among all those tested in seasons in which the strain was detected. This analysis has shown a strong correlation between estimates of the multi-season excess mortality due to an antigenic strain and its antigenic distance from previous strains for H1N1, and a weaker correlation for H3N2. Given that in recent years, most of the inter-pandemic influenza morbidity was caused by H3N2 viruses [15], we were interested in investigating the relationships between sequence and antigenic variation in the H3N2 HA and the seasonal severity (health burden) of influenza caused by this subtype. In this analysis, we aimed at reaching beyond correlations and developing a practical model for projecting the health burden of influenza H3N2.

To test the models statistically, we used two variants of a jackknife resampling procedure [28]. The first variant tests the statistical robustness of the relationships between the explanatory variables and the target variable in general. The second tests the statistical robustness of a model using a particular set of explanatory variables.

In the first variant we removed data for each single season from our dataset one-by-one and used the remaining 15 seasons to find the best-fitting model (both the set of the explanatory variables and coefficients). Then we used this model to predict the epidemiology for the remaining season (leave-one-out jackknife scheme). The sum of squares of the deviations of these predictions from the observations was accumulated over all 16 seasons to serve as the robustness indicator.

In the second variant we generated all possible bipartitions of the 16 epidemiological observations into two classes of equal size. With each bipartition, one half of the seasons data was used to train the model using the fixed set of the explanatory variables; then these model coefficients were used to predict the epidemiology for the remaining 8 seasons (leave-half-out jackknife scheme). The residual sums of squares for the control half of the data were accumulated over all bipartitions to serve as the indicator of the model robustness. We tested the top 10 models from each class and ranked them by decreasing total residual sum of squares.

In the leave-one-out jackknife test the 5-variable model family demonstrates a highly robust behavior, accounting for 0.96 of the of the variance in epi.PC1. The following 5-variable model demonstrated the best performance among all tested models:

epi.PC1’ ~ n.n2.PC1 + n.n2.PC2 + n1.n1.PC1 + a.n1.n1.PC1 + a.n1.n1.PC2 + 0

Hence this model includes genetic distance between the current season and 2 seasons before, and intra-season genetic and serological distances between isolates circulating in the previous season in the US (genetic) or North America and Northern Hemisphere (serological). When trained on the full complement of 16 seasons (Figure 4), this model explains 0.98 of the original variance (F-statistics p-value of 3x10^-9), the coefficients being -1.36, -1.54, -0.32, -1.28 and -1.69 respectively (PCA gives arbitrary signs to the principal components, so the signs of coefficients do not indicate the kind of relationship by themselves). In the leave-half-out jackknife test it explains, on average, 0.89 of the variance in epi.PC1.

To further ascertain the statistical validity of the relationship between the genetic and serological data and epidemiological severity we implemented an additional test procedure. The epidemiological data vector (epi.PC1) was permuted and all combinations of 5 out of 26 explanatory variables were tested as potential reconstruction models for the scrambled epi.PC1 data. Only 36 out of 100,000 permutations yielded a fit as good as that obtained with the unpermuted data.

The above shows that the sequence and antigenic data explain part of the variance in morbidity among seasons.  We can reject the hypothesis that the true R² equals zero with high statistical confidence, and we estimate that the true R² is 0.96. In order to obtain a statistical lower bound on the true R², we generated epidemiological severity values that mimicked various levels of unexplained variance, and determined the lowest R² value that could not be rejected at the 5% level by the R² observed for the real data (i.e., the lowest value for which the R² of the fit was as good as the observed 98% in at least 2.5% of the replicates).  To the epidemiological values predicted by the fitted model, which can be fit with R² = 1, we added independent pseudo-random numbers.  Two distributions for these were used: a normal distribution and a Laplace (double exponential) distribution.  The Laplace distribution has an excess kurtosis of 3.0, and is used because the observed residuals have excess kurtosis (approximately 2.5).  For each trial quantity of added noise, many replicates were generated, and each was independently fit to the sequence and antigenic data, with no constraint on the choice of combinations of explanatory variables.

In the second stage, we apply the standard stepwise model reduction approach to find the optimal set of original variables. We started with a model containing 8 variables based on the best-fit model identified in the 1st stage (se.n.n2, sn.n.n2, ss.n.n2, se.n1.n1, sn.n1.n1, ss.n1.n1, aa.n1.n1 and an.n1.n1) and removed one by one the variables that contribute the least to the model (i.e. with the lowest absolute t-value of the coefficient). Each derived model was compared to its parent using the anova() function of R; if the reduced model was not rejected at significance threshold of 0.05, the process continued. The final model that cannot be reduced any further without a significant drop of the goodness of fit contains the following 5 variables:

epi.PC1’ ~ se.n.n2 + sn.n.n2 + ss.n.n2 + se.n1.n1 + aa.n1.n1 + 0

When trained on all 16 seasons (Figure 5), this model explains 0.98 of the original variance (F-statistics p-value of 3x10^-9), the coefficients being 1.26, 1.78, 0.81, -0.41 and 2.08, respectively. In the leave-half-out jackknife test the model explains, on average, 0.91 of the variance in the remaining 8 seasons.

Because Hong Kong Centre for Health Protection data lacks ILI raw numbers, only %ILI for each month is available. Therefore, we calculated the total epidemic severity (%ILI*%Positive Specimens) for each month, and then totaled these numbers for a year to obtain. (As a test to show that this method gives results similar to the regular method (using cumulative raw numbers for each season to calculate %ILI and %Positive Specimens), we checked  the correlation between the epidemic severity numbers calculated using the two methods for the USA data. The correlation between the two methods is very high (r=0.97). The H3 epidemic severity (S-H3) was calculated as the total influenza epidemic severity prorated by the fraction of H3 isolates among all isolates. These data were converted to a single epidemiological index using PCA (Figure 6). Interestingly, despite the fact that Hong Kong, like USA, is a non-tropical Northern Hemisphere country, the epidemic severity measures between the two are relatively weakly correlated (r_s[epi.PC1.USA,epi.PC1.HK] = 0.36 with the p-value of 0.21). Thus the Hong Kong dataset provides a largely independent test of the model applicability.

Sequence and antigenic distances were computed in a manner similar to that for the US data (the list of Hong Kong isolates used for genetic distance calculation is in the supplementary file [21]). Again, in this analysis, a season corresponds to a calendar year. Due to a low number of serologically tested isolates from Hong Kong proper, antigenic data were computed for a combined set of Hong Kong and mainland China isolates. The seasonal epidemiology seems to be the same in Hong Kong and Southern mainland China, and only Northern mainland China follows the standard Northern Hemisphere pattern [30]. Moreover, the Northern China epidemics follow the Southern China epidemics and not vice versa [30]. Based on this observation, two possible ways of merging these data sets were implemented, producing two versions of the serological distances (ac and as in contrast to aa and an serological distances for US). The first one (ac) corresponds to the scenario where for all isolates (Hong Kong and China) the season was considered to correspond to a calendar year, following the Hong Kong pattern. The second one (as) considers Hong Kong isolates following the Hong Kong pattern and China isolates following the Northern Hemisphere pattern with seasons starting in the fall.

There are several possible ways to apply the model derived from the US data to the Hong Kong data. One would be to directly compute the model prediction using the formula:

epi.PC1' = 1.26se.n.n2 + 1.78sn.n.n2 + 0.81ss.n.n2 - 0.41se.n1.n1 + 2.08aa.n1.n1

and replacing aa.n1.n1 (US data) with either ac.n1.n1 or as.n1.n1. Using as.n1.n1 for the serological distance variable, this formula produces a prediction that is reasonably well correlated with the real epidemiology data (r_s[epi.PC1,epi.PC1’] = 0.60, with the p-value of 0.012 in a permutation test; Figure 7, "retro.direct" plot). This model allows rough prediction of the ups and downs of the H3N2 influenza epidemics in 1997-2009, but gives a relatively poor quantitative estimate.

When the model is allowed to use the actual epidemiological data from Hong Kong to adjust its coefficients (yielding 1.50, 1.46, -0.82, 0.34 and 0.64 respectively), the prediction is improved (Figure 7, "retro.adjust" plot). The adjusted model explains 0.75 of the original variance with F-statistics p-value of 2x10^-2.

Finally, the stepwise reduction of the full model containing the genetic and serological distances from the corresponding seasons (se.n.n2, sn.n.n2, ss.n.n2, se.n1.n1, sn.n1.n1, ss.n1.n1, ac.n1.n1 and as.n1.n1) leads to the following 4-parameter model (Figure 7, "retro.stepwise" plot):

epi.PC1’ ~ se.n.n2 + sn.n.n2 + ac.n1.n1 + as.n1.n1 + 0

The model coefficients are 1.25, 0.75, 0.76 and 1.10 respectively; it explains 0.78 of the original variance with F-statistics p-value of 5x10-3. When trained on 8 out of the 13 seasons (leave-5-out jackknife test scheme) , this model explains 0.42 of the original variance on average.

For the purpose of this analysis, we considered Southern Hemisphere as comprising only non-tropical Southern countries because of a different influenza epidemic pattern in the Tropics [31]. Southern Hemisphere seasons were considered to be equivalent to calendar years, and correspond to the subsequent Northern Hemisphere seasons (e.g., 2003 Southern Hemisphere season corresponds to the 2003/04 Northern Hemisphere season). There were no Southern Hemisphere sequences from 2008, and we used the only January 2009 sequence as a substitute for 2008 because January is borderline between the two Southern seasons (2008 and 2009). List of Southern Hemisphere isolates used for genetic distance calculation is in the supplementary file [22]. Genetic and antigenic distances were computed as before.

We apply the second stage retrospective model (list of explanatory variables together with the coefficients) trained on the original dataset to this modified dataset to obtain the model projections to the current season using the Southern Hemisphere data from the preceding season (Figure 9).

As a demonstration of this approach, we divided the USA isolates from the 2002-2003 season into two subgroups - 9 Panama-like (members of the larger Sydney-like class) (supplementary file [32]) and 7 Fujian-like viruses (supplementary file [33]), based on the reconstructed phylogenetic tree for HA1 (A/H3N2 influenza) (supplementary file [34]) as described in [35].  Each of the subgroups was used to substitute for the real 2003-2004 season isolates in our data. The raw sequence distances computed for these subgroups were normalized using means and standard deviations previously determined for the full dataset. Then we used the list of variables and the coefficients determined for the second stage retrospective model to predict the severity of the 2003-2004 epidemics. The two subsets correspond to the two epidemiological hypotheses: one assumes that the 2003-2004 flu season will be dominated by Panama-like H3N2 virus, whereas the other assumes the Fujian-like virus prevalence. The results of these projections are shown on Figure 10.

Using the Panama projection predicts a very mild flu season for 2003-2004; in contrast, the Fujian projection predicts a severe epidemics. In reality the 2003-2004 flu season was dominated by the Fujian-like isolates worldwide and was the strongest on record for our range of dates in the USA.

Competition between Influenza A virus subtypes has been noted in several focused studies [45],[46],[47] as well as in a broader analysis of surveillance data [48]. However, the impact of competition from H1N1 on variation of variation of H3N2 seasonal severity is unlikely to have been significant during the time interval studied here because H3N2 was dominant most seasons, and furthermore, there is little evidence of positive selection for H1N1 antigenic novelty [49],[35]. In addition, Wolf et al.[35] have shown that during this period, H1N1 tended to dominate only after mild H3N2 seasons, so the competition appeared to be one-sided.

Although our results show that the degree of hemagglutinin novelty explains most of seasonal morbidity, other viral proteins have been shown to play important roles.  For example, consider the evolutionary history of the Fujian H3N2 strain:  initially an antigenically novel minor variant, succeeded by a reassortant containing the novel Fujian HA gene, which in turn is succeeded by a later "pure" Fujian variant in which all genes are derived from the original Fujian strain.  An analysis of these events led to a proposal that deleterious mutations in genes other than HA prevented the earliest members of this clade from dominating despite their antigenic novelty but subsequent compensatory mutations ultimately led to its dominance over the reassortant strain [35].  Recent experiments on the later Fujian variant support the idea that mutations in the polymerase PA gene were responsible for the decreased fitness of the original Fujian strain [50].  

Antigenic diversity in the preceding season is an important component of our statistical model which is consistent with the above observation that an antigenically novel strain may circulate for a season as a minor variant before compensating mutations allow it to dominate the following season.  Mutations in the hemagglutinin gene itself might have effects on fitness beyond that of changing antigenicity.  For example, mutations have been identified in the HA gene that might allow the virus to evade the host's immune response but also modify binding to the host receptor leading to decreased replication fitness (cf. [51],[52],[53]).  Our results suggest that this could be a common feature of early drift variants.

Antigenic diversity in the preceding season is not the only warning of a rise in morbidity but it is difficult to extract other clear signals or rules from our analysis.  In addition to the underlying complexity of the system, the sequence and HI data we have for clinical isolates are not sufficiently comprehensive or representative.  For example, in the 2002-2003 season, according to the CDC's antigenic characterization of the 2002 - 03 U.S. Influenza season [16], 85% of the H3N2 clinical isolates were similar to A/Panama/2007/99 and 15% were similar to A/Fujian/411/2002. However, our sequence data were about half Panama and half Fujian whereas virtually all of the HI data we were able to obtain were similar to Fujian.

While more representative input data could improve the performance of our statistical model, refinements in other aspects of the method may also be useful.  For example, other sequence-derived measures of novelty might work better than our simple amino acid replacements-based approach.  The method of Gupta et al. [12] yielded equivalent results but there are additional sequence-based methods that should be evaluated as well (cf.[54]).

Because our method provides an accurate reconstruction of epidemiological severity given sequence and serological data, one can evaluate the epidemiological consequences of various assumptions about the population of viruses for the upcoming season.  As seen in Figure 9, projections using data from Southern Hemisphere isolates provide realistic estimates of severity for the Northern Hemisphere.  When there is significant diversity among co-circulating clades in a given season, one can assess different scenarios for the upcoming season, as was done using data from 2002-2003 seasonal isolates to project severity for the 2003-2004 season.  The panel advising the US Food and Drug Administration on vaccines had difficulties that season because the Fujian strain - a minor strain in 2002-2003 - was antigenically novel but was difficult to grow in chicken eggs.   Scientists at the CDC influenza branch had been able to passage the Fujian strain through dog kidney cells and obtain egg-adapted virus but concern about contamination steered the FDA advisory committee towards a decision to use the Panama strain ( cf [55]).  The ability to make accurate projections of epidemiological severity might have been helpful to further inform this decision.

Although, in some cases, projections based on co-circulating clades could be helpful prior to the availability of Southern Hemisphere isolates, thus assisting in vaccine strain decisions, in other situations, such as the 1997-1998 Sydney season, the rapid emergence of a new dominant strain makes earlier projections effectively useless.  Deeper sampling of representative viral isolates might provide an earlier warning of such novel clades, and more sophisticated models using additional data might also eventually prove helpful.

Finally, we note that our predictions of influenza morbidity based on sequence and serological data are limited to the inter-pandemic period, and that our model cannot be used to project the severity of an emerging pandemic influenza virus, such as 2009 H1N1pdm. In addition, it remains unclear whether this model rooted in empirical data collected before the 2009 pandemic will be suitable to predict the severity of influenza A/H3N2 season in the post-pandemic period. Both H3N2 and H1N1 have begun to co-circulate in 2010 [56], and further studies will need to determine whether the dynamics of these two subtypes has changed and whether the algorithms proposed here need to be revised to accommodate these changes.