Dynamical Analysis of a Stochastic Cholera Epidemic Model

Abstract

Environmental disturbances have a strong impact on cholera transmission. Stochastic differential equations are an effective tool for characterizing environmental perturbations. In this paper, a stochastic infectious disease model for cholera is established and investigated. The dynamics of the stochastic cholera model are discussed. Firstly, the existence and uniqueness of the positive solution are proven. Then, the asymptotical stability of the disease-free equilibrium of the system is investigated. Furthermore, the asymptotical stability of the endemic equilibrium of the deterministic system corresponding to the stochastic system is obtained. Then, the theoretical results are verified by some numerical simulations. Finally, the optimal problem is considered as the theoretical basis for the control of cholera. Both theoretical and numerical results indicate that the random perturbations may make the model more realistic, which provides theoretical assessment for the control of cholera transmission.

1. Introduction

Cholera is an acute diarrheal infectious disease produced by Vibrio cholerae that is present in contaminated food or water [1]. Researchers have estimated that 1.3 to 4.0 million people become infected with cholera each year, and 21,000 to 143,000 of them die worldwide [2]. Many scholars have used mathematical models to study the spread of cholera. In [3], Zhou et al. considered an epidemic model for cholera with a periodic transmission rate. The threshold dynamics are established and investigated. In [4], Kazuo Yamazaki et al. studied the global well-posedness and asymptotic behavior of solutions to a reaction–convection–diffusion cholera epidemic model. In [5], Zhang et al. used a stochastic model to study the stationary distribution of a cholera epidemic model with vaccination under regime switching.

In [6], Liao and Wang considered the following cholera model:

(1)$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle \frac{d\mathcal{S}(t)}{dt}={\mu }_1-\frac{{\beta }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}-\frac{{\beta }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}-{\mu }_1\mathcal{S}(t),}\hfill \\ {\displaystyle \frac{d\mathcal{I}(t)}{dt}=\frac{{\beta }_1\mathcal{S}(t)B(t)}{1+{\alpha }_1\mathcal{B}(t)}+\frac{{\beta }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}-(r+\alpha +{\mu }_1)\mathcal{I}(t),}\hfill \\ {\displaystyle \frac{d\mathcal{R}(t)}{dt}=\alpha \mathcal{I}(t)-{\mu }_1\mathcal{R}(t),}\hfill \\ {\displaystyle \frac{d\mathcal{B}(t)}{dt}=\eta \mathcal{I}(t)-{\mu }_2\mathcal{B}(t).}\hfill \end{array}\right.\end{array}$

In system (1), $\mathcal{S}(t)$, $\mathcal{I}(t)$, $\mathcal{R}(t)$ represent the numbers of the susceptible individuals, infected individuals and recovered individuals at time t, respectively. It is assumed that the total population density is fixed, and $\mathcal{S}(t)+\mathcal{I}(t)+\mathcal{R}(t)=1.$ $\mathcal{B}(t)$ reflects the bacterial concentration at time t. All the parameters ${\beta }_1,{\beta }_2,{\mu }_1,{\mu }_2,\alpha ,r,\eta ,{\alpha }_1,{\alpha }_2$ are non-negative. The parameters ${\beta }_1$ and ${\beta }_2$ represent the contact rates for the human–environment and human–human interactions, respectively. The parameter ${\mu }_1$ denotes natural birth and mortality rates for humans; ${\mu }_2$ is the net mortality rate of the vibrios in the environment; $\alpha$ represents the natural recovery rate of the infected class; r is the disease-related mortality; $\eta$ indicates the human contribution rate to the growth of Vibrio cholerae. The constants ${\alpha }_1$ and ${\alpha }_2$ are the saturation factors, which measure the inhibitory effect. In [6], the authors investigated the globally asymptotic stability of the endemic equilibrium by using the Lyapunov functions method combined with the Volterra–Lyapunov stable matrices theory.

Numerous studies have shown that environmental disturbances have a huge impact on the spread of disease [7,8,9,10,11,12,13]. In [8], Zhou et al. proposed a stochastic SQEIAR model incorporating media reports and asymptomatic infection. They obtained the sufficient conditions for the existence of an ergodic stationary distribution for the system. A stochastic optimal control strategy for the system was presented. In [10], Ikram et al. established a stochastic delayed epidemic model of the COVID-19 outbreak. They discussed the existence of an ergodic stationary distribution by using the stochastic Lyapunov function. In [11], Zhang et al. studied a stochastic SIQS epidemic model with L$\acute{e}$vy jumps. In [13], Agarwal et al. formulated a stochastic model for the transmission dynamics of pine wilt disease. They developed some new results about the stochastic optimal control of the stochastic epidemic model. In order to explore the effect of environmental disturbances on cholera transmission, we introduce random effects into model (1) by replacing parameters ${\beta }_1$ and ${\beta }_2$ with ${\beta }_1\to {\beta }_1+{\sigma }_1\dot{w}(t)$ and ${\beta }_2\to {\beta }_2+{\sigma }_2\dot{w}(t)$, respectively. Here, $\dot{w}(t)$ is a standard Brownian motion in complete probability space $(\Omega ,{\{{\mathcal{F}}_t\}}_{t\ge 0},\mathbb{P})$. ${\sigma }_i(i=1,2)$ are the intensity of environmental white noise. Based on the above assumptions, we can build the following cholera epidemic model with random perturbations according to the literature [6]:

(2)$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle d\mathcal{S}(t)=[{\mu }_1-\frac{{\beta }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}-\frac{{\beta }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}-{\mu }_1\mathcal{S}(t)]dt}\hfill \\ {\displaystyle -\frac{{\sigma }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_{1}B(t)}dw(t)-\frac{{\sigma }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}dw(t),}\hfill \\ {\displaystyle d\mathcal{I}(t)=[\frac{{\beta }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}+\frac{{\beta }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}-(r+\alpha +{\mu }_1)\mathcal{I}(t)]dt}\hfill \\ {\displaystyle +\frac{{\sigma }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}dw(t)+\frac{{\sigma }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}dw(t),}\hfill \\ {\displaystyle d\mathcal{R}(t)=[\alpha \mathcal{I}(t)-{\mu }_1\mathcal{R}(t)]dt,}\hfill \\ {\displaystyle d\mathcal{B}(t)=[\eta \mathcal{I}(t)-{\mu }_2\mathcal{B}(t)]dt}\hfill \end{array}\right.\end{array}$

with initial condition

(3)$\begin{array}{c}\hfill \mathcal{S}(0)={\mathcal{S}}_0>0,\mathcal{I}(0)={\mathcal{I}}_0>0,\mathcal{R}(0)={\mathcal{R}}_0>0,\mathcal{B}(0)={\mathcal{B}}_0>0.\end{array}$

The meaning of the variables and parameters is the same as in system (1). When ${\sigma }_i=0(i=1,2)$, system (2) becomes the deterministic model (1).

The first, second, and fourth equations in system (2) are independent of $\mathcal{R}(t)$; in the rest of the paper, we will mainly study the following simplified model for convenience:

(4)$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle d\mathcal{S}(t)=[{\mu }_1-\frac{{\beta }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}-\frac{{\beta }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}-{\mu }_1\mathcal{S}(t)]dt}\hfill \\ {\displaystyle -\frac{{\sigma }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}dw(t)-\frac{{\sigma }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}dw(t),}\hfill \\ {\displaystyle d\mathcal{I}(t)=[\frac{{\beta }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}+\frac{{\beta }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}-(r+\alpha +{\mu }_1)\mathcal{I}(t)]dt}\hfill \\ {\displaystyle +\frac{{\sigma }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}dw(t)+\frac{{\sigma }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}dw(t),}\hfill \\ {\displaystyle d\mathcal{B}(t)=[\eta \mathcal{I}(t)-{\mu }_2\mathcal{B}(t)]dt.}\hfill \end{array}\right.\end{array}$

Throughout this paper, unless otherwise noted, let $(\Omega ,\mathcal{F},{\{{\mathcal{F}}_i\}}_{t\ge 0},\mathbb{P})$ be a complete probability space with a filtration ${\{{\mathcal{F}}_i\}}_{t\ge 0}$ satisfying the usual conditions (namely, it is increasing and right continuous while ${\mathcal{F}}_0$ contains all $\mathbb{P}$-null sets) (see [14]). Let $W(t)=(W_1(t),W_2(t),W_3(t),W_4(t))$ be a scalar Brownian motion defined on the complete probability space. Define

${\mathbb{R}}_{+}^d=\{x=(x_1,\cdots ,x_d)\in {\mathbb{R}}^d:x_i>0,1\le i\le d\}$

and

${\overline{\mathbb{R}}}_{+}^d=\{x=(x_1,\cdots ,x_d)\in {\mathbb{R}}^d:x_i\ge 0,1\le i\le d\}.$

Let us recall some basic theories of stochastic differential equations [14]. Consider the d-dimensional It$\widehat{o}$ process described by the stochastic differential equation

(5)$\begin{array}{c}\hfill dX(t)=f(X(t),t)dt+g(X(t),t)dW(t)fort\ge t_0,\end{array}$

with initial value $X(0)=X_0\in {\mathbb{R}}^d$, where $W(t)$ denotes a d-dimensional standard Brownian motion defined on the complete probability space $(\Omega ,\mathcal{F},{\{{\mathcal{F}}_i\}}_{t\ge 0},\mathbb{P})$. Denote by $C^{2,1}({\mathbb{R}}^d\times {\overline{\mathbb{R}}}^d;{\overline{\mathbb{R}}}^d)$ the family of all non-negative functions $V(X)$ defined on ${\mathbb{R}}^d\times {\overline{\mathbb{R}}}_{+}$ such that they are continuously twice differentiable in X. The differential operator L of Equation (5) is defined by Mao [14] as

$L=\frac{\partial }{\partial t}+\sum _{i=1}^d{f}_i(X,t)\frac{\partial }{\partial X_i}+\frac{1}{2}\sum _{i,j=1}^d{[g^{\top }(X,t)g(X,t)]}_{ij}\frac{{\partial }^2}{\partial X_i\partial X_j}.$

Let L act on a function $V\in C^{2,1}({\mathbb{R}}^d\times {\overline{\mathbb{R}}}^d;{\overline{\mathbb{R}}}^d)$, then

$LV(X,t)=V_t(X,t)+V_X(X,t)f(X,t)+\frac{1}{2}{g}^{\top }(X,t)V_{XX}g(X,t).$

Here, $V_t=\frac{\partial V}{\partial t}$, $V_X=(\frac{\partial V}{\partial X_1},\frac{\partial V}{\partial X_2},\cdots ,\frac{\partial V}{\partial X_d})$, $V_{XX}={(\frac{{\partial }^{2}V}{\partial X_i\partial X_j})}_{d\times d}$. In view of It$\widehat{o}$’s formula [14], if $X(t)\in {\mathbb{R}}^d$, then

$dV(X(t),t)=LV(X(t),t)dt+V_X(X(t),t)g(X(t),t)dW(t).$

The structure of the remainder of the paper is as follows. In Section 2, we give the existence and uniqueness of the positive solutions for system (4). We obtain the asymptotical stability of the disease-free equilibrium point of system (4) in Section 3. In Section 4, the asymptotic stability of the endemic equilibrium point of the corresponding deterministic model of the stochastic system (4) is examined. In Section 5, some numerical simulations are given to illustrate the theoretical results. In Section 6, the stochastic version of the optimal control problem is discussed. The article ends with a conclusion.

2. Existence and Uniqueness of the Positive Solution

Since variables $\mathcal{S}(t),\mathcal{I}(t),\mathcal{R}(t),\mathcal{B}(t)$ describe the population sizes, they should be non-negative. For this matter, we will prove that the solution of the system (4) is global and positive with the initial condition $\mathcal{S}(0)>0,\mathcal{I}(0)>0$ and $\mathcal{B}(0)>0$.

Theorem 1.

