key: cord-0756720-iohclx2m
authors: Zhou, Baoquan; Jiang, Daqing; Dai, Yucong; Hayat, Tasawar; Alsaedi, Ahmed
title: Stationary distribution and probability density function of a stochastic SVIS epidemic model with standard incidence and vaccination strategies
date: 2020-12-24
journal: Chaos Solitons Fractals
DOI: 10.1016/j.chaos.2020.110601
sha: 9bf24d396245baa5722f9d708994a03c8b281bca
doc_id: 756720
cord_uid: iohclx2m

Considering the great effect of vaccination and the unpredictability of environmental variations in nature, a stochastic Susceptible-Vaccinated-Infected-Susceptible (SVIS) epidemic model with standard incidence and vaccination strategies is the focus of the present study. By constructing a series of appropriate Lyapunov functions, the sufficient criterion [Formula: see text] is obtained for the existence and uniqueness of the ergodic stationary distribution of the model. In epidemiology, the existence of a stationary distribution indicates that the disease will be persistent in a long term. By taking the stochasticity into account, a quasi-endemic equilibrium related to the endemic equilibrium of the deterministic system is defined. By means of the method developed in solving the general three-dimensional Fokker-Planck equation, the exact expression of the probability density function of the stochastic model around the quasi-endemic equilibrium is derived, which is the key aim of the present paper. In statistical significance, the explicit density function can reflect all dynamical properties of an epidemic system. Next, a simple result of disease extinction is obtained. In addition, several numerical simulations and parameter analyses are performed to illustrate the theoretical results. Finally, the corresponding results and conclusions are discussed at the end of the paper.

It is well established that many infectious diseases have a critical influence on global social economies and human health. More precisely, the detailed statistics reported by the World Health Organization (WHO) show that approximately one-third of all deaths worldwide are caused by various epidemics. Recently, the global outbreak of COVID-19 with high transmission has also increased awareness of the importance of preventing and controlling infectious diseases. In epidemiology, mathematical models have provided several effective approaches to describe the characteristics and spread of epidemics in the last hundred years. In 1927, by dividing the population into two clusters, which includes people susceptible to the disease and infected individuals, Kermack and McKendrick [1] initially proposed the classical susceptible-infected-susceptible (SIS) epidemic model and established the corresponding threshold theory. Since then, various realistic ordinary differential equations (ODEs) have been extended to analyze and control the transmission of diseases [2] [3] [4] [5] [6] [7] [8] . For instance, Hove-Musekwa and Nyabadza [4] developed a HIV/AIDS model with active screening of disease carriers and obtained the corresponding basic reproduction number. Considering the effect of vertical infection, Tuncer and Martcheva [6] formulated a hepatitis B model with acute infection and carriers.

With the accelerated development of science and technology, vaccination comprises a common precaution that reduces the infection rate and even immunizes against some contagious diseases, such as measles, cholera, and tuberculosis [9] . According to a 2005 WHO report, the eradication of smallpox has been considered the most spectacular success of routine vaccination. Thus, some basic epidemic models with vaccination strategies have been studied in the last several decades [9] [10] [11] [12] [13] . In [9] , Liu et al. obtained the global stability of equilibria and analyzed the effect of pulse vaccination. Gao et al. [11] proposed mixed vaccination strategies in the SIRS epidemic model with seasonal variability on infection. However, the vast majority of infection processes are caused by person-toperson contact. As Ma and Wang [2] described, the classical bilinear incidence rate is reasonably assumed by the simple massaction law. From Anderson and May [14] , Hethcote [15] , this law is a good approximation for some communicable diseases, such as dengue fever and avian influenza. However, studies (see, e.g., [16, 17] ) have shown that the underlying assumption of homogeneous mixing and homogeneous environmental for several sexually transmitted diseases, e.g., HIV/AIDS and syphilis, may be invalid. In addition, owing to the psychological effect, susceptible individuals may tend to reduce the number of contacts with the infected per unit time as the numbers of the infected individuals increase [18, 19] . As a result, the corresponding adequate incidence rate should be modified as a nonlinear form. More importantly, Anderson and May [20, 21] pointed out that various epidemic models with standard incidence are suitable for human beings and some gregarious animals.

Given the above, a SVIS epidemic system with standard incidence and vaccination is the focus of the present study.

The total population N(t ) is divided into three compartments, namely susceptible people S(t ) , infected individuals I(t ) , and vaccinees V (t) that are in the vaccination process at time t. Then, the corresponding deterministic SVIS epidemic model with standard incidence and vaccination strategies takes the form

