 Research
 Open Access
 Published:
Stochastic latticebased modelling of malaria dynamics
Malaria Journal volume 17, Article number: 250 (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.
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,4,5,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,15,16,17,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
SLIM consists of two stochastic timecontinuous, spacediscrete models: one for vector dispersal that represents the entomological life cycles of female Anopheles, and a second for SEIRS malaria dynamics that simulates the circulation of the pathogen between host and vector populations (see Fig. 1). The model is driven by the distribution of human population (density, locations), atmospheric forcing (air temperature, precipitation), and hydrological conditions (soil moisture distribution and ponding persistence in topographic depressions). The spread of malaria in time and space is modelled as follows. Let \(\mathbb {D}\) be a bounded domain in \(\mathbb {R}^2\) and let \(m,n \ge 1\) be integers. \(\mathbb {D}\) is partitioned into uniform and rectangular lattice grid indexed as \(\{\zeta := (i,j) \in \mathbb {Z}^2_{+} : i \le n, j \le m\}\). In each grid cell \(\zeta\), a vector dispersal model for mosquito population dynamics is coupled with a SEIRS formulation for malaria transmission. It is assumed that mosquitoes in aquatic stages are immobile, but adult mosquitoes can move between adjacent grid cells, spreading the pathogen.
Vector dispersal model
The deterministic vector dispersal model of Lutambi et al. [28] that serves as the basis for constructing the stochastic formulation is described first. The life cycle of the female mosquito has six distinct phases, including three aquatic stages: egg (E), larval (L), pupal (P); and three adult stages: hostseeking (\(A_{h}\)), resting (\(A_{r}\)), and oviposition site searching (\(A_{o}\)). Each phase of this model, hereafter referred to as ELPAs, is subject to fluctuations due to recruitment, mortality, and progression of survivors into the next states. \(A_{h}\) and \(A_{o}\) are further affected by the movement of mosquitoes in space. For \(t > t_0\) and \(\forall \zeta \in \mathbb {D}\), ELPAs is given as a set of dynamical equations as described below. The summary of key parameters presented in ELPAs and their ranges of values are given in Table 1.

The rate of change in egg population (\(E_{\zeta }\)) as a function of oviposition, egg mortality and hatching:
$$\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)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.

The rate of change in larvae population (\(L_{\zeta }\)) as a function of egg population, larval mortality and maturation into pupae:
$$\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)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.

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

The rate of change in population of hostseeking adults (\(A_{h,\zeta }\)) as a function of pupal emergence, oviposition, mortality, and blood feeding rates:
$$\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)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.

The rate of change in population of resting adults (\(A_{r,\zeta }\)) as a function of blood feeding, mortality, and protein digestion rates:
$$\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)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.

