Dynamical analysis of SIR model with Gamma distribution delay driven by Lévy noise and switching

1.   Introduction

As everyone knows, epidemics have a significant impact on socioeconomic conditions and human health, such as tuberculosis, measles, pertussis, and so on. Infectious diseases have great transmission power. If they are not effectively controlled, the epidemic will spread rapidly anywhere ^[1,2]. Mathematical models are typically employed in epidemiology to depict the transmission dynamics of diseases. Since Kermack and McKendrick introduced the SIR model to describe the spread of infectious diseases in ^[3], the SIR model considering the time delay has been widely adopted in the studies of epidemics ^[4,5,6,7]. Zaman et al. ^[8] developed the time-delayed SIR epidemic model incorporating control strategies to optimize the combined population of susceptible and recovered individuals. Kudryashov et al. ^[9] examined the SIR model for the propagation of coronavirus and investigated the analytical solutions of the model in conjunction with epidemiological data. Kumar et al. ^[10] analyzed the SIR model with time delay and determined the global stability of the disease-free equilibrium. To effectively regulate the dynamics of susceptible, infected, and recovered populations under time-varying vaccination control, the SEIR model with distributed delay was developed in ^[11]. This model utilizes a feedback mechanism to achieve this regulation. Krylova et al. ^[12] discussed the SIR model and found that when the average infection period is fixed, the transmission dynamics will change with the change of the distribution shape of the infection period. Many biological processes, especially disease transmission, involve the complexity of delays, so introducing the delays into epidemic models to characterize the delay of the infection state. For vector-borne diseases, the distributed time delay needs to be further considered.

Numerous research have found that environmental noise significantly influences the transmission of infectious diseases ^{[13,14,15,16]}. Ji et al. ^[17] presented the SIR model with stochastic disturbances and found that in the case of low noise intensity, there is a threshold that determines whether the epidemic will die out or persist. The SIS model of disease transmission coefficient and mortality subject to stochastic disturbance was analyzed in ^[18], which obtained sufficient conditions for thestochastic SIS model to have a unique ergodic stationary distribution. Du NH and Nhu NN ^[19] introduced the stochastic SIR model in which both natural mortality and morbidity were disturbed by white noise, which obtained sufficient and necessary conditions for the extinction and persistence of SIR. The global dynamic characteristics of the epidemic model and its corresponding stochastic differential equations were studied in ^[20], which found that the stochastic fluctuations could suppress the outbreak of disease. Li and co-authors ^[21] investigated the stochastic tumor immune model of integrated pulse therapy, which showed that the tumor status lasted from weak to disappeared as the noise intensity gradually increased. Dieu et al. ^[22] investigated the asymptotic behavior of a class of stochastic SIR models with degenerate diffusion, providing sufficient conditions approximating permanence.

The epidemic model may also be subjected to stochastic and some environmental disturbances, such as volcanic eruption, bird flu, tsunami, etc. These environmental factors will cause the mutation of the population size, which cannot be described by white noise. The introduction of Lévy noise into epidemic models by some researchers aims to more accurately reflect the impact of discontinuous, sudden environmental changes found in nature ^{[23,24,25,26]}. Zhou et al. ^[27] discussed the SIR model with Lévy noise disturbances and adopted an innovative approach to derive sufficient and almost necessary conditions for disease extinction and persistence. The stochastic SIR model with multidimensional Lévy noise is discussed in ^[28], which shows that persistence and disappearance depend not only on the variance of the process increment but also on the shape of its distribution. The dynamical behavior of the predation model driven by Lévy noise is discussed in ^[29], proving the existence and uniqueness of the understanding. Li et al. ^[30] analyzed the asymptotic behavior of the stochastic SIR model at its equilibrium point, comparing it to a deterministic SIR model.

