key: cord-0899403-1t4gl53f
authors: Zhou, Baoquan; Jiang, Daqing; Dai, Yucong; Hayat, Tasawar
title: Stationary distribution and density function expression for a stochastic SIQRS epidemic model with temporary immunity
date: 2021-06-08
journal: Nonlinear Dyn
DOI: 10.1007/s11071-020-06151-y
sha: e1054db7a217781dc84364118fe6591002a8bace
doc_id: 899403
cord_uid: 1t4gl53f

Recently, considering the temporary immunity of individuals who have recovered from certain infectious diseases, Liu et al. (Phys A Stat Mech Appl 551:124152, 2020) proposed and studied a stochastic susceptible-infected-recovered-susceptible model with logistic growth. For a more realistic situation, the effects of quarantine strategies and stochasticity should be taken into account. Hence, our paper focuses on a stochastic susceptible-infected-quarantined-recovered-susceptible epidemic model with temporary immunity. First, by means of the Khas’minskii theory and Lyapunov function approach, we construct a critical value [Formula: see text] corresponding to the basic reproduction number [Formula: see text] of the deterministic system. Moreover, we prove that there is a unique ergodic stationary distribution if [Formula: see text] . Focusing on the results of Zhou et al. (Chaos Soliton Fractals 137:109865, 2020), we develop some suitable solving theories for the general four-dimensional Fokker–Planck equation. The key aim of the present study is to obtain the explicit density function expression of the stationary distribution under [Formula: see text] . It should be noted that the existence of an ergodic stationary distribution together with the unique exact probability density function can reveal all the dynamical properties of disease persistence in both epidemiological and statistical aspects. Next, some numerical simulations together with parameter analyses are shown to support our theoretical results. Last, through comparison with other articles, results are discussed and the main conclusions are highlighted.

Over time, an increasing number of people are becoming concerned with health and the desire to improve the quality of life worldwide. However, major infectious diseases such as Ebola, avian influenza, cholera, and heptitis B are one of the biggest threats to public health [1] [2] [3] . Epidemiology, greatly supported by various mathematical models, is the study of the spread of diseases and trace factors that give rise to their occurrence. Following the classical Susceptible-Infected-Recovered (SIR) epidemic models proposed by Kermack and McKendrick [4] , some authors have developed a series of reasonable ordinary differential equations (ODEs) to describe the transmission of various epidemics [5] [6] [7] [8] [9] [10] [11] [12] [13] . In [5] , Liu et al. established an Susceptible-Vaccinated-Infected-Recovered (SVIR) epidemic model with vaccination strategies and obtained the corresponding global stability of equilibria. Hove-Musekwa and Nyabadza [9] considered a deterministic HIV/AIDS model taking into account the active screening of disease carriers and seeking of treatment. They also derived the relevant basic reproduction number. Considering the sequence diversity and highly infectious nature of some contagious diseases, we occasionally need to implement quarantine strategies to control the spread of disease. For example, without an effective vaccine, the rapid spread of coronavirus disease 2019 (COVID-19) worldwide has had a serious socioeconomic impact and imposed a potentially great threat to human safety [14] . Hence, many suitable SIQR epidemic models have been developed in the past few decades [15] [16] [17] [18] [19] . In [15] , Herbert and Ma obtained the corresponding basic reproduction number of a deterministic SIQR model with quarantineadjusted incidence. Nevertheless, recovered individuals with temporary immunity may be susceptible to the disease again in the future [20] [21] [22] [23] . Zhang et al. [20] studied the global asymptotic stability of two equilibria to a SIQS epidemic model with the nonlinear incidence rate β S I f (I ) . Following the above analyses, in this study, a deterministic SIQRS with temporary immunity is developed for further epidemiological investigation.

In our daily life, it is obvious that the spread of infection, travel of populations and the design of control strategies are critically perturbed by some environmental variations [24] . For dynamical study and simulation, by taking the effect of stochastic perturbation into consideration, some scholars considered and analyzed various stochastic differential equations (SDEs) for the spread of epidemics [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] . From [25] , Zhao and Jiang established a universal theory of extinction and persistence in mean based on a stochastic SIS epidemic model with vaccination. Han and Jiang [27] introduced a stochastic staged progression AIDS model with second-order perturbation and proved the ergodicity of the global positive solution if R H 0 > 1. Considering the delay influence, Caraballo and Fatini [29] derived the existence of stationary distribution for a stochastic SIRS epidemic model with distributed delay. For cholera epidemic, a stochastic SIQRB infectious disease model was researched by Liu and Jiang [32] . Recently, Zhou and Zhang (2020) obtained the explicit expression of the probability density function for the three-dimensional avian-only influenza model, which is described in [33] .