The rate of change in population of oviposition site searching adults (\(A_{o,\zeta }\)) as a function of emergence, oviposition, mortality, and digestion rates:
$$\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)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].
The above model (1) is wellposed in a positively invariant domain:
Unlike the original model [28], environmental variability is further incorporated to the developmental rates of aquatic mosquitoes using the temperature relationships [29]:
for \(k=\{ E, L, P\}\), where \(T_a\) is the mean air temperature (K) over the time interval \(\Delta t\); and \(r_k(T_a) : \mathbb {R} \mapsto \mathbb {R}\) is the development rate per unit time at \(T_a\). Details for the temperature dependencies of egg, larval, and pupal populations are presented in Depinay et al. [29].
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.
A continuous SDE model for the dynamical system (1) is developed in the following way (see Chapter 5, ref [30] for more details about the approach to stochastic modelling applied herein). First, it is assumed that there is demographic variability due to births, deaths, and transitions between the states in (1). Second, a discrete (Markov chain) stochastic model for (1) is constructed. For small interval \(\Delta t\), we then identify all possible changes and their corresponding transition probabilities for the discrete stochastic processes. Third, the expected changes and the covariance matrix of changes of these processes are determined. Finally, the continuous SDE model for (1), hereafter referred to as SELPAs, is inferred by similarities in the forward Kolmogorov equations between the discrete and continuous stochastic processes [30, 31]. Note that solutions of the discrete Markov chain and continuous SELPAs models approximately have the same probability distribution.
Let \(\mathcal {E}_{\zeta }(t)\), \(\mathcal {L}_{\zeta }(t)\), \(\mathcal {P}_{\zeta }(t)\), \(\mathcal {A}_{h,\zeta }(t)\), \(\mathcal {A}_{r,\zeta }(t)\), \(\mathcal {A}_{o,\zeta }(t) \in \mathbb {R}_{\ge 0}\) denote continuous random variables for the density of eggs, larvae, pupae, hostseeking adults, resting adults, and oviposition site searching adults in a grid cell \(\zeta\), respectively (see Fig. 2). The SELPAs model depends on the state variables:
where \(\mathbf{{X}}_{k}(t) = \{X_{k,\zeta }(t) \in \mathbb {R} : \zeta \in \mathbb {D} \} \subseteq \mathbb {R}^{m \times n}_{\ge 0}\) for \(k=1,\ldots ,6\) in which \(\mathbf{{X}}_{k}(t)\) has an associated probability density function \({{p}}_{k}(x, t)\):
Let \(\Delta {\mathbf{X}}(t) = {\mathbf{X}}(t+\Delta t)  {\mathbf{X}}(t)\). For small \(\Delta t\), there are 13 possible unit changes in \(\Delta X_{k,\zeta }(t)\) of the discrete stochastic model associated with different probabilities in each cell \(\zeta\) (see Table 2). The It\(\hat{o}\) stochastic differential equations for the SELPAs model are given as follow [30]:
where \({\mathbf{F}}^x : \mathbb {R}^{6\times m \times n} \rightarrow \mathbb {R}^{6\times m \times n}\); \({\mathbf{G}}^x : \mathbb {R}^{6\times 13 \times m \times n} \rightarrow \mathbb {R}^{6 \times 13 \times m \times n}\); \({\mathbf{W}}(t) \subset \mathbb {R}^{13 \times m \times n}\) is the matrix of independent Wiener processes; and \({\mathbf{p}}^x \subseteq \mathbb {R}^{m \times n}\) is the matrix of parameters that are functions of time t and climatic, anthropogenic, and entomological factors \({\mathbf{C}^x}(t) \subseteq \mathbb {R}^{m \times n}\).
The drift term \({\mathbf{F}}^x\) in (5), to order \(\Delta t\), is the expectation of all possible changes in the discrete stochastic model computed as:
for \(k = 1,\ldots,6\) and \(\forall \zeta \in \mathbb {D}\).
Additional Wiener processes are included into the stochastic systems to simplify the derivation of the diffusion term \({{\varvec{G}}}^x\) [30, 31], written as follows:
in which the covariance matrix associated with the transition probabilities to form the diffusion term is computed as:
The diffusion term \({{\varvec{G}}}^x_{\zeta }\) in the SELPAs model is obtained as:
where \({{\varvec{I}}}_{13}\) is the 13 \(\times 13\) identity.
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 lumped, deterministic SEIRS formulations [9, 10] shown below is extended to develop a stochastic malaria model. The human host population is divided into four distinct classes: susceptible (\(S_h\)), exposed (\(E_h\)), infectious (\(I_h\)), and recovered (\(R_h\)). The adult vector population is divided into three classes: susceptible (\(S_v\)), exposed (\(E_v\)), and infectious (\(I_v\)). Here, the aquatic stages of vectors are not considered in SEIRS models. For \(t > t_0\) and \(\forall \zeta \in \mathbb {D}\), the deterministic SEIRS model is given by another set of ODEs that characterize:

The rate of change in susceptible host (\(S_{h,\zeta }\)) as a function of immigration, birth, human infection, recovery from infection, and human mortality:
$$\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)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} \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}$$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.

The rate of change in exposed host (\(E_{h,\zeta }\)) as a function of new host infections, latent period, and human mortality:
$$\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)where \(\nu _h\) is per capita rate of progression of humans from exposed to infectious state.

The rate of change in infected host (\(I_{h,\zeta }\)) as a function of latent period, recovery and human mortality rates:
$$\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)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.

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:
$$\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 \(\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} \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}$$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.

