Comparison of deterministic and stochastic SIRS epidemic model with saturating incidence and immigration

The purpose of this work is to compare the stochastic and deterministic versions of an SIRS epidemic model. The SIRS models studied here include constant inflows of new susceptibles, infectives and removeds. These models also incorporate saturation incidence rate and disease-related death. First, we study the global stability of deterministic model with and without the presence of a positive flow of infectives into the population. Next, we extend the deterministic model system to a stochastic differential system by incorporating white noise. We show there is a unique positive solution to the system, and the long-time behavior of solution is studied. Mainly, we show how the solution goes around the endemic equilibrium and the disease-free equilibrium of deterministic system under different conditions on the intensities of noises and the parameters of the model. Finally, we introduce some numerical simulation graphics to illustrate our main results.


Introduction
Mathematical modeling is an essential tool in studying the spread of infectious diseases. Understanding the transmission of an infectious disease in communities, regions, and countries is a crucial issue to prevent major outbreaks of an epidemic. Mathematical models have been used in planning, comparing, implementing, evaluating, and optimizing various detection, prevention, therapy, and control programs. Deterministic models for communicable diseases have been introduced in a systematic way by Kermack and McKendrick [9]. Since then, large literature has grown for both deterministic and stochastic models [1,2,4,7]. These models typically take the form of ordinary or stochastic differential system which split the total population into distinct classes. One line of investigation classifies individuals as one of susceptible, infectious and removed with permanent acquired immunity. Such model is termed as an SIR model. More realistically, if the removed individuals lose immunity and return to the susceptible compartment, it becomes an SIRS type [15].
In the mathematical study of epidemiological problems, the incidence rate that measures the rate of new infection is considered to be a very crucial parameter. It is assumed to be, in most classical disease transmission models, of mass action type with bilinear interactions given by β S I , where β is the per capita contact rate, and S and I represent the susceptible and infected populations, respectively. But, there are several reasons for using non-linear incidence rates. For instance, during the outbreak of SARS in 2003, many protection measures and control policies were taken by the Chinese government such as closing schools, closing restaurants, postponing conferences, isolating infectives, etc [21]. These actions greatly reduced the contact number per unit time. A number of non-linear incidence rate have suggested [5,10]. Therein, Liu et al. [10] discussed possible mechanisms that could lead to non-linear incidence and demonstrated that the SIRS model with incidence β S p I q can have periodic solutions for q > 1. To prevent the unboundedness of contact rate, Capasso and Serio [5] used a saturated incidence rate of the form β S I 1+ I , where β I measures the infection force of the disease and 1 1+ I measures the inhibition effect from the behavioral change of the susceptible individuals when their number increases or from the crowding effect of the infective individuals.
The rapid increases in international travel and trade and the mass movement of populations witnessed in the last few decades mean that infectious diseases can spread from one continent to another in a matter of hours or days, whether they are conveyed by individual travelers or in the cargo holds of aircraft or ships [20]. The role of immigration can have dramatic consequences on national public health programs. For instance, in the Great Toronto Area, approximately 95 % of tuberculosis cases are reported among foreign-born persons, who account for just one-half of the city's inhabitants. The proportion of TB cases among foreign-born individuals has been on the rise in Canada since the 1970s [22]. Epidemics ignited or enhanced by immigration of infectious cases include HIV, SARS, avian influenza and measles [18].
In this paper, we construct an SIRS model with saturated incidence rate that includes the immigration of distinct compartments, that is, there are susceptible, infective, and removed individuals in the new number of immigration. We assume a constant flow of new members into the population per unit time, of which a fraction p 1 ( p 1 ≥ 0) is infective, and a fraction p 2 ( p 2 ≥ 0) is removed, so the fraction 1− p 1 − p 2 is susceptible with 0 ≤ p 1 + p 2 ≤ 1. Thereby we study the SIRS model described by the following differential system: here, we assume that the susceptible, infective and removed individuals die at different rate respectively μ, μ+c 1 and μ+c 2 . α is the recovery rate of the infective individuals and γ is the rate at which removed individuals lose immunity and return to the susceptible class. For the bilinear incidence rate (when = 0) and p 2 = 0, Brauer and Van den Driessche [3] studied the global behavior of the SIR model (γ = 0) that includes immigration of infective individuals and variable population size. Therein, the stability of the unique endemic equilibrium state is proved by reducing the system to an integro-differential equation. In the special case where there is no input of infective and removed individuals, Mena-Lorca and Hethcote [15], constructed a Lyapunov function for an SIR epidemic model. Unfortunately, they did not extend the method to the SIRS model with temporary immunity.
Recently, by combining quadratic functions and the voltera function (I − I * ln I ), Lahrouz et al. [13] have constructed a Lyapunov function to prove the global stability of an SIRS model with saturated incidence rate. As a matter of fact, there are real benefits to be gained in using stochastic models because real life is full of randomness and stochasticity. Wherefore we will perturb the deterministic model (1) by a white noise. Recently, several authors studied stochastic epidemic systems. In [2,6,14], the situation of a white noise stochastic perturbations around the endemic equilibrium state was considered. Beretta et al. [2] proved, by constructing Lyapunov functions, the stability of an SIR model with delays influenced by stochastic perturbations under certain conditions on the intensity of the white noise. Carletti [6] investigated the stability properties of a stochastic model for phage-bacteria interaction in open marine environment both analytically and numerically. Lahrouz et al. [14] studied the mean-square and asymptotic stability in probability of a mathematical model of smoking. Dalal et al. [7] introduced stochasticity into a model of AIDS and condom use via the technique of parameter perturbation which is standard in stochastic population modeling. Lahrouz et al. [13] introduced noise into an SIRS model by perturbing the contact rate. Both of them show that the models established in their paper possess non-negative solutions. They also carried out a detailed analysis on asymptotic stability both in probability one and in pth moment. In this paper, we introduce randomness into the deterministic model (1) by replacing the parameters c 1 and c 2 by their respective stochastic counterparts c 1 + σ 1 where B 1 (t) and B 2 (t) are mutually independent Brownian motions, σ 1 and σ 2 represent the intensities of the white noise. Thus the stochastic version of the deterministic System (1) is given by the following Itô equation: In reality, due to environmental fluctuations, not only c 1 and c 2 that can be affected by a random noise, but also all the parameters involved with the deterministic system exhibit random fluctuations to a greater or lesser extent. However, perturbing all the parameters in System (1) makes the analysis intractable because of the strong non-linearity of the resulting stochastic system. The rest of this paper is organized as follows. In the next section, we discuss the existence of equilibrium states for model (1). The third section will deal with the global behavior of the deterministic System (1). We prove by constructing a Lyapunov functions that the equilibria are globally asymptotically stable. In the fourth section, we show that there is a unique global and positive solution of (2). We also show that the stochastic Systems (2) will be stochastically ultimately bounded. Thereupon, we investigate the asymptotic properties of the solution. Finally, a discussion and numerical simulations are presented to confirm our mathematical findings.

