 Research
 Open Access
Stochastic latticebased modelling of malaria dynamics
 Phong V. V. Le^{1, 2},
 Praveen Kumar^{1, 3}Email authorView ORCID ID profile and
 Marilyn O. Ruiz^{4}
Received: 25 January 2018
Accepted: 22 June 2018
Published: 5 July 2018
Abstract
Background
The transmission of malaria is highly variable and depends on a range of climatic and anthropogenic factors. In addition, the dispersal of Anopheles mosquitoes is a key determinant that affects the persistence and dynamics of malaria. Simple, lumpedpopulation models of malaria prevalence have been insufficient for predicting the complex responses of malaria to environmental changes.
Methods and results
A stochastic latticebased model that couples a mosquito dispersal and a susceptibleexposedinfectedrecovered epidemics model was developed for predicting the dynamics of malaria in heterogeneous environments. The It\(\hat{o}\) approximation of stochastic integrals with respect to Brownian motion was used to derive a model of stochastic differential equations. The results show that stochastic equations that capture uncertainties in the life cycle of mosquitoes and interactions among vectors, parasites, and hosts provide a mechanism for the disruptions of malaria. Finally, model simulations for a case study in the rural area of Kilifi county, Kenya are presented.
Conclusions
A stochastic latticebased integrated malaria model has been developed. The applicability of the model for capturing the climatedriven hydrologic factors and demographic variability on malaria transmission has been demonstrated.
Keywords
 Malaria
 Climate change
 Metapopulation
 Stochastic
 Ecohydrology
Background
Malaria is a vectorborne disease with complex nonlinear dynamics [1, 2]. The disease, caused by protozoan parasites of the genus Plasmodium, is transmitted between humans by female Anopheles mosquitoes. Many factors are determinants for the transmission of malaria, including climate suitability, life cycles of pathogens and vectors, and the local capacity to control the mosquito population [3–6]. The fundamental basis for malaria risk prediction and early warning lies in the ability to estimate the effects of weather, to identify vector habitat, and to model the population dynamics of Anopheles mosquitoes and the transmission of the pathogens in the human population.
Mathematical models have been used to provide an explicit framework for understanding malaria transmission dynamics for over 100 years, starting with the pioneering work of Ross [7]. In a simple form, he used a few ordinary differential equations (ODEs) to describe quantitative changes in densities of infected human and mosquitoes. Since then, more sophisticated models have been developed that include factors such as latent periods [4, 8], vector density and human age structure [4], varying human population size and migration [9, 10], socioeconomic developments [11], temporary immunity [12], and weather effects [13]. In addition, agentbased and metapopulation based models have been used to allow simulations of heterogeneous communities subject to realistic transmission scenarios [14–18]. Over the last 60 years, much scientific research was undertaken and progress made in understanding the biology of malaria vectors and hostparasitevector interactions. Systematic reviews of mathematical modelling of malaria [19] and other mosquitoborne diseases [20, 21] indicate the need for models to address the complexities is the hostvectorparasite interactions and to incorporate heterogeneous environments.
Deterministic models with susceptibleexposedinfectedrecovered (SEIRS) patterns, often based on nonlinear ODEs, are among the standard approaches for estimating the transmission potential for a wide range of infectious diseases, including malaria [22]. These models are useful in understanding the temporal dynamics of infection cycles and in coping with different epidemiological situations including both epidemicity and endemicity. However, SEIRS models ignore the aquatic stage of the mosquito life cycle and the spatial dynamics of mosquitoes when habitats are heterogeneously available, thus limiting the ability to couple and accurately predict the link between the life cycle of Anopheles mosquitoes and malaria transmission. Furthermore, deterministic models cannot capture fluctuations dominated by the random nature of population events, environmental conditions, and variability in the controlling parameters, which inevitably occur in a real system [23]. To that end, stochastic models have proved valuable in estimating asymptotic expressions for the probability of occurrence of major outbreaks as well as stochastic extinctions [24, 25]. Nevertheless, there still exists a lack of complete predictive capability which may allow a mechanism to efficiently capture uncertainty in the temporal and spatial dynamics of malaria. This shortcoming has been attributed to the differences in stochastic modelling approaches, as agentbased models [15, 17, 18] are inefficient for largescale simulations with very large number of vector individuals involved, while lumpedpopulation models [26, 27] bypass the spatial dynamics of the diseases.
In this paper, a stochastic latticebased integrated malaria (SLIM) model for investigating the dynamics of mosquitoes and transmission of malaria in time and space is developed. More precisely, an entomological mosquito dispersal formulation is coupled with a classic SEIRStype model to capture the interaction between malaria dynamics and the life cycle of Anopheles vectors. The model is linked to a sophisticated ecohydrologic model to incorporate climatedriven hydrologic and ecologic processes as factors that determine mosquito population and malaria transmission dynamics. A metapopulation approach is incorporated to describe the movements of vectors among discrete geographic subdomains. The It\(\hat{o}\) approximation of stochastic integrals is used to derive governing stochastic differential equations (SDEs) of the model. Then, the derived SDEs are replicated across lattice grid cells in the domain to capture temporal and spatial dynamics of malaria.
The rest of this paper is organized as follows. In the next section, the mathematical formulations of the stochastic malaria model is described. Then, numerical simulations of this model using topographic and reanalysis data for a case study in Kenya are presented. Finally, the paper closes with discussion of the key points.
Methods: model description
Vector dispersal model
Name  Description  Unit  Range 

