Stochastic Dynamics of Cholera Epidemic Model: Formulation, Analysis and Numerical Simulation

Yohana Maiga Marwa^1,2*, Isambi Sailon Mbalawata³, Samuel Mwalili⁴, Wilson Mahera Charles⁵

¹Mathematics and Statistics Department, Institute for Basic Sciences, Technology and Innovation, Pan African University, Nairobi, Kenya.
²Tanzania Public Service College, Singida, Tanzania.
³African Institute for Mathematical Sciences, Dar es Salaam, Tanzania.
⁴Jomo Kenyatta University of Agriculture and Technology, Nairobi, Kenya.
⁵University of Dar es Salaam, Dar es Salaam, Tanzania.

DOI: 10.4236/jamp.2019.75074                

Abstract

In this paper, we describe the two different stochastic differential equations representing cholera dynamics. The first stochastic differential equation is formulated by introducing the stochasticity to deterministic model by parametric perturbation technique which is a standard technique in stochastic modeling and the second stochastic differential equation is formulated using transition probabilities. We analyse a stochastic model using suitable Lyapunov function and Itô formula. We state and prove the conditions for global existence, uniqueness of positive solutions, stochastic boundedness, global stability in probability, moment exponential stability, and almost sure convergence. We also carry out numerical simulation using Euler-Maruyama scheme to simulate the sample paths of stochastic differential equations. Our results show that the sample paths are continuous but not differentiable (a property of Wiener process). Also, we compare the numerical simulation results for deterministic and stochastic models. We find that the sample path of SI_sI_aR-B stochastic differential equations model fluctuates within the solution of the SI_sI_aR-B ordinary differential equation model. Furthermore, we use extended Kalman filter to estimate the model compartments (states), we find that the state estimates fit the measurements. Maximum likelihood estimation method for estimating the model parameters is also discussed.

Keywords

Stochastic Differential Equations, Stability Condition, Extended Kalman Filter, Itô Formula, Lyapunov Function, Euler-Maruyama Scheme

1. Introduction

It is well known that diseases have impacts on people’s health. Therefore, it is necessary to study the mechanism by which the disease spread, conditions for the disease to have minor or major outbreak and get the knowledge on how to control the diseases such as cholera, Ebola, etc. In this study we are mainly concern with cholera epidemic. Cholera is an infectious disease that causes severe watery diarrhea, which leads to dehydration and even death if untreated. According to WHO report [1] , about 1.4 to 4.3 million cases of cholera are reported each year worldwide and more than 140,000 deaths per year are reported due to cholera. As a result of this, several mathematical models have been developed to model cholera epidemics, through these models one can predict the behaviour of the disease and control the particular epidemic. Any epidemic of an infectious disease can be modelled by using either deterministic or stochastic models. The deterministic models are formulated as a system of ordinary different equations and are preferred by many researchers since its analysis is simple compared to stochastic models. However, the shortcomings of deterministic models are: they give less information, rely on the law of large numbers, difficulty to do estimation, when the population is very small it becomes difficult to do analysis and also, experimentally the measured trajectories do not behave as predicted due to some random effects that disturb the system [2] .

Due to these limitations of ordinary differential equations in modelling infectious diseases, stochastic modelling of infectious diseases in both heterogeneous and homogeneous population emerged as an alternative to deterministic model and alleviated some of the problems of deterministic models in modelling epidemic diseases. In reality many phenomena in nature are usually affected by stochastic noise and the ordinary differential Equation (ODEs) models ignore the stochastic effects [2] . The stochasticity can be added to the ODEs by including the random terms or elements by parametric perturbation technique. This technique introduces other parameters to the model known as noise intensities.

Many of the models that have been employed in water-borne settings have been deterministic, thus ignoring the possible effects of randomness; see, e.g. [3] [4] [5] [6] [7] . These models incorporate an environmental pathogen component that is the concentration of the vibrios into a SIR (susceptible-infected-recovered) and SI_sI_aR (susceptible-symptomatic infected-asymptomatic infected-recovered) epidemic framework. Some of these models only considered compartment I as a single compartment (individuals with severe symptoms only) (e.g. [3] [4] [5] [6] ) instead of splitting it into two heterogeneous groups that are symptomatic and asymptomatic infected individuals to study the transmission dynamics of cholera epidemic. Also, some of these models considered only direct transmission of disease i.e., human to human and other models ignored the feedback loop from infected individuals to the environment reservoir. However, in their studies no stochastic models considered, so as to keep track where the disease is at continuous time and not only at discrete time as shown in their studies.

In [7] they developed deterministic model by extending the work by [6] . The work in [6] considered infected individuals (I) as a homogeneous group that is people with severe symptoms like vomiting and diarrhoea. Therefore, [7] divided infected individuals (I) into symptomatic infected (I_s) and asymptomatic infected (I_a), in order to observe the contribution of concentration of Vibrio cholerae in the environment through excretion from each compartment and hence how it leads to the transmission dynamics of cholera epidemic.

In the stochastic modeling of cholera epidemic few papers emphasized on cholera stochastic models. [8] developed a deterministic model and further, extended it stochastic differential equations, the limitations to their paper are: no numerical analysis considered in their paper and compartment I is considered as single compartment, it could be better to split it into two groups, individuals with severe symptoms (symptomatic infected individuals) and those with mild symptoms (asymptomatic infected individuals) so as to observe the contribution of these two groups to the concentration of Vibrio cholerae through excretion.

[9] [10] developed a simple deterministic and stochastic model to discuss the spread of cholera. In their paper, they described the spread of cholera by modeling the bacteria population in contaminated water and human interaction with the bacteria in the water supply. The limitation to their paper is the absence of enough theoretical and numerical analysis. In [11] proposed a deterministic model that described the interaction among the two types of vibrios and viruses. The deterministic model proposed was further extended to include the random effects and from the stochastic model formulated it indicated that there is always a positive probability of disease extinction within the human host.

In this paper, we extend the deterministic model developed in [7] by formulating an equivalent stochastic differential Equation (SDEs). The formulated stochastic differential Equation (SDEs) models will be analyzed theoretically using suitable Lyapunov functions, Itô formula and some stochastic techniques. Numerical simulation will be carried using Euler-Maruyama scheme and extended Kalman filter. Also, maximum likelihood estimation method will be used to estimate the model parameters.

In the next section, we present a deterministic cholera epidemic model with water treatment, which is a model formulated by [7] . In Section 3, the corresponding two different stochastic models are formulated using parametric perturbation technique and transition probabilities approach. In Section 4, the formulated stochastic differential equations are analyzed by proving the existence and positivity of solutions, stochastic boundedness, global stability in probability, moment exponential stability, and almost sure convergence. In Section 5, the numerical results from the deterministic and stochastic simulations are presented, discussed and compared. In Section 6, we provide a brief discussion to support our findings. Ultimately, the last section concludes our paper.

2. Cholera Deterministic Model

From the description of the dynamics of cholera as shown in Figure 1, we the following set of ordinary differential equation system as in [7] .

$\begin{array}{l}\frac{\text{d}S\left(t\right)}{\text{d}t}=b{N}_e-\frac{\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}-\mu S\left(t\right),\\ \frac{\text{d}{I}_s\left(t\right)}{\text{d}t}=\frac{p\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}-\left(\mu +r_1+d\right)I_s\left(t\right),\\ \frac{\text{d}{I}_a\left(t\right)}{\text{d}t}=\frac{q\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}-\left(\mu +r_2\right)I_a\left(t\right),\\ \frac{\text{d}R\left(t\right)}{\text{d}t}=r_1{I}_s\left(t\right)+r_2{I}_a\left(t\right)-\mu R\left(t\right),\\ \frac{\text{d}B\left(t\right)}{\text{d}t}={\alpha }_1{I}_s\left(t\right)+{\alpha }_2{I}_a\left(t\right)-\left(\delta +\varphi \right)B\left(t\right),\end{array}$ (1)

with suitable initial conditions. The total human population is given by $N_e=S\left(t\right)+I_s\left(t\right)+I_a\left(t\right)+R\left(t\right)$.

From the model: $S\left(t\right)$, $I_s\left(t\right)$, $I_a\left(t\right)$, $R\left(t\right)$ and $B\left(t\right)$ refer to susceptible individuals, symptomatic infected individuals, asymptomatic infected individuals, recovered individuals and the pathogen concentration in the contaminated environment respectively. The model parameters are described as follows.

b is the birth rate or recruitment rate, $\beta$ is rate of exposure to contaminated water, $\kappa$ is the concentration of Vibrio cholerae in water, ${\alpha }_1$ is the contribution

Figure 1. The interaction of cholera epidemic transmission dynamics between compartments is shown.

of $I_s\left(t\right)$ to the population of Vibrio cholerae through excretion, ${\alpha }_2$ is the contribution of $I_a\left(t\right)$ to the population of Vibrio cholerae through excretion, d death rate due to cholera, $\delta$ death rate of Vibrio cholerae due to water treatment, $\varphi$ is the mortality rate for bacteria including phage degradation, $r_1$ is the recovery rate of symptomatic infected individuals, $r_2$ is the recovery rate of asymptomatic infected individuals, p is the prob. of new infected from $S\left(t\right)$ to be symptomatic, q is the prob. of new infected $S\left(t\right)$ to be asymptomatic , $\mu$ is the natural death rate.

The key parameter in epidemiology is the basic reproduction number, which is defined as the average number of secondary infectious cases transmitted by a single primary infectious cases introduced into a whole susceptible population [12] . This parameter is useful because it helps determine whether an infectious disease will spread within the population or not. To compute $R_0$, the next generation matrix approach is used as in [13] . It is obtained by taking the largest (dominant) eigenvalue value (spectral radius) of

$\left[\frac{\partial {\mathcal{F}}_i\left(E_0\right)}{\partial {\mathcal{X}}_i}\right]{\left[\frac{\partial {\mathcal{V}}_i\left(E_0\right)}{\partial {\mathcal{X}}_i}\right]}^{-1}\mathrm{,}$

where ${\mathcal{F}}_i$ is the rate of appearance of new infection in compartment i, ${\mathcal{V}}_i$ is the net transition between compartments, $E_0$ is the disease free equilibrium and ${\mathcal{X}}_i$ stand for the terms in which the infection is in progression i.e., $I_a\left(t\right)$, $I_a\left(t\right)$ and $B\left(t\right)$ in the model (1). Hence, the basic reproduction number obtained in [7] was as follows:

$R_0=\left[\frac{{\alpha }_1\beta \mu p+{\alpha }_2\beta q{r}_1+{\alpha }_1\beta p{r}_2+{\alpha }_2\beta dq+{\alpha }_2\beta \mu q}{D_n}\right]\frac{b{N}_e}{\mu },$ (2)

where

$\begin{array}{c}{D}_n=\delta \kappa \mu d+\delta \kappa {\mu }^2+\delta \kappa \mu r_1+\delta \kappa d{r}_2+\delta \kappa \mu r_2+\delta \kappa r_1{r}_2+\delta \kappa \mu \varphi \\ +\varphi \kappa {\mu }^2+\kappa \mu \varphi r_1+\kappa \varphi d{r}_2+\kappa \mu \varphi r_2+\kappa \varphi r_1{r}_2.\end{array}$

From Equation (2), when $R_0<1$, the highly infectious vibrios will not grow within human host and the environmental vibrios ingested into human body will not cause cholera infection. But when $R_0>1$ the human vibrios will grow and persist, and hence leads to human cholera.

3. Stochastic Differential Equations

In this section we provide two approaches of formulating stochastic differential equations from the deterministic model (1). These approaches are parametric perturbation Itô SDEs and Itô SDEs from the transition probabilities.

3.1. Parametric Perturbation Itô SDEs

The idea of parametric perturbation is to choose a parameter of interest from the deterministic model and change it to random variable [14] [15] . In this study we introduce the randomness into the model by replacing parameters $\beta$ and d by $\beta \mapsto \beta +{\sigma }_1\text{d}{\mathcal{W}}_1\left(t\right)$ and $d\mapsto d+{\sigma }_2\text{d}{\mathcal{W}}_2\left(t\right)$, where ${\mathcal{W}}_1$ and ${\mathcal{W}}_2$ are two independent standard Brownian motions with the following properties:

1) ${\mathcal{W}}_0=0$ $\mathbb{P}$ a.s,

2) For all times $0\le s<t$, the increment ${\mathcal{W}}_t-{\mathcal{W}}_s$ is normally distributed with zero mean and variance $t-s$ ; i.e., ${\mathcal{W}}_t-{\mathcal{W}}_s\sim N\left(\mathrm{0,}t-s\right)$,

