Forecasting the spatial transmission of influenza in the United States

Metapopulation Model

To forecast the spatial diffusion of ILI+, we developed a metapopulation model that can flexibly generate patterns of spatial transmission. We used ILI+ resolved at the state level. Within the model, the different states are connected by mobility flows (11), including movement by recurrent commuters and random visitors. The numbers of work commuters crossing state lines are available from the US census survey (23). For example, for the six southeastern US states, most commuting events occur between contiguous states, while long-range movements are quite rare (Fig. 1C and SI Appendix, Analysis of the Commuting Pattern Across States and Fig. S1). Using the commuting data, the total population in $l$ states is subdivided into $l^2$ subpopulations ${\{N_n^k\}}_{l\times l}$, where $N_n^k$ represents the fixed subpopulation living in state $k$ and commuting to work in state $n$. In constructing our metapopulation model, we aim to track the evolution of the susceptible ($S_n^k$) and infected ($I_n^k$) population in each $N_n^k$, and simulate influenza transmission using a humidity-driven susceptible–infected–recovered–susceptible (SIRS) model (24, 25) (SI Appendix, Humidity-Driven SIRS Model).

Because the population composition in each state differs during day and night time, we separate the transmission equations into two parts. During daytime, subpopulation $N_n^k$ is present in state $n$ and participates in transmission there. The total, susceptible, and infected populations in state $n$ during daytime are simply $N_n={\sum }_r{N}_n^r$, $S_n={\sum }_r{S}_n^r$, and $I_n={\sum }_r{I}_n^r$, respectively.

In addition to fixed commuting flows, we also simulate the irregular movement of visitors who travel for reasons other than work: business trips or vacations, for example. These individuals circulate among subpopulations following a Markov process, causing a population exchange in all subpopulations ${\{N_n^k\}}_{l\times l}$. Transportation records of flights, trains, and buses exist; however, this information is typically released after a long lag period and varies in availability across different locations and data sources. Due to a lack of real-time information on this movement, we assume the number of random visitors traveling from state $m$ to $n$ is proportional to the average number of commuters between them: ${\overline{N}}_n^m=(N_n^m+N_m^n)/2$. In particular, the random visitors from state $m$ to $n$ are assumed to be $\theta {\overline{N}}_n^m$, where $\theta$ is an adjustable parameter. The rationale behind this proportionality assumption is that, if two locations exchange a larger number of commuters, there should exist a stronger economic tie or other incentives to facilitate the exchange of random visitors.

During daytime, $\theta {\overline{N}}_n^m$ persons, drawn uniformly from the population present in state $m$, move to state $n$ and are randomly redistributed into the subpopulations there. Such population exchange exists for all pairs of states. This implies that the infected individuals leaving subpopulation $N_n^k$ are $\theta {\sum }_{m\ne n}{\overline{N}}_m^n\times I_n^k/N_n$, where $\theta {\sum }_{m\ne n}{\overline{N}}_m^n$ is the total number of random visitors leaving state $n$, and $I_n^k/N_n$ is the ratio of the number of infected people from subpopulation $N_n^k$ to the total population of state $n$. Similarly, the infected individuals entering $N_n^k$ can be calculated as $\theta {\sum }_{m\ne n}{\overline{N}}_n^m{I}_m/N_m\times N_n^k/N_n$, where $\theta {\sum }_{m\ne n}{\overline{N}}_n^m{I}_m/N_m$ is the total number of infected persons entering state $n$, and $N_n^k/N_n$ is the ratio of subpopulation $N_n^k$ to the total population of state $n$. The random exchange of susceptible individuals can be calculated similarly. Combining the regular commuters and random visitors, the transmission equations for $S_n^k$ and $I_n^k$ during daytime are

