Threshold dynamics of a stochastic SIRS epidemic model with transfer from infected individuals to susceptible individuals and log-normal Ornstein-Uhlenbeck process

1.   Introduction

Infectious disease has been a natural enemy of human health. There is an urgent need to explore the transmission mechanism of the disease. Mathematical model plays a central role in predicting the spread trend of infectious disease. Kermack and Mckendrick ^[1] proposed a landmark susceptible-infected-recovered (SIR) epidemic model and assumed that the total population $N(t)$ is divided into three categories at time $t$, including susceptible individuals $S(t)$, infected individuals $I(t)$, and recovered individuals $R(t)$. Recovered individuals in the SIR epidemic model are permanent immunity. In fact, instead of lifelong immunity, recovered individuals may have acquired immunity loss. For this phenomenon, numerous researchers are dedicated to developing the SIRS epidemic model, and various improved versions of the SIRS model have been proposed. For instance, Liu et al. ^[2] investigated the influence of nonlinear incidence rates on the dynamics of SIRS epidemiological models. They presented a modified general criterion for supercritical or subcritical Hopf bifurcation of two-dimensional systems. For a multi-group SIRS epidemic model with varying population sizes, Muroya et al. ^[3] established sufficient conditions for the global stability of an endemic equilibrium. Sekiguchi ^[4] considered a discrete SIRS epidemic model with time delays and obtained the sufficient condition for the permanence by the Euler method.

It is possible for the infected individuals to recover after undergoing certain treatments and go back directly to the susceptible individuals due to the presence of transient antibody. To describe this phenomenon, Bai and Mu ^[5] considered the transfer rate from infected individuals to susceptible individuals in a generalized SIRS epidemic model. Avila-Vales and Pérez ^[6] proposed a diffusive SIRS epidemic model with the same transfer rate and analyzed the dynamical characteristics of the system. Chen ^[7] constructed the following SIRS epidemic model with partial immunity, transfer from $I(t)$ to $S(t)$, and standard incidence

where $N(t) = S(t)+I(t)+R(t)$. $A$ is the natural birth rate. $d$ represents the natural mortality rate. $a$ is the disease-related death rate. $\bar{\beta}$ denotes the average incidence rate. $\gamma$ stands for the immunity loss rate of recovered individuals. $c$ and $\eta$ represent the transfer rate from infected individuals to recovered individuals and susceptible individuals, respectively. The basic reproduction number of system (1.1) is $R_0 = \frac{\bar{\beta}}{a+c+\eta+d}$. The authors in ^[7] showed the dynamical properties of system (1.1) as follows.

Theorem 1.1. If $R_0 < 1$, the disease-free equilibrium $E_0 = (\frac{A}{d}, 0, 0)$ is locally asymptotically stable. If $R_0 > 1$, the system has a positive equilibrium $E^\ast(S^*, I^*, R^*)$, where $S^* = \frac{(a+c+\eta+d)N^*}{\bar{\beta}}$, $I^* = \frac{A-dN^*}{a}$, $R^* = \frac{cI^*}{d+\gamma}$, and $N^*$ is the positive root of the following equation:

It is widely recognized that the real world is full of randomness. Deterministic models often use fixed parameter values to describe the transmission mechanism of the disease. Although this approach provides a clear framework, it does not fully take into account the random factors in the environment. As May ^[8] noted, the parameter in the model is inevitably affected by all kinds of noise in the surrounding environment, which can come from temperature, humidity, air pressure, and many other factors. Stochastic models can capture the inherent uncertainty in the disease transmission. In recent years, many researchers have introduced random perturbations into deterministic models. For instance, Privault and Wang ^[9] considered a general stochastic SIR epidemic model driven by a three-dimensional Lévy noise. They showed that the variance of the processes increments and the shapes of their distributions can have a significant influence on the dynamical behavior of the system. For a SIRS epidemic model with white noise and standard incidence rate, Zhang et al. ^[10] derived the condition for the extinction of the disease and proved that there exists stationary distribution. By replacing $\beta dt$ by $\beta dt+\sigma dB(t)$ in the classical deterministic SIRS model, Settati et al. ^[11] proposed a stochastic SIRS epidemic model and obtained that the dynamics is governed by a threshold, which characterizes the extinction and persistence of the disease. Besides, their results showed that when the threshold is equal to one, small noise will make the disease extinct and large noise promotes persistence.

