研究登革热的数学方法

rf_001

生物谷报道：Nature最新报道采用数学方法研究疾病的流行病学规律，这一研究成果将大大开创流行病学研究的新篇章。该研究成为Nature的封面文章，可见其价值。登革热是作为本期封面文章的一篇结构论文所研究的对象，每年受其影响的人口多达1亿。登革热最严重的传染形式为登革出血热，会对东南亚国家的医疗卫生体系造成沉重负担。一组流行病学家、数学家和物理学家，用一种最初由NASA开发用来分析物理材料中的波的数学方法，对登革热病例的“流行波”进行了研究。分析结果表明，泰国的登革热疫情在传播速度上表现出一个可预测的三年振荡，疫情发源于曼谷，然后向全国各地传播，传播速度可达每月148公里，相当于疾病的一个“行波”。因此，对大城市进行有效监测，有可能对疫情进行早期预警，使人们有时间采取措施消灭传播该疾病的蚊子，有时间正确分配有限的公共卫生资源。

Nature 427, 344 - 347 (22 January 2004); doi:10.1038/nature02225  


Travelling waves in the occurrence of dengue haemorrhagic fever in Thailand


DEREK A.T. CUMMINGS1,2, RAFAEL A. IRIZARRY3, NORDEN E. HUANG4, TIMOTHY P. ENDY5, ANANDA NISALAK6, KUMNUAN UNGCHUSAK7 & DONALD S. BURKE2

1 Department of Geography and Environmental Engineering, Johns Hopkins University, Baltimore, Maryland 21218, USA
2 Department of International Health, Johns Hopkins Bloomberg School of Public Health, Baltimore, Maryland 21205, USA
3 Department of Biostatistics, Johns Hopkins Bloomberg School of Public Health, Baltimore, Maryland 21205, USA
4 Laboratory for Hydrospheric Processes/Oceans and Ice Branch, NASA Goddard Space Flight Center, Greenbelt, Maryland 20771, USA
5 Virology Division, United States Army Medical Research Institute in Infectious Disease, Fort Detrick, Maryland 21702, USA
6 Department of Virology, Armed Forces Research Institute of Medical Sciences, Bangkok 10400, Thailand
7 Bureau of Epidemiology, Ministry of Public Health, Nonthaburi 11000, Thailand


Correspondence and requests for materials should be addressed to D.S.B. (dburke@jhsph.edu).


Dengue fever is a mosquito-borne virus that infects 50–100 million people each year1. Of these infections, 200,000–500,000 occur as the severe, life-threatening form of the disease, dengue haemorrhagic fever (DHF)2. Large, unanticipated epidemics of DHF often overwhelm health systems3. An understanding of the spatial–temporal pattern of DHF incidence would aid the allocation of resources to combat these epidemics. Here we examine the spatial–temporal dynamics of DHF incidence in a data set describing 850,000 infections occurring in 72 provinces of Thailand during the period 1983 to 1997. We use the method of empirical mode decomposition4 to show the existence of a spatial–temporal travelling wave in the incidence of DHF. We observe this wave in a three-year periodic component of variance, which is thought to reflect host–pathogen population dynamics5, 6. The wave emanates from Bangkok, the largest city in Thailand, moving radially at a speed of 148 km per month. This finding provides an important starting point for detecting and characterizing the key processes that contribute to the spatial–temporal dynamics of DHF in Thailand.


The incidence of DHF in Thailand varies widely from year to year, showing as much as a tenfold difference between years. Dengue is a leading cause of hospitalization and death among children in Thailand, where all four serotypes of the virus circulate7. Few tools exist to control dengue virus infection and transmission. Control efforts focus on controlling the mosquito vector of the disease, Aedes aegypti, and on effective management of cases of infections1. Reliable prediction of the location and times of high incidence would allow public health systems to allocate their limited resources more effectively.

It is difficult to predict the pattern of DHF over time and geography (that is, the spatial–temporal pattern) owing to the presence of nonstationarity and nonlinearity in incidence data. Several factors are thought to influence the pattern of DHF, including environmental and climate factors8, 9, predator–prey dynamics between the pathogen and the host population5, 6 and viral factors10, 11. Incidence patterns reflect the complex interaction of all of these factors. As a result, incidence data show strong seasonality, multiyear oscillations and changes in period over time. In the present analysis, we use empirical mode decomposition (EMD) to isolate a 3-yr periodic mode of variance. The travelling wave revealed in this periodic mode is obscured in the raw incidence data by the presence of many periodic and roughly periodic components.

An array of spatial–temporal patterns have been observed12-14 and are predicted by theory15, 16 for host–pathogen and predator–prey ecological systems. A repeating, spatial–temporal wave has not been observed, however, in a vector-borne disease of humans. Characterization of the particular spatial–temporal pattern of an ecological system can be used to identify the mechanisms most important in the dynamics of these systems17. One feature of disease systems that has been shown to produce waves in incidence is spatial heterogeneity in the host population12. Spatial temporal travelling waves in measles incidence have been attributed to the reintroduction of measles to small communities through infective sparks from larger communities12. The size of Bangkok's population and its large role in the commerce of the country suggest that, if spatial heterogeneity in the host population is important to DHF dynamics, Bangkok may have a central role.

