To compile a dataset of laboratory-confirmed and possible human cases of HPAI H5N1 worldwide, we systematically searched four databases: PubMed, Scopus, Google Scholar, and the World Health Organization Global Alert and Response compilation (WHO GAR) [2]. We created an initial dataset using cases published in the WHO GAR website, which includes cases reported since November 2003. We assumed laboratory confirmation satisfying the WHO reporting standard for all such cases, although specific clinical laboratory data were not provided.

We then reviewed articles published in PubMed, Scopus, and Google Scholar using keywords “H5N1” and “human” or “humans.” We included articles in any language published between January 1, 1997 (the year of the initial Hong Kong outbreak) and April 19, 2013. We included articles describing confirmed or possible human cases of HPAI H5N1 virus infection. We defined confirmed cases as those meeting WHO reporting criteria: isolation of HPAI H5N1 virus, a positive result by reverse transcription polymerase chain reaction (RT-PCR) testing of clinical specimens using H5-specific primers and probes, an elevated H5-specific antibody titer of ≥1∶80 (or equivalent using the WHO protocol), or at least a fourfold rise in H5N1 virus neutralization antibody titer in paired sera [13]. We defined possible cases as those lacking laboratory confirmation but exhibiting symptoms and having known exposure to a confirmed human case. Due to limited availability of data on other modes of transmission (such as via consumption of raw poultry meat or exposure to feces-contaminated environments), this definition is more restrictive than the WHO definitions of suspected and probable cases. Two independent reviewers (RP, MM) assessed each article for these inclusion criteria; a third reviewer (NK) resolved disagreements.

For articles in Japanese, Russian, French, and Spanish, we verified inclusion with native speakers. A professional translator (YX) assessed inclusion criteria for the large number of Chinese articles. For all other languages, we converted articles to text format using PDF OCR X software (version 1.9.32, Burnaby, British Columbia), then used an online translation program to translate the article to English and verify inclusion [14].

For articles meeting inclusion criteria, we added each reported case to our dataset, abstracting a predefined set of variables (Abstraction Form S1). To avoid duplicate entries, we systematically compared demographic data on each case to all others in our dataset.

For each case, we extracted data on location (country and city) and date (month and year) of illness onset. We created a season variable based on month of case occurrence (defining summer as June–August, fall as September–November, winter as December–February, and spring as March–May). Based on country, we also defined a hemisphere variable (Northern or Equatorial; no Southern Hemisphere cases have been reported to date).

We defined a climate variable characterizing the predominant climate of the patient's city or region. We used the Köppen-Geiger climate classification map, which assigns one of 31 classifications incorporating three dimensions of climate: main climate (equatorial, arid, warm temperate, snow, polar), precipitation, and temperature [15]. The map is based on extensive data from the second half of the 20^th century, collected by the Global Precipitation Climatology Centre at the German Weather Service and the Climatic Research Unit of the University of East Anglia. For a small number of cases missing city or region data, but occurring in a country with a clearly predominant climate (n = 7), we imputed the country's predominant climate.

We performed all statistical analyses in R (Version 3.0.0, Vienna, Austria) and defined statistical significance by an alpha level of 0.05. Our primary analytic goal was to assess seasonal patterns in occurrence (frequency) of human cases of HPAI H5N1.

We first computed the proportion of cases occurring in each season. To assess significance of seasonal incidence patterns and to consider hemisphere as a moderator, we modeled case incidence using negative binomial regression. While the more common Poisson regression model is a canonical choice for modeling rare event occurrence, it relies on the assumption that the mean and variance of event frequency per unit time are equal. A common violation of this distributional assumption occurs when event frequency is over-dispersed: that is, the variance greatly exceeds the mean [16].

Negative binomial regression is a generalization of the Poisson model that allows a variance term to be fit to the data independent of the mean term, thus better modeling over-dispersion [16]–[18]. We note, however, that because Model 1 includes only orthogonal categorical factors, the maximum likelihood estimation process of any generalized linear model (GLM) in this case reduces mathematically to fitting separate means to each disjoint subset (defined by seasons). Thus, here, the choice of GLM is consequential only to the fitted standard errors and resulting p-values, not to the coefficients or predicted values.