3) For all times $0<t_1<t_2<\cdots <t_n$, the increments ${\mathcal{W}}_{t_1}\mathrm{,}{\mathcal{W}}_{t_2}\mathrm{,}\cdots \mathrm{,}{\mathcal{W}}_{t_n}-{\mathcal{W}}_{t_{n-1}}$ of the process are independent random variables,

4) All samples paths $\mathcal{W}\left(\mathrm{.,}\omega \right)\mathrm{<mo>\; :\; </mo>}\left[\mathrm{0,}\infty \right)\to ℝ\mathrm{,}\omega \in \Omega$, are continuous $\mathbb{P}$ a.s.

Similarly, ${\sigma }_1$ and ${\sigma }_2$ are real constants, known as the intensities of the stochastic environment.

From the deterministic model (1) we get the following system of stochastic differential equations:

$\text{d}S\left(t\right)=\left[b{N}_e-\frac{\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}-\mu S\left(t\right)\right]\text{d}t-\frac{{\sigma }_{1}B\left(t\right)S\left(t\right)\text{d}{\mathcal{W}}_1\left(t\right)}{\kappa +B\left(t\right)}\mathrm{,}$

$\begin{array}{c}\text{d}{I}_s\left(t\right)=\left[\frac{p\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}-\left(\mu +r_1+d\right)I_s\left(t\right)\right]\text{d}t-\frac{{\sigma }_{1}pB\left(t\right)S\left(t\right)\text{d}{\mathcal{W}}_1\left(t\right)}{\kappa +B\left(t\right)}\\ -{\sigma }_2{I}_s\left(t\right)\text{d}{\mathcal{W}}_2\left(t\right)\mathrm{,}\end{array}$

$\text{d}{I}_a\left(t\right)=\left[\frac{q\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}-\left(\mu +r_2\right)I_a\left(t\right)\right]\text{d}t-\frac{{\sigma }_{1}qB\left(t\right)S\left(t\right)\text{d}{\mathcal{W}}_1\left(t\right)}{\kappa +B\left(t\right)},$ (3)

$\text{d}R\left(t\right)=\left[r_1{I}_s\left(t\right)+r_2{I}_a\left(t\right)-\mu R\left(t\right)\right]\text{d}t,$

$\text{d}B\left(t\right)=\left[{\alpha }_1{I}_s\left(t\right)+{\alpha }_2{I}_a\left(t\right)-\left(\delta +\varphi \right)B\left(t\right)\right]\text{d}t,$

with suitable initial conditions. The technique of parameter perturbation introduces another parameters ${\sigma }_1$ and ${\sigma }_2$ in a model.

3.2. Itô SDEs with Transition Probabilities

This method was proposed by [16] [17] , the stochastic differential equations follow from the diffusion process. The nature of stochastic differential equation to be formulated is in this form:

$\text{d}x\left(t\right)=f\left(t\mathrm{,}x\left(t\right)\mathrm{<mo>\; ;\; </mo>}\theta \right)\text{d}t+G\left(t\mathrm{,}x\left(t\right)\mathrm{<mo>\; ;\; </mo>}\theta \right)\text{d}\mathcal{W}\left(t\right)\mathrm{,}$ (4)

where $G\left(t\mathrm{,}x\left(t\right)\mathrm{<mo>\; ;\; </mo>}\theta \right)$ is a matrix satisfying $G{G}^{\text{T}}=\Sigma$ [16] [17] , $\Sigma$ is the co-variance to order $\Delta t$, $\Sigma \Delta t$ is the approximate covariance matrix, $\mathcal{W}\left(t\right)$ is the vector of independent Wiener process, $\theta$ is a vector of parameters and $f\left(t\mathrm{,}x\left(t\right)\mathrm{<mo>\; ;\; </mo>}\theta \right)$ is the drift part or deterministic part.

From Equation (1) we formulate an equivalent SDEs as

$\text{d}x\left(t\right)=F_t\text{d}t+G_t\text{d}\mathcal{W}\left(t\right)\mathrm{,}$ (5)

where $F_t=\frac{\mathbb{E}\left[\Delta x\right]}{\Delta t}$ and $G_t=\sqrt{\frac{\mathbb{E}\left[\Delta x{\left(\Delta x\right)}^{\text{T}}\right]}{\Delta t}}$. Clearly, $x\left(t\right)$ is an Itô process,

such that $x\left(t\right)={\left[x_1\mathrm{,}{x}_2\mathrm{,}{x}_3\mathrm{,}{x}_4\right]}^{\text{T}}$. Then, $x_1\mathrm{,}{x}_2\mathrm{,}{x}_3\mathrm{,}{x}_4$, corresponds to the number of individuals $S\left(t\right)\mathrm{,}{I}_s\left(t\right)\mathrm{,}{I}_a\left(t\right)\mathrm{,}R\left(t\right)\in \left[\mathrm{0,}{N}_e\right]$ respectively and $B\left(t\right)$ the Vibrio cholerae concentration in the environment. However, $B\left(t\right)$ is not a compartment occupancy as the other one are so its transition is not considered in the formulation of transitions.

In order to find the expectation $\mathbb{E}\left[\Delta x\right]$ and covariance matrix $\Sigma$, we need to consider the transition probabilities as stated in Table 1. The transition probabilities are formulated from Equation (1). From Table 1, the expectation and variance co-variance matrix are computed as follows:

$\begin{array}{l}E\left[\Delta x\right]={\displaystyle \underset{i=1}{\overset{10}{\sum }}}{P}_i\left(\Delta x_i\right)=P_1\Delta x_1+P_2\Delta x_2+P_3\Delta x_3+P_4\Delta x_4+P_5\Delta x_5+P_6\Delta x_6\\ +P_7\Delta x_7+P_8\Delta x_8+P_9\Delta x_9+P_{10}\Delta x_{10}.\end{array}$ (6)

On substituting the values of $P_i$, $\Delta x_i$ and ${\left[x_1\mathrm{,}{x}_2\mathrm{,}{x}_3\mathrm{,}{x}_4\right]}^{\text{T}}={\left[S\left(t\right)\mathrm{,}{I}_s\left(t\right)\mathrm{,}{I}_a\left(t\right)\mathrm{,}R\left(t\right)\right]}^{\text{T}}$ to Equation (6), we get

$E\left[\Delta x\right]=\left[\begin{array}{c}b{N}_e-\frac{\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}-\mu S\left(t\right)\\ \frac{p\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}-\left(\mu +r_1+d\right)I_s\left(t\right)\\ \frac{q\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}-\left(\mu +r_2\right)I_a\left(t\right)\\ r_1{I}_s\left(t\right)+r_2{I}_a\left(t\right)-\mu R\left(t\right)\end{array}\right]\text{d}t$

and the co-variance matrix is given by

Table 1. Possible changes in the process for the $S{I}_s{I}_{a}R-B$ model.

$\begin{array}{c}E\left[\Delta x{\left(\Delta x\right)}^{\text{T}}\right]={\displaystyle \underset{i=1}{\overset{10}{\sum }}}{P}_i\left(\Delta x_i\right){\left(\Delta x_i\right)}^{\text{T}}\\ =P_1\left(\Delta x_1\right){\left(\Delta x_1\right)}^{\text{T}}+P_2\left(\Delta x_2\right){\left(\Delta x_2\right)}^{\text{T}}+P_3\left(\Delta x_3\right){\left(\Delta x_3\right)}^{\text{T}}\\ +P_4\left(\Delta x_4\right){\left(\Delta x_4\right)}^{\text{T}}+P_5\left(\Delta x_5\right){\left(\Delta x_5\right)}^{\text{T}}+P_6\left(\Delta x_6\right){\left(\Delta x_6\right)}^{\text{T}}\\ +P_7\left(\Delta x_7\right){\left(\Delta x_7\right)}^{\text{T}}+P_8\left(\Delta x_8\right){\left(\Delta x_8\right)}^{\text{T}}+P_9\left(\Delta x_9\right){\left(\Delta x_9\right)}^{\text{T}}\\ +P_{10}\left(\Delta x_{10}\right){\left(\Delta x_{10}\right)}^{\text{T}}.\end{array}$ (7)

Also, by substituting the values of $P_i$, $\Delta x_i$ and ${\left[x_1,x_2,x_3,x_4\right]}^{\text{T}}={\left[S\left(t\right),I_s\left(t\right),I_a\left(t\right),R\left(t\right)\right]}^{\text{T}}$ to Equation (7), we get

$\begin{array}{l}E\left[\Delta x{\left(\Delta x\right)}^{\text{T}}\right]\\ =\left[\begin{array}{cccc}{\sigma }_1^2& -\frac{p\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}& -\frac{q\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}& 0\\ -\frac{p\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}& {\sigma }_2^2& 0& -r_1{I}_s\left(t\right)\\ -\frac{q\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}& 0& {\sigma }_3^2& -r_2{I}_a\left(t\right)\\ 0& -r_1{I}_s\left(t\right)& -r_2{I}_a\left(t\right)& {\sigma }_4^2\end{array}\right]\Delta t,\end{array}$

where

${\sigma }_1^2=b{N}_e+\frac{\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}+\mu S\left(t\right),$

${\sigma }_2^2=\frac{p\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}+\left(\mu +r_1+d\right)I_s\left(t\right),$

${\sigma }_3^2=\frac{q\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}+\left(\mu +r_2\right)I_a\left(t\right),$

${\sigma }_4^2=r_1{I}_s\left(t\right)+r_2{I}_a\left(t\right)+\mu R\left(t\right).$

The Itô stochastic differential equation, where global Lipschitz conditions are satisfied to ensure the existence and uniqueness of strong solution, can be written as in Equation (4) as found in [14] . Hence by Equation (5), the Itô stochastic differential equations become:

$\text{d}x\left(t\right)=\left[\begin{array}{c}b{N}_e-\frac{\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}-\mu S\left(t\right)\\ \frac{p\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}-\left(\mu +r_1+d\right)I_s\left(t\right)\\ \frac{q\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}-\left(\mu +r_2\right)I_a\left(t\right)\\ r_1{I}_s\left(t\right)+r_2{I}_a\left(t\right)-\mu R\left(t\right)\end{array}\right]\text{d}t$

$+\left[\begin{array}{cccc}\sqrt{{\sigma }_1^2}& -\frac{p\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}& -\frac{q\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}& 0\\ -\frac{p\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}& \sqrt{{\sigma }_2^2}& 0& -r_1{I}_s\left(t\right)\\ -\frac{q\beta B\left(t\right)S\left(t\right)}{\kappa +B\left(t\right)}& 0& \sqrt{{\sigma }_3^2}& -r_2{I}_a\left(t\right)\\ 0& -r_1{I}_s\left(t\right)& -r_2{I}_a\left(t\right)& \sqrt{{\sigma }_4^2}\end{array}\right]\text{d}\mathcal{W}\left(t\right).$ (8)

Note that if $I_s\left(t\right)=I_a\left(t\right)=0$ and $B\left(t\right)=0$, then cholera epidemic stops.

4. Analysis of Parametric Perturbation Itô Stochastic Differential Equation Model

In this paper, unless otherwise stated, we let $\left(\Omega \mathrm{,}\mathcal{F}\mathrm{,}\left\{{\mathcal{F}}_t\right\}\mathrm{,}\mathbb{P}\right)$ be a complete filtered probability space, with ${\left\{{\mathcal{F}}_t\right\}}_{t\ge 0}$ satisfying the usual conditions (i.e., increasing and right continuous also ${\mathcal{F}}_0$ contains all $\mathbb{P}$ -null sets). Let

$\Delta =\left\{x\in {ℝ}_{+}^4:x_1+x_2+x_3+x_4<\frac{b{N}_e}{\mu }\right\},$

where

$x_1=S\left(t\right),x_2=I_s\left(t\right),x_3=I_a\left(t\right),x_4=R\left(t\right).$

4.1. Global Existence and Uniqueness of Positive Solutions

Before proving the existence and uniqueness of positive solutions, let us state the conditions that guarantee the existence and uniqueness of solution of Equation (4).

Lemma 1. Assume that there exist two positive constants $\hat{K}$ and K such that

1) (Lipschitz condition) for all $x\mathrm{,}y\in {ℝ}^n$ and $t\in \left[t_0\mathrm{,}T\right]$

${\left|f\left(x\mathrm{,}t\right)-f\left(y\mathrm{,}t\right)\right|}^2\le \hat{K}{\left|x-y\right|}^2$

2) (Linear growth condition) for all