where Λ denotes the recruitment rate of the susceptible, β is the effective contact rate, μ depicts the natural death rate of the population, α denotes the additional death rate due to the disease, ϑ is the vaccination rate of the susceptible, γ denotes the immunity loss coefficient of the vaccinated, and δ reflects the recovery rate of the infected. In epidemiology, these biological parameters are assumed to be positive. Following similar results described by Ma and Zhou [22] , the corresponding basic reproduction number of system (1.1) takes the form

(1.2)

By defining a positive invariant set D 0 = (S, V, I) | S ≥ 0 , V ≥ 0 , I ≥ 0 , S + V + I ≤ Λ μ , two possible equilibria and their dynamical properties are then given as follows.

• Assuming that R 0 ≤ 1 , the disease-free, E 0 = (S 0 , V 0 , I 0 ) = Λ(μ+ γ ) μ(μ+ γ + ϑ ) , Λϑ μ(μ+ γ + ϑ ) , 0 , are then globally asymptotically stable in D 0 .

• If R 0 > 1 , there is a unique endemic equilibrium . Moreover, E + is globally asymptotically stable in D 0 , but E 0 is unstable.

In fact, Truscott and Gilligan [23] pointed out that the spread of infection, travel of populations, and design of control strategies are critically perturbed by some environmental variations. Therefore, it is more reasonable to construct a corresponding stochas-tic model to reveal the epidemiological characteristics of infectious diseases by comparison with the deterministic model. Notably, there are various possible approaches to simulate the random effects from biological significance and mathematical perspective [24] . For instance, making use of the fatal properties and multiplex networks, Zhu et al. [25] , Jia et al. [26] studied the SIR epidemic spreading process and analyzed individual decisionmaking behavior. In 2002, the most classical assumption that random changes always fluctuate around some average values due to continuous disturbances in nature, adopted by Mao et al. [27] , became a common way of describing environmental variations. Moreover, the above random fluctuations are all assumed to be types of white noise. Therefore, many authors have formulated the relevant stochastic differential equations (SDEs) with linear noises for the transmission of various epidemics [28] [29] [30] [31] [32] [33] [34] [35] [36] . As an example, Qi and Jiang [29] studied the impact of virus carrier screening and actively seeking treatment on the dynamical behavior of a stochastic HIV/AIDS epidemic model with bilinear incidence. In [34] , Shi and Zhang focused on the corresponding stochastic avian influenza system and investigated the existence of the unique ergodic stationary distribution. In addition, several dynamical analyses of the stochastic SIS models or epidemic systems with vaccination have been conducted [37] [38] [39] [40] . In [37] , Zhao and Jiang creatively proposed a general theory about extinction and persistence in mean based on a stochastic SIS epidemic model with vaccination. Zhang and Jiang [39] obtained sufficient conditions for a stochastic SIS system with saturated incidence and double epidemic diseases. By taking the periodicity effect into account, they still investigated a stochastic SVIR epidemic model with vaccination strategies, and derived the criteria for the existence of non-trivial positive periodic solution [40] . Given the above, in the present study it is assumed that the environmental noises are separately proportional to the compartments S, V and I. Then, the corresponding system (1.1) with the stochastic perturbations is described by

where B 1 (t) , B 2 (t) and B 3 (t) are three independent standard Brownian motions (or Wiener processes), with σ 2 i > 0 (i = 1 , 2 , 3) denoting their intensities. From the perspective of biomathematics, the existence and ergodicity of stationary distribution indicates that an infectious disease will prevail and persist in long-term development. More importantly, the corresponding probability density function of the stationary distribution can reflect all statistical properties of the individuals S, V and I. It can be regarded as a great intersection of epidemiological dynamics and statistics. It should be pointed out that there are relatively few studies devoted to the explicit expression of probability density function due to the difficulty of solving the corresponding Fokker-Planck equation. To the best of our ability, several general methods of solving the corresponding algebraic equations are developed herein that are equivalent to the Fokker-Planck equation, and the exact expression of density function is derived.

The rest of this paper is organized follows. In Section 2 , several mathematical notations and important lemmas for the dynamical analyses of system (1.3) are presented. The sufficient conditions for the existence and uniqueness of the ergodic stationary distribution of system (1.3) are obtained in Section 3 . By means of the developed approaches in solving the general three-dimensional Fokker-Planck equation, the exact expression of the probability density function for the stationary distribution is derived in Section 4 . In Section 5 , several simple criteria for the disease extinction of system (1.3) are given. In Section 6 , several numerical simulations are performed, together with parameter analyses to validate the theoretical results. Finally, the corresponding result discussions and main conclusions are shown, compared with existing articles, at the end of the paper.