In addition to Lévy noise, infectious diseases are also affected by color noise during the transmission process, causing the system to switch between different environmental states. Telegraph noise can be interpreted as a switch between two or more different states, such as temperature, seasonal turnover, or rainfall, among others ^{[31,32,33,34,35,36,37]}. For example, the stochastic SIR epidemic models with state switching are discussed, and Markov semigroup theory is used to prove the existence of stationary distributions in ^[31]. Continuous-time Markov chains are used to describe the state-switching behavior of the environment ^[33,34]. On this basis, Sereno et al. ^[37] proposed a switching nonlinear model predictive control strategy for SIR models to optimize non-drug interventions.

As mentioned above, numerous researchers have explored stochastic SIR epidemic models with time delays. Unfortunately, few researchers have considered introducing both Lévy noise and Markov switching into SIR epidemic models with distributed delays, which are more realistic than other models. However, this makes it challenging for us to theoretically analyze the asymptotic behavior of disease transmission. Hence, the primary contributions of this research are manifested in the following aspects:

$\bullet$ Based on ^[4,38], this paper describes the latency characteristics by introducing Gamma distributed delay, describes sudden environmental changes combined with Lévy noise, and adds Markov switching mechanism to reflect the dynamic changes of the environment. This multi-factor integrated modeling method improves the authenticity and reliability of the model.

$\bullet$ Inspired by Liu et al. ^[39], we construct a stochastic Lyapunov function tailored for regime switching, which allows us to prove the model possesses a unique stationary distribution. Applying this approach, we can analyze the long-term behavior of infectious diseases by uncovering the complex dynamics of disease persistence and the inherent weak stability of stochastic systems.

$\bullet$ Through theoretical analysis and numerical simulations, we explore the influence of environmental noise on disease transmission dynamics. Special attention is given to the probability density distribution of different compartments, analyzing their relationship with Markov switching and uncovering the internal mechanisms of disease spread under stochastic environments.

The paper proceeds as follows. Section 2 introduces the construction of the model. Section 3 establishes the existence of a unique global positive solution for any initial value. Section 4 derives the stationary distribution for the time-delay stochastic model. Section 5 presents numerical simulations verifying the theoretical results. The paper concludes with final remarks and future research directions.

2.   Model building

The SIR epidemic model with distributed delay, as presented in ^[4], is formulated as follows:

where $\Lambda$ is the birth rate, $\mu$ is the death rate, $\gamma$ represents the recovery rate, $\varepsilon$ represents the mortality rate caused by disease, $\beta$ is the rate of conversion from the susceptible to the infected state. In reality, the incubation period and infection period of the disease are usually not fixed values, but fluctuate within a certain range. The use of Gamma distribution delay can describe this variability well and can more accurately reflect the memory effect and time dependence of disease transmission. For instance, the quantity of infected persons is reliant upon not just the immediate transmission status but also the long-term impact of infection events. Thus, utilizing a Gamma distribution for the kernel delay $F(\cdot)$ offers a practical approach. That is $F(\tau) = \frac{{{\tau ^n}{\rho ^{n + 1}}{e^{ - \rho \tau }}}}{{n!}}, \tau > 0$. $\rho$ is the rate parameter of the gamma distribution, controlling the scale of the distribution and $\rho > 0$. In this study, we consider the weak kernel $n = 1$, we set $Z({\rm{t}}) = \int_{ - \infty }^t {\rho {e^{ - \rho (t - \tau)}}I(\tau)d\tau, }$ model (2.1) is converted as follows:

model (2.2) has the corresponding basic reproduction number and the invariant attracting set, which are given by ${R_0} = \frac{{\beta \Lambda}}{{\mu (\mu + \gamma + \varepsilon)}}$ and ${\wp _0} = \{ (S, I, R, Z)\left| {S \ge 0, } \right.I \ge 0, R \ge 0, Z \ge 0, S + I + R + Z \le \frac{\Lambda }{\mu }\}$. Moreover, two potential equilibrium states are described as follows:

$\bullet$ If ${R_0} \le 1$, the disease-free equilibrium ${E_0} = (\frac{\Lambda }{\mu }, 0, 0, 0)$, then ${E_0}$ is globally asymptotically stable in ${\wp _0}$, indicating that the disease will be eliminated from the population.

$\bullet$ If ${R_0} > 1$, the endemic equilibrium ${E^ + } = ({S^ + }, {I^ + }, {R^ + }, {Z^ + })$, where ${S^ + } = \frac{{\mu + \gamma + \varepsilon }}{\beta }$, ${I^ + } = \frac{{\beta \Lambda - \mu (\mu + \gamma + \varepsilon)}}{{(\mu + \gamma + \varepsilon)\beta }}$, ${R^ + } = \frac{{\gamma (\beta \Lambda - \mu (\mu + \gamma + \varepsilon))}}{{\beta \mu(\mu + \gamma + \varepsilon)}}$, ${Z^ + } = {I^ + }$, then ${E^ + }$ is globally asymptotically stable, but ${E_0}$ is unstable in the domain ${\wp _0}$. This suggests that the disease will continue to be prevalent for an extended period.

Inspired by the above discussion, this paper introduces the stochastic change factor to the environment, constructs the SIR epidemic system model with Lévy noise and Markov switching mechanism. The introduction of Lévy noise is used to capture the large stochastic disturbance that may occur in the real society; it is described by a Markov chain, the instantaneous dynamics of life rate switching between two or more different infectious environmental conditions, so as to build a more comprehensive and accurate disease transmission model. Based on the characteristics of Markov chains, for any finite state space $n \in N$, we construct the following switched hybrid model

where ${B_i}(t)(i = 1, 2, 3)$ are mutually independent standard Brownian motions, ${\sigma _{\rm{i}}}^2 > 0\; (i = 1, 2, 3)$ represents white noises. This paper proposes the following hypotheses:

${Assumption(A1)}$: ${\sigma _{\rm{i}}}^2(\cdot) > 0\; (i = 1, 2, 3)$ represent white noises. The process $\eta (t)$ is an irreducible and continuous-time Markov chain defined on the state space $N = \{ 1, 2, 3, ..., k\}$, $\eta (t)$ is assumed to be generated by the transition rate matrix ${({\mu _{ij}})_{k \times k}},$ so

where ${\mu _{ij}} > 0$ denotes the transition rate from state $i$ to $j$ and $i \ne j$, ${\mu _{ii}} = - \sum\nolimits_{i \ne j, i = 1}^k {{\mu _{ii}}}.$ Based on the irreducibility property of $\eta (t)$, there exists a unique stationary probability distribution $\phi = ({\phi _1}, {\phi _2}, ..., {\phi _k}) \in {\Re^{1 \times k}}$ subject to $\sum\nolimits_{n = 1}^k {{\phi _n} = 1}$ and ${\phi _n} > 0$ $(n \in N)$. It is also assumed that $\eta (t)$ and ${B_i}(t)(i = 1, 2, 3)$ are independent.

${Assumption(A2)}$: $r(u)$ is a bounded function $1 + r(u) > 0$ and $\left| {\frac{\Lambda }{\mu }r(u)} \right| \le \delta < 1, u \in Y.$

${Assumption(A3)}$: If the above assumptions are satisfied, there exists $\widetilde \chi (t, n) = {\left({S\left(t \right), I\left(t \right), R\left(t \right), Z\left(t \right), n} \right)^T} \in {R _ + }^4 \times {\rm{ }\it N}$ has $\mathop {\lim }\limits_{t \to \infty } \sup {\rm E}{\left| {\widetilde \chi (t, n)} \right|^n} \le Q(\varphi).$

3.   Existence and uniqueness of global positive solution

To investigate the dynamical behavior of the infectious disease model, the initial focus is on determining whether its solutions are globally positive. We will demonstrate that model (2.3) possesses a unique global positive solution for any initial condition, employing the technique of constructing a stochastic Lyapunov function.