b  Integer number of female eggs laid per oviposition  −  50–300 
\(\psi ^\mathsf {W}_{\zeta }\)  50% of the eggs are assumed to hatch into female mosquitoes parameter represent water availability in cell \(\zeta\)  −  0.0–1.0 
\(\psi ^\mathsf {H}_{\zeta }\)  Binary parameter represent human presence in cell \(\zeta\)  −  0–1 
\(\rho _E\)  Egg hatching rate into larvae  day\(^{1}\)  0.33–1.0 
\(\rho _L\)  Rate at which larvae develop into pupae  day\(^{1}\)  0.08–0.17 
\(\rho _P\)  Rate at which pupae develop into adult/emergence rate  day\(^{1}\)  0.33–1.0 
\(\mu _E\)  Egg mortality rate  day\(^{1}\)  0.32–0.80 
\(\mu _{L_1}\)  Natural mortality rate of larvae  day\(^{1}\)  0.30–0.58 
\(\mu _{L_2}\)  Densitydependent mortality rate of larvae  day\(^{1}\)mosq\(^{1}\)  0.0–1.0 
\(\mu _P\)  Pupae mortality rate  day\(^{1}\)  0.22–0.52 
\(\rho _{A_h}\)  Rate at which hostseeking mosquitoes enter the resting state  day\(^{1}\)  0.322–0.598 
\(\rho _{A_r}\)  Rate at which resting mosquitoes enter oviposition searching state  day\(^{1}\)  0.30–0.56 
\(\rho _{A_o}\)  Oviposition rate  day\(^{1}\)  3.0–4.0 
\(\mu _{A_h}\)  Mortality rate of mosquitoes of searching for hosts  day\(^{1}\)  0.125–0.233 
\(\mu _{A_r}\)  Mortality rate of resting mosquitoes  day\(^{1}\)  0.0034–0.01 
\(\mu _{A_o}\)  Mortality rate of mosquitoes searching for oviposition sites  day\(^{1}\)  0.41–0.56 

The rate of change in egg population (\(E_{\zeta }\)) as a function of oviposition, egg mortality and hatching:where b is the average number of eggs laid during an oviposition with 1:1 sex ratio; \(\rho _{A_o}\) is the rate at which eggs are oviposited by gravid mosquitoes; \(\mu _{E}\) is egg mortality rate; and \(\rho _{E}\) is the hatching rate into larvae. The term \(\psi ^\mathsf {W}_{\zeta }\) on the right hand side represents water availability in a particular cell \(\zeta\) and is discussed later.$$\begin{aligned} \frac{dE_{\zeta }}{dt} = b \psi ^\mathsf {W}_{\zeta } \rho _{A_o} A_{o,\zeta }  \left( \mu _{E} + \rho _{E} \right) E_{\zeta } \end{aligned}$$(1a)

The rate of change in larvae population (\(L_{\zeta }\)) as a function of egg population, larval mortality and maturation into pupae:where \(\rho _{L}\) is the progression rate from larvae to pupae; \(\mu _{L_1}\) and \(\mu _{L_2}\) represent natural and densitydependent death rates of larvae, respectively.$$\begin{aligned} \frac{dL_{\zeta }}{dt} = \rho _{E} E_{\zeta }  \left( \mu _{L_1} + \mu _{L_2} L_{\zeta } + \rho _{L} \right) L_{\zeta } \end{aligned}$$(1b)

The rate of change in pupae population (\(P_{\zeta }\)) as a function of larval maturation, pupal mortality, and emergence into adults:where \(\mu _{P}\) is the mortality rate of pupae, and \(\rho _{P}\) represent the rate of emergence from pupae into adults.$$\begin{aligned} \frac{dP_{\zeta }}{dt} = \rho _{L} L_{\zeta }  \left( \mu _{P} + \rho _{P} \right) P_{\zeta } \end{aligned}$$(1c)

The rate of change in population of hostseeking adults (\(A_{h,\zeta }\)) as a function of pupal emergence, oviposition, mortality, and blood feeding rates:where \(\mu _{A_h}\) is the death rate of hostseeking adults; and \(\rho _{A_h}\) is the rate at which they enter a resting state after blood feeding. Hostseeking adults can spread to adjacent cells for searching human host in which \(\omega ^\mathsf {H} _{\zeta _1:\zeta _2}\) represents their movement rate from cell \(\zeta _1\) to \(\zeta _2\) modelled as a decreasing exponential function of human population and average area of cell \(\zeta _1\) and \(\zeta _2\) (see [28]); \(\mathcal {N}\) denotes the neighbors of the cell under consideration. The last two terms in (1d) represent the movements of vectors from all neighbors \(\zeta '\) into cell \(\zeta\) and vice versa, respectively. After laying eggs, gravid mosquitoes return to the hostseeking state for subsequent blood feeding.$$\begin{aligned} \begin{array}{l} \frac{{d{A_{h,\zeta }}}}{{dt}} = {\rho _P}{P_\zeta } + \psi _\zeta ^W{\rho _{{A_o}}}{A_{o\zeta }}  \left( {{\mu _{{A_h}}} + \psi _\zeta ^H{\rho _{{A_h}}}} \right) {A_{h,\zeta }}\\ \qquad \quad+ \left( {\sum \limits _{\zeta ' \in {\mathcal {N}}} {\omega _{\zeta ':\zeta }^H} } \right) {A_{h,\zeta '}}  \left( {\sum \limits _{\zeta ' \in \mathcal{N}} {\omega _{\zeta :\zeta '}^H} \times {A_{h,\zeta }}} \right) \end{array} \end{aligned}$$(1d)

The rate of change in population of resting adults (\(A_{r,\zeta }\)) as a function of blood feeding, mortality, and protein digestion rates:where \(\mu _{A_r}\) is the death rate of resting adults; and \(\rho _{A_r}\) is the progression rate at which the survivors enter the oviposition site searching phase. In the resting state, female mosquitoes are usually dormant to digest protein.$$\begin{aligned} \frac{dA_{r,\zeta }}{dt} = \psi ^\mathsf {H}_{\zeta } \rho _{A_h} A_{h,\zeta }  \left( \mu _{A_r} + \rho _{A_r} \right) A_{r,\zeta } \end{aligned}$$(1e)