For any initial value${(\mathcal{S}(0),\mathcal{I}(0),\mathcal{B}(0))}^{\top }\in {\mathbb{R}}_{+}^3$, there is a positive solution$(\mathcal{S}(t),\mathcal{I}(t),\mathcal{B}(t))$of the stochastic model (4) for$t\ge 0$and the solution will be maintained in${\mathbb{R}}_{+}^3$with probability one, i.e., for$t\ge 0$, ${(\mathcal{S}(t),\mathcal{I}(t),\mathcal{B}(t))}^{\top }\in {\mathbb{R}}_{+}^3$almost surely (the abbreviation: a.s.).

Proof. 

Since the coefficients involved in the equations are locally Lipschitz continuous for the given initial population sizes ${(\mathcal{S}(0),\mathcal{I}(0),\mathcal{B}(0))}^{\top }\in {\mathbb{R}}_{+}^3$, by [15], for any initial value ${(\mathcal{S}(0),\mathcal{I}(0),\mathcal{B}(0))}^{⊺}\in {\mathbb{R}}_{+}^3$, system (4) has a unique local solution $(\mathcal{S}(t),\mathcal{I}(t),\mathcal{B}(t))$, $t\in [0,{\tau }_e]$, where ${\tau }_e$ is the explosion time. To show that actually this solution is global, we need to prove that ${\tau }_e=\infty$ a.s. Define the stopping time

${\tau }_{\rho }=inf\{t\in [0,{\tau }_{\rho }):(\mathcal{S}(t),\mathcal{I}(t),\mathcal{B}(t))\in (\frac{1}{\rho },\rho )\},$

Let $inf\varphi =\infty$ (whenever $\varphi$ denotes the empty set). Obviously, ${\tau }_k\le {\tau }_e$, if ${\tau }_k=\infty$ a.s., then ${\tau }_e=\infty$ a.s., and all the solutions of system (4) are positive for any $t\ge 0$.

If ${\tau }_k<\infty$, there exists a constant $T>0$ such that $\mathbb{P}({\tau }_k<T)>0$. Define function $V(X):{\mathbb{R}}_{+}^3\to {\mathbb{R}}_{+}$ as follows

$V={(\mathcal{S}+\mathcal{I}+\mathcal{B})}^2+\frac{1}{\mathcal{S}},$

where $X=(X_1,X_2,X_3)=(\mathcal{S},\mathcal{I},\mathcal{B}).$

For $t\in (0,{\tau }_m\wedge T)$ and $m\ge m_1$, by applying the Itô formula, we can obtain

(6)$\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{\partial V}{\partial X}f& =\frac{\partial V}{\partial \mathcal{S}}{f}_1+\frac{\partial V}{\partial \mathcal{I}}{f}_2+\frac{\partial V}{\partial \mathcal{B}}{f}_3\hfill \\ & =2(\mathcal{S}+\mathcal{I}+\mathcal{B})({\mu }_1-{\mu }_1\mathcal{S}-(r+\alpha +{\mu }_1)\mathcal{I}+\eta \mathcal{I}-{\mu }_1\mathcal{B})\hfill \\ & -\frac{1}{{\mathcal{S}}^2}({\mu }_1-\frac{{\beta }_1\mathcal{S}\mathcal{B}}{1+{\alpha }_1\mathcal{B}}-\frac{{\beta }_2\mathcal{S}\mathcal{I}}{1+{\alpha }_2\mathcal{I}}-{\mu }_1\mathcal{B}).\hfill \end{array}\end{array}$

Hence,

$\frac{1}{2}{g(t)}^{\mathrm{T}}{V}_{XX}g(t)=\frac{2{\sigma }_1^2{\mathcal{S}}^2{\mathcal{B}}^2}{{(1+{\alpha }_1\mathcal{B})}^2}+\frac{2{\sigma }_1^2{\mathcal{S}}^2{\mathcal{I}}^2}{{(1+{\alpha }_2\mathcal{I})}^2}+\frac{{\sigma }_2^2{\mathcal{S}}^2{\mathcal{I}}^2}{{\mathcal{S}}^3{(1+{\alpha }_2\mathcal{I})}^2}+\frac{2{\sigma }_1^2{\mathcal{S}}^2{B}^2}{{\mathcal{S}}^3{(1+{\alpha }_1\mathcal{B})}^2},$

and

(7)$\begin{array}{c}\hfill \begin{array}{cc}\hfill V_{X}g(t)& =(2(\mathcal{S}+\mathcal{I}+\mathcal{B})-\frac{1}{{\mathcal{S}}^2})(-\frac{{\sigma }_1\mathcal{S}\mathcal{B}}{1+{\alpha }_1\mathcal{B}})dw+(2(\mathcal{S}+\mathcal{I}+\mathcal{V})-\frac{1}{{\mathcal{S}}^2})(-\frac{{\sigma }_2\mathcal{S}\mathcal{I}}{1+{\alpha }_2\mathcal{I}})dw\hfill \\ & +2(\mathcal{S}+\mathcal{I}+\mathcal{B})\frac{{\sigma }_1\mathcal{S}\mathcal{B}}{1+{\alpha }_1\mathcal{B}}dw+2(\mathcal{S}+\mathcal{I}+\mathcal{B})\frac{{\sigma }_2\mathcal{S}\mathcal{I}}{1+{\alpha }_2\mathcal{I}}dw.\hfill \end{array}\end{array}$

From (6) and (7), we have

(8)$\begin{array}{c}\hfill \begin{array}{cc}\hfill dV(X(t))& =LV(X(t))dt+(2(\mathcal{S}+\mathcal{I}+\mathcal{B})-\frac{1}{{\mathcal{S}}^2}(-\frac{{\sigma }_1\mathcal{S}\mathcal{B}}{1+{\alpha }_1\mathcal{B}})dw\hfill \\ & +(2(\mathcal{S}+\mathcal{I}+\mathcal{B})-\frac{1}{{\mathcal{S}}^2})(-\frac{{\sigma }_2\mathcal{S}\mathcal{I}}{1+{\alpha }_2\mathcal{I}})dw\hfill \\ & +2(\mathcal{S}+\mathcal{I}+\mathcal{B})\frac{{\sigma }_1\mathcal{S}\mathcal{B}}{1+{\alpha }_1\mathcal{B}}dw+2(\mathcal{S}+\mathcal{I}+\mathcal{B})\frac{{\sigma }_2\mathcal{S}\mathcal{I}}{1+{\alpha }_2\mathcal{I}}dw,\hfill \end{array}\end{array}$

where $LV:{\mathbb{R}}_{+}^3\to {\mathbb{R}}_{+}$ is defined by

$\begin{array}{c}\hfill \begin{array}{cc}\hfill LV& =2(\mathcal{S}+\mathcal{I}+\mathcal{B})({\mu }_1-{\mu }_1\mathcal{S}-(d+\alpha +{\mu }_1)\mathcal{I}+\eta \mathcal{I}-{\mu }_1\mathcal{B})\hfill \\ & -\frac{1}{{\mathcal{S}}^2}({\mu }_1-\frac{{\beta }_1\mathcal{S}\mathcal{B}}{1+{\alpha }_1\mathcal{B}}-\frac{{\beta }_2\mathcal{S}\mathcal{I}}{1+{\alpha }_2\mathcal{I}}-{\mu }_1\mathcal{S})+\frac{2{\sigma }_1^2{\mathcal{S}}^2{\mathcal{B}}^2}{{(1+{\alpha }_1\mathcal{B})}^2}+\frac{2{\sigma }_1^2{\mathcal{S}}^2{\mathcal{I}}^2}{{(1+{\alpha }_2\mathcal{I})}^2}\hfill \\ & +\frac{{\sigma }_2^2{\mathcal{S}}^2{\mathcal{I}}^2}{{\mathcal{S}}^3{(1+{\alpha }_2\mathcal{I})}^2}+\frac{2{\sigma }_1^2{\mathcal{S}}^2{\mathcal{B}}^2}{{\mathcal{S}}^3{(1+{\alpha }_1\mathcal{B})}^2}.\hfill \end{array}\end{array}$

Because ${(a-b)}^2\ge 0$, we can obtain

$\begin{array}{c}\hfill \begin{array}{cc}\hfill LV(X)& \le {\mu }_1^2+{(\mathcal{S}+\mathcal{I}+\mathcal{B})}^2+2(\mathcal{S}+\mathcal{I}+\mathcal{B})\eta \mathcal{I}+2{\sigma }_1^2{\mathcal{S}}^2{\mathcal{B}}^2+2{\sigma }_2^2{\mathcal{S}}^2{\mathcal{I}}^2\hfill \\ & +\frac{2{\sigma }_1^2}{\mathcal{S}{\alpha }_1^2}+\frac{2{\sigma }_2^2}{\mathcal{S}{\alpha }_2^2}+\frac{{\beta }_1}{\mathcal{S}{\alpha }_1}+\frac{{\beta }_2}{\mathcal{S}{\alpha }_2}+\frac{{\mu }_1}{\mathcal{S}}\hfill \\ & ={\mu }_1^2+2{\sigma }_1^2\frac{{\eta }^2}{{\mu }_2^2}+2{\sigma }_2^2+\frac{2{\sigma }_1^2}{S{\alpha }_1^2}+\frac{2{\sigma }_2^2}{S{\alpha }_2^2}+\frac{{\beta }_1}{\mathcal{S}{\alpha }_1}+\frac{{\beta }_2}{\mathcal{S}{\alpha }_2}+{\mathcal{S}}^2+{\mathcal{I}}^2(1+2\eta )\hfill \\ & +{\mathcal{B}}^2+2\mathcal{S}\mathcal{I}(1+\eta )+2\mathcal{B}\mathcal{I}(1+\eta )+2\mathcal{B}\mathcal{S}\hfill \\ & \le {\mu }_1^2+2{\sigma }_1^2\frac{{\eta }^2}{{\mu }_2^2}+2{\sigma }_2^2+\frac{2{\sigma }_1^2}{S{\alpha }_1^2}+\frac{2{\sigma }_2^2}{S{\alpha }_2^2}+\frac{{\beta }_1}{S{\alpha }_1}+\frac{{\beta }_2}{S{\alpha }_2}+C_1{(\mathcal{S}+\mathcal{I}+\mathcal{B})}^2\hfill \\ & \le {\mu }_1^2+2{\sigma }_1^2\frac{{\eta }^2}{{\mu }_2^2}+2{\sigma }_2^2+C_2({(\mathcal{S}+\mathcal{I}+\mathcal{B})}^2+\frac{1}{\mathcal{S}})\hfill \\ & \le C(1+V),\hfill \end{array}\end{array}$

where

$C_1=max\{1,1+2\eta \}=1+2\eta ,$

$C_2=max\{\frac{2{\sigma }_1^2}{{\alpha }_1^2}+\frac{2{\sigma }_2^2}{{\alpha }_2^2}+\frac{{\beta }_1}{{\alpha }_1}+\frac{{\beta }_2}{{\alpha }_2}+{\mu }_1,C_1\},$

$C=max\{{\mu }_1^2+2{\sigma }_1^2\frac{{\eta }^2}{{\mu }_2^2}+2{\sigma }_2^2,C_2,C_1\}.$

Integrating both sides of (8) from 0 to $(t_{\rho }\wedge t)$, then taking expectation

$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathbb{E}V(X(t\wedge {\tau }_{\rho }))& =V(X(0))+\mathbb{E}{\int }_0^{t\wedge {\tau }_{\rho }}LV(s)dr\hfill \\ & \le V(X(0))+\mathbb{E}{\int }_0^{t\wedge {\tau }_{\rho }}C(1+V(X(r)))dr\hfill \\ & \le V(X(0))+CT++\mathbb{E}{\int }_0^{t\wedge {\tau }_{\rho }}CV(X(r))ds\hfill \\ & \le V(X(0))+CT+\mathbb{E}{\int }_0^{t}CV(X(r\wedge {\tau }_{\rho }))dr.\hfill \end{array}\end{array}$

By using Gronwall inequality [16], we have

$V(X(T\wedge {\tau }_{\rho }))\le V(X(0))+CT)e^{CT}.$

Obviously, for every $\omega \in {\tau }_{\rho }\le T$, this function $X({\tau }_{\rho },\omega )$ equals $\rho$ or $\frac{1}{\rho }$. Hence,

(9)$\begin{array}{c}\hfill (9{\rho }^2+\frac{1}{\rho })\wedge (\frac{9}{{\rho }^2}+\rho )\le V(X({\tau }^{\rho })).\end{array}$

