Global Spatio-temporal Patterns of Influenza in the Post-pandemic Era

We study the global spatio-temporal patterns of influenza dynamics. This is achieved by analysing and modelling weekly laboratory confirmed cases of influenza A and B from 138 countries between January 2006 and January 2015. The data were obtained from FluNet, the surveillance network compiled by the the World Health Organization. We report a pattern of skip-and-resurgence behavior between the years 2011 and 2013 for influenza H1N1pdm, the strain responsible for the 2009 pandemic, in Europe and Eastern Asia. In particular, the expected H1N1pdm epidemic outbreak in 2011/12 failed to occur (or “skipped”) in many countries across the globe, although an outbreak occurred in the following year. We also report a pattern of well-synchronized wave of H1N1pdm in early 2011 in the Northern Hemisphere countries, and a pattern of replacement of strain H1N1pre by H1N1pdm between the 2009 and 2012 influenza seasons. Using both a statistical and a mechanistic mathematical model, and through fitting the data of 108 countries, we discuss the mechanisms that are likely to generate these events taking into account the role of multi-strain dynamics. A basic understanding of these patterns has important public health implications and scientific significance.

of 2009 and early 2011. Unusually an H1N1pdm outbreak failed to appear at all in the "skip year" of 2011/12, given that the strain was very new, although the outbreak returned and resurged in 2012/13. As we will discuss shortly, this same pattern was generic to many countries across Europe, with slight differences from country to country. A visualisation of the extraordinary skip-year across 45 countries is given in Fig. 2.
Many of the features of the time series in Fig. 1 can be explained in terms of basic epidemiological theory. Briefly, when the new 2009 pandemic influenza strain confronted a large susceptible human population it was able to generate a large-scale global epidemic. This placed in motion a succession of epidemic "waves" that followed one after the other. Since infected individuals who recover from the disease gain temporary immunity, each new epidemic wave also served to build up further the level of immunity in the population. In effect, this served to reduce the number of susceptible individuals available for future infection. At some point, when the number of available susceptibles fell below a threshold level, it became impossible for a new outbreak of the pandemic strain to trigger. This explains the "skip" year in 2011/12 in which the strain was mostly absent (see Fig. 2). The H1N1pdm strain resurged in 2012/13 presumably because recovered individuals gradually lost their immunity, providing enough new susceptibles to trigger further outbreaks. A systematic theory for understanding epidemic oscillations and skips has been developed over the last decade 8 , which we will use to explain these long term dynamics.
Our detailed spatio-temporal analysis is based on time series data obtained from FluNet, a comprehensive global surveillance tool for influenza developed by the World Health Organization (WHO) in 1997 9 , in which virological data are documented in real-time and publicly available. When discussing and presenting the FluNet data it is convenient to use the following notation. We denote H1N1pre as the pre-2009-pandemic seasonal A strains, H1N1pdm as the pandemic strain (H1N1pdm09) responsible for the 2009 influenza pandemic, and H3N2 as the seasonal H3N2 strains whose original form was responsible for the 1968 influenza pandemic.
We note that FluNet has previously been employed to study the spread of influenza on global or large-scale spatio-temporal patterns in three other studies that we know of. Finkelman et al. 10 studied the pre-2009-pandemic period between January 1997 and July 2006 in 19 temperate countries in both Hemispheres. They identified large scale co-existence of influenza A and B, interhemispheric synchronized pattern for subtype A H3N2, and latitudinal gradients in the epidemic timing for seasonal influenza A. A recent study 11 that focused on the Western Pacific Region between January 2006 and December 2010, found that dominant strains of influenza A were reported earlier in Southern Asia than in other countries. Thus, status in South Asian countries may provide early warning for other countries. Bloom-Feshbach et al. 12 examined latitudinal variations in seasonal activity of influenza and respiratory syncytial virus (RSV) and applied a time series model to the seasonal influenza data from 85 countries. They found evidence of latitudinal gradients in timing, duration, seasonal amplitude and between-year Countries with α 1 < log1/10 are coded in red and countries with log1/10 < α 1 < log1/5 are in blue. The color (grey) scheme is in a square-root scale. variability of epidemics. In terms of the temporal pattern in a single region, Dorrigatti et al. 13 studied the third wave of infection by the H1N1pdm pandemic strain in England in the 2010/11 season. They found that increased transmissibility and loss-of-immunity among the population may be responsible for this unexpected wave.
However, to our knowledge, no study has been conducted focusing on the global pattern of seasonal activities of the H1N1pdm pandemic virus and interactions among different strains based on the FluNet large-scale dataset from 2010 to 2013. This is of special interest given that the surveillance scale was substantially improved since 2010. Such a study is important to aid the development of strategies for relieving the burden of seasonal influenza. Understanding the spatial pattern may be useful in a global effort to reduce the impact of a deadly influenza pandemic. The activity of H1N1pdm still causes substantial attention in the post-pandemic era and has led to substantial morbidity and mortality in most years since its appearance, including 2013/14.
On the global network of influenza transmission, it is known that China and Southeast Asia lie at the center of the global network and USA acts as the primary hub of temperate transmission 14,15 . The expansion of H1N1pdm during 2009 can be explained with data on human mobility (air travel) and viral evolution 16 .

