Stochastic lattice-based modelling of malaria dynamics.

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, lumped-population models of malaria prevalence have been insufficient for predicting the complex responses of malaria to environmental changes. METHODS AND
RESULTS: A stochastic lattice-based model that couples a mosquito dispersal and a susceptible-exposed-infected-recovered epidemics model was developed for predicting the dynamics of malaria in heterogeneous environments. The It[Formula: see text] 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 lattice-based integrated malaria model has been developed. The applicability of the model for capturing the climate-driven hydrologic factors and demographic variability on malaria transmission has been demonstrated.

Background

Malaria is a vector-borne 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], socio-economic developments [11], temporary immunity [12], and weather effects [13]. In addition, agent-based and meta-population 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 host-parasite-vector interactions. Systematic reviews of mathematical modelling of malaria [19] and other mosquito-borne diseases [20, 21] indicate the need for models to address the complexities is the host-vector-parasite interactions and to incorporate heterogeneous environments.

Deterministic models with susceptible-exposed-infected-recovered (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 agent-based models [15, 17, 18] are inefficient for large-scale simulations with very large number of vector individuals involved, while lumped-population models [26, 27] bypass the spatial dynamics of the diseases.

In this paper, a stochastic lattice-based 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 SEIRS-type 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 climate-driven hydrologic and ecologic processes as factors that determine mosquito population and malaria transmission dynamics. A meta-population approach is incorporated to describe the movements of vectors among discrete geographic sub-domains. The It 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 time-continuous, space-discrete 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 be a bounded domain in and let be integers. is partitioned into uniform and rectangular lattice grid indexed as . In each grid cell , 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.

Fig. 1

Schematic of the SLIM model that couples a vector dispersal model with a malaria dynamics formulation. a Ponding and moisture index obtained from an ecohydrologic model (Dhara) provide the habitat for gravid female Anopheles to deposit eggs. The dispersal of host-seeking mosquitoes is based on a host searching index calculated as a function of human density in each grid cell. The rates of movement of mosquitoes among adjacent grid cells are diffusive. In each cell, the changes of sub-populations (x-axes) in vector dispersal and malaria dynamics models over time period are described by transition probabilities. b Sub-population of compartmental malaria dynamics model in each grid cell. c Sub-population of vector dynamics model in each grid cell. The vector dispersal and malaria dynamics models share the adult vector population which affects both the aquatic density and malaria transmission in human hosts

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: host-seeking (), resting (), and oviposition site searching (). 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. and are further affected by the movement of mosquitoes in space. For and , 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.

Table 1

Description and values of parameters of the ELPA model.

Modified from [28]

Name	Description	Unit	Range
b	Integer number of female eggs laid per oviposition	−	50–300
ψζW	50% of the eggs are assumed to hatch into female mosquitoes parameter represent water availability in cell ζ	−	0.0–1.0
ψζH	Binary parameter represent human presence in cell ζ	−	0–1
ρ_E	Egg hatching rate into larvae	day^-1	0.33–1.0
ρ_L	Rate at which larvae develop into pupae	day^-1	0.08–0.17
ρ_P	Rate at which pupae develop into adult/emergence rate	day^-1	0.33–1.0
μ_E	Egg mortality rate	day^-1	0.32–0.80
μ_L_1	Natural mortality rate of larvae	day^-1	0.30–0.58
μ_L_2	Density-dependent mortality rate of larvae	day^-1mosq^-1	0.0–1.0
μ_P	Pupae mortality rate	day^-1	0.22–0.52
ρ_A_h	Rate at which host-seeking mosquitoes enter the resting state	day^-1	0.322–0.598
ρ_A_r	Rate at which resting mosquitoes enter oviposition searching state	day^-1	0.30–0.56
ρ_A_o	Oviposition rate	day^-1	3.0–4.0
μ_A_h	Mortality rate of mosquitoes of searching for hosts	day^-1	0.125–0.233
μ_A_r	Mortality rate of resting mosquitoes	day^-1	0.0034–0.01
μ_A_o	Mortality rate of mosquitoes searching for oviposition sites	day^-1	0.41–0.56

The above model (1) is well-posed 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 , where is the mean air temperature (K) over the time interval ; and is the development rate per unit time at . Details for the temperature dependencies of egg, larval, and pupal populations are presented in Depinay et al. [29].

The rate of change in egg population () 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; is the rate at which eggs are oviposited by gravid mosquitoes; is egg mortality rate; and is the hatching rate into larvae. The term on the right hand side represents water availability in a particular cell and is discussed later.

The rate of change in larvae population () as a function of egg population, larval mortality and maturation into pupae: where  is the progression rate from larvae to pupae; and represent natural and density-dependent death rates of larvae, respectively.

The rate of change in pupae population () as a function of larval maturation, pupal mortality, and emergence into adults: where is the mortality rate of pupae, and represent the rate of emergence from pupae into adults.

The rate of change in population of host-seeking adults () as a function of pupal emergence, oviposition, mortality, and blood feeding rates: where is the death rate of host-seeking adults; and is the rate at which they enter a resting state after blood feeding. Host-seeking adults can spread to adjacent cells for searching human host in which represents their movement rate from cell to modelled as a decreasing exponential function of human population and average area of cell and  (see [28]); denotes the neighbors of the cell under consideration. The last two terms in (1d) represent the movements of vectors from all neighbors into cell and vice versa, respectively. After laying eggs, gravid mosquitoes return to the host-seeking state for subsequent blood feeding.

The rate of change in population of resting adults () as a function of blood feeding, mortality, and protein digestion rates: where is the death rate of resting adults; and 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 () as a function of emergence, oviposition, mortality, and digestion rates: where is the death rate of gravid female mosquitoes. Vectors in ovipositing state also spread in space to find water for oviposition in which represents their movement rate from cell to modelled as a decreasing exponential function of surface soil moisture and average area of cell and  [28].

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 , 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 S-ELPAs, 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 S-ELPAs models approximately have the same probability distribution.

Let , , , , , denote continuous random variables for the density of eggs, larvae, pupae, host-seeking adults, resting adults, and oviposition site searching adults in a grid cell , respectively (see Fig. 2). The S-ELPAs model depends on the state variables:where for in which has an associated probability density function :Let . For small , there are 13 possible unit changes in of the discrete stochastic model associated with different probabilities in each cell (see Table 2). The It stochastic differential equations for the S-ELPAs model are given as follow [30]:where ; ; is the matrix of independent Wiener processes; and is the matrix of parameters that are functions of time t and climatic, anthropogenic, and entomological factors .

Fig. 2

Schematic representation of Anopheles life and feeding cycles. The first three stages are aquatic. The last three stages are adult, which are able to carry the pathogens

(Modified from Lutambi et al. [28])

Table 2

Probabilities associated with changes in ELPA model

l	Change, ΔXk,ζl(t)	Probability, pk,ζl(t)	Description
1	[1,0,0,0,0,0]^T	bψζWρ_A_oA_o,ζΔt	A new egg E is deposited by A_o
2	[-1,0,0,0,0,0]^T	μ_EE_ζΔt	An egg E dies
3	[-1,1,0,0,0,0]^T	ρ_EE_ζΔt	An egg E hatches into a larva L
4	[0,-1,0,0,0,0]^T	(μ_L_1+μ_L_2L_ζ)L_ζΔt	A larva L dies
5	[0,-1,1,0,0,0]^T	ρ_LL_ζΔt	A larva L develops into a pupa P
6	[0,0,-1,0,0,0]^T	μ_PP_ζΔt	A pupa P dies
7	[0,0,-1,1,0,0]^T	ρ_PP_ζΔt	A pupa P develops into a host-seeking adult A_h
8	[0,0,0,1,0,-1]^T	ψζWρ_A_oA_o,ζΔt	An oviposition adult A_o enters host-seeking state
9	[0,0,0,-1,0,0]^T	μ_A_hA_h,ζΔt	A host-seeking adult A_h dies
10	[0,0,0,-1,1,0]^T	ψζHρ_A_hA_h,ζΔt	A host-seeking adult A_h enters resting state
11	[0,0,0,0,-1,0]^T	μ_A_rA_r,ζΔt	A resting adult A_r dies
12	[0,0,0,0,-1,1]^T	ρ_A_rA_r,ζΔt	A resting adult A_r enters oviposition searching state
13	[0,0,0,0,0,-1]^T	μ_A_oA_o,ζΔt	An oviposition searching adult A_o dies

The drift term in (5), to order , is the expectation of all possible changes in the discrete stochastic model computed as:for and .

Additional Wiener processes are included into the stochastic systems to simplify the derivation of the diffusion term  [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 in the S-ELPAs model is obtained as:where is the 13  identity.

The S-ELPAs 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. S-ELPAs is coupled with a stochastic SEIRS formulation presented below for malaria transmission dynamics.

Malaria transmission model

The S-ELPAs 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 climate-driven 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 (), exposed (), infectious (), and recovered (). The adult vector population is divided into three classes: susceptible (), exposed (), and infectious (). Here, the aquatic stages of vectors are not considered in SEIRS models. For and , the deterministic SEIRS model is given by another set of ODEs that characterize:

The rate of change in susceptible host () as a function of immigration, birth, human infection, recovery from infection, and human mortality: where is immigration rate; and is per capita birth rate of humans; 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]. is total population size for humans in each cell ; is the human per capita death rate; and is the infection rate from mosquitoes to humans defined as: in which is total population of mosquitoes in cell ; represents the number of times one mosquito attempt to bite humans per unit time; is the maximum number of mosquito bites a human can have per unit time; and 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 () as a function of new host infections, latent period, and human mortality: where is per capita rate of progression of humans from exposed to infectious state.

The rate of change in infected host () as a function of latent period, recovery and human mortality rates: where represents per capita recovery rate for humans from infectious to recovered states; and represents per capita disease-induced death rate for humans.

The rate of change in recovered host () as a function of recovery, immunity loss, and human mortality:

The rate of change in susceptible vector () as a function of reproduction, vector infection, and vector mortality: where represents per capita birth rate of the vectors; is the per capita death rate for vectors in each cell ; and is the infection rate from humans to mosquitoes defined as: where and 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 () as a function of new vector infection, vector latent period, and vector mortality: where is per capita rate of progression of mosquitoes from the exposed state to the infectious state.

The rate of change in infected vector () as a function of latent period and mortality:

A summary of parameters associated with (10) are shown in Table 3. All parameters described are strictly positive with the exception of the disease-induced death rate, , which is non-negative [10]. The model is well-posed in a positively invariant domain:The development rate of Plasmodium parasites within humans , or the intrinsic incubation period, is assumed to be temperature independent.

Table 3

The parameters for SEIRS malaria model.

From [10]

Name	Description	Unita
Λ_h	Immigration rate of humans	H × T^-1
ψ_h	Per capita birth rate of humans	T^-1
ψ_v	Per capita birth rate of mosquitoes	T^-1
σ_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
σ_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
β_hv	Probability of transmission of infection from an infectious mosquito to a susceptible human, given that a contact between the two occurs	−
β_vh	Probability of transmission of infection from an infectious human to a susceptible mosquito, given that a contact between the two occurs	−
β~_hv	Probability of transmission of infection from a recovered (asymptomatic carrier) human to a susceptible mosquito, given that a contact between the two occurs	−
ν_h	Per capita rate of progression of humans from the exposed state to the infectious state. 1/ν_h is the average duration of the latent period	T^-1
ν_v	Per capita rate of progression of mosquitoes from the exposed state to the infectious state. 1/ν_v is the average duration of the latent period	T^-1
γ_h	Per capita recovery rate for humans from the infectious state to the recovered state. 1/γ_h is the average duration of the infectious period	T^-1
δ_h	Per capita disease-induced death rate for humans	T^-1
ρ_h	Per capita rate of loss of acquired temporary immunity for humans. 1/ρ_h is the average duration of the immune period	T^-1
μ_1h	Density-independent part of the death (and emigration) rate for humans	T^-1
μ_2h	Density-dependent part of the death (and emigration) rate for humans	H × T^-1
μ_1h	Density-independent part of the death (and emigration) rate for mosquitoes	T^-1
μ_2h	Density-dependent part of the death (and emigration) rate for mosquitoes	M × T^-1

aIn the Unit, H represents humans, M represents mosquitoes, and T represents time

However, the development rate of Plasmodium within mosquitoes , or extrinsic incubation period, is highly dependent on air temperature. Unlike the original models [9, 10], the fitted temperature-development function for the development rate of parasites within the vector [33-35] is also incorporated to the extended model:where is the progression rate of mosquitoes from exposed to infectious state and .

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 S-SEIRS, is derived. Solutions of the discrete Markov chain and continuous S-SEIRS models approximately have the same probability distribution as well.

Let , , , and denote continuous random variables for the density of susceptible, exposed, infectious for human () and vector (), and recovered human, respectively (Fig. 3). The S-SEIRS model depends on the state variables:where for in which has an associated probability density function :

Fig. 3

Schematic representation of malaria transmission. The model divides the human population into four classes: susceptible, ; exposed, ; infectious, ; and recovered (immune), . Vector population is divided into three classes: susceptible, ; exposed, ; and infectious, . Both species follow a logistic population model, with humans having additional immigration and disease-induced death. Birth, death, and migration into and out of the population are not shown in the figure

(Adapted from Chitnis et al. [10])

Let . For small , there are 15 possible unit changes in of the discrete stochastic processes associated with different probabilities in each grid cell (Table 4). The It stochastic differential equations for the S-SEIRS model are given as:in which , , is matrix of independent Wiener processes, and is the matrix of parameters which are functions of time t and climatic and socio-economic factors . The drift term is calculated as:for and . Using a similar approach for derivation of the diffusion term in S-ELPAs model, the form of diffusion term in S-SEIRS is written as:in which:Thus, the diffusion term in the S-SEIRS model is obtained as follows:where is the 15  identity.

Table 4

Probabilities associated with changes in SEIRS model

l	Change, ΔYk,ζl(t)	Probability, pk,ζl(t)	Description
1	[1,0,0,0,0,0,0]^T	(Λ_h+ψ_hN_h,ζ)Δt	A new host enters the human susceptible class
2	[1,0,0,-1,0,0,0]^T	ρ_hR_h,ζΔt	A recovered host becomes susceptible again
3	[-1,1,0,0,0,0,0]^T	σ_vσ_hβ_hvI_v,ζS_h,ζσ_vN_v,ζ+σ_hN_h,ζΔt	A susceptible host enters exposed state
4	[-1,0,0,0,0,0,0]^T	(μ_1h+μ_2hN_h,ζ)S_h,ζΔt	A susceptible host dies
5	[0,-1,1,0,0,0,0]^T	ν_hE_h,ζΔt	An exposed host enters infectious state
6	[0,-1,0,0,0,0,0]^T	(μ_1h+μ_2hN_h,ζ)E_h,ζΔt	An exposed host dies
7	[0,0,-1,1,0,0,0]^T	γ_hI_h,ζΔt	An infectious host enters recovered state
8	[0,0,-1,0,0,0,0]^T	(μ_1h+μ_2hN_h,ζ+δ_h)I_h,ζΔt	An infectious host dies
9	[0,0,0,-1,0,0,0]^T	(μ_1h+μ_2hN_h,ζ)R_h,ζΔt	A recovered host dies
10	[0,0,0,0,1,0,0]^T	ψ_vN_v,ζΔt	A new mosquito enters the vector susceptible class
11	[0,0,0,0,-1,1,0]^T	σ_vσ_hβ_hvI_v,ζS_h,ζβ_vhI_h,ζ+β~_vhR_h,ζΔt	A susceptible vector enters exposed state
12	[0,0,0,0,-1,0,0]^T	(μ_1v+μ_2vN_v,ζ)S_v,ζΔt	A susceptible vector dies
13	[0,0,0,0,0,-1,1]^T	ν_vE_v,ζΔt	An exposed vector enters infectious state
14	[0,0,0,0,0,-1,0]^T	(μ_1v+μ_2vN_v,ζ)E_v,ζΔt	An exposed vector dies
15	[0,0,0,0,0,0,-1]^T	(μ_1v+μ_2vN_v,ζ)I_v,ζΔt	An infectious vector dies

The S-SEIRS model incorporates environmental perturbation and stochastic interactions among sub-populations 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 mosquito-borne diseases.

Model couplings

S-SEIRS model is coupled with the S-ELPAs formulation through the equality constraint of adult vector population. In essence, S-SEIRS represents different states (i.e. susceptible, exposed, infected) of adult vectors in S-ELPAs. Therefore, the total populations of adult mosquitoes in the two model are the same. The equality constraints are given as: and . For spatial movement of adult vectors among adjacent cells, it is assumed that vector populations are well-mixed or the fraction of malaria classes in the adult vector population remain unchanged during movements: where and represent time before and after vector movement in every modelled time step , respectively. Furthermore, birth and mortality rates of vectors in the S-SEIRS model are excluded as these processes are already represented in S-ELPAs model. The proposed coupling approach between S-ELPAs and S-SEIRS models presented allows to simulate stochastically the spatial dynamics of both vector population and malaria transmission over time. Finally, the coupled stochastic lattice-based vector dispersal and malaria model (SLIM) can be written as: in which and are random variables described in (3) and (12). Solutions for equations (21) are obtained numerically.

Estimating moisture index

The moisture index shown in (1a) represents water availability in a particular cell 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 physically-based surface-subsurface flow model coupler (GCS-flow, see [40]) designed for capturing moisture transport on the land surface and in the below-ground systems. It incorporates vegetation acclimation to elevated CO 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 S-ELPAs and S-SEIRS) model is linked with Dhara for incorporating climate-driven 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 small-scale topographic depressions inside a particular non-saturated cell  that hydrologic modelling cannot capture. To address this scale mismatch issue, we incorporate the fractal structure found in topographic depressions to the estimation of . 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 and assume that this relationship remains unchanged at smaller scales for estimation of as below:in which is the degree of saturation of the soil surface in cell obtained from the Dhara model, is the minimum degree of saturation that is a function of soil properties, and is the negative slope of the power-law scaling relationship for the exceedance probability of the area of topographic depressions found in the domain under consideration. In this work, 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 (S-ELPAs and S-SEIRS) 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.

