Threshold dynamics and density function of a stochastic cholera transmission model

1.   Introduction

Cholera, a pervasive endemic disease, poses substantial risks to global public health, resulting in significant morbidity and mortality. The World Health Organization (WHO) reported that cholera incidence reached approximately 3.1 million cases and 95,000 fatalities in 2022, marking a $145\%$ increase relative to the average of the preceding five years. Transmission of cholera primarily occurs through the ingestion of the pathogen, characterizing it as both a waterborne and foodborne disease. The consumption of water contaminated with sewage along with the ingestion of victuals prepared in unsanitary conditions can significantly facilitate the transmission of the pathogen. Additionally, direct interpersonal contact with infected individuals can transmit the pathogen to susceptible hosts. High-risk individuals carrying the pathogen may spread cholera to their family members through close contact. Individuals who have recovered from cholera may develop a temporary acquired immunity to the pathogen. However, recent studies indicate that this acquired immunity may wane within a few months or weeks. These findings emphasize the critical necessity for ongoing research into the mechanisms of cholera transmission and the development of effective control measures. Various epidemiological frameworks and models have been extensively studied ^[1,2,3,4]. Gui ^[5] developed a four-dimensional epidemiological model for the dynamics of cholera transmission, as follows:

where $N (S + I + R = N)$ represents the total population of China. The population is separated into three groups: individuals who are susceptible $(S)$, infected $(I)$, and recovered $(R)$. In addition to these demographic groups, it is important to consider the potential implications of disinfection in controlling cholera, which is closely related to the concentration of vibrios in contaminated water (denoted by state B). In addition, Table 1 provides further information on parameters.

In system (1.1), the third equation is independent of the others, indicating that we only need to study the dynamics of the following subsystems

According to ^[5], the basic reproduction number is obtained by $R_0 = \beta_{h} \frac{\mu N}{(\mu+\nu)(\gamma+\mu)} +\bar{\beta}_e \frac{\mu N \xi }{(\mu+\nu)(\gamma+\mu)(\delta+c)\kappa}$. Moreover, the relevant threshold dynamics of system (1.2) is as follows:

$\bullet$ If $R_0 < 1,$ system (1.2) has a disease-free equilibrium $E_0 = (\frac{\mu N}{\mu + \nu}, 0, 0),$ which is globally asymptotically stable.

$\bullet$ If $R_0 > 1,$ there exists a unique endemic equilibrium $E^+ = (S^+, I^+, B^+) = (\frac{\mu N \xi -(\gamma+ \mu)(\delta+c)B^+}{(\mu+\nu) \xi}, \frac{ (\delta+ c)B^+}{\xi}, B^+)$ and it is globally asymptotically stable, where $B^+$ is the unique positive root of the following quadratic equation

and $\mu N \xi > (\gamma+ \mu)(\delta+c)B^+.$ Since it has been proven that the stochastic model can more accurately explain biological processes and infectious illnesses, there is increasing scholarly interest in examining the impact of environmental disturbances on epidemic models ^[6,7,8,9,10]. As a result, the development and research into stochastic models have intensified. For example, Jiang et al. ^[11] examined stationary distributions and extinction in non-autonomous logistic equations with random perturbations. In addition, several studies ^{[12,13,14,15]} have yielded notable conclusions in this field. One of the most important factors in the epidemic model (1.2) is $\beta_e$, which is always fluctuating around the average value $\bar{\beta}_e$ owing to the continuous spectrum of environmental noise. In this sense, $\beta_e$ should be considered a random variable. We assume that $\beta_e(t)$ is an Ornstein–Uhlenbeck process and satisfies the following form to imitate environmental noise's effect on transmission rates:

where $\alpha$ and $\sigma$ are positive constants indicating the speed of reversion and intensity of volatility, respectively, and $\bar{\beta}_e$ is a positive constant representing the long-term average transmission rate $\beta_e$. $W(t)$ is a standard Brownian motion defined on a complete probability space $(\varOmega, \mathscr {F}, \{\mathscr {F}_{t}\}_{_{t \geq 0}}, \mathbb{P })$ with a filtration $\{\mathscr {F}_{t}\}_{_{t \geq 0}}.$ According to Mao's monograph ^[16], we can obtain that