From (8) and (9), we have

$(V(X(0))+CT)e^{CT}\ge \mathbb{E}[I_{{\tau }_{\rho }\le T(\omega )}V(X({\tau }_{\rho }))].$

When $\rho \to +\infty$, we obtain

$\underset{\rho \to \infty }{lim}\mathbb{P}({\tau }_{\rho }\le T)\le \underset{\rho \to \infty }{lim}\frac{(V(X(0)+CT)e^{CT}}{(9{\rho }^2+\frac{1}{\rho })(\frac{9}{{\rho }^2}+\rho )}.$

Hence,

$\mathbb{P}({\tau }_{\infty }\le T)=0.$

Since we assume that $T\in (0,+\infty )$ and T is a positive real number, $\mathbb{P}({\tau }_{\infty }=\infty )=1$.   □

Remark 1.

By Theorem 1, for $t\le {\tau }_{\rho }$, from the first three equations of system (4), we can obtain

$\begin{array}{cc}d(\mathcal{S}(t)+\mathcal{I}(t)+\mathcal{R}(t))\hfill & =[{\mu }_1-{\mu }_1(\mathcal{S}(t)+\mathcal{I}(t)+\mathcal{R}(t))-r\mathcal{I}(t)]dt\hfill \\ & \le [{\mu }_1-{\mu }_1(\mathcal{S}(t)+\mathcal{I}(t)+\mathcal{R}(t))]dt.\hfill \end{array}$

Then,

$\begin{array}{cc}\mathcal{S}(t)+\mathcal{I}(t)+\mathcal{R}(t)\hfill & \le 1+e^{-{\mu }_{1}t}[\mathcal{S}(0)+\mathcal{I}(0)+\mathcal{R}(0)-1]\hfill \\ & \le \left\{\begin{array}{c}1,if\mathcal{S}(0)+\mathcal{I}(0)+\mathcal{R}(0)\le 1\hfill \\ \mathcal{S}(0)+\mathcal{I}(0)+\mathcal{R}(0),if\mathcal{S}(0)+\mathcal{I}(0)+\mathcal{R}(0)>1.\hfill \end{array}\right.\hfill \end{array}$

Similarly, from the last equation of system (4), we can obtain

$d\mathcal{B}(t)\le [\eta -{\mu }_2\mathcal{B}(t)]dt.$

Then,

$\mathcal{B}(t)\le max\{\frac{\eta }{{\mu }_2},\mathcal{B}(0)\}.$

Hence,

$\mathcal{S}(t)+\mathcal{I}(t)+\mathcal{R}(t)\le 1,\mathcal{B}(t)\le \frac{\eta }{{\mu }_2}.$

Thus,

$\Gamma =\{(\mathcal{S}(t),\mathcal{I}(t),\mathcal{R}(t),\mathcal{B}(t))\in {\mathbb{R}}_{+}^4|\mathcal{S}(t)+\mathcal{I}(t)+\mathcal{R}(t)\le 1,\mathcal{B}(t)\le \frac{\eta }{{\mu }_2}\}$

is the positive invariant set of stochastic system (4).

3. Stochastic Disease-Free Dynamics

In [6], Liao et al. presented the basic reproduction number ${\mathfrak{R}}_0$ of the deterministic system corresponding to system (2) as follows:

${\mathfrak{R}}_0=\frac{{\beta }_2{\mu }_2+{\beta }_1\eta }{{\mu }_2(r+\alpha +{\mu }_1)}.$

Furthermore, they obtained the local and global stability of the disease-free equilibrium and the endemic equilibrium associated with the threshold ${\mathfrak{R}}_0$, which determines whether a disease becomes extinct or not.

If ${\mathfrak{R}}_0<1$, the deterministic system (1) has a globally asymptotically stable disease-free equilibrium ${\mathcal{E}}^0=(1,0,0)$. Obviously, ${\mathcal{E}}^0=(1,0,0)$ is also the disease-free equilibrium of system (4). In the following, we will discuss the asymptotical stability of the disease-free equilibrium of system (4).

Remark 2.

Assume that $(\mathcal{S}(t),\mathcal{I}(t),\mathcal{B}(t))$ is the solution of system (4), which satisfies the initial condition $(\mathcal{S}(0),\mathcal{I}(0),\mathcal{B}(0))\in \Gamma$. If ${\displaystyle {\mathfrak{R}}_0=\frac{{\beta }_2{\mu }_2+{\beta }_1\eta }{{\mu }_2(r+\alpha +{\mu }_1)}<1}$, then

$\frac{1}{t}\underset{t\to \infty }{lim}\frac{1}{t}\mathbb{E}{\int }_0^t[(1-{\mathfrak{R}}_0)\mathcal{I}+{\mu }_1\mathcal{S}{(1-\frac{1}{\mathcal{S}})}^2]dr\le \frac{1}{2}\frac{{\sigma }_1^2}{{\alpha }_1^2}+\frac{1}{2}\frac{{\sigma }_2^2}{{\alpha }_2^2}.$

In other words, when the random perturbation intensity ${\sigma }_1$ and ${\sigma }_2$ are small, the number of infected individuals exponentially approaches zero a.s., i.e., cholera becomes extinct with probability 1.

Proof. 

In the following demonstration, we use the method of reference [17,18]. Constructing function $V:{\mathbb{R}}_{+}^3\to {\mathbb{R}}_{+}$ as follows,

$V_1(\mathcal{S},\mathcal{I},\mathcal{B})=\mathcal{S}-1-ln\mathcal{S}+\mathcal{I}+\frac{{\beta }_1}{{\mu }_2}\mathcal{B}.$

Applying the Itô formula, we arrive at

$\begin{array}{c}\hfill \begin{array}{cc}\hfill V_{1_X}f=& (1-\frac{1}{\mathcal{S}})f_1+f_2+\frac{{\beta }_1}{{\mu }_2}{f}_3\hfill \\ \hfill =& (1-\frac{1}{\mathcal{S}})({\mu }_1-\frac{{\beta }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}-\frac{{\beta }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}-{\mu }_1\mathcal{S}(t))+(\frac{{\beta }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}\hfill \\ & +\frac{{\beta }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}-(r+\alpha +{\mu }_1)\mathcal{I}(t))+\frac{{\beta }_1}{{\mu }_2}(\eta \mathcal{I}-{\mu }_2\mathcal{B})\hfill \\ \hfill =& -{\mu }_1\mathcal{S}{(1-\frac{1}{\mathcal{S}})}^2-(1-\frac{1}{\mathcal{S}})(\frac{{\beta }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}+\frac{{\beta }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)})+\frac{{\beta }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}\hfill \\ & +\frac{{\beta }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}-(r+\alpha +{\mu }_1)\mathcal{I}(t)+\frac{{\beta }_1}{{\mu }_2}\eta \mathcal{I}-\frac{{\beta }_1}{{\mu }_2}{\mu }_2\mathcal{B}\hfill \\ \hfill =& -{\mu }_1\mathcal{S}{(1-\frac{1}{\mathcal{S}})}^2+\frac{{\beta }_1\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}+\frac{{\beta }_2\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}-(r+\alpha +{\mu }_1)\mathcal{I}(t)\hfill \\ & +\frac{{\beta }_1}{{\mu }_2}\eta \mathcal{I}-\frac{{\beta }_1}{{\mu }_2}{\mu }_2\mathcal{B}\hfill \\ \hfill =& -{\mu }_1\mathcal{S}{(1-\frac{1}{\mathcal{S}})}^2+(\frac{{\beta }_1}{1+{\alpha }_1\mathcal{B}(t)}-{\beta }_1)\mathcal{B}+[{\beta }_2+\frac{\eta {\beta }_1}{{\mu }_2}-(r+\alpha +{\mu }_1)]\mathcal{I}\hfill \\ \hfill \le & -{\mu }_1\mathcal{S}{(1-\frac{1}{\mathcal{S}})}^2+(r+\alpha +{\mu }_1)[\frac{{\beta }_2{\mu }_2+\eta {\beta }_1}{{\mu }_2(r+\alpha +{\mu }_1)}-1]\mathcal{I}.\hfill \end{array}\end{array}$

Simplifying the above expression, we have

(10)$\begin{array}{c}\hfill V_{1_X}f⩽(r+\alpha +{\mu }_1)[{\mathfrak{R}}_0-1]\mathcal{I}-{\mu }_1\mathcal{S}{(1-\frac{1}{\mathcal{S}})}^2.\end{array}$

In addition,

(11)$\begin{array}{c}\hfill \frac{1}{2}\mathrm{trance}(g^{\mathrm{T}}{V}_{1_{XX}}g)=\frac{1}{2}\frac{{\sigma }_1^2{\mathcal{B}}^2}{{(1+{\alpha }_1\mathcal{B})}^2}+\frac{1}{2}\frac{{\sigma }_2^2{\mathcal{I}}^2}{{(1+{\alpha }_2\mathcal{I})}^2}\end{array}$

and

(12)$\begin{array}{cc}\hfill V_{1_X}gdw=& (1-\frac{1}{\mathcal{S}})(-\frac{{\sigma }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)})dw+(1-\frac{1}{\mathcal{S}})(-\frac{{\sigma }_1\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)})dw\hfill \\ \hfill & +\frac{{\sigma }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}dw-\frac{{\sigma }_1\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)})dw.\hfill \end{array}$

From (10)–(12), we obtain

$\begin{array}{c}\hfill \begin{array}{c}\hfill d{V}_1=(r+\alpha +{\mu }_1)[{\mathfrak{R}}_0-1]\mathcal{I}-{\mu }_1\mathcal{S}{(1-\frac{1}{\mathcal{S}})}^2+\frac{1}{2}\frac{{\sigma }_1^2{\mathcal{B}}^2}{{(1+{\alpha }_1\mathcal{B})}^2}+\frac{1}{2}\frac{{\sigma }_2^2{\mathcal{I}}^2}{{(1+{\alpha }_2\mathcal{I})}^2}.\end{array}\end{array}$

Integrating both sides of $d{V}_1$ from 0 to t, then taking the expectation, we have

$\begin{array}{c}\hfill \begin{array}{cc}\hfill 0\le & \mathbb{E}(V_1(t))-V_1(0)=\mathbb{E}{\int }_0^{t}L{V}_1(s)ds\hfill \\ & \le -\mathbb{E}{\int }_0^t[(1-{\mathfrak{R}}_0)\mathcal{I}+{\mu }_1\mathcal{S}{(1-\frac{1}{\mathcal{S}})}^2]dr+\frac{1}{2}\frac{{\sigma }_1^2{\mathcal{B}}^2}{{(1+{\alpha }_1\mathcal{B})}^2}+\frac{1}{2}\frac{{\sigma }_2^2{\mathcal{I}}^2}{{(1+{\alpha }_2\mathcal{I})}^2}.\hfill \end{array}\end{array}$

Hence,

$\underset{t\to \infty }{lim}\frac{1}{t}\mathbb{E}{\int }_0^t[(1-{\mathfrak{R}}_0)\mathcal{I}+{\mu }_1\mathcal{S}{(1-\frac{1}{\mathcal{S}})}^2]dr\le \frac{1}{2}\frac{{\sigma }_1^2}{{\alpha }_1^2}+\frac{1}{2}\frac{{\sigma }_2^2}{{\alpha }_2^2},$

which completes the proof.   □

4. Stochastic Endemic Dynamics

If ${\mathfrak{R}}_0>1$, there is a globally asymptotically stable endemic equilibrium ${\mathcal{E}}^{\ast }({\mathcal{S}}^{\ast },{\mathcal{I}}^{\ast },{\mathcal{B}}^{\ast })$ of the deterministic system (1), where the expressions of ${\mathcal{S}}^{\ast },{\mathcal{I}}^{\ast },{\mathcal{B}}^{\ast }$ can be found in [6]. Obviously, there is no endemic equilibrium in system (4), which is derived from the parametric perturbations of its deterministic model. Therefore, we will discuss the asymptotic behavior of the endemic equilibrium of the system (4) at the place of the corresponding deterministic system (1), which can disclose whether the disease is prevalent to some extent.

Theorem 2.