In the construction of stochastic models, parameter perturbation technology is a widely recognized tool. This technique allows researchers to introduce stochasticity into complex dynamic models and simulate random changes in the real-world environment, so that the model can better predict and explain the actual observed parameter variation. According to the existing research, there are two common approaches to incorporate environmental variability by modifying the parameters. One way is to assume that parameter is a linear function of Gaussian white noise ^{[11,12,13,14]}. In infectious disease models, the incidence rate refers to the frequency at which individuals in a population come into contact with one another, potentially leading to the transmission of pathogens. This parameter is particularly susceptible to random fluctuations. For example, individuals' social behaviors can vary widely due to factors such as mood, seasonality, or public health interventions, like lockdowns. These variations are often unpredictable and can lead to increases or decreases in incidence rate. Assume that $\bar{\beta}$ is a linear function of Gaussian white noise, which is expressed as

where $B(t)$ is the standard Brownian motion and is defined on a complete probability space $\{\varOmega, \mathscr{F}, \{\mathscr{F}_t\}_{t\geq0}, \mathbb{P}\}$ with a filtration $\{\mathscr{F}_t\}_{t\geq0}$ satisfying the usual conditions (i.e., it is right continuous and $\mathscr{F}_0$ contains all $\mathbb{P}$-null sets), and $\varrho^2$ is the noise intensity of $B(t)$. Let $\langle \beta(t) \rangle$ be the time average of $\beta(t)$ for any time interval $[0, t]$. By directly integrating (1.2), one can get

where $\mathbb{N}(\cdot, \cdot)$ denotes the one-dimensional Gaussian distribution. The variance of $\langle \beta(t) \rangle$ goes to infinity when $t\rightarrow 0$.

The other way is to suppose that biological parameter satisfies the Ornstein-Uhlenbeck (OU) process. The OU process, also known as the mean-reverting process, is a stochastic process that tends to move back toward a long-term equilibrium or mean value over time ^[15]. It is characterized by two main parameters: the speed of reversion, which determines how quickly the process returns to the mean, and the level of mean reversion, which represents the long-term equilibrium value toward which the process tends to revert. The OU process has applications in finance, economics, and physics, where it is used to capture the inherent tendency of certain systems to return to a stable state over time ^[16]. In recent years, with the continuous progress of statistics, increasing attention has been paid to the dynamics of biological models with the OU process ^{[17,18,19,20,21]}. This way avoids the problem of the first approach, but there is a problem that the parameter satisfying this way may take negative values. To guarantee the positive of the biological parameter, Allen ^[15] introduced environmental disturbance to deterministic system by assuming that the parameter satisfies the log-normal OU process.

Applying the log-normal OU process to infectious disease models has significant biological implications: (ⅰ) Understanding stochastic dynamics: Infectious diseases do not spread in a perfectly predictable manner. Instead, they are influenced by various unpredictable factors like environmental changes, human behavior, and random events. The log-normal OU process incorporates randomness into the model, which allows us to account for these uncertainties and better reflect real-world dynamics. (ⅱ) Mean-reverting behavior: The log-normal OU process has a natural tendency to return to an average state over time. This means that if the number of infections spikes or drops significantly, it will likely stabilize back toward an expected level. It is important for understanding diseases that become endemic constantly present in populations, as they tend to fluctuate around a certain level rather than going extinct or growing uncontrollably. (ⅲ) Logarithmic transformation: Infectious diseases often grow exponentially at first but may slow down as the population becomes immune or interventions are applied. Logarithmic transformation helps model this growth more accurately by stabilizing variations. In addition, logarithmic transformation can ensure the positivity of parameters subject to random perturbation. Based on the above background and biological significance, we assume that $\bar{\beta}$ is affected by the following log-normal OU process. It can be expressed as

where $K$ and $\sigma$ are positive constants. $K$ is the speed of reversion. $\sigma$ is the intensity of volatility. Let $z(t) = \log \beta(t)$ and $\bar{z} = \log \bar{\beta}$. Equation (1.3) becomes

By integrating (1.4), the explicit solution with initial value $z(0)$ is

For convenience, letting $z(0) = \bar{z}$, we calculate that the average value of $z(t)$ is

One can obtain

where $o (t^2)$ is the high-order infinitesimal of $t^2$. Since $z(0) = \bar{z} = \log \bar{\beta}$, we have

By virtue of the result in ^[15], for large $t$, the probability density of $\beta(t)$ approaches a stationary log-normal density with mean $\bar{\beta} e^{\frac{\sigma^2}{4K}}$ and variance $\bar{\beta}^2(e^{\frac{\sigma^2}{K}(1-e^{-2Kt})}-e^{\frac{\sigma^2}{2K}(1-e^{-2Kt})})$. One can see that with sufficient shortening of the time interval, it becomes apparent that the variance of the variability level of the incidence rate gradually approaches zero. This suggests that the incidence rate has remained stable and consistent over time. Therefore, this modeling method is more reasonable.