The expectation and variance of $\beta_e(t)$ are

As a result, the limit distribution of the Ornstein–Uhlenbeck process $\beta_e(t)$ is $\mathbb {N} (\bar{\beta}_e, \, \frac{\sigma^2}{2\alpha}).$ In other words, the probability density of the limit distribution is $\pi (x) = \frac{\sqrt{\alpha}}{\sqrt{\pi}\sigma}e^{-\frac{\alpha (x-\bar{\beta}_e)^2}{\sigma^2}}.$ Furthermore, it is straightforward to deduce that $\lim_{t\rightarrow0^+} \mathbb {E}[\beta_e(t)] = \beta_e (0)$ and $\lim_{t\rightarrow0^+} \mathrm{Var }[\beta_e(t)] = 0.$ This result suggests the Ornstein–Uhlenbeck process could be more suitable for describing random perturbations. Moreover, we let $\beta_e^+(t): = \max \{\beta(t), 0\},$ since the transmission rate coefficient should be non-negative. Thus, we obtain the following stochastic model:

In addition, we obtain

Considering the third equation of system (1.3), we have

Accordingly, region

is positively invariant set with respect to model (1.3). As a result, we take the supposition that initial values fulfill $S(0)+I(0) < N, \; \; B(0) < \frac{\xi N}{c}$, throughout the entire paper. We summarize our main contributions and innovations in comparison with the existing literature as follows: (ⅰ) To obtain the existence of a stationary distribution, which is a probability distribution with some invariant properties, we develop an innovative approach to build some stochastic Lyapunov functions. (ⅱ) During the following discussion, we explore the sufficient conditions for the existence of stationary distribution. The probability density function corresponding to the quasi-equilibrium is established, which is of great significance. (ⅲ) We also study the sufficient conditions for the disease to exterminate exponentially. In summary, the remainder of our paper is structured as follows. The existence and uniqueness of a global solution to system (1.3), as well as other relevant mathematical Lemma, are all presented in Section 2. The sufficient criteria for an ergodic stationary distribution of system (1.3) are derived in Section 3. To be complete, we build sufficient conditions for eliminating viruses and infected cells in Section 4. Lastly, numerical simulations are carried out to demonstrate our theoretical findings.

2.   Preliminaries

An important lemma is introduced in this section to obtain the accurate probability density expression.

Lemma 2.1. ^[17] For the algebraic equation $H_0^2+A_0\Sigma_0+\Sigma_0A_0^T = 0$, where $H_0 = diag(1, 0, 0, 0)$, $\Sigma_0$ is a real symmetric matrix and the standard matrix

If $\eta_1 > 0, \; \eta_3 > 0, \; \eta_4 > 0$ and $\eta_1\eta_2\eta_3-\eta_3^2-\eta_1^2\eta_4 > 0$, then $\Sigma_0$ is a positive definite matrix, where

Here, $A_0$ in this form is called the standard $R_1$ matrix.

Theorem 2.1. The system (1.3) has a unique global solution $(S(t), I(t), B(t), \beta_e(t))$ that will remain in $\Gamma$ with probability one $(a.s.)$ for any initial value of $(S(0), I(0), B(0), \beta_e(0))\in \Gamma$.

Proof. To prove the existence and uniqueness of positive solutions to stochastic models, we usually rely on the standard approach and the classical Khasminskii Lyapunov functional method, which is similar to ^[18]. Constructing the appropriate Lyapunov function is crucial. We construct a non-negativity $C^2$-function $U_0(S, I, B, \beta_e)$ as follows

Since $u-1-\ln u \geq 0$ for any $u > 0$, then the above function is non-negativity. Applying Itô's formula to $U_0$, we have

where $W_0$ is a positive constant, which is independent of $S, \; I, \; B$ and $\beta_e$.

3.   Stationary distribution

The aim of this section is to analyze whether the stochastic model has a stationary distribution that reflects infectious disease persistence. Define