Here ${\beta }_n(t)$ is the contact rate at time $t$ at state $n$, $D$ is the duration of infection, and $L$ is the duration of immunity. The transmission equations at nighttime are formulated similarly (SI Appendix, Metapopulation Model). During implementation, we assume daytime and nighttime transmissions last for 8 and 16 h, respectively. That is, the model is integrated deterministically using the daytime equations for an interval of 1/3 d, and then the nighttime equations are integrated for the subsequent 2/3 d. The number of new infections in subpopulation $N_n^k$ is calculated using the contact terms determined through model integration. Observations at location $k$, i.e., weekly incidence, are obtained by summing the new infections in subpopulations whose living location is $k$.

Parameter Inference

The metapopulation model is capable of generating abundant scenarios of spatial spread from variable ratios of work commuters and random visitors by adjusting the parameter $\theta$. This relatively parsimonious construct allows its application in conjunction with statistical filtering techniques, which, in combination, can be used to infer model state space and parameters and to generate ensemble forecasts (SI Appendix, Model–Data Assimilation Framework).

Similar model–data assimilation frameworks have been successfully used to forecast the local epidemic trajectory for a number of infectious diseases (13–16, 26). During model training, the state variables and parameters in the dynamical model are repeatedly calibrated by available observations; a forecast is then generated by integrating the optimized model into the future. In practice, an ensemble of randomly generated epidemic trajectories is simulated to produce a probabilistic distribution of forecast outcomes.

The metapopulation model possesses a high-dimensional state vector. For each subpopulation $N_n^k$, we record three variables: $S$, $I$, and weekly incidence. The parameters $R_{0max}$ (the maximal basic reproductive number), $R_{0min}$ (the minimal basic reproductive number), $L$ (duration of immunity), $D$ (length of infectious period), and $\theta$ (random movement ratio) are shared for all subpopulations. For metapopulation forecasts at a high geographical resolution, the dimension of state vectors will undergo a quadratic growth as the number of locations increases. To deal with this high-dimensional data assimilation problem, we chose an efficient filtering algorithm: the Ensemble Adjustment Kalman Filter (EAKF) (31). Unlike particle filters (32), which require a larger number of particles for high-dimensional problems and are more suitable for low-dimension systems (33), the EAKF uses a limited number of ensemble members and is conducive for use with high-dimension systems. Previous forecast work has shown that the performances of particle filters and the EAKF are similar when used to produce influenza forecasts in conjunction with a SIRS model (13). In Fig. 1B, the posterior fitting of the metapopulation model–EAKF to the southeastern US states during the 2009 pandemic is presented. The epidemic curves are well captured by the model–data assimilation system.

To validate the forecast system further, we next examined whether key model parameters can be accurately recovered from model-generated synthetic outbreaks with known true values (i.e., the “truth”). Two such outbreaks were generated using the metapopulation model for the six southeastern US states and the same initial conditions and parameters excepting a change of parameter $\theta$ (Fig. 2 A and B). To mimic observational error, we added noise to the model-generated time series of influenza incidence. The prescribed observational error variance, OEV, at week $t$ is defined as ${\text{OEV}}_t=1\times {10}^5+{(\text{average ILI}+\text{in preceding}3\text{wk})}^2/5$. We then used the metapopulation model–EAKF framework, in conjunction with these error-laden observations, to estimate the model parameters. In Fig. 2 C–F, we report the posterior mean of two of the more sensitive parameters, $R_{0max}$ (the maximal basic reproductive number) and $\theta$ (the random movement ratio), as inferred by the metapopulation model–EAKF system. For both outbreak simulations, the metapopulation model estimates of these parameters, inferred using observations from all six states, are rapidly adjusted toward the truth. Additional analyses indicate that inference accuracy is not particularly sensitive to the initial conditions of the model (SI Appendix, Inference of Parameters and Variables and Figs. S2–S4). Further, the metapopulation forecast system generates more-accurate forecasts of onset week, peak timing, and peak intensity for the synthetic truth, compared with a reference SIRS–EAKF system in which each state is modeled in isolation (13, 14) (SI Appendix, Forecasts of Synthetic Outbreaks and Fig. S5).