According to the above analysis, we formulate the stochastic version corresponding to deterministic system (1.1) as follows:

This paper makes the following main contributions and innovations: (ⅰ) To introduce environmental perturbations and keep the positive of the parameter, we consider that the incidence rate satisfies the log-normal OU process. This process enables more precise forecasting of the spread of the disease. (ⅱ) Mathematical analysis is employed to determine a critical value $R_0^*$, which serves as a threshold determining the extinction and persistence of the disease. In precise terms, if $R_0^* < 1$, then the disease will be extinct with probability one; if $R_0^* > 1$, then the disease persists by demonstrating the existence of stationary distribution. (ⅲ) By solving the four-dimensional algebraic equations, we derive the exact form of probability density function of the stationary distribution. (ⅵ) Numerical simulations are conducted. We use the variable-controlling approach to analyze the influence of some key parameters on the dynamics of the considered model.

The paper is organized as follows. Section 2 proves the existence of stationary distribution. In Section 3, we derive the exact form of covariance matrix in the probability density function for distribution. Sufficient condition for the extinction of infectious disease is obtained in Section 4. In Section 5, using some examples and numerical simulations, we confirm the theoretical results. Finally, we make a discussion in Section 6.

For convenience, there are some mathematical notations. $\mathbb{R}_{+}^n = \{(x_1\cdots x_n)\in\mathbb{R}^n|x_i > 0, 1\leq i\leq n\}$. If $P$ is a matrix, $P^T$ and $P^{-1}$ represent its transpose and inverse matrix, respectively. $P^{(n)}$ represents a matrix consisting of the first $n$ rows and first $n$ columns of matrix $P$. To proceed, a theorem that the stochastic system (1.7) has a unique positive global solution is introduced.

Theorem 1.2. For any initial value $(S(0), I(0), R(0), \beta(0))\in \mathbb{R}_{+}^4$, system (1.7) has a unique positive solution $(S(t), I(t), R(t), \beta(t))$ and the solution will remain in $\mathbb{R}_{+}^4$ almost surely for all $t > 0$.

Remark 1.3. It is easy to get the following inequation concerning $N$ from system (1.7)

which implies that there is a positive invariant set $\Gamma$ of model (1.7)

In the following sections, we always assume the initial value $\frac{A}{a+d} < S(0)+I(0)+R(0) < \frac{A}{d}$.

2.   Stationary distribution

Stationary distribution provides insights into the equilibrium states that a disease can reach, helping researchers predict the prevalence of infection over time under various conditions. This information is vital for public health planning and resource allocation, as it allows policymakers to assess the potential burden of a disease and implement appropriate control measures. In this section, we construct sufficient condition for the existence of stationary distribution, which means the persistence of the disease. Define

Consider a $k$-dimensional nonlinear stochastic differential equation with initial value $\Upsilon(0)$

where $\check{B}(t)$ is a $\ell$-dimensional Brownian motion. $g_1:\mathbb{R}^k\rightarrow \mathbb{R}^k$ and $g_2:\mathbb{R}^k\rightarrow \mathbb{R}^{k\times \ell}$ are Borel measurable.

Lemma 2.1. (^[22,23,24]). Assume that there is a bounded closed domain $\Xi\subset{\mathbb{R}}^k$ with a regular boundary $\mathcal{M}$. For any initial value $\Upsilon(0)\in\mathbb{R}^k$, if

where $\mathbb{P}(\tau, \Upsilon(0), \cdot)$ is the transition probability of $\Upsilon(\tau)$, then there exists a solution of system (2.1) and admits at least one stationary distribution on $\mathbb{R}^k$.

Theorem 2.2. Assume that $R_0^* > 1$, then system (1.7) admits at least one stationary distribution on set $\Gamma$.

Proof. Step 1. (Construction of stochastic Lyapunov functions): Construct a $C^2$-Lyapunov function $\bar{V}(S, I, R, \beta):\Gamma\rightarrow \mathbb{R}$ as follows:

where $V_1 = -\log I+\dfrac{(a+d)}{4\varepsilon A(d+\gamma)}R$, $V_2 = -\log S-\log R-\log\left(\dfrac{A}{d}-(S+I+R)\right)-\log\left((S+I+R)-\frac{A}{a+d}\right)+\beta-1-\log \beta$, $\varepsilon = \frac{(a+c+\eta+d)(R_{0}^{*}-1)}{2\hat{\beta}}$, $\hat{\beta} = \bar{\beta}^{2}e^{\frac{\sigma^2}{K}}$. $M_0$ is a sufficiently large positive number satisfying the following inequality:

in which

Letting the differential operator $\mathcal{L}$ act on $-\log I$ yields

Using the Young inequality $z\leq \xi z^{2}+\frac{1}{4\xi}(\xi > 0)$ and the lower bound of $N$, we obtain

Set $F(\beta) = (\tilde{\beta}-\beta)+\varepsilon(\beta^2-\hat{\beta}^2)$. It follows from (2.2) and (2.3) that

Therefore,

Making use of $It \widehat{o}'s$ formula to $V_2$ results in

Combining (2.4) and (2.5), we obtain

where

Since $\bar{V}$ tends to $+\infty$ as $(S, I, R, \beta)$ approaches the boundary of $\Gamma$, $\bar{V}$ has a minimum value $\bar{V}^*(\bar{S}^*, \bar{I}^*, \bar{R}^*, \bar{\beta}^*)$. Therefore, a nonnegative function $V(S, I, R, \beta)$ is defined by $V(S, I, R, \beta) = \bar{V}(S, I, R, \beta)-\bar{V}^*(\bar{S}^*, \bar{I}^*, \bar{R}^*, \bar{\beta}^*)$. Applying $It \widehat{o}'s$ formula to $V$ has

Step 2. (Construction of a compact set): Construct a bounded set $\mathbf{U}$

where $\rho$ is a sufficiently small positive number satisfying

We divide the set $\Gamma \setminus \mathbf{U}$ into seven subsets $\mathbf{U}^c_{j}$, $j = 1, 2, \ldots, 7$:

We can prove $H(S, I, R, \beta)\leq-1$ in the above seven subsets. In other words,

Analysis of the expression of $H(S, I, R, \beta)$ leads us to the conclusion that there exists a positive constant $L$ such that

Step 3. (Existence of stationary distribution): To simplify the notation, we define $\textbf{Z}(t) = (S(t), I(t), R(t), \beta(t))$ as the solution of system (1.7) with initial value $\textbf{Z}(0) = (S(0), I(0), R(0), \beta(0))$. Integrating both sides of (2.6) from $0$ to $t$, dividing by $t$, and taking mathematical expectation yields

From (2.7) and (2.8), it follows that

where $\mathbf{1}$ is an indicator function. This, together with (2.9) yields

In view of the ergodicity theorem ^[25], one can obtain that

This implies immediately that

Taking the inferior limit on both sides of (2.10) and combining with (2.12), we get

According to Lemma 2.1, one can obtain that system (1.7) has at least one stationary distribution on $\Gamma$. This completes the proof. □

3.   Probability density function

Studying the probability density function of the distribution for the infectious disease system is crucial as it provides a detailed and continuous description of the long-term behavior. We can gain insights into the likelihood of various states within the system and understand how often certain conditions are expected to occur by analyzing the probability density function. The main purpose of this section is to calculate the concrete form of the covariance matrix and prove that it is positive definite.

Lemma 3.1. ^[21] For the real algebraic equation

where $\Omega = diag(1, 0, 0, 0)$, $\Delta$ is a real symmetric matrix, and

If $\theta_1 > 0$, $\theta_3 > 0$, $\theta_4 > 0$ and $\theta_1(\theta_2\theta_3-\theta_1\theta_4)-\theta_3^2 > 0$, then $\Delta$ is positive definite.

The dynamical property of the corresponding deterministic model of (1.7) is the same as the following system about $N$.

As a matter of convenience, we study the density function of this system. If $R_0 > 1$, system (3.1) exists a positive equilibrium $\check{E}$ = ($N^*$, $I^*$, $R^*$, $\log \bar{\beta}$). Let $Z = (z_1, z_2, z_3, z_4)^T = (N-N^*, I-I^*, R-R^*, \log \beta-\log \bar{\beta})^T$. System (3.1) can be linearized at $\check{E}$

System (3.2) can be equivalent to writing in the following form:

in which

and $a_{11} = d$, $a_{12} = a$, $a_{21} = \frac{\bar{\beta}I^*(I^*+R^*)}{(N^*)^2}$, $a_{22} = (a+c+\eta+d)-\frac{\bar{\beta} (N^*-2I^*-R^*)}{N^*}$, $a_{23} = \frac{\bar{\beta}I^*}{N^*}$, $a_{24} = \frac{\bar{\beta}I^*(N^*-I^*-R^*)}{N^*}$, $a_{32} = c$, $a_{33} = d+\gamma$, $a_{44} = K$.

