Threshold dynamics of stochastic cholera epidemic model with direct transmission

1.   Introduction

Infectious diseases significantly impact human health since they can have terrible effects. The dynamical behavior of infectious diseases must be understood in order to safeguard human health and control an infection ^[1]. Cholera, a diarrhea illness caused by Vibrio Cholera and characterized by diarrhea, is one of the most severe infectious diseases ^[2]. WHO has identified cholera as a public health issue ^[3]. Numerous theoretical and clinical research have been investigated in the literature. Cholera has been documented ^[4], but it still poses a serious hazard to public health in developing nations.

Mathematical modeling is a powerful tool for understanding how disease spreads ^[5]. Additionally, it provides future outbreak predictions and epidemic control strategies ^[6]. A disease's epidemic model may be deterministic or stochastic. An easy technique to analyze a system of ordinary differential equations is to use a deterministic model ^[7]. Due to their straightforward formulation, these models are commonly utilized. Deterministic models, however, have a number of limitations, including the fact that when the population is small enough, they are less informative and more challenging to estimate and analyze. Additionally, random effects disrupt deterministic models, causing trajectories to deviate from expected noise behavior ^[8]. Deterministic models have their strengths and can be useful for systems with well-defined rules and clear cause-and-effect relationships, but they often fall short in capturing the complexity, randomness, and uncertainty present in many real-world systems. Stochastic systems offer a more flexible and realistic way to model such complexity, making them a valuable tool in various scientific, engineering, and social applications^[9,10].

Various techniques can be used to study stochastic epidemic systems ^[11]. The time and state variables show that these processes are different from one another ^[12]. The literature has a variety of stochastic models that investigate how environmental noise affects the spread of infectious diseases. The stochastic models are efficient at calculating asymptotic expressions for the likelihood of an outbreak occurring ^[13]. One of the traditional models is a stochastic $SIS$ epidemiological model. Gray et al. ^[14] presents such a model whose basic characteristics with vaccination are studied in ^[15]. By employing the threshold quantity of the deterministic system, the authors of ^[16] investigated a stochastic $SIS$ epidemic model's dynamics and produced conditions for the noise of disease persistence and extinction. Some authors studied the long-term behavior of an $SIR$ stochastic epidemic model ^[17]. They established a threshold condition for the extinction and persistence of the model and supported their theoretical predictions through numerical simulations. Song et al. ^[18] used a stochastic SIRS model for studying the noise impact of an infectious disease with a saturated incidence rate and ^[19] reported a $SIQS$ model for investigating the threshold dynamics for white noise.

Numerous scholars reported different mathematical models to study the cholera epidemic ^[20,21,22]. Depending on the environmental noises, researchers studied a cholera stochastic system with vaccination and derived significant conditions regarding the basic threshold number ^[23]. The role of the aquatic reservoir on cholera disease was reported in ^[26], where they studied the endemic and epidemic dynamics of the infection. The impact of factors like contaminated water, temperature, rainfall, etc., on the cholera outbreak was investigated through a stochastic model system ^[27]. However, this work is solely concerned with theoretical considerations, and no analysis has been performed. More recently, authors reported a fundamental deterministic model that considers the direct contact transmission mechanism and the dynamics of the cholera sickness ^[28]. The model they used for the said investigation has the following structure:

A description of the parameters in system (1.1), for $t > 0$ is given below.

The term $S$ is the compartment of the susceptible individuals, $I$ is the compartment of the infected populace, $T$ denotes the class of treated people, and $R$ is the class of the recovered class.

The birth rate of the population is represented by $\Pi$ while the natural rate of mortality is symbolized by $\zeta$. $\delta$ denotes the contact rate between the susceptible class and the infected compartment; $\mu$ is the rate of treatment while death due to infection is denoted by $\alpha$ while deaths during therapy are described by $\gamma$ and $\eta$ is the rate of recovery, the loss rate of immunity is $\Lambda$.