where $\hat{\beta}_e = (\int_0^\infty x^{\frac{1}{4}}\pi (x) dx)^4, \; \; \pi (x) = \frac{\sqrt{\alpha}}{\sqrt{\pi}\sigma}e^{-\frac{\alpha (x-\bar{\beta}_e)^2}{\sigma^2}}$, $a_1 = \frac{\beta_h \mu N}{(\mu+ \nu)^2}, \; \; b_1 = \frac{\mu N \xi \hat{\beta}_e}{(\mu+\nu)^2(\delta+c) \kappa}.$

Theorem 3.1. Suppose that $R_0^s > 1$, then the stochastic system (1.3) admits at least one ergodic stationary distribution $\eta (\cdotp)$ on $\Gamma$.

Proof. Using Itô's formula to $-\ln S, \, -\ln I \, \text{and}\, -\ln B,$ respectively, we have

Then, define

where positive constants $a_1, b_1, b_2, b_3$ will be approved later. Using Ito's formula on $U_1$ and combining (3.1), we obtain

where

Note that $\beta_e^+ = \frac{\arrowvert \beta_e\arrowvert+\beta_e}{2},$ we have

Choose

Substituting (3.3) and (3.4) into (3.2), we have

where

Then, define

Making the use of Itô's formula to $U_2$, we have

where

Next, we define

Combining (3.1) and applying Itô's formula to $U_3,$ we have

where

Then, we define

where $M_0$ is a sufficiently large constant satisfying

Then, from (3.6)–(3.8), we have

where

Next, we construct a compact set $D_\varepsilon \subset \Gamma$ as follows

such that $g_3(S, I, B, \beta_e)\leq -1\, \text{for any} \, (S, I, B, \beta_e) \in \Gamma\setminus D_\varepsilon: = D^c_\varepsilon.$ Then, let $D_\varepsilon^c = \bigcup\limits_{i = 1}^6D_{\varepsilon, i}^c,$ where

$\varepsilon$ is a small enough constant satisfying the following inequalities

Case 1. If $(S, I, B, \beta_e) \in D_{\varepsilon, 1}^c,$ from (3.10) and (3.11), we have

Case 2. If $(S, I, B, \beta_e) \in D_{\varepsilon, 2}^c,$ from (3.10) and (3.12), we obtain

Case 3. If $(S, I, B, \beta_e) \in D_{\varepsilon, 3}^c,$ from (3.10) and (3.12), we derive

Case 4. If $(S, I, B, \beta_e) \in D_{\varepsilon, 4}^c,$ from (3.10) and (3.12), we get

Case 5. If $(S, I, B, \beta_e) \in D_{\varepsilon, 5}^c,$ from (3.10) and (3.12), we have

Case 6. If $(S, I, B, \beta_e) \in D_{\varepsilon, 6}^c,$ from (3.10) and (3.12), we obtain

In summary, we have

Let

Then, we have

The function $U_4$ has the minimum value $U_4(S^0, I^0, B^0, \beta_e^0),$ since it tends to $+\infty$ as $(S, I, B, \beta_e)$ approaches the boundary of $\varGamma.$ Thus, we obtain a non-negative function

Then, applying Itô's formula to $U$, we have

For any initial value $(S(0), I(0), B(0), \beta_e(0))\in \Gamma$ and a interval $[0, t],$ using the Itô's integral and then taking mathematical expectation to $U$, we get

Given the ergodicity of $\beta_e(t)$ and the strong law of large numbers, we have

Taking the inferior limit on both sides of (3.13) and combining with (3.12), (3.14) and (3.15), we get that

Therefore, we have

Making use of Fatou's lemma ^[19], we have

