STATIONARY DISTRIBUTION OF STOCHASTIC SIRS EPIDEMIC MODEL WITH STANDARD INCIDENCE

. We study stochastic versions of a deterministic SIRS(Susceptible, Infective, Recovered, Susceptible) epidemic model with standard incidence. We study the existence of a stationary distribution of stochastic system by the theory of integral Markov semigroup. We prove the distribution densities of the solutions can converge to an invariant density in L 1 . This shows the system is ergodic. The presented results are demonstrated by numerical simulations.

1. Introduction. Infectious diseases are the second cause of leading death worldwide, after heart disease. Therefore, it is imperative to know the dynamical behavior of such diseases and to forecast what may happen. The spread of communicable diseases is generally described mathematically by compartmental models. Most epidemiological models stimulate from the classical SIR compartmental model of Kermack and McKendrick [14] which assumes a constant population. The assumption that the population is constant or asymptotically constant is often a reasonable approximation when modeling epidemics where the disease spreads quickly in the population and dies out within a short time (influenza, SARS, etc.) and the disease rarely causes deaths (West Nile virus in human or livestock). However, for endemic diseases in communities with changing populations (malaria), or diseases with high mortality rate (HIV/AIDS in poor countries), it is reasonable to assume varying population size and standard incidence rate in the mathematical theory of epidemiology (see, for instance, [1,2,3,4,17,11]).
In the SIRS models, susceptible individuals may become infected by contact with infective individuals, and after some time, recovered individuals become susceptible again. Busenberg and van den Driessche [7] have proposed the continuous SIRS epidemic model in a population with varying size: Here S, I and R denote the total numbers of susceptible, infective and recovered (removed) individuals respectively, and N is the total population size, The parameter b denotes the per capita birth rate, µ denotes the per capita disease free death rate, α denotes the disease-related per capita death rate of infected individuals. δ denotes the per capita loss of immunity rate of recovered individuals, γ denotes the per capita recovery rate of infected individuals. The force of infection is βI N with β as the effective per capita contact rate of an infective individual. All parameter values are assumed to be nonnegative. They discuss the extinction and persistence of the epidemic according to the threshold As a matter of fact, population systems are often subject to environmental noise. May [21] pointed out the fact that due to environmental fluctuations, the birth rates, death rates, carrying capacity and other parameters involved with the model system exhibit random fluctuations to a greater or lesser extent. Realistic models of population dynamics must take into account both predictable and unpredictable changes in those factors. The standard technique of parameter perturbation has been used by many authors for building stochastic epidemic models [8,9,20,15,30,28,27].
In studying epidemic dynamical system, we are always interested in when the disease will die out or prevail. In the deterministic models, the second problem is solved by showing that the endemic equilibrium of corresponding model is a global attractor or is globally asymptotically stable. But, there may be no endemic equilibrium in stochastic systems. Therefore, it is necessary to study the existence of stationary distribution to the solution and whether the solution is ergodic. These properties can reveal the disease to persist. Mao [9] showed the existence of a unique stationary distribution of the solution of the stochastic SIS epidemic model based on the theory of Hasminskii [10]. Lahrouz and Omari [16] studied a stochastic SIRS epidemic model with general incidence rate in a population of varying size. Sufficient conditions for the extinction and the existence of a unique stationary distribution are obtained. Lin and Jiang [18] discussed a stochastic classic SIR system with bilinear incidence. They showed sufficient conditions for the disease to extinct exponentially. In the case of persistence they proved the existence of a stationary distribution and found the support of the invariant density.
Recently, Zhao et al. [29] considered the stochastic versions of model (1) with varying population and standard incidence: Here B t is a standard Brownian motion with B(0) = 0, the white noise intensity is σ 2 > 0. The authors [29] obtained that, a crucial threshold is determined which indicates the differences between the deterministic and stochastic model. When the threshold is less than one or the noise intensity is large, they deduce the disease to extinct exponentially. When the threshold is more than one, sufficient conditions for persistence in the mean are established. But in the case of persistence they can not obtain the existence of stationary distribution of system (3). The aim of this paper is to fill the gap. Hence our work can be considered as the further work of Zhao et al. [29]. From system (3), the equation for the total population size isṄ (t) = (b − µ)N t − αI t . Obviously, the total population is a variable. Thus it is convenient to work with proportions of susceptibles, infectives and recovereds in the population (as in deterministic model [7]). Defining it gives x t + y t + z t = 1. Then it is sufficient to study the stochastic differential equation of x t , y t and z t . Using Itô's formula and the relation z t = 1 − x t − y t , we can omit the equation of z t and discuss the following system with any given initial value (x 0 , y 0 ) ∈ R 2 + and x 0 + y 0 < 1. In this paper, we will study the long-time behavior of system (6). The main aim of this paper is to study the existence of a stationary distribution of system (6) and its asymptotic stability. To the best of our knowledge, there are no literatures analytically concerning the existence of a stationary distribution of epidemic models with standard incidence. This is mainly because the Fokker-Planck equation (7) corresponding to system (6) is of a degenerate type. Therefore, it is not within the scope of application of Has'minskii theorem which is used in [9,19,13]. The key to our analysis approach is based on the theory of integral Markov semigroups, which was used in [23,24,18]. The specific strategy is given section 3. In the appendix, we present some auxiliary results and the main tools concerning Markov semigroups.

