Analysis of a stochastic SIB cholera model with saturation recovery rate and Ornstein-Uhlenbeck process

1.   Introduction

Cholera is an acute intestinal infectious disease caused by Vibrio cholerae, which is mainly transmitted through unclean water and food (see ^[1]). According to news reports, more than 20 countries and regions (southwestern Cameroon, Afghanistan, the Republic of Mozambique, northern Iraq, Haiti, Malawi, etc.) have experienced or are experiencing cholera epidemics in the past two years. Therefore, many scholars have studied the cholera transmission model, such as differential dynamics system ^[2,3,4,5], age-structured transmission model ^[6], reaction-convection-diffusion equations ^[7], generalized fractional model ^[8].

In ^[9], the scholars investigated the following SIB cholera model:

where $S(t)$ and $I(t)$ denote the numbers of susceptible individuals and infected individuals at time t, respectively. $B(t)$ denotes the concentration of vibrios in contaminated water at time t. Besides, The parameter $A$ is the recruitment rate, and parameter $\mu_1$ is the natural human death rate. The parameters $\beta$ is the transmission coefficients of environment-to-human pathways. The parameter $K$ is the concentration of vibrios in contaminated water that yields $50\%$ chance of catching cholera. The parameter $\gamma$ is the recovery rate of infected individuals. $c$ is the maximum recovery per unit of time, and $b$ is the infected size at $50\%$ saturation. $\alpha$ is the disease-induced human death rate. The parameter $\eta$ is the contribution rate of each infected individual to the concentration of vibrios, and $\mu_2$ is the net death rate of vibrios. All parameters are usually also assumed to be nonnegative. For system (1.1), the basic reproduction number is defined by

where $S_0 = \frac{A}{\mu_1}$. In ^[9], the scholars obtained that if $R_0 < 1$, model has only a disease-free equilibrium $E_0 = (S_0, 0, 0)$ that is locally asymptotically stable; if $R_0 > 1$, system (1.1) has a unique endemic equilibrium $E^+ = (S^+, I^+, B^+)$ that is locally asymptotically stable.

In the real environment, the spread of cholera is disturbed by random factors, so many scholars have proposed a kind of stochastic differential equation cholera spread model with random perturbation, such as ^[10,11], stochastic cholera model between communities linked by migration ^[12], stochastic model with L$\acute{e}$vy process^[13], stochastic model under regime switching ^[14]. At present, the O-U process is popular in various random interference channels (see ^{[15,16,17,18]}). Therefore, in order to reveal the impact of environmental noise on the transmission rate, we assume that it is a random variable and satisfies the following form ^{[19,20,21,22]}:

where $\bar{\gamma}, \bar{c}$ are measure the long-time mean levels of the infection rates $\gamma, c$; $\lambda_i (i = 1, 2)$ are the speeds of reversion. $B_i(t)$ are independent standard Brownian motion parameters defined on a complete probability space $(\Omega, \mathcal{F}, P)$, and parameter $\sigma_i > 0\, (i = 1, 2)$ represents the intensity of $B_i(t)$. All parameters are usually also assumed to be nonnegative. In order to discuss the need for positivity of the stochastic model, variable $\max\{\gamma(t), 0\}, \max\{c(t), 0\}$ is used instead of variable $\gamma(t), c(t)$ in ^[19]. Then, we investigated the following stochastic SIB model with O-U process:

Throughout this paper, we define $\mathbb{R}_{+}^{3} = \{(x_1, x_2, x_3): x_i > 0, i = 1, 2, 3\}$. For an integrable function $f(t)$ defined on $[0, \infty)$, define $\langle f(t)\rangle = \frac{1}{t}\int_0^tf(s) \mathrm{d} s$. For the numbers $a$ and $b$, we define $a\vee b = \max\{a, b\}, a\wedge b = \min\{a, b\}$.

In the next section, we verify the existence and uniqueness of a global solution of the model (1.2) with any initial value. In Section 3, the criterion on the ergodicity and existence of unique stationary distribution for any solution of model (1.2) is stated and proved. In Section 4, We obtain the probability density function of model (1.2) around the positive equilibrium point $E^*$. In Section 5, the numerical examples are carried out to illustrate the main theoretical results.

2.   Existence and uniqueness of a global solution

In this section, we will discuss the existence and uniqueness of a global solution for model (1.2).

