An ecohydrological model of malaria outbreaks

Malaria is a geographically widespread infectious disease that is well known to be affected by climate variability at both seasonal and interannual timescales. In an effort to identify climatic factors that impact malaria dynamics, there has been considerable research focused on the development of appropriate disease models for malaria transmission driven by climatic time series. These analyses have focused largely on variation in temperature and rainfall as direct climatic drivers of malaria dynamics. Here, we further these efforts by considering additionally the role that soil water content may play in driving malaria incidence. Specifically, we hypothesize that hydro-climatic variability should be an important factor in controlling the availability of mosquito habitats, thereby governing mosquito growth rates. To test this hypothesis, we reduce a nonlinear ecohydrological model to a simple linear model through a series of consecutive assumptions and apply this model to malaria incidence data from three South African provinces. Despite the assumptions made in the reduction of the model, we show that soil water content can account for a significant portion of malaria’s case variability beyond its seasonal patterns, whereas neither temperature nor rainfall alone can do so. Future work should therefore consider soil water content as a simple and computable variable for incorporation into climate-driven disease models of malaria and other vector-borne infectious diseases.


Introduction
The World Health Organization estimates that 250 million clinical episodes of malaria occur annually, resulting in at least one million disease-associated deaths (World Health Organization, 2008).Malaria incidence is especially high in developing countries, where it is a leading cause of morbidity and mortality, in particular among children and pregnant women.Since the pioneering work by Ross (1910) and MacDonald (1957), progress in understanding malaria dynamics has been made through the development of mathematical models and their statistical inference with incidence data (e.g.Bailey, 1982;Hay et al., 2002;Depinay et al., 2004;Zhou et al., 2004;Pascual et al., 2008;Chaves et al., 2012;Bhadra, 2011;Laneri et al., 2010, among others).A subset of these models has considered the role that external forcing plays in generating patterns of seasonal and interannual case variability.Despite these advances, early warning systems for malaria have still only limited ability (and thereby efficacy) to predict outbreaks, and the factors contributing to malaria case variability still require more thorough investigation (Pascual et al., 2008;Craig et al., 2004a,b).
Two climatic variables that have long been known to influence malaria's seasonal and interannual dynamics are temperature and rainfall.Temperature is known to affect the development time of mosquito larvae, the probability of mosquito survival, and the development time of the malaria parasite Plasmodium falciparum in infected mosquitoes (Bayoh and Lindsay, 2003;Hoshen and Morse, 2004).Rainfall is hypothesized to affect malaria incidence through the creation of mosquito breeding sites during wet periods and the reduction of their availability during droughts (Patz et al., 1998).Although a moderate level of rainfall appears to have a positive effect on mosquito recruitment, intense rainfall events may destroy mosquito habitats and thereby reduce malaria incidence shortly following their occurrence (Briet et al., 2008).

E. Montosi et al.: Ecohydrology and malaria outbreaks
Due to the multiple effects of temperature and rainfall on the malaria parasite and its mosquito vector, previous work linking malaria incidence to climate forcing has frequently focused on the direct impacts of temperature and rainfall on the risk of malaria (Craig et al., 2004a;Lindasy et al., 2000).Several studies have found a correlation between malaria incidence and either minimum, mean or maximum temperatures (Hay et al., 2002;Zhou et al., 2004;Devi and Jauhar, 2007;Gomez-Elipe et al., 2007), while others have instead considered the effects of diurnal temperature fluctuations (Paaijmans et al., 2009) and indices of the upper atmospheric circulation combined with sea-surface temperatures (Jury and Kanemba, 2007) on the disease.At regional spatial scales, malaria incidence has been correlated with rainfall amounts (Pascual et al., 2008;Craig et al., 1999;Hay et al., 2001), interestingly with time lags of several months, consistent with the presence of a physical "buffering" mechanism that dampens rainfall fluctuations.These latter findings hint at the fundamental role that soil water content dynamics may play in malaria outbreaks.Recent work has moved in this direction, with the role of surface hydrology being more fully recognized and explicitly included in mechanistic models.Some of these models have focused on water level fluctuations in small reservoirs and the effect of these fluctuations on malaria prevalence (Porphyre et al., 2005;Shaman et al., 2002).In an effort to predict the number of malaria cases, several statistical models have also considered vegetation density (which depends on soil water availability) and distances between water bodies and human populations (Gomez-Elipe et al., 2007;Kleinschmidt et al., 2001).Finally, detailed models linking hydrology with entomology have recently been proposed and tested for a semi-arid region (Bomblies et al., 2008;Yamana et al., 2011).The role of hydrologic processes in the dynamics of other vector-borne infectious diseases is also starting to become more fully recognized (Bertuzzo et al., 2008).
Motivated by these previous efforts, here we test the hypothesis that soil water content is an important driver of malaria dynamic.To test this hypothesis, we develop a minimalist ecohydrological model of coupled surface water malaria dynamics and use it to describe malaria incidence in three South African provinces (Fig. 1).Although the ecohydrological model we derive makes a series of simplifying assumptions, our results show that variability in soil water content is significantly correlated with variability in malaria incidence, whereas neither rainfall nor temperature alone shows this correlation (beyond their seasonal associations).Future work should therefore focus on soil water content as a simple and computable environmental variable that can be incorporated into more mechanistic models of malaria transmission that include internal feedbacks.