2.
Preliminaries. Now we introduce a Markov semigroup connected with the Fokker-Planck equation (7). Let X = R 2 + , Σ be the σ-algebra of Borel subsets of X, and m be the Lebesgue measure on (X, Σ). We denote P(t, x, y, A) to the transition probability function for the diffusion process (x t , y t ) of (6) with the initial condition (x 0 , y 0 ) = (x, y), i.e. P(t, x, y, A) = Prob((x t , y t ) ∈ A). If, for t > 0, the distribution of (x t , y t ) is absolutely continuous with respect to the Lebesgue measure with the density u(t, x, y), then u(t, x, y) satisfies the Fokker-Planck equation: Since P (t) is a contraction on D, it can be extended to a contraction on L 1 (X, Σ, m).
The adjoint operator of L is of the form In [29], the existence and uniqueness of positive solution of system (1.6) is given.
There is a unique solution (x(t), y(t)) of system (6) on t ≥ 0 for any initial value (x(0), y(0)) ∈ Γ, and the solution will remain in Γ with probability 1, where 3. The stationary distribution of the solution. In this section, we investigate the existence for a stationary distribution of system (6).
Theorem 3.1. Let (x t , y t ) be a solution of system (6). Then for every t > 0, the distribution of (x t , y t ) has a density u(t, x, y) which satisfies Fokker-Planck equation (7). IfR 0 > 1, then there exists a unique density u * (x, y) which is a stationary solution of (7) and In addition, Remark 1. By the support of a measurable function f we simply mean the set . This is not the customary definition of the support of a function, which is usually the closure of the set supp f , but this slightly unusual definition is more suitable our purposes. By Theorem 2.1, for any initial value (x(0), y(0)) ∈ Γ, and the solution (x(t), y(t)) of system (6) on t ≥ 0 will remain in invariant set Γ with probability 1. Therefore we consider Γ is the whole space. As a result, the support of the invariant density u * is shown in (9). This results from the fact that the Fokker-Planck equation corresponding to system (6) is of a degenerate type. Our approach comes from Markov semigroup theory, which was used in [23,24,18], the specific strategy is as follows: First, using the Hörmander condition [6] we show that the transition function of the process (x t , y t ) is absolutely continuous (see Lemma 3.1). Then, using support theorems [26,5,25] we find a set G on which the density of the transition function is positive (G is given in (9)) (see Lemma 3.2 and Lemma 3.3). Next, we show that the Markov semigroup satisfies the "Foguel alternative" ( see Appendix), i.e. it is either asymptotically stable or "sweeping" (see Lemma 3.3). Finally, we exclude sweeping by showing that there exists a Khasminskiȋ function (24) (see Lemma 3.5). In this way we prove the most difficult part of the paper to show asymptotic stability of system (6).
In the following, we realize this strategy by Lemma 3.2-3.6.
Proof. If a(x) and b(x) are vector fields on R d , then the Lie bracket [a, b] is a vector field given by ] are linearly independent on Γ. This implies that at any point (ξ, η) ∈ Γ, linearly independent vectors b(ξ, η), [a, b](ξ, η) can span R 2 . Thus the vector fields a and b satisfy the Hörmander condition [6]. Then using Hörmander Theorem we show that the transition probability function P(t, x 0 , y 0 , A) has a density k(t, x, y; x 0 , y 0 ) and k ∈ C ∞ ((0, ∞) × Γ × Γ).
Remark 2. According to Lemma 3.2, there exists a measurable function k(t, x, y; x 0 , y 0 ), such that for every density f . By the definition (see Appendix A), the semigroup {P (t)} t≥0 is an integral Markov semigroup.
Next we prove that our claim holds on G 0 . There exists a positive constant T and a differentiable function then u ψ satisfies (14) for t ∈ [0, T ]. We split the construction of the function u ψ on three intervals [0, τ ], (τ, T − τ ) and [T − τ, T ], where 0 < τ < T /2. Let and Therefore, we can construct a C 2 -function u ψ : Taking T sufficiently large we can extend the function u ψ : and therefore, the function u ψ satisfies (16) on [0, T ] in view of (17), (18). It follows that we can find a C 1 -function x ψ which satisfies (11) and finally we can determine a continuous function ψ from (11) . Therefore, our claim holds. This completes our proof.
where G is given in (9).
Proof. Let U t = x t + y t . The system (6) can be replaced by where g 1 and g 2 are as in (12). Since (x t , y t ) is a positive solution of system (6) with probability 1, from the expression of g 2 , we get Similar to (14), it is rewritten as Now we claim that for almost every ω ∈ Ω there exists t 0 = t 0 (ω) such that which completes our proof. According to the position of initial value U 0 we consider two cases: (21) is false, then we know that there exists Ω ⊂ Ω with Prob(Ω ) > 0 such that U t (ω) ∈ (0, m 1 ], ω ∈ Ω . By (20), it follows that for any ω ∈ Ω , U t (ω) is strictly increasing on [0, +∞), and therefore lim t→∞ U t (ω) = m 1 , ω ∈ Ω . From the expression of g 2 , it follows that lim t→∞ x t (ω) = 0 and lim By Itô's Formula, we get which yields Let M (t) = t 0 x r dB r , by using Strong Law of Large Numbers (Lemma A.2), we obtain lim t→∞ Case 2. U 0 ∈ (m 1 , 1). From (20), it is obvious that our claim (21) holds and lim t→∞ U t = m 1 , lim t→∞ U t = 1. By similar arguments to Case 1, we obtain lim t→∞ U t = m 1 . If lim t→∞ U t = 1, there exists Ω ⊂ Ω with Prob(Ω ) > 0 such that for any ω ∈ Ω , lim t→∞ U t (ω) = 1, lim t→∞ x t (ω) = 1 and lim t→∞ y t (ω) = 0.
Remark 3. From lemma 3.3 and 3.4, we know that if (7) has a stationary solution u * , then suppu * = G. So according to Lemma A.1, The semigroup {P (t)} t≥0 is asymptotically stable or is sweeping with respect to compact sets.
Proof. We will construct a nonnegative C 2 -function V and a closed set U ∈ Σ such that sup Such a function is called a Khasminskiȋ function in [22]. The existence of a Khasminskiȋ function implies that the semigroup is not sweeping from the set U . Consider the function where r is a positive constant satisfied The function H(x, y) satisfies From the above equations, we can deduce y = rx. Then It is not difficult to know f (x) is a monotonic increasing function and f ( m1 r+1 ) = −∞, f ( 1 r+1 ) = +∞. So H(x, y) attains its minimum value at the only stable point (θ, rθ), and H(x, y) ≥ H(θ, rθ).
We construct a Lyapunov function V : G →R + by Applying the Itô's formula, we obtain by (14), Then (26) where ε, κ > 0 are two small numbers such that and Denote ε,κ . In the following, we consider four cases: Case 1. On D 1 ε,κ , we have Then Noting that ε satisfied (28), so we can obtain By (25) and (27), we have Case 3. On D 3 ε,κ , let κ = ε 2 , from (29), then Case 4. On D 4 ε,κ , in view of (30), we have In summary, on G \ U ε,κ = D 1 According to Lemma A.1, the semigroup {P (t)} t≥0 is asymptotically stable. In another words, the semigroup has a unique stationary solution on G. This completes the proof.

