Seasonal transmission potential and activity peaks of the new influenza A(H1N1): a Monte Carlo likelihood analysis based on human mobility

Estimating the transmission potential of a newly emerging virus is crucial when planning for adequate public health interventions to mitigate its spread and impact, and to forecast the expected epidemic scenarios through sophisticate computational approaches [1–4]. With the current outbreak of the new influenza A(H1N1) strain having reached pandemic proportions, the investigation of the influenza situation worldwide might provide the key to the understanding of the transmissibility observed in different regions and to the characterization of possible seasonal behavior. During the early phase of an outbreak, this task is hampered by inaccuracies and incompleteness of available information. Reporting is constrained by the difficulties in confirming large numbers of cases through specific tests and serological analysis. The cocirculation of multiple strains, the presence of asymptomatic cases that go undetected, the impossibility to monitor mild cases that do not seek health care and the possible delays in diagnosis and reporting, all worsen the situation. Early modeling approaches and statistical analysis show that the number of confirmed cases by the Mexican authorities during the early phase was underestimated by a factor ranging from one order of magnitude [5] to almost three [6]. The Centers for Disease Control (CDC) in the US estimate a 5% to 10% case detection, similar to other countries facing large outbreaks, with expected heterogeneities due to different surveillance systems. Even within the same country, the setup of enhanced monitoring led to improved notification with respect to the earlier phase of the pandemic, later relaxed as reporting requirements changed [7].

By contrast, the effort put in place by the World Health Organization (WHO) and health protection agencies worldwide is providing an unprecedented amount of data and, at last, the possibility of following in real time the pandemic chronology on the global scale. In particular, the border controls and the enhanced surveillance aimed at detecting the first cases reaching uninfected countries appear to provide more reliable and timely information with respect to the raw count of cases as local transmission occurs, and this data has already been used for early assessment of the number of cases in Mexico [5]. Moreover, data on international passenger flows from Mexico was found to display a strong correlation with confirmed H1N1 importations from Mexico [8]. Here we present an estimate of the reproduction number, R ₀, (that is, the average number of secondary cases produced by a primary case [9]) of the current H1N1 epidemic based on knowledge of human mobility patterns. We use the GLEaM (for GLobal Epidemic and Mobility) structured metapopulation model [10] for the worldwide evolution of the pandemic and perform a maximum likelihood analysis of the parameters against the actual chronology of newly infected countries. The method is computationally intensive as it involves a Monte Carlo generation of the distribution of arrival time of the infection in each country based on the analysis of 10⁶ worldwide simulations of the pandemic evolution with the GLEaM model. The method shifts the burden of estimating the disease transmissibility from the incidence data, suffering notification/surveillance biases and dependent on country specific surveillance systems, to the more accurate data of the early case detection in newly affected countries. This is achieved through the modeling of human mobility patterns on the global level obtained from high quality databases. In other words, the chronology of the infection of new countries is determined by two factors. The first is the number of cases generated by the epidemic in the originating country. The second is the mobility of people from this country to the rest of the world. The mobility data are defined from the outset with great accuracy and we can therefore find the parameters of the disease spreading as those that provide the best fit for the time of infection of new countries. This method also allows for uncovering the presence of a seasonal signature in the observed pattern, not hindered or effectively caused by notification and reporting changes in each country's influenza monitoring. The obtained values for the reproduction numbers are larger than the early estimates [5], though aligned with later works [11–13]. The simulated geographic and temporal evolution of the pandemic based on these estimates shows the possibility of an early epidemic activity peak in the Northern hemisphere as soon as mid October. While the simulations refer to a worst-case scenario, with no intervention implemented, the present findings pertain to the timing of the vaccination campaigns as planned by many countries. For this reason we also present an analysis of scenarios in which the systematic use of antiviral drug therapy is implemented with varying effectiveness, according to the national stockpiles, and study their effect on the epidemic timeline.