Theorem 1. For any given initial value $\left({S\left(0 \right), I\left(0 \right), R\left(0 \right), Z\left(0 \right), n} \right) \in {R _ + }^4 \times {\rm{ }\it N}$, there exists a unique solution $(S, I, R, Z)$ for the model (2.3) on t $\geq$ 0, and the solution will remain in ${R _ + }^4$ with probability one, namely, the solution $\left({S\left(t \right), I\left(t \right), R\left(t \right), Z\left(t \right), n} \right) \in {R _ + }^4 \times {\rm{ }\it N}$ for all t $\geq$ 0 almost surely (a.s.).

Proof: Since the coefficients of model (2.3) satisfy the local Lipschitz condition, then for model (2.3) with initial value $\left({S\left(0\right), I\left(0 \right), R\left(0 \right), Z\left(0 \right), n} \right) \in {R _ + }^4 \times {\rm{ }\it N}$, there exists a unique local solution$\left({S\left(t \right), I\left(t \right), R\left(t \right), Z\left(t \right), n} \right)$ on $t \in \left[ {0, {\tau _e}} \right)$, where ${\tau _e}$ is an explosion time. To prove the local solution is global, we only need to prove ${\tau _e } = \infty$ a.s.. Assume that ${k_0}$ is so large that $\left({S\left(0 \right), I\left(0 \right), R\left(0 \right), Z\left(0 \right), n} \right) \in [\frac{1}{{{k_0}}}, {k_0}]$. For each integer $k > {k_0}$, define the stopping time as follows:

where $\inf\phi = \infty\; (\phi\; is\; the empty\; set)$. ${\tau _k}$ is increasing as $k \to \infty$. Assume that ${\tau _\infty } = \mathop {\lim }\limits_{k \to \infty } {\tau _k}$ and ${\tau _\infty } \le {\tau _e}$ a.s.. We need to prove that ${\tau _\infty } = \infty$ a.s., then ${\tau _e} = \infty$ a.s.. Otherwise, there exist the constants $T > 0$ and $\delta \in \left({0, 1} \right)$ such that $P({\tau _\infty } \le T) \geq \delta$. Hence there exist the integers ${k_1} \geq {k_0}$, which take

Define a fundamental function V: ${R_ + }^4 \to {R_ + }$

We have the following inequality to get $V\left({S, I, R, Z} \right)$ is non-negative