where $\mathbb{P }(t, (S, I, B, \beta_e, \mathbb{A})$ represents the transition probability of $(S(t), I(t), B(t), \beta_e(t))^T$ belonging to the set $\mathbb{A}$. According to the Lemma 2.1 in ^[20], system (1.3) has at least one stationary stochastic distribution when $R_0^s > 1.$

4.   Probability density function

The purpose of this section is to derive the explicit expression of probability density function around the quasi-endemic equilibrium based on the corresponding four-dimensional matrix equation.

The quasi-endemic equilibrium $P^* = (S^*, I^*, B^*, \beta_e^*)$ is the solution of the following equations.

By direct calculation, if $R_0 > 1,$ the solution of the above equation is unique, and it is

where $S^+, I^+\; \text{and }\; B^+$ are the same as Section 1.

Letting $(z_1, z_2, z_3, z_4)^T = (S-S^*, I-I^*, B-B^*, S-S^*, \bar{\beta}_e-\beta_e^*)^T,$ system (1.3) can be linearized as follows:

where

Letting

and

then, linearized Eq (4.2) can be expressed by matrix

Based on the continuous Markov processes in ^[21], system (1.3) has a unique probability density function $\Phi (z_1, z_2, z_3, z_4)$ surrounding the quasi-endemic equilibrium; according to Fokker–Plank equation, its form is as follows:

As a result, $\Phi (z_1, z_2, z_3, z_4)$ can be expressed by a quasi-Gaussian distribution since the diffusion matrix $G$ is constant. Therefore

where $q$ is a constant to ensure that the normalized condition is established $\int_{R^4}\Phi (z_1, z_2, z_3, z_4)\mathrm{d}z_1\mathrm{d}z_2\mathrm{d}z_3\mathrm{d}z_4 = 1.$ Then, $PG^2P+A^TP+PA = 0$ is satisfied by the real symmetric matrix $P.$ If the matrix $P$ is inverse, which denotes $P^{-1} = \Sigma,$ it can be equivalently changed into

To find the exact expression of the probability density function $\Phi (z_1, z_2, z_3, z_4),$ we will solve Eq (4.3).

Theorem 4.1. Assuming that if $R_0 > 1,$ then for any initial value $(S(0), I(0), B(0), \beta_e(0))\in \Gamma$, the stationary distribution of system (1.3) around $P^*$ has a normal density function as follows:

where $\Sigma$ is a positive definite matrix that takes the form

Here,

where $\varrho = -\sigma a_{14}m_1(a_{12}+a_{21}+a_{22}-a_{11})$, $m_1 = \frac{a_{32}(a_{21}+a_{33}-a_{11})}{a_{12}+a_{21}+a_{22}-a_{11}}$, $m_2 = m_1[a_{12}(a_{11}-a_{21}+a_{22}+a_{33})+a_{22}(a_{22}+a_{33})+a_{13}a_{32}+a_{33}^2]$, $m_3 = -a_{13}(a_{12}+a_{21}+a_{22}-a_{11})m_1-a_{33}^3$, $r_1 = -a_{14}m_1(a_{12}+a_{21}+a_{22}-a_{11})$, $r_2 = m_1(a_{12}+a_{22}+a_{21}-a_{11})(-a_{11}-a_{22}-a_{33})$, $\eta_1 = H_1+a_{44}, \, \eta_2 = H_2+H_1a_{44}, \eta_3 = H_3+H_2a_{44}$, $\eta_4 = H_3a_{44}, \, H_1 = a_{11}+a_{22}+a_{33}$, $H_2 = a_{11}(a_{22}+a_{33})+a_{12}a_{21}+(a_{22}a_{33}-a_{13}a_{32})$, $H_3 = a_{33}(a_{11}a_{22}+a_{12}a_{21})+a_{32}a_{13}(a_{21}-a_{11})$.

Proof. The first step we need is to confirm that $A$ is a Hurwitz matrix. Let $\varphi_A(\lambda)$ be the characteristic polynomial of matrix $A$, which can be written as follows:

where $H_1 = a_{11}+a_{22}+a_{33}, \, H_2 = a_{11}(a_{22}+a_{33})+a_{12}a_{21}+(a_{22}a_{33}-a_{13}a_{32})$, $H_3 = a_{33}(a_{11}a_{22}+a_{12}a_{21})+a_{32}a_{13}(a_{21}-a_{11}).$ The characteristic roots of $\varphi_A(\lambda)$ can be derived by: $\lambda_1 = -a_{44}$ and $\lambda^3+H_1\lambda^2+H_2\lambda+H_3 = 0.$ By calculations, we can obtain

which yields that $H_1 > 0, \, H_2 > 0, \, H_3 > 0,$ and $H_1H_2-H_3 > 0.$ This implies that all the roots of the characteristic equation (4.4) have negative real parts and the matrix $A$ is a Hurwitz matrix.

Next, by solving the equation $G^2+A\Sigma+\Sigma A^T = 0$, we shall find the form of $\Sigma$ and demonstrate that it is positive definite.

Let $A_1 = J_1AJ_1^{-1},$ where the ordering matrix $J_1$ is given by

Then,

Define $A_2 = J_2A_1J_2^{-1},$ where the elimination $J_2$ takes the form

then, we have

Next, let $A_3 = J_3A_2J_3^{-1},$ where

by simple calculation, then

where $m_1 = \frac{a_{32}(a_{21}+a_{33}-a_{11})}{a_{12}+a_{21}+a_{22}-a_{11}}.$

Let $A_4 = J_4A_3J_4^{-1},$ where the standardized transform matrix $J_4$ takes the form

where $r_1 = -a_{14}m_1(a_{12}+a_{21}+a_{22}-a_{11}), r_2 = m_1(a_{12}+a_{22}+a_{21}-a_{11})(-a_{11}-a_{22}-a_{33})$, $m_2 = m_1[a_{12}(a_{11}-a_{21}+a_{22}+a_{33})+a_{22}(a_{22}+a_{33})+a_{13}a_{32}+a_{33}^2]$, $m_3 = -a_{13}(a_{12}+a_{21}+a_{22}-a_{11})m_1-a_{33}^3.$

By direct calculation, we obtain

This characteristic polynomial is also invariant according to the invariant of matrix elementary transformations. Consequently, we can obtain

Additionally, Eq (4.3) is equivalently convertible into the following form

where $G_0 = diag(1, 0, 0, 0), \, \Sigma_0 = \varrho^{-2}(J_4J_3J_2J_1)\Sigma(J_4J_3J_2J_1)^T, \varrho = -\sigma a_{14}m_1(a_{12}+a_{21}+a_{22}-a_{11}).$

Using Lemma 2.1, the form of $\Sigma_0$ can be given as

More importantly, from $H_1 > 0, \, H_2 > 0, \, H_3 > 0, \, a_{44} > 0$ and $H_1H_2-H_3 > 0,$ we can deduce that $\eta_1 > 0, \, \eta_3 > 0, \, \eta_4 > 0$ and $\eta_1(\eta_2\eta_3-\eta_1\eta_4)-\eta_3^2 = (H_1+a_{44})(H_1H_2-H_3)a_{44}^2+H_3(H_1H_2-H_3)+H_2(H_1H_2-H_3)a_{44} > 0.$ So, the matrix $\Sigma_0$ is a positive definite matrix. Therefore, the matrix $\Sigma = \varrho^2(J_4J_3J_2J_1)^{-1}\Sigma_0[(J_4J_3J_2J_1)^{-1}]^T$ is also positive definite. As a result, the density function around quasi-endemic equilibrium $P^* = (S^*, I^*, B^*, \beta^*)$ is as follows

5.   Extinction

In this section, we explore the condition of disease extinction. Defining

Theorem 5.1. If $R_0^E < 1,$ then the disease of system (1.3) will exponentially die out a.s.

Proof. We define a $C^2$-function $F$ as follows

where $l_1 = (1+\frac{\nu}{\mu})R_0 \, \text{and}\, l_2 = \frac{\bar{\beta}_e N}{\kappa(\delta+c)}.$ By using the Itô's formula to $\ln F$, we obtain

By integrating (5.1) from $0$ to $t$ and dividing by $t$ on both sides, we get

Combining (3.14) and letting $t\rightarrow +\infty,$ we know that

Therefore, if $R_0^E < 1,$ then we have $\limsup _{t\rightarrow +\infty}\frac{\ln F(t)}{t}\leq \min\bigg\{\delta+c, \frac{\beta_hN\mu}{(\mu+\nu)R_0}\bigg\}(R_0^E-1) < 0,$ which indicates $\lim_{t\rightarrow +\infty}I(t) = 0$ and $\lim_{t\rightarrow +\infty}B(t) = 0.$ This suggests that system (1.3) sickness will disappear exponentially.

6.   Numerical simulation

Our analytical results are illustrated in this section through numerical simulations. Assuming $x(t) = \beta_e(t)-\bar{\beta}_e$, the discretization equation for system (1.3) is:

where the time increment $\triangle t > 0, \, \, \zeta_j$ are random variables and satisfy Gaussian distribution $N(0, 1)$ for $j = 1, 2\ldots, n.$

Example 6.1. Let $(S(0), I(0), B(0), x(0)) = (1.4008\times10^8, 3.944\times10^7, 1.088\times10^8, 1.36\times10^8),$ $\bar{\beta}_e = 5.35\times10^{-2}, \, \beta_h = 2.05\times 10^{-9}, \, \nu = 0.04/365$ and the other parameter values are presented in Table 1. Direct computation leads to that $(S^*, I^*, B^*, \beta_e^*) = (1.3704\times10^8, 3.8996\times10^7, 9.6683\times10^7, 5.35\times10^{-2})$ and

From Figure 1, it can be seen that system (1.3) has a stationary distribution.

Therefore, the solution $(S(t), I(t), B(t), \beta_e(t))$ of system (1.3) follows the normal density function $\Phi(S, I, B, \beta_e)\sim N\big((S^*, I^*, B^*, \beta_e^*)^T, \Sigma\big) = \mathbb{N }\big((1.3704\times10^8, 3.8996\times10^7, 9.6683\times10^7, 5.35\times10^{-2})^T, \Sigma\big).$

From Figure 2, we can find that the marginal density functions basically coincide with the corresponding fitting curves.

Example 6.2. Choose $\bar{\beta}_e = 3.1\times 10^{-4}, \beta_h = 1.32\times 10^{-5}, \nu = 0.31017/365, \sigma = 0.6, \alpha = 4,$ and the rest of the parameters unchanged, by direct calculation leads to $R_0^E = (1+\frac{\nu}{\mu})R_0+\frac{R_0(\mu+\nu)(\delta+c)\sigma}{\mu\bar{\beta}_e\sqrt{\pi\alpha}\min\bigg\{\delta+c, \frac{\beta_h\mu N}{(\mu+\nu)R_0}\bigg\}} = 0.8489 < 1.$ Thus, from Theorem 5.1, the disease of the system (1.3) will be extinct exponentially in a long time, which is supported by Figure 3.

7.   Conclusions

In this paper, on the basis of both biological significance and mathematically reasonable hypotheses, we illustrate that the Ornstein–Uhlenbeck process has a more stable variability than the linear perturbation and nonlinear perturbation. In this sense, we establish and analyze a stochastic cholera model with Ornstein–Uhlenbeck process to describe the propagation mechanism of cholera disease in the population. It is proved that the stochastic model has a unique global positive solution. We derive two critical conditions $R_0^s$ and $R_0^E;$ the system has at least one stationary distribution if $R_0^s > 1;$ the virus will be cleared out if $R_0^E < 1.$ The probability density function near the quasi-positive equilibrium is obtained by solving the corresponding Fokker–Planck equation. Mathematically, the existence of a stationary distribution implies the weak stability in the viewpoint of stochastic process. Biologically, the existence of a stationary distribution and probability density indicates the persistence of the disease. Finally, numerical simulations are used to verify our theoretical results.

Author contributions

Ying He: Conceptualization, Investigation, Formal analysis, Writing – review and editing. Bo Bi: Formal analysis, Writing – review and editing, Numerical simulation. All authors have read and approved the final version of the manuscript for publication.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work is supported by Hainan Provincial Natural Science Foundation of China (No. 122RC679,121RC554), Talent Program of Hainan Medical University (No. XRC2020030) and Northeast Petroleum University Special Research Team Project (No. 2022TSTD-05).

Conflict of interest

The authors declare there is no conflict of interest.