Forecasting seasonal outbreaks of influenza

Results

Before using the GFT estimates we first validated the functionality of the combined SIRS-EAKF framework, using a synthetic, model-generated record. This preliminary validation of the combined SIRS-EAKF framework with a wholly known synthetic record was performed to determine whether the EAKF appropriately constrains the SIRS model state and parameters.

To conduct the validation, a single run of the humidity-driven SIRS model was initiated with prescribed parameters and initial conditions and integrated from October 1, 1972 into the spring of 1973, using daily absolute humidity (AH) conditions for New York state. This simulation produced a January outbreak of influenza (Fig. 1A), which we defined as the “true” outbreak. Prescribed observational error was then added to this truth to create a synthetic, error-laden observational record of influenza infection to be assimilated in the combined SIRS-EAKF framework.

For this first assimilation, an ensemble of 100 simulations was initiated, each with randomly chosen parameters and initial conditions, and integrated through time using the same daily AH record for New York state used to generate the synthetic truth. The synthetic observations were then iteratively assimilated using the EAKF as the model moved forward in time. Ensemble posterior mean estimates of the true state and parameters were well constrained by the combined SIRS-EAKF framework (Fig. 1A).

The observed state variable, I, the number of infected persons, is well captured in the mean by the combined SIRS-EAKF framework. The nonobserved state variable, S, the number of susceptible persons, is also well captured, such that before the cessation of the outbreak, there is little bias between the ensemble mean and the truth. Through time the variance of the 100-member ensemble decreases continuously for S (Fig. 1B). For I, error growth and ensemble spread occur before the outbreak as individual ensemble members engender premature outbreaks, only to be constrained by the observations. During the outbreak, I undergoes more intense error growth as outbreaks are then supported by the observations.

Epidemiologically significant parameters within the SIRS model, including the mean infection period, D, and R_0max, which sets an upper bound on the transmissibility of the virus, are also constrained by the EAKF. By the end of the true outbreak, the ensemble mean estimates of these parameters are close to the true values used to generate the synthetic observations.

For comparison with the EAKF approach, we alternatively tested a particle filter method, using the same SIRS model and synthetic observations (SI Methods). Unlike the EAKF, particle filters make no assumptions on the linearity underlying the model or the mapping from state to observational space. However, particle filter methods are prone to particle degeneracy, which can lead to model divergence and poor state estimation (26, 27). The particle filter method tested performed reasonably well; however, the EAKF method provided better and more computationally efficient estimates of the unobserved state, S, and parameters (Fig. S1). This better performance may stem from the fact that the EAKF approach maintains greater spread during the preoutbreak period, which may enable better state and parameter estimation once the outbreak commences.

We have also elected to present the EAKF as the method of data assimilation because the EAKF can be readily applied to higher-dimensional systems, unlike particle filter methods (28). Whereas the current system is only dimension 6 (two state variables and four parameters), we anticipate use of higher-dimensional influenza model systems that include multiple strains, age stratification, and geographic structure in the future.

Additional simulations were performed to gauge the sensitivity of the combined SIRS-EAKF framework to changes of the ensemble size, the random sampling of initial ensemble member states, the frequency of observations (e.g., every second day, weekly, etc.), and the prescribed error variance used to generate the synthetic observations from the truth (Figs. S2–S4). These sensitivity tests indicated that the SIRS-EAKF framework provides a robust estimate of the truth provided the ensemble is of sufficient size (100 or more members), the observations are available with sufficient frequency (at least every 10 d), and the observational error variance is not too large.

We next applied the combined SIRS-EAKF framework to perform retrospective simulation of influenza infection in New York City during the 2004–2005 and 2007–2008 seasons, using GFT weekly estimates of ILI (Figs. S5–S7). A small level of multiplicative inflation was included in these simulations to avoid filter divergence (29, 30) and better replicate the GFT observations (SI Methods). This need for inflation, which is commonly used during data assimilation, in part stems from the fact that the GFT observations reflect real-world transmission dynamics, which include structures and processes, such as spatial heterogeneity and preferential mixing, that are not represented in the SIRS model and contribute to model error.