Materials and Methods
Weekly time series data of lab-confirmed cases (isolates) of influenza were obtained from FluNet for 138 countries that have non-zero cases between January 2006 and January 2015. The analysis included six different types of time series: i) total specimens processed, ii) H1N1pre strains, iii) H3N2 strains, iv) H1N1pdm strain, v) un-subtyped influenza A, and vi) influenza B (including two circulation lineages).
The number of un-subtyped influenza A is often substantial and needs to be accounted for. Following 10 , we proportionally assigned the un-subtyped influenza A to the three subtypes as follows. Let the number of lab-confirmed cases of subtypes H1N1pre, H1N1pdm, H3N2 and un-subtyped A for any country in a particular week be a 1 , a 2 , a 3 and a 0 , respectively. Then the new revised number for each of the three subtypes was taken to be: The statistical analysis was implemented in the R programming language (http://www.r-project.org/). We generally preferred to focus on regions (macro geographical continental regions and geographical sub-regions) rather than their constituent countries (see Fig. 1) because aggregated regional data is less influenced by stochasticity. We observed that nearly all countries, especially those in the temperate regions, largely followed their regional patterns.
The breakdown of countries of the eight regions used in this study may be found in the Supplementary Material §S9.
To compare between hemispheres it was convenient to redefine the initiation and termination dates of calendar years in a manner that makes influenza seasons (usually winter) line up. We therefore moved forwards the beginning and end dates that define years for Northern Hemisphere (NH) countries to stretch from the 35th week of a calendar year to the 34th week of the following calendar year, roughly overlapping the school calendar year. For countries in the Southern Hemisphere (SH), the calendar year remains the reference frame. Thus the skip-year (skip-season) of H1N1pdm is 2011/12 in NH and is 2012 in SH, see Fig. 2.
A skip-year, or simply a skip, is defined as a season with an uninitiated or minor epidemic. What constitutes a "minor" epidemic is difficult to quantify precisely. For the purposes of this study, we formulated the following practical quantitative comparative definition. If, after an epidemic year the number of influenza cases drops by more than a factor of ten, we consider this to be a skip year. We thus use the following simple measure defined for H1N1pdm: where h 10 , h 11 are the total number of lab-confirmed cases of H1N1pdm in that region during the 2010/11 and 2011/12 seasons, respectively. The index compares the ratio of the number of cases in 2011/12 season to those in the 2010/11 season in NH (or 2012 to 2011 in SH). Our criterion for a "skip" is generally that α 1 < log(1/10), i.e., an order of magnitude difference. We set k = 50 in α 1 to reduce errors magnified when case numbers are small. The merit of using skip index (which is a ratio of two years) rather than the actual number is clear. In this way, we can remove the differences in the testing effort among countries and we can also remove the effects of different population sizes among countries. We assume that the total numbers of specimens processed had not varied much from year to year, which was true from 2010 to 2012.
Similarly we define skip-indices α 2 , α 3 , α 4 , for subtype H3N2, influenza B and total specimens processed, respectively. We argue that as long as surveillance efforts and testing policies were implemented consistently in each country between 2010 and 2013, then effects due to differences in testing policies are removed by taking the ratio of confirmed cases over total cases in consecutive years. To our knowledge, there were no dramatic changes in surveillance effort in most countries from 2010 to 2013 (as observed from the total number of specimens processed).

Results
Skip-and-resurgence pattern of H1N1pdm. Figure 1 panels (a-h) show weekly aggregated lab-confirmed cases of subtypes H1N1pdm (red triangle) and H3N2 (blue circle) in the eight geographical regions having the largest case numbers in the period January 2009 and January 2015. A similar set of panels for thirty different countries (having the largest number of confirmed cases) may be found in Fig. S2 in the Supplementary Material . It should be emphasised that the graphs are scaled to highlight the trends (and also accommodate the extremely high peak in 2009) by plotting the square-root of the weekly lab-confirmed cases. The figures immediately identify a number of clear features. With regard to H1N1pdm, we summarise here: (i) All regions in Europe and Eastern Asia have identical trends and experienced skip-years in 2011/12.
In more detail these regions experienced two initial waves of H1N1pdm in 2009/10, followed by a single wave in 2010/11, a skip-year in the 2011/12 season and then a reemergence of H1N1pdm in the following 2012/13 season. The skip was more evident in Eastern/Southern Europe than Western/ Northern Europe. Although the latter experienced a minor outbreak in 2011/12, its peak was an order of magnitude lower than the previous season, and thus by our criteria could be classified as a skip. The size of the mini-outbreak is misrepresented and appears exaggerated due to the square-root scaling.  (Fig. 1b). The dynamics over these years were essentially biennial. (iv) South America shows an irregular pattern (Fig. 1c).
With regard to H3N2 dynamics, we observe: (i) All regions in Europe and Eastern Asia experienced significant H3N2 epidemics in 2011/12, which was a skip-year for H1N1pdm. Moreover, apart from 2012-2014, the H3N2 dynamics were essentially negatively correlated with H1N1pdm. (ii) In Northern America H3N2 tended to oscillate synchronously in-phase with H1N1pdm in 2010-2012.
A spatial summary of the dynamics of each geographic region has been superimposed on the world map of Fig. 1 panel (i). The regions colour coded in cyan experienced a skip year in 2011/12 and constitute a considerable proportion of the global map.
We identified 27 countries with α 1 < log(1/10) and thus skip years. Total confirmations of the skip year was one order of magnitude lower than that of the previous year. Using a higher threshold α 1 < log(1/5), the number increases to 45 countries. Namely 18 countries have α 1 between log(1/10) and log(1/5). The weekly confirmations of these latter countries are displayed in Fig. 2 (in latitude order) which gives a remarkable demonstration of the broad geographic synchrony of the epidemic skip over the globe. We note that most of the 45 countries experienced a resurgence of H1N1pdm in 2012/13. To study this relation in more detail, we tested whether these two strains are correlated across all countries in the 2011/12 season alone. That is, we asked whether countries with smaller H1N1pdm outbreaks tend to have larger H3N2 outbreaks, in 2011/12 in NH (or 2012 in SH). Using the skip-index, eqn. (1), our analysis showed that of the 108 countries which reported more than 500 cases of all strains, the correlation coefficient across all countries is r = − 0.61. (Without the threshold of 500, the correlation is r = − 0.63.)

2011/12 skip year and strain dynamics.
For a more detailed analysis of the 2010-2012 years, we considered a linear model with α 1 as the response, and having seven different predictors: α 2 (H3N2), α 3 (flu B), rank of population size in the year 2005, rank of area, rank of absolute latitude, rank of distance from Mexico and geographical region code. Regarding the distance from Mexico, we considered both Euclidean distance (defined as + x y 2 2 where x, y are in terms of longitude and latitude, respectively) and effective distance 17  where the c i are constants to be fitted. The H3N2 skip index (α 2 ) was found to be a significant predictor (p-value < 0.001) of α 1 though negatively correlated, while α 3 (influenza B) was not a significant predictor (p-value ≈ 0.452). Region code and area rank were also found to be significant predictors (p-value < 0.001 and ≈ 0.01 respectively), while all other predictors were not significant. These results parallel our observations that in broad terms, countries in the same region share a common pattern.
In this study, we have focused largely on the 2011/12 skip. However, it is evident from Fig. 1 of 10  The following assumptions were made when fitting the model: • The initial susceptible proportion of the population was taken to lie between 40% and 75% for all countries rather than 100%. This takes into account that many of the elderly population had preexisting cross-reactive antibodies 19 . In addition it was found that the cross-protection provided by the pre-pandemic vaccine was as high as 19% 20 . The model is used to estimate the actual number of initial susceptible in the population • The transmission rate β(t) was taken to be seasonal and modelled by a periodic function of time, with a period of one year. Weather variations and school terms are understood to be responsible for the seasonal variability 21,22 . We adopted a seven-node cubic spline function, and fixed the parameter of node seven to be equal to node one. The function is second-order differentiable except for the seventh node. Thus there were six free parameters in the transmission rate β(t). • The reporting rate ρ(t) of each country was modelled by a three-piece step function of the following form: • If the infection dies out in a country after the invasion in the simulation, we introduced a single infected individual. This mimics the transmission of influenza between countries, so that no country is completely isolated. • We fix σ = 365 and γ = 182.5. Thus the latent and infectious periods are approximately 1.58 and 2.54 days which lie in biologically reasonable range 25 . We used eqn. A7 in 18 to calculate the latent and infectious periods in a time discretized setting with a time step size of 1 day. Also choosing a smaller time step size (such as 0.5 day) has a negligible effect on the results 23 . • We also fit the duration of the immunity (λ −1 ) by calculating the maximum log likelihood profile for it. Namely we fixed the duration to eight discrete values spanning from 1.5 to 7 years, and maximized the performance of the model while fitting other parameters. Then from this profile we estimated λ −1 , and its 95% confidence interval 18 .
The model essentially finds the best fitting estimates of the transmission rate and reporting ratio (as defined above) to the influenza A time series data by maximizing the relevant log likelihood. The output of the model is a plot of the profile log likelihood as a function of the duration of immunity. For example, Fig. 3 shows the best fitting model to the FluNet time series data (plotted in black) from ten countries when aggregating influenza A (i.e., by combining H1N1pdm data with H3N2. The inset figure plots the likelihood profile and shows that the maximum occurs when the immunity duration (λ −1 ) is approximately two years in most countries and four years in some countries. In a recent study on H1N1pdm, Dorigatti et al. 13 found that "the half-life of the decay of prior immunity is estimated to be ~1 year" where the authors only considered the first three waves in England. Cowling et al. 26 found that "homosubtypic immunity against (influenza) infection lasted for at least 18 months" in a three-year clinical study. In summary, our estimates of ~2 years (and some around 4 years) is not far away from other recent studies.
We plot time series data of the median value of reported cases for 1000 model simulations. The median values are plotted in red, while the grey shaded regions indicate the 95% confidence interval. The median values sit close to the observed values (black lines) for all ten countries. The observed data falls well within the envelope of simulated model runs (and their 95% CIs), and therefore this signals that the stochastic model is performing as it should; we cannot reject the model fit as statistically implausible.
The fitted reporting ratio is shown in the top-left corner. The simulations match most of the observed waves for all ten countries. The estimated seasonal amplitudes in the transmission rate are small and largely consistent across countries, which suggests that the post-pandemic waves are largely due to loss-of-immunity and its associated replenishment of the susceptible pool (possibly impacted by dynamical resonance 27 ). The estimated reporting ratio is small and varied considerably across countries. (Note that in 13 , it was found that only 0.7% infected individuals actually visited a General Practitioner. ) We also attempted to use the same model for fitting single strain data including H1N1pdm alone and H3N2 alone. However, when modelling a single strain (either H1 or H3), fits were generally poor. In particular, the model was unable to predict the skip in 2011/12 for any country. Since the main shortcoming of the model used is that it is only a single strain, we conclude that a multi-strain model that includes the interaction between the H1N1pdm and H3N2 strains is needed to capture the 2011/12 skip. The model will be substantially more complicated and is beyond the scope of this paper.
We also considered using a sinusoidal function (three parameters) to replace the seven-node (six free parameter) cubic spline function. We found that the cubic spline model performs better in 9 out of ten countries. The difference in the AIC c 24 , the estimated reproductive number  0 and initial susceptible proportions are given in Table S5 (Supplementary Material). Synchrony Patterns. We examined the synchrony pattern across Northern Hemisphere (NH) countries (latitude > 29°) for the three strains, H1N1 (combined H1N1pdm and H1N1pre), H3N2, and influenza B. We focused on the period from Jan 2006 to Jan 2015.
In order to quantify synchrony, following 10  If all countries have the same MCT, so that MCT j = c is the same constant for each country-j, then the distribution of the MCT j is just a single spike indicating that the countries are highly synchronized. Obviously, the smaller the standard deviation amongst the different MCT j the more synchronized are the different countries. The synchrony analysis examines the distribution of MCT j over all countries j = 1, 2,…, N (median and standard deviation) for different strains and different seasons. Countries with no cases were removed from the analysis. The results are shown in Fig. 4 and listed in Table S1 in the Supplementary Material . We examined the NH countries and for this purpose grouped H1N1pre and H1N1pdm together as H1N1. Fig. 4 show that countries were more synchronized by H1N1pdm than by H3N2 or influenza B in 2010/11 and 2012/13. Indeed, from Table S1 in the Supplementary Material, the standard deviation of MCT for the H1N1pdm strain in the 2010/11 and 2012/13 seasons were significantly smaller than the other seasons and any of the other strains. Table S2 showed the correlation between the epidemic size (total confirmations) versus the median and SD of MCT. From Table S2, we found negative correlations, which suggests that more intense epidemics tend to initiate earlier (small median of MCT) with stronger synchronized pattern (small variance of MCT). This effect is more evident for H1N1pdm and flu B rather than H3N2. This may reflect a more efficient transmissibility of the H1N1pdm virus which allows it to spread more rapidly between countries. These findings corroborate what is observed by eye in Fig. 2, and that the MCT of H1N1 only varies by some 5 weeks across all countries. Note that the median of MCT of influenza B seems larger than the two influenza A strains (by two weeks), suggesting that the flu B epidemic lagged behind the other two flu A strains 10 . We reproduced the same tables for the other countries (latitude < 29) in Table S3 and S4 and did not notice any clear pattern.
Patterns of Strain Replacement. Table 1 lists the annual total confirmations of H1N1pre. It is interesting to note that the 2008 total (pre-pandemic) was double that in 2006, which was due to an increase in testing effort (total specimens processed). The high number in 2009 was most likely due to the extensive testing during the pandemic. The numbers decreased quickly after 2009 suggested that H1N1pre was replaced by H1N1pdm 28 The low numbers in 2012-2014 are likely to be errors, either misclassification or mis-input. For example, six cases of H1N1pre were reported in Poland in the 7th week of 2014. However, close observation revealed that there were minor epidemics of both H1N1pdm and H3N2 in that period, yet data was unexpectedly missing in these categories. But despite its absence anywhere else, six cases of H1N1pre were recorded suggesting possible misdiagnosis. No evidence for an epidemic of H1N1pre after 2011 has been found so far. It is interesting to note that the original form of the H1N1pre subtype had an unusual re-emergence in 1977 some 20 years after its disappearance 29 .

Discussion
Fundamental epidemiological principles are able to explain the skip dynamics seen in Figs. 1,2 in relatively simple terms. When the new strain H1N1pdm first appeared in March 2009, the population at large had no previous exposure to the strain. This allowed the pandemic to develop into a global-scale epidemic even though outside the normal influenza season in many countries. With the passage of time, each successive epidemic outbreak exposed the population at large further to the new H1N1pdm strain, thereby building up population immunity and reducing the number of susceptible individuals 13 . Thus by the end of 2011/12, the susceptible pool of individuals available for infection had reduced below a critical threshold level, so that the epidemic failed to trigger over the 2011 winter season. In epidemiological parlance, by "burning out" the available susceptible pool, the virus effectively reduced the effective reproductive number  below unity making it impossible for the epidemic to initiate in the 2011/12 season. This set the stage for the opportunist H3N2 virus to outcompete and replace H1N1pdm, thus accounting for the H3N2 outbreak at the end of 2011. Fig. 1 makes clear the complex interaction between the H1N1pdm and H3N2 strains as they compete for the available pool of susceptible individuals as well as offer cross-protection.
It is interesting that in Central America where H1N1pdm first appeared, the outbreak progression was different to the above pattern. Two skip-seasons were observed in 2010/11 and 2012/13 (see Mexico in Fig. S2 in the Supplementary Material). These skips were in all likelihood an outcome of the same underlying process, namely a burn-out of susceptibles from the previous waves.
Other mechanisms such as climatic variation, poor surveillance (changes in reporting ratio) and results of new births unlikely played a key role here. There were no previous studies showing that these factors favor H3N2 rather than H1N1pdm. These factors are largely the same between Europe and Northern America. Sampling bias in different age groups might have played some role, as elderly is more vulnerable to H3N2 than H1N1pdm. Unfortunately FluNet data contains no age information.
The occurrence of skips gives information regarding the loss-of-immunity (strain specific) in the population, particularly if there might be only a single viral strain, or if the viral strain is stable and evolves only at a relatively slow rate. The latter is the case for the H1N1pdm strain which is believed to have been antigenically stable since its emergence in 2009 30 . As an indication of its stability, the vaccine component against H1N1pdm recommended by the WHO and the United States Food and Drug Administration (FDA) was not updated since fall of 2009, while those vaccine components against H3N2 and influenza B have been updated more than twice already (see Table S10 in the Supplementary Material). If the H1N1pdm strain was indeed stable over these last years, and the virus evolved relatively slowly, then the main source of new susceptibles in the population was largely derived through natural loss of immunity. In these circumstances, the resurgence of H1N1pdm in 2012/13 after the skip in 2011/12 should be viewed as a consequence of the natural loss-of-immunity in the population 26 .
The differences in influenza dynamics between Northern America (no skip) and Europe (skip) given that they share many common factors with regard to economics, culture, climate and latitude are in some respects surprising. We speculate the different dynamics may be connected with the influenza vaccination coverage which was consistently higher in Northern America than in Europe (and the rest of the world). High coverage of vaccination (against H3N2) among general population could have slowed down the transmission of H3N2, thus saved H1N1pdm from skip a year in Northern America. Vaccination  H1N1pre  13268  18983  29807  37879  709  41  3  5  17   Total Specimens  355834  513430  671232  2290733  1186197  1270287  1350542 1672204 810570  Fig. S3 in the Supplementary Material). Also intriguing is that many parts of Northern America and Central America had much higher attack rates of H1N1 in 2009 and influenza-associated mortality in 2009 was almost 20-fold higher in some countries in America than in Europe (see 31 ). Additional work is still needed to understand which factors are responsible for the different spatio-temporal patterns of influenza seen in Europe and America. It should be noted that the vaccine is trivalent, and thus affords protection against H1N1pdm as well as H3N2. However, in 2011 it had most impact on the H3N2 susceptible members of the populations, since a significant proportion of the population (who received the vaccine such as school children and elderly) were previously infected, and thus naturally immunized against H1N1pdm. If vaccination has most impact on H3N2 susceptibles, then it may release H1N1pdm from competition with H3N2. Such a competitive release could favour the spreading of H1N1pdm, as was seen in the US and Canada. Our attempts to fit the time series data of aggregated influenza A confirmed cases using the same simple SEIR model are shown in Fig. 3 and are reasonably well given that they capture the different trends observed in ten different countries. That the same SEIR model reproduced the different trends suggests that the dynamics of influenza epidemics have a large degree of determinism and the model is considerably robust. Moreover this indicates that the key assumptions behind the SEIR model are largely being met. Namely, the classical mathematical concept of infection being spread by a randomly mixing population and mean field dynamics appear to apply when modelling large real human populations. The different features of each country's influenza A time series can be explained through a change in the SEIR model's parameters. It is also interesting that the model fits (likelihood profiles) indicate a reasonably fast loss of immunity in the vicinity of 2-4 years. This could explain the fast susceptible buildup required during skip years, which would be necessary to generate the resurgent epidemics observed in the following years.
Our analysis has given interesting insights into the global patterns of the invasion of H1N1pdm, which first appeared as a pandemic and then, within a few years, apparently outcompeted and completely replaced the H1N1pre seasonal flu strain. The synchrony between countries of the H1N1pdm outbreaks is striking particularly as witnessed in the highly visible skip (Fig. 2), where for a large number of countries, H1N1pdm failed to outbreak in the 2011 flu season. Moreover, the synchrony between countries in H1N1pdm outbreak years was also very strong (see Fig. 3). The FluNet data gave a comprehensive picture of these phenomenon as they evolved in time over several years, and also in space over 138 countries.