By applying $\mathrm{It\hat o's}$ formula

where

choosing $a = \frac{\beta }{\rho }, b = \frac{{\mu + \gamma + \varepsilon }}{\rho }$. which ${H_1} = - \int_Y [ \ln (1 - {r_1}(n, u)) - {r_1}(n, u)]vdu, {H_2} = - \int_Y [ \ln (1 - {r_2}(n, u)) - {r_2}(n, u)]vdu$, ${H_3} = - \int_Y [ \ln (1 - {r_3}(n, u)) - {r_3}(n, u)]vdu$. From Taylor's formula and hypothesis A2, we can get:

where $\theta \in (0, 1)$ is an arbitrary constant, ${H_2}$ and ${H_3}$, so on.

where $\widetilde K\ > 0$, therefore we can obtain

Integrate both sides of Eq (3.2) from 0 to ${\tau _k} \wedge T$ and draw the expectation then

Let ${\Omega _k} = \left\{ {{\tau _k} \le T} \right\}$, for $w \in {\Omega _\xi }$ and in view of Eq (3.2), we obtain $P({\Omega _k}) \ge \varepsilon$ such that for each $w \in {\Omega _\xi }$, $S({\tau _k}, \omega)$ or $I({\tau _k}, \omega)$ or $R({\tau _k}, \omega)$ or $Z({\tau _k}, \omega)$ equals either $k$ or $\frac{1}{k}$, so

From Eq (3.2), we can obtain

where ${1_{{\Omega _k}}}$ is the indicative function of ${\Omega _k}$, let $k \to \infty$, then

Thus ${\tau _\infty } = \infty$ holds almost everywhere. Consequently, Theorem 1 is verified.

Remark 1. In a biological sense, the population of each compartment cannot be negative. For the established model (2.3) to have practical biological significance, it is necessary to prove the existence and uniqueness of a global positive solution.

4.   The existence of ergodic stationary distribution

Next, we prove the existence of an ergodic stationary distribution of model (2.3). To facilitate this, we begin by presenting the following symbols and lemma.

Lemma 4.1. ^[38,40]For a bounded open domain $\bigcup \in R^4$ with regular boundary $\partial \bigcup$, the following statement is valid

$(i)$ There exists a positive constant $K$ such that for all $x \in \cup, \xi \in {R^4}$, $\sum\limits_{i, j = 1}^4 {{a_{ij}}(x){\xi _i}} {\xi _j} \ge K{\left| \xi \right|^2}.$

$(ii)$ There is a non-negative ${C^2}$ function ${u}(t)$ such that $L{u}(t) < 0$ for any ${R^4}\backslash \cup.$

Then Markov process $\chi(t)$ is characterized by a unique ergodic stationary distribution, which is represented by $\prod(\cdot)$.

Theorem 2. For any initial value $\left({S\left(0 \right), I\left(0 \right), R\left(0 \right), Z\left(0 \right), n} \right) \in {R _ + }^4 \times {\rm{ }\it N}$, the presence of ${R_s} > 1$ guarantees that a unique ergodic stationary distribution is present in $\in {R _ + }^4 \times {\rm{ }\it N}$.

Proof: The model (2.3) has the following diffusion matrix

From this we can get that condition $(i)$ of Lemma 4.1 is fulfilled. To confirm Lemma 4.1 criterion $(ii)$, we introduce the function ${u_1}(t)$ as follows:

where ${m_1} = 1, {m_2} = \frac{{\rho(n) \beta(n) \Lambda(n) }}{{\rho(n) [\mu(n) + \gamma(n) + \varepsilon(n) + {\varsigma _2} + \frac{1}{2}{{\sum\nolimits_{n = 1}^k {{\phi _n}{\sigma _2}(n)]} }^2}}}, {m_3} = \frac{{\rho(n) \beta(n) \Lambda(n) }}{{{\rho(n) ^2}[\mu(n) + {\varsigma _1} + \frac{1}{2}{{\sum\nolimits_{n = 1}^k {{\phi _n}{\sigma _2}(n)]} }^2}}}.$ The function $\overrightarrow \omega (n)$ represents a Markov process and is twice continuously differentiable. Its $\mathrm{It\hat o's}$ derivative takes the following form:

We obtain

Utilizing the irreducibility characteristic of the Markov process, for ${\sigma _1}^2(n) = ({\sigma _1}^2(1), {\sigma _1}^2(2), ...., {\sigma _1}^2(K)),$ there exist $\overrightarrow \omega (n) = (\overrightarrow \omega (1), \overrightarrow \omega (2), ...\overrightarrow \omega (K)),$ and $\overrightarrow \omega (n)$ is defined in Eq (4.2) such that

For ${\sigma _2}^2(n) = ({\sigma _2}^2(1), {\sigma _2}^2(2), ...., {\sigma _2}^2(K)),$ we have

where $n \in N = \left\{ {1, 2, 3, ..., K} \right\}$ and ${u_{nj}} > 0$, then it yields

where ${R_s} = \frac{{\beta(n) \Lambda(n) }}{{[\mu(n) + {\varsigma _1} + \frac{1}{2}\sum\nolimits_{n = 1}^k {{\phi _n}{\sigma _1}^2(n)][\mu(n) + \gamma(n) + \varepsilon(n) + {\varsigma _2} + \frac{1}{2}\sum\nolimits_{n = 1}^k {{\phi _n}{\sigma _2}^2(n)]} } }}$. Then, we define ${u_2}(t)$,

By $\mathrm{It\hat o's}$ formula, we obtain

In addition, we have

And then

Therefore, it can be inferred that

where ${r_{\min }}(u) = \{ {r_1}(n, u) \wedge {r_2}(n, u) \wedge {r_3}(n, u)\}$ and ${r_{\max }}(u) = \{ {r_1}(n, u) \vee {r_2}(n, u) \vee {r_3}(n, u)\}.$ Moreover, it can be deduced that

where ${\sigma _{\max }}^2(n) = \{ {\sigma _1}^2(n) \vee {\sigma _2}^2(n) \vee {\sigma _3}^2(n)\}$ and ${w_1} = \frac{1}{{2(S(t) + I(t) + R(t))}}\{ \mu(n) + \varepsilon(n) + {\sigma _{\max }}^2(n) + \int_Y [ \frac{1}{{1 + {r_{\min }}(n, u)}} - 1 + {r_{\max }}(n, u)]vdu\} + \Lambda(n).$ Then, we define ${u_3}(t)$

Using $\mathrm{It\hat o's}$ formula, we have

where ${w_2} = \frac{1}{{R(t)}}\int_Y [ \frac{1}{{1 + {r_3}(n, u)}} - 1 + {r_3}(n, u)]vdu$. Finally, $u(t)$ is defined

The function $u(t)$ is continuous. When $u(t) \to \infty$, there is a lower value ${u_*}({S_*}, {I_*}, {R_*}, {Z_*})$ at point of $({S_*}, {I_*}, {R_*}, {Z_*})$ formulate a nonnegative function.

if assumption (A3) holds, it can be obtained

When $S(t) \to {0^ + }$ or $I(t) \to {0^ + }$ or $R(t) \to {0^ + }$ or $Z(t) \to {0^ + }$, it gives

When $S(t) \to + \infty$ or $I(t) \to + \infty$ or $R(t) \to + \infty$ or $Z(t) \to + \infty$, it gives

Based on the analysis in (4.3) and (4.4), if ${R_s} > 1$, then there exist sufficiently small have $L\tilde u(t) < - 1$. This satisfies the conditions stated in Lemma 4.1, thereby concluding the proof of the theorem.

Remark 2. When ${R_s} > 1$, the model (2.3) has the ergodic stationary distribution, indicating that the epidemic will persist in the population. Additionally, by examining ${R_s}$ and ${R_0}$, it can be obtained that ${R_s}$ = ${R_0}$ when ${\sigma _1} = {\sigma _2} = {\sigma _3} = 0.$ Thus, this reveals the critical impact of stochastic fluctuations on the spread of epidemics.

5.   Numerical simulations

In this paper, we discuss the long-term behavior of a class for SIR models with distributed delays driven by Lévy noise and switching. Next, we will employ numerical simulations to validate the aforementioned theoretical findings. Define the finite state space of continuous-time Markov chain $\eta (t)$ as $N = \left\{ {1, 2} \right\}$, with its generator matrix $\Gamma$ is

By solving $\pi \Gamma = 0,$ we find that the unique stationary distribution of the aforementioned Markov chain $\eta (t)$ is $({\pi _1}, {\pi _2}) = (0.8, 0.2)$.

The parameter values presented in Table 1 have been obtained.

5.1.   Theoretical verification

To be reasonable, around each parameter value in Table 1, we chose the following parameter values for numerical simulations. The model with initial value $(S(0), I(0), R(0)) = (50, 10, 0.1)$.

By direct calculation, we have

The conditions of Theorem 2 are satisfied. Theorem 2 implies the trend observed in model (2.3), a unique stationary distribution evident in Figure 2. By looking closely at the curves representing S(t), I(t), and R(t), it is clear that the disease is still present. From a biological perspective, this indicates that the disease will persist in its present state.

Figure 3 exhibits the stationary probability density (SPD) of S(t), I(t), and R(t), using the same parametric values as those in Figure 2. Figure 3 shows that the probability density functions of S(t), I(t), and R(t) are confined to a narrow interval, each exhibiting two local peaks that correspond to the two states 1, 2 of Markov chain $\eta (t)$, where $N = \left\{ {1, 2} \right\}$. States 1 and 2 delineate two different epochs or conditions within the disease transmission process. Specifically, the concentration of the probability density function within the cell indicates that the values of S(t), I(t), and R(t) tend to cluster around specific points, reflecting the stability of the system in these states. The presence of two local peaks further indicates two states corresponding to Markov chains.

In order to further study the high transmission state and low transmission state in the process of disease transmission, we drew the time trajectory diagram of each state, as shown in Figure 4. The figure illustrates state 1 with a black line and state 2 with a green line. By observing, we can clearly see that the number of susceptible persons in state 1 is located below state 2. Based on this observation, it can be concluded that state 1 represents the stage of low or limited transmission of the disease and state 2 represents the stage of high or widespread transmission of the disease. This conclusion is also consistent with the value of the state parameter we choose.

5.2.   Influence of noise on disease

Figure 5 illustrates how a switching mechanism influences disease dynamics when subjected to Lévy noise. The selection of parameters is consistent with Figure 2; these parameters satisfy the theorem's conditions for traversal stationary distributions. By looking at Figure 5, it is clear that diseases with Lévy noise exhibit stronger volatility under the action of the switching mechanism. Using the parameters from Table 3, we examine how variations in noise intensity affect disease progression. The other parameters are the same as those in Figure 5.

As shown in the comparative analysis of Figures 2 and 6, with the increase of noise intensity, the sizes of the susceptible population and the infected population gradually tend to a stable state and show a significant attenuation trend. This dynamical behavior indicates that the enhancement of stochastic disturbances may drive the disease transmission system to shift towards an extinction phase. Specifically, when the intensity of environmental noise exceeds a certain critical threshold, the stochastic stability of the system is disrupted, resulting in the inability of the number of infected people to maintain sustainable transmission, and thereby promoting the disappearance of the disease within a local area.

6.   Conclusions

This paper primarily investigates the dynamic behavior of the SIR epidemic model with gamma-distributed delays under the influence of Lévy noise and Markovian switching. Initially, we establish the existence and uniqueness of the global positive solution for the model and subsequently explore the stationary distribution of model (2.3). Under ${R_s} > 1$, the existence of the ergodic stationary distribution for model (2.3) is derived. This indicates that the infectious disease will persist within the population for an extended period. Furthermore, from the ${R_s}$ and ${R_0}$, it can be obtained that ${R_s}$ = ${R_0}$, when ${\sigma _1} = {\sigma _2} = {\sigma _3} = 0.$

The research methodology employed in this study holds significant potential for broader application in future studies of stochastic systems. For example, it can be effectively applied to models featuring multiple delays and fractional order noise. In addition, research can be extended to other related fields, such as ecosystems, complex networks, public opinion dissemination, etc., where stochastic processes and their interactions play a crucial role. The insights derived from our research results can effectively advance the surveillance, prevention, and control of infectious diseases.

Use of AI tools declaration

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

Authors' contribution

Jing Yang: Writing - editing, data processing, software. Shaojuan Ma: Funding, writing - review, methodology. Dongmei Wei: Formal analysis, Review.

Acknowledgments

This work was supported by the grants from the National Natural Science Foundation of China (No.12362005), Key Project of Natural Science Foundation of Ningxia (No.2024AAC02033), Major Special project of North Minzu University (No.ZDZX201902) and Postgraduate innovation project of North Minzu University (YCX24073).

Conflict of interest

The authors declare there is no conflict of interest.