Two hundred-member ensembles of the SIRS-EAKF framework were run with observed daily AH for New York City as forcing. Each ensemble member was initialized with a unique sample of state variables and SIRS parameters in September before the wintertime influenza season. The ensemble was then integrated through time and GFT ILI observations were assimilated using the EAKF to generate a posterior estimate of infection rates. This process was repeated 250 times, each time with a different 200-member ensemble of initial parameter and state variable combinations. With a small amount of inflation, the average ensemble mean posterior estimates match the GFT weekly estimates of ILI well (Fig. 2). The complicated outbreak structure during the 2004–2005 influenza season, which includes three individual peaks, is represented by the SIRS-EAKF framework. The model–data assimilation system also captures the long rise and single peak of infection during 2007–2008.

A retrospective forecast for the 2007–2008 influenza season in New York City was next generated. As in preceding simulations, after the assimilation of a weekly GFT observation, the model posterior was integrated forward 7 d to the next GFT observation; however, in addition, the same model posterior was also integrated forward 300 d without any further EAKF constraint, although with perfect knowledge of the daily AH conditions. These second longer runs constitute retrospective forecasts; i.e., after assimilation of a GFT ILI estimate in near-real time and consequent adjustment of SIRS ensemble member parameters and state variables by the EAKF process, the model is integrated into the future, using those new parameters and the current state, to generate anticipated possible outbreak outcomes.

The SIRS-EAKF system can be evaluated for its predictive skill for any number of outcomes, such as outbreak magnitude, duration, or onset. Here we choose to focus on one metric: the timing of the peak of the outbreak—i.e., the week during the influenza season in which the highest GFT ILI estimate is recorded (e.g., week 25, or the week ending February 17 for the 2007–2008 season). The retrospective forecasts indicate that model predictions of influenza outbreak peak converge toward the observed peak in the weeks before the actual event (Fig. 3 A and B). In fact, the ensemble mode predicted peak is within 1 wk of being correct 5 wk in advance of the observed peak. In addition, the spread of peak predictions among the 200-member ensemble decreases as the forecast nears the observed peak. This decrease in ensemble spread can be seen in weekly histograms of predicted peak timing, as well as in the ensemble variance of these same predictions.

The percentage of ensemble members predicting the peak timing within ±1 wk increases dramatically beginning in week 18 as model spread decreases (Fig. 3 C and D). This relationship indicates that, for a real-time prediction performed in the absence of knowledge of when the actual observed peak will occur, a measure of confidence in the prediction can be made on the basis of the ensemble spread. That is, decreased ensemble spread will coincide with increased confidence in the prediction.

To further codify this finding and verify its generality across a larger set of seasonal outbreaks, we next ran repeated retrospective ensemble forecasts for the 2003–2004 through 2008–2009 influenza seasons in New York City. These forecasts were run in 200-member ensembles, using two different approaches: (i) with perfect knowledge of future daily AH conditions for only the first 5 d of the influenza forecast—these conditions would be locally available from numerical weather forecasts—and daily climatological AH conditions for the remaining 295 d of each forecast; and (ii) with perfect knowledge of future daily AH conditions for the first 5 d of the forecast followed by 295 d of AH climatology, but with only the model infected state constrained by the EAKF, i.e., the posterior parameters and other state variable are reset to an initial distribution before commencing each forecast.

The first approach is how a real-time forecast could be run for this SIRS-EAKF framework: with full use of the EAKF constraint and only limited knowledge of future AH conditions. The difference between the outcomes of the first and second approaches indicates the forecast improvement due to full EAKF constraint of all variables and parameters.

The same relationship between the accuracy of the ensemble mode predicted peak and the spread of the predictions within the ensemble holds for these more numerous retrospective forecasts (Fig. S8). A clearer picture of this predictive skill emerges when the forecasts are grouped by how many weeks into the future the ensemble mode peak is predicted (SI Methods). The two different forecast approaches have comparable predictive skill when the ensemble peak mode is predicted to be 1–3 wk in the future (Fig. 4 and Fig. S9); however, for ensemble peak modes forecast greater than 3 wk in the future, only the first forecast form, which makes full use of the EAKF state variable and parameter constraint, has demonstrated skill (Fig. 4). The second forecast form, which lacks full EAKF constraint, has no forecast skill for mode peak predictions greater than 3 wk; indeed, without full EAKF constraint mode peak predictions rarely occur more than 3 wk in the future.