The GLEaM structured metapopulation model is based on a metapopulation approach [4, 14–22] in which the world is divided into geographical regions defining a subpopulation network where connections among subpopulations represent the individual fluxes due to the transportation and mobility infrastructure. GLEaM integrates three different data layers [10]. The population layer is based on the high-resolution population database of the 'Gridded Population of the World' project of the SocioEconomic Data and Applications Center (SEDAC) [23] that estimates the population with a granularity given by a lattice of cells covering the whole planet at a resolution of 15 × 15 minutes of arc. The transportation mobility layer integrates air travel mobility obtained from the International Air Transport Association (IATA [24]) and Official Airline Guide (OAG [25]) databases that contain the list of worldwide airport pairs connected by direct flights and the number of available seats on any given connection [26]. The combination of the population and mobility layers allows the subdivision of the world into georeferenced census areas defined with a Voronoi tessellation procedure [27] around transportation hubs. These census areas define the subpopulations of the metapopulation modeling structure (see Figure 1). In particular, we identify 3,362 subpopulations centered around IATA airports in 220 different countries (see [10] and Additional file 1 for more details). GLEaM integrates short scale mobility between adjacent subpopulations by considering commuting patterns worldwide as obtained from the data collected and analyzed from more than 29 countries in 5 continents across the world [10]. Superimposed on these layers is the epidemic layer that defines the disease and population dynamics. The model simulates the mobility of individuals from one subpopulation to another by a stochastic procedure in which the number of passengers of each compartment traveling from a subpopulation j to a subpopulation l is an integer random variable defined by the actual data from the airline transportation database (see Additional file 1). Short range commuting between subpopulations is modeled with a time scale separation approach that defines the effective force of infections in connected subpopulations [10, 28, 29]. The infection dynamics takes place within each subpopulation and assumes the classic influenza-like illness compartmentalization in which each individual is classified by a discrete state such as susceptible, latent, infectious symptomatic, infectious asymptomatic or permanently recovered/removed [9, 30]. The model therefore assumes that the latent period is equivalent to the incubation period and that no secondary transmissions occur during the incubation period (see Figure 1 for a detailed description of the compartmentalization). All transitions are modeled through binomial and multinomial processes to preserve the discrete and stochastic nature of the processes (see Additional file 1 for the full description). Asymptomatic individuals are considered as a fraction p _a = 33% of the infectious individuals [31] generated in the model and assumed to infect with a relative infectiousness of r _β = 50% [5, 30, 32]. Change in traveling behavior after the onset of symptoms is modeled with the probability 1 - p _t, set to 50%, that individuals would stop traveling when ill [30]. The spreading rate of the disease is ultimately governed by the basic reproduction number R ₀. Once the disease parameters and initial conditions based on available data are defined, GLEaM allows the generation of stochastic realizations of the worldwide unfolding of the epidemic, with mobility processes entirely based on real data. The model generates in silico epidemics for which we can gather information such as prevalence, morbidity, number of secondary cases, number of imported cases and other quantities for each subpopulation and with a time resolution of 1 day. While global models are generally used to produce scenarios in which the basic disease parameters are defined from the outset, here we use the model to provide a maximum likelihood estimate of the transmission potential by finding the set of disease parameters that best fit the data on the arrival time of cases in different countries worldwide. It is important to stress that the model is not an agent-based model and does not include additional structure within a subpopulation, therefore it cannot provide detailed information at the level of households or workplaces. The projections for the winter season in the northern hemisphere are also assuming that there will be no mutation of the virus with respect to the spring/summer of 2009. Furthermore, while at the moment the novel H1N1 influenza is accounting for 75% of the influenza cases worldwide, the model does not consider the cocirculation of different influenza strains and cannot provide information on cocirculation data.

The initial conditions of the epidemic are defined by setting the onset of the outbreak near La Gloria in Mexico on 18 February 2009, as reported by official sources [33] and analogously to other works [5]. We tested different localizations of the first cases in census areas close to La Gloria without observing relevant variations with respect to the observed results. We also performed sensitivity analysis on the starting date by selecting a seeding date anticipated or delayed by 1 week with respect to the date available in official reports [33]. The arrival time of infected individuals in the countries seeded by Mexico is clearly a combination of the number of cases present in the originating country (Mexico) and the mobility network, both within Mexico and connecting Mexico with countries abroad. For this reason we integrated into our model the data on Mexico-US border commuting (see Figure 2a), which could be relevant in defining the importation of cases in the US, along with Mexican internal commuting patterns (see Figure 1) that are responsible for the diffusion of the disease from rural areas as La Gloria to transportation hubs such as Mexico City. In addition, we used a time-dependent modification of the reproductive number in Mexico as in [6] to model the control measures implemented in the country starting 24 April and ending 10 May, as those might affect the spread to other countries.