$\left(x\mathrm{,}t\right)\in {ℝ}^n\times \left[t_0\mathrm{,}T\right]{\left|f\left(x\mathrm{,}t\right)\right|}^2\le K\left(1+{\left|x\right|}^2\right)\mathrm{.}$

Theorem 1. For any initial value $\left(S\left(0\right)\mathrm{,}{I}_s\left(0\right)\mathrm{,}{I}_a\left(0\right)\mathrm{,}R\left(0\right)\mathrm{,}B\left(0\right)\right)\in {ℝ}_{+}^5$, there is a unique solution $x\left(t\right)=S\left(t\right),I_s\left(t\right),I_a\left(t\right),R\left(t\right),B\left(t\right)$ to system (3) for all $t\ge 0$, and the solution will remain in ${ℝ}_{+}^5$ with the probability 1, namely $x\left(t\right)\in {ℝ}_{+}^5$, for all $t\ge 0$ almost surely.

Proof. Since the coefficients of system (3) satisfy the local Lipschitz condition, then for initial values $\left(S\left(0\right)\mathrm{,}{I}_s\left(0\right)\mathrm{,}{I}_a\left(0\right)\mathrm{,}R\left(0\right)\mathrm{,}B\left(0\right)\right)\in {ℝ}_{+}^5$, there is unique local solution $S\left(t\right),I_s\left(t\right),I_a\left(t\right),R\left(t\right),B\left(t\right)$ on $\left[\mathrm{0,}{\nu }_{\epsilon }\right]$, where ${\nu }_{\epsilon }$ is the first hitting time or explosion time. In order to prove that the solution is global and positive, we need to show that ${\nu }_{ϵ}=\infty$ a.s. Let $k_0>0$ be sufficiently large so

that $S\left(0\right)\mathrm{,}{I}_s\left(0\right)\mathrm{,}{I}_a\left(0\right)\mathrm{,}R\left(0\right)\mathrm{,}B\left(0\right)$ remain within the interval $\left[\frac{1}{k_0}\mathrm{,}{k}_0\right]$. For

integer $k\ge k_0$. Let ${\nu }_k\mathrm{<mo>\; :\; </mo>}\Omega \mapsto \left[\mathrm{0,}\infty \right]$ be any random variable (taking eventually infinite value) on a filtered probability space $\left(\Omega \mathrm{,}\mathcal{F}\mathrm{,}\left\{{\mathcal{F}}_t\right\}\mathrm{,}\mathbb{P}\right)$. We say that ${\nu }_k$ is a stopping time for filtration ${\mathcal{F}}_t$ if at each (deterministic) time $t\ge 0$, the event $\left[{\nu }_k\le t\right]$ is known in ${\mathcal{F}}_t$. In this study we define

${\nu }_k=\inf \left\{t\in \left[0,{\nu }_k\right):x_i\left(t\right)\notin \left(\frac{1}{k},k\right)\text{for}\text{some}i,1\le i\le 5\right\}$

Let $\varnothing$ be empty set then the $\inf \varnothing =\infty$. By definition ${\nu }_k$ increases as $k\to \infty$. Let ${\nu }_{\infty }={\lim }_{k\to \infty }{\nu }_k$, then ${\nu }_{\infty }\le {\nu }_{\epsilon }$ a.s. We need to show that ${\nu }_{\infty }=\infty$ a.s. Then ${\nu }_{\epsilon }=\infty$ and $x\left(t\right)\in {ℝ}_{+}^5$. If this statement is false we say that, there exist $\epsilon \in \left(\mathrm{0,1}\right)$ and a constant $T>0$, such that $\mathbb{P}\left\{{\nu }_{\infty }\le T\right\}>\epsilon$. As a consequence to this, there exist an integer $k_1>k_0$ such that

$\mathbb{P}\left\{{\nu }_k\le T\right\}\ge \epsilon \mathrm{,}\forall k\ge k_1$ (9)

For $t\le {\nu }_k$ and at each k,

$\begin{array}{c}\text{d}\left(S+I_s+I_a+R\right)=\left[\Lambda -\left(S+I_s+I_a+R\right)\mu -d{I}_s\right]\text{d}t\\ \le \left[\Lambda -\left(S+I_s+I_a+R\right)\mu \right]\text{d}t\mathrm{,}\end{array}$ (10)

where $\Lambda =b{N}_e$. Then