The method of EMD is analogous to other methods available for processing nonstationary data, such as wavelet analysis and singular spectrum analysis (SST). Unlike wavelet analysis, however, it does not assume a basis a priori. Unlike both wavelet analysis and SST, EMD is appropriate for data describing nonlinear phenomena4. EMD uses an adaptive basis that is derived from each data set to decompose the variance of that set into a finite number of intrinsic mode functions (IMFs)4. These IMFs represent oscillations around a mean at a characteristic timescale of the data—the spacing between extrema. The IMFs are not restricted to a particular frequency or band of frequency, but can experience both amplitude and frequency modulation. An example of the EMD sifting process is shown in Fig. 1. Features suggesting nonlinearity, such as front–back asymmetry in the waveforms, can be seen in the IMFs presented. The decomposition is local, complete and, for practical purposes (but not theoretically), orthogonal4.

  
Figure 1 Example of the EMD sifting process.  Full legend

High resolution image and legend (49k)



Figure 2 shows the monthly incidence of DHF for each Thai province. Brief inspection suggests that there is temporal synchrony across the country, with peaks in incidence occurring roughly at the same time across all provinces. Figure 3 shows the IMFs of roughly 3-yr periodicity for all provinces. EMD analysis of each incidence time series yields, at most, six IMFs. The two most energetic IMFs across all provinces are the IMF associated with seasonal variance and the 3-yr periodic IMF. The 3-yr periodic modes account for 44% (95% confidence interval (C.I.) 39–48%) of the interannual variability in dengue incidence.


Figure 2 Monthly DHF incidence in each of the 72 provinces of Thailand.  Full legend

High resolution image and legend (63k)

  


Figure 3 The 3-yr periodic mode for each of the 72 provinces of Thailand.  Full legend

High resolution image and legend (88k)

  

The nonparametric covariance function18 was used to characterize the spatial synchrony of incidence fluctuations. Spatial synchrony provides a measure of the spatial dependence of temporal correlation among incidence series. Spatial synchrony of the 3-yr mode differs markedly from the synchrony of the raw data. The spatial extent, the distance for which local synchrony is statistically significantly different from the average synchrony across all data, and the global synchrony of these two sets differ greatly (Fig. 4). The spatial extent of synchrony in the raw DHF incidence is about 180 km, whereas the extent of synchrony of the 3-yr mode is about 420 km.

  
Figure 4 Spatial synchrony of DHF incidence (blue) and the 3-yr periodic mode of variance (red) across 72 provinces of Thailand with 95% C.I. envelopes (see Methods).  Full legend

High resolution image and legend (28k)



Spatial synchrony reflects both the timing and relative amplitude of incidence across provinces. By contrast, phase coherence measures only the relative timing of peaks and troughs in incidence12. Phase coherence of the 3-yr mode (see Supplementary Information), although showing a similar spatial extent, is significantly less than spatial synchrony. This suggests that the timing of peaks and troughs is synchronous during large amplitude changes in incidence, but less synchronous during periods of low incidence. The first half of the data, a period of high energy in the 3-yr mode, has both higher spatial synchrony and coherence than the second half.

To examine the role of Bangkok in this wave, cross-correlation functions (CCFs) between the 3-yr mode of incidence in Bangkok and all other provinces were calculated. Figure 5 shows a marked pattern in these functions. The lag at which each of these CCFs are at their maximum absolute value is found to be a function of the distance from Bangkok (P < 10-8; see Supplementary Information). Out of 71 provinces, 68 are either synchronous or lag behind Bangkok. These results describe a repeating, spatial–temporal travelling wave, emanating from Bangkok at a speed of 148 km per month (95% C.I., 114–209 km month-1).


Figure 5 Cross-correlation coefficients between the 3-yr oscillatory mode of DHF incidence in Bangkok and the same mode of DHF incidence in the 71 other provinces of Thailand.  Full legend

High resolution image and legend (69k)



This travelling wave is surprising in its spatial extent. The spatial extent of synchrony of this mode (420 km) and the large number of provinces that are either synchronous or lag behind Bangkok suggest that almost all of the country is affected by this wave. But several border regions do not correlate well with Bangkok. Future research will investigate these border provinces and whether other urban areas in Southeast Asia influence the incidence patterns in these provinces.

To test the robustness of the central role of Bangkok in the dynamics of DHF in Thailand, we repeated the analysis presented here for all 72 provinces. Bangkok has the largest number of provinces lagging behind its incidence pattern (65/71). Three other provinces tie with Bangkok in having the same number of provincial incidence patterns that are either synchronous or lagging behind their own incidence patterns. In a direct comparison with Bangkok, all three of these provinces lag behind Bangkok&#39;s incidence pattern.