Focusing on the temporary immunity phenomenon of infected people, quarantine strategies and random perturbations, our study aims to develop a stochastic SIQRS epidemic system with temporary immunity. As is well known, the corresponding basic reproduction number and unique endemic equilibrium can reflect the disease persistence of a deterministic system. Nevertheless, the positive equilibrium no longer exists in a stochastic model owing to the effect of unpredictable environmental noises. Hence, ergodicity theory and the existence of stationary distribution, which greatly reflect the stochastic permanence of disease, are gradually becoming more popular in the transmission of epidemics. In practical application, some statistical properties of epidemic models still need to be estimated to effectively prevent and control the spread of infectious diseases. Notably, there are relatively few studies devoted to deriving the explicit expression of probability density function due to the difficulty of solving the high-dimensional Fokker-Planck equation. To the best of our knowledge, some studies of probability density functions for stationary distributions are shown in the present study. As a result, we will concentrate on the following three aims: (i) construct a reasonable stochastic threshold R S 0 corresponding to the basic reproduction number R 0 ; (ii) investigate the disease persistence of stochastic SIQRS model under R S 0 > 1, namely, the existence of the uniqueness of an ergodic stationary distribution and the exact expression of this unique probability density function; and (iii) provide the corresponding numerical simulations and parameter analyses for our analytical results.

The rest of our study is arranged as follows. For the later dynamical investigation, Sect. 2 introduces the corresponding mathematical models, important notations and necessary lemmas. By constructing a series of suitable Lyapunov functions, a stochastic critical value R S 0 involved in the random noises is obtained. Based on the global positive solution property and Khas'minskii theory, Sect. 3 shows that there is a unique ergodic stationary distribution when R S 0 > 1. In Sect. 4, by developing some solving theories of the relevant algebraic equations, the corresponding fourdimensional Fokker-Planck equation is solved for the explicit expression of log-normal density function to the stochastic model if R S 0 > 1. Section 5 presents some empirical examples and parameter analyses to validate the above theoretical results. Finally, the relevant results are discussed and the main conclusions are introduced by comparison with the existing results in Sect. 6.

Given the above descriptions, we assume that the investigated population N (t) can be divided into susceptible S(t), infectious I (t), quarantined Q(t) and recovered R(t) individuals at time t. A deterministic SIQRS epidemic model with temporary (short-term) immunity is studied herein, which is given by

where Λ is the recruitment rate of the susceptible individuals, μ depicts the natural death rate of the population, β is the transmission rate, α 1 , α 2 represent the average disease-induced death rate of the infected and quarantined individuals, respectively, δ denotes the isolated rate of the infected individuals, γ and ε are the recovery rate of the infected and quarantined individuals, and ω denotes the immune loss rate of the recovered individuals. All the parameters are assumed to be positive constants. In the similar methods described by Ma and Zhou [21] , system (2.1) has the corresponding basic reproduction number and the invariant attracting set, which are given by

Additionally, two possible equilibria are shown as follows. (i) The disease-free equilibrium E 0 = (S 0 ,

(μ+ω)(μ+α 1 +ε) > 0, 2 = γ (μ+α 2 + ε) + εδ > 0. These two equilibria have the following dynamical properties.

• If R 0 ≤ 1, then E 0 is globally asymptotically stable in Θ, which means the disease will be eradicated in a population.

• If R 0 > 1, then E * is globally asymptotically stable, but E 0 is unstable in the domain Θ. This indicates the disease will prevail and persist long-term.

In reality, the dynamical behavior of most epidemics is inevitably affected by random factors in nature. By means of the relevant assumptions and forms of stochastic perturbations developed in [25] [26] [27] [28] [29] [30] [32] [33] [34] , in this study, we assume that these stochastic noises are directly proportional to the groups S(t), I (t), Q(t) and R(t). Hence, the corresponding stochastic SIQRS epidemic model with temporary immunity is described by

and B 4 (t) are four independent standard Brownian motions, and σ 2 i > 0 (i = 1, 2, 3, 4), respectively, denote their intensities.

Throughout this study, unless otherwise specified, let { , F , {F t } t≥0 , P} be a complete probability space with a filtration {F t } t≥0 satisfying the usual conditions (i.e., it is increasing and right continuous while F 0 contains all P-null sets). For detailed descriptions, refer to Mao [36] . Moreover, for convenience and simplicity, all stochastic approaches and theories are based on the above space.

In order to study the later dynamical behavior of the stochastic system (2.2), some common notations shall be defined in the first place. Let R n be an n-dimensional Euclidean space and

In addition, let A τ be the transposed matrix of the inverse matrix A, and A −1 be the relevant inverse matrix of A. Next, the corresponding global existence of the positive solution to the system (2.2) is introduced as follows.

2) on t ≥ 0, and the solution will remain in R 4 + with probability one (a.s.).