Theorem 2.1. For any initial value $(S(0), I(0), B(0), \gamma(0), c(0))\in\mathbb{R}_{+}^{3}\times\mathbb{R}^{2}$, model (1.2) has a unique global solution $(S(t), I(t), B(t), \gamma(t), c(t))$. That is, solution $(S(t), I(t), B(t), \gamma(t), c(t))$ is defined for all $t\geq 0$ and remains in $\mathbb{R}_{+}^{3}\times\mathbb{R}^{2}$with probability one.

Proof. Since the coefficients of model (1.2) satisfy the local Lipschitz conditions, for any initial value $(S(0), I(0), B(0), \gamma(0), c(0))\in\mathbb{R}_{+}^{3}\times\mathbb{R}^{2}$, there exists a unique local solution $(S(t), I(t), B(t), \gamma(t), c(t))$ on $t\in [0, \tau_e)$, where $\tau_e$ denotes the explosion time. To show this solution is global, we only need to prove that $\tau_e = \infty$ a.s.. To this end, let $k_0\geq 1$ be sufficiently large such that $S(0), I(0), B(0)$, $e^{\gamma(0)}$ and $e^{c(0)}$all lie within the interval $[\frac{1}{k_0}, k_0]$. For each integer $k\geq k_0$, define the stopping time

where throughout this paper, we set $\inf \emptyset = \infty$ (as usual $\emptyset$ represents the empty set). Clearly, $\tau_k$ is increasing as $k\rightarrow \infty$. Let $\tau_\infty = \lim_{k\rightarrow \infty}\tau_k$, whence $\tau_\infty \leq \tau_e$ a.s. If $\tau_\infty = \infty$ a.s. is not false, then $\tau_e = \infty$ a.s. and $(S(t), I(t), B(t), \gamma(t), c(t))\in \mathbb{R}_{+}^{3}\times\mathbb{R}^{2}$ a.s. for all $t > 0$. That is to say, if we want to finish the proof, we only need to show $\tau_\infty = \infty$ a.s. If this assertion is not true, then there is a pair of constants $T > 0$ and $\varepsilon \in (0, 1)$ such that

Consequently, there exists an integer $k_1\geq k_0$ such that

Define a Lyapunov function

By calculating, we have

By model (1.2), we have

This implies that

where $\tilde{N} = \max\{\frac{A}{\mu_1}, S(0)+I(0)\}$.

which implies that

where $\bar{N} = \max\{\frac{A\eta}{\mu_1\mu_2}, B(0)\}$. Then,

Therefore, integrating both sides of (2.4) from $0$ to $\tau_k\wedge T = \min\{\tau_k, T\}$ for any $k\geq k_1$ and then taking the expectations result in

Thus,

Set $\Omega_k = \{\tau_k\leq T\}$ for $k\geq k_1$, and in view of (2.1), we get $P(\Omega_k)\geq \varepsilon$. Notice that for every $\omega \in \Omega_k$, there exists $S(\tau_k, \omega), I(\tau_k, \omega), B(\tau_k, \omega), \gamma(\tau_k, \omega)$ or $c(\tau_k, \omega)$, which equals either $k$ or $\frac{1}{k}$. Therefore, $V(S(\tau_k, \omega), I(\tau_k, \omega), B(\tau_k, \omega), \gamma(\tau_k, \omega), c(\tau_k, \omega))$ is no less than either

Thereby, we can obtain

By (2.5), if follows that

where ${I}_{\Omega_k}$ denotes the indicator function of $\Omega_k$. Letting $k\rightarrow \infty$, then one can get that

which is a contradiction, and then we derive $\tau_\infty = \infty$. The proof is complete. Define the set $\Gamma$ as follows:

Corollary 2.2. For any initial value $x(0) = (S(0), I(0), B(0), \gamma(0), c(0))\in\mathbb{R}_{+}^{3}\times\mathbb{R}^{2}$ theglobal solution $x(t) = (S(t), I(t), B(t), \gamma(t), c(t))$ of model (1.2) ultimately enters into region $\Gamma$ with probability one as $t\to\infty$, and when $x(0)\in\Gamma$, then $x(t)\in\Gamma$ with probability one for all $t\geq 0$.

3.   Stationary distribution

In this section, we study the existence of the stationary distribution of model (1.2). Define

Theorem 3.1. Assume that $R_{0}^s > 1$, and then model (1.2) has at least one stationary distribution and ergodic property.

Proof. Define the Lyapunov function as follows:

where

By calculating, we have

where $y_1(t) = \gamma(t)-\bar{\gamma}, y_2(t) = c(t)-\bar{c}$ and

We have the following stochastic differential equation:

$y_i(t)$ has the ergodic property and density as follows:

and

Then,

where

Define the bounded closed set

where $\varepsilon$ are small enough positive constants, which will be determined later.

For convenience, we divide $\mathbb{R}_{+}^{3}\times\mathbb{R}^{2}\backslash U$ into six domains.

We will prove that $F(S, I, B, \gamma, c)\leq-1$ on $\mathbb{R}_{+}^{3}\times\mathbb{R}^{2}\backslash U$, which is equivalent to show it on the above seven domains.

Case 1. If $(S, I, B, \gamma, c)\in U_{1}$, we can obtain

where

We choose a constant $\varepsilon > 0$ small enough such that $-\frac{A}{\varepsilon}+G_{1}\leq-1$, and then it follows that

Case 2. If $(S, I, B, \gamma, c)\in U_{2}$, we can obtain

where

Choose constants $M > 0$ large enough and $\varepsilon > 0$ small enough such that $(c_1\frac{\beta}{K\mu_2}+c_3)M\eta \varepsilon+G_{2}\leq-1,$ and then it follows that

Case 3. If $(S, I, B, \gamma, c)\in U_{3}$, we can obtain

Choose a constant $\varepsilon > 0$ small enough such that $-\frac{\eta}{\varepsilon}+G_{1}\leq-1$, and then it follow that

Case 4. If $(S, I, B, \gamma, c)\in U_{4}$, we can obtain

Choose a constant $\varepsilon > 0$ small enough such that $-\frac{\alpha}{\varepsilon}+G_{1}\leq-1$, and then we have

Case 5. If $(S, I, B, \gamma, c)\in U_{5}$, we can obtain

Choose a constant $\varepsilon > 0$ small enough such that $-\frac{1}{2\varepsilon^2}+G_{1}\leq-1$, and then we get

Case 6. If $(S, I, B, \gamma, c)\in U_{6}$, we can obtain

Choose a constant $\varepsilon > 0$ small enough such that $-\frac{1}{2\varepsilon^2}+G_{1}\leq-1$, and then we get

Finally, from (3.3)-(3.7) we obtain

Define

where $(S^0, I^0, B^0, \gamma^0, c^0)$ is the point in the interior $\mathbb{R}_{+}^{3}\times\mathbb{R}^{2}$, such that $V(S, I, B, \gamma, c)$ will be minimized. By (2.2), we have

Taking the mathematical expectation and It$\hat{{\rm o}}$'s integral of (3.10), for any initial value $(S(0), I(0), B(0), \gamma(0), c(0))$, one has

By the ergodicity of $y_i (i = 1, 2)$ and the strong law of large numbers, one has

By (3.12), we have the inferior limit of (3.11),

From (3.9), we have

By (3.13) and (3.14), we have

In view of the definition of event probability and Fatou's lemma (see^[23]), the result of (3.15) is

where $P(t, (S, I, B, \gamma, c), U)$ is the transition probability of $(S(t), I(t), B(t), \gamma(t), c(t))$ belonging to set $U$. This shows that the solution $(S(t), I(t), B(t), \gamma(t), c(t))$ of model (1.2) has the Feller and ergodic property. This completes the proof.

4.   Density function

The positive equilibrium $P^* = (S^*, I^*, B^*, \gamma^*, c^*)$ of system (1.2) satisfied the following equations:

When $R_0 > 1$, we have

Let $(x_1, x_2, x_3, y_1, y_2) = (S-S^*, I-I^*, B-B^*, \gamma-\gamma^*, c-c^*)$. The linearization model of system (1.2) is:

where

Define

Let $X(t) = (x_1(t), x_2(t), x_3(t), y_1(t), y_2(t))^T, P(t) = (0, 0, 0, 0, B_1(t), B_2(t))^T$, and then the linear system (4.2) can be written as follows:

The characteristic equation of

where

and

Lemma 4.1. For the five-dimensional algebraic equation $P_0^2+D_0\Sigma_0+\Sigma_0 D_0^T = 0$, where $P_0 = diag(1, 0, 0, 0, 0)$, $\Sigma_0$ is a symmetrical matrix,

(1)

If $D_0$ is a Hurwitz matrix, that is, $a_i > 0 (i = 1, 2, 3), a_2a_1-a_3 > 0$, then $\Sigma_0$ is positive definite and has the following expression:

(2)

If $D_0$ is a Hurwitz matrix, that is, $a_i > 0 (i = 1, 2, 3, 4), a_2a_1-a_3 > 0$ and, $a_1 a_2 a_3-a_3^2-a_1^2 a_4 > 0$, then $\Sigma_0$ is positive semi-definite and has the following expression:

Theorem 4.2. If $R_0^s > 1$, the stationary solution $(S(t), I(t), B(t), \gamma(t), c(t))$ of model (1.2) around $(S^*, I^*, B^*, \gamma^*, c^*)^T$ has a unique normal density function, which takes the form

with the positive definite matrix $\Sigma = \Sigma_1+\Sigma_2.$

(1) $\mu_1\neq\mu_2$

(2) $\mu_1 = \mu_2$

$H_i(i = 1, 2, 3, 4, 5, 6, 7, 8, 9, 10)$, $\Sigma_{ij}(j, i = 1, 2)$ are given in the following proof.

Proof. Let the symmetric matrices $\Sigma = (r_{ij})_{5\times 5}$, where $r_{ji} = r_{ij}$. $\Sigma$ is determined by the following algebraic equation:

Furthermore, let $\Sigma = \Sigma_1+\Sigma_2$, $Q_1 = diag(0, 0, 0, 0, \sigma_1, 0), Q_2 = diag(0, 0, 0, 0, 0, \sigma_2)$, and then equation (4.5) can be decomposed into the following two equations: $Q_i^{2}+D\Sigma_i+\Sigma_iD^{T} = 0\; (i = 1, 2)$. We will calculate the matrices $\Sigma_i\; (i = 1, 2)$ in the following two steps.

Step (1) $Q_1^{2}+D\Sigma_1+\Sigma_1D^{T} = 0$. Let $D_1 = H_1DH_1^{-1}$, where

Let $D_2 = H_2D_1H_2^{-1}$, where

where $a_{12}-a_{22}-(a_{21}-a_{11}) = -\alpha$. Let $D_3 = H_3D_2H_3^{-1}$, where

Case 1. Let $a_{21}-a_{11}+\mu_2\neq0$, that is, $\mu_2\neq \mu_1$. Let $D_4 = H_4D_3H_4^{-1}$, where

Let $D_5 = H_5D_4H_5^{-1}$, where

where $a_{1}, a_{2}, a_{3}$ are given in (4.3).

Then, by Lemma 4.1(2), for the matrix equation

with $Q_0 = diag(1, 0, 0, 0, 0)$, we can obtain the positive definite matrix as follows:

where $H_{11} = (a_1a_2-a_3)(a_{3}+\lambda_1 a_{2}+\lambda_1^2a_1+\lambda_1^3)$. The element $a_{15}\eta(a_{21}-a_{11}+\mu_2)$ in the second row and fifth column of matrix $D_5$ does not affect the calculation results using the method of Lemma 4.1(2).

Then, we have that $Q_1^{2}+A\Sigma_1+\Sigma_1A^{T} = 0$ is equivalent to

which is further equivalent to

Therefore, we finally have that

Thus, $\Sigma_{1}$ is calculated and is positive semi-definite.

Case 2. $\mu_2 = \mu_1$. Let $D_6 = H_6D_{3}H_6^{-1}$, where

where $D^* = (a_{22}+a_{21})(a_{11}-a_{21})+a_{21}\alpha-a_{13}\eta$. Let $D_7 = H_{7}D_6H_{7}^{-1}$, where

where $D_{71} = -a_{13}\alpha(a_{11}+a_{22}+\mu_2), D_{72} = -a_{15}\alpha(a_{11}+a_{22}-\lambda_2)$. Then, by Lemma 4.1(1), for the matrix equation

we can obtain the positive define matrix as follows:

$H_{12} = (a_{11}+a_{22}+\lambda_1)a_{11}+a_{22}\lambda_1+a_{11}+a_{22}D^*$. The elements $-a_{13}\alpha, a_{15}\alpha$ in the second row, fourth column and fifth column of matrix $D_5$ do not affect the calculation results using the method of Lemma 4.1(1). Then, we have that $Q_1^{2}+D\Sigma_1+\Sigma_1D^{T} = 0$ is equivalent to

which is further equivalent to

Therefore, we finally have that

Thus, $\Sigma_{1}$ is calculated and is positive semi-definite.

Step (2) $Q_2^{2}+D\Sigma_2+\Sigma_2D^{T} = 0$. Let $D_8 = H_8DH_8^{-1}$, $D_9 = H_2D_8H_2^{-1}$, $D_{10} = H_3D_9H_3^{-1}$, where $H_2, H_3$ are given in step (1), and