Remark 4.
In the proof of lemma 3.6, we take X = Γ. To verify V is a Khasminskiȋ function, it suffices that there exist a closed set U ⊆ Σ ( which lies entirely in Γ) such that sup x / ∈U L * V (x) < 0.

4.
Simulations. Next we make numerical simulations to illustrate our results by using Milstein's Higer Order Method [12]. We assume that the unit of time is one day and the population sizes are measured in units of 1 million. The parameters in (6) are given by b = 0.2, β = 0.6, α = 0.1, γ = 0.2, δ = 0.2, σ = 0.1, t = 1.
In this case, the condition of Theorem 3.1 is satisfied. We find that these lines in Figure 1 fit very well which implies that wherever x(t) and y(t) start from, the density functions of x(t) and y(t) converge to the same functions respectively. Figure 2 indicates that there is a stationary distribution for system (6). Hence, Figure 1 and Figure 2 approve the result of theorem 3.1.
Appendix A. . Since the proof of our result is based on the theory of integral Markov semigroups, we need some auxiliary definitions and results concerning Markov semigroups (see [23], [24]). For the convenience of the reader, we present these definitions and results in the appendix. Let the triple (X, Σ, m) be a σ-finite measure space.  for every density f . A family {P (t)} t≥0 of Markov operators which satisfies conditions: (1) P (0) = Id, (2) P (t + s) = P (t)P (s) for s, t ≥ 0, (3) for each f ∈ L 1 the function t → P (t)f is continuous with respect to the L 1 norm, is called a Markov semigroup. A Markov semigroup {P (t)} t≥0 is called integral, if for each t > 0, the operator P (t) is an integral Markov operator.
We also need two definitions concerning the asymptotic behaviour of a Markov semigroup. A density f * is called invariant if P (t)f * = f * for each t > 0. The Markov semigroup {P (t)} t≥0 is called asymptotically stable if there is an invariant density f * such that lim t→∞ P (t)f − f * = 0 for f ∈ D.  We need some result concerning asymptotic stability and sweeping which can be found in [23] (see Corollary 1).
Lemma A.1. Let X be a metric space and Σ be the σ-algebra of Borel sets. Let {P (t)} t≥0 be an integral Markov semigroup with a continuous kernel k(t, x, y) for t > 0, which satisfies (34) for all y ∈ X. We assume that for every f ∈ D we have ∞ 0 P (t)f dt > 0 a.e.