
==== Front
Heliyon
Heliyon
Heliyon
2405-8440
Elsevier

S2405-8440(24)11780-X
10.1016/j.heliyon.2024.e35749
e35749
Research Article
Analysis of a stochastic SEIuIrR epidemic model incorporating the Ornstein-Uhlenbeck process
Mediani Mhammed medi.mohammed@univ-adrar.edu.dz
a
Slama Abdeldjalil aslama@univ-adrar.edu.dz
a
Boudaoui Ahmed ahmedboudaoui@univ-adrar.edu.dz
a
Abdeljawad Thabet tabdeljawad@psu.edu.sa
thabetabdeljawad@gmail.com
bcde⁎
a Laboratory of Mathematics, Modeling and Applications (LaMMA), University of Adrar, Adrar, Algeria
b Department of Mathematics and Sciences, Prince Sultan University Riyadh, Saudi Arabia
c Department of Medical Research, China Medical University, Taichung 40402, Taiwan
d Center for Applied Mathematics and Bioinformatics (CAMB), Gulf University for Science and Technology, Hawally, 32093, Kuwait
e Department of Mathematics and Applied Mathematics, Sefako Makgatho Health Sciences University, Garankuwa, Medusa 0204, South Africa
⁎ Corresponding author. tabdeljawad@psu.edu.sathabetabdeljawad@gmail.com
06 8 2024
30 8 2024
06 8 2024
10 16 e357493 3 2024
30 7 2024
2 8 2024
© 2024 The Author(s)
2024
https://creativecommons.org/licenses/by-nc-nd/4.0/ This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).
This article aims to analyze a stochastic epidemic model SEIuIrR (Susceptible-exposed-undetected infected-detected infected (reported -recovered) assuming that the transmission rate at which people undetected become detected is perturbed by the Ornstein–Uhlenbeck process. Our first objective is to prove that the stochastic model has a unique positive global solution by constructing a nonnegative Lyapunov function. Afterward, we provide a sufficient criterion to prove the existence of an ergodic stationary distribution of the mode by constructing a suitable series of Lyapunov functions. Subsequently, we establish sufficient conditions for the extinction of the disease. Finally, a series of numerical simulations are carried out to illustrate the theoretical results.

Keywords

Stochastic epidemic model
Ornstein–Uhlenbeck process
Stationary distributions
Disease extinction
==== Body
pmc1 Introduction

Mathematical modeling has the potential to play a significant role in solving the problem of the spread of epidemics. Vaccination programs, physical separation and disease eradication efforts could all benefit from mathematical analysis of epidemic models. In recent years, mathematical models have been developed to analyze infectious diseases such as Covid 19, HIV/AIDS and influenza [17], [24], [26], [38], [43], [45], [47]. In the context of epidemic modeling, the widely used deterministic models are SIR (Susceptible-Infectious-Recovered) and SEIR (Susceptible-Exposed-Infectious-Recovered) models [7], [9], [20], [27], [28], [29], [34], [36], [42], among others. These deterministic models have provided invaluable insights into infection rates, healthcare system demands, and potential intervention effectiveness. However, to capture the memory and hereditary properties of biological systems more accurately, Proportional Caputo Fractional Derivative models are increasingly recognized as necessary [1], [2], [3], [4], [5], [6], [19], [37].

However, the intricacies of real-world dynamics and the inherent variability in human behavior introduce an element of randomness that cannot be ignored. From variations in individual susceptibility and contact patterns to the uncertainty surrounding the emergence of new viral strains, randomness permeates every aspect of the pandemic's progression. This realization has led to the evolution of epidemic modeling beyond deterministic analyses, prompting researchers to embrace stochastic processes and probabilistic methods to capture the unpredictable nature of the virus's spread [10], [13], [16], [39], [50]. Zhou et al. [49] formulated a stochastic SIR epidemic model with nonlinear incidence rate and general stochastic noises. Su and Zhang [41] proposed a stochastic SEI epidemic model in which the transmission rates are general functions and satisfy the log-normal Ornstein-Uhlenbeck process. Song and Zhang [40] studied a stationary distribution and the exponential extinction of a stochastic SVEIS epidemic model. Gatyeni et al. [18] suggested and analyzed a model incorporates the vital dynamics to capture the dynamics of COVID-19 infection using the South African setting in addition to optimizing the control strategies. Based on the work done by Gatyeni et al. [18] and El hadj Moussa et al. [15], we consider an SEIuIrR epidemic model as follows:(1) {dS(t)dt=Δ−β(v1Iu+v2Ir)S−μS,dE(t)dt=β(v1Iu+v2Ir)S−(σ+μ)E,dIu(t)dt=σ(1−ρ)E−(μ+d1+δ+γIu)Iu,dIr(t)dt=σρE+δIu−(μ+d2+γIr)Ir,dR(t)dt=γIuIu+γIrIr−μR,

where S(t) denotes susceptible population, E(t) represents exposed people, Iu(t) represents undetected infected people, Ir(t) represents the detected infected people (or reported), and R(t) denotes the recovered people. Assuming all parameters to be constant and positive, their respective descriptions are provided in Table 1. It is feasible to achieve that system (1) possesses a positive invariant region [15], [18]:Γ={(S,E,Iu,Ir,R)∈R+5,0≤S+E+Iu+Ir+R≤Δμ}.

The expression for the basic reproduction number of system (1) is provided as follows:R0=R1+R2+R3,

whereR1=σβS0v1(1−ρ)(σ+μ)(δ+γIu+μ+d1),R2=σβS0ν2ρ(σ+μ)(γIr+μ+d2),R3=σβS0v2(1−ρ)δ(σ+μ)(δ+γIu+μ+d1)(γIr+μ+d2),

which plays a crucial role in determining the occurrence of the disease, in a situation where S0=Δμ.Table 1 Description of state parameters of the model (1).

Table 1Parameter	Description	
Δ	Birth rate	
μ	Natural death rate	
β	Disease transmission rate	
ρ	Proportion of E which becomes undiagnosed	
v1	Transmissibility relative to undetected people	
v2	Transmissibility relative to detected people	
δ	The transmission rate at which undetected people become detected	
d1	Disease related death rate in Iu compartment	
d2	Disease related s death rate in Ir compartment	
σ	Incubation rate Fitted	
γIu	Recovery rate of people in Iu compartment	
γIr	Recovery rate of people in Ir compartment	

Furthermore, the system's (1) relevant threshold dynamics can be described in the following manner:• If R0<1. The system (1) exhibits a disease-free equilibrium E(S0,0,0,0,0)=(Δμ,0,0,0,0), which is locally asymptotically stable within the region Γ.

• If R0≥1. The system (1) possesses a globally asymptotically stable endemic equilibrium X⁎=(S⁎,E⁎,Iu⁎,Ir⁎,R⁎) within the region Γ.

Lately, an increasing number of researchers have shown a significant interest in investigating the dynamic behavior of epidemiological models incorporating stochastic perturbations, as it has become evident that the utilization of stochastic modeling techniques provides a more accurate representation of infectious diseases [21], [23], [31], [32].

Until now, multiple methods exist for incorporating stochastic perturbations into deterministic models. Among these approaches, one of the widely adopted strategies assumes that the system's parameters follow an Itô process known as the Ornstein-Uhlenbeck process, [8], [33], [46], [51], [52].

Motivated by the aforementioned discussions, this study considers the analysis of a stochastic SEIuIrR (Susceptible-exposed-undetected infected-detected infected (reported)-recovered) epidemic model by incorporating the Ornstein-Uhlenbeck process, as outlined in the following manner:dδ(t)=ρ1[δ¯−δ(t)]dt+σ1dB(t),

where δ¯ is measure the long-time mean level of the infection rate δ, ρ1 is the speeds of reversion, B(t) is independent standard Brownian motion parameter defined on a complete probability space (Ω,F,{Ft}t≥0,P) with a filtration {Ft}t≥0 satisfying the usual conditions, and σ1 represents the intensity of B(t). Typically, it is assumed that all parameters in the stochastic model are nonnegative to examine the necessity of incorporating positivity into the discussion, variable max⁡{δ(t),0} is used instead of variable δ(t) in [31]. Consequently, the ensuing stochastic model is derived as follows:(2) {dS(t)=[Δ−β(v1Iu+v2Ir)S−μS]dt,dE(t)=[β(v1Iu+v2Ir)S−(σ+μ)E]dt,dIu(t)=[σ(1−ρ)E−(μ+d1+max⁡{δ(t),0}+γIu)Iu]dt,dIr(t)=[σρE+max⁡{δ(t),0}Iu−(μ+d2+γIr)Ir]dt,dR(t)=[γIuIu+γIrIr−μR]dt,dδ(t)=ρ1[δ¯−δ(t)]dt+σ1dB(t).

For convenience and clarity, we introduce the following notation conventions:

R+m={(x1,...,xm)|xb≥0,1≤b≤m},p1∨p2=max⁡{p1,p2} for any p1,p2∈R. Consider 1D as the indicator function for set D. If Q is a vector or matrix, we represent its transpose as QT.

The main innovations and contributions of this paper are given as follows:• This paper introduces and investigates a novel stochastic epidemic model SEIuIrR with the incorporation of Ornstein–Uhlenbeck process to perturb the transmission rate δ.

• The model stochastic considered is more realistic and biologically meaningful framework for describing the transmission rate δ at which undetected people become detected, because this rate can be subject to fluctuations and uncertainties in real-world scenarios. For example, it can depend on factors like changes in testing capabilities, public health interventions, or population behavior. By modeling δ as a stochastic process, the model can account for these variations in detection rates over time.

• By employing novel Lyapunov functions, we establish the existence of a unique global solution of the model for any initial condition, we describe the sufficient conditions for establishing an ergodic stationary distribution and we proceed to define the sufficient criteria for extinction of the disease.

The rest of the paper is arranged as follows: In Section 2, we establish the existence of a unique global solution of system (2) for any initial condition. Section 3 outlines sufficient conditions for establishing a distinctive ergodic stationary distribution through the utilization of the stochastic Lyapunov method. In Section 4, we proceed to delineate the sufficient criteria for the extinction of the disease. In addition, numerical simulations are given in Section 5 to illustrate the results of the previous analysis. The last section concludes the paper.

2 The global solution's existence and uniqueness

Before delving into the properties of the epidemic system, it is crucial to establish whether the solution it exhibits is globally valid or not, and this theorem addresses the issue of the existence and uniqueness of the global solution for system (2) under any given initial value. Theorem 2.1 If there is an initial value(S(0),E(0),Iu(0),Ir(0),R(0),δ(0))T∈R+5×R, the system(2)has a unique solution(S(t),E(t),Iu(t),Ir(t),R(t),δ(t))T. That is,(S(t),E(t),Iu(t),Ir(t),R(t),δ(t))Tis defined for∀t≥0and remains inR+5×Ralmost surely (a.s.).

Proof 2.1 Since that all the coefficients of system (2) satisfy the local Lipschitz conditions, there will exist a unique local solution (S(t),E(t),Iu(t),Ir(t),R(t),δ(t))T on the interval [0,ϱe) for any initial value (S(0),E(0),Iu(0),Ir(0),R(0),δ(0))T∈R+5×R, where ϱe denotes the time of explosion [35]. So, to prove that this solution is global, it suffices to prove that: ϱe=∞ a.s. Letting Dn=(1/n,n)×(1/n,n)×(1/n,n)×(1/n,n)×(1/n,n)×(1/n,n), for any (S(0),E(0),Iu(0),Ir(0),R(0),δ(0))T∈R+5×R, one easily obtains a sufficiently large integer j0 to satisfy (S(0),E(0),Iu(0),Ir(0),R(0),exp⁡δ(0))∈Dj0. In this sense, for every integer j≥j0, a stopping time set ϱj is defined by [35]:ϱj=inf⁡{t∈[0,ϱe):min⁡{S(t),E(t),Iu(t),Ir(t),R(t),eδ(t)}≤1jormax⁡{S(t),E(t),Iu(t),Ir(t),R(t),eδ(t)}≥j}

Where we assume in our paper that inf⁡∅=∞ (with ∅ representing the empty set) and ϱ∞=limj→∞⁡ϱj, whence ϱ∞≤ϱe a.s. If ϱ∞=∞ a.s. is true, then ϱe=∞ a.s. and (S(t),E(t),Iu(t),Ir(t),R(t),δ(t))T∈R+5×R a.s. for all t≥0. So we must prove that ϱ∞=∞ a.s. To prove the latter, we use the proof backwards, assuming that it is false. It means there is a pair of constants (ε,T1)∈((0,1),R+) such that:P{ϱ∞≤T}>ε.

Consequently, there exists an integer j1≥j0 such that(3) P{ϱj≤T}:=P{Πj}≥ε∀j≥j1.

Define the function U:R+5×R→R+ byU(S,E,Iu,Ir,R,δ)=(S−1−ln⁡S)+(E−1−ln⁡E)+(Iu−1−ln⁡Iu)+(Ir−1−ln⁡Ir)+(R−1−ln⁡R)+δ22.

Applying the Itô's formula [35], we obtain(4) dU(S,E,Iu,Ir,R,δ)=LU(S,E,Iu,Ir,R,δ)dt+σ1δdB(t).

Where LU:R+5×R→R+ is defined by(5) LU=Δ+5μ+σ+d1+d2+γIu+γIr+σ122+max⁡{δ(t),0}+β(v1Iu+v2Ir)−μS−μE−μIu−d1Iu−μIr−d2Ir−μR−ΔS−β(v1Iu+v2Ir)SE−σ(1−ρ)EIu−σρEIr−max⁡{δ(t),0}IuIr−γIuIuR−γIrIrR+ρ1δ¯δ−ρ1δ2≤Δ+5μ+σ+d1+d2+γIu+γIr+σ122+|δ|+β(v1Iu+v2Ir)+ρ1δ¯δ−ρ1δ2.

Based on model (2) we haved(S+E+Iu+Ir+R)dt=Δ−μ(S+E+Iu+Ir+R)−(d1Iu+d2Ir)≤Δ−μ(S+E+Iu+Ir+R).

Consequently, this suggests that(6) S(t)+E(t)+Iu(t)+Ir(t)+R(t)≤{S(0)+E(0)+Iu(0)+Ir(0)+R(0), if S(0)+E(0)+Iu(0)+Ir(0)+R(0)≥Δμ,Δμ, if S(0)+E(0)+Iu(0)+Ir(0)+R(0)<Δμ≤ϒ1.

Whereϒ1:=sup⁡{S(0)+E(0)+Iu(0)+Ir(0)+R(0),Δμ}.

Substituting (6) into (5) leads to that(7) LU≤Δ+5μ+σ+d1+d2+γIu+γIr+σ122+β(v1+v2)ϒ1+supδ∈R⁡{−ρ1δ2+|δ|+ρ1δ¯δ}:=ϒ2.

In this case, ϒ2 represents a positive constant which is independent of the variables S,E,Iu,Ir,R and δ.

Consequently, by inetegrating both sides of equation (4) from 0 to ϱj∧T for any j≥j1, and then taking the mathematical expectations lead to(8) 0≤EU(S(ϱj∧T),E(ϱj∧T),Iu(ϱj∧T),Ir(ϱj∧T),R(ϱj∧T),δ(ϱj∧T))=EU(S(0),E(0),Iu(0),Ir(0),R(0),δ(0))+E∫0ϱj∧TLU(S(ξ),E(ξ),Iu(ξ),Ir(ξ),R(ξ),δ(ξ))dξ≤EU(S(0),E(0),Iu(0),Ir(0),R(0),δ(0))+ϒ2T.

Similarly, based on (3), it can be deduced that for all ζ∈Πj, at least one of the variables S,E,Iu,Ir,R and δ must be either 1j or j, ensuring that U(S(ϱj,ζ),E(ϱj,ζ),Iu(ϱj,ζ),Ir(ϱj,ζ),R(ϱj,ζ),δ(ϱj,ζ)) will not be less than(j−1−ln⁡j)∧ln2⁡j2 or (1j−1+ln⁡j)∧ln2⁡j2.

Based on (3) and (8), we have(9) EU(S(0),E(0),Iu(0),Ir(0),R(0),δ(0))+ϒ2T≥EU(S(ϱj∧T),E(ϱj∧T),Iu(ϱj∧T),Ir(ϱj∧T),R(ϱj∧T),δ(ϱj∧T))≥E[1Πj(ζ)U(S(ϱj∧T),E(ϱj∧T),Iu(ϱj∧T),Ir(ϱj∧T),R(ϱj∧T),δ(ϱj∧T))]≥P(Πj(ζ))U(S(ϱj,ζ),E(ϱj,ζ),Iu(ϱj,ζ),Ir(ϱj,ζ),R(ϱj,ζ),δ(ϱj,ζ))≥ε[(j−1−ln⁡j)∧(1j−1+ln⁡j)∧ln2⁡j2].

Considering the unlimited nature of j, as j approaches positive infinity, a contradiction appears.+∞<EU(S(0),E(0),Iu(0),Ir(0),R(0),δ(0))+ϒ2T=+∞.

That is to say ϱ∞=∞ a.s. The demonstration is finished.

Remark 2.1 From the conditions provided in the proof of Theorem 2.1, we can deduce that if S(0)+E(0)+Iu(0)+Ir(0)+R(0)<ΔμΣ={(S,E,Iu,Ir,R,δ)∈R+5×R,0≤S+E+Iu+Ir+R≤Δμ}

is positively invariant for system (2).

3 Ergodic stationary distribution

In the present section, our primary focus is on establishing adequate requirements for the existence of a stationary ergodic distribution, which, in turn, indicates the significant persistence of the susceptible population, exposed individuals, undetected infected individuals, and detected infected individuals. To achieve this objective, we introduce the following theorem.

Let us consider a b-dimensions nonlinear stochastic differential equation (SDE):(10) dV(t)=h1(V(t))dt+h2(V(t))dB(t).

With the initial value V(0)∈Rb, where B(t) is a b-dimensional Brownian motion defined on the complete probability space (Ω,F,{Ft}t≥0,P). Moreover, h1:Rb→Rb and h2:Rb→Rb×n are Borel measurable.

The following lemma is useful to prove the next theorem. Lemma 3.1 (See[14], Theorem 2.2) Assume the presence of a bounded closed domainA⊂Rbwith a regular boundary Λ. For any initial valueV(0)∈Rb, ifliminft→∞1t∫0tP(s,V(s),A)ds>0a.s.

WhereP(s,V(s),.)signifies the transition probability ofV(t), then a solution exists for system(10)that possesses the Feller property. Moreover, system(10)accommodates at least one stationary distributionΘ(.)onRb.

Theorem 3.1 Let R0s=S0(σ+μ)(σβv1(1−ρ)(δ¯+γIu+μ+d1)+σβν2ρ(γIr+μ+d2)+σβv2(1−ρ)(1π∫−δ¯ρ1σ1∞(σ1ρ1y+δ¯)14e−y2dy)4(δ¯+γIu+μ+d1)(γIr+μ+d2))−(σβv1(1−ρ)Δμ(μ+d1+δ¯+γIu)2+σβv2(1−ρ)Δ(1π∫−δ¯ρ1σ1∞(σ1ρ1y+δ¯)14e−y2dy)4μ(μ+d1+δ¯+γIu)2(μ+d2+γIr))σ12πρ1(σ+μ)>1 . Suppose the initial values (S(0),E(0),Iu(0),Ir(0),R(0),δ(0))∈R+5×R then, the solution (S(t),E(t),Iu(t),Ir(t),R(t),δ(t)) of system (2) possesses a single stationary distribution Θ(.) with the ergodic property.

Proof 3.1 Establish the Lyapunov function in the following manner:W¯(S,E,Iu,Ir,R,δ)=N˜W0(S,E,Iu,Ir,R)+W1(S,E,Iu,Ir)+W2(S,E,Iu,Ir,R)+W3(δ)

WhereW0(S,E,Iu,Ir,R,δ)=−2ln⁡E−(a1+a2+a3)ln⁡S−(a4+a5)ln⁡Iu−(a6+a7)ln⁡Ir+(a1+a2+a3)βR,W1(S,E,Iu,Ir)=−ln⁡S−ln⁡E−ln⁡Iu−ln⁡Ir,W2(S,E,Iu,Ir,R)=−ln⁡(Δμ−S−E−Iu−Ir−R),W3(δ)=δ22.

The positive constants a1, a2, a3, a4, a5, a6 and a7 will be determined later, while N˜ is a suitably large positive constant satisfying the given condition(11) −N˜(σ+μ)(R0s−1)+ϒ3≤−2

and(12) ϒ3:=supδ∈R⁡{−ρ12δ2+|δ|+ρ1δ¯δ}+βΔ(v1+v2)μ+6μ+2σ+d2+d2+γIu+γIr+σ122<∞.

Indeed, W¯(S,E,Iu,Ir,R,δ) exhibits continuity, conforming toliminfm→∞,(S,E,Iu,Ir,R,δ)∈Σ﹨DmW¯(S,E,Iu,Ir,R,δ)=+∞.

Consequently, a non-negative C2-function W(S,E,Iu,Ir,R,δ) is provided as followsW(S,E,Iu,Ir,R,δ)=W¯(S,E,Iu,Ir,R,δ)−W¯(S0,E0,Iu0,Ir0,R0,δ0).

Where (S0,E0,Iu0,Ir0,R0,δ0)∈Σ is the minimum point of W¯(S,E,Iu,Ir,R,δ).

By Applying Itô's formula [35] to W0 and the arithmetic-geometric men inequalityb1+b2+...+bnn≥(b1b2...bn)1n For any bn>0and n∈N,

we obtain(13) LW0=−(a1+a2+a3S)(Δ−β(v1Iu+v2Ir)S−μS)−2E(β(v1Iu+v2Ir)S−(σ+μ)E)−(a4+a5Iu)(σ(1−ρ)E−(μ+d1+max⁡{δ(t),0}+γIu)Iu)(a6+a7Ir)(σρE+max⁡{δ(t),0}Iu−(μ+d2+γIr)Ir)−(a1+a2+a3)β(γIuIu+γIrIr−μR)≤(−a1ΔS−a4σ(1−ρ)EIu−βv1IuSE)+a1μ+a4(μ+d1+max⁡{δ(t),0}+γIu)+(−a2ΔS−a6σρEIr−βv2IrSE)+a2μ+a6(μ+d2+γIr)+(−a3ΔS−a5σ(1−ρ)EIu−a7max⁡{δ,0}IuIr−βv2IrSE)+a3μ+a5(μ+d1+max⁡{δ(t),0}+γIu)+a7(μ+d2+γIr)+(a1+a2+a3)β(v1Iu+v2Ir)+(a1+a2+a3)β(γIuIu+γIrIr)+2(σ+μ)≤−3a1a4σβv1(1−ρ)Δ3+a1μ+a4(μ+d1+δ¯+γIu)−3a2a6σβv2ρΔ3+a2μ+a6(μ+d2+γIr)−4a3a5a7σβv2(1−ρ)Δδ˜4+a3μ+a5(μ+d1+δ¯+γIu)+a7(μ+d2+γIr)+(a4+a5)A(t)∨0+2(σ+μ)+(a1+a2+a3)β(v1Iu+v2Ir)+(a1+a2+a3)β(γIuIu+γIrIr),

withδ˜=max⁡{δ(t),0}andA(t)=δ(t)−δ¯.

Concerning the equation corresponding to the sixth position in system (2), that isdδ(t)=ρ1[δ¯−δ(t)]dt+σ1dB(t).

Based on the references [11], [31], [48], [51] it can be inferred that the process δ(t) exhibits the ergodic property and it is expected to undergo weak convergence towards the invariant densityf(y)=ρ1πσ1e−ρ1(y−δ¯)2σ12,y∈R.

Having incorporated the ergodic theorem [30], the aforementioned leads us to the following conclusion(14) ∫−∞∞(y∨0)14f(y)dy=∫0∞y14f(y)dy=∫0∞y14ρ1πσ1e−ρ1(y−δ¯)2σ12dy=1π∫−δ¯ρ1σ1∞(σ1ρ1y+δ¯)14e−y2dy.

Likewise, in the case of the stochastic differential equationdA(t)=ρ1A(t)dt+σ1dB(t).

The ergodic property of A(t) and its eventual weak convergence to the invariant density can be readily derivedg(y)=ρ1πσ1e−ρ1y2σ12,y∈R.

Drawing upon the ergodic theorem [30], we arrive at the following conclusion(15) ∫−∞∞(y∨0)g(y)dy=∫0∞yg(y)dy=∫0∞yρ1πσ1e−ρ1y2σ12dy=σ12πρ1.

Substituting (14) and (15) into (13) leads to that(16) LW0≤−3a1a4σβv1(1−ρ)Δ3+a1μ+a4(μ+d1+δ¯+γIu)−3a2a6σβv2ρΔ3+a2μ+a6(μ+d2+γIr)−4a3a5a7σβv2(1−ρ)Δδˆ4+(4a3a5a7σβv2(1−ρ)Δδˆ4−4a3a5a7σβv2(1−ρ)Δδ˜4)+a3μ+a5(μ+d1+δ¯+γIu)+a7(μ+d2+γIr)+2(σ+μ)+(a4+a5)σ12πρ1+(a4+a5)(A(t)∨0−∫0∞yg(y)dy)+(a1+a2+a3)β(v1Iu+v2Ir)+(a1+a2+a3)β(γIuIu+γIrIr).

Whereδˆ=(∫−∞∞(y∨0)14f(y)dy)4=(1π∫−δρ1σ1∞(σ1ρ1y+δ¯)14e−y2dy)4.

Leta1μ=a4(μ+d1+δ¯+γIu)=σβv1(1−ρ)Δμ(μ+d1+δ¯+γIu),a2μ=a6(μ+d2+γIr)=σβv2ρΔμ(μ+d2+γIr)a3μ=a5(μ+d1+δ¯+γIu)=a7(μ+d2+γIr)=σβv2(1−ρ)Δδˆμ(μ+d1+δ¯+γIu)(μ+d2+γIr).

Consequently, we acquirea1=σβv1(1−ρ)Δμ2(μ+d1+δ¯+γIu),a4=σβv1(1−ρ)Δμ(μ+d1+δ¯+γIu)2a2=σβv2ρΔμ2(μ+d2+γIr),a6=σβv2ρΔμ(μ+d2+γIr)2a3=σβv2(1−ρ)Δδˆμ2(μ+d1+δ¯+γIu)(μ+d2+γIr),a5=σβv2(1−ρ)Δδˆμ(μ+d1+δ¯+γIu)2(μ+d2+γIr),a7=σβv2(1−ρ)Δδˆμ(μ+d1+δ¯+γIu)(μ+d2+γIr)2.

Consequently, as a result(17) LW0≤−σβv1(1−ρ)Δμ(μ+d1+δ¯+γIu)−σβv2ρΔμ(μ+d2+γIr)−σβv2(1−ρ)Δδˆμ(μ+d1+δ¯+γIu)(μ+d2+γIr)+2(σ+μ)+(a4+a5)σ12πρ1+4a3a5a7σβv2(1−ρ)Δ4(δˆ4−δ˜4)+(a4+a5)(A(t)∨0−∫0∞yg(y)dy)+(a1+a2+a3)β(v1Iu+v2Ir)+(a1+a2+a3)β(γIuIu+γIrIr)≤−(σ+μ)(R0s−1)+σ+μ+(a4+a5)(A(t)∨0−∫0∞yg(y)dy)+(a1+a2+a3)β(v1Iu+v2Ir)+(a1+a2+a3)β(γIuIu+γIrIr).

WhereR0s=S0(σ+μ)(σβv1(1−ρ)(δ¯+γIu+μ+d1)+σβν2ρ(γIr+μ+d2)+σβv2(1−ρ)δˆ(δ¯+γIu+μ+d1)(γIr+μ+d2))−(a1+a5)σ12πρ1(σ+μ).

Analogously, employing Itô's formula [35] to W1, W2, and W3, respectively, yields the following results(18) LW1=−ΔS+β(v1Iu+v2Ir)+μ−β(v1Iu+v2Ir)SE+σ+μ−σ(1−ρ)EIu+(μ+d1+max⁡{δ(t),0}+γIu)−σρEIr−max⁡{δ(t),0}IuIr+(μ+d2+γIr)≤−ΔS−β(v1Iu+v2Ir)SE+βΔ(v1+v2)μ+|δ|+4μ+σ+d1+d2+γIu+γIr

(19) LW2=1Δμ−S−E−Iu−Ir−R(Δ−μS−μIu−μIr−d1Iu−d2Ir−μR)≤μ−d1Iu+d2IrΔμ−S−E−Iu−Ir−R

and(20) LW3=ρ1δ¯δ−ρ1δ2+σ122.

By (17), (18), (19) and (20), we get(21) LW≤−N˜(σ+μ)(R0s−1)+4a3a5a7σβv2(1−ρ)Δ4(δˆ4−δ˜4)+N˜(a1+a2+a3)β(v1Iu+v2Ir)+N˜(a1+a2+a3)β(γIuIu+γIrIr)+N˜(a4+a5)(A(t)∨0−∫0∞yg(y)dy)−ΔS−β(v1Iu+v2Ir)SE−d1Iu+d2IrΔμ−S−E−Iu−Ir−R+ρ1δ¯δ−ρ12δ2−ρ12δ2+|δ|+βΔ(v1+v2)μ+6μ+2σ+d1+d2+γIu+γIr+σ122:=H(S,E,Iu,Ir,R;δ)+4a3a5a7σβv2(1−ρ)Δ4(δˆ4−δ˜4)+N˜(a4+a5)(A(t)∨0−∫0∞yg(y)dy).

Where(22) H(S,E,Iu,Ir,R;δ)=−N˜(σ+μ)(R0s−1)+N˜(a1+a2+a3)β(v1Iu+v2Ir)+N˜(a1+a2+a3)β(γIuIu+γIrIr)−ΔS−β(v1Iu+v2Ir)SE−d1Iu+d2IrΔμ−S−E−Iu−Ir−R+ρ1δ¯δ−ρ12δ2−ρ12δ2+|δ|+βΔ(v1+v2)μ+6μ+2σ+d1+d2+γIu+γIr+σ122.

Following that, we define a suitable bounded subset Eˆ is described byEˆ={(S,E,Iu,Ir,R)T∈Σ,S≥ε,Iu≥ε,Ir≥ε,E≥ε3,S+E+Iu+Ir+R≤Δμ−ε2,|δ|≤1ε}.

Let ε be a sufficiently small positive constant that fulfills the following conditions.(23) −Δε+ϒ4≤−1

(24) ε≤1N˜β(a1+a2+a3)(v1+v2+γIu+γIr)

(25) −β(v1+v2)ε+ϒ4≤−1

(26) −d1+d2ε+ϒ4≤−1

(27) −ρ12ε2+ϒ4≤−1.

Where(28) ϒ4:=supδ∈R⁡{−ρ12δ2+|δ|+ρ1δ¯δ}+N˜(a1+a2+a3)β(v1Iu+v2Ir)+N˜(a1+a2+a3)β(γIuIu+γIrIr)+βΔ(v1+v2)μ+6μ+2σ+d1+d2+γIu+γIr+σ122.

We can then partition the set Σ﹨Eˆ into the following five subsets Eˆic,i=1,2,3,4,5, whereEˆ1c={(S,E,Iu,Ir,R,δ)T∈Σ,S<ε}

Eˆ2c={(S,E,Iu,Ir,R,δ)T∈Σ,Iu<ε,Ir<ε}

Eˆ3c={(S,E,Iu,Ir,R,δ)T∈Σ,E<ε3,Iu≥ε,Ir≥ε}

Eˆ4c={(S,E,Iu,Ir,R,δ)T∈Σ,S+E+Iu+Ir+R>Δμ−ε2,Iu≥ε,Ir≥ε}

Eˆ5c={(S,E,Iu,Ir,R,δ)T∈Σ,|δ|>1ε}.

It is then obvious that,Σ﹨Eˆ=⋃i=15Eˆic.

We shall now proceed to establish the proof that H(S,E,Iu,Ir,R;δ)≤−1 on Eˆc.

Put differently, we must verify its fulfillment across the aforementioned five regions• Case 1. For any (S,E,Iu,Ir,R,δ)T∈Eˆ1c, according to (22), we get(29) H(S,E,Iu,Ir,R,δ)≤−ΔS+N˜(a1+a2+a3)β(v1Iu+v2Ir)+N˜(a1+a2+a3)β(γIuIu+γIrIr)+ρ1δ¯δ−ρ12δ2+|δ|+βΔ(v1+v2)μ+6μ+2σ+d1+d2+γIu+γIr+σ122≤−ΔS+N˜(a1+a2+a3)β(v1+v2)Δμ+N˜(a1+a2+a3)β(γIu+γIr)Δμ+ρ1δ¯δ−ρ12δ2+|δ|+βΔ(v1+v2)μ+6μ+2σ+d1+d2+γIu+γIr+σ122≤−Δε+ϒ4≤−1.

Which follows from (23) and (28)

• Case 2. For any (S,E,Iu,Ir,R,δ)T∈Eˆ2c, by (22), we obtain(30) H(S,E,Iu,Ir,R,δ)≤−N˜(σ+μ)(R0s−1)+N˜(a1+a2+a3)β(v1Iu+v2Ir)+N˜(a1+a2+a3)β(γIuIu+γIrIr)+ρ1δ¯δ−ρ12δ2+|δ|+βΔ(v1+v2)μ+6μ+2σ+d1+d2+γIu+γIr+σ122≤−N˜(σ+μ)(R0s−1)+ϒ3+N˜β(a1+a2+a3)(v1+v2+γIu+γIr)ε≤−1.

Which follows from (11) and (24).

• Case 3. For any (S,E,Iu,Ir,R,δ)T∈Eˆ3c, from (22) it follows that(31) H(S,E,Iu,Ir,R,δ)≤−β(v1Iu+v2Ir)SE+N˜(a1+a2+a3)β(v1Iu+v2Ir)+N˜(a1+a2+a3)β(γIuIu+γIrIr)+ρ1δ¯δ−ρ12δ2+|δ|+βΔ(v1+v2)μ+6μ+2σ+d1+d2+γIu+γIr+σ122≤−β(v1Iu+v2Ir)SE+N˜(a1+a2+a3)β(v1+v2)Δμ+N˜(a1+a2+a3)β(γIu+γIr)Δμ+ρ1δ¯δ−ρ12δ2+|δ|+βΔ(v1+v2)μ+6μ+2σ+d1+d2+γIu+γIr+σ122≤−(v1+v2)ε2ε3+ϒ4≤−(v1+v2)ε+ϒ4≤−1.

Which results from (25).

• Case 4. For any (S,E,Iu,Ir,R,δ)T∈Eˆ4c, in view of (22), we obtain(32) H(S,E,Iu,Ir,R,δ)≤−d1Iu+d2IrΔμ−S−E−Iu−Ir−R+N˜(a1+a2+a3)β(v1Iu+v2Ir)+N˜(a1+a2+a3)β(γIuIu+γIrIr)+ρ1δ¯δ−ρ12δ2+|δ|+βΔ(v1+v2)μ+6μ+2σ+d1+d2+γIu+γIr+σ122≤−d1Iu+d2IrΔμ−S−E−Iu−Ir−R+N˜(a1+a2+a3)β(v1+v2)Δμ+N˜(a1+a2+a3)β(γIu+γIr)Δμ+ρ1δ¯δ−ρ12δ2+|δ|+βΔ(v1+v2)μ+6μ+2σ+d1+d2+γIu+γIr+σ122≤−(d1+d2)εε2+ϒ4≤−(d1+d2)ε+ϒ4≤−1.

Which results from (26).

• Case 5. For any (S,E,Iu,Ir,R,δ)T∈Eˆ5c, in view of (22), we have(33) H(S,E,Iu,Ir,R,δ)≤−ρ12δ2+N˜(a1+a2+a3)β(v1Iu+v2Ir)+N˜(a1+a2+a3)β(γIuIu+γIrIr)+ρ1δ¯δ−ρ12δ2+|δ|+βΔ(v1+v2)μ+6μ+2σ+d1+d2+γIu+γIr+σ122≤−ρ12δ2+N˜(a1+a2+a3)β(v1+v2)Δμ+N˜(a1+a2+a3)β(γIu+γIr)Δμ+ρ1δ¯δ−ρ12δ2+|δ|+βΔ(v1+v2)μ+6μ+2σ+d1+d2+γIu+γIr+σ122≤−ρ12ε2+ϒ4≤−1.

Which follows from (27).

Drawing from the evidence presented in inequalities (29), (30), (31), (32) and (33), a straightforward conclusion can be reached, establishing the existence of a sufficiently small ε that satisfies The following condition(34) H(S,E,Iu,Ir,R,δ)≤−1for any(S,E,Iu,Ir,R,δ)T∈Eˆc.

So we have(35) H(S,E,Iu,Ir,R,δ)≤ϒ5<∞for any(S,E,Iu,Ir,R,δ)T∈R+5×R.

Whereϒ5:=sup(S,E,Iu,Ir,R,δ)∈R+5×R⁡{−N˜(σ+μ)(R0s−1)+N˜(a1+a2+a3)β(v1Iu+v2Ir)+N˜(a1+a2+a3)β(γIuIu+γIrIr)−ΔS−β(v1Iu+v2Ir)SE−d1Iu+d2IrΔμ−S−E−Iu−Ir−R+ρ1δ¯δρ12δ2−ρ12δ2+δ+βΔ(v1+v2)μ+6μ+2σ+d1+d2+γIu+γIr+σ122}.

By integrating equation (21) over the interval [0,t] for any initial values (S(0),E(0),Iu(0),Ir(0),R(0),δ(0))∈Σ and subsequently calculating the mathematical expectation, the following expression is obtained(36) 0≤EW(S(t),E(t),Iu(t),Ir(t),R(t),δ(t))t=EW(S(0),E(0),Iu(0),Ir(0),R(0),δ(0))t+1t∫0tE(LW(S(s),E(s),Iu(s),Ir(s),R(s),δ(s)))ds≤EW(S(0),E(0),Iu(0),Ir(0),R(0),δ(0))t+1t∫0tE(H(S(s),E(s),Iu(s),Ir(s),R(s),δ(s)))ds+4N˜a3a5a7σβv2(1−ρ)Δ4E[∫−∞∞(y∨0)14f(y)dy−1t∫0t(δ(s)∨0)14ds]+N˜E[1t∫0t(A(s)∨0)ds−∫0∞yg(y)dy].

Utilizing the ergodic properties of both δ(t) and A(t), along with the strong law of large numbers as demonstrated in [35], we derive the followinglimt→∞⁡E[∫−∞∞(y∨0)14f(y)dy−1t∫0t(δ(s)∨0)14ds]=0a.s.

andlimt→∞⁡E[1t∫0t(A(s)∨0)ds−∫0∞yg(y)dy]=0a.s.

Therefore, by applying the inferior limit to both sides of (36), we deduce the subsequent outcome(37) 0≤liminft→∞EW(S(0),E(0),Iu(0),Ir(0),R(0),δ(0))t+liminft→∞1t∫0tE(H(S(s),E(s),Iu(s),Ir(s),R(s),δ(s)))ds=liminft→∞1t∫0tE(H(S(s),E(s),Iu(s),Ir(s),R(s),δ(s)))dsa.s.

Furthermore, in accordance with to (34) and (35), it follows that we obtain(38) liminft→∞1t∫0tE(H(S(s),E(s),Iu(s),Ir(s),R(s),δ(s)))ds=liminft→∞1t∫0tE(H(S(s),E(s),Iu(s),Ir(s),R(s),δ(s)))1(S(s),E(s),Iu(s),Ir(s),R(s),δ(s))T∈Eˆds+liminft→∞1t∫0tE(H(S(s),E(s),Iu(s),Ir(s),R(s),δ(s)))1(S(s),E(s),Iu(s),Ir(s),R(s),δ(s))T∈(Σ﹨Eˆ)ds≤ϒ5liminft→∞1t∫0t1(S(s),E(s),Iu(s),Ir(s),R(s),δ(s))T∈Eˆds−liminft→∞1t∫0t1(S(s),E(s),Iu(s),Ir(s),R(s),δ(s))T∈(Σ﹨Eˆ)ds≤−1+(ϒ5+1)liminft→∞1t∫0t1(S(s),E(s),Iu(s),Ir(s),R(s),δ(s))T∈Eˆds.

Consequently, from (37) and (38) lead to the conclusion that(39) liminft→∞1t∫0t1(S(s),E(s),Iu(s),Ir(s),R(s),δ(s))T∈Eˆds≥1ϒ5+1>0a.s.

Considering the event probability definition and the application of Fatou's lemma [12], [31], [51]. The result (39) is equivalent to(40) liminft→∞1t∫0tP(s,S(s),E(s),Iu(s),Ir(s),R(s),δ(s),Eˆ)ds≥1ϒ5+1>0a.s.

Where P(s,S(s),E(s),Iu(s),Ir(s),R(s),δ(s),Eˆ) is the transition probability of (S(s),E(s),Iu(s),Ir(s),R(s),δ(s))T belonging to set Eˆ. Thus, we have fulfilled the conditions of Lemma 3.1 and thus the system (2) has at least one stationary distribution Θ(.) on R+5×R which has the Feller property. This completes the proof.

4 Extinction of the disease

Within this portion, we will set forth the sufficient conditions for the complete eradication of the disease. Theorem 4.1 Consider the solution(S(t),E(t),Iu(t),Ir(t),R(t),δ(t))to system(2)with initial value(S(0),E(0),Iu(0),Ir(0),R(0),δ(0))∈R+5×R. IfR0EX=R32+(R32)2+(R1+R23)33+R32−(R32)2+(R1+R23)33<1

andψ=min⁡{σ+μ,μ+d1+δ¯+γIu,μ+d2+γIr}(R0EX−1)+(v1(μ+d1+δ¯+γIu)(μ+d2+γIr)REXR0EX(μ+d2+γIr)S0v1+δ¯S0v2+(μ+d2+γIr)R0EXS0)∫0∞∫−∞∞|l−s0|Θ(l,δ)dldδ+(v2(μ+d1δ¯+γIu)R0EXv1(μ+d2+γIr)R0EX+δ¯v2+1)σ1πρ1is negative.

Thenlimt→∞⁡E(t)=0,limt→∞⁡Iu(t)=0,limt→∞⁡Ir(t)=0a.s.

Specifically, the disease undergoes exponential extinction with a almost surely.

Proof 4.1 Based on the first equation of system (2), we have:dS(t)≤(Δ−μS)dt

Let the following auxiliary logistical equation be(41) dL(t)=(Δ−μL)dt.

Assume that L(t) represents the solution of (41) with the initial condition L(0)=S(0)>0. Referring to the theorem shown in [25], we obtain S(t)≤L(t), for any t≥0 a.s.

Furthermore, B it is readily apparent that a three-dimensional matrix possesses a non-negative eigenvector on the left and R0EX(ϑ1,ϑ2,ϑ3)=(ϑ1,ϑ2,ϑ3)B andB=(0βS0v1σ+μβS0v2σ+μσ(1−ρ)μ+d1+δ¯+γIu00σρμ+d2+γIrδ¯μ+d2+γIr0)

Define a C2-lyapunov function Ģ(E,Iu,Ir) byĢ(E,Iu,Ir)=α1E+α2Iu+α3Ir.

Whereα1=ϑ1σ+μ,α2=ϑ2μ+d1+δ¯+γIu,α3=ϑ3μ+d2+γIr.

Which implies(42) L(ln⁡Ģ)=1Ģ[α1(β(v1Iu+v2Ir)S−(σ+μ)E)+α2(σ(1−ρ)E−(μ+d1+δ˜+γIu)Iu)+α3(σρE+δ˜Iu−(μ+d2+γIr)Ir)]=1Ģ[α1(β(v1Iu+v2Ir)S0−(σ+μ)E)+α2(σ(1−ρ)E−(μ+d1+δ¯+γIu)Iu)+α3(σρE+δ¯Iu−(μ+d2+γIr)Ir)]+α1(β(v1Iu+v2Ir)(S−S0)Ģ+α2(δ¯−δ˜)IuĢ+α3(δ˜−δ¯)IuĢ≤1Ģ[ϑ1σ+μ(β(v1Iu+v2Ir)S0−(σ+μ)E)+ϑ2μ+d1+δ¯+γIu(σ(1−ρ)E−(μ+d1+δ¯+γIu)Iu)+ϑ3μ+d2+γIr(σρE+δ¯Iu−(μ+d2+γIr)Ir)]+α1(β(v1Iu+v2Ir)(L−S0)Ģ+α2|δ˜−δ¯|IuĢ+α3|δ˜−δ¯|IuĢ≤1Ģ(ϑ1,ϑ2,ϑ3)(B(E,Iu,Ir)T−(E,Iu,Ir)T)+α1(β(v1Iu+v2Ir)(L−S0)Ģ+(α2+α3)|δ˜−δ¯|IuĢ≤1Ģ(R0EX−1)(ϑ1E+ϑ2Iu+ϑ3Ir)+α1β(α3v1+α2v2)α2α3|L−S0|+(α2+α3)α2|δ˜−δ¯|≤1Ģ(R0EX−1)(α1(σ+μ)E+α2(μ+d1δ¯+γIu)Iu+α3(μ+d2+γIr)Ir)+α1β(α3v1+α2v2)α2α3|L−S0|+(α2+α3)α2|δ˜−δ¯|≤min⁡{σ+μ,μ+d1+δ¯+γIu,μ+d2+γIr}(R0EX−1)+χ1|L−S0|+χ2|δ˜−δ¯|.

Whereχ1:=v1(μ+d1+δ¯+γIu)(μ+d2+γIr)REXREX(μ+d2+γIr)S0v1+δ¯S0v2+(μ+d2+γIr)R0EXS0

χ2=v2(μ+d1δ¯+γIu)R0EXv1(μ+d2+γIr)R0EX+δ¯v2+1.

Upon integrating both sides of (42) over the interval [0,t] and subsequently dividing by t, it follows that(43) ln⁡Ģ(E(t),Iu(t),Ir(t))t≤ln⁡Ģ(E(0),Iu(0),Ir(0))t+min⁡{σ+μ,μ+d1+δ¯+γIu,μ+d2+γIr}(R0EX−1)+χ11t∫0t|L(s)−S0|ds+χ21t∫0t|δ˜(s)−δ¯|ds.

Considering Theorem 3.1 along with the strong law of large numbers [30], [31], [35], [51], we can conclude that the process (L(t),δ(t)) possesses a distinct stationary distribution Θ(.,.) and exhibits the property of ergodicity so(44) limt→∞⁡1t∫0t|L(s)−S0|ds=∫0∞∫−∞∞|l−S0|Θ(l,δ)dldδ,

and(45) limt→∞⁡1t∫0t|δ˜(s)−δ¯|ds=∫−∞∞|δ˜(y)−δ¯|f(y)dy=σ1πρ1.

Applying the upper limit to both sides of (43) and consolidating it with (44) and (45), we arrive atlimsupt→∞ln⁡Ģ(E(t),Iu(t),Ir(t))t≤min⁡{σ+μ,μ+d1+δ¯+γIu,μ+d2+γIr}(R0EX−1)+χ1∫0∞∫−∞∞|l−S0|Θ(l,δ)dldδ+χ2σ1πρ1:=ψa.s.

and if ψ is negative, we can deduce from this:limsupt→∞ln⁡E(t)t<0,limsupt→∞ln⁡Iu(t)t<0andlimsupt→∞ln⁡Ir(t)t<0a.s.

This entails thatlimt→∞⁡E(t)=0,limt→∞⁡Iu(t)=0andlimt→∞⁡Ir(t)=0a.s.

This finalizes the theorem's proof.

5 Numerical results

Using the advanced technique outlined by Milstein in reference [22], the discretized equation corresponding to system (2) can be derived(46) {S[j+1]=S[j]+(Δ−β(v1Iu[j]+v2Ir[j])S[j]−μS[j])Δt,E[j+1]=E[j]+(β(v1Iu[j]+v2Ir[j])S[j]−(σ+μ)E[j])Δt,Iu[j+1]=Iu[j]+(σ(1−ρ)E[j]−(μ+d1+max⁡(δ[j],0)+γIu)Iu[j])Δt,Ir[j+1]=Ir[j]+(σρE[j]+max⁡(δ[j],0)Iu[j]−(μ+d2+γIr)Ir[j])Δt,R[j+1]=R[j]+(γIuIu[j]+γIrIr[j]−μR[j])Δt,δ[j+1]=δ[j]+ρ1(δ¯−δ[j])Δt+σ1Δtιj.

In accordance with the jth iteration of Equation (46), denoted as (S[j],E[j],Iu[j],Ir[j],R[j],δ[j]), where Δt represents the positive time increment, ιj signifies independent Gaussian random variables adhering to the N(0,1) distribution for j=1,...,n. Realistic parameter values are selected from established sources, and the biological parameters are comprehensively listed in Table 2.Table 2 Values of parameters of stochastic model (2).

Table 2Parameters	values	Source	
Δ	0.1	Assumed	
μ	0.0399	Assumed	
β	0.6594	[44]	
ρ	0.2929	[44]	
v1	0.3958	[18]	
v2	0.4941	[18]	
δ¯	0.1	Estimated	
d1	0.0290	[44]	
d2	0.4897	[44]	
σ	0.5732	[18]	
γIu	0.0458	[44]	
γIr	0.0806	[44]	

In this section, our primary focus is on confirming the validity of the following two outcomes:1. The condition for R0s>1 leads to the existence of a distinctive ergodic stationary distribution.

2. Model (2) undergoes exponential extinction when R0EX<1 and ψ<0.

Example 5.1 Given the values β=0.6594,ρ1=0.1 and σ1=0.02, a calculation yields

R0s=S0(σ+μ)(σβv1(1−ρ)(δ¯+γIu+μ+d1)+σβν2ρ(γIr+μ+d2)+σβv2(1−ρ)(1π∫−δ¯ρ1σ1∞(σ1ρ1y+δ¯)14e−y2dy)4(δ¯+γIu+μ+d1)(γIr+μ+d2))−(σβv1(1−ρ)Δμ(μ+d1+δ¯+γIu)2+σβv2(1−ρ)Δ(1π∫−δ¯ρ1σ1∞(σ1ρ1y+δ¯)14e−y2dy)4μ(μ+d1+δ¯+γIu)2(μ+d2+γIr))σ12πρ1(σ+μ)≈2.54697>1.

According to Theorem 3.1, we deduce the existence of a singular stationary distribution Θ(.) with ergodic properties, as depicted in Fig. 1.Figure 1 The left illustration depicts the simulated evolution of the solution (S(t),E(t),Iu(t),Ir(t),R(t)) for the models (1) and (2) using the given parameter values β = 0.6594,ρ1 = 0.1, and σ1 = 0.02. The right displays the frequency histogram and probability densities for S, E, Iu, Ir and R of model (2).

Figure 1

It's clear that the Fig. 1 illustrates a stationary distribution of a solution of the system, (S(t),E(t), Iu(t),Ir(t),R(t),δ(t))T, which means that the disease will last for a long time.

Example 5.2 Assume that β=0.2,ρ1=4 and σ1=0.08, then we similarly compute thatR0EX=R32+(R32)2+(R1+R23)33+R32−(R32)2+(R1+R23)33≈0.16<1

andψ=min⁡{σ+μ,μ+d1+δ¯+γIu,μ+d2+γIr}(R0EX−1)+(v1(μ+d1+δ¯+γIu)(μ+d2+γIr)REXREX(μ+d2+γIr)S0v1+δ¯S0v2+(μ+d2+γIr)R0EXS0)∫0∞∫−∞∞|l−S0|Θ(l,δ)dldδ+(v2(μ+d1δ¯+γIu)R0EXv1(μ+d2+γIr)R0EX+δ¯v2+1)σ1πρ1≈−0.003.

The implications of Theorem 4.1 readily demonstrate that the infected population will undergo exponential extinction, thereby guaranteeing the eradication of the disease outbreak. The graphs in Fig. 2 representing the curves of the populations S(t), E(t), Iu(t) and Ir(t) confirm this result.Figure 2 The simulations of the solution (S(t),E(t),Iu(t),Ir(t),R(t)) of the deterministic model (1) and the stochastic model (2) with the parameter values β = 0.2,ρ1 = 4, and σ1 = 0.08.

Figure 2

The Fig. 2 presents the extinction dynamics of the disease, showing a steady and continuous decline in the number of exposed E, undetected infected Iu, and infected Ir individuals over time. This trend ultimately leads to the disease's eradication.

6 Conclusion

In this paper we analyzed a novel stochastic SEIuIrR model with the Ornstein–Uhlenbeck process to describe the transmission rate from undetected to detected individuals. We proved the theoretical results by constructing a series of suitable Lyapunov functions. In the first, we gave the theoretical result that the stochastic SEIuIrR system (2) has a unique global positive solution and proved it. Then, we established sufficient criteria for the existence of stationary distribution and exposed the effects of the Ornstein–Uhlenbeck process on the existence of stationary distribution. Specifically, if R0s>1 and the parameters δ of the Ornstein–Uhlenbeck process meet certain conditions, the system (2) exists with a stable distribution. In addition, we derived the sufficient conditions when R0EX<1 for the extinction of the disease. We used numerical simulation to simulate and verify the theoretical results in the paper.

CRediT authorship contribution statement

Mhammed Mediani: Writing – original draft, Visualization, Methodology, Investigation. Abdeldjalil Slama: Writing – original draft, Validation, Conceptualization. Ahmed Boudaoui: Writing – review & editing, Writing – original draft, Validation, Conceptualization. Thabet Abdeljawad: Writing – review & editing, Supervision, Investigation, Formal analysis, Conceptualization.

Declaration of Competing Interest

The authors declare that there is no conflict of interest in this work.

Data availability

The data that supports the findings of this study are available within the article.

Acknowledgements

The author T. Abdeljawad would like to thank 10.13039/501100012639 Prince Sultan University for the support through TAS research lab.
==== Refs
References

1 Abbas S. Ahmad M. Nazar M. Ahmad Z. Amjad M. Garalleh H.A. Jan A.Z. Soret effect on mhd Casson fluid over an accelerated plate with the help of constant proportional Caputo fractional derivative ACS Omega 2024
2 Abbas S. Ahmad M. Nazar M. Amjad M. Ali H. Jan A.Z. Heat and mass transfer through a vertical channel for the Brinkman fluid using Prabhakar fractional derivative Appl. Therm. Eng. 232 2023 121065
3 Abbas S. Ahmad M. Rahimzai A.A. Nazar M. Khan I. Active and Passive Control of Mhd Jeffrey Nanofluid over a Vertical Plate with Constant Proportional Caputo Fractional Derivative 2023
4 Abbas S. Gilani S.F.F. Nazar M. Fatima M. Ahmad M. Nisa Z.U. Bio-convection flow of fractionalized second grade fluid through a vertical channel with Fourier's and Fick's laws Mod. Phys. Lett. B 37 23 2023 2350069
5 Abbas S. Nazar M. Nisa Z.U. Amjad M. Din S.M.E. Alanzi A.M. Heat and mass transfer analysis of mhd Jeffrey fluid over a vertical plate with cpc fractional derivative Symmetry 14 12 2022 2491
6 Abbas S. Nisa Z.U. Nazar M. Amjad M. Ali H. Jan A.Z. Application of heat and mass transfer to convective flow of Casson fluids in a microchannel with Caputo–Fabrizio derivative approach Arab. J. Sci. Eng. 49 1 2024 1275 1286
7 Adnan Thirthar A. Stability and bifurcation of an sis epidemic model with saturated incidence rate and treatment function Iran. J. Math. Sci. Inform. 15 2 2020 129 146
8 Allen E. Environmental variability and mean-reverting processes Discrete Contin. Dyn. Syst., Ser. B 21 7 2016 2073 2089
9 Atangana A. Modelling the spread of covid-19 with new fractal-fractional operators: can the lockdown save mankind before vaccination? Chaos Solitons Fractals 136 2020 109860
10 Britton T. Stochastic epidemic models: a survey Math. Biosci. 225 1 2010 24 35 20102724
11 Cai Y. Jiao J. Gui Z. Liu Y. Wang W. Environmental variability in a stochastic epidemic model Appl. Math. Comput. 329 2018 210 226
12 Chen S.-S. Cheng C.-Y. Takeuchi Y. Stability analysis in delayed within-host viral dynamics with both viral and cellular infections J. Math. Anal. Appl. 442 2 2016 642 672
13 Chu Y.-M. Zafar Z.U.A. Inc M. Javeed S. Ali A.S. Numerical modeling of a novel stochastic coronavirus Fractals 30 08 2022 2240211
14 Du N.H. Nguyen D.H. Yin G.G. Conditions for permanence and ergodicity of certain stochastic predator–prey models J. Appl. Probab. 53 1 2016 187 202
15 El hadj Moussa Y. Boudaoui A. Ullah S. Muzammil K. Riaz M.B. Application of fractional optimal control theory for the mitigating of novel coronavirus in Algeria Results Phys. 39 2022 105651
16 Ez-Zetouni A. Khyar O. Allali K. Akdim K. Zahid M. Stochastic and Deterministic Analysis of a covid-19 Pandemic Model Under Vaccination Strategy: Real Cases Application 2022
17 Garnett G.P. An introduction to mathematical models in sexually transmitted disease epidemiology Sex. Transm. Infect. 78 1 2002 7 12 11872850
18 Gatyeni S. Chukwu C. Chirove F. Nyabadza F. Application of optimal control to the dynamics of covid-19 disease in South Africa Sci. Afr. 16 2022 e01268
19 Gul H. Alrabaiah H. Ali S. Shah K. Muhammad S. Computation of solution to fractional order partial reaction diffusion equations J. Adv. Res. 25 2020 31 38 32922971
20 Guo Y. Li T. Modeling the competitive transmission of the omicron strain and delta strain of covid-19 J. Math. Anal. Appl. 526 2 2023 127283
21 Han B. Jiang D. Zhou B. Hayat T. Alsaedi A. Stationary distribution and probability density function of a stochastic SIRSI epidemic model with saturation incidence rate and logistic growth Chaos Solitons Fractals 142 2021 110519
22 Higham D.J. An algorithmic introduction to numerical simulation of stochastic differential equations SIAM Rev. 43 3 2001 525 546
23 Hussain G. Khan A. Zahri M. Zaman G. Ergodic stationary distribution of stochastic epidemic model for hbv with double saturated incidence rates and vaccination Chaos Solitons Fractals 160 2022 112195
24 Hyman J.M. Stanley E.A. Using mathematical models to understand the aids epidemic Math. Biosci. 90 1–2 1988 415 473
25 Ikeda N. Watanabe S. Stochastic Differential Equations and Diffusion Processes 1981 North-Holland Publishing Company New York
26 Jawad S. Winter M. Rahman Z.-A.S. Al-Yasir Y.I. Zeb A. Dynamical behavior of a cancer growth model with chemotherapy and boosting of the immune system Mathematics 11 2 2023 406
27 Khan M.A. Atangana A. Modeling the dynamics of novel coronavirus (2019-ncov) with fractional derivative Alex. Eng. J. 59 4 2020 2379 2389
28 Khan S.A. Shah K. Kumam P. Seadawy A. Zaman G. Shah Z. Study of mathematical model of hepatitis b under Caputo-Fabrizo derivative AIMS Math. 6 1 2021 195 209
29 Kifle Z.S. Lemecha Obsu L. Optimal control analysis of a covid-19 model Appl. Math. Sci. Eng. 31 1 2023 2173188
30 Kutoyants A.Y. Statistical Inference for Ergodic Diffusion Processes 2003 Springer London
31 Liu Q. Stationary distribution and extinction of a stochastic hliv model with viral production and Ornstein–Uhlenbeck process Commun. Nonlinear Sci. Numer. Simul. 119 2023 107111
32 Liu Q. Jiang D. Stationary distribution and extinction of a stochastic sir model with nonlinear perturbation Appl. Math. Lett. 73 2017 8 15
33 Liu Q. Jiang D. Analysis of a stochastic logistic model with diffusion and Ornstein–Uhlenbeck process J. Math. Phys. 63 5 2022
34 Liu X. Ullah S. Alshehri A. Altanji M. Mathematical assessment of the dynamics of novel coronavirus infection with treatment: a fractional study Chaos Solitons Fractals 153 2021 111534
35 Mao X. Stochastic Differential Equations and Their Applications 1997 Horwood Publishing Chichester
36 Ndaïrou F. Area I. Nieto J.J. Torres D.F. Mathematical modeling of covid-19 transmission dynamics with a case study of Wuhan Chaos Solitons Fractals 135 2020 109846
37 Sher M. Shah K. Rassias J. On qualitative theory of fractional order delay evolution equation via the prior estimate method Math. Methods Appl. Sci. 43 10 2020 6464 6475
38 Siettos C.I. Russo L. Mathematical modeling of infectious disease dynamics Virulence 4 4 2013 295 306 23552814
39 Sk N. Mondal B. Thirthar A.A. Alqudah M.A. Abdeljawad T. Bistability and tristability in a deterministic prey–predator model: transitions and emergent patterns in its stochastic counterpart Chaos Solitons Fractals 176 2023 114073
40 Song Y. Zhang X. Stationary distribution and extinction of a stochastic SVEIS epidemic model incorporating Ornstein–Uhlenbeck process Appl. Math. Lett. 133 2022 108284
41 Su T. Zhang X. Stationary distribution and extinction of a stochastic generalized SEI epidemic model with Ornstein-Uhlenbeck process Appl. Math. Lett. 143 2023 108690
42 Thirthar A.A. Abboubakar H. Khan A. Abdeljawad T. Mathematical modeling of the covid-19 epidemic with fear impact AIMS Math. 8 3 2023 6447 6465
43 Thirthar A.A. Naji R.K. Bozkurt F. Yousef A. Modeling and analysis of an si1i2r epidemic model with nonlinear incidence and general recovery functions of i1 Chaos Solitons Fractals 145 2021 110746
44 Ullah S. Khan M.A. Modeling the impact of non-pharmaceutical interventions on the dynamics of novel coronavirus with optimal control analysis with a case study Chaos Solitons Fractals 139 2020 110075
45 Virgin H.W. IV Speckt S.H. Unraveling immunity to γ-herpesviruses: a new model for understanding the role of immunity in chronic virus infection Curr. Opin. Immunol. 11 4 1999 371 379 10448140
46 Wen B. Liu B. Cui Q. Analysis of a stochastic SIB cholera model with saturation recovery rate and Ornstein-Uhlenbeck process Math. Biosci. Eng. 20 7 2023 11644 11655 37501413
47 Zhai X. Li W. Wei F. Mao X. Dynamics of an hiv/aids transmission model with protection awareness and fluctuations Chaos Solitons Fractals 169 2023 113224
48 Zhang X. Yuan R. A stochastic chemostat model with mean-reverting Ornstein-Uhlenbeck process and Monod-Haldane response function Appl. Math. Comput. 394 2021 125833
49 Zhou B. Han B. Jiang D. Ergodic property, extinction and density function of a stochastic sir epidemic model with nonlinear incidence and general stochastic perturbations Chaos Solitons Fractals 152 2021 111338
50 Zhou B. Jiang D. Dai Y. Hayat T. Threshold dynamics and probability density function of a stochastic avian influenza epidemic model with nonlinear incidence rate and psychological effect J. Nonlinear Sci. 33 2 2023 29
51 Zhou B. Jiang D. Han B. Hayat T. Threshold dynamics and density function of a stochastic epidemic model with media coverage and mean-reverting Ornstein–Uhlenbeck process Math. Comput. Simul. 196 2022 15 44
52 Zhou B. Jiang D. Hayat T. Analysis of a stochastic population model with mean-reverting Ornstein–Uhlenbeck process and Allee effects Commun. Nonlinear Sci. Numer. Simul. 111 2022 106450