In order to ascertain the effect of seasonality on the observed pattern, we explored different seasonality schemes. The seasonality is modeled by a standard forcing that rescales the value of the basic reproductive number into a seasonally rescaled reproductive number, R(t), depending on time. The seasonal rescaling is time and location dependent by means of a scaling multiplicative factor generated by a sinusoidal function with a total period of 12 months oscillating in the range α_min to α_max, with α_max = 1.1 days (sensitivity analysis in the range 1.0 to 1.1) and α_min a free parameter to be estimated [17]. The rescaling function is in opposition in the Northern and Southern hemispheres (see Additional file 1 for details). No rescaling is assumed in the Tropics. The value of R ₀ reported in the Tables and the definition of the baseline is the reference value in the Tropics. In each subpopulation the R(t) relative to the corresponding geographical location and time of the year is used in the simulations.

The major problem in the case of projections on an extended time horizon is the seasonality effect that in the long run is crucial in determining the peak of the epidemic. In order to quantify the degree of seasonality observed in the current epidemic, we estimate the minimum seasonality scaling factor α_min of the sinusoidal forcing by extending the chronology under study and analyzing the whole data set composed of the arrival dates of the first infected case in the 93 countries affected by the outbreak as of 18 June. We studied the correlation between the simulated arrival time by country and its corresponding empirical value, by measuring the regression coefficient between the two datasets. Given the extended time frame under observation, the arrival times considered in this case are expected to provide a signature of the presence of seasonality. They included the seeding of new countries from outbreaks taking place in regions where seasonal effects might occur, as for example in the US or in the UK. For the simulated arrival times we have considered the median and 95% confidence interval (CI) emerging from the 2 × 10³ stochastic runs. The regression coefficient is found to be sensitive to variations in the seasonality scaling factor, allowing discrimination of the α_min value that best fits the real epidemic. A detailed presentation of this analysis is reported in Additional file 1. The full exploration of the phase space of epidemic parameters and seasonality scenarios reported in Additional file 1 required data from 10⁶ simulations; the equivalent of 2 million minutes of PowerPC 970 2.5 GHz CPU time.

Table 1 reports the results of the maximum likelihood procedure and of the correlation analysis on the arrival times for the estimation of α_min. In the following we consider as the baseline case the set of parameters defined by the best estimates: G _t = 3.6 days, μ^-1 = 2.5 days, R ₀ = 1.75.

The empirical arrival time data used for the likelihood analysis are necessarily an overestimation of the actual date of the importation of cases as cases could go undetected. If we assume a shift of 7 days earlier for all arrival times available from official reports, the resulting maximum likelihood is increasing the best estimate for R ₀ to 1.87 (95% CI 1.73 to 2.01), as expected since earlier case importation necessitates a larger growth rate of the epidemic. The official timeline used here therefore provides, all other parameters being equal, a lower estimate of the transmission potential. We have also explored the use of a subset of the 12 countries, always generating results within the confidence interval of the best estimate.

The best estimates reported in Table 1 do not show any observable dependence on the assumption about the seasonality scenario (as reported in Additional file 1). The analysis is restricted to the first countries seeded from Mexico to preserve the conditional independence of the variables and it is natural to see the lack of any seasonal signature since these countries receive the disease from a single country, mostly found in the tropical region where no seasonal effects are expected.

In order to find the minimum seasonality scaling factor α_min that best fits the empirical data, we performed a statistical correlation analysis of the arrival time of the infection in the 93 countries infected as of 18 June, as detailed in the Methods section and Additional file 1. By considering a larger number of countries and a longer period for the unfolding of the epidemic worldwide as seasons change, the correlation analysis for the baseline scenario provides clear statistical indications for a minimum rescaling factor in the interval 0.6 < α_min < 0.7. In the full range of epidemic parameters explored, the correlation analysis yields values for α_min in the range 0.4 to 0.9. This evidence for a mild seasonality rescaling is consistent with the activity observed in the months of June and July in Europe and the US where the epidemic progression has not stopped and the number of cases keeps increasing considerably (see also Table 2 for the corresponding values of R(t) in those regions during summer months).