S-ELPAs model

S-ELPAs 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 partitioned into uniform and lattice grid . The size of each grid cell is . In addition, periodic boundary conditions are applied. The model is run with four different initial conditions and sizes of vector population uniformly distributed over to analyse the effects of random processes on mosquito population dynamics. As S-ELPAs is run independently, homogeneous moisture index and human distribution in are assumed for simplicity. In addition, to isolate the effects of stochastic noise on the dynamics of vector population, all parameters in S-ELPAs are assumed to be constants over time. In other words, model parameters’ dependences on hydro-climatic factors are excluded in the simulations. S-ELPAs simulations are conducted in 800 days using daily time step.

Figure 4a shows the variations in log-scale of adult mosquito density averaged over under the effects of random processes. It can be seen that S-ELPAs 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 meta-population and lattice-based models is thus important to study the dynamics of vector density in large areas, where populations at various sizes are interconnected.

Fig. 4

The dynamics of S-ELPAs and S-SEIRS models with white noise. a Simulation of total adult mosquitoes at different sizes of initial population. The graph shows how the oscillatory behavior becomes disrupted by noise in smaller populations, whereas large populations conform close to the equilibria. b Comparison of the malaria infected cases in humans between deterministic and stochastic simulations. The red curve shows the mean, and the gray shaded region shows the range for simulations of stochastic SEIRS model. The blue curve is from the original deterministic SEIRS model. While deterministic simulation tends to an endemic equilibrium, stochastic simulations show possible extinctions of the disease, as expected. The agreement between deterministic and mean stochastic simulations implies that a small fraction of the stochastic trajectory go to extinction in the simulations

