Environ. Sci. Technol. 2008, 42, 2643–2649
Coupling Hydrologic and Infectious Disease Models To Explain Regional Differences in Schistosomiasis Transmission in Southwestern China JUSTIN REMAIS,* SONG LIANG, AND ROBERT C. SPEAR Center for Occupational and Environmental Health, School of Public Health, 50 University Hall, University of California, Berkeley, Berkeley, California 94720-7360
Received May 4, 2007. Revised manuscript received December 7, 2007. Accepted December 11, 2007.
Rainfall-runoff models have become essential tools for conceptualizing and predicting the response of hydrologic processes to changing environments, but they have rarely been applied to challenges facing health scientists. Yet with their efficient parameterization and modest data requirements, they hold great promise for epidemiological application. A modeling analysis incorporating simple hydrologic constraints on transmission of the human parasite Schistosoma japonicum in southwestern China was conducted by coupling a lumped parameter rainfall-runoff model (IHACRES) with a delay-differential equation schistosomiasis transmission model modified to account for channel flows and on-field egg inactivation. Model predictions of prevalence and infection timing agree with observations in the region, which indicate that hydrological differences between sites can lead to pronounced differences in transmission. Channel flows are shown to be important in determining infection intensity and timing in modeled village populations. In the periodic absence of flow, overall transmission intensity is reduced among all modeled risk groups. However, the influence of hydrologic variability was greater on the cercarial stage of the parasite than the miracidial stage, due to the parasite ova’s ability to survive dormant on fields between rain events. The predictive power gained from including hydrological data in epidemiological models can improve risk assessments for environmentally mediated diseases, under both long-term climate change scenarios and near-term weather fluctuations.
Introduction Hydrological processes have been shown to strongly influence transmission of water-related infectious diseases. The interannual transmission of malaria has been associated with precipitation and river height anomalies in Peru, Guyana, and Colombia (1). Analyses of cholera cases and remotely sensed data have shown that increases in sea surface heights (SSH) correlate positively with cholera incidence, possibly due to inland intrusion of contaminated water during SSH increases (2). Above-average precipitation accompanying a 1991–1992 El Niño event is believed to have initiated an outbreak of Hantavirus pulmonary syndrome in the U.S (3). * Corresponding author phone: +1 510 642 9016; fax: +1 510 642 5815; e-mail:
[email protected]. 10.1021/es071052s CCC: $40.75
Published on Web 02/20/2008
2008 American Chemical Society
These empirical relationships serve as useful starting points; however, a mechanistic understanding of the links between transmission and intervening hydrological processes would aid in estimating the risk of water-related disease in space and through time. While rare, several studies have used hydrologic models to reveal underlying mechanisms linking hydrology and disease. In a series of papers, Shaman et al. (4–6) construct a topographically based hydrologic model combining a module simulating the vertical movement of water and heat between the soil, vegetation, and the atmosphere with a TOPMODEL (7) component tracking the horizontal movement of shallow groundwater along elevation gradients. Modeled surface wetness was a significant predictor of mosquito abundance using both statistical (6) and mechanistic (5) models and was empirically correlated with human St. Louis encephalitis cases (4). The hydrologic model was not dynamically coupled to a transmission model, however; thus, predictions of disease incidence were not presented. Patz et al. (8) estimated soil moisture using a dynamic hydrologic model and then associated these estimates with key malaria transmission parameters, including Anopheles human-biting rate (HBR) and entomological inoculation rate (EIR). Statistical models including soil moisture conditions, estimated using a simple water balance model, considerably improved prediction of HBR and EIR over rainfall data alone. The authors did not estimate the parameters’ influence on disease incidence. Such studies ultimately aspire to combine their hydrologic models and associated parameters with epidemiologic models (5) but have been hindered by the complexities involved in coupling disease dynamics, hydrological processes, vector habitat, and abundance. Schistosomiasis, a disease long associated with irrigation and dam projects, offers a unique opportunity to explore this coupling, due to its straightforward dependence on surface water for transmission. Schistosomes are water-based parasites that spend part of their life cycle in aquatic snails, part as free-living aquatic larvae, and part in human and other mammalian hosts. Schistosomiasis is the second most prevalent tropical disease after malaria, with 600 million people at risk of infection globally and 200 million infected (9). Schistosoma japonicum (S. japonicum), the species causing the disease in east and southeast Asia, threatens nearly 30 million people in the tropical and subtropical zones of the People’s Republic of China (PRC; 10). In the mountainous region of Sichuan Province, PRC, infected humans shed S. japonicum ova in their excreta, transmitting the parasite to snails through use of manurebased fertilizers. Infected snails shed a larval stage of the parasite, cercariae, which penetrates the intact skin of humans and other mammals that contact contaminated waters. Irrigation canals, the principal snail habitat in this region, serve as the dominant site of water contact activities (11). Timing and intensity of infections is controlled by a number of time-varying factors termed gating functions, described elsewhere (12). These factors, including weather variables, hydrological dynamics, seasonal water-contact patterns, and irrigation patterns, limit incidence to the period roughly between April and October. Our work in southwestern China exploring the environmental determinants of schistosomiasis transmission has highlighted the importance of hydrology to transmission. Maszle et al. (13) estimated the hydrological parameters (e.g., the Chezy coefficient) of typical irrigation systems in the Xichang VOL. 42, NO. 7, 2008 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2643
TABLE 1. Calibration Parameters and Calibration and Validation Model Performance
CQa XCb
rs
βs
βq
τs
υs
υq
τw
f
NSEcc
Bc
ARPE
NSEvc
Bvd
-0.997 -0.995
0.002 0.004
0.449 0.139
303.91 193.60
0.551 0.861
0.449 0.139
146.9 15
4.0 0.1
0.85 0.26
3.1 × 10-3 -1.3 × 10-2
0.007 0.037
0.18 0.09
2.5 × 10-3 -1.8 × 10-4
a CQ, Changqiu. validation.
b
XC, Xichang.
c
NSE subscripts: c, calibration; v, validation.
area and the travel time of the free-swimming forms of the parasite due to advective transport. In a follow-up study, Lowe et al. (14) showed empirically that there is significant transport of viable parasite larvae within irrigation channels, with viable organisms detectable as far as 400 m from source snails. Using these parameters, the impact of larval transport on endemic disease transmission was assessed using a spatial-temporal model of networked villages (15), showing that diffusion of larvae via the surface water pathway influences not just the intensity of transmission but also the effectiveness of interventions. The motivation for the work reported here arose from dissimilarities in epidemiological field data between the Xichang study sites and villages in a separate endemic area some 400 km distant, in Sichuan’s Changqiu Mountains. There are minimal differences in climate, topology, and agriculture between the two sites. While historical data indicate at least some periods when disease prevalence has been similar (16), environmental measures of parasite presence in recent years have been strikingly different. A mouse bioassay procedure has shown cercarial density, and consequent infection risk, in irrigation systems in Changqiu villages to be episodic over the infection season, in contrast with the essentially continuous infection risk, with a higher annual average, found in Xichang villages. Infection was present in just 4% (n ) 404) of sentinel mice sampled in Changqiu throughout the 2004–2005 infection seasons (unpublished data), compared with 14% (n ) 807) observed in Xichang (11). Annual average infection intensity among positive mice differed by 2 orders of magnitude between sites. Likewise, infected snails were located in nearly 3% (n ) 4101) of Xichang snail sites (11), compared with 0.1% (n ) 1363) of sites in Changqiu villages (unpublished data). Field observations (discussed below) suggested that surface water dynamics underlie this difference in risk. Mathematical models accounting for the dominant mechanistic aspects of transmission have served as valuable tools for exploring the influence of environmental parameters on disease (17). Here, we explore the interaction between local hydrological dynamics and schistosomiasis transmission using a transmission model modified to include a lumped parameter hydrologic module. We first develop a predictive runoff model that generates the surface water flow behaviors of the two different study regions, with rainfall and temperature data serving as inputs. Flow predictions are then used as inputs to a mechanistic schistosomiasis transmission model, which explores the degree to which hydrological differences explain observed differences in infection risk between the two locations.
Materials and Methods Study Sites. The Xichang study area (E102°18′ N27°52′) is an approximately 12 × 12 km endemic region in southern Sichuan, characterized by a subtropical climate, with an annual average temperature of 18 °C and annual rainfall of about 1000 mm. Nearly 400 km north of Xichang and 100 km southwest of the capital city of Chengdu, Changqiu (E103°36′ N30°12′) is a 15 × 15 km endemic region with a somewhat cooler subtropical climate, an annual average temperature of 16.4 °C, and annual rainfall of about 1000 mm. Both regions 2644
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 42, NO. 7, 2008
d
B (m3 s-1) subscripts: c, calibration; v,
are mountainous and dominated by irrigated agricultural land. Networks of irrigation channels crisscross their heavily terraced terrain, providing key habitat for the snail host (11). The availability of water in irrigation channels differs dramatically between the regions. Changqiu has low or no baseflow, punctuated by high peak events during the wet season. In Xichang, where local water managers pump water from nearby Qionghai Lake, baseflow is consistent yearround, with occasional wet season surges. Changqiu relies mainly on rainwater and small reservoirs, which are tapped once or twice yearly to flood local paddies. These anecdotal observations led us to pursue quantification of channel flow differences and explore how these might impact infection dynamics. Hydrologic Model Development. A conceptual rainfallrunoff model (IHACRES; 18, 19) was used to estimate daily channel flow in representative channels in the two sites. The IHACRES model, having been applied to diverse catchments worldwide including small and ephemeral catchments (20–23), is particularly attractive in this application because it is parametrically efficient. It consists of a nonlinear rainfall loss module, which converts observed rainfall (mm) into effective or excess rainfall, and a linear module, which converts excess rainfall into observed streamflow (m3 s-1; Supporting Information Figure S1). The nonlinear module permits linearization of the rainfall-runoff relationship by considering only that part of the observed rainfall that contributes to flow. The IHACRES implementation used here, as well as the goodness of fit criteria used to evaluate model performance, are fully detailed in the Supporting Information. IHACRES Calibration and Validation. Two irrigation channels (one each in Changqiu and Xichang; both approximately 2 m wide at high water) were selected for stage monitoring on the basis of the presence of active transmission. Continuously logging pressure transducers (Hobo Onset U20-001-01), compensated for temperature drift and barometric pressure, were installed at the bottom of each channel to measure stage. Flow measurements were made at multiple flow volumes, and a simple rating curve was constructed to estimate flow (m3 s-1) from the stage. One year of daily aggregated flow estimates was available for the analysis described herein. Daily precipitation and temperature data were collected from local rain gauges and thermometers. The calibration procedure, following Jakeman et al. (19), involves a fine grid sampling of nonlinear model parameters, each set producing a candidate effective rainfall series. A simple refined instrumental variable (SRIV) technique was used to estimate linear module parameters. The preferred model is identified by the set of parameters that define the global maximum of Nash-Sutcliffe efficiency (NSE), while minimizing bias (B) and average relative parameter error (ARPE). Calibration of the IHACRES model (version 2.1; 24) was undertaken using a daily time step with rainfall and flow data spanning the late spring through fall of 2005 in Xichang and Changqiu (Table 1). The calibration interval was selected to represent the range of flows experienced at the sites. Validation of the best fitting models was undertaken with rainfall and flow data in the spring of 2006. The authors acknowledge the limitations of calibrating and validating over
FIGURE 1. Modeled and observed flow, the residual error and observed rainfall during the calibration interval in Changqiu.
FIGURE 2. Modeled and observed flow, the residual error and observed rainfall during the calibration interval in Xichang. short time intervals, but short periods have been used successfully in the past (25). Further, our aim of identifying a dynamic model that reproduces qualitative surface water behavior at the sites, rather than meeting quantitative modeling requirements (such as hydrograph volume prediction), is compatible with the limited availability of data. NSE values show model performance in Changqiu exceeded that in Xichang. Low ARPEs for both models indicate that linear module parameter identification during calibration was successful. While validation fits are qualitatively reasonable, NSE for both validation models suffered from limited available data. Figure 1 and Figure 2 show the modeled and observed flow, residual error, and observed rainfall during the calibration interval in Changqiu and Xichang, respectively; Supporting Information Figure S2 depicts their respective performance during model validation. The IHACRES model occasionally missed some flow events; others were indicated in the simulation but not observed. Spatial variability in precipitation not captured by rain gauges could have led to these anomalies. Bias in validation runs was tolerable (less than 5% average flow) for both models. Xichang’s flow pattern resembles that of other small catchments with persistent baseflow (26). Changqiu’s flow pattern resembles that of a more ephemeral catchment (27); in fact, reservoir releases may account for the minimal baseflow observed there. As anticipated, the Changqiu flow regime as captured by the model is characterized by strong quickflow patterns, accounting for nearly half (45%) the total flow volume. In contrast, the Xichang flow regime is dominated by baseflow, accounting for about 86% of the total flow volume. Both sites exhibit a long recession time for baseflow (304 days for Changqiu; 194 for Xichang), characteristic of agricultural water supplies with extensive terraces and contour cropping (28). This may indicate a minimum flow allowance introduced by reservoir storage, particularly during the dry season. Considerable variance
in the Xichang flows in early spring coincides with the planting season, which may indicate pulsed reservoir releases for row crops (Supporting Information Figure S2). Spring crop planting takes place later in Changqiu, which may contribute to the absence of early spring flow variance there. Rice paddy flooding, occurring in late spring in Xichang, may account for dips observed in late April flows. Changqiu’s row crop and orchard dominated agriculture may preclude these large extractions (Supporting Information Figure S2). Integration into Schistosomiasis Transmission Model. A schistosomiasis transmission model following the Anderson and May framework (29) was modified to incorporate local environmental characteristics vital to disease transmission in Sichuan. It is a nonlinear, differential-equation model with state variables representing the parasite’s mean village infection intensity in three human risk groups (farmer, student, and others, which includes domestic workers, teachers, and administrators, etc.), and the village’s infected snail population. Structure and parametrization of the model are described elsewhere (30); the basic framework is summarized in the Supporting Information. The calibrated IHACRES model was coupled to the transmission model as follows. Flow modifies C(t) (Supporting Information eq S9), the mean spatial cercarial density, and M(t) (Supporting Information eq S11), the mean miracidial density, in the irrigation system at time t. As a modifier of C(t), flow predictions (qt) can provide estimates of rc(t) in Supporting Information eq S9, the precipitation-dependent modulation of the average daily cercarial production entering the aquatic environment, where, at time t, we propose the following simple binary formulation: rc(t) )
{
0 qt < Ωc 1 q t g Ωc
(1)
If flow falls below the threshold for cercarial release, Ωc, then rc(t) ) 0, effectively prohibiting cercarial penetration of hosts. When the flow threshold is met or exceeded, rc(t) ) 1 and transmission proceeds unimpeded. Interpretation of this formulation of rc(t) is straightforward: when water is unavailable in irrigation channel environments, exposure to cercariae, and therefore new infections, cannot occur. During and after rain events, when flowing water is available as predicted by the IHACRES submodel, cercarial dispersion and penetration can occur. This formulation is consistent with the ecology of Oncomelanian snails, which reside above the waterline but are submerged and shed cercariae when channel flows rise (31). A similar formulation is suggested for egg/miracidial dynamics, except that, unlike cercariae, ova are resilient under environmental stress. Average cercariae lifetime is approximately 10–12 h (32); schistosome eggs can persist for days in the environment. In animal feces, S. japonicum eggs can survive up to 4 months, dependent on temperature (33–35). Desiccation, however, reduces egg survival to a week or less (36). A composite parameter representing all inactivation processes was calculated from previous experimental data, resulting in a first-order decay constant of 0.2 day-1. A simple first-order inactivation process was implemented in which the miracidia term, M(t) (Supporting Information eq S11), is a function of the sum of decaying eggs contributed since the last flow event: T
E *(t) )
∑ E(t)e
-d(T-t)
(2)
T-τ
where E*(t) is the sum of viable eggs shed by infected humans since the last flow event; d is the decay constant governing inactivation of eggs lying dormant on fields between flow events; and T - τE is the time since the last flow event. VOL. 42, NO. 7, 2008 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
2645
FIGURE 3. Predicted flows in Changqiu (top) and Xichang (bottom) normalized as a fraction of the range (maximum-minimum) of the series.
FIGURE 4. Predicted worm burden, w, over a 4 year period with the year-round presence of flow in channels (Xichang, left) and accounting for the absence of flow in channels (Changqiu, right) using a binary rc(t) function allowing approximately 60 transmission days/year.
Analogous to the cercarial equation, at time t: re(t) )
{
0 qt < Ωe 1 q t g Ωe
(3)
Here, if water flow falls below the threshold for egg hatching, Ωe, then re(t) ) 0 and eggs lie dormant. When the flow threshold is met or exceeded, re(t) ) 1, and viable eggs on fields are washed into the irrigation system, where they hatch and can infect snails. Parameterization of the transmission model followed a Monte Carlo procedure described in the Supporting Information. Simulations in the current work rely on the calibration reported by Liang et al. (37). Calibration runs were conducted for the village Shian 5 for the interval of June 2000-October 2002, and outputs were classified using field data collected during that period. Rainfall and temperature data from 2002 to 2004 were used to predict channel flows using the calibrated IHACRES models. Flow predictions were normalized as qk′ )
qk - qmin qmax - qmin
(4)
where 0 > qk′ > 1. Normalized predictions for this interval are depicted in Figure 3 and agree with our field observations at the sites. Thresholds Ωc and Ωe were set low to disallow transmission under nearly standing-water conditions. Predicted daily flows were classified as insufficient (