Equilibria
Equilibria for System (1) can be found by setting the right sides of the three differential equations of (1) equal to zero, giving the algebraic system From the third equation of (3) we have Using the first equation of (3) and (4) we get Substituting (5) into the second equation of (3). This gives the quadratic equation where If p 1 = 0, one root of (6) is I 0 0 = 0 then (1) where there is no new infective members of immigration. The second root of (6) is and is the endemic equilibrium state of the model (1) without new infective members of immigration which exists provided that the reproduction number R 0 > 1 [19].
If p 1 > 0, the quadratic Eq. (6) has one positive and one negative root. The disease-free equilibrium state E 0 0 that occurs when p 1 = 0 now becomes negative (not biologically feasible). The positive root is Thus, System (1) has unique endemic equilibrium state E * p 1 (S * , I * , R * ) in the presence of new infective members of immigration, where S * and R * are given by (5) and (4), respectively.

Analysis of the deterministic model
It is easy to show that System (1) is well posed, in the sense that if S(0), I (0) and R(0) are positive, then there exists a unique solution and S(t), I (t) and R(t) are positives for all t. By summing all equations of System (1) we find that the total population size verifies the equation, It is convenient to use I , R and N as variables and replace S by N − I − R. This gives the following system: Throughout this paper, we use the following notation 3.1 Global stability of the endemic equilibrium states In [13], the global stability of the endemic equilibrium was proved only for the case of c 2 = p 1 = p 2 = 0.
Here, we obtain the global stability of the endemic equilibrium for c 2 , p 1 , p 2 = 0 by constructing another Lyapunov function.
Proof Changing the variables such that For constructing a Lyapunov function for the trivial solution of System (3.1), we need the following elementary functions: Thus, the derivative of these functions along the solution of (4) is given bẏ Consider now the following function where A 1 , A 2 and A 3 are positive constants to be determined in the rest of this proof. From (10)(11)(12) and (13) we haveV It is easy to see that System (15) admits the unique solution given by V is positive definite andV is negative definite. Therefore the function V is a Lyapunov function for System (4) and consequently, by Lyapunov asymptotic stability theorem [11], the unique endemic equilibrium state E * p 1 is globally asymptotically stable as same as E * 0 if R 0 > 1.

