Stochastic modeling of a mosquito-borne disease

*Correspondence: pwitbooi@uwc.ac.za 1Department of Mathematics and Applied Mathematics, University of the Western Cape, Private Bag X17, Bellville, 7535, South Africa Full list of author information is available at the end of the article Abstract We present and analyze a stochastic differential equation (SDE) model for the population dynamics of a mosquito-borne infectious disease. We prove the solutions to be almost surely positive and global. We introduce a numerical invariantR of the model withR < 1 being a condition guaranteeing the almost sure stability of the disease-free equilibrium. We show that stochastic perturbations enhance the stability of the disease-free equilibrium of the underlying deterministic model. We illustrate the main stability theorem through simulations and show how to obtain interval estimates when making forward projections. We consulted a wide range of literature to find relevant numerical parameter values.


Introduction
A variety of mosquito-borne infectious diseases are the cause of millions of illnesses and deaths. In particular, the malaria disease is responsible for millions of fatalities in Africa each year and is a serious burden of the disease worldwide. The report [27] by Southern African Development Community (SADC) gives a summary of malaria statistics in Southern Africa. A map of the region clearly shows the intensity of the malaria problem in different areas. In Mozambique, for instance, malaria is the leading cause of death [12], whereas the biggest part of South Africa is malaria-free [27,Malaria Map]. Mathematical modeling is useful in the planning of interventions to curb the spread of malaria and other diseases. Indeed, different models have been proposed for different situations. A popular type of models is the compartmental model in terms of ODEs. Recent models of this type include, for instance, [4,22,23,25].
Malaria prevalence numbers have been proved to be influenced by climatic factors; see, for example, [1,2,6,14,28]. Malaria population dynamics is also correlated with other variables such as altitude and topography, land use, land cover, human behavior, and living conditions. Some regions have partial protection against malaria through indoor residual spraying or bednets [5]. In some regions, people may use traditional plant remedies that are effective against malaria; see, for example, [8]. Consequently, modeling malaria pop-ulation dynamics can become quite complex, [7]. Over time, many different approaches and methodologies were introduced into compartmental modeling of disease dynamics, and these innovations are able to deal with different complexities we experience in real life. So, for instance, we have seen the use of networks, multigroups, age-structured models, and stochastic differential equation (SDE) models [13,15,17,[31][32][33]. In a region where the link between malaria prevalence and the relevant factors is not well understood, it may be wise to introduce some randomness into an ODE compartmental model to compensate for uncertainty. In this way, we obtain an SDE model. For further motivation for SDE models, we refer the reader to Higham [11]. Such SDE models have already been proposed for various diseases, for instance, HIV [24], TB [20], and vector-borne diseases [15]; see also [16,18,21,37]. In this paper, we propose an SDE model for the population dynamics of a disease such as malaria or dengue fever. The underlying deterministic model of the present paper is the same as that in [31] and is also similar to that of [17]. However, the stochastic perturbation of the present paper differs from those in [31] and [17]. In the current paper the stochasticity is similar to perturbation of the force of infection, whereas in [17] and [31] the mortality rates are perturbed. Also, the methodologies are quite different from both mathematical analysis and simulation sides. In particular, we observe that the stochastic perturbation in our model enhances the stability of the disease-free equilibrium, similarly as in [10,18,34], and a few other cases. Another feature of the model, simply due to the stochastic perturbations, is that when we make a future projection, we are able to specify a confidence interval for the estimate.
The remainder of this paper is organized as follows. In Sect. 2, we describe the stochastic model. In Sect. 3, we show the existence and uniqueness of a global positive solution of the model. In Sect. 4, we investigate the asymptotic behavior of solutions to the stochastic model around the trivial equilibrium point. In Sect. 5, we run simulations to make future projections of the state of the population relatively to the disease and illustrate the longtime behavior of the model. The parameter values were obtained from a wide range of literature and are generally applicable to Southern Africa. Finally, in Sect. 6, we provide a few concluding remarks and suggestions for future research.

The stochastic model
In this paper, we let (Ω, F, {F t } t≥0 , P) be a filtered complete probability space with the filtration satisfying the usual conditions (i.e., the filtration is right-continuous, and F 0 contains all events of zero probability). We consider a pair of independent Wiener processes W (t) and Z(t) on this probability space. We use the notation R n We consider a stochastic disease model based on the deterministic model in [30]. As in [30], the host population at time t ≥ 0 is of size N(t) and is subdivided into three compartments, and the vector disease population of size M(t) is subdivided into two compartments or classes. We shall ambiguously use the same symbol for a population class and its magnitude. The first compartment in the human population consists of the individuals susceptible to be infected with the pathogen. We denote this class by S(t). The second class, denoted by I(t), consists of all the individuals infected with the pathogen. The third class consists of all the human individuals recovered from the infection and having temporary immunity against the pathogen. This class of individuals is denoted by R(t). The vector compartment is subdivided into two classes, susceptible V and infective J classes.
Due to the short life-span of vectors, we assume that vectors do not recover from the infection and their infective period ends with their death [30]. For the total N(t) of the human population and the total M(t) of the mosquito population, we have the identities The human birth rate is denoted by A. There is no vertical transmission, and all the newborns are susceptible. The per capita death rate (excluding death due to malaria) is assumed to be the same constant μ for all human compartments, and the rate of mortalities due to malaria is denoted by δ. The mosquito population has B and θ as the natural birth rate and per capita mortality rate, respectively. The vectors bite humans at rate a. The fraction of the bites that successfully infect humans is b, whereas c represents the fraction of bites that infect vectors when they bite infected humans. The rate of infection of the human host in class S by infected vectors in J depends on the total number of humans available per infected vector. The per capita rate of transfer from the I-class to the R-class is k. The human immune individuals lose their immunity at a (per capita) rate h.
There are SDE epidemic models in the literature (e.g., [10]) in which the transmission parameter is stochastically perturbed. In such cases, transmission is usually from susceptibles S to the infectious class I, with no latent class and without vector. In particular, it involves a complementary pair of perturbations, which results in the total population size to be a deterministic function of time. Also, each perturbation in such a pair is proportional to the product SI and can be viewed as perturbations on the class sizes S and I. In the current model, we introduce such a complementary pair of perturbations, each proportional to SI, on the class sizes S and I. We place a similar complementary pair of perturbations on the class sizes V and J.
After introducing stochastic perturbations, with nonnegative constants σ and ζ , we obtain the following system of SDEs: Note that stochastically perturbing the transmission rate from the S-class to the I-class would amount to a complementary pair of perturbations proportional to SJ rather than to SI. The difficulty with perturbations proportional to SJ is that the standard proof of positivity of solutions does not apply. Nevertheless, in this regard, not all is lost, since in general there are lagged correlations between J and I due to the incubation periods of the pathogen in the bodies of the host and vector.
We further find it convenient to use the short-hand notation, introducing the symbol A solution of this system over a time interval D is the set of points The stochastic model has a unique disease-free equilibrium point In the particular case σ = ζ = 0, we refer to system (1) as the underlying deterministic model. The basic reproduction number, of the underlying deterministic model has been calculated in [30] and is an indicator of local asymptotic stability of the disease-free equilibrium. In this paper, we are interested in the global asymptotic stability of the disease-free equilibrium of the SDE model. To this end, in Sect. 4, we introduce an invariant R, which serves as an indicator of global asymptotic stability of the disease-free equilibrium for the stochastic model.

Existence and uniqueness of a positive solution
In this section, we show that the solution of system (1) is almost surely (a.s.) global and positive. If the coefficients of the equations satisfy the local Lipschitz condition and the linear growth condition, then the system of equations has a unique global solution (i.e., no explosion in finite time) for any given initial value (see, e.g., [21]). However, the coefficients of (1) do not satisfy the linear growth condition. In this section, using Lyapunov analysis (as mentioned in [21] and popularly applied in the literature), we show that the solution of (1) is (a.s.) global and positive (Theorem 3.2).

Proposition 3.1 Given any t
Proof Given any local solution with This proves (a). The proof of (b) is similar. Theorem 3.2 For any given initial value X(0) ∈ R 5 + , there exists a unique positive solution X(t) of (1) for t ≥ 0 such that the solution remains in R 5 + with probability 1, that is, X(t) ∈ R 5 + for all t ≥ 0 almost surely.
Sketch of proof Motivated by the work of Mao et al. [21] and similar proofs in [16,17,36], and elsewhere, we present a sketch of a proof. Note that the coefficients of system (1) are locally Lipschitz continuous. Thus there exists a unique local solution on t ∈ [0, τ en ), where τ en is the explosion time. We need to show that this solution is global almost surely, that is, τ en = ∞ (a.s.). This is commonly done using the Lyapunov method. Essentially, we need a boundedness result, holding for t < τ en , which we obtain further. We denote by L the infinitesimal generator of the stochastic process in Eq. (1). Let us define the stochastic process Y 1 as .
Note that each of the five terms are nonnegative, when S, I, R, V , and J are strictly positive. Applying Itô's formula, we obtain where LY 1 (X(t)) is calculated as (recall that μ 1 = μ + k + δ). From the latter expression we obtain the following inequality by removing some of the negative terms on the right-hand side: The rest of the proof is similar to that in [16,21,36], or [17], and we skip the details.
A stability theorem for the stochastic model is formulated in terms of an indicator R(z), which is similar to the basic reproduction number R 0 of the UDM. In fact, R(z) is of a form similar to R 0 , and as the intensities of the perturbations tend to zero, then R(z) tends to R 0 . The difference between R(z) and R 0 is determined by the term E 0 in Eq. (4). Somehow we must find a lower bound for E 0 . Hence in Item 3.3 further, we briefly derive such a lower bound. Similar techniques appear in [20,24], and [35].

Remark 3.3 (A function H(x))
In the stability analysis, we encounter a function for which we require a lower bound. We introduce the function We find that, indeed, H has such a minimum value at x = x * with Noting that S ≤ A/μ and V ≤ B/θ for S and V as in the model, we obtain the inequality for all x ∈ [0, 1].
We find it useful to write this inequality in the following form:

Stability of the disease-free equilibrium
For stochastic systems, there are many different versions of the concept of stability and very sophisticated methods of stability analysis; see, for example, [26]. In this paper, we focus on almost sure exponential stability, which is also studied in [26].
Remark 4.1 Let us fix some nonnegative constants q, p 1 , p 2 , and p 3 . We define the following stochastic processes u(t) and Y 2 (t): and Y 2 (t) = ln u(t).
Note that Y 2 (t) is well defined since by Theorem 3.2 u(t) > 0 for all t > 0 a.s.
We can calculate LY 2 (t) as where E(t) is given by This can be expressed as We are interested in lim sup t→∞ Y 2 (t), and we introduce the necessary notation for this analysis. For a stochastic process ψ(t), we write Note that for every w ∈ Ω, there exists an unbounded increasing sequence of positive numbers (t n ) such that lim n→∞ Y 2 (t n ) = lim sup t→∞ Y 2 (t) and such that the following sequences are convergent: I/u t n , J/u t n , R/u t n , (A -μS)/u t n , (Bθ V )/u t n . We denote their five limits by i, j, r, s # , v # , respectively, and we write Γ = lim sup t→∞ Y 2 t . We note that these values depend on the 4-tuple (q, p 1 , p 2 , p 3 ).   Consequently, by the strong law of large numbers for local martingales (see, e.g., [21]) we can deduce that lim t→∞ 1 t M(t) = 0 (a.s.).
Further, note that This means that lim sup The main stability theorem focuses on the convergence of I(t) and J(t). For the particular case p 1 = p 2 = p 3 = 0, we can calculate LY 2 (t) as follows, where we denote Y 2 by Y 0 : where E 0 (t) is given by This can be expressed as We require a suitable lower bound for E 0 . In particular, we further use the identity Thus we proceed as follows: with H(x) as in Remark 3.3. Therefore, noting the lower bound for H(x) established in Remark 3.3, we can deduce the inequality Again, using the identity I/u + qJ/u = 1, we obtain Now we can deduce the following inequality: If R(z) < 1, then I(t) and J(t) converge exponentially to 0 (a.s.).
Proof Let R(z) < 1, which is equivalent to the inequality We can find a number 0 < < 1 such that Now let For the given value of q, we now consider Y 0 . To prove our theorem, it suffices to prove that u(t) converges exponentially to zero (a.s.). The proof is concluded by showing that the Lyapunov exponent lim t→∞ 1 t Y 0 (t) of u(t) is negative almost surely. It suffices to prove that (see Remark 4.2) Recall inequality (4). Now performing the operationt and taking limits, we obtain A routine calculation reveals that If we substitute the value of q into the expression for C 1 , then from inequality (5) it follows that C 2 < 0. The coefficients of i and j on the right-hand side of the inequality are negative and constant. Note that Therefore at least one of the quantities i or j must be nonzero. Thus lim sup t→∞ < LY 0 > t < 0 (a.s.), and the proof is complete.
Remark 4.4 (a) As a particular case of Theorem 4.3, we have that if R(0) < 1, then I(t) and J(t) converge exponentially to 0 (a.s.). So if R(0) < 1, then starting from any initial value, the solutions of the stochastic system (1) almost surely converges to the diseasefree equilibrium.
(b) The theorem gives a version of local asymptotic stability. It can be interpreted as follows: if R(z) < 1 for some z > 0, then (a.s.) either the equilibrium X * is exponentially stable, or Therefore the theorem shows that for parameter values with R 0 slightly greater than 1, the stability of the disease-free equilibrium (i.e., extinction of the disease) is enhanced by stochastic perturbations. exponentially to 0 (a.s.), then the disease-free equilibrium X * is almost surely exponentially stable.

Theorem 4.5 If I(t) and J(t) converge
Proof In defining u(t), let us choose q = p 1 = p 2 = p 3 = 1. Suppose on the contrary that the stochastic process u(t) does not exponentially converge to zero. Then i = j = 0. Consequently, from Eq. (3) we can deduce the inequality This is a contradiction, which completes the proof.

Numerical simulations
We apply our model to malaria disease over a suitable region in Southern Africa. For sample simulations, we find numerical values for parameters from the literature. Such parameter values depend on the particular situation or region of application. There are parameters affected by human lifestyle, and others vary according to the different species of vector or species of the pathogen. We do not pick a specific small region of application, but our choice of parameters is relevant to malaria disease in Southern Africa. Our values are taken from the given sources either directly or possibly adjusted. The values of the perturbation parameters σ and ζ are declared for each simulation. The parameter values as per Table 1 yield R 0 = 1.201. Depends on the region. We assume a population of size 10,000 when disease-free. B Total birth rate of vector mosquitos 240,000θ Depends on the region. We assume a population of size 0.03 million. a The probability of a specific human getting bitten by a mosquito during a one-day period 6.417 × 10 -6 Chosen (this also depends on the region) b The probability that a bite by an infected mosquito will lead to a (new) human infection 0.075 [23,25] c The probability that a bite on an infected human will lead to a (new) mosquito infection 0.0375 cf. [23,25] k Transfer rate from I-class to R-class (recovery rate) 1 180 per day [5,9,25] h Transfer rate from R-class to S-class (rate of loss of temporary immunity) 1 2×365 per day [5,9,25]

The long-run mean
The initial values for our first simulation are

Projections with percentiles
In Fig. 2, we present some graphs obtained in a sample simulation. The initial values are chosen as The parameter values that we use here are strictly as in Table 1. The intensities of the perturbations are declared with the graph. We present the mean values of I(t) over n = 4000 sample paths together with, for every t value in the discretization, the first and the ninth deciles (respectively, denoted by I ten and I nty ) of I(t) for this collection of sample paths. The deciles form, for each t, an approximation of an 80% confidence interval. This interval estimation of future projection is an advantage of SDE modeling compared to ODE modeling. The graphs show, in particular, the behavior of the deciles in the long term.

Illustrating the stability theorem
We present numerical simulations to illustrate the results of Theorem 4.3 with parameter values given in Table 1 In order to showcase the stability theorem, we need to compare different situations of the model, having different values of the basic reproduction number. It has been found in [29] that the presence of livestock living near humans is correlated with a lower prevalence of malaria in humans. An explanation for this may be that the biting rate of mosquitos per human individual decreases. In the following simulations, we utilize the parameter values of Table 1, except that we replace the value of the parameter a with a smaller value, 5.866 × 10 -6 (as compared to the value 6.417 × 10 -6 calculated from Table 1), which arises from the assumption that the population now has livestock, kept in such a way as to reduce means that the deterministic model yields the instability of the disease-free equilibrium and stochastic perturbations cause the convergence (almost surely) to disease-free equilibrium the biting rate of vectors on humans. This leads to the situation R 0 = 1.002. Since R 0 > 1, for the underlying deterministic model, the disease-free equilibrium is not stable. However, for the stochastic model, in Fig. 3, we observe that the mean of I(t) seems to converge. This agrees with Theorem 4.3. Let us agree that R 0 is relatively close to 1. Then we expect for the deterministic case that the limiting value S * of S(t) will be very close to A/μ, and likewise, V * should be relatively close to B/θ . Let us make quite a conservative guess that for the stochastic case, lim inf(μS/A) × (θ V /B) will be above the value (0.95) × (0.95). This choice yields R(0.8145) = 0.998 < 1, for σ = ζ = 0.02. Indeed, then we do expect the convergence that we observe.

Conclusion
We have presented an SDE model for the population dynamics of a mosquito-borne disease, which has solutions that are almost surely positive and global. We have introduced an invariant R which is not higher than R 0 , with R 0 < 1 being a condition that guarantees the global stability of the disease-free equilibrium in the underlying deterministic model. With almost sure exponential stability holding for R < 1 (together with other requirements), it follows that the stochastic perturbation enhances the stability of the disease-free equilibrium. Simulations suggest (as expected) that, in general, the expectation of I(t)-values in the stochastic model is lower than its counterpart in the underlying deterministic model. Furthermore, we can make point estimates of forward projections of the compartment sizes, and also we can give estimates of confidence intervals (i.e., percentiles for a large number of simulations). We have consulted a variety of literature to find relevant numerical parameter values.
We mentioned in Remark 4.4 that for parameter values that yield R 0 slightly greater than 1, extinction of the disease is enhanced by the stochastic perturbations. Understanding of the asymptotic probability distributions of S and V will improve the applicability of Theorems 4.3 and 4.5. This is a problem that can be pursued in the future. As a next step, it will be interesting to study a similar model, but with perturbation on the force of infection, that is, with SI and VJ replaced by SJ and VI, respectively. Also, future work toward improving the model could address, in particular, the inclusion of classes of latent infection or time delays to account for incubation of the pathogen both in the host and in the vector.