The characteristic polynomial of matrix $A^{(3)}$ takes the form

where

and

In view of the Routh-Hurwitz criterion ^[26], $A^{(3)}$ has all negative real-part eigenvalues.

Theorem 3.2. If $R^*_0 > 1$, then the stationary distribution of system (3.1) around the equilibrium point $\check{E}$ has a normal density function as follows:

The concrete form of covariance matrix $\varSigma$ is presented in the following proof.

Proof. Similar to the proof in paper ^[21], our goal is to solve the following equation.

There is an ordering matrix

Calculation yields

Let the elimination matrix $J_2$ be

It is easy to get

Define

where $\rho = a_{12}a_{24}(a_{11}-a_{33})$, $\varpi = a_{12}a_{33}(a_{22}+a_{33})-a_{11}a_{12}(a_{11}+a_{22})$, $\upsilon = a_{12}a_{23}(a_{33}-a_{11})+\frac{a_{21}a_{12}^2(a_{33}-a_{11})+a_{12}(a_{11}^3-a_{33}^3)}{a_{32}}$, $\omega = a_{12}a_{21}(a_{11}-a_{33})-a_{11}^3$. Then, we get

where $\xi_1 = a_{11}+a_{22}+a_{33}+a_{44}$, $\xi_2 = a_{11}(a_{12}+a_{33}+a_{44})+a_{22}(a_{33}+a_{44})+a_{12}a_{21}+a_{23}a_{32}+a_{33}a_{44}$, $\xi_3 = a_{11}(a_{22}a_{33}+a_{22}a_{44}+a_{33}a_{44})+a_{22}a_{33}a_{44}+a_{12}a_{21}(a_{33}+a_{44})+a_{23}a_{32}(a_{11}+a_{44})$ and $\xi_4 = a_{44}(a_{11}a_{22}a_{33}+a_{12}a_{21}a_{33}+a_{11}a_{23}a_{32})$. Through above similar transformation, the equivalent equation of (3.4) can be written as

The characteristic polynomial of matrix $A_{3}$ takes the form

Obviously, $A_{3}$ has all negative real-part eigenvalues. In view of the Routh-Hurwitz criterion ^[26], we have

From Lemma 3.1, we derive that the matrix $\bar{\varSigma}: = \frac{1}{(\rho\sigma)^2}(J_3J_2J_1)\varSigma(J_3J_2J_1)^T$ is positive definite and by solving (3.5), the exact form of $\bar{\varSigma}$ is

Therefore, it is easy to derive the exact expression of $\varSigma$. That is,

which is also a positive definite matrix. □

4.   Extinction

Examining extinction offers valuable knowledge about the likelihood of a disease disappearing from a population. Understanding the factors that influence extinction can guide effective intervention strategies aimed at eradicating or controlling outbreaks. In this part, we will give sufficient condition for the extinction of infectious disease.

Theorem 4.1. If $R_0^* < 1,$ then

which means the disease of system (1.7) will die out almost surely.

Proof. Considering the second equation of system (1.7) yields

Integrating (4.1) from 0 to $t$ and dividing by $t$ on both sides, we have

Taking the superior limit on both sides of (4.2) and combining (2.11) leads to

It implies $\lim\limits_{t \to \infty}I(t) = 0$ a.s. That is to say, the disease will be extinct in a long term. □

Remark 4.2. Theorems 2.1 and 4.1 provide sufficient conditions for persistence and extinction of the disease, indicating threshold value $R_{0}^{*}$ plays a critical role in determining the dynamics of disease transmission. If $R_{0}^{*} < 1$, each infected individual leads to less than one new infection on average, resulting in a decline in the number of cases over time. This condition is crucial for disease extinction. If $R_{0}^{*} > 1$, this indicates that each infected individual produces more than one new infection on average, allowing the disease to persist in the population. The balance between new infections and recoveries determines whether the disease remains endemic.

Remark 4.3. We can observe that $R_{0}^{*}$ coincides with $R_{0}$ of deterministic system (1.1) when there is no noise. This implies we generalize the results of the deterministic system. In addition, it can be easily observed that $R_{0}^{*} > R_{0}$ when there exists noise. There is a situation where $R_{0}^{*} > 1$, when $R_{0} < 1$, that is, in the case of disease extinction predicted by the deterministic system, the disease in the stochastic system is likely to be persistent. Therefore, the influence of noise on system dynamics cannot be ignored.