In these synthetic tests, the metapopulation model was run deterministically. For infectious diseases for which only sporadic case numbers are observed, a deterministic model could introduce significant errors by rounding the real number estimates of infection cases to integer values, and a stochastic model may be preferable; however, the deterministic model is suitable in this study of seasonal influenza, as large numbers of infections exist. Indeed, we also ran a stochastic version of the metapopulation model in synthetic tests and found no significant change in forecast accuracy (SI Appendix, Fig. S6).

Retrospective Forecasts

We next performed retrospective forecasts for 35 states in the continental United States for the 2008–2009 through 2012–2013 influenza seasons, using the AFHSB type A ILI+ data. The remaining 15 states were excluded, as they were not well represented in the AFHSB dataset due to low numbers of ILI records (SI Appendix, Fig. S7). The remaining AFHSB ILI+ data provide reasonable time series representations of influenza incidence at national and state levels (Fig. 3A and SI Appendix, Fig. S8). We also performed synthetic tests using the full 35-state system, which showed that key parameters, such as $R_{0max}$ and $\theta$, can still be accurately inferred (SI Appendix, Fig. S9).

For each season, we generated 100 independent ensemble forecasts for each week during the flu season using both the metapopulation model and the isolated SIRS model (SI Appendix, Retrospective Forecasts of Historical Outbreaks and Figs. S10–S12). As the actual onset and peak weeks are unknown in real time, we evaluated forecast accuracy relative to the predicted onset or peak week. Specifically, forecast accuracy for onset and peak week (peak intensity) is the fraction of ensemble mean predictions within $\pm$1 wk ($\pm$25%) from the observed. Fig. 3B shows average forecast accuracy for onset week, peak week, and peak intensity over all seasons and states. Here, a positive predicted lead indicates that onset or peak week is predicted to occur in the future, while a negative value implies that the onset or peak week is predicted to have already passed. Onset week was predicted up to 6 wk in advance by the metapopulation system, whereas the isolated forecasts failed to predict most onset events until a predicted lead of 2 wk. In Fig. 3C, we summarize the performance of the metapopulation forecasts of onset week compared with the isolated forecasts for the individual states. Onset predictions for a few states were degraded (particularly at shorter leads); however, the majority of states had improved onset forecast accuracy with the metapopulation model. Sensitivity analysis showed that this improvement held for alternate definitions of onset (SI Appendix, Fig. S13). In addition, the metapopulation forecast system was more accurate in predicting the peak week and peak intensity for all prospective forecasts (i.e., lead > 0 wk). Forecast accuracy improvement for onset week, peak week, and peak intensity is summarized in Table 1. All improvements were found to be statistically significant. In particular, the improvements for onset week and peak week are more pronounced, up to 35% and 31%, respectively, several weeks before onset/peak. For peak intensity, which is a more challenging forecast target, the metapopulation forecasts still outperformed the isolated predictions, albeit to a lesser extent.

Our results indicate that the benefits of using a metapopulation system to predict onset are two-fold: (i) The number of state-season instances where an onset was predicted (as opposed to a forecast of no onset) was larger with the metapopulation system (SI Appendix, Table S2), and (ii) for instances where both the metapopulation and isolated systems predicted onset, the metapopulation system was more often correct.

To further validate the effectiveness of the proposed forecast framework, we applied it to a smaller geographical scale at the county level. In particular, seven contiguous counties in southeast Virginia with sufficient AFHSB ILI+ data were selected. Similar improvements were observed at the county level and are reported (SI Appendix, Retrospective Forecasts at the County Level and Fig. S14).