The rate of change in population of oviposition site searching adults (\(A_{o,\zeta }\)) as a function of emergence, oviposition, mortality, and digestion rates:where \(\mu _{A_o}\) is the death rate of gravid female mosquitoes. Vectors in ovipositing state also spread in space to find water for oviposition in which \(\omega ^\mathsf {W} _{\zeta _1:\zeta _2}\) represents their movement rate from cell \(\zeta _1\) to \(\zeta _2\) modelled as a decreasing exponential function of surface soil moisture and average area of cell \(\zeta _1\) and \(\zeta _2\) [28].$$\begin{aligned} \frac{dA_{o,\zeta }}{dt}&= \rho _{r} A_{r,\zeta }  (\mu _{A_o} + \psi ^\mathsf {W}_{\zeta } \rho _{A_o}) A_{o,\zeta } \nonumber \\&\quad + \left( \sum _{\zeta ' \in \mathcal {N}} \omega ^{\mathsf {W}}_{\zeta ':\zeta } \right) A_{o,\zeta '}  \left( \sum _{\zeta ' \in \mathcal {N}} \omega _{\zeta :\zeta '}^{\mathsf {W}} \times A_{o,\zeta } \right) \end{aligned}$$(1f)
Although the ELPAs model offers a simple approach to incorporate the effects of vector mobility on spatial population dynamics of vectors in heterogeneous environments, it does not account for fluctuations dominated by randomness in population events and environmental conditions. In some cases, such fluctuations may result in critical state transitions of vector population dynamics. This fact highlights the need to develop stochastic tools to address the complexities arising from vector population dynamics.
Probabilities associated with changes in ELPA\(_{s}\) model
l  Change, \(\Delta X^{l}_{k,\zeta }(t)\)  Probability, \(p^{l}_{k,\zeta }(t)\)  Description 

1  \([1, 0, 0, 0, 0, 0]^T\)  \(b \psi ^{\mathsf {W}}_{\zeta } \rho _{A_o} \mathcal {A}_{o,\zeta } \Delta t\)  A new egg E is deposited by \(A_o\) 
2  \([1, 0, 0, 0, 0, 0]^T\)  \(\mu _E \mathcal {E}_{\zeta } \Delta t\)  An egg E dies 
3  \([1, 1, 0, 0, 0, 0]^T\)  \(\rho _E \mathcal {E}_{\zeta } \Delta t\)  An egg E hatches into a larva L 
4  \([0, 1, 0, 0, 0, 0]^T\)  \((\mu _{L_1} + \mu _{L_2} \mathcal {L}_{\zeta } ) \mathcal {L}_{\zeta } \Delta t\)  A larva L dies 
5  \([0, 1, 1, 0, 0, 0]^T\)  \(\rho _L \mathcal {L}_{\zeta } \Delta t\)  A larva L develops into a pupa P 
6  \([0, 0, 1, 0, 0, 0]^T\)  \(\mu _P \mathcal {P}_{\zeta } \Delta t\)  A pupa P dies 
7  \([0, 0, 1, 1, 0, 0]^T\)  \(\rho _P \mathcal {P}_{\zeta } \Delta t\)  A pupa P develops into a hostseeking adult \(A_h\) 
8  \([0, 0, 0, 1, 0, 1]^T\)  \(\psi ^\mathsf {W}_{\zeta } \rho _{A_o} \mathcal {A}_{o,\zeta } \Delta t\)  An oviposition adult \(A_o\) enters hostseeking state 
9  \([0, 0, 0, 1, 0, 0]^T\)  \(\mu _{A_h} \mathcal {A}_{h,\zeta } \Delta t\)  A hostseeking adult \(A_h\) dies 
10  \([0, 0, 0, 1, 1, 0]^T\)  \(\psi ^\mathsf {H}_{\zeta } \rho _{A_h} \mathcal {A}_{h,\zeta } \Delta t\)  A hostseeking adult \(A_h\) enters resting state 
11  \([0, 0, 0, 0, 1, 0]^T\)  \(\mu _{A_r} \mathcal {A}_{r,\zeta } \Delta t\)  A resting adult \(A_r\) dies 
12  \([0, 0, 0, 0, 1, 1]^T\)  \(\rho _{A_r} \mathcal {A}_{r,\zeta } \Delta t\)  A resting adult \(A_r\) enters oviposition searching state 
13  \([0, 0, 0, 0, 0, 1]^T\)  \(\mu _{A_o} \mathcal {A}_{o,\zeta } \Delta t\)  An oviposition searching adult \(A_o\) dies 
The SELPAs model provides the basis to capture spatial variation of mosquito population dynamics. It incorporates random processes and heterogeneity in both densities of human hosts and breeding sites for the feeding and life cycles of the vectors. SELPAs is coupled with a stochastic SEIRS formulation presented below for malaria transmission dynamics.
Malaria transmission model
The SELPAs model described above extends the deterministic ELPAs model to incorporate stochastic variability associated with population and spatial dynamics of the mosquitoes. Similarly, a stochastic version of SEIRS formulations is developed to capture the variability associated with the circulations of malaria parasites between human and vector populations. The aim is then to combine and link them to an ecohydrological model that explicitly considers climatedriven hydrologic factors for simulating mosquito population dynamics and malaria transmission.

The rate of change in susceptible host (\(S_{h,\zeta }\)) as a function of immigration, birth, human infection, recovery from infection, and human mortality:where \(\Lambda _h\) is immigration rate; and \(\psi _h\) is per capita birth rate of humans; \(\rho _h\) is per capita rate of losing acquired temporary immunity. Acquired temporary immunity represents the enhancement of the defense mechanism of the human host as a result of a previous encounter with the pathogen [32]. \(N_{h,\zeta } = S_{h,\zeta } + E_{h,\zeta } + I_{h,\zeta } + R_{h,\zeta }\) is total population size for humans in each cell \(\zeta\); \(f_h(N_{h,\zeta }) = \mu _{1h} + \mu _{2h} N_{h,\zeta }\) is the human per capita death rate; and \(\lambda _{h,\zeta }\) is the infection rate from mosquitoes to humans defined as:$$\begin{aligned} \frac{dS_{h,\zeta }}{dt} =\,& \Lambda _h + \psi _h N_{h,\zeta } + \rho _h R_{h,\zeta }\\&  \lambda _{h,\zeta }(t) S_{h,\zeta }  f_h(N_{h,\zeta }) S_{h,\zeta } \end{aligned}$$(10a)in which \(N_{v,\zeta } = S_{v,\zeta } + E_{v,\zeta } + I_{v,\zeta }\) is total population of mosquitoes in cell \(\zeta\); \(\sigma _v\) represents the number of times one mosquito attempt to bite humans per unit time; \(\sigma _h\) is the maximum number of mosquito bites a human can have per unit time; and \(\beta _{hv}\) is the probability of infection transmission from an infectious mosquito to a susceptible human, given that a contact between the two occurs.$$\begin{aligned} \lambda _{h,\zeta } = \frac{\sigma _v \sigma _h}{\sigma _v N_{v,\zeta } + \sigma _h N_{h,\zeta }} \times \beta _{hv} I_{v,\zeta } \end{aligned}$$

The rate of change in exposed host (\(E_{h,\zeta }\)) as a function of new host infections, latent period, and human mortality:where \(\nu _h\) is per capita rate of progression of humans from exposed to infectious state.$$\begin{aligned} \frac{dE_{h,\zeta }}{dt} = \lambda _{h,\zeta }(t) S_{h,\zeta }  \nu _h E_{h,\zeta }  f_h(N_{h,\zeta }) E_{h,\zeta } \end{aligned}$$(10b)

The rate of change in infected host (\(I_{h,\zeta }\)) as a function of latent period, recovery and human mortality rates:where \(\gamma _h\) represents per capita recovery rate for humans from infectious to recovered states; and \(\delta _h\) represents per capita diseaseinduced death rate for humans.$$\begin{aligned} \frac{dI_{h,\zeta }}{dt} = \nu _h E_{h,\zeta }  \gamma _h I_{h,\zeta }  f_h(N_{h,\zeta }) I_{h,\zeta }  \delta _h I_{h,\zeta } \end{aligned}$$(10c)

The rate of change in recovered host (\(R_{h,\zeta }\)) as a function of recovery, immunity loss, and human mortality:$$\begin{aligned} \frac{dR_{h,\zeta }}{dt} = \gamma _h I_{h,\zeta }  \rho _h R_{h,\zeta }  f_h(N_{h,\zeta }) R_{h,\zeta } \end{aligned}$$(10d)

The rate of change in susceptible vector (\(S_{v,\zeta }\)) as a function of reproduction, vector infection, and vector mortality:where \(\psi _v\) represents per capita birth rate of the vectors; \(f_v(N_{v,\zeta }) = \mu _{1v} + \mu _{2v} N_{v,\zeta }\) is the per capita death rate for vectors in each cell \(\zeta\); and \(\lambda _{v,\zeta }\) is the infection rate from humans to mosquitoes defined as:$$\begin{aligned} \frac{dS_{v,\zeta }}{dt} = \psi _v N_{v,\zeta }  \lambda _{v,\zeta }(t) S_{v,\zeta }  f_v(N_{v,\zeta }) S_{v,\zeta } \end{aligned}$$(10e)where \(\beta _{vh}\) and \(\tilde{\beta }_{hv}\) represent the transmission probability of infection from an infectious and a recovered human, respectively, to a susceptible mosquito, given that a contact between them occurs.$$\begin{aligned} \lambda _{v,\zeta } = \frac{\sigma _v \sigma _h}{\sigma _v N_{v,\zeta } + \sigma _h N_{h,\zeta }} \times \left( \beta _{vh} I_{h,\zeta } + \tilde{\beta _{vh}} R_{h,\zeta } \right) \end{aligned}$$

The rate of change in exposed vector (\(E_{v,\zeta }\)) as a function of new vector infection, vector latent period, and vector mortality:where \(\nu _v\) is per capita rate of progression of mosquitoes from the exposed state to the infectious state.$$\begin{aligned} \frac{dE_{v,\zeta }}{dt} = \lambda _{v,\zeta }(t) S_{v,\zeta }  \nu _v E_{v,\zeta }  f_v(N_{v,\zeta }) E_{v,\zeta } \end{aligned}$$(10f)

The rate of change in infected vector (\(I_{v,\zeta }\)) as a function of latent period and mortality:$$\begin{aligned} \frac{dI_{v,\zeta }}{dt} = \nu _v E_{v,\zeta }  f_v(N_{v,\zeta }) I_{v,\zeta } \end{aligned}$$(10g)
Name  Description  Unit^{a} 

\(\Lambda _h\)  Immigration rate of humans  H \(\times\) T\(^{1}\) 
\(\psi _h\)  Per capita birth rate of humans  T\(^{1}\) 
\(\psi _v\)  Per capita birth rate of mosquitoes  T\(^{1}\) 
\(\sigma _v\)  Number of times one mosquito would want to bite humans per unit time, if humans were freely available. This is a function of the mosquito’s gonotrophic cycle (the amount of time a mosquito requires to produce eggs) and its anthropophilic rate (its preference for human blood)  T\(^{1}\) 
\(\sigma _h\)  The maximum number of mosquito bites a human can have per unit time. This is a function of the human’s exposed surface area  T\(^{1}\) 
\(\beta _{hv}\)  Probability of transmission of infection from an infectious mosquito to a susceptible human, given that a contact between the two occurs  − 
\(\beta _{vh}\)  Probability of transmission of infection from an infectious human to a susceptible mosquito, given that a contact between the two occurs  − 
\(\tilde{\beta }_{hv}\)  Probability of transmission of infection from a recovered (asymptomatic carrier) human to a susceptible mosquito, given that a contact between the two occurs  − 
\(\nu _h\)  Per capita rate of progression of humans from the exposed state to the infectious state. \(1/\nu _h\) is the average duration of the latent period  T\(^{1}\) 
\(\nu _v\)  Per capita rate of progression of mosquitoes from the exposed state to the infectious state. \(1/\nu _v\) is the average duration of the latent period  T\(^{1}\) 
\(\gamma _h\)  Per capita recovery rate for humans from the infectious state to the recovered state. \(1/\gamma _h\) is the average duration of the infectious period  T\(^{1}\) 
\(\delta _h\)  Per capita diseaseinduced death rate for humans  T\(^{1}\) 
\(\rho _h\)  Per capita rate of loss of acquired temporary immunity for humans. \(1/\rho _h\) is the average duration of the immune period  T\(^{1}\) 
\(\mu _{1h}\)  Densityindependent part of the death (and emigration) rate for humans  T\(^{1}\) 
\(\mu _{2h}\)  Densitydependent part of the death (and emigration) rate for humans  H \(\times\) T\(^{1}\) 
\(\mu _{1h}\)  Densityindependent part of the death (and emigration) rate for mosquitoes  T\(^{1}\) 
\(\mu _{2h}\)  Densitydependent part of the death (and emigration) rate for mosquitoes  M \(\times\) T\(^{1}\) 
Probabilities associated with changes in SEIRS model
l  Change, \(\Delta Y^{l}_{k,\zeta }(t)\)  Probability, \(p^{l}_{k,\zeta }(t)\)  Description 

1  \([1, 0, 0, 0, 0, 0, 0]^T\)  \((\Lambda _h + \psi _h N_{h,\zeta }) \Delta t\)  A new host enters the human susceptible class 
2  \([1, 0, 0, 1, 0, 0, 0]^T\)  \(\rho _h R_{h,\zeta } \Delta t\)  A recovered host becomes susceptible again 
3  \([1, 1, 0, 0, 0, 0, 0]^T\)  \(\frac{\sigma _v \sigma _h \beta _{hv}I_{v,\zeta }S_{h,\zeta }}{\sigma _v N_{v,\zeta } + \sigma _h N_{h,\zeta }}\Delta t\)  A susceptible host enters exposed state 
4  \([1, 0, 0, 0, 0, 0, 0]^T\)  \((\mu _{1h} + \mu _{2h} N_{h,\zeta })S_{h,\zeta } \Delta t\)  A susceptible host dies 
5  \([0, 1, 1, 0, 0, 0, 0]^T\)  \(\nu _h E_{h,\zeta } \Delta t\)  An exposed host enters infectious state 
6  \([0, 1, 0, 0, 0, 0, 0]^T\)  \((\mu _{1h} + \mu _{2h} N_{h,\zeta })E_{h,\zeta } \Delta t\)  An exposed host dies 
7  \([0, 0, 1, 1, 0, 0, 0]^T\)  \(\gamma _h I_{h,\zeta } \Delta t\)  An infectious host enters recovered state 
8  \([0, 0, 1, 0, 0, 0, 0]^T\)  \((\mu _{1h} + \mu _{2h} N_{h,\zeta } + \delta _h) I_{h,\zeta } \Delta t\)  An infectious host dies 
9  \([0, 0, 0, 1, 0, 0, 0]^T\)  \((\mu _{1h} + \mu _{2h} N_{h,\zeta }) R_{h,\zeta } \Delta t\)  A recovered host dies 
10  \([0, 0, 0, 0, 1, 0, 0]^T\)  \(\psi _v N_{v,\zeta } \Delta t\)  A new mosquito enters the vector susceptible class 
11  \([0, 0, 0, 0, 1, 1, 0]^T\)  \(\frac{\sigma _v \sigma _h \beta _{hv}I_{v,\zeta } S_{h,\zeta }}{\beta _{vh} I_{h,\zeta } + \tilde{\beta }_{vh} R_{h,\zeta }} \Delta t\)  A susceptible vector enters exposed state 
12  \([0, 0, 0, 0, 1, 0, 0]^T\)  \((\mu _{1v} + \mu _{2v} N_{v,\zeta }) S_{v,\zeta } \Delta t\)  A susceptible vector dies 
13  \([0, 0, 0, 0, 0, 1, 1]^T\)  \(\nu _v E_{v,\zeta } \Delta t\)  An exposed vector enters infectious state 
14  \([0, 0, 0, 0, 0, 1, 0]^T\)  \((\mu _{1v} + \mu _{2v} N_{v,\zeta }) E_{v,\zeta } \Delta t\)  An exposed vector dies 
15  \([0, 0, 0, 0, 0, 0, 1]^T\)  \((\mu _{1v} + \mu _{2v} N_{v,\zeta }) I_{v,\zeta } \Delta t\)  An infectious vector dies 
The SSEIRS model incorporates environmental perturbation and stochastic interactions among subpopulations of human hosts and vectors in different states. It provides a stochastic framework to study uncertainty and dynamics of malaria transmission as well as other mosquitoborne diseases.
Model couplings
Estimating moisture index
The moisture index \(\psi ^\mathsf {W}_{\zeta }\) shown in (1a) represents water availability in a particular cell \(\zeta\) and is estimated using a sophisticated ecohydrologic modelling framework (Dhara, see [36]). The Dhara framework includes a collection of canopy process models (MLCan, see [37–39]) and a physicallybased surfacesubsurface flow model coupler (GCSflow, see [40]) designed for capturing moisture transport on the land surface and in the belowground systems. It incorporates vegetation acclimation to elevated CO\(_2\) and the retention of moisture flow dynamics associated with topographic variability. This integration provides predictive capability to capture the impacts of environmental changes on the formation and persistence of breeding habitat. The SLIM (coupled SELPAs and SSEIRS) model is linked with Dhara for incorporating climatedriven hydrologic factors to vector population and malaria transmission dynamics.
Model performance
In order to evaluate the performance of SLIM, each of its components (SELPAs and SSEIRS) is first analysed independently. Then simulations of the fully coupled SLIM model using observed meteorological and topographic data for a case study in the rural area in Kilifi county, Kenya are presented.
SELPAs model
Figure 4a shows the variations in logscale of adult mosquito density averaged over \(\mathbb {D}_1\) under the effects of random processes. It can be seen that SELPAs simulations with larger populations are affected slightly by stochastic noises, and the dynamics tend to be close to equilibria predicted by the deterministic systems. In contrast, smaller population sizes experience proportionally more noise and their behaviors tend to be further from the deterministic systems, highlighting the different dynamics of vector population at different sizes. The variations of small vector population can also be disrupted significantly by sudden changes induced by stochastic noises in the system. Linking these similarities and differences between stochastic and deterministic systems in metapopulation and latticebased models is thus important to study the dynamics of vector density in large areas, where populations at various sizes are interconnected.
SSEIRS model
A large number (\(\sim\)100) of independent SSEIRS simulations are conducted and compared with the deterministic SEIRS model to investigate the modelled dynamics of malaria transmission in noisedominated systems. All the stochastic and deterministic simulations have the same initial conditions chosen randomly. Moreover, periodic boundary conditions are applied for all simulations. Parameters for the two models are the same and selected from a previous published study [10]. Further, to separate the effects of stochastic noise on the dynamics of malaria, the dependences of parameters on hydroclimatic factors are also excluded in both SSEIRS and SEIRS simulations. As the spatial movement of mosquitoes is not considered in SSEIRS model, only simulations for a single grid cell are implemented. These simulations are conducted in 800 days with daily time step as well.
The variation of infected human malaria cases (\(\mathcal {I}_h\)) obtained from stochastic SSEIRS (shaded area) and deterministic SEIRS (blue solid line) models are shown in Fig. 4b. Unlike the deterministic approach, SSEIRS simulations provide a range of possibilities for malaria transmission given the same initial conditions and model parameters. SSEIRS simulations shown in Fig. 4b highlight the random nature of malaria transmission that inevitably occur in real systems. Although the mean values of stochastic simulations (red solid line) are found to be close to those in the deterministic simulation, variability obtained from SSEIRS implies that there are probabilities that (i) local outbreaks can be disrupted by random noises or stochastic extinctions may occur and (ii) the intensity of the outbreaks can be larger than the theoretical endemic equilibrium point shown in deterministic models. Note that periodic boundary conditions applied in the simulations imply isolated systems. This may allow shorter persistence time to extinction of malaria than in nonisolated systems. In our simple tests, the extinction of malaria in a specific region (i.e. entire domain) is defined as events in which the number of infected human population in this region is smaller than 0.5. Moreover, the probability of resurgence of the diseases in isolated systems is low. Capturing the range of variability and possible stochastic extinctions plays a fundamental role in understanding and breaking the circulation of the pathogen. This information, usually ignored in deterministic approaches, is important for preparing malaria control in reality.
Case study
Figure 6a–d presents the variation of total aquatic and adult phases of mosquito populations obtained from the SELPA component in SLIM. The results reveal that the variation of the mosquito population in both aquatic and adult stages are highly dependent on climatic factors. Specifically, positive correlations are found between monthly averaged mosquito populations modelled in SELPAs with observed monthly air temperature (\(R^2 = 0.75  0.86\)) and rainfall (\(R^2 = 0.69  0.77\)), respectively. The largest and smallest total mosquito population during the years are found correspondingly with the highest and lowest air temperature and rainy seasons with several days of time lag (10 and 18 days, respectively). In aquatic stages, the sensitivity of larvae development to air temperature change is found to be much lower than of eggs and pupae which was also shown in previous studies [29]. The population of adult Anopheles mosquitoes is also sensitive to climatic conditions (see Fig. 6d). The result shows that the fraction of hostseeking mosquitoes (\(A_h\)) in adult stage is high, consisting of \(\sim 7080\%\) of the total adult population. The subpopulation of oviposition site searching mosquitoes (\(A_o\)) are usually \(23\) times larger than the resting mosquitoes (\(A_r\)). The high number of egg deposited by female Anopheles during reproduction is likely a key factor that explains the high density of vectors in aquatic environment, thus mosquito population. The total number of adult mosquitoes is equal to the number of adults in the SEIRS model and plays a key role in malaria transmission.
The dynamics of malaria in human hosts and mosquito populations in the study area are presented in Fig. 6e–f. Positive correlations are also found between the monthly average malaria incidence (\(\mathcal {E}_h\), \(\mathcal {I}_h\)) with air temperature (\(R^2 = 0.58  0.69\)) and rainfall (\(R^2 = 0.53  0.67\)), respectively. The results show that, similarly to the vector population, the variation of malaria incidences, including both exposed and infected cases, in the region is sensitive to climatic factors as it is directly dependent on vector density. The largest values of exposed human cases (\(E_h\)) are usually found after the rainy seasons start and air temperature was high. The peaks of \(E_h\) are also followed by the largest values of infected human cases (\(I_h\)) in several days (Fig. 6e). During the peaks and troughs of the season, the rates of infected cases are about 2.5 and 1.0% of total population, respectively.
Conclusions
In summary, we have presented a stochastic latticebased integrated malaria (SLIM) model that consists of a vector dispersal (SELPAs) and a malaria dynamics (SSEIRS) component. SLIM is developed to predict mosquito population dynamics and malaria transmission in response to heterogeneity and variability in the environment. It is well known that climatic and hydrological conditions strongly control Anopheles mosquito populations and thus influence malaria incidence, and indeed the associations have been demonstrated repeatedly. The details of malariaenvironment interactions are highly nonlinear and uncertain both in time and space, thus the optimal predictive ability arises from complex models that involve processes from hydroclimatology, ecology, and entomology.
The stochastic coupled model is constructed based on deterministic systems [9, 10, 28]. The model is also link with a an ecohydrologic model (Dhara) used to capture soil moisture dynamics on the ground. This integration provides the capability to incorporate climatedriven hydrologic and ecologic processes with the dynamics of vector population and malaria transmission. In this manner, the presented SLIM model augments existing models by explicitly simulating all of the aforementioned complexities and incorporating a range of possible outcomes to the dynamics of vector population and transmission.
Declarations
Authors' contributions
PVVL, PK designed research. PVVL, PK, MOR carried out the modelling. PVVL analysed data. PVVL, PK, MOR wrote the paper. All authors read and approved the final manuscript.
Acknowledgements
We would like to thank people in Kumar research group at Ven Te Chow Hydrosystem Laboratory for their help and support with this study. The work also used the ROGER supercomputer, which is supported by NSF grant number ACI 1429699.
Competing interests
The authors declare that they have no competing interests.
Availability of data and materials
The SLIM model is publicly available at https://github.com/HydroComplexity.
Consent for publication
Not applicable.
Ethics approval and consent to participate
Not applicable.
Funding
PVVL received support from Computational Science and Engineering (CSE) fellowship. PK received support from NSF (CBET1209402, ACI 1261582, EAR 1331906, EAR 1417444).
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Authors’ Affiliations
References
 Miller LH, Baruch DI, Marsh K, Doumbo OK. The pathogenic basis of malaria. Nature. 2002;415:673–9.View ArticlePubMedGoogle Scholar
 Smith DL, McKenzie EF. Statics and dynamics of malaria infection in Anopheles mosquitoes. Malar J. 2004;3:1–14.View ArticleGoogle Scholar
 Anderson RM. The population dynamics of infectious diseases: theory and applications. Population and community biology series. London: Chapman & Hall Ltd.; 1982.View ArticleGoogle Scholar
 Anderson RM, May RM. Infectious diseases of humans: dynamics and control. Dynamics and control. Oxford: Oxford University Press; 1992.Google Scholar
 Paaijmans KP, Thomas MB. Health: wealth versus warming. Nat Clim Change. 2011;1:349–50.View ArticleGoogle Scholar
 Caminade C, Kovats S, Rocklov J, Tompkins AM, Morse AP, ColónGonzález FJ, et al. Impact of climate change on global malaria distribution. Proc Natl Acad Sci USA. 2014;111:3286–91.View ArticlePubMedGoogle Scholar
 Ross R. The prevention of malaria. 2nd ed. Dutton; 1910.Google Scholar
 MacDonald G. The Epidemiology and Control of Malaria. Oxford Medical Publications. Oxford, UK: Oxford University Press; 1957.Google Scholar
 Ngwa GA, Shu WS. A mathematical model for endemic malaria with variable human and mosquito populations. Math Comput Model. 2000;32:747–63.View ArticleGoogle Scholar
 Chitnis N, Cushing J, Hyman J. Bifurcation analysis of a mathematical model for malaria transmission. SIAM J Appl Math. 2006;67:24–45.View ArticleGoogle Scholar
 Yang HM. Malaria transmission model for different levels of acquired immunity and temperaturedependent parameters (vector). Rev Saude Publica. 2000;34:223–31.View ArticlePubMedGoogle Scholar
 Filipe JAN, Riley EM, Drakeley CJ, Sutherland CJ, Ghani AC. Determination of the processes driving the acquisition of immunity to malaria using a mathematical transmission model. PLoS Comput Biol. 2007;3:e255.View ArticlePubMedPubMed CentralGoogle Scholar
 Parham PE, Michael E. Modeling the effects of weather and climate change on malaria transmission. Environ Health Perspect. 2010;118:620–6.View ArticlePubMedGoogle Scholar
 Ariey F, Duchemin JB, Robert V. Metapopulation concepts applied to falciparum malaria and their impacts on the emergence and spread of chloroquine resistance. Infect Genet Evol. 2003;2:185–92.View ArticlePubMedGoogle Scholar
 Bomblies A, Duchemin JB, Eltahir EAB. Hydrology of malaria: model development and application to a Sahelian village. Water Resour Res. 2008;44:W12445.View ArticleGoogle Scholar
 Gu W, Novak RJ. Agentbased modelling of mosquito foraging behaviour for malaria control. Trans R Soc Trop Med Hyg. 2009;103:1105–12.View ArticlePubMedPubMed CentralGoogle Scholar
 Arifin SN, Zhou Y, Davis GJ, Gentile JE, Madey GR, Collins FH. An agentbased model of the population dynamics of Anopheles gambiae. Malar J. 2014;13:1–20.View ArticleGoogle Scholar
 Pizzitutti F, Pan W, Barbieri A, Miranda JJ, Feingold B, Guedes GR, et al. A validated agentbased model to study the spatial and temporal heterogeneities of malaria incidence in the rainforest environment. Malar J. 2015;14:1–19.View ArticleGoogle Scholar
 Mandal S, Sarkar R, Sinha S. Mathematical models of malaria—a review. Malar J. 2011;10:202.View ArticlePubMedPubMed CentralGoogle Scholar
 Reiner RC, Perkins TA, Barker CM, et al. A systematic review of mathematical models of mosquitoborne pathogen transmission: 1970–2010. J R Soc Interface. 2013;10:20120921.View ArticlePubMedPubMed CentralGoogle Scholar
 Smith DL, Perkins TA, Reiner RC, Barker CM, Niu T, Chaves LF, et al. Recasting the theory of mosquitoborne pathogen transmission dynamics and control. Trans R Soc Trop Med Hyg. 2014;108:185–97.View ArticlePubMedPubMed CentralGoogle Scholar
 Keeling MJ, Rohani P. Modeling infectious diseases in humans and animals. Princeton: Princeton University Press; 2008.Google Scholar
 Azaele S, Maritan A, Bertuzzo E, RodriguezIturbe I, Rinaldo A. Stochastic dynamics of cholera epidemics. Phys Rev E. 2010;81:051901.View ArticleGoogle Scholar
 Herwaarden OA, Grasman J. Stochastic epidemics: major outbreaks and the duration of the endemic period. J Math Biol. 1995;33:581–601.View ArticlePubMedGoogle Scholar
 van Herwaarden AO. Stochastic epidemics: the probability of extinction of an infectious disease at the end of a major outbreak. J Math Biol. 1997;35:793–813.View ArticlePubMedGoogle Scholar
 Britton T. Stochastic epidemic models: a survey. Math Biosci. 2010;225:24–35.View ArticlePubMedGoogle Scholar
 Krstic M. The effect of stochastic perturbation on a nonlinear delay malaria epidemic model. Math Comput Simul. 2011;82:558–69.View ArticleGoogle Scholar
 Lutambi AM, Penny MA, Smith T, Chitnis N. Mathematical modelling of mosquito dispersal in a heterogeneous environment. Math Biosci. 2013;241:198–216.View ArticlePubMedGoogle Scholar
 Depinay JM, Mbogo C, Killeen G, Knols B, Beier J, Carlson J, et al. A simulation model of African Anopheles ecology and population dynamics for the analysis of malaria transmission. Malar J. 2004;3:29.View ArticlePubMedPubMed CentralGoogle Scholar
 Allen E. Modeling with Itô Stochastic differential equations. Mathematical modelling: theory and applications. Heidelberg, Germany: Springer Berlin Heidelberg; 2007.Google Scholar
 Allen LJS. An introduction to stochastic processes with applications to biology. 2nd ed. Florida: CRC Press; 2010.Google Scholar
 Doolan DL, Dobaño C, Baird JK. Acquired immunity to malaria. Clin Microbiol Rev. 2009;22:13–36.View ArticlePubMedPubMed CentralGoogle Scholar
 Detinova TS. Agegrouping methods in Diptera of medical importance with special reference to some vectors of malaria. WHO Monograph series. 1962;47:13–191.PubMedGoogle Scholar
 Briere JF, Pracros P, Le Roux AY, Pierre JS. A novel rate model of temperaturedependent development for arthropods. Environ Entomol. 1999;28:22–9.View ArticleGoogle Scholar
 Paaijmans KP, Read AF, Thomas MB. Understanding the link between malaria risk and climate. Proc Natl Acad Sci USA. 2009;106:13844–9.View ArticlePubMedGoogle Scholar
 Le PVV, Kumar P. Interaction between ecohydrologic dynamics and microtopographic variability under climate change. Water Resour Res. 2017;53:8383–403.View ArticleGoogle Scholar
 Drewry DT, Kumar P, Long S, Bernacchi C, Liang XZ, Sivapalan M. Ecohydrological responses of dense canopies to environmental variability: 1. Interplay between vertical structure and photosynthetic pathway. J Geophys Res. 2010;115:G04022.Google Scholar
 Le PVV, Kumar P, Drewry DT, Quijano JC. A graphical user interface for numerical modeling of acclimation responses of vegetation to climate change. Comput Geosci. 2012;49:91–101.View ArticleGoogle Scholar
 Le PVV, Kumar P, Drewry DT. Implications for the hydrologic cycle under climate change due to the expansion of bioenergy crops in the Midwestern United States. Proc Natl Acad Sci USA. 2011;108:15085–90.View ArticlePubMedGoogle Scholar
 Le PVV, Kumar P, Valocchi AJ, Dang HV. GPUbased highperformance computing for integrated surfacesubsurface flow modeling. Environ Modell Softw. 2015;73:1–13.View ArticleGoogle Scholar
 Le PVV, Kumar P. Power law scaling of topographic depressions and their hydrologic connectivity. Geophys Res Lett. 2014;41:1553–9.View ArticleGoogle Scholar
 Nyakeriga AM, TroyeBlomberg M, Chemtai AK, Marsh K, Williams TN. Malaria and nutritional status in children living on the coast of Kenya. Am J Clin Nutr. 2004;80:1604–10.View ArticlePubMedGoogle Scholar
 Snow RW, Kibuchi E, Karuri SW, Sang G, Gitonga CW, Mwandawiro C, et al. Changing malaria prevalence on the Kenyan Coast since 1974: climate, drugs and vector control. PLoS ONE. 2015;10:1–14.Google Scholar
 Mogeni P, Williams TN, Fegan G, Nyundo C, Bauni E, Mwai K, et al. Age, spatial, and temporal variations in hospital admissions with malaria in Kilifi County, Kenya: a 25year longitudinal observational study. PLoS Med. 2016;13:1–17.View ArticleGoogle Scholar
 Le PVV, Kumar P, Ruiz MO, Mbogo C, Muturi JE. Predicting the direct and indirect impacts of climate change on malaria in coastal Kenya. PLOS (under review). 2018.Google Scholar
 Tatem AJ, Noor AM, von Hagen C, Di Gregorio A, Hay SI. High resolution population maps for low income nations: combining land cover and census in East Africa. PLoS ONE. 2007;2:1–8.View ArticleGoogle Scholar
Comments
By submitting a comment you agree to abide by our Terms and Community Guidelines. If you find something abusive or that does not comply with our terms or guidelines please flag it as inappropriate. Please note that comments may be removed without notice if they are flagged by another user or do not comply with our community guidelines.