Case 1. Let $\mu_2\neq\mu_1$. Let $D_{11} = H_4D_{10}H_4^{-1}$, where $H_4$ is given in step (1), and

where $a_{1}, a_{2}, a_{3}$ are given in (4.3). Let $D_{12} = H_9D_{11}H_9^{-1}$, where

Then, by Lemma 4.1(2), for the matrix equation:

we can obtain the positive define matrix as follows:

where $H_{21} = (a_1a_2-a_3)(a_{3}+\lambda_2 a_{2}+\lambda_2^2a_1+\lambda_2^3)$.

Then, we have that $Q_2^{2}+D\Sigma_2+\Sigma_2D^{T} = 0$ is equivalent to

which is further equivalent to

Therefore, we finally have that

Thus, $\Sigma_{2}$ is calculated and is positive semi-definite.

Case 2. $\mu_1 = \mu_2$. Let $D_{13} = H_6D_{10}H_6^{-1}$, where $H_6$ is given in step (1), and

Let $D_{14} = H_{10}D_{13}H_{10}^{-1}$, where

where $D_{73} = -a_{14}\alpha(a_{11}+a_{22}-\lambda_1)$. Then, by Lemma 4.1(1), for the matrix equation

we can obtain the positive define matrix as follows:

where $H_{22} = (a_{11}+a_{22}+\lambda_2)(a_{11}+a_{22})\lambda_2+(a_{11}+a_{22})D^*$. Then, we have that $Q_2^{2}+D\Sigma_2+\Sigma_2D^{T} = 0$ is equivalent to

which is further equivalent to

Therefore, we finally have that

Thus, $\Sigma_{2}$ is calculated and is positive semi-definite. Then, $\Sigma = \Sigma_1+\Sigma_2$ is positive define. Therefore, the expression of a normal density function around the quasi-endemic equilibrium of model (1.2) is obtained by

This completes the proof.

5.   Numerical examples

In model (1.2), we take the parameters $A = 2$, $\beta = 0.5$, $K = 0.4$, $\mu_1 = 0.2$, $\mu_2 = 0.5$, $\alpha = 0.2$, $\eta = 0.2$, $\theta_1 = \theta_2 = 0.2$, $\bar{c} = 0.06$, $\bar{\gamma} = 0.1$, $b = 1$ and $\sigma_1 = \sigma_2 = 0.02$. The numerical simulation of solution $(S(t), I(t), B(t))$ is done with initial value $(S(0), I(0), B(0)) = (3, 1.5, 1.5)$. We obtain $R_{0}^{S} = 8.5436 > 1$. The positive equilibrium of model (1.2) is $E^* = (S^*, I^*, B^*, \gamma^*, c^*)\approx(4.08, 2.96, 1.18, 0.1, 0.06)$. (See Fig. 1) The conditions in Theorem 3.1 and Theorem 4.2 are satisfied. Therefore, there exists the stationary distribution of model (1.2) (See Fig. 2), and there is a unique normal density function near equilibrium $E^*$. Then, the positive definite matrix $\Sigma$ is calculated as

and density function

Then, $S$, $I$ and $B$ have the following marginal probability densities:

The numerical simulation of the marginal probability densities for $S$, $I$ and $B$ is given in Figure 3.

6.   Conclusion

In this paper, we investigated a stochastic SIBS cholera model with saturation recovery rate and Ornstein-Uhlenbeck process. First, we proved that there is a unique global solution for any initial value of model. Secondly, the sufficient criterion of the stationary distribution of the model was obtained by constructing a suitable Lyapunov function, and the expression of the probability density function was calculated by the same condition. Finally, the correctness of the theoretical results is verified by numerical simulation, and the specific expression of the marginal probability density function is obtained.

This paper analyzes the interference of random disturbance in the process of recovering the infected person to the susceptible person during the transmission of cholera, but in fact, other parameters will also be more or less affected by random disturbance. In addition, people in many countries and regions have been vaccinated against cholera. In this process, cholera transmission is also affected by random factors. Therefore, there is still much work worthy of further study.

Acknowledgments

This research is supported by the National Natural Science Foundation of China (Grant No. 12201274, 12171004, 12001305), the Natural Science Foundation of Liaoning Province of China (Grant No. 2022-BS-287), the Basic Scientific Research Projects of Department of Education of Liaoning Province of China (Grant No. LJKQZ20222424) and the Scientific Research Foundation of Higher Education Institutions of Ningxia (Grant No. NGY2020005).

Conflict of interest

The authors declare there is no conflict of interest.