Antigenic drift and subtype interference shape A(H3N2) epidemic dynamics in the United States

Influenza viruses continually evolve new antigenic variants, through mutations in epitopes of their major surface proteins, hemagglutinin (HA) and neuraminidase (NA). Antigenic drift potentiates the reinfection of previously infected individuals, but the contribution of this process to variability in annual epidemics is not well understood. Here we link influenza A(H3N2) virus evolution to regional epidemic dynamics in the United States during 1997—2019. We integrate phenotypic measures of HA antigenic drift and sequence-based measures of HA and NA fitness to infer antigenic and genetic distances between viruses circulating in successive seasons. We estimate the magnitude, severity, timing, transmission rate, age-specific patterns, and subtype dominance of each regional outbreak and find that genetic distance based on broad sets of epitope sites is the strongest evolutionary predictor of A(H3N2) virus epidemiology. Increased HA and NA epitope distance between seasons correlates with larger, more intense epidemics, higher transmission, greater A(H3N2) subtype dominance, and a greater proportion of cases in adults relative to children, consistent with increased population susceptibility. Based on random forest models, A(H1N1) incidence impacts A(H3N2) epidemics to a greater extent than viral evolution, suggesting that subtype interference is a major driver of influenza A virus infection dynamics, presumably via heterosubtypic cross-immunity.