Results from the first forecast form, which uses climatological AH conditions and full EAKF constraint of model variables and parameters, can now be used to explore how a real-time influenza forecast would be interpreted (Fig. 4, red lines). On the basis of the ensemble spread and how many weeks in the future the outbreak peak is predicted, confidence in a given prediction, or, conversely, the forecast uncertainty, can be assigned on the basis of the probabilities shown in Fig. 4. For instance, an ensemble mode peak predicted for 7 wk in the future with a log-transformed ensemble variance of 2.5 wk squared is expected to be accurate within ±1 wk about 40% of the time.

Discussion

The results presented here demonstrate that weekly, local predictions of influenza risk, with estimates of forecast certainty, can be made in real time. These forecasts can be performed using GFT ILI data at municipal, state, and country levels and used to help inform public health decisions including vaccine allocation and antiviral drug distribution.

In the future, as more years of GFT ILI estimate data and more locations are entrained into this analysis, a more robust estimate of the relationship between predictions of peak timing and ensemble variance will be developed. This will further codify the relationships presented in Figs. 3 and 4 and improve prediction skill. Forecasts can also be developed for other outcome metrics, such as peak magnitude and outbreak duration. Predictions of large local maxima, rather than the one seasonal peak we have used here, might also prove useful; such a metric may improve prediction skill, as some years (e.g., 2005 in New York City) have no distinct, substantive peak. Additional information, such as influenza strain or school calendar, might also be used to further optimize predictions.

Information on influenza strain may be particularly important. The three peaks during 2004–2005 (Fig. 2A) stem principally from differently timed outbreaks of three distinct influenza subtypes in circulation that season. The SIRS model, as run here with only one strain, is not equipped to produce separate parameter and state estimates for each subtype, which would be a truer representation of actual seasonal dynamics. A model that accounts for and assimilates infection data by strain type would likely improve prediction skill and uncertainty estimates for each strain and the season as a whole; however, at present, such partitioned data are not available in real time. Alternatively, a segregated analysis of predictions from seasons with one predominant strain, vs. those with more than one strain, might improve prediction and uncertainty estimates for those single-strain seasons, simply by aggregating those “cleaner” years.

In weather and climate, forecasts made with a suite of models are generally more robust than forecasts made with any single model (31, 32). This finding motivates the development and use of additional model–data assimilation influenza forecast frameworks to be used in conjunction with the approach presented here. Other models might be spatially distributed, discrete, or age stratified or include stochasticity. Alternate assimilation approaches should also be tested further, including additional particle filtering methods (18, 26, 27) or alternate ensemble Kalman filtering forms (33, 34). Indeed, no one model or assimilation method need be used to the exclusion of all others, and it is not our intention to advocate exclusive use of either the SIRS model or the EAKF.

The general model–data assimilation framework presented here is flexible and adaptive and could be used to develop predictions for other seasonally recurring respiratory diseases, such as respiratory syncytial virus and rhinovirus, should real-time estimates of these pathogens become available. The forecasts developed here indicate that we will soon reach an era when reliable forecasts of some infectious pathogens are as commonplace as weather predictions.

Methods

Description of the SIRS Model.

The model used for this study is a perfectly mixed, absolute humidity-driven SIRS construct. The form is similar to that of the model described in Shaman et al. (3). The SIRS model equations are

where S is the number of susceptible people in the population, t is time in years, N is the population size, I is the number of infectious people, N − S − I is the number of resistant individuals, β(t) is the contact rate at time t, L is the average duration of immunity, D is the mean infectious period, and α is the rate of travel-related import of influenza virus into the model domain. The basic reproductive number, which is the number of secondary infections the average infectious person would produce in a fully susceptible population, at time t is related to the contact rate through the expression .

AH modulates transmission rates within this model by altering R₀(t) through an exponential relationship similar to how AH has been shown to affect both influenza virus survival and transmission in laboratory experiments (2),