The detailed proof of Lemma 2.1 is mostly similar to that in Theorem 3.1 of Liu and Jiang [37] , so we omit it here. By means of the Khas'minskii theory [38] , considering the following stochastic differential equation defined in the space R n ,

where the diffusion matrix F(X ) = (ā i j (X )), and a i j (X ) = n k=1 σ i k (X )σ j k (X ). Furthermore, the relevant existence theory of the unique ergodic stationary distribution is shown by the following Lemma 2.2. [38] ) The Markov process X (t) has a unique ergodic stationary distribution (·) if there exists a bounded domain D ⊂ R n with a regular boundary Γ and (A 1 ). There is a positive number κ 0 such that

Then for all x ∈ R n and integral function φ(·) with respect to the measure φ(·), it follows that P lim

Now, by the relevant definitions described in Zhou and Zhang [33] , we will develop some solving theories for the corresponding four-dimensional algebraic equations, which are described by the following Lemmas 2.3-2.5.

Assuming that a 1 > 0, a 3 > 0, a 4 > 0, and a 1 (a 2 a 3 − a 1 a 4 ) − a 2 3 > 0, then θ 0 is a positive definite matrix.

Let θ 1 be a symmetric matrix in the fourdimensional algebraic equation

Lemma 2.5 Let θ 2 be a symmetric matrix in the fourdimensional algebraic equation G 2 0 +C 0 θ 2 +θ 2 C τ 0 = 0, where G 0 = diag(1, 0, 0, 0), and

If c 1 > 0 and c 2 > 0, then θ 2 is semi-positive definite.

Remark 2.6 For convenience, A 0 , B 0 and C 0 are, respectively, called standard R 1 , R 2 , R 3 matrices in this study. The corresponding proofs of Lemmas 2.3-2.5 are separately shown in subsections (I), (II) and (III) of "Appendix A".

In this section, by Lemmas 2.1 and 2.2, we are devoted to obtain the sufficient conditions for ergodicity of the global positive solution and the existence of stationary distribution. First, we define

Theorem 3.1 Assuming that R S 0 > 1, for any initial value (S(0), I (0), Q(0), R(0)) ∈ R 4 + , then the solution (S(t), I (t), Q(t), R(t)) of the system (2.2) is ergodic and has a unique stationary distribution (·).

Proof By Lemma 2.1, we derive that there is a unique global positive solution (S(t), I (t), Q(t), R(t)) ∈ R 4 + . Hence, the proof of Theorem 3.1 is divided into the following three steps: (i) construct a series of Lyapunov functions to derive a suitable non-negative C 2 -function V (S, I, Q, R) and a stochastic critical value R s 0 related to R 0 ; (ii) establish a reasonable bounded domain D and prove the assumption (A 2 ) of Lemma 2.2; and (iii) validate the condition (A 1 ) of Lemma 2.2.

Step 1 Define an important C 2 -function V (S, I, Q, R) by 

For simplicity, we let

By means of Itô's formula which is shown in "Appendix C", the function V 1 satisfies

Employing Itô's formula to V 2 , one has

Similarly, by the definition of λ, we have Consequently, we can construct a suitable non-negative C 2 -function V (S, I, Q, R) in the following form

Combining (3.1)-(3.4) and the definition ofλ, one can see that

(3.5)

Step 2 Consider the following bounded set

where > 0 is a sufficiently small constant such that the following inequalities hold.

(3.8)

For simplicity, let X (t) = (S(t), I (t), Q(t), R(t)) τ . Consider the following seven subsets of R 4

Clearly, R 4 + \ D = 7 k=1 D k, . By (3.6)-(3.9), we can derive

Given the above, we can therefore obtain a pair of sufficiently small > 0 and closed domain D such that

Hence, the assumption (A 2 ) of Lemma 2.2 holds.

Step 3 System (2.2) has the corresponding diffusion matrix

Obviously, for any (S, I, Q, R) ∈ D , F is a positive definite matrix. In other words, we can determine a positive constant

Then the condition (A 1 ) of Lemma 2.2 also holds. Given the above three steps, system (2.2) admits a unique ergodic stationary distribution (·) with respect to the solution (S(t), I (t), Q(t), R(t)). The proof of Theorem 3.1 is confirmed.

Remark 3.2 From the expressions of R 0 and R S 0 , we can easily obtain that R S 0 ≤ R 0 . As we know, the existence of an ergodic stationary distribution denotes the stochastic positive equilibrium state. Hence, R S 0 > 1 can be regarded as the unified criterion which guarantees the disease persistence of the deterministic model (2.1) and stochastic system (2.2). Furthermore, R S 0 = R 0 while σ 1 = σ 2 = 0. This means that the disease persistence is critically affected by the random fluctuations of susceptible and infected individuals rather than quarantined and recovered individuals.