Malaria data
Monthly malaria incidence data over the period July 1996-March 2007 were obtained from the South African Department of Health (http://www.health-e.org.za/resources/statistics.php,last access: February 2012) for three South African provinces: Limpopo, Mpumalanga and KwaZulu-Natal (Figs. 2-3).Implementation of malaria control strategies, mainly due to the Lumbobo Spatial Development Initiative, resulted in a steady decrease in the number of reported malaria cases in the provinces of Mpumalanga and Limpopo starting in 2005 and in the province of KwaZulu-Natal starting in 2002.To neglect this transient and consider a time series driven only by climate variability, we therefore limited our analysis to incidence data before June 2005 for Mpumalanga and Limpopo and before July 2001 for KwaZulu-Natal (Figs. 2-3).

Climate data
The meteorological data consisted of daily rainfall and daily minimum and maximum temperature records collected from weather stations managed by the South African Weather Service.We restricted our analysis to datasets with less than 20 % missing meteorological data that were obtained from stations in areas with high population densities and intermediate to high malaria risk (Fig. 1).For each province, spatially integrated daily time series were first obtained by averaging across the selected weather stations.Finally, to have the same temporal resolution for the malaria and climate datasets, we computed monthly meteorological data from these spatially integrated daily time series (Figs.2-3).A preliminary inspection of monthly malaria data alongside the meteorological data highlighted different responses to temperature and rainfall.We observed that months characterized by one or more days with daily maximum temperatures above 39 • C tended to be followed by a decrease in the number of malaria cases, most likely resulting from the effect of heat-stress on mosquitoes (Craig et al., 1999).Periods of moderate to high precipitation had a positive and delayed effect on malaria incidence, presumably because of the increased availability of mosquito breeding sites (Wyse et al., 2007).The few instances of anomalously high daily precipitation were followed by a rapid decrease in malaria incidence, presumably a result of mosquito habitat destruction (Briet et al., 2008).Moreover, Limpopo province shows a shorter delay in malaria occurrence with respect to the provinces of Mpumalanga and KwaZulu-Natal (insets of Figs.2-3).

Model description
We assessed the association between climatic drivers and malaria cases in three ways: (i) standard linear regression between climate anomalies and malaria case anomalies, at time lags of zero and one month; (ii) an ecohydrological model of malaria dynamics; and (iii) a transfer function model accounting for delayed climate effects on malaria dynamics.The first approach (i) provided a baseline to quantify the improvement of models (ii) and (iii).We chose not to consider the monthly maximum temperature data as an explanatory variable for malaria cases, because it did not have significant explanatory power in preliminary analysis.Instead, we focused our analysis on monthly data of both minimum temperature and precipitation anomalies.
The full ecohydrological model explicitly coupled a model of malaria transmission dynamics with a hydrologic model describing soil water content.Through a series of assumptions, detailed below, this full ecohydrological model was reduced to a minimal, linear model that describes the expected relationship between precipitation levels in previous months and current malaria incidence.
The full model of malaria transmission is given by the following equations, where M is the population size of the Anopheles mosquito vector, M I is the population size of infected mosquitoes, H S is the population size of susceptible individuals, and H I is the population size of infected individuals (Ross, 1910;Porphyre et al., 2005;Kermack and McKendrick, 1927;Smith and McKenzie, 2004;McCallum et al., 2001): The number of individuals recovered and immune to malaria infection at time t is not explicitly modeled, but given by H TOT − H S − H I where H TOT is the constant total host population size.
The dynamics of the mosquito population, given by Eq. ( 1), are modeled by a delay differential equation where 0 is the mosquito growth rate (a function of soil moisture, w, and temperature, T ), 1 δ (∼10 days) is the average lifespan of a mosquito, and τ M (∼10 days) is the time delay between oviposition and mosquito emergence (Chitnis et al., 2008).The number of infected mosquitoes, governed by Eq. (2), increases through feeding on infected hosts at a rate α (whose unit is (months • number of infected individuals) −1 as we assumed a temporal resolution of 1 month), with τ I being a time delay (in months) representing the incubation period for malaria parasites, and decreaseing with background mortality at a rate δ.The number of susceptible hosts increases with births (µH TOT ) and loss of immunity (γ (H TOT − H S − H I ), where γ is dimensionless) and decreases with background mortality (µH S ) and malaria transmission via infected mosquitoes ( η 0 M I H S H TOT , Eq. 3).Finally, the number of infected hosts increases with transmission and decreases with background mortality and recovery from infection (υH I ).µ, η 0 and υ are all quantities expressed in terms of (months) −1 .
Soil water content, w, is described by a soil-moisture and surface-water balance equation (Wyse et al., 2007;Porporato et al., 2004): where P (t) is rainfall at time t, and mw(t) is a linearized soil water loss function accounting for plant transpiration, surface evaporation, and deep infiltration, where m is a lumped soil water loss rate (with units of months −1 ).The use of this function is justified by the large spatial scale (Fig. 1) and relatively low temporal resolution (∼1 month) we are considering.In contrast to classical models of soil water content, we chose a form for Eq. ( 5) that does not saturate, so as to account for both soil and surface water storage.By surface water storage, we mean primarily soil moisture, but also seasonal ponds and rivers that might be breeding sites for mosquitoes.In the following, the assumptions employed to reduce the full ecohydrological model (i.e., Eqs.1-4) for malaria dynamics and Eq. ( 5) for surface hydrology to the linear model are described.
-Assumption 1: Because the number of infected is much smaller than the total population in the considered areas, as a first approximation we can assume that the susceptible fraction of the host population, H S H TOT , is relatively stable through time.With this assumption, the dynamics in the number of infected individuals are simplified to where η = η 0 H S H TOT is expressed in (months) −1 .-Assumption 2: Consistent with the available malaria data, we focus on the monthly time scale and we assume that both the time delay between mosquito oviposition and emergence, and the incubation period, can be neglected.These assumptions yield -Assumption 3a: when climate is the only limiting factor to mosquito emergence, the effect of M(t) on the growth rate can be neglected, effectively approximating the exponential growth resulting from Eq. ( 7) with a linear growth with moisture-and temperature-controlled rate, which can be expressed as -Assumption 3b: alternatively, in the case that the mosquito population exhibits logistic growth, with climate determining its carrying capacity, we can re-write Eq. ( 7) as where r is the maximum growth rate and K the carrying capacity of the mosquito population.
-Assumption 4: we assume that the total mosquito density M is approximately in equilibrium at the monthly time scale dM(t) dt = 0, on the grounds that both the soil water content (with a mean transit time of 1-3 months) and the size of the infected population (transit time ∼10 months) fluctuate at longer time scales than mosquito density (transit time 1 month) (Chitnis et al., 2008).With Assumption 3a in place, the equilibrium mosquito population, M, is found as Alternatively, with Assumption 3b in place, -Assumption 5: since we are considering a large geographic area where spatial heterogeneity likely weakens the nonlinear nature of the interaction between infected humans and mosquitoes (as shown for other systems; see e.g.Katul et al., 2007b), we also assume that the density of infected mosquitoes is proportional to the total mosquito density: M I (t) ∼ = γ M(t).With these assumptions, M I (t) is given by and from Eqs. ( 11) and ( 12), respectively.Equations ( 13) and ( 14), together with Eq. ( 5) and ( 6), lead to a system of two coupled linear equations driven by rainfall P , water content w, and temperature T : where -Assumption 6: the effects of environmental conditions on mosquito growth rate are assumed to be linear, i.e., A(t) = a + bw(t) + cT (t).
Under these assumptions, the minimalist model given by Eqs. ( 15) and ( 16) is linear in both the state variables and the climatic factors, allowing us to remove the seasonal component of these dynamics, leaving the dynamics themselves in terms of anomalies from the seasonal averages.This solution is appropriate because the disease and climate time series are both strongly seasonal, therefore hindering the identification of how climate variability drives disease dynamics on an interannual time scale (Hay et al., 2000;Briet et al., 2008).Seasonal averages for malaria case incidence and for precipitation levels are shown in the insets of Figs.2-3 for each of the three provinces.
Indicating the monthly anomalies with prime signs and the seasonal monthly averages with an overbar, we have P (t) = P + P (t), A(t) = A + A (t) and H (t) = H + H (t). Because the monthly temperature anomalies are negligible in comparison to the anomalies in the hydrologic variables (the maximum standard deviation of the anomalies of temperature is 0.84 • C, recorded for the province of Limpopo, and the maximum standard deviation of the temperature seasonal average is 4.19 • C in the province of KwaZulu-Natal), the effect of temperature anomalies on A can be neglected.Accordingly, Eqs. ( 15) and ( 16) can be re-written as with Given monthly malaria case incidence anomalies A and monthly precipitation anomalies P , we sought to fit parameters m and b.To this end, Eq. ( 17) was integrated numerically with a given value of m through a finite difference approach to yield soil water content levels over time, w (t), and parameter b was then estimated through linear regression (Eq.19).
Our expectation is that b assumes positive values (anomalously high water soil content should have a positive effect on mosquito growth rate) and that the rate of soil water loss m is on the order of 0.3-0.5 months −1 (Katul et al., 2007a).To yield predictions of monthly malaria cases (Fig. 4), malaria case anomalies A were first estimated (given the best fit parameter values b and m) and then added to the seasonal averages in malaria case numbers A. By looking at the monthly series of malaria and precipitation, we realized that there is a high seasonality both in the data per se and in their standard deviations.Therefore, we developed the aforementioned ecohydrological model to describe these high seasonal dynamics mainly looking at the deseasonalized data by means of a physically based approach.However, it is interesting to assess the impact of the seasonality of the standard deviation of the data as well.To this aim, we developed a stochastic model, our third approach, to study the association between climatic drivers and malaria cases.Specifically, we used a transfer function model, i.e. a linear filter of the climatic variables to predict malaria cases.The model was conceived by mimicking the relationship suggested by the ecohydrological model (Eqs.17 and 19 above).
Given that the time series of both monthly malaria cases (A) and monthly rainfall rates (P ) showed a seasonal pattern in their standard deviations, making the time series heteroscedastic and limiting robust parameter estimation, we adopt a transfer function model of the following form: where double prime signs indicate that the variables P , w , and A are normalized with respect to their seasonal standard deviation.Equation ( 20) describes soil moisture anomaly at month i as rainfall anomaly plus a portion of the previous month anomaly, while Eq. ( 21) is analogous to Eq. ( 19).Parameters d and n are conceptually similar to b and m.Again, we expect that d is positive.To obtain predictions of monthly malaria cases using this transfer function model (Fig. 4), malaria case anomalies A were estimated (based on the best fit parameter values d and n), denormalized with respect to standard deviations and, as before, added to the seasonal averages in malaria case numbers A. To check statistical consistency, transfer function model residuals were tested against Gaussianity, homoscedasticity, and independence, using the Kolmogorov-Smirnov, Bartlett and Portmanteau test, respectively.Calibration for both models was performed by using the differential evolution adaptive Metropolis (DREAM) algorithm (Vrugt et al., 2008), yielding best fit values for the parameters along with their posterior probability density functions.DREAM algorithm convergence is monitored through the R-statistic of Gelman and Rubin (1992) and with a visual inspection of the likelihood behavior of each individual chain.To approximate the posterior probability density functions of the model parameters, a total of 80 000 model evaluations are performed with DREAM; the first 40 % of the samples of each chain are discarded and considered as burn-in, in order to evaluate the marginal densities of the parameters from the samples whose likelihood is close to convergence.Figure 5 shows the posterior probability distributions for the parameters of both models applied to the Mpumalanga data.It can be seen that calibration converged to well-defined optimal values.17) and ( 19)).Best fit parameters are reported in Table 2.