Global stability of the free-diseases equilibrium state
Theorem 3.2 The disease-free equilibrium state E 0 0 is globally asymptotically stable in R 3 + , whenever R 0 ≤ 1.
Proof Let p 1 = 0. Define a C 2 -function W : R 3 + → R + by where N 0 0 = S 0 0 + R 0 0 and A 1 , A 2 , A 3 are the constants defined by (16). We get the derivative of W along the solution of (8) using β S 0 0 = (μ + c 1 + α)R 0 and the Eq. (16) verified by the constants A 1 , A 2 and A 3 we obtaiṅ Then, under the condition R 0 ≤ 1,Ẇ will be negative definite which leads to W being a Lyapunov function with respect to the disease-free equilibrium E 0 0 . This completes the proof.

Stochastic model analysis
Throughout the reset of this paper, let ( , F, 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).

Existence of the global and positive solution
To investigate the dynamical behavior of System (2), the first concern is whether the solution exists globally. Moreover, for a model of epidemic population dynamics, whether the value is positive is also requested. The following theorem shows that the solution of System (2) is global and positive.

Theorem 4.1 For any initial value X
, R(t)) to stochastic differential System (2) on t ≥ 0 and the solution will remain in R 3 + with probability 1, namely X (t) ∈ R 3 + for all t ≥ 0 almost surely.
Proof Since the coefficients of System (2) are locally Lipschitz continuous, for any initial value X 0 ∈ R 3 + , there is a unique local solution on [0, τ e ) where τ e is the explosion time [16]. To show this solution is global, we need to have τ e = ∞ almost surely (briefly a.s.). Define the stopping time τ + by with the traditional setting inf ∅ = ∞, where ∅ denotes the empty set.

Stochastically ultimate boundedness
Theorem 4.1 shows that the solutions to System (2) will remain in R 3 + . This property makes us continue to discuss how the solution varies in R 3 + in more detail. We first present the definition of stochastic ultimate boundedness (see, e.g. [12])

Definition 4.2
The solution X = (S, I, R) of System (2) is said to be stochastically ultimately bounded, if for any ∈ (0, 1), there is a positive constant η = η( ), such that for any initial value X 0 ∈ R 3 + , the solution X to (2) has the property that lim sup t→∞ P( X > η) ≤ .

Theorem 4.3 For any initial value X
Proof Let k 0 > 0 be sufficiently large such that every component of X 0 is contained within the interval ( 1 k 0 , k 0 ). For each integer k ≥ k 0 , define the stopping time which by Theorem 4.1 has the properties that, τ k → ∞ almost surely as k → ∞. We have Applying Itô's formula to e μt N gives By integrating this inequality and then taking expectations on both sides, one can see that Letting k → ∞ shows that Thus for any given > 0, let η = μ , by virtue of Markov's inequality, we can derive that lim sup