where a = −180, , R_0max is the maximum daily basic reproductive number, R_0min is the minimum daily basic reproductive number, and q(t) is the time-varying specific humidity, a measure of AH. The value of a is estimated from the laboratory regression of influenza virus survival upon AH (2).

This SIRS model, as run previously, reproduces both the observed 1972–2002 seasonal cycle of excess pneumonia and influenza mortality within the United States and the association between negative AH anomalies and the onset of local influenza outbreaks (3). It is a strong starting point for the proposed project because it has been validated against datasets of historical influenza outbreaks, and its inclusion of environmental forcing (i.e., AH conditions) may help the model better simulate outbreak dynamics.

For this study, the model was run continuously and without stochasticity. That is, fractional persons were generated within model state classes, and transitions between model states were calculated deterministically directly from Eqs. 1 and 2. The complete temporal evolution of the SIRS model is fully described by the model equations (Eqs. 1–3), initial conditions, and AH boundary forcing. Simulations were performed with an influenza importation rate, α, of 0.1 infections per day (1 infection every 10 d). Only one influenza strain was simulated per season, and partial and cross-immunities were not represented.

Further description of the model initialization and forcing for individual SIRS-EAKF experiments is provided below.

Description of the EAKF.

Generally, sequential ensemble filtering can be viewed as the problem of estimating the probability of the system state at a given time (Z_t) conditional on observations (y_t) taken up until and including time t. In the case of the SIRS model described above, the state is composed of the number of individuals in a population that are susceptible and infectious at a given time and the set of independent model parameters; that is, Z_t = (S_t, I_t, R_0max, R_0min, L, D) and y_t is a time series representing observations (e.g., the Google Flu Trends ILI estimates).

Bayes’ rule provides a target for the update of the system state given an observation:

Here the first term on the right-hand side is the likelihood of observing the data given the state and the second term is the prior distribution of the system state. The updated distribution (the left-hand side of Eq. 4) is called the posterior. Generally, Kalman filters fall into a class of filters that assumes normality of both the likelihood and prior distributions during an update. This assumption allows parameterization of these distributions in terms of the first two moments only (mean and covariance). However, in ensemble filtering, one need never form the covariance of the prior or posterior; instead we have a finite ensemble of states that are samples from these distributions. The model mean and covariance are computed directly from the ensemble.

Knowledge of these prior moments, as well as observations and their error, allows for the computation of the mean and covariance of the posterior. The method of transforming the ensemble members of the prior into ensemble members of the posterior is what distinguishes different types of ensemble filters. The EAKF [see Anderson (20) for algorithmic details] adjusts the prior ensemble members such that their new moments match the target moments of the posterior predicted by Bayes’ theorem. The method is sequential in that given samples from the posterior, the SIRS model can be used to integrate each ensemble member forward in time to the point at which new observations become available. The update is then repeated.

In models that consist of multiple prognostic variables—in this instance, the state variables S and I, as well as the model parameters—covariant relationships between variables arise naturally from the dynamics of the system. Mathematically rigorous methods for assimilating real-world data into numerical models rely on knowledge of these intervariable relationships. This information allows “balanced” adjustments to the entire model solution, even when some of the prognostic variables, such as S, and the parameters are not directly observable quantities. Ensemble filters store all of the information about the state variable and model parameter interrelationships in the form of multiple model solutions (the ensemble), all of which are possible realizations given past measurements.

In Kalman filtering these intervariable relationships are assumed linear (i.e., jointly distributed by a multivariate Gaussian). The EAKF inherits this assumption, but the adjustment operates only on the first two moments of the prior distribution, leaving the higher-order moments unchanged. In the case of a system in which intervariable relationships are nearly linear (such as the SIRS model), the ensemble filter can have high utility even when it is not strictly optimal. Ensemble Kalman filtering (like particle filtering) also has the attractive quality that the algorithm can be implemented independently of the dynamic model. In contrast to variational methods, which require adjoints or parametric covariance forms, ensemble methods use intervariable relationships that are statistically diagnosed from the ensemble of model solutions.

Supplementary Material