Epidemic deterministic models assume even mixing in large populations, but this fails at outbreak start due to few infected people and stochastic transmission. Homogeneous mixing is not accurate for small initial infections, making deterministic models unsuitable. To address this, we introduce random parameter variation using a parametric perturbation approach, reflecting real-world heterogeneity and contact patterns. This enhances our understanding of disease dynamics during the early stages of epidemics. In this work, we extend the notion of the cholera deterministic model presented in ^[28] into a stochastic problem by taking direct human-to-human transmission into account. For this purpose, we apply the idea of parametric perturbation by choosing parameters from model (1.1) and transforming it to a random variable ^[9,18,19]. For more results, we refer to the following collection of articles in the book Disease Prevention and Health Promotion in Developing Countries ^[29].

We outline the environmental changes that affect the parameter $\delta$ such that $\delta \rightarrow \delta + \varepsilon d W(t)$. Here, the standard Brownian motion with intensity $\sigma^2 > 0$ is denoted by the symbol $W(t)$. Suppose $(\Omega, \mathcal{F}, \mathcal{P})$ represents an entire probability space with filtration so that $\{ \mathcal{F} \}_{t\geq0}$ meets the prerequisites (i.e., $\mathcal{F}_{0}$, including all null sets, is assumed to be continuous from the right and growing). Our stochastic model extends the deterministic one presented in ^[24] and up to the random excitations, it keeps similar extensions as in ^[23,30]. We study the following stochastic model with identic stochastic perturbation for the classes $S$ and $I$:

The manuscript is organized as follows: Section 2 of the manuscript is devoted to the basic preliminaries, including concepts and presentation of important formulas. In Section 3, the dynamic behavior of a positive solution is examined within a global context, utilize the Lyapunov analysis. Section 4 deals with exploration of the conditions under which the disease is guaranteed to become extinct with a probability of one. In Section 5, we establish the necessary conditions to demonstrate weak permanence and mean permanence with a probability one for the infection. Section 6 is dedicated to the execution of numerical simulations aimed at illustrating the core theoretical findings. The derivation of the basic threshold number is presented in Section 11. The calculated threshold number provides sufficient conditions for the infection to persist or eventually dies out. Finally we conclude our work in Section 7.

2.   Preliminaries

Introduce the notation $\mathcal{R}_{+}^{4} = \{(y_{1}, y_{2}, y_{3}, y_{4})| y_{i} > 0, i = 1, 2, 3, 4 \}$. Suppose that $(\mathcal{F}, \Omega, \mathcal{P})$ illustrates a probability measure space and the filtration $\{ \mathcal{F} \}_{t\geq0}$ observing the necessary relations (that is, all null sets include in $\mathcal{F}_{0}$, and is continuously increasing). Take into account the underlying $m$-dimensional stochastic system

such that $\mathfrak{f}(t, \chi)$ is $m$-dimensional vector function in $\mathcal{R}^{m}$ define in $[t_{0}, \infty]\times \mathcal{R}^{m}$, $\mathfrak{g}(t, \chi)$ is $m\times d$ matrix, $\mathfrak{f}$ and $\mathfrak{g}$ are locally Lipschitz functions in $\chi$. The differential operator $\mathcal{L}$ define by ^[12,31] associated with (2.1) as follow

Once the operator $\mathcal{L}$ is applied to the function $\mathscr{N} \in \mathcal{C}^{2, 1} ([t_{0}, \infty]\times \mathcal{R}^{m})$, one may consequently obtain

such that $\mathscr{N}_{t} = \frac{\partial \mathscr{N}}{\partial t}$, $\mathscr{N}_{y} = \sum_{i = 1}^{m}\frac{\partial \mathscr{N}}{\partial y_{i}}$, $V_{yy} = \sum_{i, j = 1}^{m}\frac{\partial^{2} \mathscr{N}}{\partial y_{i}\partial y_{j}}$. So we can define $It\hat{o}$ formula for $y(t) \in \mathcal{R}^{m}$ as