satisfying the usual conditions (i.e., it is increasing and right continuous, while Γ 0 contains all P -null sets. Assuming that A m ×n is a real matrix, let A τ be the transpose matrix of A . If m = n, A −1 depicts the inverse matrix of A . The reader is referred to Mao [41] for detailed descriptions. For convenience, let R n be the n -dimensional Euclidean space, and R n

Clearly, the values S, V and I that satisfy system (1.3) are required to be positive for the corresponding dynamical behavior. To this end, the existence of uniqueness of the global positive solution of system (1.3) is described by the following Lemma 2.1 .

and the solution will remain in R 3 + with probability 1 (a.s.).

The detailed proof is almost the same as those in Zhou et al. [28] , and thus it is omitted here. Next, let X (t) be a homogeneous Markov process defined on R n that satisfies the following SDE, 

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

By Zhou et al. [28] , two important lemmas of solving the special algebraic equations are given as follows. 

If p 1 > 0 , p 3 > 0 and p 1 p 2 − p 3 > 0 , then ϒ 1 is positive definite, which follows [28] ) Consider the three-dimensional algebraic equation G 2

Assuming that q 1 > 0 and q 2 > 0 , then ϒ 2 is positive semi-definite, which takes the form

Finally, combining Lemmas 2.3 -2.4 and the Routh-Hurwitz criterion [43] , two general theories are developed in solving the similar algebraic equations, i.e., Lemmas 2.5 and 2.6 . 

. Assume that 0 is a symmetric matrix, for the three-dimensional algebraic equation

By defining the characteristic polynomials of A as ψ A (λ) = λ 3 + r 1 λ 2 + r 2 λ + r 3 , if A has all negative real part eigenvalues -that is,

Remark 1. From Zhou et al. [28] , A 0 and B 0 are called standard R 1 and R 2 matrices, respectively. Similarly, it is assumed that C 0 is a standard R 3 matrix. In addition, subsection (I) of Appendix A gives the detailed proof of Lemma 2.5 . The corresponding proof of Lemma 2.6 and the special form of 0 are shown in subsection (II) of Appendix A.

In this section, the focus is on the sufficient conditions for the existence and ergodicity of stationary distribution for system (1.3) . Moreover, one must guarantee that the results have no difference from those in the deterministic system (1.1) . Define

(3.1) Theorem 3.1. Assuming that R s 0 > 1 , then the solution (S(t) , V (t) , I(t )) of system (1.3) with any initial value (S(0) , V (0) , I(0)) ∈ R 3 + is ergodic and has a unique stationary distribution (·) .

Proof. By Lemma 2.1 , for any initial value (S(0) , V (0) , I(0)) ∈ R 3 + , system (1.3) has a unique global positive solution (S(t) ,

Then, the proof of Theorem 3.1 is divided into the following two steps: (i) construct a pair of a C 2 -Lyapunov function U (S, V , I) and bounded domain D such that L U ≤ −1 for any (S, V, I) ∈ R 3 + \ D , and (ii) validate the condition (A 1 ) of Lemma 2.2 .

Step 1. Consider a suitable C 2 -function W (S, V, I) in the form

where a 1 , a 2 , a 3 are all positive constant and are determined in (3.4) , and M 0 = Λ+2 μ+ γ + ϑ+

Define, for simplicity,

Applying It ˆ o 's formula to W 1 , which is shown in subsection (III) of Appendix B, obtains

Choosing a 1 , a 2 and a 3 such that , one can obtain

Additionally, note that W (S, V, I) is a continuous function satisfying

Hence, W (S, V, I) has a minimum value W 0 . Defining a non-negative C 2 -function U (S, V, I) : R 3

and combining (3.2) and (3.5) -(3.7) , it can be shown that

(3.8)

Next, the corresponding compact subset D is contructed by

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

(3.10)

For convenience, consider the following four subsets of R 3

Case 3. Assuming that (S, V, I) ∈ D 3 , , by (3.8) -(3.9) it can be seen that 

Notably, R 3

Hence, one equivalently obtains

Therefore, the condition (A 2 ) of Lemma 2.2 holds. The corresponding diffusion matrix is given by

Clearly, H is a positive-definite matrix. Then, the assumption (A 1 ) of Lemma 2.2 also holds. Given the above, the global positive solution (S(t) , V (t) , I(t )) of system (1.3) follows a unique ergodic stationary distribution (·) . This completes the proof of Theorem 3.1 .

Under R s 0 > 1 , the existence of the ergodic stationary distribution for system (1.3) is derived. This reveals that the contagious disease will prevail and persist in a population. Furthermore, from the expressions of R s 0 and R 0 , it can be obtained that R s 0 ≤ R 0 , and the sign holds if and only if σ 1 = σ 2 = σ 3 = 0 . Consequently, not only does this reveal that random fluctuations have a critical effect on the spread of epidemic, but it also indicates that R s 0 is a unified threshold of the disease persistence of systems (1.1) and (1.3) .

By Theorem 3.1 , one obtains that the global solution (S(t) , V (t) , I(t )) of system (1.3) follows a unique ergodic stationary distribution (·) . This section is devoted to deriving the explicit expression of the density function of the distribution (·) while R s 0 > 1 . In fact, the result will present a wide range of possibilities for the further development of epidemiological dynamics. Before this, two necessary transformations of system (1.3) should be first introduced.

(4.1)

By taking the random effect into consideration, another critical value is defined:

+ is determined by the following Eq. (4.2) :

3 ) ; thus, the corresponding linearized system of (4.1) takes the form