If ${\mathfrak{R}}_0>1$, there exists an appropriate positive constant ρ, for any initial value $(\mathcal{S}(0),\mathcal{I}(0),\mathcal{B}(0))\in {\mathbb{R}}_{+}^3$ and $\mathcal{S}(0)+\mathcal{I}(0)\le 1$, and the solution $(\mathcal{S}(t),\mathcal{I}(t),\mathcal{B}(t))$ of system (4) has the following property

$\begin{array}{c}\hfill \begin{array}{cc}\hfill \underset{t\to \infty }{lim}\frac{1}{t}\mathbb{E}& {\int }_0^t[(c_1{\mu }_1+{\mu }_2+\frac{\rho }{2}({\mu }_1+{\mu }_2)){({\mathcal{S}}^{\ast }-\mathcal{S})}^2+(c_1\frac{k}{{\mathcal{I}}^{\ast }}-\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}+\frac{k}{{\mathcal{I}}^{\ast }}){({\mathcal{I}}^{\ast }-\mathcal{I})}^2\hfill \\ & +(c_2{\mu }_2+{\mu }_2+\frac{1}{2\rho }({\mu }_1+{\mu }_2)){(\mathcal{B}-{\mathcal{B}}^{\ast })}^2)]dr\hfill \\ & \le \frac{{\sigma }_1^2}{{\alpha }_1^2}(1+c_1)+\frac{{\sigma }_2^2}{{\alpha }_2^2}(1+c_2),\hfill \end{array}\end{array}$

where

$c_1=\frac{\frac{{\mu }_1{\mathcal{B}}^{\ast }}{I^{\ast }}-\frac{k}{{\mathcal{I}}^{\ast }-{\mu }_1}}{\frac{k}{{\mathcal{I}}^{\ast }}+{\mu }_1},c_2=\frac{\frac{{\mu }_2{\mathcal{B}}^{\ast }}{I^{\ast }}-\frac{k}{{\mathcal{I}}^{\ast }-{\mu }_2}}{\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}},k=\frac{{\beta }_1{\mathcal{S}}^{\ast }{\mathcal{B}}^{\ast }}{1+{\alpha }_1{\mathcal{B}}^{\ast }}+\frac{{\beta }_2{\mathcal{S}}^{\ast }{\mathcal{I}}^{\ast }}{1+{\alpha }_2{\mathcal{I}}^{\ast }}.$

Proof. 

At equilibrium ${\mathcal{E}}^{\ast }({\mathcal{S}}^{\ast },{\mathcal{I}}^{\ast },{\mathcal{B}}^{\ast })$, the following equations are satisfied:

$\begin{array}{c}\hfill \begin{array}{cc}& {\mu }_1=\frac{{\beta }_1{\mathcal{S}}^{\ast }{\mathcal{B}}^{\ast }}{1+{\alpha }_1{\mathcal{B}}^{\ast }}+\frac{{\beta }_2{\mathcal{S}}^{\ast }{\mathcal{I}}^{\ast }}{1+{\alpha }_2{\mathcal{I}}^{\ast }}+{\mu }_1{\mathcal{S}}^{\ast },\hfill \\ & r+\alpha +{\mu }_1=\frac{{\beta }_1{\mathcal{S}}^{\ast }{\mathcal{B}}^{\ast }}{(1+{\alpha }_1{\mathcal{B}}^{\ast }){\mathcal{I}}^{\ast }}+\frac{{\beta }_2{\mathcal{S}}^{\ast }}{1+{\alpha }_2{\mathcal{I}}^{\ast }},\hfill \\ & \eta =\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}.\hfill \end{array}\end{array}$

Define a $C^2$-function:

$V_2(\mathcal{S},\mathcal{I},\mathcal{B})=\frac{c_1}{2}{V}_{21}+\frac{c_2}{2}{V}_{22}+V_{23},$

where

$V_{21}=\frac{{(\mathcal{S}-{\mathcal{S}}^{\ast }+\mathcal{I}-{\mathcal{I}}^{\ast })}^2}{2},V_{22}={(\mathcal{B}-{\mathcal{B}}^{\ast })}^2,V_{23}=\frac{{(\mathcal{S}-{\mathcal{S}}^{\ast }+\mathcal{I}-{\mathcal{I}}^{\ast }+\mathcal{B}-{\mathcal{B}}^{\ast })}^2}{2}.$

Hence,

$\begin{array}{c}\hfill \begin{array}{cc}\hfill L{V}_{21}=& (\mathcal{S}-{\mathcal{S}}^{\ast }+\mathcal{I}-{\mathcal{I}}^{\ast })(\frac{{\beta }_1{\mathcal{S}}^{\ast }{\mathcal{B}}^{\ast }}{1+{\alpha }_1{\mathcal{B}}^{\ast }}+\frac{{\beta }_2{\mathcal{S}}^{\ast }{\mathcal{I}}^{\ast }}{1+{\alpha }_2{\mathcal{I}}^{\ast }}+{\mu }_1{\mathcal{S}}^{\ast }-\frac{{\beta }_1\mathcal{S}\mathcal{B}}{1+{\alpha }_1\mathcal{B}}-\frac{{\beta }_2\mathcal{S}\mathcal{I}}{1+{\alpha }_2\mathcal{I}}-{\mu }_1\mathcal{S})\hfill \\ & +(\mathcal{S}-{\mathcal{S}}^{\ast }+\mathcal{I}-{\mathcal{I}}^{\ast })(\frac{{\beta }_1\mathcal{S}\mathcal{B}}{1+{\alpha }_1\mathcal{B}}+\frac{{\beta }_2\mathcal{S}\mathcal{I}}{1+{\alpha }_2\mathcal{I}}-\frac{{\beta }_1{\mathcal{S}}^{\ast }{\mathcal{B}}^{\ast }}{1+{\alpha }_1{\mathcal{B}}^{\ast }}\frac{\mathcal{I}}{{\mathcal{I}}^{\ast }}-\frac{{\beta }_2{\mathcal{S}}^{\ast }{\mathcal{I}}^{\ast }}{1+{\alpha }_2{\mathcal{I}}^{\ast }}\frac{\mathcal{I}}{{\mathcal{I}}^{\ast }})\hfill \\ & +\frac{{\sigma }_1^2{\mathcal{S}}^2{\mathcal{B}}^2)}{{(1+{\alpha }_1\mathcal{B})}^2}+\frac{{\sigma }_2^2{\mathcal{S}}^2{\mathcal{I}}^2)}{{(1+{\alpha }_2\mathcal{I})}^2}\hfill \\ \hfill \le & (\frac{{\beta }_1{\mathcal{S}}^{\ast }{\mathcal{B}}^{\ast }}{1+{\alpha }_1{\mathcal{B}}^{\ast }}+\frac{{\beta }_2{\mathcal{S}}^{\ast }{\mathcal{I}}^{\ast }}{1+{\alpha }_2{\mathcal{I}}^{\ast }})(\mathcal{S}-{\mathcal{S}}^{\ast }+\mathcal{I}-{\mathcal{I}}^{\ast })(1-\frac{\mathcal{I}}{{\mathcal{I}}^{\ast }})\hfill \\ & -{\mu }_1(\mathcal{S}-{\mathcal{S}}^{\ast })(\mathcal{I}-{\mathcal{I}}^{\ast })+\frac{{\sigma }_1^2}{{\alpha }_1^2}+\frac{{\sigma }_2^2}{{\alpha }_2^2}.\hfill \end{array}\end{array}$

Let ${\displaystyle k=\frac{{\beta }_1{\mathcal{S}}^{\ast }{\mathcal{B}}^{\ast }}{1+{\alpha }_1{\mathcal{B}}^{\ast }}+\frac{{\beta }_2{\mathcal{S}}^{\ast }{\mathcal{I}}^{\ast }}{1+{\alpha }_2{\mathcal{I}}^{\ast }}}$, then

$\begin{array}{c}\hfill \begin{array}{cc}\hfill L{V}_{21}\le & -\frac{k}{{\mathcal{I}}^{\ast }}{(\mathcal{I}-{\mathcal{I}}^{\ast })}^2-\frac{k}{{\mathcal{I}}^{\ast }}(\mathcal{S}-{\mathcal{S}}^{\ast })(\mathcal{I}-{\mathcal{I}}^{\ast })-{\mu }_1{(\mathcal{S}-{\mathcal{S}}^{\ast })}^2\hfill \\ & -{\mu }_1(\mathcal{S}-{\mathcal{S}}^{\ast })(\mathcal{I}-{\mathcal{I}}^{\ast })+\frac{{\sigma }_1^2}{{\alpha }_1^2}+\frac{{\sigma }_2^2}{{\alpha }_2^2}\hfill \\ \hfill =& -\frac{k}{{\mathcal{I}}^{\ast }}{(\mathcal{I}-{\mathcal{I}}^{\ast })}^2-{\mu }_1{(\mathcal{S}-{\mathcal{S}}^{\ast })}^2-(\frac{k}{{\mathcal{I}}^{\ast }}+{\mu }_1)(\mathcal{S}-{\mathcal{S}}^{\ast })(\mathcal{I}-{\mathcal{I}}^{\ast })+\frac{{\sigma }_1^2}{{\alpha }_1^2}+\frac{{\sigma }_2^2}{{\alpha }_2^2},\hfill \end{array}\end{array}$

$\begin{array}{c}\hfill \begin{array}{cc}\hfill L{V}_{22}=& (\mathcal{B}-{\mathcal{B}}^{\ast })(\frac{{\mu }_2{\mathcal{B}}^{\ast }\mathcal{I}}{{\mathcal{I}}^{\ast }}-{\mu }_2\mathcal{B})\hfill \\ \hfill =& {\mu }_2{\mathcal{B}}^{\ast }(\mathcal{B}-{\mathcal{B}}^{\ast })(\mathcal{I}-{\mathcal{I}}^{\ast })(\frac{\mathcal{I}}{{\mathcal{I}}^{\ast }}-\frac{\mathcal{B}}{{\mathcal{B}}^{\ast }})\hfill \\ \hfill =& {\mu }_2{\mathcal{B}}^{\ast }(\mathcal{B}-{\mathcal{B}}^{\ast })(\mathcal{I}-{\mathcal{I}}^{\ast })(\frac{\mathcal{I}}{{\mathcal{I}}^{\ast }}-1+1-\frac{\mathcal{B}}{{\mathcal{B}}^{\ast }})\hfill \\ \hfill =& {\mu }_2{\mathcal{B}}^{\ast }\frac{\mathcal{I}-{\mathcal{I}}^{\ast }}{{\mathcal{I}}^{\ast }}+{\mu }_2{\mathcal{B}}^{\ast }\frac{{\mathcal{B}}^{\ast }-\mathcal{B}}{{\mathcal{B}}^{\ast }}\hfill \\ \hfill =& \frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}(\mathcal{B}-{\mathcal{B}}^{\ast })(\mathcal{I}-{\mathcal{I}}^{\ast })-{\mu }_2{(\mathcal{B}-{\mathcal{B}}^{\ast })}^2,\hfill \end{array}\end{array}$