3.   Qualitative analysis of the system (1.2)

The first question to ask when examining dynamic behavior is whether the solution exists globally. Additionally, a populace dynamics system takes into account whether the solution is non-negative. So, in this part, we begin by demonstrating that the solution of model (1.2) is positive as well as global. Locally Lipschitz coefficients ensure a unique local solution for an ordinary differential equation on a given interval by guaranteeing a bounded rate of change for the solution within a neighborhood around each point. This controlled growth prevents solutions from diverging and overlapping, leading to a distinct trajectory for each initial condition, thereby ensuring uniqueness on the interval ^[25]. For this reason, the coefficients of the equation typically need to meet the local Lipschitz condition and the linear growth condition in order to produce a stochastic differential equation that has a distinct global solution for a given initial value ^[10]. Even if they are locally Lipschitz continuous, the coefficients of the model (1.2) do not meet the linear growth condition, hence the solution of the system (1.2) may blow up at a certain point in time. This part of the paper demonstrates that the solution of the system (1.2) is positive and global utilizing the Lyapunov analysis approach.

Theorem 1. The distinct solution of the stochastic model (1.2) on $t$ greater than or equal to zero is unique and for a starting point $(S_{0}, T_{0}, I_{0}, R_{0})\in \mathcal{R}^{4}$, with probability one will stay in $\mathcal{R}^{4}$, specifically $(S, T, I, R) \in \Omega$ for positive $t$ almost surely (a.s).

Proof. Since the parameters of the system (1.2) meet the local Lipschitz conditions, then there exists a distinct local solution $(I, T, S, R)$ for an assigned starting value given as $(S_{0}, T(_{0}, I_{0}, R_{0})\in R^{4}$on $[0, \tau_{e})$ a.s., with $\tau_{e}$ indicates the time of outbursts. Furthermore, to ascertain that solution is global, we require to verify $\tau_{e} = \infty$ as sure. Assume that $j_{0} > 0$ is sufficient large satisfies $(I_{0}, T_{0}, S_{0}, R_{0})$ remain in the interval $[\frac{1}{j_{0}}, j_{0}]$. For each integer $j\geq j_{0}$ we describe the time of stopping

In this paper, $\emptyset$ represents the empty set such that $\inf \emptyset = \infty$. By definition $\tau_{j}$ is rising as $j\rightarrow \infty$. Take $\tau_{\infty} = \lim_{j \rightarrow \infty} \tau_{j}$, from which $\tau_{e} \geq \tau_{\infty}$ a.s. Now we require to prove that $\tau_{\infty} = \infty$ a.s, therefore $\tau_{e} = \infty$ and $(I, T, S, R) \in R^{4}$ a.s. Furthermore, to accomplish this task, we must prove that $\tau_{\infty} = \infty$ a.s. Suppose that this assumption is false, then one may find some constants $B > 0$ and any $0 < \epsilon < 1$ for which $\mathcal{P} \{ \tau_{\infty} \leq B \} > \epsilon.$ As a significance of this, there exist an integer $j_{1}\geq j_{0}$ in such a manner

Clearly when $t \leq \tau_{j}$, we can write

which can be further written as

Now, let us define a Lyapunov function $V$ from $R^{4}$ to $\bar{R}$ as follows:

This function is obviously positive as $-\ln u-1+u \geq 0, \forall\; u > 0$. Suppose that $j \geq j_{0}$ and for arbitrary $B > 0$, applying the $It\hat{o}$ Integral formula, one may arrive at

According to the definition of the operator $\mathcal{L}$, we have

After further simplification, one may write

Inserting this into (3.1) leads to

Integrating (3.3) from 0 to $\tau_{j}\wedge B$ one may arrive at

where $\tau_{j}\wedge B = \min (\tau_{n}, t)$. After taking expectations on both sides, it becomes

Put $\Omega_{j} = \{\tau_{j}\leq B \}$ for $j \geq j_{1}$ and by (3.1), $P \{ \Omega_{k} \} \geq \epsilon$. It can be noticed that for each $\omega\in \Omega_{j}$, there corresponds at least one of $S(\tau_{j}, \omega), I(\tau_{j}, \omega), T(\tau_{j}, \omega), R(\tau_{j}, \omega)$ and that is equal to $j$ or $\frac{1}{j}$, and thus $V (S(\tau_{j}, \omega), I(\tau_{j}, \omega), T(\tau_{j}, \omega), R(\tau_{j}, \omega))$ which must not minimum than

Subsequently one may arrive at

Next, it is easy to write from Eqs (3.1) and (3.4) that

Note that $1_{\Omega_{k}}$ is the indicator function of $\Omega_{k}$. Suppose that $j \rightarrow \infty$ induces to the contradiction

Thus accordingly, it is necessary that

4.   Extinction of disease

For the dynamical behavior in epidemiology, our major concern is how to regulate the illness spread so that the infectious illness vanishes in a long term. For this purpose, the investigation of the sufficient conditions of disease extinction is performed in this section.

To establish sufficient conditions for the cholera disease extinction, we introduce the new notation in this paper: let $\langle x(t)\rangle = \frac{1}{t} \int x(r)dr$.

For our planned stochastic model (1.2), the basic reproduction value $\tilde{R}_{0}$ is described as

Lemma 1. Suppose that a continuous real-valued martingale denoted by $M = \{M_{t}\}_{t\geq0}$ is disappearing at positive $t$. So $\lim\limits_{t\rightarrow0} \langle M, M\rangle_{t} = \infty$ a.s., infer that $\lim\limits_{t\rightarrow0}\frac{M_{t}}{\langle M, M\rangle_{t}} = 0$ a.s. Moreover, $\lim\limits_{t\rightarrow0} \frac{\langle M, M\rangle_{t}}{t} < \infty$ almost surely, interprets that $\lim\limits_{t\rightarrow0} \frac{M_{t}}{t} = 0$ almost surely.

Theorem 2. Let $(I, S, T, R)$ be the solution of the planned system (1.2) with starting value $(I_{0}, T_{0}, S_{0}, R_{0})\in \eta^{*}$. If

a) $\varepsilon^{2} > \frac{\delta^{2}}{2(\zeta+\mu+\alpha)}$, therefore $\lim\limits_{t\rightarrow0} \sup \frac{\ln I}{t} \leq - (\zeta+\mu+\alpha)+\frac{\delta^{2}}{2\varepsilon^{2}} < 0$ very nearly surely;