(·) Theorem 4.1. For any initial value (S(0) , V (0) , I(0)) ∈ R 3 + , if R s 0 > 1 ,

(·) around E * follows a unique lognormal probability density function Φ(S, V, I) , which is given by

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

(1) . If m 1 = 0 , m 2 = 0 and a 13 = 0 , then

(2) . If m 1 = 0 , m 2 = 0 and a 13 = 0 , then

(3) . If m 1 = 0 , m 2 = 0 and a 13 = 0 , then

(4) . If m 1 = 0 and m 2 = a 13 = 0 , then

Proof. For convenience and simplicity, let

−a 21 0 a 32 + a 33 −a 32 −a 33 .

Hence, system (4.3) can be rewritten as d X = AXd t + Md B (t) . By the theory of Gardiner [44] , the unique density function Φ(X ) around the quasi-endemic equilibrium E * satisfies the following Fokker-Plauck equation:

(4.5)

Since the diffusion matrix M is a constant matrix, Roozen [45] pointed out that Φ(X ) can be described by a quasi-Gaussian distribution, i.e., Φ(X ) = c 0 e − 1 2 XQX τ , where c 0 > 0 is determined by the normalized condition R 3 + Φ(X ) dX = 1 and Q is a symmetric matrix.

Substituting these results into (4.5) , one can obtain that Q obeys the algebraic equation QM 2 Q + A τ Q + QA = 0 . If Q is a inverse matrix, by letting = Q −1 , an equivalent equation is given by

Next, it will be proved that A has all negative real-part eigenvalues. The characteristic polynomial of A is defined as ψ (λ) = λ 3 + p 1 λ 2 + p 2 λ + p 3 , where r 1 = a 11 + a 21 + a 33 , r 2 = a 21 (a 11 − a 12 + a 33 ) + [ a 11 a 33 − a 13 (a 32 + a 33 )] , r 3 = a 21 a 33 (a 11 − a 12 − a 13 ) .

By the expressions of S * , V * , I * and N * , it can be shown that

Consequently, it follows from (i)-(iv) that (1) . r 1 = a 11 + a 21 + a 33 > 0 , r 3 = a 21 a 33 (a 11 − a 12 − a 13 ) > 0 , (2) . r 2 = a 21 (a 11 − a 12 + a 33 ) + [ a 11 a 33 − a 13 (a 32 + a 33 )] > (a 12 + a 13 ) a 33 − a 13 (a 32 + a 33 )

= a 12 a 33 − a 13 a 32 > 0 .

(3) . r 1 r 2 − r 3 = (a 11 + a 21 + a 33 ) { a 21 (a 11 − a 12 + a 33 ) + [ a 11 a 33 − a 13 (a 32 + a 33 )] } − a 21 a 33 (a 11 − a 12 − a 13 ) = a 11 r 2 + a 21 [(a 11 + a 33 ) a 33 − a 13 a 32 ] + a 33 [ a 11 a 33 − a 13 (a 32 + a 33 )] = a 11 r 2 + a 21 a 33 (a 11 − a 12 + a 33 ) + a 2 33 (a 11 − a 12 − a 13 ) + (a 21 + a 33 )(a 12 a 33 − a 13 a 32 ) > a 11 r 2 > 0 .

