Original Research ARTICLE
A Space-Time Study of Hemorrhagic Fever with Renal Syndrome (HFRS) and Its Climatic Associations in Heilongjiang Province, China
- 1Institute of Islands and Coastal Ecosystems, Ocean College, Zhejiang University, Zhoushan, China
- 2Department of Geography, San Diego State University, San Diego, CA, United States
- 3Institute of Disease Control and Prevention, Academy of Military Medical Science, Beijing, China
Background: Hemorrhagic fever with renal syndrome (HFRS) is highly endemic in China, especially in Heilongjiang province (90% of all reported HFRS cases worldwide occur in China). The dynamic identification of high HFRS incidence spatiotemporal regions and the quantitative assessment of HFRS associations with climate change in Heilongjiang province can provide valuable guidance for HFRS monitoring, preventing and control. Yet, so far there exist very few and of limited scope quantitative studies of the spatiotemporal HFRS spread and its climatic associations in Heilongjiang province. Making up for this lack of quantitative studies is the reason for the development of the present work.
Method: To address this need, the well-known Bayesian maximum entropy (BME) method of space-time modeling and mapping together with its recently proposed variant, the projected BME (P-BME) method, were employed in this work to perform a composite space-time analysis and mapping of HFRS incidence in Heilongjiang province during the years 2005–2013. Also, using multivariate El Niño-Southern Oscillation index as a proxy, we proposed a combination of Hilbert-Huang transform and wavelet analysis to study the “HFRS incidence-climate change” associations.
Results: The main results of this work were two-fold: (1) three core areas were identified with high HFRS incidences that were spatially distributed and exhibited distinct biomodal temporal patterns in the eastern, western, and southern parts of Heilongjiang province; and (2) there exists a considerable association between HFRS incidence and climate change, particularly, an ~6 months period coherency was clearly detected.
Conclusions: The combination of modern space-time modeling and mapping techniques (P-BME theory, Hilbert-Huang spectrum analysis, and wavelet analysis) used in this work led to valuable quantitative findings concerning the spatiotemporal spread of HFRS incidence in Heilongjiang province and its association with climate change. Our essential findings include the identification of three core areas with high HFRS incidences in Heilongjiang province, and considerable evidence that HFRS incidence is closely related to climate change.
Hemorrhagic fever with renal syndrome (HFRS) is a rodent-borne zoonosis caused by Hantavirus (belonging to the Bunyaviridae family). In China, the Hantaan and Seoul viruses dominate HFRS infection, the leading rodent hosts of which are Apodemus agrarius and Rattus norvegicus, respectively [1, 2]. The virus is transmitted from rodents to humans via inhalation of aerosols contaminated by rodents' urine, saliva, excreta, and dung, possibly through ingestion of contaminated food and by direct contact of contaminated materials with broken skin or mucous membranes, or by rodent bites [3, 4]. Clinical HFRS manifestations include fever, headache, nausea, and abdominal pain. Complications, like adverse kidney effects and subsequent pulmonary edema, shock, renal insufficiency, encephalopathy, hemorrhages, and cardiac complications, can cause death [5, 6]. The disease goes through five stages: febrile, hypotensive shock, oliguric, polyuric, and convalescent, which last, respectively, 1–7 days, 1–3 days, 2–6 days, 2 weeks, and 3–6 months .
Historically, numerous HFRS-like cases have occurred in China going back to the tenth century AD. Currently, about 20,000–50,000 cases/year are reported in mainland China, which account for 90% of all reported cases worldwide [5, 8, 9]. The foci of HFRS often locates in rural areas, which constitute more than 70% of the total number of cases, because of poor housing conditions and abundant rodent hosts [10, 11]. According to the data of the National Health and Family Planning Commission of China, the HFRS death rate was 2.89% during 1950–2014 . Rattus norvegicus, which hosts the Seoul virus, is regarded as one of the most damaging invasive species around the world, and it is closely associated with humans, particularly in large metropolitan areas [12, 13]. As HFRS remains a severe public health problem, it is necessary to study historical HFRS evidence to provide rigorous scientific support to current disease monitoring and control procedures. Yet, so far there exists a very limited number of quantitative studies regarding the spatiotemporal HFRS distribution and spread in Heilongjiang province. In fact, the reason for the development of the present work is to make up the lack of such quantitative studies.
Based on province-level data, a number of nationwide studies have been carried out in China. Their results indicated that the geographical distribution of HFRS incidences was clustered, particularly in the northeastern, central and eastern parts of China. The observed hotspots shifted and expanded from year to year, whereas most HFRS cases were mainly reported during the spring and the autumn-winter seasons [14–16]. In these studies, global indicators of spatial autocorrelation (GISA), local indicators of spatial association (LISA) and Kulldorff's scan statistic were employed to characterize the spatial variation of HFRS incidence during several time periods. However, these studies suffered from certain drawbacks: (a) they considered neither the temporal nor the combined space-time (spatiotemporal) correlation of HFRS incidences, and (b) a fine temporal resolution (i.e., monthly data) was assumed, whereas the spatial resolution used for mapping purposes was rather coarse (i.e., at the province-level).
For more than two decades, the Bayesian Maximum Entropy theory of space-time data analysis and mapping (BME, [17–19]) has been proven to provide efficient and cost-effective methods for characterizing, predicting and mapping disease attributes (such as disease incidence) in a composite space-time domain under conditions of in-situ uncertainty . For example, Law et al.  used BME to qualitatively and quantitatively detect core areas with high syphilis incidence density in the city of Baltimore (USA) between the years 1994 and 2002; also, based on age-adjusted influenza mortality data at the county-level, Choi et al.  used BME to represent the space-time disease dependence structure of the disease, to map the influenza mortality rates and to assess the associated disease risk in the state of California (USA). In this work, we will use both the original BME and a BME variant (projected BME, P-BME) to study the spatiotemporal HFRS dependence pattern in Heilongjiang province, including the identification of particular disease features and the detection of high incidence areas.
Since HFRS is transmitted by reservoir hosts (especially rodents), it is expected that climatic factors (such as precipitation, temperature, humidity, and global climate pattern) should influence human HFRS morbidity by affecting the reproduction and abundance of rodents [3, 23–27]. To some extent, understanding climate change can offer a preliminarily assessment or an early warning concerning the epidemic situation. Few studies have investigated the intrinsic HFRS period and its inherent relationships with climate attributes and factors. In recent years, the Hilbert–Huang transformation (HHT, an adaptive method combining empirical mode decomposition and spectral analysis) has been developed for analyzing nonlinear and temporally non-stationary data, in general [28–30]. With this method, the intrinsic mode and intensity of a disease attribute is obtained that can provide useful insight regarding the temporal regulation of disease variation. Moreover, cross-wavelet transforms and wavelet coherence can explore the relationship between two series in the time-frequency domain . As a matter of fact, the wavelet method has been employed in the past to study disease incidence and climate change. For example, Thai et al.  have found a strong non-stationary association between El Niño-Southern Oscillation (ENSO) indices and climate variables and the corresponding dengue incidence in the Binh Thuan province (Vietnam) during a 2–3 years period; also, Chowell et al.  have suggested that in Peru the dengue incidence is significantly linked with the seasonal cycle timing of mean temperature variations.
In view of the above considerations, the objective of this work is two-fold: (1) to investigate the characteristics and spatiotemporal distribution of HFRS incidence in Heilongjiang province (China) using the BME and P-BME methods; and (2) to assess the intrinsic mode and coherence similarities between the HFRS incidence spread and the time series of key climatic factors using the Hilbert-Huang spectrum method and wavelet analysis.
Materials and Methods
Study Area and Data Collection
The Heilongjiang province is located in northeastern China, with an area of ~473,000 Km2 and a population of 38.35 million people. Remarkably, the Heilongjiang province is one of the highest HFRS morbidity regions in China . The Heilongjiang basin includes four major river systems: the Heilong, Songhua, Wusuli, and Suifen rivers (Figure 1). Monthly data of HFRS cases (21,383 cases in total) were collected at 130 counties and districts during the January 2005–December 2013 period by the China Information System for Disease Control and Prevention (CISDCP). Demographic data for each county were obtained from the National Bureau of Statistics of China. Subsequently, the HFRS cases were population-standardized and used in the present work.
The multivariate ENSO Index (MEI) integrates 6 oceanic and meteorologic variables over the tropical Pacific region to represent global climatic cycles (i.e., El Niño-Southern Oscillation): sea-level pressure, zonal and meridional components of surface wind, sea surface temperature, surface air temperature, and total cloudiness fraction of sky . The MEI has been used in the scientific literature to diagnose ENSO phenomena that can cause global climate variability, including world-wide correlations with temperature and precipitation data . In the present study, the MEI is employed as a proxy for representing climatic factors in order to explore their effects on HFRS incidence (MEI data is available at https://www.esrl.noaa.gov/psd/enso/mei/table.html).
Spatiotemporal Analysis and Mapping
Methodologically, the standardized HFRS incidence was regarded as a spatiotemporal random field (S/TRF, [19, 37]), denoted as X(p), with arguments p = (s,t) ∈ R2 × T in a composite space-time domain, where s denote the centroid coordinates of each administration unit and t ∈ T denotes the time argument. In this quantitative modeling setting, space (s) represents the order of co-existence and time (t) represents the order of successive existence of HFRS incidence distribution in Heilongjiang province. In-situ uncertainty manifests itself as an ensemble of possible HFRS realizations x regarding the space-time X(p) distribution, where the likelihood that each one of these possible realizations occurs is expressed by the corresponding HFRS probability density function f. In the BME method, two main knowledge bases are considered: the general or core knowledge base (G-KB), and the site-specific knowledge base (S-KB). The G-KB includes theoretical models of the space-time HFRS mean and covariance (correlation), whereas the S-KB consists of hard (exact) and soft (uncertain) HFRS data . The BME method uses the general knowledge available to generate the G–based (prior) probability density function, fG, of HFRS incidence distribution. Subsequently, the S-KB is incorporated to generate the combined G- and S-based probability density function 
where K = G+S denotes the total KB (core G and site-specific S), x = (xh, xs, xk) are HFRS realizations at the hard data points (ph), the soft data points (ps), and the unsampled (prediction) points (pk), and A is a normalization constant. After the probability density function fK has been derived at all prediction points pk, various HFRS incidence, X(p), estimates at these points are readily available, like the mean, the mode and the median X(p) values at each pk.
It is widely-recognized that the joint spatial-temporal covariance of a disease attribute is rather difficult to calculate experimentally (primarily due to the limited number of sample points) and to fit to a theoretical covariance model [39, 40]. Hence, it would be very useful to develop an improved method of space-time covariance fitting. Responding to this need, Christakos et al.  presented a BME variant, the projected BME (P-BME) method, which projects the disease incidence distribution (HFRS incidence in our case) from the original space-time disease domain R2 × T onto a lower dimensionality traveling space domain R2 by means of the simple set of equations
where denotes the spatial coordinates of HFRS incidence values, denotes the spatial lags (separation distances) between HFRS incidence values, and υ = (υ1, υ2) is the HFRS traveling vector. Then, are the space-time coordinates in the original R2 × T domain, which are matched one-to-one with the coordinates in the traveling R2 domain according to the P-BME method, i.e., the following transformation of HFRS incidence domains is used,
Accordingly, cX(h, τ) and are the HFRS incidence covariances in the space-time R2 × T and the traveling (spatial) R2 domains, respectively. The fact that these covariance functions are related by Equation (2b) provides a practical way to calculate υ in-situ. Significant advantages of the P-BME method is that, after projection, the covariance is located in the spatial (R2) domain, where (a) it is easier to fit a theoretical HFRS covariance model to the data, (b) more choices of theoretical HRFS covariance models are allowed, and (c) HSRF incidence mapping is more accurate and computationally efficient than in the space-time (R2 × T) domain.
By way of a summary, the P-BME method combines the BME Equation (1) with the traveling Equations (2a–c) to derive HFRS incidence predictions across space-time. The mean absolute HFRS incidence prediction error can be used to evaluate the accuracy of model cross-validation prediction. For comparison purposes, the 10-fold cross-validation method will be used below to test the performance of the direct BME method and the P-BME method in the R2 × T and the R2 domains, respectively. More theoretical and technical details regarding BME and P-BME can be found in the cited literature.
The HHT can help discover certain characteristics of the cumulative HFRS incidence in the entire Heilongjiang province, together with the MEIs and the underlying rules of their variation. Basically, HHT consists of two steps: empirical mode decomposition (EMD), and Hilbert transformation (HT).
More specifically, EMD is used to extract the intrinsic mode functions (IMFs), each of which is independent of the others, from the raw series
where cm represents the IMF component, and rn denotes residuals representing the raw series trend. We start by fitting the local minima and maxima values of the raw time series, Y(t), using cubic splines, and the mean time series, m1(t), is defined (i.e., the spline mean). We also define the standard deviation (SD)
where k = 1, 2, …. denotes the number of times the process is repeated (by convention, h10(t) = Y(t)). The first difference h11(t) = Y(t) − m1(t) is considered as the first IMF if it satisfies the criterion that the SD is between 0.2 and 0.3 , otherwise, the process is repeated using h11(t) as a raw series until the criterion are met. After the first IMF is extracted, the difference between the first IMF and raw series is used to identify subsequent IMFs. Then, a monotonic series is defined as the residual of the raw series trend mentioned earlier.
The HT can be applied for each IMF above to obtain the analytical function and its polar form
respectively, where H denotes the Hilbert transform, and ai(t) and θi(t) are the IMF amplitude and phase, respectively. Each IMF is the real part of its corresponding analytical function zi(t). In view of the IMF summation, see Equation (3), the raw series can be calculated by
where ωi(t) is the IMF instantaneous frequency. The Hilbert spectrum, H(ω, t), is subsequently defined as a frequency-time distribution with the amplitude of the raw series of Equation (6). And the marginal spectrum can be calculated by
Considering that EMD may cause mode mixing , in this work it was replaced by ensemble EMD (also denoted as EEMD). The process is as follows: add a white noise series into the raw series; use EMD to decompose the series with white noise into IMFs; repeat the two steps above (say, m times in total) with various white noise series; calculate the mean values of the corresponding m IMFs, which will be the final IMFs. The HHT can provide an amplified way of gaining insight into raw time series in the time-frequency domain, whereas the marginal spectrum offers an easy way to identify the cumulative amplitude distribution across different frequencies in a probabilistic sense . In this work, the entire standardized HFRS incidences in Heilongjiang province and the MEI, 108 months in total, were analyzed by HHT for comparison. The software code used for this purpose can be downloaded from http://rcada.ncu.edu.tw/research1_clip_program.htm.
The coherency between standardized HFRS incidence and MEI expresses the association between these two variables, in which case the coherent time-frequency areas can be easily discovered in terms of wavelet analysis [31, 44].
Specifically, coherency can measure the covariation intensity between two series. For this purpose, the time series, say Y1 and Y2, are transformed by the continuous wavelet equation
where ψ(t) is the mother wavelet, and * denotes the complex conjugate. Subsequently, the two transformed series are cross-wavelet transformed by
Finally, the wavelet transform coherency across time-frequency space can be calculated by
where α and τ are the scale factor and time shift, respectively, and S is a smoothing operator. As in the case of the BME method, more theoretical and technical details regarding HHT and wavelet analysis can be found in the cited scientific literature.
A total of 14,040 space-time records were included in the spatiotemporal analysis of this work. In the R2 × T domain, the HFRS incidence variability in Heilongjiang province was measured by means of the isotropic covariance plotted in Figure 2A. This space-time HFRS covariance combined two theoretical models, an exponential and a Gaussian (squared exponential) model, i.e.,
where 72 × 103 (m) and 2.6 (months) are, respectively, the spatial and temporal correlation ranges of the HFRS incidence distribution. The interpretation of the spatiotemporal covariance plots of Figure 2A implies that the distribution of HFRS cases during the period Jan 2005–Dec 2013 were controlled by spatial and temporal dependences. In quantitative terms, the covariance value [around 0.3 (cases/105)2] at the time lag τ = 4 (months) indicates a rather strong temporal dependence among the HFRS incidence values. On the other hand, the spatial neighborhood effect is about 200 (km), which also indicates a significant spatial dependence among incidence values.
Figure 2. Plots of the spatiotemporal empirical covariances and the fitted theoretical models in (A) the R2 × T domain, and (B) the R2 domain.
By inserting Equation (11) into Equation (2b), the traveling coefficient was calculated to be υ = |υ| = −10650.89τ. Then, the HFRS incidence data points were projected from the R2 × T domain onto the reduced dimensionality R2 domain (Figure S1). Following the projection process, the empirical HFRS covariance and the fitted covariance model in the R2 domain are plotted in Figure 2B, where the HFRS covariance model is analytically given by
This covariance is also of reduced dimensionality compared to that of Equation (11). In addition, the theoretical covariance model of Equation (12) provided a better fit to the corresponding empirical covariance than the covariance model of Equation (11).
HFRS Incidence Mapping
The HFRS incidence maps were firstly generated using the P-BME technique. For illustration, the monthly HFRS incidence maps for the year 2006 (January–December) are shown in Figure 3. Additional HFRS incidence maps for all the years considered can be found in the section of “Supporting Information.”
As regards spatial variation, three areas with considerable HFRS incidence are identified in the maps of Figure 3, particularly, in the eastern, the western and the southern parts of the Heilongjiang province. Among them, the eastern part shows high HFRS incidence values over a larger area than in the other two parts. As is noted in the section of Discussion, this happens because there exists a corresponding large area of croplands and rivers in the eastern part that are linked to increasing HFRS incidence. As regards the temporal HFRS variation in Heilongjiang province, the HFRS incidence variation exhibited two outbreaks within a year's time. Specifically, as is shown in Figure 3 and Figure S12 (Supporting Information), the HFRS incidence begins to increase in April of 2006, then a peak is reached in June of 2006, and the HFRS incidence reduces significantly in September of 2006. The next outbreak is observed during September 2006 to February 2007 with the peak occurring in November of 2006. Interestingly, the number of HFRS cases during the autumn-winter period is much larger than those during the spring-summer period. As noted in the section of Discussion, the interpretation of this phenomenon is that it is probably due to the fact that the autumn-winter period coincides with the rice harvest season, i.e., after harvest the soil condition switches from a flood state to a dry state, which leads to the rodents dispersal or migration causing a higher number of infected cases).
In the western part of the Heilongjiang province, low HFRS incidence values were observed during the months February–September of the period 2006–2013 (but not for the year 2005). The HFRS incidence at the southern part of the Heilongjiang province remained high during the months of June, October, November, and December of each year considered. Apparently, HFRS is transmitted to the southern part of the province from its eastern part. Overall, a declining trend of HFRS incidence is observed in the maps of Figures S3–S11 (Supporting Information) for the period 2005–2013, which may be due, at least in part to the development of medical condition and disease prevention.
A cross-validation analysis of P-BME mapping technique vs. direct BME mapping is plotted in Figure S2 during the years 2005–2013: specifically, the P-BME mapping was more accurate in predicting HFRS at low incidence points, whereas the direct BME mapping was more accurate in predicting HFRS incidence at high incidence points during 2005 and 2007 (see, also, the months 1–36 of the time series of HFRS incidence in Figure 4A). Overall, the P-BME was on average a better predictor of the space-time HFRS incidence distribution in the Heilongjiang province than the direct BME: the mean absolute prediction error for BME over the entire domain was 0.524 cases/100,000 individuals, whereas that of P-BME was 0.459 cases/100,000 individuals. An explanation for the above results is given in the section of Discussion below.
Figure 4. Ensemble EMD of (A) HFRS series and (B) MEI series. S is the original series and res. denotes the residual series.
Ensemble Empirical Mode Decomposition
Using the ensemble EMD method, five IMFs and one residual component were technically extracted from the HFRS and MEI series. The results are shown in Figure 4. In the first 40 months, the HFRS incidence series experienced much higher peaks than during the remaining months (Figure 4A), which is consistent with the findings of Figures S3–S11 (Supporting Information) above. The IMF frequency decreases from IMF1 to IMF5, whereas the corresponding IMF period increases. As it can be also seen, the long trends (residual component) of HFRS and MEI are decreasing.
Each IMF expresses a different fluctuation period (Table 1). For HFRS, the IMF1 to IMF5 represent incidence periods lasting 5.959, 8.936, 18.024, 39.002, and 64.004 months, respectively. The main HFRS inherent periods are 5.959, 8.936, and 18.024 months, according to the contribution percentage of IMF's variance. The corresponding MEI periods are 5.138, 17.662, 34.629, 50.129, and 108.466 months, respectively, whereas the main periods are 17.662, 34.629, and 50.129 months. By comparing the periods of each IMF, it was found that the IMF1 of HFRS has the same (6 months) period as the IMF1 of MEI. Similarly, the IMF3 of HFRS has the same (18 months) period with the IMF2 of MEI.
Hilbert–Huang and Marginal HFRS Incidence Spectra
The reason why the time-frequency technique (i.e., Hilbert–Huang transformation) is used at this study stage is that we seek to detect the similarities of the two series, in which case it can be concluded that the two series are inter-related. The Hilbert-Huang spectra of the HFRS and MEI series displayed in Figures 5A,B represent the time-frequency-energy distributions of the original series. As can be seen in Figures 5A,B, continuous instantaneous frequencies are detected in the low frequency region of the spectra (<0.15 cycle/month). The energy is also high in the low frequency region, especially for MEI, although some amounts of energy are detected in the middle and high frequency region of the spectra (>0.15 cycle/month).
Figure 5. Hilbert-Huang spectra of (A) the HFRS series and (B) the MEI series. The color bar ranging from dark blue to yellow indicates energy variation from minimum to maximum.
Moreover, the energy associated with HFRS incidence is more discretely distributed than that associated with MEI. As is shown in Figure 6, the marginal spectra of HFRS and MEI have a similar peak occurring at ~0.015–0.025 cycle/month. Also, we notice the peaks at 0.0694 and 0.0787 cycle/month. From 0 to 0.15 cycle/month, both spectra exhibit a decreasing trend. However, the MEI spectrum decreases rapidly to a small level and the energy is concentrated in this frequency range. Compared to the marginal MEI spectrum, the marginal HFRS spectrum presents a more complex fluctuation pattern across the entire frequency domain, and the amplitude (or energy) remains above the 0.15 cycle/month threshold.
Wavelet Coherency Analysis
As a result of the wavelet coherency analysis we obtained the coherency wavelet spectrum between the HFRS and MEI series shown in Figure 7A. This figure indicates that there exists a strong coherency between the two series with a periodicity of 6-months band, particularly, during the sampled months 22–31, 41–52, 75–80, and 85–94. A weaker coherency was also detected around the 12-month band during the sampled months 18~28, 42~52, and 61~83.
Figure 7. The HFRS-MEI association. (A) Wavelet coherency between HFRS and MEI series. The colors depict coherency values from black (0) to white (1); black line represents the cone of influence that delimits the region that is not influenced by edge effects; the dash dot line shows a = 5% significance level computed based on 500 simple bootstrap. (B) The oscillation series during the period 5–7 months band that was reconstructed from the real part of the wavelet transformed series. (C) The oscillation series during the period 8–16 months band reconstructed from the real part of the wavelet transformed series. Dash lines represent oscillation of the HFRS series, while black line represents oscillation of the MEI series.
The HFRS and MEI series reconstructed by means of the wavelet transformation are plotted in Figures 7B,C. These plots may be interpreted as providing an interesting demonstration of the HFRS and MEI series oscillations, which is useful for comparison purposes. A strong coherency is easily detected with 2 period month bands (i.e., 5–7 and 8–16 month band) during the sampled months mentioned above. This explains why the same oscillation pattern is found in Figures 7B,C.
Public health scientists and epidemiologists are increasingly in need of gaining insight about the space-time HFRS incidence variation, and about how climate change affects HFRS incidence dynamics. This is particularly true in Heilongjiang province, which is one of the most HFRS affected areas in China. To the best of our knowledge, very few studies have used analytical methods to describe the space-time HFRS spread in Heilongjiang province. On the other hand, the associations between HFRS incidence and climate factors have been always assessed in terms of numerical modeling, for example, autoregressive integrated moving average models (ARIMA), seasonal ARIMA (SARIMA), ecological niche models (ENM), Poisson regression models, multiple regression, conditional logistic regression, and principal components regression (PCR) models [25, 34, 45–48]. Interestingly, none of these studies explored the association between HFRS incidence and climatic factors in the context of their co-variation.
Responding to the above need, the present work is a collaborative effort between the Zhejiang University (Zhoushan, China), the Institute of Disease Control and Prevention (Beijing, China), and the San Diego State University (California, USA). This collaboration led to the introduction of a combination of modern space-time modeling and mapping techniques from BME theory, Hilbert–Huang spectrum analysis and wavelet analysis in the study of the spatiotemporal HFRS incidence distribution in Heilongjiang province, and its association with climate change.
In particular, one of the main elements of this study is the implementation of the P-BME method  to analyze the space-time HFRS incidence spread in Heilongjiang province. Monthly HFRS incidence data were analyzed and processed across the Heilongjiang province during the period 2005–2013. Monthly HFRS data are used here because they lead to more accurate predictions than annual data (see below) and, also, they can serve better our goals to detect the temporal HFRS incidence pattern (which can be explained in ecological terms) and assess the association between the HFRS incidence pattern and climate change (MEI).
A key feature of the P-BME method is that technically it transfers the study of HFRS incidence spread from the original 3-D (i.e., two space dimensions plus time, R2 × T) domain onto a reduced dimensionality 2-D (i.e., two space dimensions, R2) domain. In this way, the difficult to determine space-time distance (metric) is reduced to a much easier to define spatial distance, which means that the empirical space-time covariance of the monthly HFRS incidence distribution is accordingly transformed into a spatial covariance (see Figure 2). As a result, it is technically much easier to fit a theoretical model to the spatial than to the spatiotemporal empirical covariance of monthly HFRS data.
Next, for comparison purposes, the empirical covariances and the fitted theoretical models (in the R2 × T and the R2 domains) for the annual HFRS data are shown in Figure S13. This comparison shows that the temporal dependence of the annual HFRS data is much stronger than that of the monthly HFRS data (e.g., by comparing Figure 2 and Figure S13, it is seen that the annual HFRS covariance value at time lag τ = 4 is about 0.7 (cases/105)2 compared to 0.3 (cases/105)2 for the monthly covariance).
As regards the space-time mapping of HFRS spread in Heilongjiang province, it was found that the weaker the temporal dependence of the HFRS empirical covariance is, the better is the P-BME performance compared to that of the direct BME method. This improved performance of P-BME in this case is explained by the fact that since the time argument is technically imbedded within the projected coordinates of the reduced dimensionality domain, the temporal points with the stronger dependence (compared to the spatial ones) are not explicitly taken into account, instead, only traveling spatial points are considered in HFRS prediction and mapping, and the HFRS prediction error would be larger in this case. Compared to the original BME method, the P-BME was found to provide on average more accurate HFRS incidence predictions in the case of monthly HFRS incidence data (the mean absolute prediction error for P-BME is 0.459 cases/100,000 individuals vs. 0.524 cases/100,000 individuals for direct BME), whereas the opposite was the case for annual HFRS incidence data (the mean absolute prediction error for direct BME was 3.34 cases/100,000 individuals and for P-BME it was 4.56 cases/100,000 individuals).
Other findings of this study included the following. Three core areas were observed in the HFRS incidence distribution maps we obtained for the period Jan 2005–Dec 2013, particularly, the eastern, western and southern parts of Heilongjiang province. As the drainage map of Heilongjiang province shows (Figure S14), the Wusuli and Songhua rivers, as well as parts of the Heilong river belong to the eastern part, and the Neng and Mudan rivers belong to western and southern parts, respectively, of the Heilongjiang Province. Regarded as major water sources, these river basins provide a suitable environment for rodent hosts and their reproduction. The presence of a river or a pond can be a risk for human HFRS infection . Bao et al.  found that HFRS incidence had a strong correlation with distance to rivers (in particular, HFRS incidence had a quadratic relationship with distance to rivers, and R2 = 0.999, p = 0.000). On the other hand, severe droughts can significantly decrease HFRS incidence [2, 51]. As rivers can provide sufficient water for vegetation irrigation purposes, the croplands are always located near rivers, and Heilongjiang province is no exception.
In the Heilongjiang province, croplands, mixed forests, and cropland/natural vegetation mosaic account for 38.98, 26.29, and 16.92% of the territory, respectively (Figure S15). Specifically, croplands are largely distributed in the eastern and western parts of Heilongjiang province, and there are also some croplands in the southern part of Heilongjiang province, which correspond to the three core HFRS areas mentioned above. It has been found that crop production is highly correlated with HFRS incidence with a correlation coefficient r = 0.96 (p = 0.005) due to the fact that crop can directly or indirectly serve as food for rodent hosts . Increasing food availability contributes to the growth of rodent host population. This can raise the infection probability of humans, especially farmers, who have a higher likelihood to come in contact with these animals [52, 53]. Farmers usually don't have steady jobs other than farming, and they also have higher mobility compared to other professions. Therefore, they may carry and spread hantaviruses to wider areas than the croplands during their traveling after the farming season. Under these conditions, HFRS may spread rapidly. Being aware of the above high HFRS spread likelihood, it is necessary to implement public health interventions in the core areas of the Heilongjiang province to avoid HFRS outbreaks and spread (these interventions include, e.g., the extermination of potential rodent hosts and vaccination).
In addition, June and November HFRS incidence peaks were found during the years 2005–2013 (Figure S12). Interestingly, this bimodal temporal pattern was also found in Hubei province after 1995, as result of the Seoul virus-related HFRS spring outbreaks and the Hantaan virus-related HFRS winter outbreaks . Another study found that HFRS associated with wild and house rodents occur during different seasons . In view of the fact that land-use can affect virus occurrence in hosts by influencing movement and contact rate [56, 57], we notice that Heilongjiang is located in a high latitude area, where paddy rice only grows once a year (during May and October). At the beginning of sowing season in May, soils are irrigated to remain in a flood state for paddy rice growth, and the soil conditions change from drying to flooding. After harvesting, the soil will return to drying conditions in October. These switches between different land environmental conditions may cause a cascade of factors contributing to infectious disease emergence, especially invasive alien species due to the fact that they are invaders over their natural range in the new environment [58–60]. More specifically, an agricultural ecosystem is particularly vulnerable to invasive alien species and anthropogenic activities that can initiate or accelerate the introduction or invasion of alien species . Sufficient food availability during June contributes to rodent reproduction, whereas insufficient food availability or the dry conditions of November will result in spur sudden dispersal or migration events, both of which can increase HFRS infection . Whats more, the edges of paddy fields may involve ecotones as habitats with infectious disease and animal reservoir hosts being abundant in wildlife . Following a month-long period of environmental change, HFRS outbreaks may occur.
Generally, the study of real world phenomena is highly complex and interdisciplinary, in particular the study of large spatial scale climate and disease variation (which can be affected by biological, social, geographic, economic, medical factors etc.; ). Hence, the observation series may contain a large amount of direct and indirect information. With such complex information, it can be really hard for scientists to collect data from various disciplines and explore the relationship among them through hypothesis- and equation-driven methods.
In view of the above considerations, time-frequency analysis methods, regarded as data-mining tools, constitutes another major component of the present study. The results can reveal the intrinsic variation patterns of HFRS incidence and MEI series, as well as the dynamic characteristics of the HFRS and MEI cycles in the time-frequency domain. Understanding the association between HFRS incidence and climate change (using MEI as a proxy, measuring coupled oceanic-atmospheric character of ENSO event) provides a potential auxiliary way to assess the public health effects of global climate change, since climate variability has important effects on wildlife population dynamics [64, 65]. The Hilbert-Huang transformation is a powerful tool for solving mode-mixing problems and can be also used as a filter for decomposing raw HFRS incidence series into several independent series with disparate modes, i.e., IMFs . Different component series were obtained that describe various inherent disease characteristics that cannot be detected in the raw series. Our results showed that both HFRS and MEI series have six types of characteristic components. A monotonic declining trend is shown in Figure 4, and both series are characterized by 6- and 18-months periods, approximately (Table 1), indicating that similar patterns are hidden in the variation features of the two series. For further analysis, the Hilbert–Huang and the marginal spectra were used to assess the strength of series variation in the combined time-frequency domain. For both series, stable and consistent variations were observed at low frequency regions, although certain discrete fluctuations can be found in the HFRS spectra that are not observed in the MEI spectra. Such differences may not be explained by climate change (MEI) but rather in terms of non-climatic factors, e.g., population immunity, public health condition, and socio-economic factors [10, 54, 67].
Moreover, wavelet coherency analysis showed that the HFRS incidence has a strong association with MEI during a 6-month period (Figure 7). These results suggest that the HFRS incidence dynamics are interrelated with climate change and the MEI can serve as a potential predictor of HFRS occurrence. We notice that similar results have been found for diseases like dengue fever, dengue hemorrhagic fever, hantavirus cardiopulmonary syndrome, and malaria [68–71]. Moreover, regional precipitation is known to be influenced by ENSO, showing the strongest interrelation with climate variability around the Globe . Increasing precipitation provides sufficient soil moisture for improving ecosystem productivity . As a result, the number of rodents grows rapidly, leading to increasing contact rates between rodents and between rodents and humans . The HFRS infection rate increases under the above ecological changes. A deeper understanding of the association between climate change and HFRS incidence can provide a potential tool of early HFRS outbreak warning, especially concerning short-term effects.
Certain limitations of the present work should be acknowledged, which are rather typical for this kind of quantitative studies. The first one is data limitation, i.e., the HFRS dataset used is an aggregated set that does not distinguish the infectious HFRS types, e.g., Hantaan virus or Seoul virus, for which the infectious dynamics may be different. Second, the impact of climate change on human HFRS incidence can be twofold: impact from rodents to rodents and infection from rodents to humans. Therefore, distinct studies of these two possibilities would offer a better understanding of HFRS transmission. Third, some other impact factors could also be included in HFRS pattern analysis, such as population immunity and socio-economic factors. Future work should focus on a more detailed analysis of spatiotemporal intensity differences of the various environmental factors impacting HFRS incidence.
In sum, the present work provides a quantitative study of the HFRS incidence spread in Heilongjiang province (China). Three core areas with high HFRS incidences were identified. In addition, time-frequency analysis provides evidence that HFRS incidence is closely related to climate change.
All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The authors acknowledge with appreciation the Matlab codes (processing wavelet analysis) provided by Professor Bernard Cazelles. The work was supported by grants from the National Natural Science Foundation of China (No. 41671399 and No. 11501339).
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/article/10.3389/fams.2017.00016/full#supplementary-material
1. Yan L, Fang LQ, Huang HG, Zhang LQ, Feng D, Zhao WJ, et al. Landscape elements and Hantaan virus–related hemorrhagic fever with renal syndrome, People's Republic of China. Emerg Infect Dis. (2007) 13:1301–6. doi: 10.3201/eid1309.061481
2. Tian H, Yu P, Bjørnstad ON, Cazelles B, Yang J, Tan H, et al. Anthropogenically driven environmental changes shift the ecological dynamics of hemorrhagic fever with renal syndrome. PLoS Pathog. (2017) 13:e1006198. doi: 10.1371/journal.ppat.1006198
3. Liu J, Xue FZ, Wang JZ, Liu QY. Association of haemorrhagic fever with renal syndrome and weather factors in Junan County, China: a case-crossover study. Epidemiol Infect. (2013) 141:697–705. doi: 10.1017/S0950268812001434
4. Chinikar S, Javadi A, Hajiannia A, Ataei B, Jalali T, Khakifirouz S, et al. First evidence of Hantavirus in central Iran as an emerging viral disease. Adv Infect Dis. (2014) 4:173–7. doi: 10.4236/aid.2014.44024
5. Simmons JH, Riley LK. Hantaviruses: an overview. Compar Med. (2002) 52:97–110. Available online at: http://www.ingentaconnect.com/content/aalas/cm/2002/00000052/00000002/art00003
14. Li S, Ren H, Hu W, Lu L, Xu X, Zhuang D, et al. Spatiotemporal heterogeneity analysis of hemorrhagic Fever with renal syndrome in china using geographically weighted regression models. Int J Environ Res Public Health (2014) 11:12129–47. doi: 10.3390/ijerph111212129
15. Zhang S, Wang S, Yin W, Liang M, Li J, Zhang Q, et al. Epidemic characteristics of hemorrhagic fever with renal syndrome in China, 2006–2012. BMC Infect Dis. (2014) 14:384. doi: 10.1186/1471-2334-14-384
16. Zhang WY, Wang LY, Liu YX, Yin WW, Hu WB, Magalhaes RJ, et al. Spatiotemporal transmission dynamics of hemorrhagic fever with renal syndrome in China, 2005–2012. PLoS Negl Trop Dis. (2014) 8:e3344. doi: 10.1371/journal.pntd.0003344
21. Law DCG, Bernstein KT, Serre ML, Schumacher CM, Leone PA, Rompalo AM, et al. Modeling a syphilis outbreak through space and time using the Bayesian maximum entropy approach. Ann Epidemiol. (2006) 16:797–804. doi: 10.1016/j.annepidem.2006.05.003
22. Choi KM, Yu HL, Wilson ML. Spatiotemporal statistical analysis of influenza mortality risk in the State of California during the period 1997–2001. Stoch Environ Res Risk Assess. (2008) 22:15–25. doi: 10.1007/s00477-007-0168-4
23. Clement J, Vercauteren J, Verstraeten WW, Ducoffre G, Barrios JM, Van Ranst A-M, et al. Relating increasing hantavirus incidences to the changing climate: the mast connection. Int J Health Geogr. (2009) 8:1. doi: 10.1186/1476-072X-8-1
24. Zhang WY, Guo WD, Fang LQ, Li CP, Bi P, Glass GE, et al. Climate variability and hemorrhagic fever with renal syndrome transmission in Northeastern China. Environ Health Perspect. (2010) 118:915–20. doi: 10.1289/ehp.0901504
25. Liu X, Jiang B, Gu W, Liu Q. Temporal trend and climate factors of hemorrhagic fever with renal syndrome epidemic in Shenyang City, China. BMC Infect Dis. (2011) 11:331. doi: 10.1186/1471-2334-11-331
26. Xiao H, Tian HY, Cazelles B, Li XJ, Tong SL, Gao LD, et al. Atmospheric moisture variability and transmission of hemorrhagic fever with renal syndrome in Changsha City, Mainland China, 1991–2010. PLoS Negl Trop Dis. (2013) 7:e2260. doi: 10.1371/journal.pntd.0002260
27. Tian HY, Yu PB, Luis AD, Bi P, Cazelles B, Laine M, et al. Changes in rodent abundance and weather conditions potentially drive hemorrhagic fever with renal syndrome outbreaks in Xi'an, China, 2005–2012. PLoS Negl Trop Dis. (2015) 9:e0003530. doi: 10.1371/journal.pntd.0003530
28. Huang NE, Shen Z, Long SR, Wu MC, Shih HH, Zheng Q, et al. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc R Soc A (1998) 454:903–95. doi: 10.1098/rspa.1998.0193
29. Huang NE, Wu ML, Qu W, Long SR, Shen SS. Applications of Hilbert–Huang transform to non-stationary financial time series analysis. Appl Stoch Models Bus Indus. (2003) 19:245–68. doi: 10.1002/asmb.501
31. Grinsted A, Moore JC, Jevrejeva S. Application of the cross wavelet transform and wavelet coherence to geophysical time series. Nonlin Process Geophys. (2004) 11:561–6. doi: 10.5194/npg-11-561-2004
32. Thai KT, Cazelles B, Van Nguyen N, Vo LT, Boni MF, Farrar J, et al. Dengue dynamics in Binh Thuan province, southern Vietnam: periodicity, synchronicity and climate variability. PLoS Negl Trop Dis. (2010) 4:e747. doi: 10.1371/journal.pntd.0000747
33. Chowell G, Cazelles B, Broutin H, Munayco CV. The influence of geographic and climate factors on the timing of dengue epidemics in Perú, 1994-2008. BMC Infect Dis. (2011) 11:164. doi: 10.1186/1471-2334-11-164
34. Li CP, Cui Z, Li SL, Magalhaes RJS, Wang BL, Zhang C, et al. Association between hemorrhagic fever with renal syndrome epidemic and climate factors in Heilongjiang Province, China. Am J Trop Med Hyg. (2013) 89:1006–12. doi: 10.4269/ajtmh.12-0473
35. Curran LM, Caniago I, Paoli GD, Astianti D, Kusneti, M, Haeruman H. Impact of El niño and logging on canopy tree recruitment in borneo. Science (1999) 286:2184–8. doi: 10.1126/science.286.5447.2184
43. Ju K, Zhou D, Zhou P, Wu J. Macroeconomic effects of oil price shocks in China: an empirical study based on Hilbert–Huang transform and event study. Appl Energy (2014) 136:1053–66. doi: 10.1016/j.apenergy.2014.08.037
45. Bi P, Wu X, Zhang F, Parton KA, Tong S. Seasonal rainfall variability, the incidence of hemorrhagic fever with renal syndrome, and prediction of the disease in low-lying areas of China. Am J Epidemiol. (1998) 148:276–81. doi: 10.1093/oxfordjournals.aje.a009636
46. Xiao H, Gao L, Li X, Lin X, Dai X, Zhu P, et al. Environmental variability and the transmission of haemorrhagic fever with renal syndrome in Changsha, People's Republic of China. Epidemiol Infect. (2013) 141:1867–75. doi: 10.1017/S0950268812002555
48. Xiao H, Liu HN, Gao LD, Huang CR, Li Z, Lin XL, et al. Investigating the effects of food available and climatic variables on the animal host density of hemorrhagic fever with renal syndrome in Changsha, China. PLoS ONE (2013) 8:e61536. doi: 10.1371/journal.pone.0061536
49. Li J, Chen Z, Hou T, Xing Y, Cai Z, Yang Y. Study on the risk factors of hemorrhagic fever with renal syndrome in Xi'an City (in Chinese). Chin J Dis Control Prevent. (2013) 17:564–6. Available online at: http://www.cqvip.com/qk/80677a/201307/46406708.html
50. Bao C, Liu W, Zhu Y, Liu W, Hu J, Liang Q, et al. The spatial analysis on hemorrhagic fever with renal syndrome in Jiangsu Province, China based on geographic information system. PLoS ONE (2014) 9:e83848. doi: 10.1371/journal.pone.0083848
51. Yang H, Xia X, Wu Y. Analysis of trend and impact factors of hemorrhagic fever with renal syndrome in Xi'an City, 1999-2000 (in Chinese). In: First Annual Academic Conference of Chinese Preventive Medicine Association (Jinan) (2002), p. 330–1. Available online at: http://cpfd.cnki.com.cn/Article/CPFDTOTAL-ZHYF200208001678.htm
52. Ruo SL, Li YL, Tong Z, Ma QR, Liu ZL, Fisher-Hoch SP. Retrospective and prospective studies of hemorrhagic fever with renal syndrome in rural China. J Infect Dis. (1994) 170:527–34. doi: 10.1093/infdis/170.3.527
53. Engelthaler DM, Mosley DG, Cheek JE, Levy CE, Komatsu KK, Frampton JW. Climatic and environmental patterns associated with hantavirus pulmonary syndrome, Four Corners region, United States. Emerg Infect Dis. (1999) 5:87. doi: 10.3201/eid0501.990110
54. Zhang YH, Ge L, Liu L, Huo XX, Xiong HR, Liu YY, et al. The epidemic characteristics and changing trend of hemorrhagic fever with renal syndrome in Hubei Province, China. PLoS ONE (2014) 9:e92700. doi: 10.1371/journal.pone.0092700
57. Zhang W-Y, Fang L-Q, Jiang J-F, Hui F-M, Glass GE, Yan L, et al. Predicting the risk of hantavirus infection in Beijing, People's Republic of China. Am J Trop Med Hyg. (2009) 80:678–83. doi: 10.4269/ajtmh.2009.80.678
59. Patz JA, Daszak P, Tabor GM, Aguirre AA, Pearl M, Epstein J, et al. Unhealthy landscapes: policy recommendations on land use change and infectious disease emergence. Environ Health Perspect. (2004) 112:1092–8. doi: 10.1289/ehp.6877
60. Tian H, Huang S, Zhou S, Bi P, Yang Z, Li X, et al. Surface water areas significantly impacted 2014 dengue outbreaks in Guangzhou, China. Environ Res. (2016) 150:299–305. doi: 10.1016/j.envres.2016.05.039
61. Parmenter RR, Yadav EP, Parmenter CA, Ettestad P, Gage KL. Incidence of plague associated with increased winter-spring precipitation in New Mexico. Am J Trop Med Hyg. (1999) 61:814–21. doi: 10.4269/ajtmh.1999.61.814
66. Tsai FF, Fan SZ, Lin YS, Huang NE, Yeh JR. Investigating power density and the degree of nonlinearity in intrinsic components of anesthesia EEG by the hilbert-huang transform: an example using ketamine and alfentanil. PLoS ONE (2016) 11:e0168108. doi: 10.1371/journal.pone.0168108
69. Hjelle B, Glass GE. Outbreak of hantavirus infection in the four corners region of the United States in the wake of the 1997–1998 El Nino—Southern oscillation. J Infect Dis. (2000) 181:1569–73. doi: 10.1086/315467
70. Zhou G, Minakawa N, Githeko AK, Yan G. Association between climate variability and malaria epidemics in the East African highlands. Proc Natl Acad Sci USA. (2004) 101:2375–80. doi: 10.1073/pnas.0308714100
71. Tipayamongkholgul M, Fang CT, Klinchan S, Liu CM, King CC. Effects of the El Niño-Southern Oscillation on dengue epidemics in Thailand, 1996-2005. BMC Public Health (2009) 9:422. doi: 10.1186/1471-2458-9-422
Keywords: hemorrhagic fever with renal syndrome, spatiotemporal, mapping, Bayesian maximum entropy, Hilbert-Huang transformation, wavelet analysis
Citation: He J, Christakos G, Zhang W and Wang Y (2017) A Space-Time Study of Hemorrhagic Fever with Renal Syndrome (HFRS) and Its Climatic Associations in Heilongjiang Province, China. Front. Appl. Math. Stat. 3:16. doi: 10.3389/fams.2017.00016
Received: 27 March 2017; Accepted: 18 July 2017;
Published: 07 August 2017.
Edited by:Aristides (Aris) Moustakas, Queen Mary University of London, United Kingdom
Reviewed by:Petros Damos, Aristotle University of Thessaloniki, Greece
Alex Hansen, Norwegian University of Science and Technology, Norway
Copyright © 2017 He, Christakos, Zhang and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.