b) $1 > \tilde{R}_{0}$ and $\varepsilon^{2}\leq \frac{\delta \zeta}{\Pi}$, therefore $\lim\limits_{t\rightarrow0} \sup \frac{\ln I}{t} \leq (\zeta+\mu+\alpha)(-1+\tilde{R}_{0}) < 0$ very nearly surely.

Proof. By using the $It\hat{o}$ integral formula to our proposed model (1.2), then

Integrating (4.2) from zero to $t$ and then take ratio of $t$, we obtain

where the non-discontinuous local martingale with property $M(0)$ is equal to zero, represented by $\; M(t) = \varepsilon \int_{0}^{t} S d W(t)$ ^[25]. Furthermore

Consider Lemma 1, which can be cast into

By applying the limit superior on both sides of (4.3), it follows that

If condition (a), namely,

holds, we obtain $\lim\limits_{t\rightarrow0} \sup\frac{\ln I}{t},$ {a.s.} Consequently, one may deduce that $\lim_{t\rightarrow0} I$ is equal to zero very nearly surely.

Integrating system (1.2) gives the underlying system

Further simplification leads to

After some manipulation, the last equation may written in the form

Consequently, we arrive at

Note that the value of $\varphi (t)$ is given by

It is apparent that as $t\rightarrow \infty$ so as $\varphi (t) \rightarrow 0$ a.s.,

Now we use the $It\hat{o}$ integral formula ^[25] to the model (1.2), then it attain the form