$\begin{array}{c}\hfill \begin{array}{cc}\hfill L{V}_{23}& =(\mathcal{S}-{\mathcal{S}}^{\ast }+\mathcal{I}-{\mathcal{I}}^{\ast }+\mathcal{B}-{\mathcal{B}}^{\ast })(\frac{{\beta }_1{\mathcal{S}}^{\ast }{\mathcal{B}}^{\ast }}{1+{\alpha }_1{\mathcal{B}}^{\ast }}+\frac{{\beta }_2{\mathcal{S}}^{\ast }{\mathcal{I}}^{\ast }}{1+{\alpha }_2{\mathcal{I}}^{\ast }}+{\mu }_1{\mathcal{S}}^{\ast }-\frac{{\beta }_1\mathcal{S}\mathcal{B}}{1+{\alpha }_1\mathcal{B}}\hfill \\ & -\frac{{\beta }_2\mathcal{S}\mathcal{I}}{1+{\alpha }_2\mathcal{I}}-{\mu }_1\mathcal{S})+(\mathcal{S}-{\mathcal{S}}^{\ast }+\mathcal{I}-{\mathcal{I}}^{\ast }+\mathcal{B}-{\mathcal{B}}^{\ast })(\frac{{\beta }_1\mathcal{S}\mathcal{B}}{1+{\alpha }_1\mathcal{B}}+\frac{{\beta }_2\mathcal{S}\mathcal{I}}{1+{\alpha }_2\mathcal{I}}\hfill \\ & -\frac{{\beta }_1{\mathcal{S}}^{\ast }{\mathcal{B}}^{\ast }}{1+{\alpha }_1{\mathcal{B}}^{\ast }}\frac{\mathcal{I}}{{\mathcal{I}}^{\ast }}-\frac{{\beta }_2{\mathcal{S}}^{\ast }{\mathcal{I}}^{\ast }}{1+{\alpha }_2{\mathcal{I}}^{\ast }}\frac{\mathcal{I}}{{\mathcal{I}}^{\ast }})+(\mathcal{S}-{\mathcal{S}}^{\ast }+\mathcal{I}-{\mathcal{I}}^{\ast }+\mathcal{B}-{\mathcal{B}}^{\ast })(\frac{{\mu }_2{\mathcal{B}}^{\ast }\mathcal{I}}{{\mathcal{I}}^{\ast }}-{\mu }_2\mathcal{B})\hfill \end{array}\end{array}$

$\begin{array}{c}\hfill \begin{array}{cc}& +\frac{{\sigma }_1^2{\mathcal{S}}^2(t){\mathcal{B}}^2)}{{(1+{\alpha }_1\mathcal{B})}^2}+\frac{{\sigma }_2^2{\mathcal{S}}^2{\mathcal{I}}^2)}{{(1+{\alpha }_2\mathcal{I})}^2}\hfill \\ \hfill \le & (\frac{{\beta }_1{\mathcal{S}}^{\ast }{\mathcal{B}}^{\ast }}{1+{\alpha }_1{\mathcal{B}}^{\ast }}+\frac{{\beta }_2{\mathcal{S}}^{\ast }{\mathcal{I}}^{\ast }}{1+{\alpha }_2{\mathcal{I}}^{\ast }})(\mathcal{S}-{\mathcal{S}}^{\ast }+\mathcal{I}-{\mathcal{I}}^{\ast }+\mathcal{B}-{\mathcal{B}}^{\ast })(1-\frac{\mathcal{I}}{{\mathcal{I}}^{\ast }})\hfill \\ & -{\mu }_1(\mathcal{S}-{\mathcal{S}}^{\ast })(\mathcal{S}-{\mathcal{S}}^{\ast }+\mathcal{I}-{\mathcal{I}}^{\ast }+\mathcal{B}-{\mathcal{B}}^{\ast })+{\mu }_2{\mathcal{B}}^{\ast }(\mathcal{S}-{\mathcal{S}}^{\ast }+\mathcal{I}-{\mathcal{I}}^{\ast }\hfill \\ & +\mathcal{B}-{\mathcal{B}}^{\ast })(\frac{I}{{\mathcal{I}}^{\ast }}-1)+{\mu }_2{\mathcal{B}}^{\ast }(\mathcal{S}-{\mathcal{S}}^{\ast }+\mathcal{I}-{\mathcal{I}}^{\ast }+\mathcal{B}-{\mathcal{B}}^{\ast })(1-\frac{\mathcal{B}}{{\mathcal{B}}^{\ast }})+\frac{{\sigma }_1^2}{{\alpha }_1^2}+\frac{{\sigma }_2^2}{{\alpha }_2^2}\hfill \\ \hfill =& \frac{-k}{{\mathcal{I}}^{\ast }}{(\mathcal{I}-{\mathcal{I}}^{\ast })}^2-\frac{k}{{\mathcal{I}}^{\ast }}(\mathcal{I}-{\mathcal{I}}^{\ast })(\mathcal{S}-{\mathcal{S}}^{\ast })-\frac{k_1}{{\mathcal{I}}^{\ast }}(\mathcal{I}-{\mathcal{I}}^{\ast })(\mathcal{B}-{\mathcal{B}}^{\ast })-{\mu }_1{(\mathcal{S}-{\mathcal{S}}^{\ast })}^2\hfill \end{array}\end{array}$

$\begin{array}{c}\hfill \begin{array}{cc}& -{\mu }_1(\mathcal{S}-{\mathcal{S}}^{\ast })(\mathcal{I}-{\mathcal{I}}^{\ast })-{\mu }_1(\mathcal{S}-{\mathcal{S}}^{\ast })(\mathcal{B}-{\mathcal{B}}^{\ast })+\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}{(\mathcal{I}-{\mathcal{I}}^{\ast })}^2\hfill \\ & +\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}(\mathcal{I}-{\mathcal{I}}^{\ast })(\mathcal{S}-{\mathcal{S}}^{\ast })+\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}(\mathcal{I}-{\mathcal{I}}^{\ast })(\mathcal{B}-{\mathcal{B}}^{\ast })-{\mu }_2{(\mathcal{B}-{\mathcal{B}}^{\ast })}^2\hfill \\ & -{\mu }_2(\mathcal{B}-{\mathcal{B}}^{\ast })(\mathcal{B}-{\mathcal{B}}^{\ast })-{\mu }_2(\mathcal{I}-{\mathcal{I}}^{\ast })(\mathcal{B}-{\mathcal{B}}^{\ast })+\frac{{\sigma }_1^2}{{\alpha }_1^2}+\frac{{\sigma }_2^2}{{\alpha }_2^2}\hfill \\ \hfill =& (\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}-\frac{k}{{\mathcal{I}}^{\ast }}){(\mathcal{I}-{\mathcal{I}}^{\ast })}^2-{\mu }_1{(\mathcal{S}-{\mathcal{S}}^{\ast })}^2-{\mu }_2{(\mathcal{B}-{\mathcal{B}}^{\ast })}^2+(\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}\hfill \\ & -\frac{k}{{\mathcal{I}}^{\ast }}-{\mu }_1)(\mathcal{I}-{\mathcal{I}}^{\ast })(\mathcal{S}-{\mathcal{S}}^{\ast })+(\frac{{\mu }_2{B}^{\ast }}{{\mathcal{I}}^{\ast }}-\frac{k}{{\mathcal{I}}^{\ast }}-{\mu }_2)(\mathcal{B}-{\mathcal{B}}^{\ast })(\mathcal{I}-{\mathcal{I}}^{\ast })\hfill \\ & -({\mu }_1+{\mu }_2)(\mathcal{S}-{\mathcal{S}}^{\ast })(\mathcal{B}-{\mathcal{B}}^{\ast })+\frac{{\sigma }_1^2}{{\alpha }_1^2}+\frac{{\sigma }_2^2}{{\alpha }_2^2}\hfill \\ \hfill \le & (\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}-\frac{k}{{\mathcal{I}}^{\ast }}){(\mathcal{I}-{\mathcal{I}}^{\ast })}^2-{\mu }_1{(\mathcal{S}-{\mathcal{S}}^{\ast })}^2-{\mu }_2{(\mathcal{B}-{\mathcal{B}}^{\ast })}^2+(\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}\hfill \\ & -\frac{k}{{\mathcal{I}}^{\ast }}-{\mu }_1)(\mathcal{I}-{\mathcal{I}}^{\ast })(\mathcal{S}-{\mathcal{S}}^{\ast })+(\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}-\frac{k}{{\mathcal{I}}^{\ast }}-{\mu }_2)(\mathcal{B}-{\mathcal{B}}^{\ast })(\mathcal{I}-{\mathcal{I}}^{\ast })\hfill \\ & +\frac{\rho }{2}({\mu }_1+{\mu }_2){(\mathcal{S}-{\mathcal{S}}^{\ast })}^2+\frac{1}{2\rho }({\mu }_1+{\mu }_2){(\mathcal{I}-{\mathcal{I}}^{\ast })}^2+\frac{{\sigma }_1^2}{{\alpha }_1^2}+\frac{{\sigma }_2^2}{{\alpha }_2^2}.\hfill \end{array}\end{array}$

Hence,

$\begin{array}{c}\hfill \begin{array}{cc}\hfill L{V}_2\le & c_1(-\frac{k}{{\mathcal{I}}^{\ast }}{(\mathcal{I}-{\mathcal{I}}^{\ast })}^2-{\mu }_1{(\mathcal{S}-{\mathcal{S}}^{\ast })}^2-(\frac{k}{{\mathcal{I}}^{\ast }}+{\mu }_1)(\mathcal{S}-{\mathcal{S}}^{\ast })(\mathcal{I}-{\mathcal{I}}^{\ast })+\frac{{\sigma }_1^2}{{\alpha }_1^2}+\frac{{\sigma }_2^2}{{\alpha }_2^2})\hfill \\ & +c_2(\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}(\mathcal{B}-{\mathcal{B}}^{\ast })(\mathcal{I}-{\mathcal{I}}^{\ast })-{\mu }_2{(\mathcal{B}-{\mathcal{B}}^{\ast })}^2)\hfill \\ & +(\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}-\frac{k}{{\mathcal{I}}^{\ast }}){(\mathcal{I}-{\mathcal{I}}^{\ast })}^2-{\mu }_1{(\mathcal{S}-{\mathcal{S}}^{\ast })}^2-{\mu }_2{(\mathcal{B}-{\mathcal{B}}^{\ast })}^2\hfill \\ & +(\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}-\frac{k}{{\mathcal{I}}^{\ast }}-{\mu }_1)(\mathcal{I}-{\mathcal{I}}^{\ast })(\mathcal{S}-{\mathcal{S}}^{\ast })+(\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}-\frac{k}{{\mathcal{I}}^{\ast }}-{\mu }_2)(\mathcal{B}-{\mathcal{B}}^{\ast })(\mathcal{I}-{\mathcal{I}}^{\ast })\hfill \\ & +\frac{\rho }{2}({\mu }_1+{\mu }_2){(\mathcal{S}-{\mathcal{S}}^{\ast })}^2+\frac{1}{2\rho }({\mu }_1+{\mu }_2){(\mathcal{I}-{\mathcal{I}}^{\ast })}^2+\frac{{\sigma }_1^2}{{\alpha }_1^2}+\frac{{\sigma }_2^2}{{\alpha }_2^2}\hfill \\ & -[c_1{\mu }_1+{\mu }_2+\frac{\rho }{2}({\mu }_1+{\mu }_2)]{(\mathcal{S}-{\mathcal{S}}^{\ast })}^2-(c_1\frac{k}{{\mathcal{I}}^{\ast }}-\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}+\frac{k}{{\mathcal{I}}^{\ast }}){(\mathcal{I}-{\mathcal{I}}^{\ast })}^2\hfill \\ & -[c_2{\mu }_2+{\mu }_2+\frac{1}{2\rho }({\mu }_1+{\mu }_2)]{(\mathcal{B}-{\mathcal{B}}^{\ast })}^2+[-c_1(\frac{k}{{\mathcal{I}}^{\ast }}+{\mu }_1)+\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}-\frac{k}{{\mathcal{I}}^{\ast }}\hfill \\ & -{\mu }_1](\mathcal{I}-{\mathcal{I}}^{\ast })(\mathcal{S}-{\mathcal{S}}^{\ast })+(\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}-\frac{k}{{\mathcal{I}}^{\ast }}-{\mu }_2-c_2\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }})(\mathcal{B}-{\mathcal{B}}^{\ast })(\mathcal{I}-{\mathcal{I}}^{\ast })\hfill \\ & +\frac{{\sigma }_1^2}{{\alpha }_1^2}(1+c_1)+\frac{{\sigma }_2^2}{{\alpha }_2^2}(1+c_2)\hfill \end{array}\end{array}$

If $c_1$, $c_2$ satisfy the following equations

$-c_1(\frac{k}{{\mathcal{I}}^{\ast }}+{\mu }_1)+\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}-\frac{k}{{\mathcal{I}}^{\ast }}-{\mu }_1=0,$

$\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}-\frac{k}{{\mathcal{I}}^{\ast }}-{\mu }_2-c_2\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}=0,$

then