By Theorem 3.1, we obtain that system (2.3) has a unique stationary distribution which has ergodic property if R S 0 > 1. For further development of infectious disease dynamics, in this section, we will study the corresponding probability density function of the distribution (·) to derive all statistical properties of system (2.3). Before this, two equivalent differential equations of system (2.2) should be firstly introduced.

, and x 4 = ln(R). By means of Itô's formula, system (2.2) can be transformed into the following equation:

Assuming that R S 0 > 1 and following the description of Zhou and Zhang [33] , we similarly define a quasi-endemic equilibrium E *

, e x * 4 ) ∈ R 4 + , which is determined by the following algebraic equations:

By detailed calculation, we obtain that S *

In fact, the quasi-endemic equilibrium E * + is the same as E * if there is no stochastic perturbation. This is the reason why the quasi-endemic equilibrium is defined.

then the corresponding linearized differential equation of system (4.1) is given by

By the definition of E * + , we easily obtain that all the parameters in (4.4) are positive constants. Next, we will study the corresponding probability density function around the quasi-endemic equilibrium E * + .

Theorem 4.1 Assuming that R S 0 > 1, for any initial value (S(0), I (0), Q(0), R(0)) ∈ R 4 + , the solution (S(t), I (t), Q(t), R(t)) of the system (2.2) follows the unique log-normal probability density function Φ(S, I, Q, R) around the quasi-endemic equilibrium E * + , which is described by

where is a positive definite matrix, and the special form of is given as follows.

(3) If w 1 = 0 and w 2 = 0, then

Moreover, the standardized transformation matrices Proof For the sake of convenience, by letting

then the linearized system (4.3) can be simplified to dY = AY dt + Gd B(t). According to the relevant theory of Gardiner [39] , there exists a unique probability density function Φ(y 1 , y 2 , y 3 , y 4 ) around the quasiendemic equilibrium E * + , which is determined by the following Fokker-Planck equation

Considering the diffusion matrix G is a constant matrix, from the results of Roozen [40] , Φ(y 1 , y 2 , y 3 , y 4 ) can be described by a Gaussian distribution. In other words,

where φ 0 is a positive constant, which is obtained by the normalized condition R 4 Φ(Y )dy 1 dy 2 dy 3 dy 4 = 1.

The real symmetric matrix Q satisfies

Then it is equivalent to the following equation

By the finite independent superposition principle, we only need to consider the corresponding solutions of the following four algebraic sub-equations:

Finally, we can obtain that = 1 + 2 + 3 + 4 and G 2 = G 2 1 + G 2 2 + G 2 3 + G 2 4 . Before proving the positive definiteness of , we firstly verify that all the eigenvalues of A have negative real parts. The corresponding characteristic polynomial of A is defined by

where (i) a 1 = a 11 + a 32 + a 42 + a 43 > 0, a 2 = (a 11 + a 32 )(a 42 + a 43 ) + a 11 a 32 + a 12 a 21 > 0, (ii) a 3 = (a 42 +a 43 )(a 11 a 32 +a 12 a 21 )+a 12 a 21 a 32 − a 14 a 21 a 42 , a 4 = a 21 a 32 (a 12 − a 14 )(a 42 + a 43 ). From the expressions of a 11 , a 12 and a 14 , one can see that

Moreover, by the corresponding simplicity, we have 

In view of the Routh-Hurwitz stability criterion [41] , we can obtain that A has all negative real-part eigenvalues. According to the matrix similar transformation theory, we realize that ϕ A (λ) is a similarity invariant, which means that a 1 , a 2 , a 3 and a 4 are also similarity invariants.

Next, the corresponding proof for the positive definiteness of is divided into four steps. More precisely, we will show that 3 , 4 are both positive definite, and 1 , 2 are both at least semi-positive definite.

Step 1 Consider the algebraic equation

For the following elimination matrix H 1 , by letting . Using the value of w 1 , we analyze the following two cases.

Case 1 If w 1 = 0, based on the method introduced in the subsection (I) of "Appendix B", let

with m 1 = w 1 [(a 32 + a 42 + a 43 )(a 42 + a 43 ) + a 2 32 ]. By direct calculation, we derive

Moreover, Eq. (4.9) can be equivalently transformed into

where ρ 1 = a 21 a 32 w 1 σ 1 , it can be simplified as

Since A has all negative real-part eigenvalues, B 1 is a standard R 1 matrix. According to Lemma 2.3, 0 is positive definite. In the similar results described in subsection (I) of "Appendix A", the form of 0 is given as follows.

(4.12)