Combining the above (1)-(3) and Lemma 2.6 , that of Eq. (4.6) is positive definite can be derived. However, following the corresponding proof of Lemma 2.6 , which is shown in subsection (II) of Appendix A, the exact expression of is given. First, by the finite independent superposition principle, (4.6) can be equivalently transformed into the sum of solution to the following algebraic sub-equations, 0 , σ 3 ) , and the symmetric matrices k (k = 1 , 2 , 3) are their respective solutions. Clearly, = 1 + 2 + 3 . Now, the special expression of are derived by the following three steps. Step 1. For the algebraic equation

where the elimination matrix J 1 and A 1 are derived by

. By the value of w 1 , the relevant discussion is divided into two subcases:

Case (i) . If m 1 = 0 , in view of the method introduced in Zhou et al. [28] , it is assumed that

By direct calculation, one obtains

where r 1 , r 2 and r 3 are the same as above. Furthermore, one can equivalently transform Eq. (4.7) into

Noting that A has all negative real-part eigenvalues, then B 1 is a standard R 1 matrix. By Lemma 2.3 , this means that 0 is positive definite, which takes the form

(4.9) Therefore, 1 = 2 

where w 1 = a 11 + a 21 , w 2 = a 21 (a 11 − a 12 − a 13 ) , and ξ 1 is abbreviation. Obviously, B 1 is a standard R 2 matrix. Additionally, (4.7) can be equivalently transformed into

By letting Θ 1 = (a 21 σ 1 ) −2 ( H 1 J 1 ) 1 ( H 1 J 1 ) τ , it can be simplified as

In view of Lemma 2.4 , Θ 1 is described by

(4.11)

Step 2. Consider the algebraic equation In fact, one still derives

which means that B 2 is also a standard R 1 matrix. By letting Θ 2 = −2 2 (H 2 J 2 ) 2 (H 2 J 2 ) τ , where 2 = −a 32 m 2 σ 2 , (4.12) is then equivalent to the following equation:

In other words, 2 , and ξ 2 is also shorthand. Similarly, B 2 is a standard R 2 matrix. By defining Θ 3 = (a 32 σ 2 ) −2 ( H 2 J 2 ) 2 ( H 2 J 2 ) τ , (4.12) is then equivalent to

According to Lemma 2.4 , Θ 1 takes the form

(4.14)

Step 3. For the following algebraic equation, Hence, (4.15) can be equivalently transformed into .

Hence, an equivalent algebraic equation of (4.16) is described as follows:

where 3 = a 13 a 21 σ 3 , the last equation can be also simplified as

Similarly, one obtains 

Based on a 33 > 0 , then Θ 4 is a positive semi-definite matrix, which means that 3 

In summary, the special form of is divided into eight cases by the different values of m 1 , m 2 and a 13 , which is shown in Theorem 4.1 . Finally, in view of the relation of systems (4.1) and (4.3) , the stationary distribution (·) around E * then has a unique log-normal probability density function

Therefore, this completes the proof.

If R s 0 > 1 , Theorem 4.1 shows that the stationary distribution (·) around E * has the unique log-normal density function Φ(S, V, I) . This reflects the stochastic permanence of system (1.3) from one side. In addition, that R s