$\begin{array}{c}\hfill \begin{array}{cc}\hfill d{V}_2\le & -[c_1{\mu }_1+{\mu }_2+\frac{\rho }{2}({\mu }_1+{\mu }_2)]{(\mathcal{S}-{\mathcal{S}}^{\ast })}^2-(c_1\frac{k}{{\mathcal{I}}^{\ast }}-\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}+\frac{k}{{\mathcal{I}}^{\ast }}){(\mathcal{I}-{\mathcal{I}}^{\ast })}^2\hfill \\ & -[c_2{\mu }_2+{\mu }_2+\frac{1}{2\rho }({\mu }_1+{\mu }_2)]{(\mathcal{B}-{\mathcal{B}}^{\ast })}^2\hfill \\ & +\frac{{\sigma }_1^2}{{\alpha }_1^2}(1+c_1)+\frac{{\sigma }_2^2}{{\alpha }_2^2}(1+c_2).\hfill \end{array}\end{array}$

Integrating both sides of $d{V}_2$ from 0 to t, then taking the expectation, we have

$\begin{array}{c}\hfill \begin{array}{cc}\hfill 0& \le \mathbb{E}(V_2(t))-V(0)\hfill \\ & =\mathbb{E}{\int }_0^{t}L{V}_2(s)dt\hfill \\ & \le -\mathbb{E}{\int }_0^t[(c_1{\mu }_1+{\mu }_2+\frac{\rho }{2}({\mu }_1+{\mu }_2)){({\mathcal{S}}^{\ast }-\mathcal{S})}^2+(c_1\frac{k}{{\mathcal{I}}^{\ast }}-\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}+\frac{k}{{\mathcal{I}}^{\ast }}){({\mathcal{I}}^{\ast }-\mathcal{I})}^2\hfill \\ & +(c_2{\mu }_2+{\mu }_2+\frac{1}{2\rho }({\mu }_1+{\mu }_2)){(\mathcal{B}-{\mathcal{B}}^{\ast })}^2)]dr+[\frac{{\sigma }_1^2}{{\alpha }_1^2}(1+c_1)+\frac{{\sigma }_2^2}{{\alpha }_2^2}(1+c_2)]t.\hfill \end{array}\end{array}$

Hence,

$\begin{array}{c}\hfill \begin{array}{cc}\hfill \underset{t\to \infty }{lim}& \frac{1}{t}\mathbb{E}{\int }_0^t[(c_1{\mu }_1+{\mu }_2+\frac{\rho }{2}({\mu }_1+{\mu }_2)){({\mathcal{S}}^{\ast }-\mathcal{S})}^2+(c_1\frac{k}{{\mathcal{I}}^{\ast }}-\frac{{\mu }_2{\mathcal{B}}^{\ast }}{{\mathcal{I}}^{\ast }}+\frac{k}{{\mathcal{I}}^{\ast }}){({\mathcal{I}}^{\ast }-\mathcal{I})}^2\hfill \\ & +(c_2{\mu }_2+{\mu }_2+\frac{1}{2\rho }({\mu }_1+{\mu }_2)){(\mathcal{B}-{\mathcal{B}}^{\ast })}^2)]dr\hfill \\ & \le \frac{{\sigma }_1^2}{{\alpha }_1^2}(1+c_1)+\frac{{\sigma }_2^2}{{\alpha }_2^2}(1+c_2),\hfill \end{array}\end{array}$

which completes the proof.   □

Remark 3.

Theorem 3 shows that if the stochastic perturbation intensities ${\sigma }_1$ and ${\sigma }_2$ are small, the mean of the solution of system (4) converges to the endemic equilibrium of its corresponding deterministic system (1).

5. Simulations

Numerical results play an important role in studying the epidemic model. In this section, we shall present the numerical simulations to illustrate the extinction and persistence of cholera. Moreover, we will compare the solutions of the corresponding deterministic system (1) numerically.

As we know, the Euler–Maruyama Method (EM) is commonly used to solve the numerical solution of stochastic differential equations; EM has a strong order of convergence 1/2, whereas the underlying deterministic Euler method converges with classical order 1. It is possible to raise the strong order of EM to 1 by adding a correction to the stochastic increment, using Milstein’s method [19]. A so-called It$\widehat{o}$–Taylor expansion can be formed by applying It$\widehat{o}$’s result, which is a fundamental tool of stochastic calculus. Milstein’s method of stochastic differential Equation (5) can be obtained by truncating the It$\widehat{o}$–Taylor expansion at an appropriate point in the following form:

$\begin{array}{c}\hfill \begin{array}{c}{X}_j=X_{j-1}+\Delta tf(X_{j-1},t)+g(X_{j-1},t)(W({\tau }_j)-W({\tau }_{j-i}))+\frac{1}{2}g(X_{j-1},t)\hfill \\ g^{\prime }(X_{j-1},t){(W({\tau }_j)-W({\tau }_{j-1}))}^2-\Delta t),j=1,2,\cdots ,L,\hfill \end{array}\end{array}$

where $W(t)$ is a scalar standard Brownian motion on $[0,T]$, $\Delta t=T/L$ for some positive integer L and ${\tau }_j=j\Delta t.$

By using Milstein’s method, we have derived the following discretization equations from system (2):

$\begin{array}{cc}& {\mathcal{S}}_{k+1}={\mathcal{S}}_k+[{\mu }_1-\frac{{\beta }_1{\mathcal{S}}_k{\mathcal{B}}_k}{1+{\alpha }_1{\mathcal{B}}_k}-\frac{{\beta }_2{\mathcal{S}}_k{\mathcal{I}}_k}{1+{\alpha }_2{\mathcal{I}}_k}-{\mu }_1{\mathcal{S}}_k]\Delta t-{\sigma }_1\frac{{\mathcal{S}}_k{\mathcal{B}}_k}{1+{\alpha }_1{\mathcal{B}}_k}\sqrt{\Delta t}{\xi }_k-{\sigma }_2\frac{{\mathcal{S}}_k{\mathcal{I}}_k}{1+{\alpha }_2{\mathcal{I}}_k}\sqrt{\Delta t}{\xi }_k\hfill \\ & +\frac{{\sigma }_1^2}{2}\frac{{\mathcal{S}}_k{\mathcal{B}}_k}{1+{\alpha }_1{\mathcal{B}}_k}({\xi }_k^2-1)\Delta t+\frac{{\sigma }_2^2}{2}\frac{{\mathcal{S}}_k{\mathcal{I}}_k}{1+{\alpha }_2{\mathcal{I}}_k}({\xi }_k^2-1)\Delta t,\hfill \\ & {\mathcal{I}}_{k+1}={\mathcal{I}}_k+[\frac{{\beta }_1{\mathcal{S}}_k{\mathcal{B}}_k}{1+{\alpha }_1{\mathcal{B}}_k}+\frac{{\beta }_2{\mathcal{S}}_k{\mathcal{I}}_k}{1+{\alpha }_2{\mathcal{I}}_k}-(r+\alpha +{\mu }_1){\mathcal{I}}_k]\Delta t+{\sigma }_1\frac{{\mathcal{S}}_k{\mathcal{B}}_k}{1+{\alpha }_1{\mathcal{B}}_k}\sqrt{\Delta t}{\xi }_k\hfill \\ & +{\sigma }_2\frac{{\mathcal{S}}_k{\mathcal{I}}_k}{1+{\alpha }_2{\mathcal{I}}_k}\sqrt{\Delta t}{\xi }_k+\frac{{\sigma }_1^2}{2}\frac{{\mathcal{S}}_k{\mathcal{B}}_k}{1+{\alpha }_{1}B{\mathcal{B}}_k}({\xi }_k^2-1)\Delta t+\frac{{\sigma }_2^2}{2}\frac{{\mathcal{S}}_k{\mathcal{I}}_k}{1+{\alpha }_2{\mathcal{I}}_k}({\xi }_k^2-1)\Delta t,\hfill \\ & {\mathcal{R}}_{k+1}={\mathcal{R}}_k+[\alpha {\mathcal{I}}_k-{\mu }_1{\mathcal{R}}_k]\Delta t,\hfill \\ & {\mathcal{B}}_{k+1}={\mathcal{B}}_k+[\eta {\mathcal{I}}_k-{\mu }_2{\mathcal{B}}_k]\Delta t,\hfill \end{array}$

where $\Delta t$ is a positive time increment parameter and ${\xi }_k$ ($k=1,2,3,\cdots$) are independent Gaussian random variables $N(0,1)$. We scale up the total population from 1 to 10,000.

Case 1. We take the initial value $\mathcal{S}(0)=8500$, $\mathcal{I}(0)=500$, $\mathcal{R}(0)=1000$ and $\mathcal{B}(0)=30$, and the parameters’ values are listed in Case 1 of Table 1. We can obtain ${\mathfrak{R}}_0=0.7394<1$ and the conditions of Theorem 2 are satisfied. The difference between Figure 1 and Figure 2 is the intensity of environmental changes. We choose ${\sigma }_1={\sigma }_2=0.1$, ${\sigma }_1={\sigma }_2=0.05$ in Figure 1 and Figure 2, respectively. According to Theorem 2, the disease will proceed to extinction.

Case 2. We take the initial value $\mathcal{S}(0)=8500$, $\mathcal{I}(0)=500$, $\mathcal{R}(0)=1000$ and $\mathcal{B}(0)=30$, and the parameters’ values are listed in Case 2 of Table 1. We can obtain ${\mathfrak{R}}_0=30.2482>1$ and the conditions of Theorem 3 are satisfied. The difference between Figure 3 and Figure 4 is the intensity of environmental changes. We choose ${\sigma }_1={\sigma }_2=0.2$, ${\sigma }_1={\sigma }_2=0.1$ in Figure 3 and Figure 4, respectively. According to Theorem 3, the disease will persist. We can obtain that the solution curves of the stochastic system (2) always oscillate with respect to the endemic equilibrium of the deterministic system (1). From Figure 3 and Figure 4, we can also find that the amplitude of the solution curves of the stochastic model (2) increases with the increase in noise intensity.

6. Optimal Control Strategy for Cholera Disease

In this section, we shall focus on the optimal control problem of system (2). In order to do so, we introduce ${𝓊}_1(t)$, ${𝓊}_2(t)$ as two control measures that can affect cholera. For simplicity, the force of infection for both human–environment and human–human interactions is reduced by $1-{𝓊}_1(t)$, where precautionary measures are denoted by ${𝓊}_1(t)$. The second control variable ${𝓊}_2(t)$ represents disinfection to kill Vibrio cholerae and reduce the birth rate of Vibrio cholerae. System (4) is modified for optimal control as follows:

(13)$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle d\mathcal{S}(t)=[{\mu }_1-(1-{𝓊}_1(t))\frac{{\beta }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}-(1-{𝓊}_1(t))\frac{{\beta }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}-{\mu }_1\mathcal{S}(t)]dt}\hfill \\ \\ {\displaystyle -\frac{{\sigma }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}dw(t)-\frac{{\sigma }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}dw(t),}\hfill \\ {\displaystyle d\mathcal{I}(t)=[(1-{𝓊}_1(t))\frac{{\beta }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}+(1-{𝓊}_1(t))\frac{{\beta }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}-(r+\alpha +{\mu }_1)\mathcal{I}(t)]dt}\hfill \\ \\ {\displaystyle +\frac{{\sigma }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}dw(t)+\frac{{\sigma }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}dw(t),}\hfill \\ {\displaystyle d\mathcal{R}(t)=[\alpha \mathcal{I}(t)-{\mu }_1\mathcal{R}(t)]dt,}\hfill \\ {\displaystyle d\mathcal{B}(t)=[\eta \mathcal{I}(t)-{\mu }_2\mathcal{B}(t)-{𝓊}_2(t)\mathcal{B}(t)]dt,}\hfill \end{array}\right.\end{array}$

where

$\mathcal{S}(0)>0,\mathcal{I}(0)>0,\mathcal{R}(0)>0,\mathcal{B}(0)>0.$

Let us define the vectors for convenience:

(14)$\begin{array}{c}\hfill \begin{array}{c}x(t)={[x_1(t),x_2(t),x_3(t),x_4(t)]}^{\prime },𝓊(t)={[{𝓊}_1(t),{𝓊}_2(t)]}^{\prime }.\hfill \end{array}\end{array}$

and