The rate of change in exposed vector (\(E_{v,\zeta }\)) as a function of new vector infection, vector latent period, and vector mortality:
$$\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)where \(\nu _v\) is per capita rate of progression of mosquitoes from the exposed state to the infectious state.

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)
A summary of parameters associated with (10) are shown in Table 3. All parameters described are strictly positive with the exception of the diseaseinduced death rate, \(\delta _h\), which is nonnegative [10]. The model is wellposed in a positively invariant domain:
The development rate of Plasmodium parasites within humans \(\nu _h\), or the intrinsic incubation period, is assumed to be temperature independent.
However, the development rate of Plasmodium within mosquitoes \(\nu _v\), or extrinsic incubation period, is highly dependent on air temperature. Unlike the original models [9, 10], the fitted temperaturedevelopment function for the development rate of parasites within the vector [33,34,35] is also incorporated to the extended model:
where \(\nu _v(T_a) : \mathbb {R} \mapsto \mathbb {R}\) is the progression rate of mosquitoes from exposed to infectious state and \(T_a \le 35\).
Next, the SDE model from the ODE systems described above is derived to incorporate random fluctuations of malaria transmission. It is again assumed that there is variability described by random noise in the transitions between states in (10). Similarly, a discrete (Markov chain) stochastic model is first developed and the expected changes and covariance matrix of changes of the discrete stochastic processes are identified. Then, the continuous stochastic SEIRS model, hereafter referred to as SSEIRS, is derived. Solutions of the discrete Markov chain and continuous SSEIRS models approximately have the same probability distribution as well.
Let \(\mathcal {S}_{i,\zeta }(t)\), \(\mathcal {E}_{i,\zeta }(t)\), \(\mathcal {I}_{i,\zeta }(t)\), and \(\mathcal {R}_{h,\zeta }(t)\) denote continuous random variables for the density of susceptible, exposed, infectious for human (\(i=h\)) and vector (\(i=v\)), and recovered human, respectively (Fig. 3). The SSEIRS model depends on the state variables:
where \({{\varvec{Y}}}_{k}(t) = \{Y_{k,\zeta }(t) \in \mathbb {R} : \zeta \in \mathbb {D} \} \subseteq \mathbb {R}^{m \times n}_{\ge 0}\) for \(k=1,\ldots ,7\) in which \({{\varvec{Y}}}_{k}(t)\) has an associated probability density function \({{\varvec{p}}}_{k}(y, t)\):
Let \(\Delta {\varvec{Y}}(t) = {{\varvec{Y}}}(t+\Delta t)  {{\varvec{Y}}}(t)\). For small \(\Delta t\), there are 15 possible unit changes in \(\Delta Y_{k,\zeta }(t)\) of the discrete stochastic processes associated with different probabilities in each grid cell \(\zeta\) (Table 4). The It\(\hat{o}\) stochastic differential equations for the SSEIRS model are given as:
in which \({{\varvec{F}}}^y : \mathbb {R}^{7\times m \times n} \rightarrow \mathbb {R}^{7\times m \times n}\), \({\varvec{G}}^y : \mathbb {R}^{7\times 15 \times m \times n} \rightarrow \mathbb {R}^{7 \times 15 \times m \times n}\), \({{\varvec{W}}}(t) \subset \mathbb {R}^{15 \times m \times n}\) is matrix of independent Wiener processes, and \({{\varvec{p}}}^y \subseteq \mathbb {R}^{m \times n}\) is the matrix of parameters which are functions of time t and climatic and socioeconomic factors \({{\varvec{C}}^y}(t) \subseteq \mathbb {R}^{m \times n}\). The drift term \({{\varvec{F}}}^y\) is calculated as:
for \(k = 1,\ldots,7\) and \(\zeta \in \mathbb {D}\). Using a similar approach for derivation of the diffusion term in SELPAs model, the form of diffusion term \({{\varvec{G}}}^y\) in SSEIRS is written as:
in which:
Thus, the diffusion term in the SSEIRS model is obtained as follows:
where \({{\varvec{I}}}_{15}\) is the 15 \(\times 15\) identity.
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
SSEIRS model is coupled with the SELPAs formulation through the equality constraint of adult vector population. In essence, SSEIRS represents different states (i.e. susceptible, exposed, infected) of adult vectors in SELPAs. Therefore, the total populations of adult mosquitoes in the two model are the same. The equality constraints are given as:
\(\forall \zeta \in \mathbb {D}\) and \(t > t_0\). For spatial movement of adult vectors among adjacent cells, it is assumed that vector populations are wellmixed or the fraction of malaria classes in the adult vector population remain unchanged during movements:
where \(t_b\) and \(t_a\) represent time before and after vector movement in every modelled time step \(\Delta t\), respectively. Furthermore, birth and mortality rates of vectors in the SSEIRS model are excluded as these processes are already represented in SELPAs model. The proposed coupling approach between SELPAs and SSEIRS models presented allows to simulate stochastically the spatial dynamics of both vector population and malaria transmission over time. Finally, the coupled stochastic latticebased vector dispersal and malaria model (SLIM) can be written as:
in which \({{\varvec{X}}}(t)\) and \({{\varvec{Y}}}(t)\) are random variables described in (3) and (12). Solutions for equations (21) are obtained numerically.
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,38,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.
Female mosquitoes deposit eggs in breeding habitat of various sizes. However, a large fraction of the breeding habitat is at scales that are not detectable by currently available topographic data. As a result, there is always a probability that ponding exists in smallscale topographic depressions inside a particular nonsaturated cell \(\zeta\) that hydrologic modelling cannot capture. To address this scale mismatch issue, we incorporate the fractal structure found in topographic depressions to the estimation of \(\psi ^\mathsf {W}_{\zeta }\). Specifically, the hypothesis is that topographic depressions exist at all sizes on the landscape following the power scaling law [41]. Available topographic data is used to find the scaling relationship of topographic depressions in bounded domain \(\mathbb {D}\) and assume that this relationship remains unchanged at smaller scales for estimation of \(\psi ^\mathsf {W}_{\zeta }\) as below:
in which \(\Theta _{\zeta } \in [\Theta _{min},1]\) is the degree of saturation of the soil surface in cell \(\zeta\) obtained from the Dhara model, \(\Theta _{min}\) is the minimum degree of saturation that is a function of soil properties, and \(\alpha\) is the negative slope of the powerlaw scaling relationship for the exceedance probability of the area of topographic depressions found in the domain under consideration. In this work, \(\alpha\) is identified separately using a topographic depression identification (TDI) algorithm [41].
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
SELPAs model simulations are implemented using similar parameter sets shown in a previous study [28] (Table 2). The model is tested in a simple rectangular domain \(\mathbb {D}_1\) partitioned into uniform and lattice grid \(\{\zeta := (i,j) \in \mathbb {Z}^2_{+} : i \le 5, j \le 5\}\). The size of each grid cell is \(100\text{m} \times 100\text{m}\). In addition, periodic boundary conditions are applied. The model is run with four different initial conditions and sizes of vector population uniformly distributed over \(\mathbb {D}_1\) to analyse the effects of random processes on mosquito population dynamics. As SELPAs is run independently, homogeneous moisture index and human distribution in \(\mathbb {D}_1\) are assumed for simplicity. In addition, to isolate the effects of stochastic noise on the dynamics of vector population, all parameters in SELPAs are assumed to be constants over time. In other words, model parameters’ dependences on hydroclimatic factors are excluded in the simulations. SELPAs simulations are conducted in 800 days using daily time step.
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
Next, the applicability of the SLIM model for largescale simulations of malaria is demonstrated in a rural area in Kilifi county, Kenya. This region has high levels of malnutrition as well as a high incidence rate of Plasmodium falciparum parasites for which the Anopheles gambiae is the main vector [42]. Indeed, the intensity of malaria parasite transmission in Kilifi county is complex, subject to long and shortterm cycles of variation driven by climate, changes in human land use and the efficacy and coverage of interventions that target the parasite and vector [43]. Previous studies showed a decline in malaria transmission during the 1998–2009 period. However, there was a steady and marked increase of malaria transmission from 2009 to 2014 [43, 44]. Here, the primary objective is to present the capabilities and advantages of SLIM model for capturing the spatial and temporal variations of factors that drive malaria transmission. Therefore, the model is not validated for the case study. Model validation using observed malaria prevalence data and the impacts of climate change on malaria in coastal Kenya is presented in another work [45].
The domain of simulation is approximately 440 km\(^2\) (22 km north to south and 20 km east to west) with medium to high percentage cover of vegetation (Fig. 5). It is assumed that natural water on the ground surface is the main habitat of Anopheles mosquitoes. Topographic data at 30 m \(\times\) 30 m resolution from Advanced Spaceborne Thermal Emission and Reflection radiometer (ASTER) global digital elevation model is used for modelling surface runoff and belowground soil moisture dynamics (Fig. 5a). Global reanalysis meteorological data at 3 h interval by European Centre for Mediumrange Weather Forecasts (ECMWF) from 2005 to 2014 is used to drive both Dhara and SLIM models. Human population census data at 100 m spatial resolution is obtained from the population maps for low income nations [46]. Topographic depressions in the study area are found using TDI algorithm for estimating a moisture index as presented in "Estimating moisture index" (see [41]). The distributions of topographic depressions in the study domain are shown in Fig. 5b. The spatial heterogeneities of vector habitat and human hosts are the main drivers for mosquitoes movement and spread of malaria parasites. Model parameter sets similar to those used in previous studies [10, 28] are used for model simulations.
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.
References
Miller LH, Baruch DI, Marsh K, Doumbo OK. The pathogenic basis of malaria. Nature. 2002;415:673–9.
Smith DL, McKenzie EF. Statics and dynamics of malaria infection in Anopheles mosquitoes. Malar J. 2004;3:1–14.
Anderson RM. The population dynamics of infectious diseases: theory and applications. Population and community biology series. London: Chapman & Hall Ltd.; 1982.
Anderson RM, May RM. Infectious diseases of humans: dynamics and control. Dynamics and control. Oxford: Oxford University Press; 1992.
Paaijmans KP, Thomas MB. Health: wealth versus warming. Nat Clim Change. 2011;1:349–50.
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.
Ross R. The prevention of malaria. 2nd ed. Dutton; 1910.
MacDonald G. The Epidemiology and Control of Malaria. Oxford Medical Publications. Oxford, UK: Oxford University Press; 1957.
Ngwa GA, Shu WS. A mathematical model for endemic malaria with variable human and mosquito populations. Math Comput Model. 2000;32:747–63.
Chitnis N, Cushing J, Hyman J. Bifurcation analysis of a mathematical model for malaria transmission. SIAM J Appl Math. 2006;67:24–45.
Yang HM. Malaria transmission model for different levels of acquired immunity and temperaturedependent parameters (vector). Rev Saude Publica. 2000;34:223–31.
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.
Parham PE, Michael E. Modeling the effects of weather and climate change on malaria transmission. Environ Health Perspect. 2010;118:620–6.
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.
Bomblies A, Duchemin JB, Eltahir EAB. Hydrology of malaria: model development and application to a Sahelian village. Water Resour Res. 2008;44:W12445.
Gu W, Novak RJ. Agentbased modelling of mosquito foraging behaviour for malaria control. Trans R Soc Trop Med Hyg. 2009;103:1105–12.
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.
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.
Mandal S, Sarkar R, Sinha S. Mathematical models of malaria—a review. Malar J. 2011;10:202.
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.
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.
Keeling MJ, Rohani P. Modeling infectious diseases in humans and animals. Princeton: Princeton University Press; 2008.
Azaele S, Maritan A, Bertuzzo E, RodriguezIturbe I, Rinaldo A. Stochastic dynamics of cholera epidemics. Phys Rev E. 2010;81:051901.
Herwaarden OA, Grasman J. Stochastic epidemics: major outbreaks and the duration of the endemic period. J Math Biol. 1995;33:581–601.
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.
Britton T. Stochastic epidemic models: a survey. Math Biosci. 2010;225:24–35.
Krstic M. The effect of stochastic perturbation on a nonlinear delay malaria epidemic model. Math Comput Simul. 2011;82:558–69.
Lutambi AM, Penny MA, Smith T, Chitnis N. Mathematical modelling of mosquito dispersal in a heterogeneous environment. Math Biosci. 2013;241:198–216.
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.
Allen E. Modeling with Itô Stochastic differential equations. Mathematical modelling: theory and applications. Heidelberg, Germany: Springer Berlin Heidelberg; 2007.
Allen LJS. An introduction to stochastic processes with applications to biology. 2nd ed. Florida: CRC Press; 2010.
Doolan DL, Dobaño C, Baird JK. Acquired immunity to malaria. Clin Microbiol Rev. 2009;22:13–36.
Detinova TS. Agegrouping methods in Diptera of medical importance with special reference to some vectors of malaria. WHO Monograph series. 1962;47:13–191.
Briere JF, Pracros P, Le Roux AY, Pierre JS. A novel rate model of temperaturedependent development for arthropods. Environ Entomol. 1999;28:22–9.
Paaijmans KP, Read AF, Thomas MB. Understanding the link between malaria risk and climate. Proc Natl Acad Sci USA. 2009;106:13844–9.
Le PVV, Kumar P. Interaction between ecohydrologic dynamics and microtopographic variability under climate change. Water Resour Res. 2017;53:8383–403.
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.
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.
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.
Le PVV, Kumar P, Valocchi AJ, Dang HV. GPUbased highperformance computing for integrated surfacesubsurface flow modeling. Environ Modell Softw. 2015;73:1–13.
Le PVV, Kumar P. Power law scaling of topographic depressions and their hydrologic connectivity. Geophys Res Lett. 2014;41:1553–9.
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.
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.
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.
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.
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.
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.
Author information
Authors and Affiliations
Corresponding author
Rights and permissions
Open Access This 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.
About this article
Cite this article
Le, P.V.V., Kumar, P. & Ruiz, M.O. Stochastic latticebased modelling of malaria dynamics. Malar J 17, 250 (2018). https://doi.org/10.1186/s129360182397z
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s129360182397z
Keywords
 Malaria
 Climate change
 Metapopulation
 Stochastic
 Ecohydrology