$S\left(t\right)+I_s\left(t\right)+I_a\left(t\right)+R\left(t\right)\le (\begin{array}{l}\frac{\Lambda }{\mu }\text{if}{N}_0\le \frac{\Lambda }{\mu }\\ N_0\text{if}{N}_0\ge \frac{\Lambda }{\mu }\end{array}\mathrm{<mo>\; :\; </mo>}=\mathbb{M}\mathrm{,}$ (11)

where

$N_0=S\left(0\right)+I_s\left(0\right)+I_a\left(0\right)+R\left(0\right),$

and

$\mathbb{M}=\max \left[\frac{\Lambda }{\mu },N_0\right].$

Define now ${ℂ}^{\mathrm{1,2}}$ the Lyapunov function $V\mathrm{<mo>\; :\; </mo>}{ℝ}_{+}^5\to {\overline{ℝ}}_{+}$ by

$\begin{array}{c}V\left(S\mathrm{,}{I}_s\mathrm{,}{I}_a\mathrm{,}R\mathrm{,}B\right)=\left(S-1-\text{In}S\right)+\left(I_s-1-\text{In}{I}_s\right)+\left(I_a-1-\text{In}{I}_a\right)\\ +\left(R-1-\text{In}R\right)+\left(B-1-\text{In}B\right)\mathrm{,}\end{array}$ (12)

and $S\left(t\right)\mathrm{,}{I}_s\left(t\right)\mathrm{,}{I}_a\left(t\right)\mathrm{,}R\left(t\right)\mathrm{,}B\left(t\right)\in \Omega$. Then $V\left(x\left(t\right)\mathrm{,}t\right)$ is an Itô stochastic process with the SDEs given by

$\begin{array}{l}\text{d}V\left(x\left(t\right)\mathrm{,}t\right)\\ =\left[V_t\left(x\left(t\right)\mathrm{,}t\right)+V_x\left(x\left(t\right)\mathrm{,}t\right)f\left(t\right)+\frac{1}{2}\text{trace}\left(g^{\text{T}}\left(t\right)V_{xx}\left(x\left(t\right)\mathrm{,}t\right)g\left(t\right)\right)\right]\text{d}t\\ +V_x\left(x\left(t\right)\mathrm{,}t\right)g\left(t\right)\text{d}\mathcal{W}\left(t\right)\mathrm{,}\end{array}$ (13)

which becomes

$\text{d}V\left(x\left(t\right),t\right)=LV\text{d}t+V_x\left(x\left(t\right)\mathrm{,}t\right)g\left(t\right)\text{d}\mathcal{W}\left(t\right),$ (14)

where

$\begin{array}{c}LV=\left(1-\frac{1}{S}\right)\Lambda -\frac{\beta BS}{\kappa +B}-\mu S+\left(1-\frac{1}{I_s}\right)\frac{p\beta BS}{\kappa +B}-\left(\mu +r_1+d\right)I_s\\ +\left(1-\frac{1}{I_a}\right)\frac{q\beta BS}{\kappa +B}-\left(\mu +r_2\right)I_a+\left(1-\frac{1}{R}\right)r_1{I}_s+r_2{I}_a-\mu R\\ +\left(1-\frac{1}{B}\right){\alpha }_1{I}_s+{\alpha }_2{I}_a-\left(\delta +\varphi \right)B+\frac{1}{2S^2}\left[\frac{{\sigma }_1^2{B}^2{S}^2}{{\left(\kappa +B\right)}^2}\right]\text{d}t\\ +\frac{1}{2I_s^2}\left[\frac{{\sigma }_1^2{p}^2{B}^2{S}^2}{{\left(\kappa +B\right)}^2}+{\sigma }_2^2{I}_s^2\right]\text{d}t+\frac{1}{2I_a^2}\left[\frac{{\sigma }_1^2{q}^2{B}^2{S}^2}{{\left(\kappa +B\right)}^2}\right]\text{d}t\end{array}$

Hence

$\begin{array}{c}LV\le \Lambda +\frac{\beta B}{\kappa +B}+\mu +\frac{p\beta BS}{\kappa +B}+\mu +r_1+d+\frac{q\beta BS}{\kappa +B}+\mu +r_2+r_1{I}_s\\ +r_2{I}_a+\mu +{\alpha }_1{I}_s+{\alpha }_2{I}_a+\delta +\varphi +\frac{{\sigma }_1^2{B}^2}{2{\left(\kappa +B\right)}^2}+\frac{{\sigma }_1^2{p}^2{B}^2{S}^2}{2I_s^2{\left(\kappa +B\right)}^2}\\ +\frac{{\sigma }_2^2}{2}+\frac{S^2}{2I_a^2}\frac{{\sigma }_1^2{q}^2{B}^2}{{\left(\kappa +B\right)}^2},\end{array}$ (15)

From $\mathbb{M}=\max \left[\frac{\Lambda }{\mu },N_0\right]$, it follows that

$\begin{array}{c}LV\le \Lambda +\left(\frac{p\beta B}{\kappa +B}+\frac{q\beta B}{\kappa +B}+r_1+r_2+{\alpha }_1+{\alpha }_2\right)\mathbb{M}+\frac{\beta B}{\kappa +B}+\delta +\varphi \\ +\frac{{\sigma }_1^2{B}^2}{2\left(\kappa +B\right)}+\frac{{\sigma }_1^2{p}^2{B}^2{\mathbb{M}}^2}{2I_s^2{\left(\kappa +B\right)}^2}+\frac{{\sigma }_2^2}{2}+4\mu +r_2+r_1+d+\frac{{\mathbb{M}}^2}{2I_a^2}\frac{{\sigma }_1^2{q}^2{B}^2}{{\left(\kappa +B\right)}^2}\\ :=\mathbb{K},\end{array}$ (16)

where $\mathbb{K}$ is a positive constant independent of variables $S\mathrm{,}{I}_s\mathrm{,}{I}_a\mathrm{,}R\mathrm{,}B$ and time t. Then we obtain

$\text{d}V\le \mathbb{K}\text{d}t-\left(\frac{S}{I_s}-\frac{1}{p}\right)\frac{{\sigma }_{1}pB\text{d}{\mathcal{W}}_1\left(t\right)}{\kappa +B}-\frac{S}{I_a}\frac{{\sigma }_{1}qB\text{d}{\mathcal{W}}_1\left(t\right)}{\kappa +B}-{\sigma }_2\text{d}{\mathcal{W}}_2\left(t\right)\mathrm{,}$ (17)

let ${\nu }_k\wedge T=\min \left\{{\nu }_k,T\right\}$. Introduce integral to Equation (17)

$\begin{array}{l}{\displaystyle {\int }_0^{{\nu }_k\wedge T}}\text{d}V\left(S\left(s\right),I_s\left(s\right),I_a\left(s\right),R\left(s\right),B\left(s\right)\right)\\ ={\displaystyle {\int }_0^{{\nu }_k\wedge T}}\mathbb{K}\text{d}s-{\displaystyle {\int }_0^{{\nu }_k\wedge T}}\left(\frac{S}{I_s}-\frac{1}{p}\right)\frac{{\sigma }_{1}pB\text{d}{\mathcal{W}}_1\left(s\right)}{\kappa +B}\\ -{\displaystyle {\int }_0^{{\nu }_n\wedge T}}{\sigma }_2\text{d}{\mathcal{W}}_2\left(s\right)-{\displaystyle {\int }_0^{{\nu }_n\wedge T}}\frac{S}{I_a}\frac{{\sigma }_{1}qB\text{d}{\mathcal{W}}_1\left(s\right)}{\kappa +B}.\end{array}$ (18)

Introducing the expectation to Equation (18) and by the Gronwall inequality, leads to

$\mathbb{E}\left[V\left(S\left({\nu }_k\wedge T\right)\mathrm{,}{I}_s\left({\nu }_k\wedge T\right)\mathrm{,}{I}_a\left({\nu }_k\wedge T\right)\mathrm{,}R\left({\nu }_k\wedge T\right)\mathrm{,}B\left({\nu }_k\wedge T\right)\right)\right]\le C+\mathbb{K}T\mathrm{,}$ (19)

where $C=V\left(S\left(0\right),I_s\left(0\right),I_a\left(0\right),R\left(0\right),B\left(0\right)\right)$. Since the mean of stochastic integral equals zero [18] . Then let ${\Omega }_k=\left\{{\nu }_k\le T\right\},\forall k\ge k_1$ from inequality (9), we have $\mathbb{P}\left({\Omega }_k\right)\ge ϵ$. Then for all $\omega \in {\Omega }_k$, there must exist at least one of

$S\left({\nu }_k\mathrm{,}T\right)\mathrm{,}{I}_s\left({\nu }_k\mathrm{,}T\right)\mathrm{,}{I}_a\left({\nu }_k\mathrm{,}T\right)\mathrm{,}R\left({\nu }_k\mathrm{,}T\right)\mathrm{,}B\left({\nu }_k\mathrm{,}T\right)\mathrm{,}$

which can take either $\frac{1}{k}$ or k. Hence we have

$\begin{array}{l}V\left(S\left({\nu }_k\mathrm{,}T\right)\mathrm{,}{I}_s\left({\nu }_k\mathrm{,}T\right)\mathrm{,}{I}_a\left({\nu }_k\mathrm{,}T\right)\mathrm{,}R\left({\nu }_k\mathrm{,}T\right)\mathrm{,}B\left({\nu }_k\mathrm{,}T\right)\right)\\ \ge \left[k-1-\log \left(k\right)\right]\wedge \left[\frac{1}{k}-1-\log \left(\frac{1}{k}\right)\right],\end{array}$ (20)

again define the indicator function on ${\Omega }_k$ denoted by $⊮$ be defined for all $\omega \in {\Omega }_k$ by

From Equation (19) and by Bienayme Chebyshev inequality we have

(21)

when $k\to +\infty$, this leads to contradiction as

$\infty >V\left(S\left(0\right),I_s\left(0\right),I_a\left(0\right),R\left(0\right),B\left(0\right)\right)+\mathbb{K}T=\infty$

Therefore, we conclude that for the solution $S\left(t\right),I_s\left(t\right),I_a\left(t\right),R\left(t\right),B\left(t\right)$ of cholera model not to explode at finite time with probability 1, we need to have ${\nu }_{\infty }=\infty$ a.s. This completes the proof. □

Theorem 1 shows that the solution of model (3) will remain in ${ℝ}_{+}^5$. Next, we give the definition of stochastic ultimate boundedness [19] , which is very important in population dynamics.

Definition 1. The solution $x\left(t\right)=\left(S\left(t\right),I_s\left(t\right),I_a\left(t\right),R\left(t\right),B\left(t\right)\right)$ of model (3) are said to be stochastically ultimately bounded, if for any $ϵ\in \left(\mathrm{0,1}\right)$, there is a positive constant $\delta =\delta \left(ϵ\right)$, such that for any initial value $x_0=\left(S\left(0\right),I_s\left(0\right),I_a\left(0\right),R\left(0\right),B\left(0\right)\right)\in {ℝ}_{+}^5$, the solution $x\left(t\right)$ to Equation (3) has a property that

$\underset{t\to \infty }{\lim }\sup \mathbb{P}\left\{\left|\sqrt{x\left(t\right)}\right|>\delta \right\}<ϵ.$

Theorem 2. The solutions of system (3) are stochastically ultimately bounded for any initial value $x_0\in {ℝ}_{+}^5$.

Proof. From Theorem 1, we showed that $S\left(t\right),I_s\left(t\right),I_a\left(t\right),R\left(t\right),B\left(t\right)$ remain in ${ℝ}_{+}^5$, for all $t\ge 0$ almost surely.

Let $V_1$, $V_2$ and $V_3$ be lyapunov functions stated as

$V_1={\text{e}}^t{S}^{\theta },V_2={\text{e}}^t{I}_s^{\theta },V_3={\text{e}}^t{I}_a^{\theta }\text{where}\theta >1.$ (22)

Taking their derivatives with respect to $\theta$ we get

${V^{\prime }}_1\left(\theta \right)=\frac{\theta {\text{e}}^t{S}^{\theta }}{S},{V^{\prime }}_2\left(\theta \right)=\frac{\theta {\text{e}}^t{I}_s^{\theta }}{I_s}\text{and}{V^{\prime }}_3\left(\theta \right)=\frac{\theta {\text{e}}^t{I}_a^{\theta }}{I_a}.$ (23)

Using Itô formula we get

$\text{d}{V}_1=L{V}_1\text{d}t-\frac{{\sigma }_1\theta B{\text{e}}^t{S}^{\theta }\text{d}{\mathcal{W}}_1\left(t\right)}{\kappa +B},$

$\text{d}{V}_2=L{V}_2\text{d}t-\frac{{\sigma }_{1}p\theta B{\text{e}}^t{S}^{\theta }\text{d}{\mathcal{W}}_1\left(t\right)}{\kappa +B}-{\sigma }_2\theta {\text{e}}^t{I}_s^{\theta }\text{d}{\mathcal{W}}_2\left(t\right),$ (24)

$\text{d}{V}_3=L{V}_3\text{d}t-\frac{{\sigma }_{1}q\theta B{\text{e}}^t{S}^{\theta }\text{d}{\mathcal{W}}_1\left(t\right)}{\kappa +B},$

but

$L{V}_1=\frac{\partial V_1}{\partial t}+V_{1\theta }f\left(t\right)+\frac{1}{2}\text{trace}\left(g_1^{\text{T}}\left(t\right)V_{1\theta \theta }{g}_1\left(t\right)\right)\mathrm{,}$

$L{V}_2=\frac{\partial V_2}{\partial t}+V_{2\theta }f\left(t\right)+\frac{1}{2}\text{trace}\left(g_2^{\text{T}}\left(t\right)V_{2\theta \theta }{g}_2\left(t\right)\right)\mathrm{,}$ (25)

$L{V}_3=\frac{\partial V_3}{\partial t}+V_{3\theta }f\left(t\right)+\frac{1}{2}\text{trace}\left(g_3^{\text{T}}\left(t\right)V_{3\theta \theta }{g}_3\left(t\right)\right)\mathrm{,}$

where

$V_{1\theta \theta }=\frac{\theta \left(\theta -1\right){\text{e}}^t{S}^{\theta }}{S^2},V_{2\theta \theta }=\frac{\theta \left(\theta -1\right){\text{e}}^t{I}_s^{\theta }}{I_s^2},V_{3\theta \theta }=\frac{\theta \left(\theta -1\right){\text{e}}^t{I}_a^{\theta }}{I_a^2},$ (26)

$g_1\left(t\right)=\left[\frac{-{\sigma }_{1}BS}{\kappa +B}\right],g_2\left(t\right)=\left[\frac{{\sigma }_{1}pBS}{\kappa +B}-{\sigma }_2{I}_s\right],g_3\left(t\right)=\left[\frac{{\sigma }_{1}qBS}{\kappa +B}\right].$ (27)

Therefore

$L{V}_1={\text{e}}^t{S}^{\theta }\left[1+\theta \left(\frac{\Lambda }{S}-\frac{\beta B}{\kappa +B}-\mu \right)\right]+\frac{1}{2}\frac{\theta \left(\theta -1\right){\text{e}}^t{S}^{\theta }{\sigma }_1^2{B}^2}{{\left(\kappa +B\right)}^2},$

$\begin{array}{c}L{V}_2={\text{e}}^t{I}_s^{\theta }\left[1+\theta \left(\frac{p\beta BS}{I_s\left(\kappa +B\right)}-\mu -r_1-d\right)\right]+\frac{1}{2}\frac{{\sigma }_1^2{p}^2{B}^2{S}^2\theta \left(\theta -1\right){\text{e}}^t{I}_s^{\theta }}{I_s^2{\left(\kappa +B\right)}^2}\\ +\frac{1}{2}{\sigma }_2^2\theta \left(\theta -1\right){\text{e}}^t{I}_s^{\theta },\end{array}$ (28)

$L{V}_3={\text{e}}^t{I}_a^{\theta }\left[1+\theta \left(\frac{q\beta BS}{I_a\left(\kappa +B\right)}-\mu -r_2\right)\right]+\frac{1}{2}\frac{{\sigma }_1^2{q}^2{B}^2{S}^2\theta \left(\theta -1\right){\text{e}}^t{I}_a^{\theta }}{I_a^2{\left(\kappa +B\right)}^2}.$

There exists positive constants ${\mathbb{K}}_1$, ${\mathbb{K}}_2$ and ${\mathbb{K}}_3$ such that

$L{V}_1<{\mathbb{K}}_1{\text{e}}^t,L{V}_2<{\mathbb{K}}_2{\text{e}}^t\text{and}L{V}_3<{\mathbb{K}}_3{\text{e}}^t.$ (29)

It follows that

$\begin{array}{l}\mathbb{E}\left(S^{\theta }\left(t\right)\right)-\mathbb{E}\left(S^{\theta }\left(0\right)\right)\le {\mathbb{K}}_1,\mathbb{E}\left(I_s^{\theta }\left(t\right)\right)-\mathbb{E}\left(I_s^{\theta }\left(0\right)\right)\le {\mathbb{K}}_2\\ \text{and}\mathbb{E}\left(I_a^{\theta }\left(t\right)\right)-\mathbb{E}\left(I_a^{\theta }\left(0\right)\right)\le {\mathbb{K}}_3.\end{array}$ (30)

Then

$\begin{array}{l}\underset{t\to \infty }{\lim }\sup \mathbb{E}\left[S^{\theta }\left(t\right)\right]\le {\mathbb{K}}_1<\infty ,\underset{t\to \infty }{\lim }\sup \mathbb{E}\left[I_s^{\theta }\left(t\right)\right]\le {\mathbb{K}}_2<\infty \\ \text{and}\underset{t\to \infty }{\lim }\sup \mathbb{E}\left[I_a^{\theta }\left(t\right)\right]\le {\mathbb{K}}_3<\infty .\end{array}$ (31)

For $x\left(t\right)=\left(S\left(t\right),I_s\left(t\right),I_a\left(t\right)\right)\in {ℝ}_{+}^3$ and consider the case of $0<\theta <3$. Then from the holders inequality we have

$\begin{array}{c}{\left|x\left(t\right)\right|}^{\theta }={\left(S^2\left(t\right)+I_s^2\left(t\right)+I_a^2\left(t\right)\right)}^{\frac{\theta }{2}}\le 3^{\frac{\theta }{2}}\max \left\{S^{\theta }\left(t\right),I_s^{\theta }\left(t\right),I_a^{\theta }\left(t\right)\right\}\\ \le 3^{\frac{\theta }{2}}\left\{S^{\theta }\left(t\right)+I_s^{\theta }\left(t\right)+I_a^{\theta }\left(t\right)\right\}\end{array}$ (32)

consequently we have

$\underset{t\to \infty }{\lim }\sup \mathbb{E}\left|x\left(t\right)\right|\le 3^{\frac{\theta }{2}}{\mathbb{K}}_4\mathrm{,}$ (33)

where ${\mathbb{K}}_4>0$ is a suitable constant.

There exists a positive constants ${\delta }_1$ such that

$\underset{t\to \infty }{\lim }\sup \mathbb{E}\left|\sqrt{x\left(t\right)}\right|\le {\delta }_1\mathrm{.}$

Given $ϵ>0$ and let $\delta =\frac{{\delta }_1^2}{{ϵ}^2}$, then by Bienayme Chebyshev inequality

$\mathbb{P}\left\{\left|x\left(t\right)\right|>\delta \right\}\le \frac{\mathbb{E}\left|\sqrt{x\left(t\right)}\right|}{\sqrt{\delta }}.$ (34)

Hence

$\underset{t\to \infty }{\lim }\sup \mathbb{P}\left\{\left|x\left(t\right)\right|>\delta \right\}\le \frac{{\delta }_1}{\sqrt{\delta }}=ϵ.$

This completes the proof as required. □

4.2. Moment Exponential Stability

The moment exponential stability of the equilibrium solutions of stochastic differential Equation (3) are established from the idea of Lyapunov function [20] .

Lemma 2. Consider a function ${ℂ}^{\mathrm{1,2}}\left({ℝ}^n\times \left[t\mathrm{,}\infty \right]\mathrm{<mo>\; ;\; </mo>}{ℝ}_{+}\right)$ satisfying these inequalities:

$M_1{\left|x\right|}^p\le V\left(x\mathrm{,}t\right)\le M_2{\left|x\right|}^p$

and

$LV\left(x\mathrm{,}t\right)\le -M_3{\left|x\right|}^p$, $t\ge \mathrm{0,}\forall \left(x\mathrm{,}t\right)\in {ℝ}^n\times \left[t\mathrm{,}\infty \right]$

where $M_1\mathrm{,}{M}_2\mathrm{,}{M}_3$ and p are positive constants. Then the equilibrium solution of Equation (3) is pth moment exponentially stable. When $p=2$, it is usually said to be exponentially stable in mean square and the DFE is globally asymptotically stable.

Lemma 3. When $p\ge 2$ and $x,y,ϵ>0$. Then

$x^{p-1}y\le \frac{\left(p-1\right)ϵx^p}{p}+\frac{1}{p}{ϵ}^{1-p}{y}^p\text{and}{x}^{p-2}{y}^2\le \frac{\left(p-2\right)ϵx^p}{p}+\frac{2}{p}{y}^p{\text{e}}^{\left(2-p\right)/2}\mathrm{.}$

Theorem 3. If $R_0<1$ and $p\ge 2$, the disease free equilibrium of the model system (3) is pth moment exponentially stable in $\Delta$.

Proof. Given that $\left(S\left(0\right)\mathrm{,}{I}_s\left(0\right)\mathrm{,}{I}_a\left(0\right)\mathrm{,}R\left(0\right)\mathrm{,}B\left(0\right)\right)\in \Delta$, as from Theorem 1 the solution of the system remains in $\Delta$. Set $p\ge 2$, then consider Lyapunov

function $V={\lambda }_1{\left(\frac{\Lambda }{\mu }-S\right)}^p+\frac{I_s^p}{p}+{\lambda }_3{I}_a^p+{\lambda }_4{R}^p+{\lambda }_5{B}^p$, we have

$\begin{array}{c}LV=-{\lambda }_{1}p\Lambda {\left(\frac{\Lambda }{\mu }-S\right)}^{p-1}+\frac{{\lambda }_{1}p\beta BS}{\kappa +B}{\left(\frac{\Lambda }{\mu }-S\right)}^{p-1}+{\lambda }_{1}p{\left(\frac{\Lambda }{\mu }-S\right)}^{p-1}\mu S\\ +\frac{1}{2}{\lambda }_{1}p\left(p-1\right){\left(\frac{\Lambda }{\mu }-S\right)}^{p-2}\frac{{\sigma }_1^2{B}^2{S}^2}{{\left(\kappa +B\right)}^2}+I_s^{p-1}\left(\frac{p_s\beta BS}{\kappa +B}-\left(\mu +r_1+d\right)I_s\right)\\ +\frac{1}{2}\left(p-1\right)I_s^{p-2}\left(\frac{p_s^2{\sigma }_1^2{B}^2{S}^2}{{\left(\kappa +B\right)}^2}+{\sigma }^2{I}_s^2\right)+{\lambda }_{3}p{I}_a^{p-1}\left(\frac{q\beta BS}{\kappa +B}-\left(\mu +r_2\right)I_a\right)\\ +{\lambda }_{4}p{R}^{p-1}\left(r_1{I}_s+r_2{I}_a-\mu R\right)+{\lambda }_{5}p{B}^{p-1}\left({\alpha }_1{I}_s+{\alpha }_2{I}_a-\delta B-\varphi B\right)\\ +\frac{1}{2}{\lambda }_{3}p\left(p-1\right)I_a^{p-2}\frac{{\sigma }_1^2{q}^2{B}^2{S}^2}{{\left(\kappa +B\right)}^2}.\end{array}$ (35)

In $\Delta$, we have $S\left(0\right)\mathrm{,}{I}_s\left(0\right)\mathrm{,}{I}_a\left(0\right)\mathrm{,}R\left(0\right)\mathrm{,}B\left(0\right)\in \left(\frac{\Lambda }{\mu }\mathrm{,0,0,0,0}\right)$, such that $S\le \frac{\Lambda }{\mu }$, then from the fact that $\frac{B}{\kappa +B}\le 1$, as $B\in \left[\mathrm{0,1}\right]$, there exists $\eta \in \left[\mathrm{0,1}\right]$ such that $\frac{\eta {\Lambda }^2}{{\mu }^2}\left(x^2+y^2\right)\le S^2{x}^{n-1}{y}^2$, Then LV becomes

$\begin{array}{c}LV\le -{\lambda }_{1}p\Lambda {\left(\frac{\Lambda }{\mu }-S\right)}^{p-1}+\left(\beta +\mu \right){\lambda }_{1}p\frac{\Lambda }{\mu }{\left(\frac{\Lambda }{\mu }-S\right)}^{p-1}\\ +\frac{1}{2}{\lambda }_{1}p\left(p-1\right)\eta {\sigma }_1^2\frac{{\Lambda }^2}{{\mu }^2}{\left(\frac{\Lambda }{\mu }-S\right)}^{p-2}+p_s\beta \frac{\Lambda }{\mu }{I}_s^{p-1}-\left(\mu +r_1+d\right)I_s^p\\ +\frac{1}{2}\left(p-1\right)\eta p_s^2{\sigma }_1^2\frac{{\Lambda }^2}{{\mu }^2}{I}_s^{p-2}+\frac{1}{2}{\sigma }_2^2\left(p-1\right)I_s^p+{\lambda }_{3}qp\beta \frac{\Lambda }{\mu }{I}_a^{p-1}\\ -{\lambda }_{3}p\left(\mu +r_2\right)I_a^p+{\lambda }_{4}p{r}_1{I}_s{R}^{p-1}+{\lambda }_{4}p{r}_2{I}_a{R}^{p-1}-{\lambda }_4\mu p{R}^p+{\lambda }_5{\alpha }_{1}p{I}_s{B}^{p-1}\\ +{\lambda }_5{\alpha }_{2}p{I}_a{B}^{p-1}-{\lambda }_{5}p\left(\delta +\varphi \right)B^p+\frac{1}{2}{\lambda }_{3}p\left(p-1\right)\eta {\sigma }_1^2{q}^2\frac{{\Lambda }^2}{{\mu }^2}{I}_a^{p-2}.\end{array}$ (36)

Applying the idea of Lemma 3, we can formulate some of the inequalities as follows

$I_s{R}^{p-1}\le \frac{p-1}{p}ϵR^p+\frac{1}{p}{ϵ}^{1-p}{I}_s^p\mathrm{,}$

$I_s^2{\left(\frac{\Lambda }{\mu }-S\right)}^{p-2}\le \frac{p-2}{p}ϵ{\left(\frac{\Lambda }{\mu }-S\right)}^p+\frac{2}{p}{ϵ}^{\frac{2-p}{2}}{I}_s^p\mathrm{,}$ (37)

$N_e{\left(\frac{\Lambda }{\mu }-S\right)}^{p-1}\le \frac{p-1}{p}ϵ{\left(\frac{\Lambda }{\mu }-S\right)}^p+\frac{1}{p}{ϵ}^{1-p}{N}_e^p\mathrm{.}$

On substituting these inequalities in Equation (19), we get

$\begin{array}{l}LV\le -{\lambda }_{1}b\left(p-1\right)ϵ{\left(\frac{\Lambda }{\mu }-S\right)}^p+\left(\beta +\mu \right)\frac{b}{\mu }{\lambda }_1\left(p-1\right)ϵ{\left(\frac{\Lambda }{\mu }-S\right)}^p\\ +\frac{1}{2}{\lambda }_1\left(p-1\right)\left(p-2\right)\eta {\sigma }_1^2\frac{b^2}{{\mu }^2}ϵ{\left(\frac{\Lambda }{\mu }-S\right)}^p+\frac{p_s^2}{2p}\left(p-1\right)\left(p-2\right)\eta {\sigma }_1^2\frac{b^2}{\mu }ϵ{\left(\frac{\Lambda }{\mu }-S\right)}^P\\ -\left(\mu +r_1+d-\frac{1}{2}\left(p-1\right){\sigma }_2^2\right)I_s^p+{\lambda }_4\left(p-1\right)r_1ϵR^p+{\lambda }_4\left(p-1\right)r_2ϵR^p\\ +{\lambda }_5{\alpha }_2\left(p-1\right)ϵB^p+{\lambda }_5{\alpha }_1\left(p-1\right)ϵB^p+\frac{p_s}{p}\left(p-1\right)\beta \frac{\Lambda }{\mu }ϵI_s^p\end{array}$

$\begin{array}{c}+{\lambda }_{3}q\left(p-1\right)\beta \frac{\Lambda }{\mu }ϵI_a^p-{\lambda }_{3}p\left(\mu +r_2\right)I_a^p-{\lambda }_4\mu p{R}^p\\ -{\lambda }_{5}p\left(\delta +\varphi \right)B^p+\frac{1}{2}{\lambda }_3\left(p-1\right)\left(p-2\right)\eta {\sigma }_1^2{q}^2\frac{b^2}{{\mu }^2}\eta {\left(\frac{\Lambda }{\mu }-S\right)}^p\\ +\frac{1}{2p}\left(p-1\right)\left(p-2\right)\eta p_s^2{\sigma }_1^2\frac{b^2}{{\mu }^2}ϵI_s^p+{\lambda }_4{r}_1{ϵ}^{1-p}{I}_s^p+{\lambda }_4{r}_2{ϵ}^{1-p}{I}_a^p\\ +{\lambda }_5{\alpha }_1{ϵ}^{1-p}{I}_s^p+{\lambda }_5{\alpha }_2{ϵ}^{1-p}{I}_a^p.\end{array}$ (38)

Hence

$\begin{array}{c}LV\le -[\mu +r_1+d-\frac{1}{2}\left(p-1\right){\sigma }_s^2+\frac{p_s}{p}\left(p-1\right)\beta \frac{\Lambda }{\mu }ϵ\\ -\frac{1}{2p}\left(p-1\right)\left(p-2\right)\eta p_s^2{\sigma }_1^2\frac{b^2}{{\mu }^2}ϵ-{\lambda }_4{r}_1{ϵ}^{1-p}-{\lambda }_5{\alpha }_1{ϵ}^{1-p}]I_s^p\\ -\left[r_1\left(p-1\right)+r_2\left(p-1\right)+\mu p\right]{\lambda }_4ϵR^p\\ -\left[p\left(\delta +\varphi \right)+{\alpha }_1\left(p-1\right)ϵ+{\alpha }_2\left(p-1\right)ϵ\right]{\lambda }_5{B}^p\\ -\left[\left(\left(p-1\right)\beta bϵ\frac{q}{\mu }+p\left(\mu +r_2\right)\right){\lambda }_3-{\lambda }_4{r}_2{ϵ}^{1-p}-{\lambda }_5{\alpha }_2{ϵ}^{1-p}\right]I_a^p\end{array}$

$\begin{array}{l}-\left[b-\left(\mu +\beta \right)\frac{b}{\mu }-\frac{1}{2}\left(p-2\right)\eta {\sigma }_1^2\frac{b^2}{{\mu }^2}-F_1\right]{\lambda }_1\left(p-1\right)ϵ{\left(\frac{\Lambda }{\mu }-S\right)}^p\\ -[{\lambda }_{1}b{ϵ}^{1-p}-\left(\beta +\mu \right)\frac{b}{\mu }{ϵ}^{1-p}{\lambda }_1-\left(p-1\right)\eta {\sigma }_1^2\frac{b^2}{{\mu }^2}{\eta }^{\frac{2-p}{2}}{\lambda }_1\\ -\frac{p_s}{p}\beta \frac{b}{\mu }{ϵ}^{1-p}-\frac{p-1}{p}\eta p_s^2{\sigma }_1^2\frac{b^2}{{\mu }^2}{ϵ}^{\frac{2-p}{2}}-F_2]N_e^p\mathrm{,}\end{array}$ (39)

where $F_1=\frac{p_s^2}{2p{\mu }^2{\lambda }_1}\left(p-2\right)\eta {\sigma }_1^2{b}^2+\frac{{\lambda }_3\left(p-2\right)\eta {\sigma }_1^2{q}^2}{2{\lambda }_1}\frac{b^2}{{\mu }^2}$ and $F_2={\lambda }_{3}q\beta \frac{b}{\mu }{ϵ}^{1-p}+{\lambda }_3\left(p-1\right)\eta {\sigma }_1^2{q}^2\frac{b^2}{{\mu }^2}{ϵ}^{\frac{2-p}{2}}$. If we choose $ϵ$ sufficiently very small, then we can choose ${\lambda }_1\mathrm{,}{\lambda }_2\mathrm{,}{\lambda }_3\mathrm{,}{\lambda }_4\mathrm{,}{\lambda }_5$ positive such that the coefficients ${\left(\frac{\Lambda }{\mu }-S\right)}^p\mathrm{,}{I}_s^p\mathrm{,}{I}_a^p\mathrm{,}{R}^p\mathrm{,}{B}^p\mathrm{,}{N}_e^p$ become negative.

Hence, according to Lemma 2 the proof is complete. □

4.3. Almost Sure Convergence

Theorem 4 If $R_0<1$, then $I_s\left(t\right)\mathrm{,}{I}_a\left(t\right)\mathrm{,}R\left(t\right)$ converge almost surely exponentially to (0, 0, 0).

Proof. Let $\left(S\left(0\right)\mathrm{,}{I}_s\left(0\right)\mathrm{,}{I}_a\left(0\right)\mathrm{,}R\left(0\right)\right)\in \Delta$. As $R_0<1$, define the lyapunov function as $V\left(I_s\mathrm{,}{I}_a\mathrm{,}R\right)=\text{In}\left(I_s\left(t\right)+I_a\left(t\right)+wR\left(t\right)\right)$, such that $w>0$. Using Itô formula we get

$\begin{array}{c}\text{d}V=\frac{1}{I_s\left(t\right)+I_a\left(t\right)+wR\left(t\right)}[\frac{p\beta BS}{\kappa +B}-\left(\mu +r_1+d\right)I_s\\ +\frac{q\beta BS}{\kappa +B}-\left(\mu +r_2\right)I_a+w\left(r_1{I}_s+r_2{I}_a-\mu R\right)]\text{d}t\\ -\frac{1}{2{\left(I_s\left(t\right)+I_a\left(t\right)+wR\left(t\right)\right)}^2}\left(\frac{{\sigma }_1^2{p}^2{B}^2{S}^2}{{\left(\kappa +B\right)}^2}+{\sigma }_2^2{I}_s^2+\frac{{\sigma }_1^2{q}^2{B}^2{S}^2}{{\left(\kappa +B\right)}^2}\right)\text{d}t\\ +\frac{{\sigma }_{1}pBS\text{d}{\mathcal{W}}_1\left(t\right)}{I_s\left(t\right)+I_a\left(t\right)+wR\left(t\right)\left(\kappa +B\right)}+\frac{{\sigma }_2{I}_s\text{d}{\mathcal{W}}_2\left(t\right)}{I_s\left(t\right)+I_a\left(t\right)+wR\left(t\right)}\\ +\frac{{\sigma }_{1}qBS\text{d}{\mathcal{W}}_1\left(t\right)}{\left(\kappa +B\right)I_s\left(t\right)+I_a\left(t\right)+wR\left(t\right)}\mathrm{.}\end{array}$ (40)

But $p+q=1$, also by dropping some of the negative terms in Equation (40) we obtain

$\begin{array}{c}\text{d}V\le \frac{1}{I_s\left(t\right)+I_a\left(t\right)+wR\left(t\right)}[\frac{\beta BS}{\kappa +B}-\left(\mu +r_1+d\right)I_s-\left(\mu +r_2\right)I_a\\ \underset{}{\overset{}{}}+w\left(r_1{I}_s+r_2{I}_a-\mu R\right)]\text{d}t+\frac{{\sigma }_{1}pBS\text{d}{\mathcal{W}}_1}{\left(\kappa +B\right)I_s\left(t\right)+I_a\left(t\right)+wR\left(t\right)}\\ +\frac{{\sigma }_2{I}_s\text{d}{\mathcal{W}}_2}{I_s\left(t\right)+I_a\left(t\right)+wR\left(t\right)}+\frac{{\sigma }_{1}qBS\text{d}{\mathcal{W}}_1\left(t\right)}{\left(\kappa +B\right)I_s\left(t\right)+I_a\left(t\right)+wR\left(t\right)}\mathrm{.}\end{array}$ (41)

$\begin{array}{c}\text{d}V\le \frac{1}{I_s\left(t\right)+I_a\left(t\right)+wR\left(t\right)}\left[\frac{\beta BS}{\kappa +B}-\left(\mu +r_1+d-w{r}_1\right)I_s\right]\text{d}t\\ +\frac{1}{I_s\left(t\right)+I_a\left(t\right)+wR\left(t\right)}\left[-\left(\mu +r_2-w{r}_2\right)I_a-w\mu R\right]\text{d}t\\ +\frac{{\sigma }_{1}pBS\text{d}{\mathcal{W}}_1}{\left(\kappa +B\right)I_s\left(t\right)+I_a\left(t\right)+wR\left(t\right)}+\frac{{\sigma }_2{I}_s\text{d}{\mathcal{W}}_2}{I_s\left(t\right)+I_a\left(t\right)+wR\left(t\right)}\\ +\frac{{\sigma }_{1}qBS\text{d}{\mathcal{W}}_1\left(t\right)}{\left(\kappa +B\right)I_s\left(t\right)+I_a\left(t\right)+wR\left(t\right)}\mathrm{.}\end{array}$ (42)

When $R_0<1$, then disease die out and hence $B=0$. Therefore

$\begin{array}{c}\text{d}V\le \frac{1}{I_s\left(t\right)+I_a\left(t\right)+wR\left(t\right)}\left[-\left(\mu +r_1+d-w{r}_1\right)I_s\right]\text{d}t\\ +\frac{1}{I_s\left(t\right)+I_a\left(t\right)+wR\left(t\right)}\left[-\left(\mu +r_2-w{r}_2\right)I_a-w\mu R\right]\text{d}t\\ +\frac{{\sigma }_2{I}_s\text{d}{\mathcal{W}}_2}{I_s\left(t\right)+I_a\left(t\right)+wR\left(t\right)}\mathrm{.}\end{array}$ (43)

Let $\theta =\min \left(\mu +r_1+d-w{r}_1,\mu +r_2-w{r}_2,w\mu \right)$. It follows that

$\text{d}V\le -\theta \text{d}t+\frac{{\sigma }_2{I}_s\text{d}{\mathcal{W}}_2}{I_s\left(t\right)+I_a\left(t\right)+wR\left(t\right)}\mathrm{.}$ (44)

Introducing integral both sides from 0 to t to Equation (44) we get

$\begin{array}{l}\text{In}\left(I_s\left(t\right)+I_a\left(t\right)+wR\left(t\right)\right)\\ \le \text{In}\left(I_s\left(0\right)+I_a\left(0\right)+wR\left(0\right)\right)-\theta t+{\displaystyle {\int }_0^t}\frac{{\sigma }_2{I}_s\left(s\right)\text{d}{\mathcal{W}}_2\left(s\right)}{I_s\left(s\right)+I_a\left(s\right)+wR\left(s\right)}\mathrm{,}\end{array}$ (45)

Then from the strong law of large number for local martingale we have

$\underset{t\to \infty }{\lim }\frac{1}{t}{\displaystyle {\int }_0^t}\frac{{\sigma }_2{I}_s\left(s\right)\text{d}{\mathcal{W}}_2\left(s\right)}{I_s\left(s\right)+I_a\left(s\right)+wR\left(s\right)}=0\text{a}\text{.s}\text{.}$ (46)

Hence, from Equation (45) and Equation (46) we have that

$\underset{t\to \infty }{\lim }\sup \frac{1}{t}\text{In}\left(I_s\left(t\right)+I_a\left(t\right)+wR\left(t\right)\right)\le -\theta <0.$

This completes the proof. □

Theorem 5. If $\frac{1}{2}{\sigma }^2>\frac{{\beta }^2}{\mu }$, then the disease free equilibrium $E^f$ is almost surely stable in $\Delta$.

Proof. Define a Lyapunov function as follows

$V=\text{In}\left[\left(\frac{\Lambda }{\mu }-S\right)+I_a+R+B\right]\mathrm{.}$

Applying the Itô formula to function V above we get

$\begin{array}{c}\text{d}V=\frac{1}{\left(\frac{\Lambda }{\mu }-S\right)+I_a+R+B}(-\Lambda +\frac{\beta BS}{\kappa +B}+\mu S+\frac{q\beta BS}{\kappa +B}-\left(\mu +r_2\right)I_a\\ \underset{}{\overset{}{}}+r_1{I}_s+r_2{I}_a-\mu R)+\frac{1}{\left(\frac{\Lambda }{\mu }-S\right)+I_a+R+B}\left({\alpha }_1{I}_s+{\alpha }_2{I}_a-\left(\delta +\varphi \right)\right)\text{d}t\\ -\frac{1}{2{\left[\left(\frac{\Lambda }{\mu }-S\right)+I_a+R+B\right]}^2}\left(\frac{{\sigma }_1^2{B}^2{S}^2}{{\left(\kappa +B\right)}^2}+\frac{{\sigma }_1^2{q}^2{B}^2{S}^2}{{\left(\kappa +B\right)}^2}\right)\text{d}t\\ +\frac{{\sigma }_{1}BS\text{d}{\mathcal{W}}_1\left(t\right)}{\left(\kappa +B\right)\left[\left(\frac{\Lambda }{\mu }-S\right)+I_a+R+B\right]}+\frac{{\sigma }_{1}qBS\text{d}{\mathcal{W}}_1\left(t\right)}{\left(\kappa +B\right)\left[\left(\frac{\Lambda }{\mu }-S\right)+I_a+R+B\right]},\end{array}$ (47)

From Equation (47) let $q=1$ and $p=0$, then

$\begin{array}{c}\text{d}V=\frac{1}{\left(\frac{\Lambda }{\mu }-S\right)+I_a+R+B}(-\Lambda +\frac{2\beta BS}{\kappa +B}-\mu \left(S-I_a-R\right)+r_1{I}_s+{\alpha }_1{I}_s\\ \underset{}{\overset{}{}}+{\alpha }_2{I}_a-\left(\delta +\varphi \right)B)\text{d}t-\frac{1}{{\left[\left(\frac{\Lambda }{\mu }-S\right)+I_a+R+B\right]}^2}\frac{2{\sigma }_1^2{B}^2{S}^2}{{\left(\kappa +B\right)}^2}\text{d}t\\ +\frac{2{\sigma }_{1}BS\text{d}{\mathcal{W}}_1\left(t\right)}{\left(\kappa +B\right)\left[\left(\frac{\Lambda }{\mu }-S\right)+I_a+R+B\right]}\mathrm{,}\end{array}$ (48)

$\begin{array}{c}\text{d}V=\frac{1}{\left(\frac{\Lambda }{\mu }-S\right)+I_a+R+B}(\frac{2\beta BS}{\kappa +B}-\Lambda -\mu \left(S-I_a-R\right)+r_1{I}_s+{\alpha }_1{I}_s\\ \underset{}{\overset{}{}}+{\alpha }_2{I}_a-\left(\delta +\varphi \right)B)\text{d}t-\frac{2{\sigma }_1^2{B}^2{S}^2}{{\left(\kappa +B\right)}^2{\left[\left(\frac{\Lambda }{\mu }-S\right)+I_a+R+B\right]}^2}\text{d}t\\ +\frac{2{\sigma }_{1}BS\text{d}{\mathcal{W}}_1\left(t\right)}{\left(\kappa +B\right)\left[\left(\frac{\Lambda }{\mu }-S\right)+I_a+R+B\right]}\mathrm{,}\end{array}$ (49)

let

$Y=\frac{BS}{\left(\kappa +B\right)\left[\left(\frac{\Lambda }{\mu }-S\right)+I_a+R+B\right]}.$

Then, $\text{d}V$ becomes

$\begin{array}{c}\text{d}V=\left(-2{\sigma }_1^2{Y}^2+2\beta Y-\frac{\Lambda +\mu S-\mu I_a-\mu R-r_1{I}_s-{\alpha }_1{I}_s-{\alpha }_2{I}_a+\left(\delta +\varphi \right)B}{\left(\frac{\Lambda }{\mu }-S\right)+I_a+R+B}\right)\text{d}t\\ +2{\sigma }_{1}Y\text{d}{\mathcal{W}}_1\left(t\right)\mathrm{,}\end{array}$ (50)

Then, $\text{d}V$ deduce to

$\text{d}V\le \left(-2{\sigma }_1^2{Y}^2+2\beta Y-\mu \right)\text{d}t+2{\sigma }_{1}Y\text{d}{\mathcal{W}}_1\left(t\right)\mathrm{,}$ (51)

since

$-2{\sigma }_1^2{Y}^2+2\beta Y-\mu =-2{\sigma }_1^2{\left(Y-\frac{\beta }{{\sigma }_1^2}\right)}^2+\frac{2{\beta }^2-\mu {\sigma }_1^2}{{\sigma }_1^2},$

we have

$\text{d}V\le \frac{2{\beta }^2-\mu {\sigma }_1^2}{{\sigma }_1^2}+2{\sigma }_{1}Y\text{d}{\mathcal{W}}_1\left(t\right)\mathrm{.}$ (52)

Applying integration from 0 to t both sides to Equation (52), we get

$\begin{array}{l}\text{In}\left(\left(\frac{\Lambda }{\mu }-S\left(t\right)\right)+I_a\left(t\right)+R\left(t\right)+B\left(t\right)\right)\\ \le \text{In}\left(\left(\frac{\Lambda }{\mu }-S\left(0\right)\right)+I_a\left(0\right)+R\left(0\right)+B\left(0\right)\right)\\ +\left(\frac{2{\beta }^2-\mu {\sigma }_1^2}{{\sigma }_1^2}\right)t+{\displaystyle {\int }_0^t}2{\sigma }_{1}Y\left(s\right)\text{d}{\mathcal{W}}_1\left(s\right)\mathrm{.}\end{array}$ (53)

Then, for all time $t\in \left[\mathrm{0,}T\right]$, the quadratic variation of the Itô integral process $I{\left(Y\right)}_t$ is deterministic integral over $\left[\mathrm{0,}t\right]$ of $Y_s^2$ i.e.,

${\left[{\displaystyle {\int }_0}Y\left(s\right)\text{d}{\mathcal{W}}_1\left(s\right)\right]}_t={\displaystyle {\int }_0^t}Y\left(s\right)\text{d}s\le Ct\mathrm{,}$ (54)

where C is a constant. Then by strong law of large number for local martingales, we have

$\underset{t\to \infty }{\lim }\frac{1}{t}{\displaystyle {\int }_0^t}Y\left(s\right)\text{d}{\mathcal{W}}_1\left(s\right)=0\text{a}\text{.s}\text{.}$ (55)

From Equation (54) and (55), we conclude that

$\underset{t\to \infty }{\lim }\sup \text{In}\left[\left(\frac{\Lambda }{\mu }-S\left(t\right)\right)+I_a\left(t\right)+R\left(t\right)+B\left(t\right)\right]\le \frac{2{\beta }^2-\mu {\sigma }_1^2}{{\sigma }_1^2}<0.$

This completes the proof 5. □

Theorem 6. If $R_0<1$, then $S\left(t\right),I_s\left(t\right),I_a\left(t\right),R\left(t\right),B\left(t\right)$ converge almost surely to

$\left(\frac{b{N}_e}{\mu }\mathrm{,0,0,0,0}\right)\text{in}\Delta \mathrm{.}$

Proof. We are required to show that

$\underset{t\to \infty }{\lim }\left(\frac{b{N}_e}{\mu }-S\left(t\right)\right)=0\text{a}\text{.s}\text{.}$

Consider the first equation of model system (3), let $\Lambda =b{N}_e$. Using Itô formula we have

$\text{d}\left(\frac{\Lambda }{\mu }-S\right)=\left(-\mu \left(\frac{\Lambda }{\mu }-S\right)+\frac{\beta BS}{\kappa +B}\right)\text{d}t+\frac{{\sigma }_{1}BS\text{d}{\mathcal{W}}_1\left(t\right)}{\kappa +B}\mathrm{,}$ (56)

Applying integration both sides from 0 to t to Equation (56), we get

$\frac{\Lambda }{\mu }-S\left(t\right)=\frac{\Lambda }{\mu }-S\left(0\right)+{\displaystyle {\int }_0^t}\frac{\beta BS\left(s\right)\text{d}s}{\kappa +B}-{\displaystyle {\int }_0^t}\mu \left(\frac{\Lambda }{\mu }-S\left(s\right)\right)\text{d}s+{\displaystyle {\int }_0^t}\frac{{\sigma }_1\beta BS\left(s\right)\text{d}{\mathcal{W}}_1\left(s\right)}{\kappa +B}\mathrm{,}$

From the fact that $\frac{B}{\kappa +B}\le 1$ and from Theorem 4, we have in $\Delta$ that

$\underset{t\to \infty }{\lim }{\displaystyle {\int }_0^t}\frac{\beta BS\left(s\right)\text{d}s}{\kappa +B}\le \underset{t\to \infty }{\lim }{\displaystyle {\int }_0^t}\beta S\left(s\right)\text{d}s\le \beta D_1{\text{e}}^{-D_{2}s}<\infty \mathrm{.}$ (57)

Hence

$\underset{t\to \infty }{\lim }\left(\frac{\Lambda }{\mu }-S\left(t\right)\right)<\infty \text{a}\text{.s}\text{.}$ (58)

Also

$\underset{t\to \infty }{\lim }{\displaystyle {\int }_0^t}\mu \left(\frac{\Lambda }{\mu }-S\left(s\right)\right)\text{d}s<\infty \text{a}\text{.s}\text{.}$ (59)

Combining Equation (57) and (59), we get

$\underset{t\to \infty }{\lim }{\displaystyle {\int }_0^t}\left(\frac{\Lambda }{\mu }-S\left(s\right)\right)\text{d}s={\displaystyle {\int }_0^{\infty }}\left(\frac{\Lambda }{\mu }-S\left(s\right)\right)\text{d}s<\infty \text{a}\text{.s}\text{.}$ (60)

If $S\left(t\right)$ does not converge almost surely to $\frac{\Lambda }{\mu }$, there is ${\Omega }_1\in \Omega$ such that the $\mathbb{P}\left({\Omega }_1\right)>0$, then for all $\omega \in \Omega$, it follows that

$\underset{t\to \infty }{\lim }\inf \left(\frac{\Lambda }{\mu }-S\left(t,\omega \right)\right)>0\text{a}\text{.s}\text{.}$

Then, there exists, $T>0$ such that

$\left(\frac{\Lambda }{\mu }-S\left(t,\omega \right)\right)>\frac{1}{2}\rho \left(\omega \right),\forall t\ge T.$

We have

$\begin{array}{c}{\displaystyle {\int }_0^{\infty }}\left(\frac{\Lambda }{\mu }-S\left(s,\omega \right)\right)\text{d}s={\displaystyle {\int }_0^T}\left(\frac{\Lambda }{\mu }-S\left(s,w\right)\right)\text{d}s+{\displaystyle {\int }_T^{\infty }}\left(\frac{\Lambda }{\mu }-S\left(s,w\right)\right)\text{d}s\\ \ge {\displaystyle {\int }_T^{\infty }}\left(\frac{\Lambda }{\mu }-S\left(s,w\right)\right)\text{d}s=\infty .\end{array}$ (61)

Then, as ${\Omega }_1\subset {\Omega }_2$, where

${\Omega }_2=\left\{\omega \mathrm{,}{\displaystyle {\int }_0^{\infty }}\left(\frac{\Lambda }{\mu }-S\left(s\mathrm{,}\omega \right)\right)\text{d}s=\infty \right\}\mathrm{.}$

It implies that $\mathbb{P}\left(\Omega \right)>0$. But from Equation (57), we see that $\mathbb{P}\left(\Omega \right)=0$. This leads to contradiction.

Hence,

$\underset{t\to \infty }{\lim }\left(\frac{\Lambda }{\mu }-S\left(t\right)\right)=0\text{a}\text{.s}\text{.}$

This completes the proof of Theorem 6. □

5. Numerical Results

In this section, we numerically solve the ordinary differential Equation (ODEs) model (1) and stochastic differential Equation (SDEs) model (3) and (8) by fourth order Runge-Kutta and Euler-Maruyama scheme respectively. The behaviour of individuals sample path of the stochastic differential equation models are compared to the deterministic solution. Initially, one infective is introduced into a population. The time step is $\Delta t=0.01$ and the time axis is the number of time steps, e.g., Time = 100 means 100 time steps and thus, an actual total time of $100\Delta t=1$. The Euler-Maruyama scheme is one of the numerical methods for computing the sample paths of SDEs (3) and (8) [21] . This method is a finite difference approximation

$x\left(t+\Delta t\right)=x\left(t\right)+f\left(x\left(t\right),t\right)\Delta t+G\left(x\left(t\right),t\right)\rho \sqrt{\Delta t},$ (62)

From Equation (62), the vector $\rho =\left({\rho }_1,\cdots ,{\rho }_j\right)$ is j independent standard normal random numbers ${\rho }_k\in N\left(\mathrm{0,1}\right)$. Also, $t=0,\Delta t,2\Delta t,\cdots$, and $\Delta t$ is chosen sufficiently very small to ensure good convergence. The sample path of the $S{I}_s{I}_{a}R-B$ stochastic differential equations is shown in Figure 6. It is observed that the sample path is continuous but not differentiable (a Wiener process property).

Maximum likelihood estimation method is used to estimate the unknown parameters $\theta$ of a SDE (8) by maximizing the likelihood [16] .

Let $p\left(t_k\mathrm{,}{x}_k\mathrm{<mo>\; |\; </mo>}{t}_{k-1}\mathrm{,}{x}_{k-1}\mathrm{<mo>\; ;\; </mo>}\theta \right)$ be the transition probability density of $\left(t_{k-1}\mathrm{,}{x}_{k-1}\right)$ given vector $\theta$. Suppose that the density of the initial state is $p_0\left(x_0\mathrm{<mo>\; |\; </mo>}\theta \right)$, in maximum likelihood estimation of $\theta$ [16] , the joint density

$D\left(\theta \right)=p_0\left(x_0|\theta \right){\displaystyle \underset{k=1}{\overset{N}{\prod }}}p\left(t_k,x_k|t_{k-1},x_{k-1};\theta \right),$ (63)

is maximized over $\theta \in {ℝ}^m$. The value of $\theta$ that maximizes $D\left(\theta \right)$ will be denoted by ${\theta }_{\mathrm{<mo>\; *\; </mo>}}\in {ℝ}^m$. For simplicity, it is more convenient to minimize the function

$L\left(\theta \right)=-\ln \left(p_0\left(x_0|\theta \right)\right)-{\displaystyle \underset{k=1}{\overset{N}{\sum }}}\ln \left(p\left(t_k,x_k|t_{k-1},x_{k-1};\theta \right)\right),$ (64)

where, $p\left(t_k\mathrm{,}{x}_k\mathrm{<mo>\; |\; </mo>}{t}_{k-1}\mathrm{,}{x}_{k-1}\mathrm{<mo>\; ;\; </mo>}\theta \right)$ is computed recursively by extended Kalman filter and its approximation is

$\eta \left(\theta \right)={\displaystyle \underset{j=1}{\overset{N}{\sum }}}\frac{1}{2}\ln \left|2\text{π}{S}_j\right|+\frac{1}{2}{\displaystyle \underset{j=1}{\overset{N}{\sum }}}{\left(y_j-h\left(m_j^{-},t\right)\right)}^{\text{T}}{S}_j^{-1}\left(y_j-h\left(m_j^{-},t\right)\right)-\ln p\left(\theta \right),$ (65)

where

$S_j=H_x\left(m_j^{-}\right)p_j^{-}{H}_x^{\text{T}}\left(m_j^{-}\right)+Z_j,$

$K_j=p_j^{-}{H}_x^{\text{T}}\left(m_j^{-},t\right)S_j^{-1},$ (66)

$m_j=m_j^{-}+K_j\left(y_j-h\left(m_j^{-}\right)\right),$

$p_j=p_j^{-}-K_j{S}_j{K}_j^{\text{T}},$

$H_x\left(x\mathrm{,}t\right)$ is the Jacobian matrix of $h\left(x\mathrm{,}t\right)$, $y_j=h\left(x\left(t_j\right)\right)+z_j$, $y_j$ is the measurement at time $t_j$, $x\left(t_j\right)$ is the state at time $t_j$, $\theta$ is the vector of parameter to estimated, $Z_j$ is the covariance matrix of the measurement error at $t_j$. However, at $t_0$ the state is assumed to have the prior distribution. For more details on EKF (see [22] ).

The estimated parameters by MLE are shown in Table 2. It is observed that the estimates are close to the true parameter values. Figures 2-5 show the extended Kalman filter estimates together with the true solution from ODE, measurements and states. Here black line, blue line and green line represents ODE solution, measurements and states realizations respectively. It is observed that the extended Kalman filter fits the states, which means states $\left(S\left(t\right),I_s\left(t\right),I_a\left(t\right),R\left(t\right)\right)$ can easily be estimated. The stochastic model (3) can be re-written in the following discretization equations:

$S_{k+1}=S_k+\left[b{N}_e-\frac{\beta B_k{S}_k}{\kappa +B_k}-\mu S_k\right]\Delta t-\frac{{\sigma }_1{B}_k{S}_k\sqrt{\Delta t}{Z}_k}{\kappa +B_k},$

$I_{sk+1}=I_{sk}+\left[\frac{p\beta B_k{S}_k}{\kappa +B_k}-\left(\mu +r_1+d\right)I_{sk}\right]\Delta t-\frac{{\sigma }_{1}p{B}_k{S}_k\sqrt{\Delta t}{Z}_k}{\kappa +B_k}-{\sigma }_2{I}_{sk}\sqrt{\Delta t}{Z}_k,$

$I_{ak+1}=I_{ak}+\left[\frac{q\beta B_k{S}_k}{\kappa +B_k}-\left(\mu +r_2\right)I_{ak}\right]\Delta t-\frac{{\sigma }_{1}q{B}_k{S}_k\sqrt{\Delta t}{Z}_k}{\kappa +B_k},$ (67)

Figure 2. The sample path of $S\left(t\right)$ for SDE (8) and deterministic solution for $S\left(t\right)$ graphed with extended Kalman filter estimates.

Figure 3. The sample path of $R\left(t\right)$ for SDE (8) and deterministic solution for $R\left(t\right)$ graphed with extended Kalman filter estimates.

Figure 4. The sample path of $I_a\left(t\right)$ for SDE (8) and deterministic solution for $I_a\left(t\right)$ graphed with extended Kalman filter estimates.

Figure 5. The sample path of $I_s\left(t\right)$ for SDE (8) and deterministic solution for $I_s\left(t\right)$ graphed with extended Kalman filter estimates.

$R_{k+1}=R_k+\left(r_1{I}_{sk}+r_2{I}_{ak}-\mu R_k\right)\Delta t,$

$B_{k+1}=B_k+\left({\alpha }_1{I}_{sk}+{\alpha }_2{I}_{ak}-\left(\delta +\varphi \right)B_k\right)\Delta t,$

where $Z_k\left(k=1,2,\cdots ,n\right)$ is the Gaussian random variables $N\left(\mathrm{0,1}\right)$. The Euler-Maruyama scheme is used to simulate the sample paths of stochastic differential Equation (3) and the result is graphed in Figure 6. From Figure 6, we observe that the susceptible proportion eventually converges to zero; the entire population becomes infected, and later they recover from the disease. Also, the sample path of $S{I}_s{I}_{a}R-B$ stochastic differential equations is continuous but not differentiable (a property of Wiener process).

The sample path of $S{I}_s{I}_{a}R-B$ stochastic differential equations model together with the solutions of ordinary differential equations is graphed in Figure 7. From Figure 7, we find that the sample path of $S{I}_s{I}_{a}R-B$ stochastic differential equations model fluctuates within the solution of the $S{I}_s{I}_{a}R-B$ ordinary differential equation model.

6. Discussion

We have proposed a new modeling framework for the dynamics of cholera using both deterministic and stochastic models. Our focus is on the interaction of environmental vibrios to human (which causes the transformation from the environmental vibrios to human) and the infected individuals shedding bacteria into the environment. For deterministic model, we derived the basic reproduction number $R_0$. The basic reproduction number is a critical parameter for disease dynamics. In the deterministic model, the value of the basic reproduction number $R_0$ determines the persistence or extinction of the disease. If $R_0<1$, the disease is eliminated, whereas if $R_0>1$, the disease persists in the population. From the deterministic model we have formulated two stochastic differential equations using parametric perturbation and transition probabilities methods.

Figure 6. The sample path of $S\left(t\right),I_s\left(t\right),I_a\left(t\right),R\left(t\right)$ for SDE model (3) when $\Delta t=0.01$, ${\sigma }_1=0.2$ and ${\sigma }_2=0.1$.

Figure 7. The sample path of $S\left(t\right),I_s\left(t\right),I_a\left(t\right),R\left(t\right)$ for SDE model (3) are graphed with the deterministic solutions (dashed curve) when $\Delta t=0.01$, ${\sigma }_1=0.2$ and ${\sigma }_2=0.1$.

We have proved the existence and uniqueness of positive solution, we showed that the solution of stochastic model are stochastically ultimately bounded, we derived that when $R_0<1$, then the infected compartments and bacteria goes to extinction. We carried out numerical simulation using Euler-Maruyama scheme to simulate the sample paths of stochastic differential Equation (3). Our results show that, the sample paths are continuous but not differentiable (a property of Wiener process) Also, we carefully compared the numerical simulation results for deterministic and stochastic models. We find that, the sample path of $S{I}_s{I}_{a}R-B$ stochastic differential equations model fluctuates within the solution of the $S{I}_s{I}_{a}R-B$ ordinary differential equation model as seen in Figure 7. However, the model parameters of SDEs are estimated by maximum likelihood estimation method. It is shown that the estimates are close to the true parameter values as seen in Table 2. Also, we used extended Kalman filter to estimate the states (compartments) of stochastic model (8) by recursively computing the transition probability density. It is observed that the state estimates fit the measurements as seen in Figures 2-5. Hence, we find that both models that are deterministic and stochastic models are very useful in understanding the dynamics of cholera epidemic. Nevertheless, Stochastic differential equation models are more important than deterministic models since they incorporate random effects such as environmental stochasticity and this enables us to model different quantities such as probability of extinction, probability of distributions and variances which cannot be captured in deterministic models.

7. Conclusions

In this paper, two stochastic differential equations models are formulated from the deterministic model using two different approaches: parametric perturbation and Transition probabilities. For deterministic model, the basic reproduction number $R_0$ determines whether the disease is eliminated or persists in the given population.

For stochastic model, the perturbed stochastic differential equation is first analyzed by proving the existence and positivity of the solutions. Secondly, we looked at the stability aspect of the model; we proved that the number of symptomatic infected, asymptomatic infected and bacteria tends to asymptotically to zero exponentially almost surely. Also, we showed that the equilibrium solution of the SDEs is pth moment exponentially stable and it is usually said to be exponentially stable in mean square. Numerical simulations are carried to simulate the sample paths of stochastic models by Euler-Maruyama scheme and the

Table 2. Estimated model parameters for cholera epidemic by MLE.

solutions of deterministic model by fourth order Runge-Kutta method. It is observed that the sample paths are continuous but not differentiable (a property of Wiener process) and the sample paths of stochastic differential equations models fluctuate within the solutions of deterministic model. However, the model parameters of SDEs are estimated by maximum likelihood estimation method. It is shown that the estimates are close to the true parameter values. Also, extended Kalman filter is used to estimate the states of stochastic model by recursively computing the transition probability density. It is observed that the state estimates fit the measurements. So, we can say that cholera transmission dynamics can be modeled using stochastic differential equations. It is clear that real world problems such as disease are not deterministic in nature so including random effects to the model gives us a more realistic way of modeling cholera epidemics and other epidemic diseases. For example, using stochastic differential equation model we managed to examine the limiting asymptotic distribution of the number of symptomatic infected, asymptomatic infected and bacteria.

Acknowledgements

The authors wish to thank Pan African University and Tanzania Public Service College for their financial support and reviewers for the insightful comments that made this manuscript better.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	World Health Organization (2016) World Health Statistics. Monitoring Health for the SDGs Sustainable Development Goals.
[2]	Roberts, M., Andreasen, V., Lloyd, A. and Pellis, L. (2015) Nine Challenges for Deterministic Epidemic Models. Epidemics, 10, 49-53. https://doi.org/10.1016/j.epidem.2014.09.006
[3]	Edward, S. and Nyerere, N. (2015) A Mathematical Model for the Dynamics of Cholera with Control Measures. Applied and Computational Mathematics, 4, 53-63. https://doi.org/10.11648/j.acm.20150402.14
[4]	King, A.A., Ionides, E.L., Pascual, M. and Bouma, M.J. (2008) Inapparent Infections and Cholera Dynamics. Nature, 454, 877-880. https://doi.org/10.1038/nature07084
[5]	Koelle, K., Rod, X., Pascual, M., Yunus, M. and Mostafa, G. (2005) Refractory Periods and Climate Forcing in Cholera Dynamics. Nature, 436, 696-700. https://doi.org/10.1038/nature03820
[6]	Obeng-Denteh, W., Andam, E.A., Obiri-Apraku, L. and Agyeil, W. (2015) Modeling Cholera Dynamics with a Control Strategy in Ghana. British Journal of Research, 2, 30-41.
[7]	Marwa, Y.M., Mwalili, S. and Mbalawata, I.S. (2018) Markov Chain Monte Carlo Analysis of Cholera Epidemic. Journal of Mathematical and Computational Science, 8, 584-610.
[8]	Dimi, J.L., Bissila, J.F. and Mbaya, T. (2016) Some Stochastic Properties of Cholera Model. Nonlinear Analysis and Differential Equations, 4, 779-792. https://doi.org/10.12988/nade.2016.61189
[9]	Gani, J. and Swift, R.J. (2009) Deterministic and Stochastic Models for the Spread of Cholera. The ANZIAM Journal, 51, 234-240. https://doi.org/10.1017/S1446181110000027
[10]	Sandro, A., Maritan, A., Bertuzzo, E., Rodriguez-Iturbe, I. and Rinaldo, A. (2010) Stochastic Dynamics of Cholera Epidemics. Physical Review E, 81, Article ID: 051901. https://doi.org/10.1103/PhysRevE.81.051901
[11]	Wang, X. and Wang, J. (2017) Modeling the Within-Host Dynamics of Cholera: Bacterial Viral Interaction. Journal of Biological Dynamics, 11, 484-501. https://doi.org/10.1080/17513758.2016.1269957
[12]	Cohn, J.N., Johnson, G., Ziesche, S., Cobb, F., Francis, G., Tristani, F. and Bhat, G. (1991) A Comparison of Enalapril with Hydralazine-Isosorbide Dinitrate in the Treatment of Chronic Congestive Heart Failure. New England Journal of Medicine, 325, 303-310. https://doi.org/10.1056/NEJM199108013250502
[13]	Van den Driessche, P. and Watmough, J. (2002) Reproduction Numbers and Sub-Threshold Endemic Equilibria for Compartmental Models of Disease Transmission. Mathematical Biosciences, 180, 29-48. https://doi.org/10.1016/S0025-5564(02)00108-6
[14]	Oksendal, B. (2003) Stochastic Differential Equations. Springer, Berlin, 65-84. https://doi.org/10.1007/978-3-642-14394-6_5
[15]	Rezaeyan, R. and Farnoosh, R. (2010) Stochastic Differential Equations and Application of the Kalman-Bucy Filter in the Modeling of RC Circuit. Applied Mathematical Sciences, 4, 1119-1127.
[16]	Allen, E. (2007) Modeling with Ito Stochastic Differential Equations. Springer Science, Business Media, Berlin, 22.
[17]	Allen, L.J. (2008) An Introduction to Stochastic Epidemic Models. In: Brauer, F., van den Driessche, P. and Wu, J., Eds., Mathematical Epidemiology, Springer, Berlin, 81-130. https://doi.org/10.1007/978-3-540-78911-6_3
[18]	Borodin, A.N. (2017) Stochastic Processes. Springer, Berlin. https://doi.org/10.1007/978-3-319-62310-8
[19]	Bahar, A. and Mao, X. (2004) Stochastic Delay Lotka-Volterra Model. Journal of Mathematical Analysis and Applications, 292, 364-380. https://doi.org/10.1016/j.jmaa.2003.12.004
[20]	Shaikhet, L. (1996) Stability of Stochastic Hereditary Systems with Markov Switching. Theory of Stochastic Processes, 2, 180-184.
[21]	Burrage, K., Burrage, P.M. and Tian, T. (2004) Numerical Methods for Strong Solutions of Stochastic Differential Equations: An Overview. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences, 460, 373-402. https://doi.org/10.1098/rspa.2003.1247
[22]	Mbalawata, I.S. (2014) Adaptive Markov Chain Monte Carlo and Bayesian Filtering for State Space Models. Doctoral Dissertation, Acta Universitatis Lappeenrantaensis.
[23]	Akman, O., Corby, M.R. and Schaefer, E. (2016) Examination of Models for Cholera: Insights into Model Comparison Methods. Letters in Biomathematics, 3, 93-118. https://doi.org/10.1080/23737867.2016.1211495
[24]	Khan, M.A., Ali, A., Dennis, L.C. and Gui, T. (2015) Dynamical Behavior of Cholera Epidemic Model with Non-Linear Incidence Rate. Applied Mathematical Sciences, 9, 989-1002. https://doi.org/10.12988/ams.2015.41166
[25]	Neilan, R.L., Schaefer, E., Gaff, H., Fister, K.R. and Lenhart, S. (2010) Modeling Optimal Intervention Strategies for Cholera. Bulletin of Mathematical Biology, 72, 2004-2018. https://doi.org/10.1007/s11538-010-9521-8
[26]	Njagarah, J.B. and Nyabadza, F. (2015) Modelling Optimal Control of Cholera in Communities Linked by Migration. Computational and Mathematical Methods in Medicine, 2015, Article ID: 898264. https://doi.org/10.1155/2015/898264