(15)$\begin{array}{c}\hfill \begin{array}{c}dx(t)=f(x(t),𝓊(t))dt+g(x(t))dw(t).\hfill \end{array}\end{array}$

Moreover, the initial approximation is given by the following formula:

(16)$\begin{array}{c}\hfill \begin{array}{c}x(0)={[x_1(0),x_2(0),x_3(0),x_4(0)]}^{\prime }=x_0.\hfill \end{array}\end{array}$

Here, f and g are vectors and comprise the following components:

(17)$\begin{array}{c}\hfill \left\{\begin{array}{c}{\displaystyle f_1(x(t),𝓊(t))=[{\mu }_1-(1-{𝓊}_1(t))\frac{{\beta }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}-(1-{𝓊}_1(t))\frac{{\beta }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}-{\mu }_1\mathcal{S}(t)]dt}\hfill \\ \\ {\displaystyle -\frac{{\sigma }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}dw(t)-\frac{{\sigma }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}dw(t),}\hfill \\ {\displaystyle f_2(x(t),𝓊(t))=[(1-{𝓊}_1(t))\frac{{\beta }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}+(1-{𝓊}_1(t))\frac{{\beta }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}-(r+\alpha }\hfill \\ \\ {\displaystyle +{\mu }_1)\mathcal{I}(t)]dt+\frac{{\sigma }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}dw(t)+\frac{{\sigma }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}dw(t),}\hfill \\ {\displaystyle f_3(x(t),𝓊(t))=[\alpha \mathcal{I}(t)-{\mu }_1\mathcal{R}(t)]dt,}\hfill \\ {\displaystyle f_4(x(t),𝓊(t))=[\eta \mathcal{I}(t)-{\mu }_2\mathcal{B}(t)-{𝓊}_2(t)\mathcal{B}(t)]dt,}\hfill \end{array}\right.\end{array}$

$g_1=-\frac{{\sigma }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}-\frac{{\sigma }_2\mathcal{S}(t)I\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)},$

$g_2=-\frac{{\sigma }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}-\frac{{\sigma }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)},$

$g_3=g_4=0.$

Given the cost function, which is quadratic,

(18)$\begin{array}{c}\hfill \mathfrak{J}(𝓊)=\frac{1}{2}\mathbb{E}\{{\int }_0^{t^f}(A_1\mathcal{I}+A_2\mathcal{B}+\frac{B_1}{2}{𝓊}_1^2+\frac{B_2}{2}{𝓊}_2^2)dt\}.\end{array}$

Here, $A_1,A_2,B_1,B_2$ are four positive constants. Our goal is to seek out a control vector ${𝓊}^{\ast }(t)=({𝓊}_1^{\ast }(t),{𝓊}_2^{\ast }(t))$ such that

$\mathfrak{J}(u^{\ast })\le \mathfrak{J}(𝓊),\forall 𝓊\in \mathcal{U},$

where $\mathcal{U}$ is admissible control defined as

$\mathcal{U}=\{{𝓊}_i(t):{𝓊}_i(t)\in [0,{𝓊}_i^{max}],\forall {𝓊}_i\in L^2[0,t_f],t\in (0,t_f],i=1,2\}.$

Here, ${𝓊}_i^{max}$ ($i=1,2$) are positive real constants. In order to take advantage of the stochastic version of the maximum principle, firstly, we define the Hamiltonian function ${\mathcal{H}}_m(x,𝓊,p,q)$ in such a way that

$\mathcal{H}(x,𝓊,p,q)=⟨f(x,𝓊),p⟩+l(x,𝓊)+⟨g(x),q⟩.$

Here, the notation $⟨·,·⟩$ represents the Euclidean inner product, and $p={[p_1,p_2,p_3]}^{\prime }$ and $q={[q_1,q_2,q_3]}^{\prime }$ are the adjoint vectors. From the standard maximum principle, we obtain

(19)$d{x}^{\ast }(t)=\frac{\partial \mathcal{H}(x^{\ast },{𝓊}^{\ast },p,q)}{\partial p}dt+g(x^{\ast }(t))dw(t),$

(20)$d{p}^{\ast }(t)=-\frac{\partial \mathcal{H}(x^{\ast },{𝓊}^{\ast },p,q)}{\partial x}dt+q(x^{\ast }(t))dw(t),$

(21)${\mathcal{H}}_m(x^{\ast },{𝓊}^{\ast },p,q)=\underset{u\in \mathcal{U}}{max}{\mathcal{H}}_m(x^{\ast },{𝓊}^{\ast },p,q),$

where $x^{\ast }$ denotes an optimal trajectory of $x(t)$. The initial condition of Equation (19) is

(22)$x^{\ast }(0)=x_0.$

Moreover, the terminal condition of Equation (20) is

(23)$p(t_f)=-\frac{\partial h(x^{\ast }(t_f))}{\partial x}.$

Equation (21) implies that the optimal control ${𝓊}^{\ast }$ is a function of $p(t),q(t)$ and $x^{\ast }(t)$. Therefore, we have

${𝓊}^{\ast }(t)=\varphi (x^{\ast },p,q),$

where $\varphi$ is an unknown function and determined by Equation (23). Thus, Equations (19) and (20) can be expressed as

$d{x}^{\ast }(t)=\frac{\partial \mathcal{H}(x^{\ast },{𝓊}^{\ast },p,q)}{\partial p}dt+g(x^{\ast }(t))dw(t),$

$d{p}^{\ast }(t)=-\frac{\partial \mathcal{H}(x^{\ast },{𝓊}^{\ast },p,q)}{\partial p}dt+q(t)dw(t).$

Hence, the associated Hamiltonian is defined by

(24)$\begin{array}{cc}\mathcal{H}=\hfill & (A_1\mathcal{I}(t)+A_2\mathcal{B}(t)+\frac{B_1}{2}{𝓊}_1^2+\frac{B_2}{2}{𝓊}_2^2)\hfill \\ \hfill & +p_1({\mu }_1-(1-{𝓊}_1)\frac{{\beta }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}-(1-{𝓊}_1)\frac{{\beta }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}-{\mu }_1\mathcal{S}(t))\hfill \\ \hfill & +p_2((1-{𝓊}_1)\frac{{\beta }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}+(1-{𝓊}_1)\frac{{\beta }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}-(r+\alpha +{\mu }_1)\mathcal{I}(t))\hfill \\ \hfill & +p_3(\alpha \mathcal{I}(t)-{\mu }_1\mathcal{R}(t))+p_4(\eta \mathcal{I}(t)-{\mu }_2\mathcal{B}(t)-{𝓊}_2\mathcal{B}(t))\hfill \\ \hfill & +q_1(-\frac{{\sigma }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}-\frac{{\sigma }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)})+q_2(\frac{{\sigma }_1\mathcal{S}(t)\mathcal{B}(t)}{1+{\alpha }_1\mathcal{B}(t)}+\frac{{\sigma }_2\mathcal{S}(t)\mathcal{I}(t)}{1+{\alpha }_2\mathcal{I}(t)}).\hfill \end{array}$

From the stochastic maximum principle, we obtain

$dp(t)=-\frac{\partial \mathcal{H}(x^{\ast },{𝓊}^{\ast },p,q)}{\partial p}dt+q(t)dw(t).$

In other words,

(25)$\begin{array}{cc}\frac{d{p}_1}{dt}=\hfill & (p_1-p_2)(1-{𝓊}_1^{\ast }(t))(\frac{{\beta }_1{\mathcal{B}}^{\ast }(t)}{1+{\alpha }_1{\mathcal{B}}^{\ast }(t)}+\frac{{\beta }_2{\mathcal{I}}^{\ast }(t)}{1+{\alpha }_2{\mathcal{I}}^{\ast }(t)})+(q_1-q_2)\frac{{\sigma }_2{\mathcal{I}}^{\ast }(t)}{1+{\alpha }_2{\mathcal{I}}^{\ast }(t)}+{\mu }_1{p}_1,\hfill \\ \frac{d{p}_2}{dt}=\hfill & (p_1-p_2)(1-{𝓊}_1^{\ast }(t))\frac{{\beta }_2{\mathcal{S}}^{\ast }(t)}{{(1+{\alpha }_2{\mathcal{I}}^{\ast }(t))}^2}+p_2(r+\alpha +{\mu }_1)-p_3\alpha -p_4\eta \hfill \\ \hfill & +(q_1-q_2)\frac{{\sigma }_2{\mathcal{S}}^{\ast }(t)}{{(1+{\alpha }_2{\mathcal{I}}^{\ast }(t))}^2}-A_1,\hfill \\ \frac{d{p}_3}{dt}=\hfill & -p_3{\mu }_1,\hfill \\ \frac{d{p}_4}{dt}=\hfill & (p_1-p_2)(1-{𝓊}_1^{\ast }(t))\frac{{\beta }_1{\mathcal{S}}^{\ast }(t)}{{(1+{\alpha }_1\mathcal{B}{(t)}^{\ast })}^2}+p_4{𝓊}_2^{\ast }(t)+(q_1-q_2)\frac{{\sigma }_1{\mathcal{S}}^{\ast }(t)}{{(1+{\alpha }_1{\mathcal{B}}^{\ast }(t))}^2}-A_2,\hfill \end{array}$

along with the auxiliary initial and terminal conditions provided as

(26)${\mathcal{S}}^{\ast }(0)=\widehat{\mathcal{S}},{\mathcal{I}}^{\ast }(0)=\widehat{\mathcal{I}},{\mathcal{R}}^{\ast }(0)=\widehat{\mathcal{R}},{\mathcal{B}}^{\ast }(0)=\widehat{\mathcal{B}},p(t_f)=0.$

Now, by differentiating the Hamiltonian equations about $u_1$ and $u_2$, we receive optimal controls ${𝓊}_1^{\ast }$ and ${𝓊}_2^{\ast }$ as follows:

(27)${𝓊}_1^{\ast }=max\{min\{{𝓊}_1^{max},\frac{1}{B_1}(p_1-p_2)(\frac{{\beta }_1{\mathcal{B}}^{\ast }(t){\mathcal{S}}^{\ast }(t)}{1+{\alpha }_1{\mathcal{B}}^{\ast }(t)}+\frac{{\beta }_2{\mathcal{I}}^{\ast }(t){\mathcal{S}}^{\ast }(t)}{1+{\alpha }_2{\mathcal{I}}^{\ast }(t)})\},0\},$

(28)${𝓊}_2^{\ast }=max\{min\{{𝓊}_2^{max},\frac{p_4{\mathcal{B}}^{\ast }(t)}{B_2}\},0\}.$

Therefore, the optimal control problem of (13) is finished.

As part of the optimal control numerical simulation, we propose a new method, i.e., the stochastic Runge–Kutta technique method, which is formed by combining the fourth-order Runge–Kutta method and Milstein’s method [13]. The implementation steps of this method are as follows. The optimal strategy is obtained by solving the state system, adjoint system, and transverse conditions. Initially, we start by solving the state system (17) using the Milstein scheme and making initial guesses about the control of the simulation time. Then, the associated adjoint equation in system (25) is solved by a backward method with the transversality condition (26) using the existing iterative method of state equations. The control is then updated by using a convex combination of the previous control and the characterizations (27) and (28). This process is repeated and iterations stopped if the values of the unknowns at the previous iterations are very close to the ones at the present iterations.

By using the fourth-order Runge–Kutta and Milstein methods, the optimal control solutions of system (13) are solved, which are shown in Figure 5. Here, ${\mu }_1=0.0008$, ${\beta }_1=0.00214$, ${\beta }_2=0.2$, $r=0.013$, $\alpha =0.2$, ${\alpha }_1=0.4$, ${\alpha }_2=0.3$, $\eta =10$, ${\mu }_2=0.33$, ${\sigma }_1=0.5$, ${\sigma }_2=0.6.$

7. Conclusions

In this paper, the dynamics of a stochastic cholera epidemic model were investigated. We verified the existence and uniqueness of the positive global solution of the proposed model (4). We discussed the asymptotical stability of the disease-free equilibrium of model (4). The asymptotic stability of the equilibrium specific to the corresponding deterministic system (1) was given. The optimal control problem was considered and solved.

