To further assess the hypothesis that AH is a fundamental driver of influenza seasonality, we examined whether a population-level model of influenza transmission forced by AH conditions could reproduce the observed seasonal patterns of P&I mortality. We simulated influenza transmission for 5 states representative of different climates within the US: Arizona, Florida, Illinois, New York and Washington. The model considers 3 disease classes: susceptible, infected, recovered; to integrate the impact of waning immunity following antigenic drift, we allow individuals to go back to the susceptible class at a defined rate (SIRS model). Observed 1972-2002 daily AH conditions within each state are used to modulate the basic reproductive number,
, of the influenza virus, i.e. the per generation transmission rate in a fully susceptible population. These daily fluctuations of
alter the transmission probability per contact within the SIRS model and thus affect influenza transmission dynamics. The SIRS model contains 4 free parameters: 2 (
) that define the range of
, 1 for the duration of immunity (D
), and 1 for the duration of infectiousness (L
If absolute humidity controls influenza seasonality, best-fit simulations with the AH-driven transmission model should meet the following criteria: 1) the mean annual model cycle of infection should match observations in each state; 2) these simulations should converge to similar parameter values, i.e. the virus response to AH should be consistent among states; and 3) AH modulation of transmission rates (
) within the model must match the large range implied by the laboratory data (Figure 1).
Multiple 31-year (1972-2002) simulations were run at each of the 5 states with randomly chosen parameter combinations. We then compared the mean annual cycle of daily infection from each simulation with a similar average of 1972-2002 observed excess P&I mortality rates ,
. Best-fit model simulations at each site capture the observed seasonal cycle of influenza (Figure 3). These simulations not only produce the late-year rise in transmission and infection, but also the wintertime peak during early January, typically followed by a secondary peak during late February/early March. In both models and observations, the dual winter peaks are not typically seen in individual years; rather these epidemic trajectories reflect the averaging of individual wintertime outbreaks that peak anytime between December and April (Figure S5).
Figure 3. Mean annual cycles for the best-fit SIRS model simulations at the 5 state sites. Here, best-fit simulations were selected individually for each state based on RMS error after scaling the 31-year mean daily infection number to the 31-year mean observed daily excess P&I mortality rate. Thick blue line shows the best-fit simulation; thinner green lines show the next 9 best simulations.
We also searched for the best-fit parameter combinations for all 5 sites evaluated together. The parameter combinations of these best ‘combined fits’ are characterized by high
(generally >2.8), high
(>1) and low mean infectious period (2<D
<4.2 days) (Table 2, Figure S6). Best-fit simulations at each of the 5 sites individually occupy a similar parameter space (Table S2, Figures S7-S11). In particular, these simulations converge to high
, which indicates a similar response to AH variability (see Supporting Information).
Parameter combinations for the 10 best-fit simulations at the Arizona, Florida, Illinois New York, and Washington state sites. 5000 simulations were performed at each site with the parameters
randomly chosen from within specified ranges. Best-fit simulations were selected for the 5 sites in aggregate based on RMS error after scaling the 31-year mean daily infection number to the 31-year mean observed daily excess P&I mortality rate at each site. The scaling factor itself, representing mortality per infection is also shown.
There is some correlation among SIRS model parameter values in simulations that fit the observed excess P&I mortality well. For instance, among better-fit simulations, L
tend to be inversely related (Figures S6-S11). In addition, broad regions of parameter space appear capable of producing high-quality, low-RMS-error simulations (Figure S6). The stochastic components of the SIRS model may contribute in part to this behavior. The flat goodness-of-fit within model parameter space indicates that no one parameter combination is strictly ‘best’; rather, a range of parameter combinations may produce good simulations of influenza transmission. These parameter ranges are: L
= 3-8 years, D
= 2-3.75 days,
= 2.6-4, and
= 1.05-1.30. We re-ran the SIRS model repeatedly sampling this approximate subset range of parameter space. Best-fit simulations from this subset range of parameter space (Table S3) were of similar quality and exhibited the same flat goodness-of-fit within model parameter space as the best-fit simulations presented in Table 2.
Because the SIRS model simulates only influenza-related infections, not deaths, a scaling factor is needed to compare model-simulated rates of infection with the observed excess P&I mortality rates. This scaling factor can be understood as the case fatality ratio, i.e. the probability of mortality given infection. Reassuringly, all best-fit simulations produce a scaling factor of the same order of magnitude and roughly consistent with the expected value of the case fatality ratio for P&I-related deaths (see Supporting Information).
The model also explains regional variations in influenza dynamics. Due to the modeled nonlinear relationship between
and AH (Figure 1c), the seasonal cycle of
is sensitive to both AH seasonal cycle amplitude and mean AH levels (Figures 1d-e). In Florida, mean AH levels are higher than for the other 4 states, but the seasonal AH cycle remains large and produces a seasonal
cycle of sufficient amplitude to generate an effective reproductive number,
, greater than 1 (Figure 1f) and organize influenza epidemics preferentially during winter. Outbreak dynamics reinforce this phase organization in that wintertime epidemics confer immunity to a large proportion of the model population, which then reduces population-level susceptibility during the following summer when
is low. In Arizona and Washington state, the seasonal AH cycle is less than for the other 3 states, but average AH levels are low, where laboratory findings indicate sensitivity to variation in AH is greater; consequently
retains a sizeable seasonal cycle (Figure 1e). For all 5 states the AH-driven seasonal variation of
is large enough that
is strongly modulated by AH conditions and exceeds 1 during winter as outbreaks develop (Figure 1f).
The humidity-driven SIRS simulations satisfy our 3 criteria for supporting the hypothesis that AH controls influenza seasonality in temperate regions. The simulations produce a consistent response in the 5 climatologically diverse US states using similar parameter values. The large sensitivity of simulated influenza transmission to AH is consistent with the analysis of laboratory experiments that show large changes in influenza virus survival and transmission in response to AH variability (Figures 1a,b).
Cross Validation of the Model Findings
To further validate the SIRS model findings, we determined whether the best-fit simulations derived from the 5 selected states could reproduce the seasonal cycles of influenza elsewhere in the US. The 10 best combined-fit parameter combinations (Table 2) were used to perform 31-year (1972-2002) SIRS simulations at each of the contiguous 48 states plus DC.
The results of this cross-validation demonstrate good simulations of observed excess P&I mortality for a majority of states (average r
> 0.7, minimum r
> 0.5, see Methods and Table 3). Some states, particularly the sparsely populated western states perform less well. These states often have low workflow 
, which may reduce the rate of introduction of the virus each winter. In addition, heterogeneous AH fields across some states (particularly large ones) create some error due to simulation with a single average statewide AH value. 13 states in the continental US, including Arizona, possess low workflow rates 
. 6 of these 13 states are among the 10 worst cross-validation performers; such a clustering is unlikely to occur by chance alone (p
< 0.005). In addition, 7 of the 10 worst performers are states with the 10 lowest population densities (p
Table 3. Correlation coefficients for the contiguous US and District of Columbia of SIRS-simulated 1972-2002 influenza incidence with 1972-2002 observed excess P&I mortality. The 10 best common-fit parameter combinations (Table 1) were used for these hindcast projections. Results are ordered based on best average correlation (among the 10 simulations for each state). States marked with an asterisk are low workflow states. The 10 states with lowest 1972-2002 population density are shown in bold.
Overall, the cross validation shows that the best combined-fit parameter combinations can simulate influenza seasonality throughout the country. Future use of higher resolution AH and observed P&I data that better represent local conditions may improve these model results.
Additional SIRS Model Results
We also used SIRS model simulations to provide additional support for the association between negative
and epidemic onset (Figure 2). Best-fit SIRS model runs reveal a comparable effect in which large negative
develop about 2 weeks prior to onset as defined by SIRS model infection rates (Figure 4). The 1-week difference in lag between this analysis with model infection rates (2 weeks) and the analysis with observed excess P&I mortality rates (3 weeks) roughly corresponds to the mediantime from infection to mortality,,
. The broader peak of negative
seen in Figure 2 is likely due to other, real-world factors that affect onset response and are not represented in the SIRS model (see Supporting Information).
Plots of average
associated with wintertime influenza onset for 10 best-fit SIRS model simulations at the 5 state sites (Arizona, Florida, Illinois, New York and Washington). The onset dates are defined as the date on which wintertime infection rates have been at or above a prescribed level for two continuous weeks (e.g. 50 infections/100,000 people). Each solid line is the averaged
associated with influenza onset as defined by a different threshold infection rate. The dashed line shows
. bottom) Plots of
anomalies using the
anomalies are calculated using the parameters
from each best-fit simulation (Table S3). The dashed line shows
Finally, we examined whether the school calendar, which alters person-to-person contact rates, could provide a better simulation of seasonal influenza than AH. School holidays have been estimated to lead to changes of ~25% in influenza transmission 
and occur during summer in the US, as well as at the end of the calendar year and again in spring. A number of SIRS model simulations were performed that included a step-wise increase of
during the school year (see Supporting Information). Simulations in which school closure was the only modulation of
were able to reproduce the seasonal cycle of influenza; however, these simulations did not reproduce observed excess P&I mortality as well as those with AH alone (see Supporting Information, Table S5, Figures S14-S15). In addition, a 40-90% change in influenza transmission (
) was needed to effect this seasonality (Table S5). This range of
changes is slightly larger than the previously estimated modulation of ~25%; however, these previous estimates were derived from an age-structured population model, so direct comparison is difficult.