THE THRESHOLD OF A STOCHASTIC SIRS EPIDEMIC MODEL IN A POPULATION WITH VARYING SIZE

. In this paper, a stochastic susceptible-infected-removed-susceptible (SIRS) epidemic model in a population with varying size is discussed. A new threshold ˜ R 0 is identiﬁed which determines the outcome of the disease. When the noise is small, if ˜ R 0 < 1, the infected proportion of the population dis-appears, so the disease dies out, whereas if ˜ R 0 > 1, the infected proportion persists in the mean and we derive that the disease is endemic. Furthermore, when R 0 > 1 and subject to a condition on some of the model parameters, we show that the solution of the stochastic model oscillates around the endemic equilibrium of the corresponding deterministic system with threshold R 0 , and the intensity of ﬂuctuation is proportional to that of the white noise. On the other hand, when the noise is large, we ﬁnd that a large noise intensity has the eﬀect of suppressing the epidemic, so that it dies out. These results are illustrated by computer simulations.


1.
Introduction. Studies of epidemic models that incorporate disease-caused death and a varying total population have become one of the important areas in the mathematical theory of epidemiology, largely inspired by the works of Anderson and May (see [1], [2]).

YANAN ZHAO, DAQING JIANG, XUERONG MAO AND ALISON GRAY
Much of the classical work on epidemiological models has been restricted to situations where the affected population is of constant size. This assumption is relatively valid for diseases of short duration with limited effects on mortality. However, it clearly fails to hold for diseases that are endemic in communities with changing populations, and for diseases which raise the mortality rate substantially. Wellknown examples of such diseases are malaria in developing countries with growing populations, and the current AIDS pandemic. In such situations, the effects of the disease-induced mortality and of the change in population size are far from negligible, and in fact, may have a crucial influence on whether or not the disease can reach epidemic levels.
In recent years there have been a number of studies on disease transmission models in populations of varying size, and some of these have given a complete global analysis of the model equations. We briefly note the work on such models that is related to ours as including the papers [3,5,7,10,15,16,17,19].
Busenberg and van den Driessche (1990) [6] discussed a SIRS epidemic model in a population with varying size. The system has the following form: Here S, I and R denote the total numbers of susceptible, infective and recovered (removed) individuals respectively, and N is the population size. We set N (t) = S(t) + I(t) + R(t) at time t, andṠ(t),İ(t) andṘ(t) are the corresponding rates of change with respect to time t. We use the following additional notation: b: per capita birth rate, µ: per capita disease free death rate, β: effective per capita contact rate of infective individuals, δ: per capita loss of immunity rate of recovered individuals, α: excess per capita death rate of infected individuals, γ: per capita recovery rate of infected individuals.
All parameter values are assumed to be nonnegative and b, γ > 0. Note that the system (1) reduces to the SIR model when δ = 0, so there is no loss of immunity. The equation for the total population size isṄ (t) = (b − µ)N − αI. Busenberg and van den Driessche (1990) [6] consider the proportions of individuals in the three epidemiological classes, namely It is easy to verify that x, y and z satisfy the system of differential equations     ẋ (t) = b − bx − βxy + δz + αxy, y(t) = βxy − (b + α + γ)y + αy 2 , z(t) = γy − (b + δ)z + αyz. ( The feasibility region becomes Γ = {x ≥ 0, y ≥ 0, z ≥ 0, x + y + z = 1} and Γ 0 = Γ − {(1, 0, 0)}. There are two distinct ways of considering a disease as being brought under control in a population of increasing or decreasing total size. The stricter way requires that the total number of infectives I(t) → 0, while a weaker requirement is that the proportion y(t) → 0. This distinction is discussed in some detail in Busenberg et al. (1991) [7]. Thus, they showed the conditions for the existence and stability of the endemic proportion steady state (x * , y * , z * ) and for the stability of the disease-free steady state (x, y, z) = (1, 0, 0). The threshold parameter of the system is which determines the extinction and persistence of the epidemic. According to the theory in Busenberg and van den Driessche (1990) [6], then (a) The disease free equilibrium proportion (x, y, z) = (1, 0, 0) always exists, and is globally asymptotically stable in the feasibility region Γ whenever R 0 ≤ 1, and it is unstable when R 0 > 1.
However the deterministic approach has some limitations in the mathematical modelling transmission of an infectious disease. Stochastic differential equation (SDE) models play a significant role in various branches of applied sciences including infectious dynamics, as they provide some additional degree of realism compared to their deterministic counterpart. Recently, Many authors have introduced parameter perturbation into epidemic models and have studied their dynamics [4,8,21,22].
Using stochastic differential equation models is one way to do this. However, compared to deterministic systems, it is extremely difficult to give the threshold of stochastic systems. Recently, Gray et al. [11] discussed the following stochastic SIS epidemic model and examined its dynamics. Here B(t) is a standard Brownian motion with B(0) = 0, the white noise intensity is σ 2 > 0, and µ, β and γ are as above. In this model S and I denote the numbers of susceptible and infected individuals in the population, respectively. As S(t) + I(t) ≡ N , the authors simplified system (5) into a single equation, and showed that • If R S 0 < 1 and σ 2 ≤ β/N , then the disease will die out with probability one. • If R S 0 > 1, then the disease will be persistent, where R S 0 = R D 0 −σ 2 N 2 /2(µ + γ) and R D 0 = βN /µ + γ is the threshold of the corresponding deterministic model. The quantity R S 0 can be considered as the threshold of system (5), which is less than the value of R D 0 . The results were illustrated by computer simulations. Tornatore, Buccellato and Vetro [18] discussed a stochastic SIR system, given by They proved that 0 < β < min{λ + µ − σ 2 /2, 2µ} is a sufficient condition for the asymptotic stability of the disease-free equilibrium (DFE). Their computer simulations for the SDE SIR model agree well with the analytical results and show that the introduction of noise into the system raises the threshold to µ + γ + σ 2 /2, so if then the DFE E 0 = (S(0), I(0), R(0)) = (1, 0, 0) is stable and the disease does not occur, whereas if β > µ + γ + σ 2 /2, then the DFE is unstable. Ji, Jiang and Shi [13] discussed system (6) further, and obtained the results that • If β < γ + µ − σ 2 /2, then the disease-free equilibrium (1, 0, 0) of system (6) is stochastically asymptotically stable in the large and is exponentially mean-square stable.
• If β > γ + µ, then the solution of system (6) fluctuates around (S * , I * , R * ), which is the endemic equilibrium of the corresponding deterministic system. They consider that the disease will prevail if the white noise is small, but they did not give an accurate threshold [13,18].
Recently Yang and Mao [20] studied a class of multi-group SEIR epidemic models with stochastic perturbations. By the method of stochastic Lyapunov functions, they studied the extinction and recurrence in terms of the intensity of the stochastic perturbations and the reproductive number R 0 .
In this paper, we suppose that some stochastic environmental factor acts simultaneously on each individual in the population. In this case β changes to a random variableβ, as in [11]. More precisely each infected individual makes potentially infectious contacts with each other individual in [t, t+dt). Here dB(t) = B(t + dt) − B(t) is the increment of a standard Brownian motion. Thus the number of potentially infectious contacts that a single infected individual makes with another individual in [t, t + dt) is normally distributed with mean βdt and variance σ 2 dt. Hence E(βdt) = βdt and var(βdt) = σ 2 dt. As var(βdt) → 0 as dt → 0 this is a biologically reasonable model. Indeed this is a well-established way of introducing stochastic environmental noise into biologically realistic population dynamic models.
The stochastic model corresponding to the deterministic SIRS model (1) then takes the following form, as an SDE model: The equation for the total population N size is obtained from (7) as: as in model (1). When the proportions x, y, z given by (2) are used as the independent variables, the system (7) becomes We will try to give a threshold of the system (8), which can easily determine the extinction and persistence of the disease. Using the relation x = 1 − y − z, we can omit analysis of the first equation of system (8), and so we discuss the following with any given initial value (y(0), z(0) ∈ R 2 + and y(0) + z(0) < 1. This paper is organized as follows. In Section 2, we show there is a unique positive solution of system (7) using the method mentioned in [9,11]. In Section 3, we deduce the conditions which will cause the disease to die out. The condition for the disease being persistent is given in Section 4. In Section 5, when R 0 > 1 and α 2 is sufficiently small, we derive that the solution of system (9) oscillates around the endemic proportion equilibrium P * (y * , z * ), and the intensity of fluctuation is proportional to the size of the white noise. The key to our analysis is choosing an appropriate Lyapunov function. Section 6 gives a short conclusion. Throughout the paper, outcomes of numerical simulations are shown in support of our analytical results.
Throughout this paper, unless otherwise specified, let (Ω, {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 right continuous and F 0 contains all P-null sets), and let B(t) be a scalar Brownian motion defined on this probability space.

Existence and uniqueness of positive solution.
To investigate the dynamical behaviour of a population model, the first concern is whether the solution of the model is positive and global. As we know, in order to get a stochastic differential equation for which a unique global solution exists (i.e. there is no explosion within a finite time) for any initial value, the coefficients of the equation should satisfy the linear growth condition and the local Lipschitz condition (cf. Mao [14]). However, the coefficients of system (7) do not satisfy the linear growth condition (as the incidence is nonlinear), so the solution of system (7) may explode at a finite time. In this section, using the Lyapunov analysis method (as mentioned in [9,11]), we show that the solution of system (7) is positive and global.
Moreover, substituting this into (12), we obtain By the Itô formula, for any t ∈ [0, T ] and k ≥ k 1 , we have where E is the expectation operator and a ∧ b denotes the minimum of a and b.
Remark 1. Theorem 2.1 shows that for any initial value (S(0), I(0), R(0)) ∈ R 3 + , there is a unique global solution (S(t), I(t), R(t)) ∈ R 3 + almost surely of system (7). Hence by (2), the region is a positively invariant set of system (8), which is similar to the feasibility region Γ of the deterministic system (3). Moreover we can state the following theorem.
3. Extinction. When studying dynamical systems for epidemics, two of the most important issues are extinction and persistence. We will discuss the extinction of the system (9) in this section but leave its persistence to the next section.
Case 3. IfR 0 < 1 and σ 2 ≤ β − α, then whereR and R 0 is defined in (4). That is, the proportion of infectives y(t) tends to zero exponentially a.s. In other words, the disease will die out with probability one. In addition, for the proportion of recovered individuals z(t), Proof. Applying the Itô formula to system (9) leads to d(log y) where f : (0, 1) → R is defined by Integrating this from 0 to t and dividing by t on both sides, we have Case 1. Compute which is negative by the condition σ 2 > max{(β − α) 2 /2(b + γ), β −α}. It therefore follows from (24) and (25) where M 1 (t) := σ t 0 (1 − y − z)dB(r). By the large number theorem for martingales (see e.g. [14]), Taking the limit superior of both sides, that we obtain the desired assertion (18) lim sup Case 2. We use the same notation as in the proof of Case 1. If β − α ≤ 0, by (23), it is easy to see that f (x) ≤ f (0) = −(b + γ). It then follows from (24) that In the same way as in the proof of Case 1, we have Case 3. We consider the quadratic function f : (0, 1) → R defined by (23). By (25) and the condition σ 2 ≤ β − α, it is easy to see that x = β − α σ 2 ≥ 1, and then f (x) takes its maximum value whereR 0 is defined by (21). It follows from (24) and (27) that
Next, let us prove the last assertion (22). According to (28), the last equation of the system (9) is an asymptotically differential system with limit system So we obtain Therefore, lim t→∞ z(t) = 0 a.s.
This finishes the proof of Theorem 3.1.
Remark 2. The condition of Case 3 states that the disease will die out ifR 0 < 1 and the white noise is not large, while if the white noise is large enough such that the condition of case 1 is satisfied, then the disease will also die out. Moreover, we notice the expression ofR 0 in (21). There is a difference compared with the threshold R 0 of system (3). In other words, the conditions for y(t) to become extinct in the SDE model (9)  Noting that with any initial value (y(0), z(0)) ∈ (0, 1) × (0, 1). That is y(t) will tend to zero exponentially with probability one and z(t) will tend to zero a.s. On the other hand, for the corresponding deterministic model, so the endemic proportion equilibrium (y * , z * ) is globally asymptotically stable in Γ 0 and the disease persists. Using the Euler-Maruyama (EM) method (see [12]), we give the simulations shown in Figure 1 to support our results. Example 2. We keep the system parameters the same as in Example 1 but let σ = 0.8. It is easy to verify that the system parameters obey the condition of Case 1 of Theorem 3.1, as We can therefore conclude that for any initial value (y(0), z(0)) ∈ (0, 1) × (0, 1), the solution of (9) obeys lim sup That is, y(t) will tend to zero exponentially with probability one. The computer simulation shown in Figure 2 clearly supports this result, showing extinction of the disease.
By the Itô formula, we have Integrating this from 0 to t and dividing by t on both sides, we have where (30) and the Schwarz inequality are used. Therefore, rearranging gives Noting that −∞ < log y(t) < 0 (as 0 < y(t) < 1) and using (26) and (32), we see that lim t→∞ φ(t) = 0 a.s.
:=ỹ * a.s., Finally, according to the last equality of (29), we get This finishes the proof of Theorem 4.2.
Remark 3. From Theorems 3.1 and 4.2, we can see that when the noise is sufficiently small that σ 2 ≤ β − α, if the value ofR 0 is above 1 or below 1 then the disease will persist or die out respectively. Therefore, we considerR 0 as the threshold of system (9). That is, we keep all the parameters the same as in Example 1, except that σ is reduced to 0.1 from 0.49. Note that so by Theorem 4.2, for any initial value (y(0), z(0)) ∈ (0, 1) × (0, 1), we conclude that the solution (y(t), z(t)) of system (9)  That is to say, the disease persists. The simulations shown in Figure 3 illustrate this.

5.
Asymptotic behavior around the endemic proportion equilibrium. In the deterministic model, the endemic proportion equilibrium P * (y * , z * ) is globally asymptotically stable. In this section, we show the effect of stochastic fluctuations of the environment on the endemic proportion equilibrium of the corresponding deterministic system.
Proof. Since P * (y * , z * ) is the endemic equilibrium of the corresponding deterministic model of system (9), we have that b + γ + α = β(1 − y * − z * ) + αy * , γy * + αy * z * = (b + δ)z * . Define a C 2 -function V : (0, 1) × (0, 1) → R + by where V 1 = y − y * − y * log(y/y * ) and V 2 = (z − z * ) 2 /2. The nonnegativity of this function can be observed from u − 1 − log u ≥ 0 on u > 0. Let L be the generating operator of system (9). Then we get where C is a positive constant. Therefore, although the solution of system (9) does not have stability like the deterministic system, there is approximate stability, provided that σ 2 is sufficiently small. That is, we keep all the parameters the same as in Example 1, except for σ which is now much smaller. Note that Then by Theorem 5.1, for any initial value (y(0), z(0)) ∈ (0, 1) × (0, 1), we conclude that the difference between the perturbed solution (y(t), z(t)) of system (9) and P * (y * , z * ) is related only to the level of the white noise, under the conditions R 0 > 1 and α 2 < 4γ 2 y * (β − α)/βz * . Using the EM method mentioned in [12], we show simulations to support our result. As expected, the solution oscillates around the endemic equilibrium P * for a long time (see Figure 4). The parameters used for the first two plots in Figure 4 are all the same but with different intensities of the white noise σ, i.e. in the first one σ = 0.05 and in the second one σ = 0.01. We observe clearly that as the white noise gets weaker, the fluctuation around P * gets smaller, which supports the result of Theorem 5.1.  6. Conclusions. In this paper, we have looked at an SDE version of the classical SIRS epidemic model in a population with varying size, with noise introduced into the disease transmission term. We showed that the SDE has a unique positive global solution, and we established conditions for the extinction and persistence of the disease. A key parameter is the thresholdR 0 , which is less than the corresponding deterministic version R 0 . Theorem 3.1 shows that ifR 0 < 1, under mild extra conditions the disease will die out. Theorem 4.2 shows that ifR 0 > 1, then the disease will persist. We also showed that if the conditions of Theorem 5.1 are satisfied, then the SDE solution oscillates around the endemic equilibrium of the deterministic system, and the intensity of fluctuation is proportional to the white noise intensity. Throughout the paper we have illustrated our theoretical results with computer simulations.