In this work, the stochastic system (2) was built based on the corresponding deterministic model. We obtained that the threshold value ${\mathfrak{R}}_0$ of extinction for the stochastic system (2) is equal to the deterministic model [6]. However, in fact, the conclusions of the stochastic system (2) are obtained in the case of small random disturbance intensity. If the random disturbance intensity becomes greater, we cannot obtain the same expression of ${\mathfrak{R}}_0$. Comparing Figure 1 and Figure 2, as well as Figure 3 and Figure 4, we find that even when the random disturbance is small, the disturbance intensity also reflects the intensity of the equilibrium shock of the corresponding deterministic model of the system unwinding. The larger the disturbance intensity, the larger the deviation between the solution of the stochastic system and the equilibrium of the deterministic model. It is also noted that the deterministic model is less pragmatic as compared to the stochastic model.

Cholera often occurs in economically underdeveloped countries; the cost of controlling the spread of the disease is a problem that the government needs to face, so the stochastic optimal control strategy is more practical. Compared with the analysis results of the general stochastic epidemic model [25,26,27], the stochastic optimal control strategy proposed by us is more advantageous for disease control. The general infectious disease model mostly obtains the threshold of disease extinction and disappearance—that is, the basic reproduction number. Most researchers give suggestions for disease control by controlling the parameters affecting the threshold of the basic reproduction number. However, the proposed method gives suggestions for disease control on the basis of the minimum cost. It also can be seen from numerical simulation that this control strategy can control the disease in a short time, as shown in Figure 5a.

In the numerical simulation of stochastic system solutions, we use the Milstein method, which further complements the principle of this method, and we point out that the Milstein method improves the order of strong convergence of the EM method. The stochastic Milstein technique meets the factual equilibria, and the deterministic outcome is its mean. The stochastic Runge–Kutta technique approaches the true equilibria of the model. However, they are both time-dependent, and the given technique may be divergent with an increase in time step size [28].

All of these theoretical findings are validated through the use of numerical simulations. This work shows that by using a stochastic epidemiological system, another option can be given to model epidemic dynamics.

Literature [28] compared the stochastic non-standard finite difference technique (SNSFD) [29,30] with the EM method, as well as the Milstein method, and the authors showed that SNSFD is the most appropriate technique to tackle all complex stochastic models, since it leads to realistic results, and this method has a number of favorable properties, such as the preservation of positivity and correct long-term behavior. The other, more complicated stochastic compartment models can be studied by applying the SNSFD technique in future work.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

I sincerely thank Xiangyun Shi and Yimeng Cao for their guidance and help in revising the paper. I thank the reviewers for their valuable comments.

Conflicts of Interest

The authors declare no conflict of interest.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Table

Figure 1. Paths of the deterministic model (1) and the corresponding stochastic model (2) with [Forumla omitted. See PDF.].

Figure 2. Paths of the deterministic model (1) and the corresponding stochastic model (2) with [Forumla omitted. See PDF.].

Figure 3. Paths of the deterministic model (1) and the corresponding stochastic model (2) with [Forumla omitted. See PDF.].

Figure 4. Paths of the deterministic model (1) and the corresponding stochastic model (2) with [Forumla omitted. See PDF.].

Figure 5. The plot represents [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.] with and without control. (a) [Forumla omitted. See PDF.]-with and without control in case of stochastic model. (b) [Forumla omitted. See PDF.]-with and without control in case of stochastic model. (c) [Forumla omitted. See PDF.]-with and without control in case of stochastic model. (d) [Forumla omitted. See PDF.]-with and without control in case of stochastic model.

Table 1

Parameter values.

Notation	Parameter Description	Case 1	Case 2
${\mu }_1$	Natural birth or death rate of humans	$5\times {10}^{-5}$/day [20]	$9.13\times {10}^{-5}$/day [21]
${\beta }_1$	Environment-to-human transmission rate	0.00011/day [Assumed]	0.213/day [21]
${\beta }_2$	Human-to-human transmission rate	0.0075/day [Assumed]	0.75/day [21]
${\alpha }_1$	The saturation constant	14.005 [Assumed]	0.3 [Assumed]
${\alpha }_2$	The saturation constant	16.05 [Assumed]	0.4 [Assumed]
$\alpha$	Recovery rate	0.2/day [21]	0.2/day [21]
${\mu }_2$	Net death rate of Vibrio cholerae	0.33/day [22]	0.33/day [22]
$\eta$	Rate of human contribution to Vibrio cholerae	10/day [22]	10/day [22]
r	Death rate due to cholera	0.015/day [23]	0.013/day [24]

References

1. Kaper, J.B.; Morris, J.G.; Levine, M.M. Cholera. Clin. Microbiol. Rev.; 1995; 8, pp. 48-86. [DOI: https://dx.doi.org/10.1128/CMR.8.1.48] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/7704895]

2. WHO. Cholera. 2021; Available online: https://www.who.int/news-room/fact-sheets/detail/cholera (accessed on 30 March 2022).

3. Zhou, X.Y.; Cui, J.A. Threshold dynamics for a cholera epidemic model with periodic transmission rate. Appl. Math. Model.; 2013; 37, pp. 3093-3101. [DOI: https://dx.doi.org/10.1016/j.apm.2012.07.044]

4. Yamazaki, K.; Wang, X. Global well-posedness and asymptotic behavior of solutions to a reaction-convection-diffusion cholera epidemic model. Discret. Contin. Dyn.-Syst.-Ser.; 2017; 21, pp. 1297-1316. [DOI: https://dx.doi.org/10.3934/dcdsb.2016.21.1297]

5. Zhang, X.H.; Peng, H. Stationary distribution of a stochastic cholera epidemic model with vaccination under regime switching. Appl. Math. Lett.; 2020; 102, 106095. [DOI: https://dx.doi.org/10.1016/j.aml.2019.106095]

6. Liao, S.; Wang, J. Global stability analysis of epidemiological models based on Volterra-Lyapunov stable matrices. Chaos Solitons Fractals; 2012; 45, pp. 966-977. [DOI: https://dx.doi.org/10.1016/j.chaos.2012.03.009]

7. Shi, X.Y.; Cao, Y.M.; Zhou, X.Y. Dynamics for a stochastic delayed SIRS epidemic model. Nonlinear Anal. Model. Control.; 2020; 25, pp. 705-725. [DOI: https://dx.doi.org/10.15388/namc.2020.25.17804]

8. Zhou, X.Y.; Shi, X.Y.; Wei, M. Dynamical behavior and optimal control of a stochastic mathematical model for cholera. Chaos Solitons Fractals; 2022; 156, 111854. [DOI: https://dx.doi.org/10.1016/j.chaos.2022.111854]

9. Zhou, X.Y.; Gao, X.W.; Shi, X.Y. Analysis of an SQEIAR stochastic epidemic model with media coverage and asymptomatic infection. Int. J. Biomath.; 2022; [DOI: https://dx.doi.org/10.1142/S1793524522500838]

10. Ikram, R.; Khan, A.; Zahri, M.; Saeed, A.; Yavuz, M.; Kumam, P. Extinction and stationary distribution of a stochastic COVID-19 epidemic model with time-delay. Comput. Biol. Med.; 2022; 141, 105115. [DOI: https://dx.doi.org/10.1016/j.compbiomed.2021.105115]

11. Zhang, X.B.; Shi, Q.H.; Ma, S.H.; Huo, H.F.; Li, D.G. Dynamic behavior of a stochastic SIQS epidemic model with Lévy jumps. Nonlinear. Dynam.; 2018; 93, pp. 1481-1493. [DOI: https://dx.doi.org/10.1007/s11071-018-4272-4]

12. Fatini, M.E.; Sekkak, I. Lévy noise impact on a stochastic delayed epidemic model with Crowly-Martin incidence and crowding effect. Phys. Stat. Mech. Its Appl.; 2020; 541, 123315. [DOI: https://dx.doi.org/10.1016/j.physa.2019.123315]

13. Garwal, R.P.; Badshah, Q.; ur Rahman, G.; Islam, S. Optimal control & dynamical aspects of a stochastic pine wilt disease model. J. Frankl. Inst.; 2019; 356, pp. 3991-4025.

14. Mao, X. Stochastic Differential Equations and Applications; Horwood Publishing: Chichester, UK, 1997.

15. Arnold, L. Stochastic Differential Equations; Wiley: New York, NY, USA, 1974.

16. Gronwall, T.H. Note on the derivatives with respect to a parameter of the solutions of a system of differential equations. Ann. Math.; 1919; 40, pp. 292-296. [DOI: https://dx.doi.org/10.2307/1967124]

17. Pang, Y.N.; Han, Y.C.; Li, W.J. The threshold of a stochastic SIQS epidemic model. Adv. Differ. Equ.; 2014; 2014, 320. [DOI: https://dx.doi.org/10.1186/1687-1847-2014-320]

18. Khasminskii, R. Ergodic properties of solutions of stochastic equations. Stoch. Model. Appl. Probab.; 2012; 4, pp. 99-144.

19. Higham, D.J. An algorithmic introduction to numerical simulation of stochastic differential equations. SIAM Rev.; 2001; 43, pp. 525-546. [DOI: https://dx.doi.org/10.1137/S0036144500378302]

20. Sanches, R.; Ferreira, C.; Kraenkel, R. The role of immunity and seasonality in cholera epidemics. Bull. Math. Biol.; 2011; 73, pp. 2916-2931. [DOI: https://dx.doi.org/10.1007/s11538-011-9652-6]

21. Hartley, D.M.; Morris, J.G., Jr.; Smith, D.L. Hyperinfectivity: A critical element in the ability of V. cholerae to cause epidemics. PLoS Med.; 2006; 3, pp. 63-69. [DOI: https://dx.doi.org/10.1371/journal.pmed.0030007]

22. Codeço, T.C. Endemic and epidemic dynamics of cholera: The role of the aquatic reservoir. BMC Infect. Dis.; 2001; 1, 1. [DOI: https://dx.doi.org/10.1186/1471-2334-1-1]

23. Mwasa, A.; Tchuenche, J. Mathematical analysis of a cholera model with public health interventions. Biosystems; 2011; 105, pp. 190-200. [DOI: https://dx.doi.org/10.1016/j.biosystems.2011.04.001]

24. Cui, J.A.; Wu, Z.M.; Zhou, X.Y. Mathematical analysis of a cholera model with vaccination. J. Appl. Math.; 2014; 2014, 324767. [DOI: https://dx.doi.org/10.1155/2014/324767]

25. Phan, T.A.; Tian, J.P.; Wang, B. Dynamics of cholera epidemic models in fluctuating environments. Stochastics Dyn.; 2021; 21, 2150011. [DOI: https://dx.doi.org/10.1142/S0219493721500118]

26. Hussain, G.; Khan, T.; Khan, A.; Inc, M.; Zaman, G.; Nisar, K.S.; Akgül, A. Modeling the dynamics of novel coronavirus (COVID-19) via stochastic epidemic model. Alex. Eng. J.; 2021; 60, pp. 4121-4130. [DOI: https://dx.doi.org/10.1016/j.aej.2021.02.036]

27. Din, A.; Li, Y. Stationary distribution extinction and optimal control for the stochastic hepatitis B epidemic model with partial immunity. Phys. Scr.; 2021; 96, 074005. [DOI: https://dx.doi.org/10.1088/1402-4896/abfacc]

28. Abodayeh, K.; Arif, M.S.; Raza, A.; Rafiq, M.; Bibi, M.; Nazeer, A. Numerical techniques for stochastic foot and mouth disease epidemic model with the impact of vaccination. Adv. Differ. Equ.; 2020; 2020, 34. [DOI: https://dx.doi.org/10.1186/s13662-020-2503-8]

29. Zafar, Z.U.; Ali, N.; Younas, S.; Abdelwahab, S.F.; Nisar, K.S. Numerical investigations of stochastic HIV/AIDS infection model. Alex. Eng. J.; 2021; 60, pp. 5341-5363.

30. Raza, A.; Arif, M.S.; Rafiq, M. A reliable numerical analysis for stochastic dengue epidemic model with incubation period of virus. Adv. Differ. Equ.; 2019; 2019, 32. [DOI: https://dx.doi.org/10.1186/s13662-019-1958-y]