5.   Numerical simulations

This section presents comprehensive numerical simulations to verify the theoretical results proved in this paper. The numerical analysis employs higher-order numerical methods of Milstein ^[27], which offers a power tool for capturing the complex dynamics of disease transmission. The findings will assist policymakers in comprehending and predicting the future trend of the epidemic, thereby facilitating the development of more efficacious prevention strategies. The corresponding discretized equation of (1.7) is acquired as follows

where $z = \log \beta$, $(S^i, I^i, R^i, z^i)^T$ is the value of the $i$-th iteration of the discretization equation (5.1). The time increment $\Delta t > 0$. $\delta_{i}$ is a random variable which obeys the Gaussian distribution $\mathbb{N}(0, 1)$ for $i = 1, 2, \ldots, n$. Let the initial value be $(S(0), I(0), R(0), z(0)) = (0.5, 0.5, 0.5, 0.5)$.

We can observe and analyze the results under different conditions by slightly adjusting a few key parameters. From the expressions of the threshold value $R_0^*$, it can be noted that $\bar{\beta}$, $a$, $\eta$, $\sigma$, and $K$ have a significant effect on the dynamical behavior. According to the variable-controlling approach, we will separately study the impact of these parameters on threshold dynamics of model (1.7). In particular, we will concentrate on two scenarios of the disease: persistence and extinction.

5.1.   The impact of $\bar{\beta}$

The aim of this part is to study how $\bar{\beta}$ affects the spread of the disease in system (1.7). The numerical simulations use the following parameters:

Example 5.1.1 (Persistence) The trend of $R_0^*$ is outlined in Figure 1. It is easy to see that system (1.7) will exist at least one stationary distribution and the solution follows a normal density function when $\bar{\beta}\in(0.2587, 0.5]$.

Consider the following three cases: (ⅰ) $\bar{\beta} = 0.3$; (ⅱ) $\bar{\beta} = 0.4$; (ⅲ) $\bar{\beta} = 0.5$. Figure 2 reflects the sample paths of $S(t)$, $I(t)$ and $R(t)$ of system (1.7) with different $\bar{\beta}$. When $\bar{\beta} = 0.3$, the top of Figure 3 depicts the solutions $S(t)$, $I(t)$ and $R(t)$ of stochastic system (1.7) and the corresponding deterministic system.

We can observe that stochastic solution fluctuates around the deterministic solution. The frequency histograms of $S(t)$, $I(t)$ and $R(t)$ with 30 rectangular bars are displayed in the bottom of Figure 3. The frequency histograms can clearly show the frequency distribution of different values, which provides a more comprehensive understanding of the distribution characteristics.

Example 5.1.2 (Extinction) We consider the situation of $\bar{\beta}\in[0, 0.2587)$. For the following three cases: (ⅰ) $\bar{\beta} = 0.05$; (ⅱ) $\bar{\beta} = 0.15$; (ⅲ) $\bar{\beta} = 0.25$, it is evident from Figure 4 that the disease will eventually disappear with the passage of time.

In brief, from Figure 2, the numbers of infectious individuals $I$ and recovered individuals $R$ will decrease when $\bar{\beta}$ decreases. We also find that the rate of disease extinction will be faster with small $\bar{\beta}$ from Figure 4. The small $\bar{\beta}$ plays an effective role in preventing the spread of the disease. To reduce incidence rate, it is necessary to formulate reasonable prevention and control strategies. Furthermore, there is a need to reinforce personal protection and health education.

5.2.   The impact of $a$

In this part, the impact of the disease-related death rate $a$ on the dynamical behavior of system (1.7) is analyzed. Let $a$ be a variable, and the other parameters are given fixed values as follows:

Example 5.2.1 (Persistence) The focus of this example is on the persistence of the disease, that is, the case where $a$ is situated within the interval $[0, 0.2005)$. This is illustrated in Figure 5.

From Figure 6, we can get information about the numbers of $S$, $I$ and $R$ with $a = 0.05, 0.1, 0.15$ and when $a$ gradually increases, the number of infected individuals and recovered individuals will decrease. According to Theorem 3.2, the solution $(N(t), I(t), R(t), \log\beta(t))$ follows a normal density function $\Phi_1\sim\mathbb{N}_4((1.6667, 0.3333, 0.0833, \log 0.4), \tilde{\varSigma}^1_*)$ when $a = 0.1$, where

Our attention is directed toward the marginal density functions of $N$, $I$, and $R$, which is expressed as follows:

Figure 7 plots the frequency histograms fitting curves and marginal density functions curves. It indicates that the marginal density functions curves and the corresponding frequency histogram fitting curves exhibit coincidence.

Example 5.2.2 (Extinction) As illustrated in Figure 5, the disease becomes extinct when $a\in(0.2005, 0.5]$. Figure 8 presents the solutions $S(t)$, $I(t)$ and $R(t)$ of system (1.7) with $a = 0.3$. The results show that the numbers of infected individuals and recovered individuals converge to zero.

5.3.   The impact of $\eta$

This part will explore the influence of the transfer rate $\eta$ from the infected individuals to the susceptible individuals on disease transmission of system (1.7). Assume that the parameters are given by

The trend of $R_0^*$ with different $\eta\in[0, 0.5]$ is presented in Figure 9. Specifically, if $\eta\in[0, 0.2913)$, system (1.7) will admit at least one stationary distribution and the solution follows a normal density function; if $\eta\in(0.2913, 0.5]$, the disease will be extinct. Four cases are considered: (ⅰ) $\eta = 0.1$; (ⅱ) $\eta = 0.2$; (ⅲ) $\eta = 0.3$; (ⅳ) $\eta = 0.4$. Figure 10 depicts the corresponding computer simulations for the stochastic solutions $S(t)$, $I(t)$, and $R(t)$ with different $\eta$. It implies that the disease of system (1.7) will persist under cases (ⅰ), (ⅱ) and will be exponentially extinct under cases (ⅲ), (ⅳ). One draws the conclusion that large $\eta$ plays a significant role in curbing the spread of the disease. In case (ⅰ), there exists a density function $\Phi_2\sim\mathbb{N}_4((1.4685, 0.5315, 0.0266, \log 0.5), \tilde{\varSigma}^2_*)$, where

We derive the following marginal density functions of $N$, $I$ and $R$

Figure 11 exhibits the frequency histograms and marginal density functions curves of $N$, $I$, and $R$ in system (1.7).

5.4.   The impact of $K$

In this part, we study the impact of reversion speed $K$ on the long-term behavior of system (1.7). To this end, the parameters of the system are shown as follows:

Example 5.4.1 (Persistence) To guarantee the persistence of the disease, we set $(\bar{\beta}, \sigma) = (0.3, 0.5)$. Three cases of reversion speed are selected for analysis: (ⅰ) $K = 1$; (ⅱ) $K = 2$; (ⅲ) $K = 3$. By calculating the value of $R_0^*$ in these cases, we find that the results are 1.9668, 1.9448 and 1.9376, which all meet the important condition, that is, $R_0^* > 1$. Figure 12 demonstrates the time evolutions of $S(t)$, $I(t)$, and $R(t)$ of deterministic system (1.1) and stochastic system (1.7) under the above cases. It is possible to observe that the solution of the stochastic model approaches that of the deterministic model as the speed reversion increases. According to Theorem 3.2, when $K = 1$, the solution follows the normal density function $\Phi_4\sim\mathbb{N}_4((1.3725, 0.6275, 0.0314, \log 0.3), \tilde{\varSigma}^4_*)$, where

Clearly, it has the following three marginal density functions

We run the numerical simulation 200,000 times and plot the marginal density functions and frequency histograms of $N$, $I$ and $R$. The results are shown in Figure 13.

Example 5.4.2 (Extinction) We rechoose $(\bar{\beta}, \sigma) = (0.8, 0.2)$. When $K = 1.5, 2.5, 3.5$, the simple calculation obtains the value of $R_0^* = 0.8558, 0.8201, 0.8052$. It implies that the disease of system (1.7) will die out in a long term. Figure 14 presents computer simulations for the numbers of susceptible individuals $S$, infected individuals $I$, and recovered individuals $R$ with $K = 1.5, 2.5, 3.5$.

One can see that stochastic fluctuation will go up with small $K$ from Figures 12 and 14. It means that speed reversion has a great destabilizing influence on the spread of the disease.

5.5.   The impact of $\sigma$

In this subsection, we will delve into how volatility intensity $\sigma$ affects the long-term stability of system (1.7). The parameter values for numerical simulations are shown as follows:

Example 5.5.1 (Persistence) Let $\bar{\beta} = 0.5$. We choose three cases of the volatility intensity: (ⅰ) $\sigma = 0.1$; (ⅱ) $\sigma = 0.2$; (ⅲ) $\sigma = 0.3$. In each of the aforementioned cases, the value of $R_0^*$ is calculated to be 1.9279, 1.9424, 1.9668, which all satisfy $R_0^* > 1$. Figure 15 presents the solutions of $S(t)$, $I(t)$, and $R(t)$ in stochastic system (1.7) and the corresponding deterministic system. It can be seen that the solution of the stochastic model closes to that of the deterministic model when volatility intensity is small. In addition, according to Theorem 3.2, we can investigate the probability density function near the positive equilibrium. When $\sigma = 0.1$, there exists a density function $\Phi_3\sim\mathbb{N}_4((1.3725, 0.6275, 0.0314, \log 0.5), \tilde{\varSigma}^3_*)$, where

In order to further deepen our understanding of the properties of density function, we calculate three two-dimensional joint density functions. Figure 16 provides an intuitive view of $\bar{\Phi}_3(N, I)$, $\tilde{\Phi}_3(N, R)$ and $\check{\Phi}_3(I, R)$. The specific expressions of joint density functions are shown as follows:

Example 5.5.2 (Extinction) To visually demonstrate extinction phenomenon of the disease, we choose $\bar{\beta} = 0.2$ and calculate directly the value of $R_0^*$. When $\sigma = 0.3, 0.6, 0.9$, the results are 0.7867, 0.8417, 0.9419. According to Theorem 4.1, the numbers of infectious individuals $I$ and recovered individuals $R$ will eventually be extinct in the long term, which is supported by Figure 17. In summary of this part, taking a closer look at Figures 15 and 17, we can observe an interesting phenomenon: the larger the volatility intensity is, the larger the fluctuation of the stochastic solution in system (1.7) is.

From Subsections 5.4 and 5.5, we find that $R_0^*\rightarrow R_0$ when $K\rightarrow \infty$ or $\sigma\rightarrow 0$. That is, the stochastic model and the corresponding deterministic model have similar properties if reversion speed $K$ approaches infinity or volatility intensity $\sigma$ tends to zero.

5.6.   The color phase diagram of $R_0^*$

The purpose of this example is to explore the influence of reversion speed and volatility intensity on the threshold $R_0^*$. Let $A = 0.2, d = 0.1, \eta = 0.1$, $\gamma = 0.1, a = 0.1, c = 0.01, \beta = 0.2, K = 1, \sigma = 0.3.$ Figure 18 plots two-dimensional diagram of $R_0^*$ in $(K, \sigma)\in[0.02, 0.05]\times[0.02, 0.05]$. It is easy to find that $R_0^*$ will rise with small reversion speed or large volatility intensity.

6.   Conclusions

In the course of a comprehensive investigation into the dynamics of epidemic models, it is recognized that the impact of environmental noise cannot be ignored. The intricate environmental conditions often lead to fluctuations in incidence rate, which can exert a considerable influence on the velocity and extent of disease transmission. In light of this, we assume that the incidence rate satisfies the log-normal OU process. This paper proposes a stochastic SIRS epidemic model with standard incidence rate and transfer from infected individuals to susceptible individuals. The main findings of this paper are listed as follows:

● The sufficient conditions for the existence of stationary distribution and extinction of the disease are established. There exists a threshold between extinction and persistence of the disease

If $R_0^* > 1$, then system (1.7) has at least one stationary distribution. It means that the disease will persist. If $R_0^* < 1$, then the disease will be extinct in a long time.

● If $R_0^* > 1$, then the solution of system (1.7) near $\check{E}$ = ($N^*$, $I^*$, $R^*$, $\log \bar{\beta}$) has a normal density function. Specifically,

● Numerical simulations show some conclusions: (a) Larger $\bar{\beta}$ makes the numbers of infectious individuals and recovered individuals increase, which is not conducive to the extinction of the disease. (b) Small disease-related death rate $a$ and the transfer rate $\eta$ from the infected individuals to the susceptible individuals have a negative effect on prevention and control of an infectious disease. (c) Reversion speed $K$ and volatility intensity $\sigma$ are the primary factors influencing the stability of the stochastic epidemic model. In particular, the dynamic behavior of the stochastic model is consistent with that of the corresponding deterministic model when the reversion speed is considerable or the volatility intensity is sufficiently minimal.

Some interesting topics deserve further consideration. On the one hand, it is well-known that time delays can be used to describe the incubation period of the infectious disease, the period of infection of a patient with the disease and the period of immunization of a recovering person against the disease. We can further study a stochastic epidemic model with time delays and the log-normal OU process. On the other hand, it is also worth it to study whether or not the method used in this paper can be applied to other stochastic epidemic models.

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 would like to thank the anonymous referee for his/her useful suggestions.

Conflict of interest

The authors declare there is no conflicts of interest.