Dynamics of Stochastically Perturbed SIS Epidemic Model with Vaccination

We introduce stochasticity into an SIS epidemic model with vaccination. The stochasticity in the model is a standard technique in stochasticpopulationmodeling. Inthedeterministicmodels, thebasicreproductionnumber 𝑅 0 isathresholdwhichdeterminesthe persistence or extinction of the disease. When the perturbation and the disease-related death rate are small, we carry out a detailed analysis on the dynamical behavior of the stochastic model, also regarding of the value of 𝑅 0 . If 𝑅 0 ≤ 1 , the solution of the model is oscillating around a steady state, which is the disease-free equilibrium of the corresponding deterministic model, whereas, if 𝑅 0 > 1 , there is a stationary distribution, which means that the disease will prevail. The results are illustrated by computer simulations.


Introduction
Epidemiology is the study of the spread of diseases with the objective to trace factors that are responsible for or contribute to their occurrence. Significant progress has been made in the theory and application of epidemiology modeling by mathematical research. Controlling infectious diseases has been an increasingly complex issue in recent years. Vaccination is an important strategy for the elimination of infectious diseases [1][2][3]. The vaccination enables the vaccinated to acquire a permanent or temporary immunity. When the immunity is temporary, the immunity can be lost after a period of time. This temporary immunity is used in many references [4][5][6][7] that assume the process of losing immunity is in the exponential form.
In [7], Li and Ma discussed an SIS model with vaccination. The system has following form:  The parameters in the model are summarized in the following list: : a constant input of new members into the population per unit time; : a fraction of vaccinated for newborns; : transmission coefficient between compartments and ; : natural death rate of , , and compartments; : the proportional coefficient of vaccinated for the susceptible; : recovery rate of infectious individuals; : the rate of losing their immunity for vaccinated individuals; : disease-caused death rate of infectious individuals.

Abstract and Applied Analysis
All parameter values are assumed to be nonnegative and , > 0.
If 0 > 1, then the disease will be persistent, where 0 = 0 − ( 2 2 /2( + )) and 0 = /( + ) is the threshold of the corresponding deterministic system. 0 can be considered as the threshold of system (3), which is less than the value of 0 .
On the other hand, white noise stochastic perturbations around the positive endemic equilibrium of epidemic models were considered in [19,25]. Beretta et al. [25] considered the following system differential equation: where * , * , and * are the positive points of equilibrium for the corresponding deterministic system and are constants, ( ) are independent from other standard Wiener processes ( = 1, 2, 3). In their paper, they proved the above stochastic system was stable in probability by using the general method of Lyapunov functionals construction and obtained the stability condition immediately in terms of the parameters of the system under consideration.
Another different approach to include stochastic perturbations in a biological model was considered by Imhof and Walcher in [11]. They introduced and analyzed a variant of the deterministic single-substrate chemostat model. Also they found all species were persistent in the multiple species by using a comparison principle. Subsequently, they turned to modeling the influence of random fluctuations on the deterministic chemostat model. They set up the following stochastic chemostat model: where 0 , 1 , 2 ≥ 0 and ( )( = 1, 2, 3) are independent Brownian motions. By analyzing this stochastic differential equation model, they proved that the stochastic model led to extinction even though the deterministic counterpart predicts persistence. In this paper, our approach to include stochastic perturbation is analogous to that of Imhof and Walcher [11]. We are devoted to studying the following stochastic system: where ( )( = 1, 2, 3) are independent Brownian motions and ( = 1, 2, 3) are their intensities. This paper is organized as follows. In Section 2, we show there is a unique positive solution of system (6) by the way mentioned in [20,23]. In Section 3, when 0 ≤ 1, we derive that the solution of the system (6) oscillates around the disease-free equilibrium 0 of system (1). In Section 4, when

Existence and Uniqueness of Positive Solution
To investigate the dynamical behavior, the first concern is whether the solution has a global existence. Moreover, for a population dynamics model, whether the value is nonnegative is also considered. Hence in this section, we first show that the solution of system (6) is global and nonnegative.
Proof. Our proof is motivated by the works of Mao et al. [17]. Since the coefficients of (6) are locally Lipschitz continuous for any given initial value ( (0), (0), (0)) ∈ R 3 + , there is a unique local solution ( ( ), ( ), ( )) on ∈ [0, ), where is the explosion time (see [20]). To show that this solution is global, we need to show that = ∞ a.s. Let 0 ≥ 0 be sufficiently large so that (0), (0), and (0) all lie within the interval [1/ 0 , 0 ]. For each integer ⩾ 0 , define the stopping time: where throughout this paper, we set inf 0 = ∞ (as usual 0 denotes the empty set). According to the definition, is increasing as → ∞. Set ∞ = lim → ∞ , whence ∞ ⩽ a.s. If we can show that ∞ = ∞ a.s., then = ∞ and ( ( ), ( ), ( )) ∈ R 3 + a.s. for all ⩾ 0. In other words, to complete the proof, all we need to show is that ∞ = ∞ a.s. If this statement is false, then there exist a pair of constants > 0 and ∈ (0, 1) such that Hence there is an integer 1 ⩾ 0 such that Define a 2 -function : R 3 + → R (see the appendix) by The nonnegativity of this function can be seen from − 1 − log ⩾ 0, for all > 0. Let ⩾ 0 and > 0 be arbitrary. Applying the Itô formula, we obtain where The remainder of the proof follows that in Mao et al. [17].