S-SEIRS model

A large number (100) of independent S-SEIRS simulations are conducted and compared with the deterministic SEIRS model to investigate the modelled dynamics of malaria transmission in noise-dominated 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 hydro-climatic factors are also excluded in both S-SEIRS and SEIRS simulations. As the spatial movement of mosquitoes is not considered in S-SEIRS 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 () obtained from stochastic S-SEIRS (shaded area) and deterministic SEIRS (blue solid line) models are shown in Fig. 4b. Unlike the deterministic approach, S-SEIRS simulations provide a range of possibilities for malaria transmission given the same initial conditions and model parameters. S-SEIRS 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 S-SEIRS 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 non-isolated 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 large-scale 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 short-term 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 (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 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 Medium-range 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.

Fig. 5

Domain of simulations in the case study at Kilifi county, Kenya. a Variation in topographic elevation. b Map of topographic depression (red polygons) identified from ASTER digital elevation model. The gray background represents hillshaded topography

Figure  6a–d presents the variation of total aquatic and adult phases of mosquito populations obtained from the S-ELPA 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 S-ELPAs with observed monthly air temperature () and rainfall (), 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 host-seeking mosquitoes () in adult stage is high, consisting of of the total adult population. The sub-population of oviposition site searching mosquitoes () are usually times larger than the resting mosquitoes (). 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.

Fig. 6

Illustration of variation of mosquitoes in different phases of their life cycle and malaria as predicted by SLIM model in response to meteorological driver for the study domain shown in Fig.  5. a Daily precipitation. b Mean daily air temperature. c Population dynamics of mosquitoes in aquatic phase; E represents egg population, L represents larvae population, and P represents pupae population, respectively. d Population dynamics of mosquitoes in adult phase; represents host-seeking mosquitoes, represents resting mosquitoes, and represents oviposition searching mosquitoes, respectively. e Dynamics of malaria within human hosts; represents exposed cases and represents infected cases, respectively. f Dynamics of malaria within human hosts; is susceptible vector, and represents exposed vector, is infected vector, and is the total vector, respectively.

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 (, ) with air temperature () and rainfall (), 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 () are usually found after the rainy seasons start and air temperature was high. The peaks of are also followed by the largest values of infected human cases () 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 lattice-based integrated malaria (SLIM) model that consists of a vector dispersal (S-ELPAs) and a malaria dynamics (S-SEIRS) 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 malaria-environment 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 climate-driven 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.