This analysis allows us to provide a comparison with the epidemic activity observed so far, and most importantly an early assessment of the future unfolding of the epidemics. For each set of parameters the model generates quantities of interest such as the profile of the epidemic behavior in each subpopulation or the number of imported cases. Each simulation generates a stochastic realization of the process and the curves are the statistical aggregate of at least 2 × 10³ realizations. In the following we report the median profiles and where indicated the 95% CI. For the sake of clarity data are aggregated at the level of country or geographical region. Additional file 1 reports a detailed comparison of the simulated number of cases in Australia, US, UK with the reported cases from official sources in the period May to July. Results are in good agreement with the reported temporal evolution of the epidemic and highlight a progressive decrease of the monitoring activity caused by the increasing number of cases, as expected [7]. The same information is also available for each single subpopulation defined in the model. We have therefore tested the model results in four territories of Australia. Interestingly, the model is able to recover the different timing observed in the four territories. A detailed discussion of this comparison is reported in Additional file 1.

The number of hospitalized at peak times in the selected countries range between 2 and 40 per 100,000 persons, for a hospitalization rate typical of seasonal influenza and for an assumed 1% rate, respectively, yielding a quantitative indication of the potential burden that the health care systems will likely face at the peak of the epidemic activity in the next few months. It is worth noting that the present analysis considers a worst-case scenario in which no effective containment measures are introduced. This is surely not the case in that pandemic plans and mitigation strategies are considered at the national and international level. Guidelines aimed at increasing social distancing and the isolation of cases will be crucial in trying to mitigate and delay the spread in the community, thus reducing the overwhelming requests on the hospital systems. Most importantly, the mass vaccination of a large fraction of the population would strongly alter the presented picture. By contrast, any mass vaccination campaign is unlikely to start before the middle of October [45, 46]. The potential for an early activity peak of the pandemic in October/November puts at risk the effectiveness of any mass vaccination program that might take place too late with respect to the pandemic wave in the Northern hemisphere. In this case it is natural to imagine the use of other mitigation strategies aimed at delaying the activity peak so that the maximum benefit can be gained with the vaccination program. As an example, we studied the implementation of systematic antiviral (AV) treatment and its effect in delaying the activity peak [19, 30, 32, 39, 47–50]. The resulting effects are clearly country specific in that each country will experience a different timing for the epidemic peak (with a local transmissibility increasing in value as we approach the winter months) and will count on antiviral stockpiles of different sizes. Here we consider the implementation of the AV treatment in all countries in the world that have drugs stockpiles available (source data from [51, 52] and national agencies), until the exhaustion of their stockpiles [4]. We have modeled this mitigation policy with a conservative therapeutic successful use of drugs for 30% of symptomatic infectious individuals. The efficacy of the AV is accounted in the model by a 62% reduction in the transmissibility of the disease of an infected person under AV treatment when AV drugs are administered in a timely fashion [30, 32]. We assume that the drugs are administered within 1 day of the onset of symptoms. We also consider that the AV treatment reduces the infectious period by 1 day [30, 32]. In Figure 3 we show the delay obtained with the implementation of the AV treatment protocol in a subset of countries with available stockpiles. As an example, we also show the incidence profiles for the cases of Spain and Germany, where it is possible to achieve a delay of about 4 weeks with the use of 5 million and 10 million courses of AV, respectively. The results of this mitigation might be extremely valuable in providing the necessary time for the implementation of the mass vaccination program.

We have defined a Monte Carlo likelihood analysis for the assessment of the seasonal transmission potential of the new A(H1N1) influenza based on the analysis of the chronology of case detection in affected countries at the early stage of the epidemic. This method allows the use of data coming from the border controls and the enhanced surveillance aimed at detecting the first cases reaching uninfected countries. This data is, in principle, more reliable than the raw count of cases provided by countries during the evolution of the epidemic. The procedure provides the necessary input to the large-scale computational model for the analysis of the unfolding of the pandemic in the future months. The analysis shows the potential for an early activity peak that strongly emphasizes the need for detailed planning for additional intervention measures, such as social distancing and antiviral drugs use, to delay the epidemic activity peak and thus increase the effectiveness of the subsequent vaccination effort.