Results and discussion
The seasonal component of malaria dynamics explained on average 41 % of the variability in cases (Table 1).None of the linear regressions between deseasonalized climate forcing and deseasonalized malaria cases explained significant variability in malaria case anomalies, resulting in at most a Nash efficiency (Nash et al., 1970) of 14 % for the rainfall-malaria correlation with a zero month time lag in the province of KwaZulu-Natal (Table 1).Introducing time lags improved the results only for the case of Mpumalanga province when one considers the rainfall-malaria correlation with one month time lag.These results indicate that a direct and linear association between climatic variables and malaria dynamics is not statistically supported, beyond seasonal effects.We therefore considered whether soil water content dynamics may be associated with malaria case anomalies.The ecohydrological model was able to explain a significant fraction of the variation in malaria cases in all three provinces (Fig. 4 and Table 1).Furthermore, the parameter estimates of the ecohydrological model, including the slope of the regression between soil water content anomalies and malaria case anomalies, are biologically and physically interpretable (Table 2).Estimates of the parameter m, the rate at which surface water is lost, yielded values between 0.1 and 0.44 months −1 , although the related 95 % confidence limits, also shown in Table 2, show the presence of significant parameter uncertainty.The optimal values of parameter m translate into soil water transit times of 2-10 months, which highlight the persistence of perturbations induced in malaria cases by higher than usual precipitation.95 % confidence limits of the transit times do not allow one rejecting the hypothesis that they are identical in all the three locations.Values of the soil water loss rate m between 0.3 and 0.5 months −1 are comparable with estimates from other systems (Katul et al., 2007a), while the smaller estimate of this parameter for Limpopo suggests the possible presence of a biological delay in the disease dynamics not explicitly considered here (Pascual et al., 2008;Hay et al., 2000) and might be further imputable to Limpopo's different climatic conditions.As a whole, from Table 1.Nash efficiencies for the seasonal component of malaria dynamics, the preliminary analysis, the ecohydrological model, and the transfer function model for the three South African provinces of Mpumalanga, Limpopo and KwaZulu-Natal.In the preliminary analysis, we examined the direct link between monthly data of malaria and rainfall, as well as malaria and minimum temperature applying linear regressions to the data before the removal of the seasonality denoted by "Non-deseasonalized" in the table and after the deseasonalization, considering only the "Anomalies" and taking into account the seasonal dynamics as well (Seasonal + anomalies).a physical perspective, these results suggest that unusually wet periods result in anomalously high surface and soil water storage, and thereby higher mosquito densities, ultimately yielding anomalously high malaria case numbers.
In our ecohydrological model, all the delays between favorable environmental conditions and malaria outbreaks are assumed to be induced by soil water dynamics.This is consistent with observed biological delays (in both the mosquito population and in disease dynamics) that are shorter than a month (e.g.Chitnis et al., 2008), while the memory induced by surface water storage is generally longer (e.g.Katul et al., 2007a).The proposed simplified model neglects compound delays and nonlinear interactions that might result in long delays between climatic forcing and malaria cases (Pascual et al., 2008;Hay et al., 2000).Our results are therefore likely to underestimate surface water loss rates m to some degree, since we attributed all physical and biological delays to soil water content dynamics alone.
In the transfer function statistical model, we expect different parameter values, because deseasonalization of standard deviation was carried out.However, the resulting transit times are consistent with those obtained with the ecohydrological model.The results show that the hypotheses of Gaussianity and homoscedasticity of the residuals cannot be re-jected at the 95 % confidence level in the three time series, therefore providing support for the assumption of model linearity.On the other hand, all the residuals exhibited a slight but still statistically significant autocorrelation (at the 95 % confidence level).This outcome is due to the inability of the model and the selected climatic determinants to fully account for the persistence of the malaria case anomalies.Such persistence could presumably be accounted for by other factors; however, this would introduce more parameters and thus uncertainties in their estimation.Alternatively, the residual autocorrelation could be eliminated introducing a further autoregressive component in the regression model to account for previous malaria cases, which would summarize the effects of the additional inputs above (see, for instance, the approach adopted by Zhou et al., 2004).However, this solution was not used here because it induces equifinality, therefore hindering the efficient identification of the causal relationship between rainfall/soil water content and malaria.Looking at the statistical model performances (Table 1), one can see that they are comparable with those of the ecohydrological model, with some differences depending on location.
It is not straightforward to provide a physical interpretation of the spatial variability for model parameter values.Additional information on the topography and climate   2.
of the study area can be obtained from the maps provided by the South African Department of Environmental Affairs and Tourism (see: http://www.environment.gov.za/?q= content/maps-graphics, last access: May 2012).The areas prone to malaria risk are located mainly in the low-lying parts of the three provinces where the climate is quite variable, as it is possible to assess from the maps of mean annual precipitation provided by the aforementioned department.
It is moreover possible to see from the insets of Figs.2-3 that the seasonal incidence of malaria in the province of KwaZulu-Natal is almost double with respect to the other two provinces.This high incidence of malaria cases is in line with the higher population counted in this province by the 2001 census (Statistics South Africa, 2003) and is also reflected by the results we obtained.In fact, parameter m is relatively stable across the three provinces (Table 2), indicating relatively mild changes in soil and vegetation properties at the large scale we are considering.In contrast, the parameter b varies remarkably.In particular, a higher value is obtained for the KwaZulu-Natal province.Because b quantifies the sensitivity of malaria cases to rainfall, the higher values in KwaZulu-Natal province may be related to its coastal proximity, which implies the presence of a more humid climate.Moreover, wetlands and water bodies are characteristically more frequent in the land cover map of this province, if compared to what is depicted in Limpopo and Mpumalanga maps.A previous work by Kleinschmidt et al. (2001) carried out in two districts located in the northern side of KwaZulu-Natal province confirms the relevant role of wetlands and water bodies.In line with our results, they stressed the importance played by proximity to permanent water bodies, which are often surrounded by ephemeral mosquito habitats, as an additional risk factor for malaria incidence.

Conclusions
Using time series of malaria cases from three South African provinces, we showed through two sets of analyses -an ecohydrological model and a transfer function model -that soil water content is an important driver of malaria dynamics.These analyses required a series of consecutive assumptions to be made in order to interface these models with available climate and malaria data.Nevertheless, we found a statistical association between modeled soil water content and malaria cases in all three provinces, with parameter estimates that were biologically and physically interpretable.Future work, with more extensive time series, should therefore focus on coupling soil water content dynamics to full epidemiological models (Eqs.1-4), which require a larger number of parameters, but that maintain the non-linear feedbacks that are known to be important for malaria dynamics in endemic regions (Koelle and Pascual, 2004;Pascual et al., 2008;Laneri et al., 2010).Further developments should consider satellite soil moisture data as a possible alternative to direct estimates of soil moisture.Using such data would provide high resolution maps without invoking a water balance model dedicated to this purpose as Brugger et al. (2011) already proposed for the study of the seasonal cycle of Culicoides spp.abundance.

Fig. 1 .
Fig. 1.Map of the investigated provinces in South Africa, colored by level of malaria risk.The locations of the weather stations for each province are marked with black triangles (see: http: //www.malaria.org.za/MalariaRisk/Risk Maps/risk maps.html).

Fig. 3 .
Fig. 3. Monthly malaria cases (solid gray line) and minimum temperature (dashed black line) for (A) Mpumalanga, (B) Limpopo, and (C) KwaZulu-Natal provinces.Insets describe seasonal averages of malaria cases and minimum temperature (time proceeds according to the arrows).

Fig. 4 .
Fig. 4. Observed and modeled monthly malaria incidence for (A) Mpumalanga, (B) Limpopo and (C) KwaZulu-Natal provinces.Observed incidence is shown alongside seasonal averages, case estimates from the ecohydrological model, and case estimates from the transfer function model.Insets are scatter plots of soil water content anomalies against malaria case anomalies for each province (dots, observations; solid lines, modeled relationships using Eqs.(17) and (19)).Best fit parameters are reported in Table2.

Fig. 5 .
Fig. 5. Posterior parameter distributions for (A) ecohydrological and (B) transfer function models at Mpumalanga.95 % confidence bands for the parameters are given inTable 2.