Integrating from zero to $t$, we reach the expression

Dividing both sides by $t$ to obtain

This can be further written as

where $M(t) = \varepsilon \int_{0}^{t}I(r) S(r) dW(r)$ is a continuous and local martingale which satisfies the condition $M(0) = 0$, see, e.g., ^[25]. By putting (4.5) in (4.8), one may obtain

Note that the term $\theta (t)$ in (4.9) is given by

Further

with the help of Lemma 1 can be written as

Next, by taking limit superior on both sides of (4.9) leads to

With the help of the above defined condition (b), inequality (4.11) gives

This gives $\lim\limits_{t\rightarrow0} \sup \frac{\ln I } {t} < 0$ ${\text{a.s}}.$, as a result of which we deduce that $\lim\limits_{t\rightarrow0} I = 0$ ${\text{ a.s}}.$

5.   Persistence of disease

This section of our manuscript is devoted for the construction of important relations for the disease dispersed in the populace. We follow the techniques discussed in ^[14,16,17]. The planned stochastic system (1.2) stated that persistence in means if

To proceed further, we consider it necessary to state the following lemmas:

Lemma 2. Assume that $f \in C [[0, \infty)\times(0, \infty)]$ and $F(t) \in C [[0, \infty)\times \Omega, R]$, in case one can find non-negative constants $\pi_{0}, \pi$ and $M$ satisfies

along with $\lim\limits_{t\rightarrow0} \frac{F(t)}{t} = 0$ almost surely. As a result $\frac{\pi}{\pi_{0}} \leq \lim\limits_{t\rightarrow0} \sup \frac{1}{t} \int_{0}^{t} f(s) ds$ almost surely.

Lemma 3. Assume that $f \in C [[0, \infty)\times(0, \infty)]$ and $F(t) \in C [[0, \infty)\times \Omega, R]$ if there exist non-negative constants $\pi_{0}, \pi$ and $M$ such that

along with $\lim\limits_{t\rightarrow0} \frac{F(t)}{t} = 0$ a.e., then $\frac{\pi}{\pi_{0}} \geq \lim\limits_{t\rightarrow0} \inf \frac{1}{t} \int_{0}^{t} f(s) ds$ a.e.

Theorem 3. If $\tilde{R_0} > 1$ and $\varepsilon^{2} < \frac{\delta \zeta}{\Pi}$, so for any starting value $(S_{0}, T_{0}, I_{0}, R_{0})\in \eta^{*}$, the solution $(S, T, I, R)$ of the planned cholera disease system (1.2) satisfies the property

very nearly surely, where

Proof. From the last inequality of (4.9), we may write

Using the value of $\tilde{ R}$ from (4.1), the last result takes the form:

Upon some algebraic manipulation, (5.1) yields:

Taking limit superior on both sides and using Lemma 2 together with (4.6), we arrive at:

On the other hand, from (4.5) and (4.7), one may write

Thus, we have

With the help of (4.1), the last in-equality gives

After some algebra, one may observe that (5.4) leads to the following equation:

Applying the inferior limit to each side and utilize Lemma 3 with (4.9) we can write

Thus from (4.5) and (4.9), we have

6.   Numerical simulation

In order to validate our theoretical results, we perform three numerical tests for the present $SITR$ stochastic system (1.2). This coupled stochastic model is derived from the classic $SITR$ model (1.1). Thus, we use a random excitation of the population classes $S$ and $I$. For computational simplicity, we use white noise as a random perturbation. These noises are independent. System (1.2) is discretized using the first-order stochastic Runge Kutta scheme. For such a stochastic framework, we refer to ^[32]. The derivation of this scheme is given as follows:

Where $\Delta t_n = t_{n+1}-t_n$ represents the time step and $\Delta W_{t_{n}} = W_{t_{n+1}}-W_{t_{n}}$ stands for the independent increments of the Gaussian Brownian motion generating the white noise. For simplicity, we implement a fixed time step $\Delta t_n = \Delta t$ for the evolution of all classes. We subdivide the temporal interval into 1, 000 equally spaced subintervals. Subject to various random initial conditions that satisfy our theoretical requirements, we numerically solve the $SITR$ system (1.2). We run our code in order to generate six simulations. The parameters of these six tests are given in Table 1. For tests 4–6, we use the same technique for the parameters $\sigma$ as implemented in ^[33,34]. Moreover, it should be stressed that the choice of our parameters proved in Theorems 2 and 3 is relevant for checking the extinctions and persistence scenarios. In tests 1–3, we simulate the extinction case and in tests 4 and 5 we simulate the persistence case. Obviously, all these parameters follow the conditions of our proven results. For every simulation, we show two samples of numerical realization and the associated numerical mean solution of a large number of solutions. Namely, we generate 1, 000 solutions for each test.

The initial values of all classes are represented by the values $S (0)$, $I (0)$, $T (0)$, and $R (0)$. The choice of these values is generated randomly in $[0, 1]$. This choice is updated by the following normalization at each time step:

In the following figures we provide two sample realizations of the stochastic numerical solution of the model under consideration. Test 1 based on simulations, we observe that all results obey the results of Theorems 1 and 2. For any on $t \geq 0$, the solution $(S, I, T, R)$ exits in $\mathbb{R}_+^4$. Also, all tests demonstrate a good accuracy and stability of the proposed model (1.2). Moreover, extinction and persistence scenarios are clearly demonstrated in the following simulations.

In Figures 1–3, we see the extinction of the class $I$ for the choice of $\widetilde{R}_0 < 1.$ The third column of these figures represents the asymptotical behavior of all population classes as a mean of 1000 random simulations using the same entries in the table above. In Figures 4–6, we show the persistence case, where the class of infected individuals $I$ remain persistent. This class will never vanish if $\widetilde{R}_0 > 1.$ Finally, it should be stressed that both out theoretical and computational results show similar conclusions as deduced in ^[23,24,30].

7.   Conclusions

Due to the nonregular and random happening, ecological phenomenon influences in the real world, such as absolute moisture, temperature, and rainfall, considerably affect the infection strength of diseases such as cholera. Thus, joining stochastic effects into the deterministic model provides a more realistic technique for modeling epidemic systems. We have evaluated a stochastic $SITR$ cholera system, which comprises variability in the direct transmission and functional our theoretical consequences to the dynamic of cholera based on realistic parameter values. To begin with, we extend the work ^[28] to the stochastic cholera model with random perturbations directly proportional to $T, I, S$, and $R$. Initially, we established the conditions for the persistence and extinction of the cholera infection. Furthermore, computer simulations confirm and indicate that white noise significantly affects the disease extinction and persistence of the cholera epidemic. Generally, we got similar results as in ^[23,24,30]. Similarly, the persistence of cholera infection is minimized with the noise strength increasing. We brought the numerical conclusions using the Runge-Kutta stochastic scheme, which supports our analytical results.

Use of AI tools declaration

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

Acknowledgments

The authors are thankful to Princess Nourah bint Abdulrahman University Researchers Supporting Project number (PNURSP2023R8), Princess Nourah bint Abdulrahman University, Riyadh, Saudi Arabia for financial support.

Conflict of interest

The authors certify that there are no reported conflicts of interest that may have appeared to have influenced the research presented in this study.

Appendix

Consider the infected compartment of stochastic system (1.2),

Let $f(t, I) = \log I.$ Using $It\hat{o}$ formula, one can write

By chain rule, we have

Now for next generation matrix, assume that

Again suppose that the Jacobian matrix of $f$ and $v$ w.r.t., I are $F$ and $V$ respectively, which can be write at the no-infection equilibrium point $E^{0} = (\frac{\Pi}{\zeta}, 0, 0, 0)$ as,

From the last two expressions, it is easy to write

The eigenvalue can be found by the Characteristics equation

After solving, we get the reproduction number of stochastic model as,