We fit three nested models:(1)(2)(3)where Y is case frequency, b₀ is the intercept, x_s is season (summer, fall, winter, spring) and x_h is hemisphere (Northern, Equatorial). Thus, Model (1) allows us to assess the existence of seasonal patterns. We computed the incidence rate ratio (IRR) for fall, winter, and spring relative to summer by exponentiating their regression coefficients. Analogously to relative risk, the IRR represents the predicted multiplicative increase in event frequency for a covariate versus the reference level. That is, an IRR of 2 represents a twofold increase in event rate. As a secondary, post-hoc analysis, we used Akaike's Information Criterion (AIC) [19], a measure of fit-complexity tradeoff, and the likelihood-ratio test to compare Model (1) to Models (2) and (3), allowing us to assess the presence of a main effect of hemisphere and interaction of season with hemisphere. AIC is an entropy-based measure that decreases as the model becomes simpler (fewer parameters) and better-fitting (higher log-likelihood): thus, models with lower AIC are preferable.

Because cases have occurred in only 15 countries, with 55% of cases occurring in Egypt and Indonesia combined, inherent power limitations made quantitative assessment of climate type effects unfeasible. We instead present graphical representations as a preliminary analysis tool. We also manually synthesized the data we had compiled on HPAI H5N1 case occurrence with monthly weather data to visually assess the relationship between disease occurrence and local temperature, humidity, and season.

Our literature search identified 617 human cases of HPAI H5N1; we excluded from analysis six cases with missing data on season, yielding a final sample size of 611 (Figure 1). Cases occurred in 12 Northern Hemisphere countries and three equatorial countries (Table 1). Of note, all but two of the cases in equatorial countries were in Indonesia; the remaining two were in Nigeria and Djibouti.

We found n = 591 cases with data for both climate and season. Cases occurred in regions representing 11 climate types, with Arid/desert/hot arid (n = 170) and Equatorial/fully humid (n = 144) being the most frequent climate types. Fewer than 90 cases occurred in each of the other nine climate groups (Table 2, Figure 2).

Figure 3a displays case incidence by season and hemisphere. Across all countries, 76% of all cases occurred in winter and spring. In the winter and spring periods combined (n = 462), 3.1 times as many cases occurred as in the summer and fall months combined (n = 149). Among non-equatorial countries, 81% of cases occurred in winter and spring periods combined (n = 355), representing a 4.3-fold increase over summer and fall months combined (n = 83).

As expected, a classical Poisson model was overdispersed (residual deviance  = 1221.98, df = 60, p<0.001; also see Figure S1) [20] and yielded artificially low p-values. The negative binomial model appeared to resolve the poor fit due to overdispersion (residual deviance  = 65.40, df = 60; p = 0.29). Table 3 displays regression results for Models (1)–(3). Based on Model (1), case occurrence is higher in winter vs. summer (IRR = 4.45, p = 0.03) and marginally higher in spring vs. summer (IRR = 2.77, p = 0.13), but not in fall vs. summer (IRR = 1.33, p = 0.68).

Model (2), incorporating a main effect of hemisphere, and Model (3), incorporating an interaction of season with hemisphere, were not significantly better-fitting than model (1) (likelihood-ratio test: chi-square  = 2.16, df = 1, p = 0.14 and chi-square  = 2.82, df = 4, p = 0.59, respectively). Model (1) had comparable AIC (367.90) to Model (2) (367.74) and lower AIC than Model (3) (373.08). Thus, the season-only model adequately described seasonal variation; incorporating a main effect or interaction with hemisphere did not significantly improve fit. Model fit is summarized in Figure 4.

Figure 3b displays case incidence by season and climate. By visual inspection, the winter peak in case occurrence relative to summer and fall is robust across all six climate types. Case frequency in the spring appears more variable by climate group: three climate types (Arid/desert/hot summer, Equatorial/fully humid, and Equatorial/monsoonal) show some increase in frequency in spring relative to summer and fall, while the remaining three (Warm temperate/winter dry/hot summer, Equatorial/winter dry, and Other) do not show a spring peak in case occurrence.

Interactive Map S1 visually depicts the relationship between cases and temperature and humidity over the course of several years in Egypt and Indonesia, the two countries with the highest frequency of cases.