Asymptotic Behavior around 0
It is clear that 0 = ( 0 , 0, 0 ) is the disease-free equilibrium of system (1). If 0 ≤ 1, then 0 is globally stable, which means the disease will die out after some period of time.
Hence, it is interesting to study the disease-free equilibrium for controlling infectious disease. But there is no disease-free equilibrium in system (6). In this section, we show the average oscillation around 0 in time to exhibit whether the disease will die out.
Remark 3. From Theorem 2, we show the solution of system (6) will oscillate around the disease-free equilibrium of system (1) under some conditions and the disturbance intensity is proportional to the intensity of the white noise. In a biological view, as the intensity of stochastic perturbations is small, we consider the disease will die out. Besides, if 1 = 0 and 3 = 0, then 0 is also the diseasefree equilibrium of system (6). From the proof of Theorem 2, we get Thus, the solution of system (6) is stochastically asymptotically stable in the large if ( + 2)(1 + 2 / ) − /4 > ( 2 + 1) 2 2 . (6) In studying epidemic dynamical system, we are also interested in when the disease will prevail and persist in a population. In the deterministic models, the problem is solved by showing that the endemic equilibrium of corresponding model is a global attractor or is globally asymptotically stable. But there is no endemic equilibrium in system (6). In this section, based on the theory of Has'minskiȋ [26] (see the appendix), we show that there is a stationary distribution subjected to some conditions on 0 and the parameters of the model which reveals that the disease will prevail also. System (6) can be written as a form of system (A.8) (see the appendix). Consider    [26], we obtain that the solution ( ) is a homogeneous Markov process in 3 + .

Numerical Simulations
In order to confirm the results above, we numerically simulate the solution of system (6) with given initial value and the parameters. Using Milstein's higher order method mentioned in [27], we get the discretization equation of system (6): where 1, , 2, , and 3, , = 1, 2, . . . , , are the independent Gaussian random variables (0, 1). Choosing suitable parameters in the system, by Matlab we get the simulation figures with initial value ( (0), (0), (0)) = (0.8, 0.4, 0.5) and time step Δ = 0.001. Example 6. Throughout the paper we will assume that the unit of time is one day, and the population sizes are measured in units of 1 million. Choose the parameters in system (6) and 1 , 2 , and 3 satisfy condition (14), then by Theorem 2, we show the solution of system (6) will oscillate around the disease-free equilibrium 0 ( 0 , 0 , 0 ) of system (1) in time which is globally asymptotically stable. Besides Theorem 2 tells us that the difference between the perturbed solution and 0 is only related with white noises 1 and 3 . Numerical simulations are shown in Figure 1.  (1) is globally asymptotically stable. But the white noise may make system (1) appear as different phenomena. In detail, the conditions (30) and (31) are also satisfied. Therefore, as Theorem 5 stated, there is a stationary distribution (see the histograms in the middle in Figure 2). In addition, the left pictures in Figure 2 show that the solution of system (6) has fluctuating in a small neighborhood. Moreover, from Figure 3, we find that 95% or more of the population distribution lies within a neighborhood, which can be imagined a circular or elliptic region centered at * ( * , * , * ) (see the red point in Figure 3). All of these imply system (6) has stochastic stability.
Testing these data for normality, all tests used were highly significant, conclusively rejecting normality in all cases. This is not surprising in view of the very large sample sizes (10,000), as even moderate deviations from the tested distribution will be significant; however the normal quantilequantile plots (see the QQ plots in the right in Figure 2) suggest that these data are not far from being normally distributed for smaller values of ( = 1, 2, 3).

Conclusions
As most real world problems are not deterministic, incorporating stochastic effects into the model give us a more realistic way of modeling epidemic models. In this paper, we have considered stochastic SIS epidemic models with vaccination. We first proved the positivity of the solutions. When the perturbation and the disease-related death rate are small, we illustrated the dynamical behavior of the SDE model according to 0 ≤ 1 or 0 > 1. When 0 ≤ 1, we proved that the asymptotic behavior around the disease-free equilibrium of the deterministic model. When 0 > 1, we also proved that the SDE model has the ergodic property as the fluctuation is small, where the positive solution converges weakly to the unique stationary distribution. Simulations are also carried out to verify our analytical results.

Appendix
Here, we give some basic theory in stochastic differential equations (see [18]).
Throughout this paper, unless otherwise specified, let (Ω, { } ≥0 , ) be a complete probability space with a filtration { } ≥0 satisfying the usual conditions (i.e., it is right continuous, and 0 contains all P-null sets). Denote  Next, there is some theory about stationary distributions (see Has'minskiȋ [26], 1980). Let ( ) be a homogeneous Markov process in ( denotes Euclidean -space) described by (B.2) If ∈ \ , the mean time at which a path issuing from reaches the set is finite and sup ∈ < ∞ for every compact subset ⊂ . Lemma A.2 (see [26]). If (B) holds, then the Markov process ( ) has a stationary distribution (⋅). Let (⋅) be a function integrable with respect to the measure . Then Remark A. 3. The proof is given in Has'minskiȋ [26]. The existence of a stationary distribution with a density is given in Theorem 4.1, page 119, and Lemma 9.4, page 138. The weak convergence and the ergodicity are obtained in Theorem 5.1, page 121, and Theorem 7.1, page 130.
Lemma A.4. Let ( ) be a regular temporally homogeneous Markov process in . If ( ) is recurrent relative to some bounded domain , then it is recurrent relative to any nonempty domain in .