which is similarly obtained by the method in subsection (II) of "Appendix B". Then B 1w 1 is

where b i (i = 1, 2, 3, 4) are shorthands, and we are only concerned with their signs. Meanwhile, (4.9) can be equivalently transformed into the following equation:

where ρ 1w 1 = a 21 a 32 σ 1 . We similarly obtain G 2 0 + B 1 0 + 0 B τ 1 = 0. By means of the uniqueness of ϕ A (λ), one can see that

= a 11 a 32 + a 12 a 21 > 0, Furthermore, we can compute that Based on (i)-(iv) and Lemma 2.4, B 1w 1 is a standard R 2 matrix, and 0 is a semi-positive definite matrix. Considering the relevant results introduced in subsection (II) of "Appendix A", 0 will be

(4.14)

Therefore, given the above cases, we get that 1 is at least semi-positive definite.

Step 2 For the following algebraic equation

consider the corresponding order matrix J 2 and elimination matrix P 2 :

Then we can obtain that

where w 1 , w 2 are the same as those in Theorem 4.1.

Similarly, based on the values of w 1 and w 2 , the following four sub-cases shall be analyzed.

Case (I 1 ) For the following elimination matrix H 2 , let C 2 = H 2 B 2 H −1 2 :

where w 3 = a 14 + (a 42 +a 43 −a 11 )w 2 w 1 . Based on the value of w 3 , we will discuss the following two sub-cases of Case (I 1 ) • If w 3 = 0, in the similar method as Case 1 of Step 1 described, we construct the following standardized transformation matrix:

where m 2 = w 3 [(a 11 + a 42 + a 43 )(a 42 + a 43 ) + a 2 11 ]. Considering the similar invariant properties of a 1 , a 2 , a 3 , a 4 , we can clearly derive that the standard R 1 matrix of A is unique. From the results shown in subsection (I) of "Appendix B", let D 2 = M 2 C 2 M −1 2 . Then we obtain that D 2 = B 1 , and (4.15) can be transformed into

• If w 3 = 0, the relevant standardized transformation matrix is given as follows:

By simple computation, we have

Similarly, b i (i = 1, 2, 3, 4) are shorthands, and we only focus on their signs. According to the method presented in Case 2 of Step 1, the characteristic polynomial A is described by Hence, D 2w 3 is a standard R 2 matrix. Meanwhile, we can transform (4.15) into the following equation:

where 0 = ρ −2 2w 3 (M 2w 3 H 2 P 2 J 2 ) 2 (M 2w 3 H 2 P 2 J 2 ) −1 with ρ 2w 1 = a 32 w 1 σ 2 . By means of Lemma 2.4 and the above results (i)-(iv), it is found that 0 is semi-positive definite. Following the detailed results described in subsection (II) of "Appendix A", we derive

(4.21)

Thus,

Case (I 2 ) If w 1 = 0 and w 2 = 0, consider the following standardized transformation matrix: where m 3 = a 14 [(a 11 + a 42 + a 43 )(a 42 + a 43 ) + a 2 11 ]. Let D 2w 2 = M 2w 2 B 2 M −1 2w 2 . From the descriptions in subsection (I) of "Appendix B", then D 2w 2 is also a standard R 1 matrix. Based on the uniqueness of the standard R 1 matrix of A, D 2w 2 = B 1 . Thus, (4.15) can be converted to the equivalent equation

Consequently, 2 is a positive-definite matrix, and

Case (I 3 ) If w 1 = 0 and w 2 = 0, for the following order matrix J 1 , by letting C 2w 1 = J 1 B 2 J −1 1 , we find that We can then define D 2w 1 = M 2w 1 C 2w 1 M −1 2w 1 , where the standardized transformation matrix M 2w 1 is given by

(4.24)

Based on subsection (II) of "Appendix B", D 2w 1 is also a standard R 2 matrix, which satisfies

which is the same as in Case 2 of Step 1. Thus, (4.15) can be transformed into the following equation:

By the property of 0 , then 2 With Lemma 2.5, we obtain that D 2w 12 is a standard R 3 matrix. Additionally, by letting¯ 0 = ρ −2 2w 12 (M 2w 12 P 2 J 2 ) 2 (M 2w 12 P 2 J 2 ) −1 with ρ 2w 12 = a 32 σ 2 , then (4.15) can be equivalently transformed into the following equation: (4.28)

Consequently,

Given the above analyses, 2 is at least a semi-positive definite matrix.

Step 3 Consider the following algebraic equation:

where the relevant order matrix J 3 is the same as that described in Theorem 4.1, and standardized transformation matrix M 3 is given by Similarly, we get that B 3 is a standard R 1 matrix. Given the uniqueness, then B 3 = B 1 . Furthermore, Eq. (4.29) equivalently converts to the following equation:

By the property of 0 , then 3 

Step 4 Consider the following algebraic equation: Obviously, B 4 is also a standard R 1 matrix, which means that B 4 = B 1 . Similarly, (4.32) is equivalent to the following algebraic equation: = 1 + 2 + 3 + 4 is a positive definite matrix, and the corresponding special form of can be determined by the above steps. By the positive definiteness of , we can compute that φ 0 = (2π) −2 | | − 1 2 . Thus, there exists a unique exact normal density function around the quasi-endemic equilibrium E * + while R S 0 > 1, which is given by Φ(y 1 , y 2 , y 3 ,

Considering the transformation of (y 1 , y 2 , y 3 , y 4 ) and (S, I, Q, R), then the unique ergodic stationary distribution (·) of system (2.2) approximately follows a unique log-normal probability density function This completes the proof of Theorem 4.1.

we can obtain that the unique ergodic stationary distribution (·) of system (2.2) admits the corresponding probability density function Φ(S, I, Q, R). Hence, R S 0 > 1 can be considered as a reasonable stochastic criterion for the disease persistence, and the exact expression of the density function of system (2.2) can provide an effective method to prevent and control many epidemics.

In this section, we will introduce some examples and simulations to illustrate the above theoretical results. Making use of the common higher-order method developed by Milstein [41] , the corresponding discretization equation of system (2.2) is given by

where the time increment is t > 0, and ξ k , η k , ζ k , ν k are the independent Gaussian random variables which follow the Gaussian distribution N (0, 1) for k = 1, 2, . . . , n. From the realistic statistics described by Hethcote [15] , Qi [34] and the Central Statistical Office of Zimbabwe (CSZ), the corresponding biological parameters and initial value of system (2.2) are shown in Table 1 .

Next, based on Table 1 , we will perform some empirical examples to focus on the following four aspects: (i) The ergodicity property and the existence of a unique stationary distribution if R S 0 > 1. (iv) The effects of the main parameters of system (2.2) on the disease dynamics, such as the recruitment rate and transmission rate. In addition, we still give the exact expression of a unique log-normal density function Φ(S, I, Q, R) with respect to the distribution (·).

Example 5.1 According to Table 1 , let the main parameters (Λ, β, μ, α 1 , α 2 , δ, γ, ε, ω) = (25000, 4 × 10 −6 , 0.0143, 0.02, 0.01, 0.2, 0.2, 0.15, 0.25) and let the stochastic perturbations (σ 1 , σ 2 , σ 3 , σ 4 ) = (0.008, 0.004, 0.005, 0.005). We then calculate that By Theorem 3.1, we can obtain that system (2.3) has a unique stationary distribution (·), which means the epidemic will be persistent long-term. This is supported by the left column of Fig. 1 . Furthermore, by detailed calculation, we derive that

It follows from Theorem 4.1 that the stationary distribution (·) obeys a log-normal density function Φ(S, I, Q, R). More precisely, 

The curves of (i)-(iv) are shown in the right column of Fig. 1 . Clearly, it verifies Theorem 4.1 from the side.

It follows from Remark 4.1 that random fluctuations of susceptible and infected individuals (i.e., σ 1 , σ 2 ) have a critical influence on the disease persistence. Therefore, by the method of controlling variables, we are devoted to studying the corresponding dynamical effects of σ 1 and σ 2 in the following Example 5.2.

Example 5.2 First, the parameters are chosen considering the following three subcases of random perturbations:

In fact, the above three subconditions (a 1 )-(a 3 ) all guarantee the existence of the ergodic stationary distribution of system (2.2). For subcases (a 1 ) and (a 2 ), i.e., sub- figure (2-1) , by only increasing the perturbation intensity of susceptible individuals, the corresponding numbers of quarantined and recovered individuals decrease apparently. In contrast, by only increasing the perturbation intensity of infected individuals, the population From the expression of R S 0 (or R 0 ), the disease persistence of system (2.2) (or system (2.1)) is critically affected by the recruitment rate Λ, transmission rate β and quarantined coefficient δ. Hence, the following Examples 5.3-5.5 will reveal these effects. Fig. 3 , respectively. Clearly, the epidemic infection will decrease as the recruitment rate decreases. Moreover, a small recruitment rate, such as Λ < 10,000, can effectively lead to disease extinction (see Fig. 3 ). Consider the subcases of transmission rate β = 3 × 10 −7 , 5 × 10 −7 , 7 × 10 −7 and 9 × 10 −7 , the corresponding population numbers of infected and isolated individuals are described in Fig.  4 , respectively. Obviously, smaller contact rate can be conducive to the reduction in infection even lead to elimination of disease. For instance, by Fig. 4 , we can take some reasonable measures to guarantee the result of β < 3 × 10 −7 to eliminate the disease. Then we can compute

According to Theorem 4.1, the existence and uniqueness of the ergodic stationary distribution of system (2.2) is unknown. Figure 3 indicates that the disease will go to extinction in the long-term. Furthermore, when R 0 > 1, the deterministic system (2.1) has an endemic equilibrium E * which is globally asymptotically stable (Fig. 6 ).

In this subsection, we draw conclusion from the theoretical results of this paper.

• By means of Theorem 3.1, the existence of the ergodic stationary distribution (·) of the solution (S(t), I (t), Q(t), R(t)) to system (2.2) is proved under

• By taking the effect of stochasticity into account, the quasi-endemic equilibrium E * + related to E * is defined. Assuming that R S 0 > 1, we determine that the stationary distribution (·) around E * + admits a log-normal density function in the following form: where the special form of is shown in Theorem 4.1. Clearly, the above results indicate that the infectious disease will prevail and persist for long-term development if R S 0 > 1. In epidemiology, the first concern is whether an epidemic will occur. Therefore, by means of the above numerical simulations and parameter analyses, we are devoted to providing some effective measures to reduce the threat of infectious diseases to human life and safety, and eventually lead to the eradication of disease. Based on the above numerical simulations, we can conclude the following three points:

(1) To provide effective treatment and a wide range of isolation measures, see Fig. 5 .

(2) To control the activities of the susceptible individuals in highly pathogenic areas and provide some effective vaccination strategies for susceptible individuals, that is to say, β → 0 + , see Fig. 4 .

(3) To implement some reasonable policies to reduce population mobility in differential risk epidemic regions, which means the value of Λ is sufficiently small, see Fig. 3 . For example, various joint prevention and control greatly inhibited the spread of COVID-19 in China.

In this paper, to the best of our knowledge, the disease persistence of an SIQRS epidemic model, which includes the existence of ergodic stationary distribution and the exact expression of probability density function, is studied. Through comparison with the existing results, our main contributions will be introduced in the following two aspects in detail.

(i) By means of the linear random disturbance shown in [15, [20] [21] [22] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] , we focus on a stochastic SIQRS epidemic model with temporary immunity in the present study. Next, we construct some suitable Lyapunov functions to derive a stochastic critical value R S 0 . By the Khas'minskii ergodicity theory, we obtain that system (2.2) admits a unique ergodic stationary distribution while R S 0 > 1. From the similar expressions of R S 0 and the basic reproduction number R 0 , it greatly reveals that the stochastic positive equilibrium state (i.e., disease persistence) is only determined by the dynamical behavior of the susceptible individuals and the infectious people. More precisely, the corresponding random fluctuations are σ 1 , σ 2 . In view of the method of controlling variables and numerical simulations, the key measure to the prevention of infectious disease lies in the control and quarantine of infected individuals. Furthermore, by the corresponding parameter analyses, we further provide some reasonable measures to reduce the transmission of epidemic.

(ii) The existence of the ergodic stationary distribution makes it difficult to determine the exact statistical properties of disease persistence. For further dynamical investigation in epidemiology, based on Zhou and Zhang [33] , we develop some solving theories of algebraic equations with respect to the four-dimensional probability density function, which are described in Lemmas 2.3-2.5. In fact, focusing on the previous studies [27] [28] [29] [30] [31] [32] , the corresponding persistence is only obtained by the existence theory of the unique stationary distribution with ergodicity. By taking the effect of stochasticity into consideration, the quasi-endemic equilibrium E * + is constructed. By means of the equivalence of system (4.1) and the corresponding linearized system (4.3), we derive the exact expression of the log-normal four-dimensional density function Φ(S, I, Q, R). In addition, the covariance matrix is solved by the algebraic equation G 2 + A + A τ = 0, that is, Eq. (4.6). Following the existing results, the corresponding stability theory of zero solution of the general linear equation, described in [42] , can validate the positive definiteness of . But the specific form of is hard to obtain. In the current study, we develop the corresponding standard R 1 , R 2 , R 3 matrices shown in Lemmas 2.3-2.5. By means of the general solving theories, we can verify that is positive definite and obtain the special form of as shown in the detailed discussions. Furthermore, compared to what the existing results cannot obtain the general expression of , it is important to highlight that our methods and theories can be used to prove that is positive definite even if the diffusion matrix G is semi-positive definite, such as in delay stochastic differential equations [29, [43] [44] [45] .

Finally, some important topics that should be further studied are noted here. First, due to the limitation of our mathematical approaches to an epidemic model with temporary immunity, the sufficient conditions for disease extinction are difficult to establish. Consequently, for a comprehensive discussion, we only plot the relevant simulation of the solution (S(t), I (t), Q(t), R(t)) while R S 0 ≤ 1. Second, by taking the effect of telegraph noises into account [31, 35, 46] , the corresponding SIQRS epidemic model with temporary immunity and regime switching should be studied. These problems are expected to be considered and solved in our future work. where

, θ 13 = −θ 22 , θ 11 = b 2 θ 22 , where ϑ 11 = 1 2c 1 , ϑ 22 = 1 2c 1 c 2 . If c 1 > 0 and c 2 > 0, then θ 2 is a semi-positive definite matrix. This completes the proof.

By means of the invertible linear transformations, we will derive the corresponding standardized transformation matrices of standard R 1 , R 2 , and R 3 matrices. Obviously, we obtain the corresponding standard R 1 matrix M AM −1 := A 0 , which refers to (2.3). Let ρ 1 = a 21 a 32 a 43 σ and θ 0 = ρ −2 1 M M τ . Then the above equation can be equivalently transformed into the following equation: 

The transmission dynamic of different hepatitis B-infected individuals with the effect of hospitalization

Impact of drainage networks on cholera outbreaks in Lusaka, Zambia

A discrete model of avian influenza with seasonal reproduction and transmission

A contribution to the mathematical theory of epidemics

SVIR epidemic models with vaccination strategies

Stability and bifurcation analysis of an SIR epidemic model with logistic growth and saturation treatment

Analysis and modeling of tuberculosis transmission dynamics

Global dynamics of an SEIR epidemic model with vertical transmission

The dynamics of an HIV/AIDS model with screened disease carriers

Avian-human influenza epidemic model

Analysis of an HIV/AIDS treatment model with a nonlinear incidence

A generalization of the Kermack-McKendrick deterministic epidemic model

The first six weekssetting up a UK urgent dental care centre during the COVID-19 pandemic

Global stability of an SIRS epidemic model with transport-related infection

Effect of quarantine in six endemic models for infectious diseases

Global dynamics of an SIQR model with vaccination and elimination hybrid strategies

Global dynamics of an SIQR epidemic model with saturated incidence rate

Recurrent outbreaks of childhood diseases revisited: the impact of isolation

Homoclinic bifurcation in an SIQR model for childhood diseases

Dynamics of the deterministic and stochastic SIQS epidemic model with nonlinear incidence

Modeling and Dynamic of Infectious Disease

Cholera models with hyperinfectivity and temporary immunity

Sufficient and necessary conditions of stochastic permanence and extinction for stochastic logistic populations under regime switching

Analysis of a delayed vaccinated SIR epidemic model with temporary immunity and Lévy jumps

A stochastic epidemic model incorporating media coverage

The threshold of a stochastic SIS epidemic model with vaccination

The extinction and persistence of the stochastic hepatitis B epidemic model

Stationary distribution and extinction of a stochastic staged progression AIDS model with staged treatment and second-order perturbation

Global dynamics of a stochastic avian-human influenza epidemic model with logistic growth for avian population

Analysis of a stochastic distributed delay epidemic model with relapse and Gamma distribution kernel

Stationary distribution of an HIV model with general nonlinear incidence rate and stochastic perturbations

Nontrivial periodic solution for a stochastic brucellosis model with application to Xinjiang, China

Dynamical behavior of a stochastic epidemic model for cholera

Dynamics and density function analysis of a stochastic SVI epidemic model with half saturated incidence rate

The impact of virus carrier screening and actively seeking treatment on dynamical behavior of a stochastic HIV/AIDS infection model

Stationary distribution of stochastic SIS epidemic model with vaccination under regime switching

Stochastic Differential Equations and Applications

Stationary distribution and extinction of a stochastic SEIR epidemic model with standard incidence

Handbook of Stochastic Methods for Physics. Chemistry and the Natural Sciences

An asymptotic solution to a two-dimensional exit problem arising in population dynamics

An algorithmic introduction to numerical simulation of stochastic differential equations

Qualitative and Stability Methods for Ordinary Differential Equations

Long-time behaviour of a stochastic chemostat model with distributed delay

Global stability of multigroup epidemic models with distributed delays

Asymptotic behavior of stochastic multi-group epidemic models with distributed delays

A stochastic SIRS epidemic model with logistic growth and general nonlinear incidence rate

Conflict of interest The authors declare that they have no conflict of interest.

(I)Proof of Lemma 2.3: Consider the algebraic equationwhere θ 0 is a symmetric matrix. By direct calculation, we have

Consider the following k-dimensional stochastic differential equationdepicts a k-dimensional standard Brownian motion defined in the above complete probability space. The common differential operator L is described byLet the operator L act on a functiondV (X (t), t) = L V (X (t), t)dt +V X (X (t), t)g(X (t), t)d B(t).