The mechanisms underlying this wave are not understood at present. Several classes of mechanism have been implicated in inducing spatial synchrony in other systems, including dispersal12, 19, activator–inhibitor dynamics20 and wide-scale correlation of environmental factors21. For the system studied here, DHF in Thailand, we speculate that immune interactions between the four serotypes may give rise to complex local dynamics in Bangkok5; it has been reported that all four serotypes are continuously present but vary in proportion from season to season22. Similar to the observations for measles in the United Kingdom12, we speculate that smaller communities elsewhere in Thailand may experience periods of no incidence of particular serotypes during some periods (&#39;stochastic fade-out&#39;), and that these serotypes are reintroduced by the migration of viruses from Bangkok. A preliminary analysis has found that the number of months of zero incidence that each province experiences is associated with population size and distance from Bangkok. A previous study detected no role for weather or climate in longer term variances in DHF incidence in Bangkok6; thus, it is unlikely that changes in weather contribute to the country-wide travelling wave.

The development of spatial transmission models of DHF to explore the role of these mechanisms for dengue is an area for future research. Although only narrowing the field of theoretical models17, the spatial–temporal pattern described here provides a much needed crucial test for models of the spatial transmission of dengue. The 3-yr periodic mode seems to modulate both in frequency and in amplitude. It is possible that this modulation has subtle influences on the wave signatures that we have not captured in this analysis.

This 3-yr periodic mode may prove important to predicting periods of high incidence of disease. The low periodicity of this mode facilitates prediction on annual timescales, which will be useful in health system preparations. In addition, surveillance in Bangkok may prove useful for preparations in the surrounding regions, as epidemics are heralded by as much as ten months in some areas. Urban centres have been thought to be important in the genesis of dengue epidemics elsewhere in Asia23. Our results suggest that high priority should be placed on surveillance systems in urban areas of Southeast Asia.

In this analysis, time series decomposition revealed a phenomenon that was not apparent in the raw incidence data or in other modes of variance in the incidence data. Although the isolation of particular modes of variance is a simplification of complex and dynamically interacting disease processes, these techniques can aid the formation of hypotheses and larger models. To our knowledge, this is the first application of EMD to the analysis of epidemiological data. We consider that the ability of the EMD to decompose nonstationary and nonlinear data makes it the most appropriate technique for this analysis.

Methods
Data Numbers of DHF cases have been collected by the Ministry of Health of Thailand since 1972. Cases are diagnosed on the basis of criteria of the World Health Organization (WHO). Serological confirmation is conducted where feasible, but logistically is impossible to do on every case report. Here we consider only the data that were available to us, that is, the monthly incidence of DHF in each of the 72 provinces of Thailand from 1983 to 1997. These data are available on the Johns Hopkins Center for Immunization Research website (http://www.jhsph.edu/cir/dengue.html). We are currently working to obtain incidence data from 1997 to the present.

EMD The method of EMD decomposes a time series into IMFs by means of a sifting process4. The sifting process begins with the identification of local minima and maxima of a raw time series, X(t). Two cubic splines are fit: one connecting the local maxima and one connecting the local minima. The time series of means of these two splines, m1(t), is calculated. The difference between the raw time series and the mean series, X(t) - m1(t), is designated h1(t). The series h1(t) is the first IMF of the data if it satisfies two admission criteria: the number of extrema and zero-crossings must not differ by more than 1, and the mean series between cubic splines connecting the extrema of h1(t) must be 0 at all times. If h1(t) does not satisfy these criteria, the algorithm is repeated using h1(t) as the raw series. The first IMF is subtracted from the raw data series, and the algorithm is repeated on this difference to identify subsequent IMFs of the data.

Log-transformations of incidence time series from each province normalized to have zero mean and unit standard deviation were decomposed24. IMFs with an approximate period (measured by Fourier analysis) of 3–4 yr (the third IMF component for all but one province) were used in subsequent analyses. The Hilbert–Huang Transformation Toolbox for Matlab (Princeton Satellite and NASA, 2001) was used for the EMD analysis.

Spatial synchrony, phase coherence and CCFs Algorithms in the NCF library for R/S-plus (available at http://asi23.ent.psu.edu/) were used to estimate the spatial correlation functions and phase coherence functions. A detailed discussion of these techniques is given in refs 12, 18. Spatial synchrony of two DHF incidence time series was quantified as the Pearson correlation coefficient of those series17. Spatial correlation functions were estimated by the nonparametric spline covariance function18. A cubic B-spline of nine equivalent degrees of freedom was used in this estimation (the square root of the number of provinces was used as a guide18. We calculated C.I.s using 500 bootstrap iterations.

The phase coherence of two DHF incidence time series was calculated as the Pearson correlation coefficient of phase angles of each series12. Again, the nonparametric spline covariance function was used to estimate the phase coherency function12. We calculated CCFs using Pearson correlation coefficients. The lag at which each CCF was at its maximum absolute value (for lags between -12 and 12) was identified for each province.

Supplementary information accompanies this paper.

Received 18 August 2003;accepted 20 November 2003