Introduction
(linear models, LMs: H3, adjusted R 2 = 0.37, P = 0.03; N2: R 2 = 0.26, P = 0.08) and peak incidence (LMs,188 H3: R 2 = 0.4, P = 0.02; N2: R 2 = 0.33, P = 0.04) and higher effective Rt (generalized linear models, GLMs: 189 H3, R 2 = 0.38, P = 0.05; N2, R 2 = 0.32, P = 0.03) (regression results: Figure 3, Spearman's correlations: 190 Figure S7). HI log2 titer distance (t -2) exhibited positive but non-significant associations with different 191 measures of epidemic impact (Figure 3, Figure S7). Seasonal diversity in the growth rates of circulating 192 lineages in the current t or prior season (t -1) had strong negative correlations with effective Rt (GLMs,193 H3 (t -1): R 2 = 0.49, P = 0.009; N2, t: R 2 = 0.46, P = 0.006) and epidemic intensity (Beta GLMs, 1): R 2 = 0.45, P = 0.003; N2, t: R 2 = 0.51, P = 0.001) ( Figures S7-S8). Seasonal mean LBI exhibited 195 similar but slightly weaker correlations with effective Rt and epidemic intensity. Pneumonia and influenza 196 excess mortality attributable to A(H3N2) also increased with H3 epitope distance, though this relationship 197 was not statistically significant ( Figure S9). The remaining indicators of viral evolution, including H3 and 198 N2 non-epitope distance (mutational load), H3 RBS distance, and H3 stalk footprint distance had weak, 199 non-significant correlations with the different measures of epidemic impact ( Figure S7). 200 We explored whether evolutionary changes in A(H3N2) may predispose this subtype to dominate 201 influenza virus circulation in a given season. A(H3N2) subtype dominance -the proportion of influenza 202 positive samples typed as A(H3N2) -increased with H3 epitope distance (t -2) and N2 epitope distance 203 (t -1) (Beta GLMs, H3: R 2 = 0.32, P = 0.05; N2: R 2 = 0.34, P = 0.03; Figure 4, Figure S7). Figure 4 204 illustrates this relationship at the regional level across two seasons in which A(H3N2) was nationally 205 dominant, but where antigenic change differed. In 2003-2004, we observed widespread dominance of 7 diversity in the current and prior season, due to weaker or non-significant correlations between the other 285 evolutionary metrics and epidemic burden ( Figure S7). To account for potential type or subtype 286 interference, we included A(H1N1) epidemic size (A(H1N1) or A(H1N1)pdm09) and B epidemic size in the 287 current and prior season and the dominant IAV subtype in the prior season. We included A(H3N2) 288 epidemic size in the prior season as a proxy of natural prior immunity to A(H3N2). To account for vaccine-289 induced immunity, we considered four categories of predictors and included estimates for the current and between currently circulating strains and the US vaccine reference strain. We could not include a 294 predictor for vaccination coverage in children or consider clade-specific VE estimates, because data were 295 not available for most seasons in our study. We did not predict excess mortality attributable to A(H3N2), 296 due to data limitation (one national estimate per season) and omitted models predicting epidemic timing, 297 due to weak or non-significant correlations between timing-related measures and most indicators of viral 298 evolution ( Figure S11). Lastly, we could not separate our analysis into pre-and post-2009 pandemic 299 periods due to small sample sizes. 300 Based on variable importance scores, A(H1N1) epidemic size in the current season was the most 301 informative predictor of A(H3N2) epidemic size and peak incidence, followed by H3 epitope distance, and 302 the dominant IAV subtype in the previous season or N2 epitope distance ( Figure 6). For A(H3N2) subtype 303 dominance, the highest ranked predictors were H3 epitope distance, N2 epitope distance, and the 304 dominant IAV subtype in the previous season ( Figure 6). We note that we did not include A(H1N1) 305 epidemic size as a predictor in this model, due to its confounding with the target variable. For models of 306 A(H3N2) effective Rt and epidemic intensity, we observed less discernable differences in variable 307 importance scores across the set of candidate predictors ( Figure 6). For the model of effective Rt, N2 LBI 308 diversity in the current season, A(H1N1) epidemic size in the current season, and N2 epitope distance 309 between circulating strains and the vaccine strain were the highest ranked variables, while the most 310 important predictors of epidemic intensity were H3 and N2 LBI diversity in the current season and adult 311 vaccination coverage in the current and prior season. Variable importance rankings from LASSO (least 312 absolute shrinkage and selection operator) regression models were qualitatively similar to those from 313 random forest models, with A(H1N1) epidemic size in the current season, H3 and N2 epitope distance, 314 and the dominant IAV subtype in the prior season consistently retained across the best-tuned models of 315 epidemic size, peak incidence, and subtype dominance ( Figure S20). Vaccine-related parameters and H3 316 antigenic drift (either H3 epitope distance or HI log2 titer distance) were retained in the best-tuned LASSO 317 models of effective Rt and epidemic intensity ( Figure S20). 318 We measured correlations between observed values and model-predicted values at the HHS region level. 319 Among our various epidemic metrics, random forest models produced the most accurate predictions of 320 A(H3N2) subtype dominance ( = 0.94, regional range = 0.8 -0.98), peak incidence (Spearman's = 321 0.91, regional range = 0.73 -0.95), and epidemic size ( = 0.9, regional range = 0.73 -0.94), while 322 predictions of effective Rt and epidemic intensity were less accurate ( = 0.8, regional range = 0.65 -323 0.91; = 0.78, regional range = 0.63 -0.91, respectively) ( Figure 7). Random forest models tended to 324 underpredict most epidemic targets in seasons with substantial H3 antigenic transitions, in particular the 325 SY97 cluster seasons (1998)(1999)(1999)(2000) and the FU02 cluster season (2003)(2004) (Figure 7). 326 For epidemic size and peak incidence, seasonal predictive error -root-mean-square error (RMSE) across 327 all regional predictions in a season -increased with H3 epitope distance (size, Spearman's = 0.51, P = 328 0.02; peak, = 0.61, P = 0.007) and N2 epitope distance (size, = 0.43, P = 0.06; peak, = 0.46, P = 329 0.04). For models of epidemic intensity, seasonal RMSE increased with N2 epitope distance ( = 0.62, P 330 = 0.006) but not H3 epitope distance ( = 0.07, P = 0.8) (Figures S21-S22 This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted October 3, 2023. ;https://doi.org/10.1101https://doi.org/10. /2023 To further refine our set of informative predictors, we performed multivariable regression with the top 10 334 ranked predictors from each random forest model and used Bayesian Information Criterion (BIC) to select 335 the best fit model for each epidemic metric, allowing each metric's regression model to include up to three 336 independent variables. This additional step of variable selection demonstrated that models with few 337 predictors fit the observed data relatively well (epidemic size, adjusted R 2 = 0.69; peak incidence, adj. R 2 338 = 0.63; effective Rt, adj. R 2 = 0.65; epidemic intensity, adj. R 2 = 0.75), except for subtype dominance (adj.

339
R 2 = 0.48) ( Table 3). The set of variables retained after model selection were similar to those with high 340 importance rankings in random forest models and LASSO regression models, with the exception that HI 341 log2 titer distance, rather than H3 epitope distance, was included in the minimal models of effective Rt and 342 epidemic intensity.

345
Antigenic drift between currently circulating influenza viruses and the previous season's viruses is 346 expected to confer increased viral fitness, leading to earlier, larger, or more severe epidemics. However, 347 prior evidence for the impact of antigenic drift on seasonal influenza outbreaks is mixed. Here, we 348 systematically compare experimental and sequence-based measures of A(H3N2) evolution in predicting 349 regional epidemic dynamics in the United States across 22 seasons, from 1997 to 2019. We also 350 consider the effects of other co-circulating influenza viruses, prior immunity, and vaccine-related 351 parameters, such as coverage and effectiveness, on A(H3N2) incidence. Our findings indicate that 352 evolution in both major surface proteins -hemagglutinin (HA) and neuraminidase (NA) -contributes to 353 variability in epidemic magnitude across seasons, though viral fitness appears to be secondary to subtype 354 interference in shaping annual outbreaks.

356
The first question of this study sought to determine which metrics of viral fitness have the strongest 357 relationships with A(H3N2) epidemic burden and timing. Among our set of candidate evolutionary 358 predictors, genetic distances based on broad sets of epitope sites (HA = 129 sites; NA = 223 epitope 359 sites) had the strongest, most consistent associations with A(H3N2) epidemic size, transmission rate, 360 severity, subtype dominance, and age-specific patterns. Increased epitope distance in both H3 and N2 361 correlated with larger epidemics and increased transmissibility, with univariate analyses finding H3 362 distance more strongly correlated with epidemic size, peak incidence, transmissibility, and excess 363 mortality, and N2 distance more strongly correlated with epidemic intensity (i.e., the "sharpness" of the 364 epidemic curve) and subtype dominance patterns. However, we note that minor differences in correlative 365 strength between H3 and N2 epitope distance are not necessarily biologically relevant and could be 366 attributed to noise in epidemiological or virological data or the limited number of influenza seasons in our 367 study. The fraction of ILI cases in children relative to adults was negatively correlated with N2 epitope 368 distance, consistent with the expectation that cases are more restricted to immunologically naïve children 369 in seasons with low antigenic novelty [7,78]. Regarding epidemic timing, the number of days from 370 epidemic onset to peak (a proxy for epidemic speed) decreased with N2 epitope distance, but other Positive associations between H3 antigenic drift and population-level epidemic burden are consistent with 385 previous observations from theoretical models [25,26,83]. For example, phylodynamic models of 386 punctuated antigenic evolution have reproduced key features of A(H3N2) phylogenetic patterns and case 387 dynamics, such as the sequential replacement of antigenic clusters, the limited standing diversity in HA 388 after a cluster transition, and higher incidence and attack rates in cluster transition years [25,26,83]. Our 389 results also corroborate empirical analyses of surveillance data [6,27,28,66] and forecasting models of 390 annual epidemics [29,30] that found direct, quantitative links between HA antigenic novelty and the 391 number of influenza cases or deaths in a season. Moving beyond the paradigm of antigenic clusters, Wolf 392 et al., Bedford et al., 2014 demonstrated that smaller, year-to-year changes in H3 antigenic drift 393 also correlate with seasonal severity and incidence [6,27]. A more recent study did not detect an 394 association between antigenic drift and city-level epidemic size in Australia [31], though the authors used 395 a binary indicator to signify seasons with major HA antigenic transitions and did not consider smaller, 396 more gradual changes in antigenicity. While Lam and colleagues did not observe a consistent effect of 397 antigenic change on epidemic magnitude, they found a negative relationship between the cumulative prior 398 incidence of an antigenic variant and its probability of successful epidemic initiation in a city [31].

400
We did not observe a clear relationship between H3 receptor binding site (RBS) distance and epidemic 401 burden, even though single substitutions at these seven amino acid positions are implicated in major 402 antigenic transitions [68,84]. The outperformance of the RBS distance metric by a broader set of epitope 403 sites could be attributed to the tempo of antigenic cluster changes. A(H3N2) viruses are characterized by 404 both continuous and punctuated antigenic evolution, with transitions between antigenic clusters occurring 405 every 2 to 8 years [6,26,32,33,36,37,67,68,85] predictions. First, epitope distances may capture properties that affect viral fitness (and in turn outbreak 442 intensity) but are unrelated to immune escape, such as intrinsic transmissibility, ability to replicate, or 443 epistatic interactions. A second set of factors concern methodological issues associated with HI assays.

444
The reference anti-sera for HI assays are routinely produced in ferrets recovering from their first influenza   This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted October 3, 2023. ; https://doi.org/10.1101/2023.10.02.23296453 doi: medRxiv preprint considered the effects of influenza A(H1N1) incidence and B incidence on A(H3N2) virus circulation 491 within a season. We detected strong negative associations between A(H1N1) incidence and A(H3N2) 492 epidemic size, peak incidence, transmissibility, and excess mortality, consistent with previous animal, 493 epidemiological, phylodynamic, and theoretical studies that found evidence for cross-immunity between 494 IAV subtypes [53][54][55]57,59,102]. For example, individuals recently infected with seasonal influenza A 495 viruses are less likely to become infected during subsequent pandemic waves [52,53,[102][103][104], and the 496 early circulation of one influenza virus type or subtype is associated with a reduced total incidence of the 497 other type/subtypes within a season [31,57]. Due to the shared evolutionary history of their internal genes In our study, univariate correlations between A(H1N1) and A(H3N2) incidence were more pronounced 506 than those observed between A(H3N2) incidence and evolutionary indicators, and A(H1N1) epidemic size 507 was the highest ranked feature by random forest models predicting epidemic size and peak incidence.

508
Consequently, interference between the two influenza A subtypes may be more impactful than viral 509 evolution in determining the size of annual A(H3N2) outbreaks. Concerning epidemic timing, we did not 510 detect a relationship between A(H3N2) antigenic change and the relative timing of A(H3N2) and A(H1N1) 511 cases; specifically, A(H3N2) incidence did not consistently lead A(H1N1) incidence in seasons with 512 greater H3 or N2 antigenic change. Overall, we did not find any indication that influenza B incidence 513 affects A(H3N2) epidemic burden or timing, which is not unexpected, given that few T and B cell epitopes 514 are shared between the two virus types [106].

516
Lastly, we used random forest models and multivariable linear regression models to assess the relative 517 importance of viral evolution, prior population immunity, co-circulation of other influenza viruses, and 518 vaccine-related parameters in predicting regional A(H3N2) epidemic dynamics. We chose conditional 519 inference random forest models as our primary method of variable selection because several covariates 520 were collinear, relationships between some predictors and target variables were nonlinear, and our goal 521 was inferential rather than predictive. We performed leave-one-season-out cross-validation to tune each 522 model, but, due to the limited number of seasons in our dataset, we were not able to test predictive 523 performance on an independent test set. With the caveat that models were likely overfit to historical data, 524 random forest models produced accurate predictions of regional epidemic size, peak incidence, and 525 subtype dominance patterns, while predictions of epidemic intensity and transmission rates were less 526 exact. The latter two measures could be more closely tied to climatic factors, the timing of influenza case 527 importations from abroad, or mobility patterns [7,13,14,16] or they may be inherently more difficult to 528 predict because their values are more constrained. Random forest models tended to underpredict 529 epidemic burden in seasons with major antigenic transitions, particularly the SY97 seasons (1998-1999, 530 1999-2000) and the FU02 season (2003)(2004), potentially because antigenic jumps of these magnitudes 531 were infrequent during our 22-season study period. An additional step of post-hoc model selection 532 demonstrated that models with only three covariates could also produce accurate fits to observed 533 epidemiological data.

535
Our study is subject to several limitations, specifically regarding geographic resolution and data 536 availability. First, our analysis is limited to one country with a temperate climate and its findings 537 concerning interactions between A(H3N2), A(H1N1), and type B viruses may not be applicable to tropical 538 or subtropical countries, which experience sporadic epidemics of all three viruses throughout the year 539 [107]. Second, our measure of population-level influenza incidence is derived from regional CDC This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted October 3, 2023. ; https://doi.org/10.1101/2023.10.02.23296453 doi: medRxiv preprint electronic health records are accessible in theory but not in the public domain. Access to ILI cases 543 aggregated at the state or city level, collected over the course of decades, would increase statistical 544 power and enable us to add more location-specific variables to our analysis, such as climatic and 545 environmental factors. A third limitation is that we measured influenza incidence by multiplying the rate of 546 influenza-like illness by the percentage of tests positive for influenza, which does not completely eliminate 547 the possibility of capturing the activity of other co-circulating respiratory pathogens [11]. Surveillance data 548 based on more specific diagnosis codes would ensure the exclusion of patients with non-influenza 549 respiratory conditions. Fourth, our data on the age distribution of influenza cases were derived from ILI 550 encounters across four broad age groups and did not include test positivity status, virus type/subtype, or 551 denominator information. Despite the coarseness of these data, we found statistically significant 552 correlations in the expected directions between N2 antigenic change and the fraction of cases in children 553 relative to adults. Lastly, a serological assay exists for NA, but NA titer measurements are not widely 554 available because the assay is labor-intensive and inter-lab variability is high. Thus, we could not test the 555 performance of NA antigenic phenotype in predicting epidemic dynamics.  In conclusion, relationships between A(H3N2) antigenic drift, epidemic impact, and age dynamics are 572 moderate, with genetic distances based on broad sets of H3 and N2 epitope sites having greater 573 predictive power than serology-based antigenic distances for the timeframe analyzed. Influenza 574 epidemiological patterns are consistent with increased population susceptibility in seasons with high 575 antigenic novelty, and our study is the first to link NA antigenic drift to epidemic burden, timing, and the 576 age distribution of cases. It is well established that anti-HA and anti-NA antibodies are independent 577 correlates of immunity [45,[47][48][49][116][117][118], and the influenza research community has advocated for NA-578 based vaccines [39,119]. The connection between NA drift and seasonal incidence further highlights the 579 importance of monitoring evolution in both HA and NA to inform vaccine strain selection and epidemic 580 forecasting efforts. Although antigenic change in both HA and NA was correlated with epidemic dynamics, 581 ecological interactions between influenza A subtypes appear to be more influential than viral evolution in 582 determining the intensity of annual A(H3N2) epidemics. The aim of our study was to retrospectively 583 assess the potential drivers of annual A(H3N2) epidemics, yet we cautiously suggest that one could 584 project the size or intensity of future epidemics based on sequence data and A(H1N1)pdm09 incidence This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. week 16 due to the emergence of the A(H1N1)pdm09 virus [57]. 595 We extracted syndromic surveillance data for the ten HHS regions from the U.S. Outpatient Influenza-like 596 Illness Surveillance Network (ILINet) [120]. ILINet consists of approximately 3,200 sentinel outpatient 597 healthcare providers throughout the United States that report the total number of consultations for any 598 reason and the number of consultations for influenza-like illness (ILI) every week. ILI is defined as fever 599 (temperature of 100°F [37.8°C] or greater) and a cough and/or a sore throat. The indicator is based on 600 the weekly proportion of outpatient consultations for influenza-like illness and is available weighted or 601 unweighted by regional population size. The number of ILI encounters by age group are also provided (0-602 4, 5-24, 25-64, and ≥65), but these data are not weighted by total encounters or population size. 33. For each region, we scaled pre-pandemic incidences by the ratio of mean weekly ILI + (for all influenza 623 type/subtypes combined) in the post-pandemic period to that of the pre-pandemic period. Incidences for 624 HHS Region 10 were not adjusted for pre-and post-pandemic reporting because surveillance data for this 625 region were not available before 2009. To account for differences in reporting rates across HHS regions, 626 we next scaled each region's type/subtype incidences by its mean weekly ILI + for the entire study period.

627
Scaled incidences were used in all downstream analyses of epidemic burden and timing. 629 Epidemic burden: We considered three complementary indicators of epidemic burden, separately for 630 each influenza type/subtype, HHS region, and season. We defined peak incidence as the maximum 631 weekly scaled incidence and epidemic size as the cumulative weekly scaled incidence. We also 632 estimated epidemic intensity based on a method previously developed to study variation in the shape 633 (i.e., sharpness) of influenza epidemics across US cities [123]. Epidemic intensity was based on the 634 inverse Shannon entropy of the weekly incidence distribution. Epidemic intensity increases when 635 incidence is more concentrated in particular weeks and decreases when incidence is more evenly spread This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

Epidemic burden and timing
The copyright holder for this preprint this version posted October 3, 2023.

696
The 95% CI for the weighted mean VE was calculated as: Correlations among epidemic metrics 699 We used Spearman's correlation coefficients to measure pairwise relationships between A(H3N2) 700 epidemiological indictors. We adjusted P-values for multiple testing using the Benjamini and Hochberg 701 method [168]. 703 We considered multiple indicators of influenza evolution based on genetic and phenotypic (serologic)

Indicators of influenza A(H3N2) evolution
wherein pi is the proportion of LBI values belonging to the ith bin.

762
Antigenic and genetic distance relative to prior seasons 763 We estimated genetic and antigenic distances between influenza viruses circulating in consecutive 764 seasons by calculating the mean distance between viruses circulating in the current season t and viruses 765 circulating during the prior season (t -1 year; one season lag) or two prior seasons ago (t -2 years; two 766 season lag). Seasonal genetic and antigenic distances are greater when currently circulating strains are 767 more antigenically distinct from previously circulating strains. We used Spearman's correlation 768 coefficients to measure pairwise relationships between scaled H3 and N2 evolutionary indictors. We We measured univariate associations between national indicators of A(H3N2) viral fitness and regional 773 A(H3N2) epidemic parameters -peak incidence, epidemic size, effective Rt, epidemic intensity, subtype 774 dominance, excess P&I deaths, onset timing, peak timing, spatiotemporal synchrony, the number of 775 weeks from onset to peak, and seasonal duration. We first measured Spearman correlation coefficients 776 between pairs of scaled fitness indicators and epidemic metrics using 1000 bootstrap replicates of the 777 original dataset (1000 samples with replacement).

778
Next, we fit regression models with different distribution families (Gaussian or Gamma) and link functions 779 (identity, log, or inverse) to observed data and used Bayesian information criterion (BIC) to select the best 780 fit model, with lower BIC values indicating a better fit to the data. For subtype dominance, epidemic 781 intensity, and age-specific proportions of ILI cases, we fit Beta regression models with logit links. For 782 each epidemic metric, we fit the best-performing regression model to the resampled dataset. To measure 783 the effects of sub(type) interference on A(H3N2) epidemics, the same approach was applied to measure 784 the univariate relationships between A(H1N1) or B epidemic size and A(H3N2) peak incidence, epidemic 785 size, effective Rt, epidemic intensity, and excess mortality. As a sensitivity analysis, we tested univariate Next, we explored multivariable approaches that would shed light on the potential mechanisms driving 794 annual epidemic impact. Considering that we had many predictors and relatively few observations (22 795 seasons x 9-10 HHS regions), several covariates were collinear, and our goal was explicative rather than 796 predictive, we settled on methods that tend to select few covariates.

798
We first used conditional inference random forest models to select relevant predictors of A(H3N2) 799 epidemic size, peak incidence, effective Rt, epidemic intensity, and subtype dominance (party and caret This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted October 3, 2023. ; https://doi.org/10.1101/2023.10.02.23296453 doi: medRxiv preprint the current and previous season, A(H3N2) vaccine effectiveness in the current and previous season, and 806 H3 and N2 epitope distances between circulating A(H3N2) viruses in the United States and the A(H3N2) 807 vaccine strain in the same season. We did not conduct variable selection analysis for excess A(H3N2) 808 mortality due to data limitations (one national estimate per season). Metrics related to epidemic timing 809 were also excluded from this analysis because we found weak or non-statistically significant associations 810 with most of the candidate evolutionary predictors in univariate analyses.

812
We created each forest by generating 3,000 regression trees from 10 repeats of a leave-one-season-out 813 (jackknife) cross-validated sample of the data. Due to the small size of our dataset, evaluating the 814 predictive accuracy of random forest models on a quasi-independent test set produced unstable 815 estimates. Consequently, we included all data in the training set and report root mean squared error 816 (RMSE) and R 2 values from the best tuned model. We used permutation importance (N = 50 817 permutations) to estimate the relative importance of each predictor in determining model outcomes.

818
Permutation importance is the decrease in prediction accuracy when a single feature (predictor) is 819 randomly permuted, with larger values indicating more important variables. Because our features were 820 collinear, we used conditional permutation importance to compute feature importance scores, rather than 821 the standard marginal procedure [177,178,180,181].

823
As an alternative method for variable selection, we performed LASSO regression on the same cross-824 validated dataset and report RMSE and R 2 values from the best tuned model (glmnet and caret R 825 packages) [179,182]. Unlike random forest models, this approach assumes linear relationships between 826 predictors and the target variable. LASSO models (L1 penalty) are more restrictive than ridge models (L2 827 penalty) and elastic net models (combination of L1 and L2 penalties) and will arbitrarily select one 828 variable from a set of collinear variables.

830
To further reduce the set of predictors for each epidemic metric, we performed model selection with linear 831 regression models that considered all combinations of the top 10 ranked predictors from conditional 832 inference random forest models. Candidate models were limited to three independent variables, and 833 models were compared using BIC. We did not include HHS region or season as fixed or random effects in 834 these models because these variables either did not improve model fit (region) or caused convergence 835 issues (season).

837
All predictors were centered and scaled prior to fitting random forest or regression models.    This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.  This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.    This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted October 3, 2023. ; https://doi.org/10.1101/2023.10.02.23296453 doi: medRxiv preprint Figure 6. Variable importance rankings from conditional inference random forest models 1499 predicting A(H3N2) epidemic dynamics. Ranking of variables in predicting regional A(H3N2) A. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.  This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted October 3, 2023. ; https://doi.org/10.1101/2023.10.02.23296453 doi: medRxiv preprint 1 Candidate models were limited to 3 independent variables and considered all combinations of the top 10 1535 ranked predictors from conditional inference random forest models ( Figure 6).

1536
for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted October 3, 2023.  1997-1998, 1998-1999, 1999-2000) and the FU02 cluster season (2003)(2004). This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted October 3, 2023. ; https://doi.org/10.1101/2023.10.02.23296453 doi: medRxiv preprint Figure S2. Pairwise correlations between H3 and N2 evolutionary indicators (one season lags). We 1552 measured Spearman's correlations between seasonal measures of H3 and N2 evolution, including H3 1553 RBS distance, H3 epitope distance, H3 non-epitope distance, H3 stalk footprint distance, HI titer distance, 1554 N2 epitope distance based on 223 or 53 epitope sites, N2 non-epitope distance, mean clade growth of H3  This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.

1567
Seasonal distances were estimated as the mean distance between strains circulating in the current 1568 season t and those circulating in the prior season (t -1). The Benjamini and Hochberg method was used  This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted October 3, 2023. ;https://doi.org/10.1101https://doi.org/10. /2023  This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted October 3, 2023. ;https://doi.org/10.1101https://doi.org/10. /2023  This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted October 3, 2023. ;https://doi.org/10.1101https://doi.org/10. /2023  This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted October 3, 2023. ;https://doi.org/10.1101https://doi.org/10. /2023  This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted October 3, 2023. ; https://doi.org/10.1101/2023.10.02.23296453 doi: medRxiv preprint Figure S10. Regional patterns of influenza type and subtype incidence from seasons 1997-1998 to This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.  This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted October 3, 2023. ; https://doi.org/10.1101/2023.10.02.23296453 doi: medRxiv preprint Figure S13. Epidemic speed increases with N2 antigenic drift. N2 epitope distance correlates with 1660 fewer days from epidemic onset to peak (A), while the relationship between H3 epitope distance and 1661 epidemic speed is less apparent (B). Seasonal epitope distance is the mean distance between strains 1662 circulating in season t and those circulating in the prior season (t -1) or two seasons ago (t -2). This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.  This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.  orange: A(H3N2)/A(H1N1)pdm09 co-dominant), and vertical bars are 95% confidence intervals of 1695 regional age distribution estimates. The fraction of cases in each age group were fit as a function of 1696 seasonal H3 or N2 epitope distance using Beta GLMs (logit link) with 1000 bootstrap resamples.

1697
for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.  This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted October 3, 2023. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted October 3, 2023. ; https://doi.org/10.1101/2023.10.02.23296453 doi: medRxiv preprint Figure S19. Wavelet analysis of influenza A and B epidemic timing. A. A(H3N2) incidence precedes 1719 A(H1N1) incidence in most seasons. Although A(H1N1) incidence sometimes leads or is in phase with 1720 A(H3N2) incidence (negative or zero phase lag), the direction of seasonal phase lags is not clearly This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted October 3, 2023. ; https://doi.org/10.1101/2023.10.02.23296453 doi: medRxiv preprint Figure S20. Variable importance rankings from LASSO models predicting A(H3N2) epidemic 1733 dynamics. Ranking of variables in predicting seasonal A(H3N2) A. epidemic size, B. peak incidence, C. 1734 transmissibility (effective reproduction number, Rt), D. epidemic intensity (inverse Shannon entropy), and 1735 E. subtype dominance. Models were tuned using a repeated leave-one-season-out cross-validated 1736 sample of the data. Variables are ranked by their coefficient estimates, with differences in prediction 1737 accuracy scaled by the total (null model) error. Abbreviations: HI titer = hemagglutination inhibition log2 1738 titer distance, t -1 = one-season lag, t -2 = two-season lag, LBI = local branching index, peak = peak 1739 incidence, distance to vaccine = epitope distance between currently circulating strains and the 1740 recommended vaccine strain, VE = vaccine effectiveness.

1741
for use under a CC0 license. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted October 3, 2023. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint this version posted October 3, 2023. This article is a US Government work. It is not subject to copyright under 17 USC 105 and is also made available (which was not certified by peer review) is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity.