Dynamic interactions of influenza viruses in Hong Kong during 1998-2018

Influenza epidemics cause substantial morbidity and mortality every year worldwide. Currently, two influenza A subtypes, A(H1N1) and A(H3N2), and type B viruses co-circulate in humans and infection with one type/subtype could provide cross-protection against the others. However, it remains unclear how such ecologic competition via cross-immunity and antigenic mutations that allow immune escape impact influenza epidemic dynamics at the population level. Here we develop a comprehensive model-inference system and apply it to study the evolutionary and epidemiological dynamics of the three influenza types/subtypes in Hong Kong, a city of global public health significance for influenza epidemic and pandemic control. Utilizing long-term influenza surveillance data since 1998, we are able to estimate the strength of cross-immunity between each virus-pairs, the timing and frequency of punctuated changes in population immunity in response to antigenic mutations in influenza viruses, and key epidemiological parameters over the last 20 years including the 2009 pandemic. We find evidence of cross-immunity in all types/subtypes, with strongest cross-immunity from A(H1N1) against A(H3N2). Our results also suggest that A(H3N2) may undergo antigenic mutations in both summers and winters and thus monitoring the virus in both seasons may be important for vaccine development. Overall, our study reveals intricate epidemiological interactions and underscores the importance of simultaneous monitoring of population immunity, incidence rates, and viral genetic and antigenic changes.


44
Influenza epidemics recur annually in many regions of the world and cause significant 45 morbidity and mortality, leading to 3 to 5 million severe infections and 291,000 to 646,000 46 deaths worldwide each year (1, 2). These recurrent epidemics are a combined outcome of 47 viral antigenic changes, interactions among co-circulating influenza viruses, transmission 48 and host immune response. Influenza viruses undergo continual genetic mutations and 49 punctuated antigenic changes, which allow them to escape prior immunity and cause new 50 epidemics (3-6). In addition, currently two influenza A subtypes, A(H1N1) and A(H3N2), 51 and influenza B co-circulate in humans; and infection of a certain type/subtype can result 52 in partial immunity to influenza strains of the same type/subtype as well as strains across 53 subtypes (7-15) and types (16-18). The latter cross-reactive immunity-termed cross-54 immunity-leads to epidemiological interactions among influenza viruses that, along with 55 viral antigenic changes, shape influenza phylodynamics and epidemic dynamics (19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30)(31). 56 Improving understanding of interactions among influenza viruses and the resulting impact 57 on epidemic dynamics will thus provide insights to aid public health efforts to mitigate the 58 burden of influenza. 59 60 As observed epidemic dynamics are the result of the aforementioned key 61 evolutionary and epidemiological processes, in this study, we utilized a long-term influenza 62 incidence surveillance dataset collected in Hong Kong since 1998 to estimate the timing of 63 punctuated antigenic changes, strength of cross-immunity, duration of immunity, as well as 64 other key epidemiological parameters for each of the three types/subtypes. Hong Kong is a 65 subtropical city located in Southeast Asia. Unlike the highly seasonal epidemics in 66 system was able to accurately estimate the infection rates during both epidemic and non-136 epidemic periods for all types/subtypes. In particular, over the 20-year study period, 137 epidemics of A(H1N1) were less frequent compared to the other two types/subtypes. For 138 instance, few A(H1N1) cases were reported during Jan 1998 -Dec 1999 (the entire 2-year 139 period) and May 2002 -Dec 2004 (the entire 2.5-year period). While R0 followed the same 140 seasonal cycle as years with epidemics (Fig 3A), the estimated infection rates remained at 141 near-zero over these extended time-periods when A(H1N1) was not in circulation (Fig 2A). 142

143
In addition to running the model continuously without re-initiation to correct for 144 model errors, we included a continuous seeding at a constant rate of 1 case per 100,000 145 population per 10 days for all types/subtypes. This epidemic seeding represents travel-146 related importation of cases into local population, since Hong Kong is a transportation hub 147 highly connected to places worldwide. This seeding is also designed to test the accuracy of 148 estimates of key model state variables and parameters, as they are intricately linked to the 149 infection rate. Foremost, overestimation of population susceptibility would lead to 150 overestimation of cases during non-epidemic periods, as the seeded cases would lead to an 151 epidemic should the effective reproductive number Re (computed as R0 times population 152 susceptibility) be above unity and sustain so. Conversely, underestimation of susceptibility 153 would lead to underestimation of cases during epidemic periods or require unrealistically 154 frequent antigenic changes prior to an epidemic. Further, population susceptibility is also 155 linked to the duration of immunity and the strength of cross-immunity. Loss of immunity 156 (determined by the immunity period) replenishes the susceptible pool whereas gaining of whereas scenario 3 (i.e., no cross-immunity) largely overestimated the frequency of co-228 epidemics caused by multiple types/subtypes. 229 230 Timing and frequency of punctuated antigenic changes. Analyses of the antigenic 231 evolution of influenza have shown that major antigenic changes occur in a punctuated 232 fashion, allowing substantial increases in population susceptibility and in turn surge of 233 influenza epidemic (4,26,50). To capture such punctuated antigenic changes, our model-234 inference system allowed abrupt increases in population susceptibility in response to 235 unexpected epidemic surges that are potentially due to major antigenic evolutions in the 236 influenza virus. We then used the major increases in population susceptibility (i.e., >15% 237 increase relative to the preceding week; results using 10% as cutoff are consistent; see 238 Table S1) to identify potential punctuated antigenic changes. Over the 20-year study 239 period, we detected such major changes in population susceptibility 7 times for A(H1N1), 240 including two due to the 2009 pandemic (Fig 2A), 12 times for A(H3N2) (Fig 2B), and 4 241 times for B ( Fig 2C). Table 2 shows the estimated timing of each potential antigenic 242 innovation. Based on these estimates, punctuated antigenic changes occurred every 4.6 243 (SD=4.6; exponential distribution) years in A(H1N1) (excluding the 2009 pandemic), 1.8 244 (SD=1.2; gamma distribution) years in A(H3N2), and 6.2 (SD=6.2; exponential distribution) 245 years in B during Jan 1998 -July 2018. These results are consistent with reported 246 estimates based on antigenic testing of influenza virus isolates (4, 50). Further, 247 interestingly, for A(H3N2), such potential antigenic innovations occurred with comparable 248 frequencies in both summers (7 times in May, June or July) and winters (5 times in Jan or 249 Feb) during the 20-year study period. In contrast, it occurred mostly in winters for 250 A(H1N1)  seasons. Here we estimate that R0 was similar for A(H1N1) and B (range: 1.1-1.8 for both) 272 and was slightly higher for A(H3N2) (range: 1.3-2.0). Estimated simultaneously with R0, 273 the infectious period, or generation time, was slightly longer for A(H3N2) and B than 274 A(H1N1); these estimates were very close to those based on viral shedding curves-3.0 vs 275 3.1 days reported for A(H3N2), 3.1 vs 3.4 days reported for B, and 2.6 vs 2.3 days reported 276 for A(H1N1) (52). Together, the higher transmissibility of A(H3N2) (i.e., larger R0 and 277 longer infectious period) in part explains its larger seasonal epidemics and resulting lower 278 average population susceptibility (mean: 60.8% for A(H3N2) vs. 66.6% for A(H1N1) and 279 68.6% for B). Further, these consistent estimates indicate that our strategy using the 280 empirical R0 seasonal cycle in the model-inference system is effective in controlling for the 281 elusive influenza seasonality in Hong Kong and allows for estimation of other key 282 epidemiological features. Future studies could apply the same strategy to study the diverse 283 influenza dynamics in other subtropical/tropical regions or other infections with similar 284 challenges. 285 286 Cross-immunity and co-circulation pattern. Multiple lines of evidence have suggested 287 cross protection offered by infection of influenza viruses of the same type/subtype as well 288 as across types/subtypes (7-18). However, some studies have reported weak or no 289 heterosubtypic and/or heterotypic cross-immunity (47, 53). As noted by authors of those 290 studies, the cross-sectional serological surveys used therein may miss the window of cross-291 immunity and/or immunity mechanisms not related to antibodies against the 292 hemagglutinin surface protein (e.g., from T cells mediated immunity). In this study, based 293 on the long-term type/subtype-specific incidence data, we estimate that there were 294 relatively strong cross-immunity against A(H3N2) by infection of A(H1N1) (~40% of the 295 are consistent with estimates from phylogenetic analyses. For instance, the two surges in 342 susceptibility to A(H1N1) during the 2009 pandemic (i.e., around 9/13/09 and 1/5/10; 343 Table 2) roughly matched with the estimated timings of amino acid mutations from 344 NextStrain (55): 1) HA1: S185T and HA2: S124N around 10/1/09 (CI: 7/26/09-1/16/10) 345 and 2) HA1: D97N round 10/30/09 (CI: 09/03/09-02/07/10). In addition, previous 346 studies estimated that punctuated antigenic innovations occur roughly every 2 to 8 years in 347 A(H3N2) (26, 50), 4 to 9 years in A(H1N1), and 3 to 14 years in B (4), based on antigenic 348 analysis of influenza virus isolates. Consistently, we estimate that antigenic innovations 349 occurred very 4.6±4.6 years in A(H1N1), every 1.8±1.2 years in A(H3N2), and every 6.2±6.2 350 years in B during Jan 1998-July 2018. These consistent results suggest incidence 351 surveillance data could be an inexpensive alternative to support detection of evolutionary 352 changes in influenza viruses, especially since these estimates reflect changes in population 353 immunity in response to antigenic innovations and thus the potential of epidemic surge. 354

355
In addition, we note the different seasonal timing of antigenic innovations for 356 A(H3N2). While the identified antigenic innovations predominantly occurred in winters for 357 A(H1N1) and B, we find that A(H3N2) could undergo major antigenic changes in both 358 summers and winters. Further investigation is warranted, as, if confirmed, this seasonality 359 suggests that sampling A(H3N2) strains in the summer may be essential for selection of 360 vaccine strains. 361

362
Our study has several limitations. First, as B lineage data were only available since 363 2014 in Hong Kong, in this study, we did not differentiate the two B lineages (i.e. Yamagata 364 and Victoria). Future work will extend the model-inference system to include and test the 365 epidemiological interactions with the two B lineages. Second, while using the empirical 366 seasonal cycle to account for the epidemic seasonality in Hong Kong has proven critical in 367 this study and may be applied to other subtropical and tropical regions with similar diverse 368 epidemics, better delineation of influenza seasonality in these climates is needed. Third, 369 we assumed that the probability of seeking medical attention for influenza is comparable to 370 other illnesses and, for simplicity, we did not consider the differences in health seeking 371 behavior among patients infected by different influenza (sub)types. Based on this 372 assumption, the estimated annual attack rates ranged from 13% in 2013 to 42% in 2005. 373 There is large uncertainty in estimated influenza attack rates. Recent estimates ranged 374 from 3-11% symptomatic infections among vaccinated and unvaccinated US residents 375 during 2010-2016 (56)to 32% confirmed infection among unvaccinated individuals during 376 2015 in New Zealand (57). In comparison, our estimates for Hong Kong are relatively high 377 but plausible given the longer epidemic duration and lower vaccination rates in Hong Kong 378 (~10% (58, 59) v. ~45% in the US (60)). It could also be due to overestimation of A(H3N2) 379 incidence rates, as patients infected with A(H3N2) may be more likely to seek treatment 380 due to more severe symptoms compared to A(H1N1) and B. Further, while we calibrated 381 the incidence rates of A(H1N1) during the 2009 pandemic to the estimated attack rate from 382 serological data (see Methods), our model-inference system did not include public health 383 interventions implemented particularly at the beginning of the pandemic. As a result, our 384 model-inference system-entirely based on incidence data-likely underestimated the 385 increase in susceptibility at the onset of the pandemic. Nevertheless, the estimated timing 386 of antigenic innovations (6 July 2019) was very close to the observed onset of the pandemic 387 in Hong Kong (30 June 2019 (61)). Finally, due to a lack of age-specific incidence data and 388 for simplicity, we did not include age structure in our model. Recent studies have reported 389 differences in immune imprinting in children compared to adults (47,48 and also a small number of samples submitted by the sentinel outpatient physicians. We 414 obtained weekly records of ILI consultation rate and concurrent (sub)type-specific 415 influenza detection rate from the week ending 04 Jan 1998 to the week ending 14 July 416 2018. To obtain a more specific measure of the incidence of influenza virus infections in the 417 community, we multiplied the weekly ILI rate by the concurrent viral isolation rate for 418 A(H1N1), A(H3N2) and B, respectively (42, 49, 65, 66). We refer to these combined 419 measures as A(H1N1)+, A(H3N2)+, and B+. 420 421

Converting ILI consultation rate to infection rate 422
The ILI consultation rate is the ratio of patients presenting ILI to all patient-visits and does 423 not measure true influenza incidence rate in the catchment population. Previously, we 424 showed that conversion to per-capita infection rate can be made based on Bayes' rule (65, 425 67). Briefly, ILI+ (here, ILI+ = A(H1N1)+, A(H3N2)+, or B+) estimates the probability that a 426 person seeking medical attention, m, has influenza, i.e., p(i|m). By Bayes' rule, the incidence 427 rate, or probability of a person contracting influenza during a given week, p(i), is: 428 This formula can thus convert ILI+ to the incidence rate p(i) if ( = ) * ) $ " is known. 430 Here to reduce model uncertainty and complexity, we estimated γ separately based on 433 independent estimates of attack rate during the pandemic. Based on serological data, 434 estimated attack rates for the 2009 pandemic ranged from 16% (68) to ~20% (61, 69). 435 Multi-strain SIRS model. A number of mathematical models have been developed (19-29, 456 70-73) to study the interplay of influenza viruses. The central idea behind these models is 457 that cross-immunity can either reduce the chance of reinfection and transmission in any 458 individual with prior infection by a similar strain (i.e. leaky immunity) or protect a portion 459 of individuals with prior infection (i.e. polarized immunity). In particular, Gog and Grenfell 460 (24) developed a parsimonious status based multi-strain model by assuming polarized 461 immunity and cross-immunity acted to reduce transmission risk. This model is able to 462 generate evolutionary dynamics similar to that observed for influenza and is commonly 463 used as the backbone for more complex phylodynamic models (26,28,73,74). The multi-464 strain SIRS model in this study took a construct similar to the Gog & Grenfell model: 465 where N is the population size; Si and Ii, are, respectively, the numbers susceptible and 467 infected, to virus-i (here A(H1N1), A(H3N2), or B); βi, Di, and Li are, respectively, the 468 transmission rate, mean infectious period, and mean immunity period, for virus-i; and cij 469 measures the strength of cross-immunity to virus-i conferred by infection of virus-j (e.g., 470 close to 0 if it is weak and cii=1 for strain-specific immunity). We assumed the total 471 population size is stable over the study period such that the birth rate and death rate are 472 equal and both are represented by parameter μ. The parameter α represents travel-related 473 importation of cases and for simplicity was set to 0.1 per 100,000 per day throughout the 474 study period. The model was run continuously in daily time-step from 1 Jan 1998 to 14 July 475 2018 in conjunction with the SR-IF2 filter ( week t of the year was then computed as R0,i(t)/Di(t). For the immunity period Li, given the 510 wide range reported in the literature (from months to ~8 years (76, 77)), we tested prior 511 ranges from 1 to 9 years, divided into four 2-year segments (i.e., [1,3], [3,5], [5,7], and [7, 512 9] years) for each type/subtype; this resulted in 4 3 =64 combinations for the three 513 type/subtypes. Cross-immunity has been observed for influenza strains of the same 514 subtype as well as across types/subtypes (7-15) and types (16-18); however, the strength 515 has rarely been reported. We thus tested prior ranges from 0% to 80% of the strength of 516 specific immunity, divided into two levels: [0%, 40%] (low) and [40%, 80%] (high); this 517 resulted in 2 6 =64 combinations for 6 virus-pairs (cij, i=H1N1, H3N2, or B and j i). To 518 reduce the number of combinations, a priori, we further assumed that heterotypic cross-519 immunity is low especially when the immunity period (Li) is long; thus, for Li>5 years, we 520 only included combinations with low A(H1N1) (or H3N2) and B interactions. Note, 521 however, with perturbations in the IF2 procedure, the posterior estimates can go outside 522 the prior range. Together, we tested a total of 736 combinations of Li and cij priors in the 523 first round of search. We then selected the top ~1% runs with the best model-fits (highest 524 likelihood and fewest space re-probing needed) used their posterior estimates to compute 525 the prior range for subsequent round of search. We repeated this process until the 526 parameters converge. It took two rounds for the synthetic data and three rounds for the 527 observed data. 528

529
For each SR-IF2 model-inference run, we performed 50 iterations of filtering with 530 5000 particles. Each run took round 8 to 10 hours to complete on a high-performance 531 cluster (https://www.mailman.columbia.edu/information-for/information-532 technology/high-performance-computing-hpc). To account for stochasticity in model 533 initiation, we performed 5 runs for each prior in the first round and 50 runs in the second 534 and third round (i.e., final round for the observed data). From the final round, we selected 535 ~10 runs with the best model-fits (highest likelihood and fewest space re-probing needed) 536 and computed the posterior estimates across those runs per Rubin's rule (78,79). 537 538 Estimating the timing and frequency of punctuated antigenic change. To identify 539 potential punctuated antigenic changes for each (sub)type, we computed the changes in 540 susceptibility for each week and considered weeks with a relative increase of 15% due to 541 punctuated antigenic changes. We also tested a cutoff of 10%, for which the identified 542 weeks were the same except for 2 additional weeks for A(H3N2) (1/10/1998 and 543 author/funder, who has granted medRxiv a license to display the preprint in perpetuity. The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/19008987 doi: medRxiv preprint Table Captions:  Table 1. Estimates of key parameters. Posterior estimates from all runs and weeks are pooled to obtain the mean and full range (in parentheses). Table 2. Estimated timing and frequency of potential punctuated antigenic changes. Dates (mm/dd/yy) of potential punctuated antigenic changes are identified based on above 15% increases in susceptibility; the standard deviations are computed based on estimates from all model-inference runs indicating punctuated antigenic changes around the same time. The time intervals between subsequent antigenic changes are estimated based on four different distributions; the best-fit estimates are shown in bold. For A(H1N1), numbers in the parentheses show results when dates during the pandemic are included.    (B). Blue lines and areas show the mean and 95% CI estimates. For comparison, the observed (black dots) and estimated (black lines) incidence rates are superimposed (right y-axis).

Supplementary Material Captions:
Supplemental testing and results. Table S1. Estimated timing of punctuated antigenic changes using a cutoff of 10% increase in susceptibility, instead of 15% in Table 2.  author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10. 1101/19008987 doi: medRxiv preprint to the true values. In (A) and (B), dots show the true incidence rates (in black; right y-axis), susceptibility (in blue; left y-axis), and R0 (in blue; left y-axis); the lines and surrounding areas show the mean and 95% credible interval (CI) estimates. In (A), vertical dashed black lines show the true weeks with punctuated antigenic changes and red lines show modelestimates. In (C)-(E), red stars show true parameter values; box plots show the median, 75%, and 95% CIs of posterior estimates by the model-inference system; orange dashed lines show the prior ranges tested.

Fig. S3.
Model validation using the second synthetic dataset. Posterior estimates for the incidence rate and population susceptibility (A), basic reproductive number R0 (B), infectious period (C), immunity period (D), and strength of cross-immunity (E), compared to the true values. In (A) and (B), dots show the true incidence rates (in black; right y-axis), susceptibility (in blue; left y-axis), and R0 (in blue; left y-axis); the lines and surrounding areas show the mean and 95% credible interval (CI) estimates. In (A), vertical dashed black lines show the true weeks with punctuated antigenic changes and red lines show modelestimates. In (C)-(E), red stars show true parameter values; box plots show the median, 75%, and 95% CIs of posterior estimates by the model-inference system; orange dashed lines show the prior ranges tested.        author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/19008987 doi: medRxiv preprint 38 Tables: Table 1. Estimates of key parameters. Posterior estimates from all runs and weeks are pooled to obtain the mean and full range (in parentheses). All rights reserved. No reuse allowed without permission.

Strain
author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/19008987 doi: medRxiv preprint author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/19008987 doi: medRxiv preprint

Fig. 2. Estimates of weekly infection rate and population susceptibility for the three 10 types/subtypes: A(H1N1) (A), A(H3N2) (B) and B (C). Black dots show observed incidence 11
rates (left y-axis); the lines run through the dots show model estimates; surrounding grey 12 areas show the 95% credible intervals (CI). Blue lines and surrounding areas show the 13 mean and 95% CI estimates of the population susceptibility. Vertical red lines and dates 14 show identified timing of punctuated antigenic changes. 15 16 0 1 2 3 4 5 6 7 0.0 0.5 1.0 1.5 2.0 0 1 2 3 4 5 6 7    author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/19008987 doi: medRxiv preprint This supplemental document includes 1) Validation of the model-inference system using 5 model-generated mock epidemics; 2) Simulations to test the impact of cross-immunity on 6 long-term epidemic pattern; and supplemental table and figures. 7 8

SI Appendix
1. Validation of the model-inference system 9

Method 10
We first tested our multi-strain SIRS model-SR-IF2 inference system using model-11 generated mock epidemics (i.e. synthetic data). Specifically, we tested its ability to 12 estimate three main epidemic features of interest: 1) The population susceptibility for each 13 virus over time, including timing of major susceptibility increase due to epochal antigenic 14 change. For the latter, we tested three levels of magnitude-small, medium, and large 15 change-with 10%, 15%, and 25% increase in susceptibility, respectively. 2) The strength 16 of cross-immunity between each virus-pair. For this, we tested two combinations of 17 heterosubtypic immunity (i.e., low H3N2→H1N1 and high H1N1→H3N2 cross-immunity for 18 the first synthetic dataset; high H3N2→H1N1 and median H1N1→H3N2 cross-immunity for 19 the second). and 3) The duration of immunity. Reported immunity period for influenza 20 varied substantially, ranging from months to ~8 years (1-4). Thus, we tested prior ranges 21 from 1 to 9 years.

23
To generate the mock epidemics for the above testing, we first generated 100,000 24 random combinations of initial conditions and model parameters except R0, using Latin 25 Hypercube sampling. Weekly R0 values were randomly sampled from a normal distribution 26 with the mean set to the empirical estimate for that week (Fig 1 in the main text) and 27 standard deviation set to 5% of the mean. To mimic punctuated antigenic changes, for each 28 (sub)type, we first identified the observed major epidemic peak weeks (peak ILI+ in the 29 upper ~75%, 60%, and 75% percentiles for the three (sub)types, respectively); for each 30 selected peak week, we increased the susceptibility at the 8 th preceding week (i.e., roughly 31 the onset) by 15% for A(H1N1) (except for the pandemic) and A(H3N2), by 10% for B, and 32 by 25% for the A(H1N1) pandemic (see the selected weeks in Figs S2 and S3). We then ran 33 the multi-strain SIRS model (Eqn 1 in the main text) for each setting and selected 2 34 simulations that most closely resembled the observed epidemics in Hong Kong but had 35 different levels of cross-immunity between the two A subtypes: the first simulated time 36 series had high cross-immunity from A(H1N1) against A(H3N2) but low cross-immunity in 37 the reverse direction; and the second had moderate cross-immunity from A(H1N1) against 38 A(H3N2) and high cross-immunity from A(H3N2) against A(H1N1). We then added Poisson 39 All rights reserved. No reuse allowed without permission.
author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/19008987 doi: medRxiv preprint random noise to the two simulated time series to mimic observational errors and used 40 these as mock data (termed synthetic 'truth') to test the model-inference system.

42
For each of the two synthetic truths, we performed two rounds of model-inference 43 using the same procedure and prior distributions for the Hong Kong dataset as described in susceptibility for the three types/subtypes) and parameters (the basic reproductive 50 number R0, infectious period, immunity period and strength of cross-immunity) for the two 51 tests, respectively, compared with the true values. The relative root-mean-square-error 52 was 0.16 and the correlation was 0.83 combining all model variables/parameters from 53 both tests. These results suggest that our model-inference system and parameter 54 estimation strategy is able to accurately estimate the model state variables and parameters. 55 However, due to large fluctuations of R0 in the synthetic datasets (e.g., Fig S2B) and the 56 collinearity between R0 and susceptibility, the filter tended to compensate sudden 57 increases in infections (due to increases in susceptibility to mimic punctuated antigenic 58 changes) with increases in R0 rather than susceptibility. As a result, it failed to identify 59 some of the punctuated antigenic changes, in particular, for the two A subtypes with larger 60 variances of R0. 61 62 2. Simulating the impact of cross-immunity on long-term epidemic pattern 63

Method 64
To test the impact of cross-immunity, we simulated epidemics under three different 65 scenarios: 1) As estimated in this study, by setting all cross-immunity terms to the 66 posterior mean estimates; 2) Stronger cross-immunity from A(H3N2) against A(H1N1), by 67 setting cH1←H3 to 0.5, i.e., slightly above the mean estimate for cH3←H1 (i.e. 0.4); and 3) No 68 interactions, by setting all cross-immunity terms to 0. For each simulation, we initiated the 69 multi-strain SIRS model (Eqn 1 in the main text) using the posterior estimates at Week 4 70 (i.e., after discarding the first three weeks of filter spin-off) and ran the model 71 stochastically from Jan 1998 to July 2018. Weekly values of model parameters (i.e., the 72 infectious period, immunity period, R0, and cross-immunity terms) not specified in the 73 above scenarios were set to the mean posterior estimates from the Hong Kong dataset. 74 Weekly susceptibilities to the three types/subtypes were simulated by the model according 75 to model-simulated infections, loss of immunity, cross-immunity, replenishment from 76 newborns, and deaths; however, to simulate the impact from antigenic innovations, we 77 reset the susceptibility to the posterior mean estimate from the Hong Kong dataset for 78 weeks identified to experience punctuated antigenic changes ( Table 1). As in the model-79 All rights reserved. No reuse allowed without permission. author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/19008987 doi: medRxiv preprint inference runs, we set the epidemic seeding (α in Eqn 1) to 0.1 per 100,000 per day for the 80 all simulations. To account for model stochasticity, we simulated each scenario for 1000 81 times.

83
To compute the number of epidemics over the study period, as in (5), we defined the 84 epidemic baseline, for each type/subtype, as the 40% quantile of the non-zero 85 corresponding ILI+ records over the ~20 year study period. We then identified 1) all 86 epidemics including small ones, defined as periods with ≥3 consecutive weeks with ILI+ 87 above the baseline and at least one of the weeks with ILI+ ≥3 times of the baseline; and 2) 88 all large epidemics, defined as periods with ≥6 consecutive weeks with ILI+ above the 89 baseline and at least one of the weeks with ILI+ ≥6 times of the baseline. If an identified 90 epidemic lasted for >1 year, we divided it into multiple epidemics at the year division(s). To compare the co-circulation patterns, we categorized the weeks into 8 possible 98 co-circulation types: 1) none in circulation, i.e., none of the three types/subtypes had ILI+ 99 above their baselines; 2) A(H1N1) alone, i.e., only A(H1N1) had ILI+ above its baseline; and deviated from the observations. However, overall, simulations with cross-immunity in 112 place were able to better reproduce the observed epidemic pattern (Figs. S8 and S9 vs. Fig.  113 S10). In particular, those simulations were able to reproduce the long periods with near-114 zero A(H1N1) activities (e.g., [2002][2003][2004][2005] whereas simulations with no cross-immunity 115 (Fig. S10) overestimated the circulation of all types/subtypes, especially for A(H1N1). 116 Comparing the first two scenarios (i.e. both with cross-immunity), simulations with 117 stronger cross-immunity from A(H3N2) against A(H1N1) underestimated the circulation of 118 A(H1N1) and were only able to generate A(H1N1) epidemics following the simulated 119 All rights reserved. No reuse allowed without permission. author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/19008987 doi: medRxiv preprint punctuated antigenic changes when the susceptibility was reset (Fig. S9 vs. Fig. S8 as estimated from the data).

122
Tally over the entire study period, the first two scenarios (i.e., with cross-immunity) 123 appeared to underestimate the total number of the epidemics (Fig. S11 A and B). However, 124 this was due to the challenge in generating clear bimodal epidemics within a year as those 125 observed in Hong Kong. For instance, there were two consecutive A(H1N1) epidemics in 126 2008; and similarly, A(H3N2) tended to cause two separate epidemics within a year, e.g. , 127 1999, 2000, 2003, and 2014. The underestimation of number of A(H1N1) epidemics was 128 severer for scenario 2 when stronger cross-immunity from A(H3N2) against A(H1N1) was 129 implemented in the model. Simulations with no cross-immunity were able to generate 130 more epidemics due to overestimation of circulation during periods without observed 131 influenza activity (Fig. S10).

133
Comparing the 8 co-circulation types (Fig. 11C), scenario 1 using the posterior mean 134 estimates most closely reproduced the observed pattern. In contrast, scenario 2 (i.e., 135 stronger cross-immunity from A(H3N2) against A(H1N1)) largely underestimated the 136 frequency of A(H1N1) epidemics, whereas scenario 3 (i.e., no cross-immunity) largely 137 overestimated co-epidemics caused by multiple types/subtypes.  author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/19008987 doi: medRxiv preprint Supplementary Table  Table S1. Estimated timing of punctuated antigenic changes using a cutoff of 10% increase in susceptibility, instead of 15% in Table 2 author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/19008987 doi: medRxiv preprint Fig. S1. Example model fits and estimates of key parameters using a particle filter. Model fits compared to the observations are shown for A(H1N1) (A), A(H3N2) (B), and B (C). Estimates of R0 are shown in D-F for the three influenza types/subtypes and strength of cross-immunity shown in G-L.

Supplementary figures
Infection rate (%) (B) Infection rate: A(H3N2) author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/19008987 doi: medRxiv preprint Fig. S2. Model validation using the first synthetic dataset. Posterior estimates for the incidence rate and population susceptibility (A), basic reproductive number R0 (B), infectious period (C), immunity period (D), and strength of cross-immunity (E), compared to the true values. In (A) and (B), dots show the true incidence rates (in black; right y-axis), susceptibility (in blue; left y-axis), and R0 (in blue; left y-axis); the lines and surrounding areas show the mean and 95% credible interval (CI) estimates. In (A), vertical dashed black lines show the true weeks with epochal antigenic changes and red lines show modelestimates. In (C)-(E), red stars show true parameter values; box plots show the median, 75%, and 95% CIs of posterior estimates by the model-inference system; orange dashed lines show the prior ranges tested.
Strain 2 0 20 40 60 80 100  Immunity period (day) 1 y 2 y 3 y 4 y 5 y 6 y 7 y 8 y 9 y * * * (D) Immunity period Model validation using the second synthetic dataset. Posterior estimates for the incidence rate and population susceptibility (A), basic reproductive number R0 (B), infectious period (C), immunity period (D), and strength of cross-immunity (E), compared to the true values. In (A) and (B), dots show the true incidence rates (in black; right y-axis), susceptibility (in blue; left y-axis), and R0 (in blue; left y-axis); the lines and surrounding areas show the mean and 95% credible interval (CI) estimates. In (A), vertical dashed black lines show the true weeks with epochal antigenic changes and red lines show modelestimates. In (C)-(E), red stars show true parameter values; box plots show the median, 75%, and 95% CIs of posterior estimates by the model-inference system; orange dashed lines show the prior ranges tested. Immunity period (day) 1 y 2 y 3 y 4 y 5 y 6 y 7 y 8 y 9 y * * * (D) Immunity period author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/19008987 doi: medRxiv preprint  author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/19008987 doi: medRxiv preprint author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/19008987 doi: medRxiv preprint All rights reserved. No reuse allowed without permission. author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/19008987 doi: medRxiv preprint Fig. S9. Simulations assuming stronger cross-immunity from A(H3N2) against A(H1N1). The cross-immunity parameter cH1←H3 was set to 0.5 and other parameters set to the posterior mean estimates. Colored lines show simulated weekly infection rates from 1000 individual stochastic model runs (A(H1N1) in orange, A(H3N2) in red, and B in green); dashed black lines show the weekly mean infection rates across 1000 simulations and solid black lines show the weekly observations for comparison. All rights reserved. No reuse allowed without permission. author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/19008987 doi: medRxiv preprint Fig. S10. Simulations assuming no cross-immunity. All cross-immunity terms were set to 0 and other parameters set to the posterior mean estimates. Colored lines show simulated weekly infection rates from 1000 individual stochastic model runs (A(H1N1) in orange, A(H3N2) in red, and B in green); dashed black lines show the weekly mean infection rates across 1000 simulations and solid black lines show the weekly observations for comparison.
All rights reserved. No reuse allowed without permission. author/funder, who has granted medRxiv a license to display the preprint in perpetuity.
The copyright holder for this preprint (which was not peer-reviewed) is the . https://doi.org/10.1101/19008987 doi: medRxiv preprint