As is known, all of the properties of disease persistence of system (1. lim sup

which means that the epidemic of system (1.3) will go to extinction with probability 1 (a.s.).

o 's formula to ln I(t ) , one obtains

Integrating from 0 to t and dividing by t on both sides of (5.1) , it can be seen that

Next, by the strong law of large numbers [1] , one derives

Taking the superior limit of t → + ∞ on both sides of (5.3) , the assertion (5.1) can then be obtained by (5.4) . Moreover, from the expressions of R s 0 and R s 0 , one can obtain that R d 0 ≤ R s 0 . Consequently, the proof of Theorem 5.1 is confirmed. 

In this section, by means of the well-known higher-order method developed by Milstein [46] , the corresponding discretization equation

where the time increment t > 0 , and ξ k , η k , andζ k are three independent Gaussian random variables that follow the distribution N(0 , 1) for k = 1 , 2 , ..., n . Furthermore, (S k , V k , I k ) is the corresponding value of the k th iteration of the discretization equation. From AI-Darabsah [13] , Zhao and Jiang [37] , Liu et al. [38] , Zhang and Jiang [39] , Arino et al. [47] , and the detailed data of the Central Statistical Office of Zimbabwe (CSZ), the corresponding realistic statistics of system (1.3) are shown in Table 1 

It follows from Theorem 3.1 that system (1.3) admits a unique ergodic stationary distribution (·) . The left-hand column of Fig. 1 can be seen to validate it. By Theorem 4.1 , the stationary distribution (·) around the quasi-endemic equilibrium E * has a unique log-normal density function Φ(S, V, I) . Moreover, it is calculated that The curves of (i)-(iii) are shown in the right-hand column of Fig. 1 . Obviously, this greatly illustrates Theorem 4.1 from the side.

Combining Remarks 3.1 -4.1 and Theorem 5.1 , one can derive that all random perturbations σ 1 , σ 2 , and σ 3 have a critical influence on the dynamical behavior of system (1.3) . Therefore, the corresponding parameter analyses of the above three white noises are shown by Example 6.2 . First, it should be pointed out that the above four subcases (i) -(iv) all guarantee the existence of a stationary distribution, which has an ergodicity property. For convenience and simplicity, only the population intensities of susceptible and infected individuals are focused on, which are presented in subfigures (2-1) and (2-2) of Fig. 2 , respectively. By only increasing the perturbation intensities of the vaccinated individuals (or infected individuals), i.e., the larger σ 2 (or σ 3 ), then the disease infection will be effectively inhibited. In contrast, by only increasing the perturbation intensity of the susceptible individuals, a great destabilizing influence on the population numbers of S and I manifests. Next, by Zhu et al. [25] , Jia et al. [26] , the impact of the main parameters of system (1.3) on the individual decision-making behavior is studied. From the expressions of R s 0 and R d 0 , the disease persistence and extinction of system (1.3) are critically affected by the transmission rate β and vaccination rate ϑ. Thus, the following Examples 6.3 and 6.4 will reveal the impact. In addition, the corresponding effect of the recruitment rate Λ on the dynamical behavior of system (1.3) is also shown in Example 6.5 . Fig. 4 . Similarly, a small vaccination rate can control the disease infection more effectively than a large one. 

By Theorem 3.1 , one cannot derive the existence of the ergodic stationary distribution of system (1.3) . In contrast, it follows from Theorem 5.1 that the disease of stochastic system (1.3) will be extinct in a long term. In addition, the deterministic model (1.1) has a globally asymptotically stable endemic equilibrium E + . On the one hand, these results validate the fact that large white noises lead to disease elimination from the side. On the other hand, the large random fluctuation σ 3 (i.e., σ 3

= 78 >> 1) indicates that it is necessary to isolate and control the infected individuals during the outbreak of an epidemic. These results are verified by Fig. 6 .

For epidemiological study, combining the above numerical simulations and parameter analyses, several reasonable and effective measures to reduce the threat of infectious diseases to human life, and even eliminate the epidemic, are provided. The special approaches are the following.

(i) Several reasonable policies of joint prevention and control are implemented to reduce the population mobility in differential risk epidemic areas. Then, the small recruitment rate Λ may lead to the elimination of disease (see Fig. 5 ).

(ii) Controlling the activities of the susceptible individuals in highly pathogenic areas to decrease the contact rate of population. Hence, β → 0 + can be guaranteed, which means R s 0 < 1 and R d 0 < 1 (see Fig. 3 ). (iii) Developing several effective vaccines and carrying out other prophylactic measures to improve the immune rate of disease (see Fig. 4 ).

The corresponding theoretical results of this paper are the following. (i) By Theorem 3.1 , system (1.3) admits a unique ergodic stationary distribution

(ii) By taking the effect of random perturbations into account, a quasi-endemic equilibrium E * related to E + is defined while R c 0 = β(μ+ γ + 

In this paper, combining the great effect of vaccination and the unpredictability of environmental fluctuations in the real world, a stochastic SVIS infectious disease model with vaccination and standard incidence is the object of concern. Adopting the descriptions in [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] , linear perturbation, which is the most intuitive assumption of a random effect, is similarly taken into consideration in this paper. Subsequently, several dynamical properties of stochastic system (1.3) are analyzed, such as the existence and uniqueness of a global positive solution, existence and ergodicity of a stationary distribution, and disease elimination. By comparison with the existing results ( [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] ), several highlights developed in the present study are detailed in the following two points.

• As is known, the endemic equilibrium and basic reproduction number, two important results of a deterministic epidemic, reflect disease permanence and elimination. Similarly, for the corresponding stochastic model, the existence of stationary distribution indicates the stochastic positive equilibrium state. In this paper, it is first proved that stochastic system (1.3) admits a unique ergodic stationary distribution under the critical value R s 0 > 1 . It should be pointed out that R s 0 > 1 is a unified threshold for the disease persistence of systems (1.1) and (1.3) . Moreover, the sufficient condition R d 0 < 1 is obtained for the disease extinction of system (1.3) . Both R s 0 > 1 and R d 0 < 1 reveal that the dynamical behavior of system (1.3) is critically affected by the random fluctuations, i.e., σ 1 , σ 2 , and σ 3 . In view of the method of controlling variables and numerical simulations, this means that a large white noise leads to the disease eradication, while a small one guarantees stochastic permanence. In addition, by the main parameter analyses, several effective measures to stop the spread of an epidemic are provided.

• It is generally accepted that the existence of an ergodic stationary distribution incurs difficulty in studying more exact statistical properties. Hence, this paper is devoted to obtaining the corresponding probability density function for further dynamical investigation. The results of Zhou et al. [28] are further perfected and general solving theories of algebraic equations with respect to the three-dimensional Fokker-Planck equation are developed, which are described in Lemmas 2.5 and 2.6 . By taking the effect of stochasticity into account again, the quasi-endemic equilibrium E * corresponding to the endemic equilibrium E + is defined. For practical application, the exact expression of the log-normal three-dimensional density function Φ(S, V, I) of system (1.3) is given. Furthermore, it is worth mentioning that the methods and theories developed herein are still suitable for the case of the diffusion matrix M being positive semi-definite, such as delay stochastic differential equations [32, 48] .

Several remaining issues are now proposed and analyzed. First, by virtue of the limitation of the present mathematical approaches for epidemiological dynamics, a value gap exists between R s 0 and R d 0 , and it is unfortunate that difficulty is encountered in obtaining the most precise threshold for disease extinction and persistence. Second, the impact of telegraph noises and periodicity on the dynamics of system (1.3) should also be studied; the reader is referred to [30, 34, 45, 49, 50] . These problems are expected to be studied and solved as planned future work.

We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted. 

Obviously, one can obtain

Before proving the positive definiteness of 0 , the following two theories of matrix algebra should be described first.

(T 1 ) . The positive definiteness of the matrix is not affected by the inverse congruence transformation.

(T 2 ) . The similarity transformation does not change the eigenvalues of the matrix.

For convenience and simplicity, an important notation is introduced as follows. For the same dimensional square matrix A and B, define A B : A − B is at least a positive semi-definite matrix .

Given the above, it is easily derived that A is also positive definite if B is a positive definite matrix.

First, consider the following algebraic equation,

after which the relevant proof can be divided into the following two conditions:

(B 1 ) . a 21 = a 31 = 0 , (B 2 ) . a 21 = 0 or a 31 = 0 .

Next, one must demonstrate that the elements a 21 

it can be noticed that (i) . Π 1 = Π 1 , (ii) . 1 and 1 have the same positive definiteness .

In addition, a 21 = a 31 , a 31 = a 21 . Therefore, the validation is completed. Namely, one must only discuss the following two cases, which are equivalent to conditions (B 1 ) and (B 2 ) , respectively: Since A has all negative real part eigenvalues, by the similarity invariance of ψ A (λ) , it indicates that

Consequently, ϕ A (λ) has an eigenvalue λ 1 = a 11 , which has a negative real part. By a 11 ∈ R , then a 11 < 0 . In other words, Δ 11 is positive semi-definite. Moreover, 

(A.5) where the parameters ξ k (k = 1 , 2 , 3) can be obtained by (A.5) . Because their sign is the only object of concern, they are omitted here. Furthermore, the characteristic polynomial ψ A (λ) follows from A that ψ A (λ) = λ − a 33 + a 23 a 31 a 21

(λ 2 + ξ 1 λ + ξ 2 ) .

By the condition that A has all negative real part eigenvalues, it thus means that the equation λ 2 + ξ 1 λ + ξ 2 = 0 has two negative real part roots. By the Routh-Hurwitz stability criterion [43] , it can be shown that ξ 1 > 0 , ξ 2 > 0 . 

(A.7)

By means of ξ 1 > 0 and ξ 2 > 0 , it is derived that Δ 12 , L 2 and (F 3 F 2 ) −1 L 2 [(F 3 F 2 ) −1 ] τ are all positive semi-definite. It is then implied that 1 Δ 12 .

(A.8) ⎞ ⎠ := L 3 + L 4 , one can therefore derive, by a similar method as that described in (A.7) , In addition, for the following two algebraic equations, (i) . Π 2 + A 2 + 2 A τ = 0 , (ii) . Π 3 + A 3 + 3 A τ = 0 , and letting Π 2 = F 5 Π 2 F τ 5 , Π 3 = F 6 Π 3 F τ 6 , 2 = F 5 2 F τ 5 , 3 = F 6 3 F τ 6 , andA 2 = F 5 AF −1 5 , A 3 = F 6 AF −1 6 , where Given the above definitions and discussions, 0 is a positive-definite matrix. This completes the proof.

: For the above complete probability space { , F , { F t } t≥0 , P } , it is assumed that B (t) is an n -dimensional standard Brownian motion defined on it. Consider the following n -dimensional SDE, dX (t) = f (X (t ) , t ) dt + g(X (t ) , t ) dB (t ) , for t ≥ t 0 , with the initial value X (t 0 ) = X 0 ∈ R n . A common differential operator L is given by

Letting the operator L act on a function V ∈ C 2 , 1 (R n × [ t 0 , + ∞ ] ; R 1 + ) , one has

) n ×n . If X (t) ∈ R n , one has dV (X (t ) , t ) = L V (X (t ) , t ) dt + V x (X (t ) , t ) g(X (t ) , t ) dB (t ) .

A contribution to the mathematical theory of epidemics

A discrete model of avian influenza with seasonal reproduction and transmission

Global behaviour of an SIR epidemic model with time delay

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

Transmission dynamics and control strategies assessment of avian influenza a (H5N6) in the philippines

Modeling seasonality in avian influenza H5N1

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

Global stability of an SEI epidemic model with general contact rate

SVIR epidemic models with vaccination strategies

Qualitative analyses of SIS epidemic model with vaccination and varying total population size

Mixed vaccination strategy in SIRS epidmeic model with seasonal variability on infection

A susceptible-infected epidemic model with voluntary vaccinations

Threshold dynamics of a time-delayed epidemic model for continuous imperfect-vaccine with a generalized nonmonotone incidence rate

Infectious diseases in humans: dynamics and control

The mathematics of infectious diseases

The ergodicity and extinction of stochastically perturbed SIR and SEIR epidemic models with saturated incidence

Viral dynamics of an HIV model with latent infection incorporating antiretroviral therapy

Global dynamics of avian influenza epidemic models with psychological effect

Global analysis of an epidemic model with nonmonotone incidence rate

Infectious diseases of human: dynamics and control

Population biology of infectious diseases

Modeling and dynamic of infectious disease

Response of a deterministic epidemiological system to a stochastically varying environment

Dynamics of positive solutions to SIR and SEIR epidemic models with saturated incidence rates

Investigation of epidemic spreading process on multiplex networks by incorporating fatal properties

Evolutionary dynamics drives role specialization in a community of players

Environmental Brownian noise suppresses explosions in population dynamics

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

The extinction and persistence of the stochastic hepatitis b epidemic model

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

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

A stochastic epidemic model incorporating media coverage

Dynamics of an avian influenza model with half-saturated incidence

The ergodicity and extinction of stochastically perturbed SIR and SEIR epidemic models with saturated incidence

A stochastic epidemic model incorporating media coverage

The threshold of a stochastic SIS epidemic model with vaccination

Threshold behavior in a stochastic delayed SIS epidemic model with vaccination and double diseases

Dynamics of a stochastic SIS model with double epidemic diseases driven by L é vy jumps

Dynamical behavior of a stochastic SVIR epidemic model with vaccination

Stochastic differential equations and applications

Stochastic stability of differential equations. The Netherlands;: Sijthoff Noordhoff

Qualitative and stability methods for ordinary differential equations

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

Global results for an epidemic model with vaccination that exhibits backward bifurcation

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

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

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

This work was supported by the National Natural Science Foundation of China (Grant No. 11871473 ) and Shandong Provincial Natural Science Foundation (Grant Nos. ZR2019MA010 and ZR2019MA006 ).

where ϒ 3 is a symmetric matrix. Letting ϒ 3 := (κ i j ) 3 ×3 , by direct calculation one has