The above validation metrics are defined based on observations. However, these observations in real-world surveillance systems intrinsically contain a certain level of stochasticity. Another approach to validating the forecast system is to regard each observation as a single realization of a stochastic process. From this point of view, we can explore how the observations distribute with respect to the predicted values. If the forecasts are providing a good estimate of the truth, the observations, i.e., single realizations of a stochastic process, should be normally distributed around the predictions with zero mean, that is, an accurate forecast. If, on the other hand, the mean error is nonzero, there remains bias/inaccuracy in the prediction. By looking at the scatter of observations around many predictions, we can obtain an estimate of that inaccuracy. In Fig. 4, we show the distributions of the distance of the observations (onset week, peak week, and peak intensity) from predicted values across all 35 states and five seasons. In particular, we first grouped the forecasts according to their predicted lead to onset (peak) for onset week (peak week and peak intensity) predictions from 6 wk to 1 wk, and then plotted the distribution of the discrepancy of observations from corresponding predictions within each category, for both metapopulation and isolated predictions.

Unlike the metapopulation forecast system, the isolated forecast appears more sensitive to observation noise. For onset week, low incidence produces no-onset predictions due to a lack of signal, whereas, for the metapopulation model, early incidence signal is picked up from neighboring localities and inferred random movement levels. In addition, for the isolated forecasts predicting an onset, initial low-incidence weeks can drive the forecast system to erroneously predict a late onset, as indicated by the negative bias in onset predictions (Fig. 4, Top). For peak week, the isolated forecast tends to predict an outbreak peak too early in the future, and fails to capture the growing trend during low-incidence weeks (the positive bias observed in Fig. 4, Middle). For peak intensity, the isolated forecast trajectories tend to undershoot observed incidence (SI Appendix, Fig. S15B), which produces a positive bias in peak intensity prediction (Fig. 4, Bottom). These biases in the isolated forecasts are substantially reduced by the metapopulation forecast for all three targets.

The accuracy of onset week predictions can be further discriminated by the spread of the ensemble forecasts (13, 14). Metapopulation forecasts for the five seasons and 35 states were grouped into three categories according to their predicted onset lead (i.e., how far in the future onset is forecast to occur). Within each lead category, we plot forecast accuracy as a function of the ensemble variance of predicted onset log transformed, $\log ({\sigma }_{onset}^2)$, a measure of the within-ensemble spread of predictions (Fig. 5, Top). Forecasts with a smaller ensemble spread are found to be more reliable. This also holds for peak week and peak intensity, calibrated by the ensemble predicted peak variance log transformed, $\log ({\sigma }_{peak}^2)$ (Fig. 5, Middle and Bottom). This finding indicates that the likelihood a particular metapopulation model forecast is accurate can be estimated in real time.

Discussions

In this work, we have developed and validated an operational forecast system that can successfully predict the spatial transmission of influenza in the United States at both the state and county levels. By assimilating surveillance data from multiple locations, forecast accuracy for onset week, peak week, and peak intensity is enhanced by the metapopulation forecast system up to 35%, 31%, and 13%, respectively, compared with baseline forecasts made at each location in isolation. The proposed framework is flexible because its implementation is independent of specific disease transmission dynamics. As a result, it has potential to be adapted for use in the forecast of other infectious diseases.

The last few decades have witnessed the rapid emergence and worldwide spread of multiple novel infectious diseases, including severe acute respiratory syndrome, Middle East respiratory syndrome, Zika, Chikungunya, Ebola, H5N1 avian influenza, and 2009 pandemic H1N1 influenza. In an increasingly interconnected world, accurate forecast of the spatial progression of infectious disease can provide useful evidence-based information for policy-making and containment coordination. The forecast system presented here, which combines an intermediate-complexity metapopulation model, an efficient data assimilation technique, publicly available commuting and humidity data, and state-level surveillance data, can be run operationally in real time to generate predictions of future disease incidence and spread.

In the upcoming 2017–2018 season, the CDC will launch a real-time ILI forecast challenge for participating US teams (34). The proposed framework will be tested and evaluated in real time for this challenge. We will also incorporate the metapopulation framework into the real-time forecasts of ILI+ we generate and publish at state and municipal scales (35). These real-time efforts will facilitate the transformation of infectious disease forecasting from research to routine operation and its integration into decision-making.

Materials and Methods

Supplementary Material