Asymptotic behavior around the equilibria of the deterministic system
In the deterministic models, the biological significant of the global stability of the disease-free equilibrium state is that the disease will be extinct while the global stability of the endemic equilibrium state signifies that the disease will persist in a population. Whereas, the stochastic System (2) has no equilibrium states. So, how do we measure whether the disease will die out or persist?. In this paragraph, following Imhof and Walcher [8], we show the difference between the solution of System (2) and E * p 1 is small if the intensities of noises are low to reflect that the disease is prevalent. Besides, if the reproduction number R 0 ≤ 1, we estimate the oscillation around E 0 0 to reflect whether the disease will die out. and where A 1 , A 2 and A 3 are defined by (16).
Proof As in the analysis of deterministic System (1), utilizing the variables I , R and N and changing the variables such that x = I − I * , y = R − R * and z = N − N * , System (2) can be written as By applying Dynkin's formula (see, e.g. [17]) to the Lyapunov function V defined by (14) we obtain where τ k is the stoping time defined by (21) and L is the differential operator associated to System (22) (see also [17]). Since τ k → ∞, so we get by letting k → ∞ where V 0 = V (x(0), y(0), z(0)). Let us now compute LV (x, y, z) which can be simplified to Applying the basic inequality (a + b) 2 ≤ 2a 2 + 2b 2 to (x + I * ) 2 and (y + R * ) 2 and utilizing the expression (17) ofV . It then follows from (24) after simplification that Combining (23) and (25) gives It then follows from the positivity of V that lim sup The proof is complete.
Remark If p 1 = 0, we change the equilibrium state (I * , R * , N * ) in the result of Theorem 4.4 by the equilibrium state (I * 0 , R * 0 , N * 0 ) provided that the reproduction number R 0 > 1.
Proof The differential operator L acting on the Lyapunov function W defined by (18) gives It follows from the expression (19) ofẆ that By Dynkin's formula, W (I, R, N ) satisfies which implies together with the utilization of the positivity of W and (26) that Dividing both sides by t and then letting t → ∞ we get the required result.

Discussion and numerical simulations
This work studies a deterministic and a stochastic epidemic model of SIRS type which includes constant inflows of new susceptibles, infectives, removeds and saturation incidence rate. If the constant inflow of the infectives is null in the deterministic model (1), we find the basic reproduction number R 0 . It determines completely the dynamical behavior of an epidemic describing by (1). If R 0 ≤ 1 the disease will die out (see Fig. 1). If not (i.e., R 0 > 1), the disease will persist at the endemic equilibrium state level (see Fig. 2). If there is a positive flow of infectives into the population, the model has unique endemic equilibrium state approached by all solutions. Therefore, the infection cannot be eliminated from the population (see Fig. 3). It is worth mentioning that from the expression (7), I * is a decreasing function of p 1 and we have lim Therefore, we can control an epidemic and reduce the force of infection by reducing the reproduction number R 0 and the positive constant of the new infective individuals p 1 . Noting that R 0 can be reduced by reducing the contact rate β, increasing the fraction of the new removeds in the population or increasing the duration of lose of immunity 1 γ , all this by taking effective measures such as quarantine, isolation, vaccination, etc.
On the other hand, population dynamics is inevitably subjected to environmental noise. So, it is important to examine the inclusion of stochastic effects into deterministic models. In this paper a stochastic model is developed based on the deterministic formulation by perturbing randomly two parameters, namely the diseaserelated death rates of I and R. We show that there is a unique global positive solution to the stochastic model (2) for any positive initial value. Moreover, we show that this solution is stochastically bounded. Next we investigated the behavior of the stochastic solution. In contrast to the deterministic solutions, the stochastic solutions do not converge to one of the equilibria of the deterministic System (1). However, Theorems 4.3 and 4.4 relate the behavior of the stochastic system to the asymptotic deterministic behavior. Indeed, if the intensities of noise are sufficiently small, the stochastic solution can be expected to remain close to the diseasefree or the endemic equilibrium state (see Figs. 1, 2, 3). However, in the case when R 0 ≤ 1 and there is no immigration of infected individuals, that is p 1 = 0, it seems likely that the number of infectives will tend to 0 almost surely as t goes to infinity. We leave the study of this statement for future investigation.
In summary, the theoretical analysis and numerical simulations in this study show that adding noise to deterministic model affects the stability and able to change the dynamics of the model system from stable situation to unstable one. Moreover, as illustrated in Figs. 4, 5, 6, the strong noise may give a divergence between stochastic and deterministic behaviors.  (1) and (2) for the same parameters of Fig. 3 except σ 1 = 4, σ 2 = 5