Dynamics of an SLIR model with nonmonotone incidence rate and stochastic perturbation

Abstract: In this paper we study an SLIR epidemic model with nonmonotonic incidence rate, which describes the psychological effect of certain serious diseases on the community when the number of infectives is getting larger. By carrying out a global analysis of the model and studying the stability of the disease-free equilibrium and the endemic equilibrium, we show that either the number of infective individuals tends to zero or the disease persists as time evolves. For the stochastic model, we prove the existence, uniqueness and positivity of the solution of the model. Then, we investigate the stability of the model and we prove that the infective tends asymptotically to zero exponentially almost surely as R0 < 1. We also proved that the SLIR model has the ergodic property as the fluctuation is small, where the positive solution converges weakly to the unique stationary distribution.


Introduction
The incidence of a disease is the number of individuals who become infected per unit of time and plays an important role in the study of mathematical epidemiology. Generally, incidence rate depends on both the susceptible and infective classes. In many epidemic models, the bilinear incidence rate is frequently used [1]. However, in recent years, researchers have taken into account oscillations in incidence rates and proposed many nonlinear incidence rates. Let S (t) be the number of susceptible individuals, I(t) be the number of infective individuals, and R(t) be the number of removed individuals at time t, respectively. After studying the cholera epidemic spread in Bari in 1973, Capasso and Serio [2] introduced the saturated incidence g(I)S into epidemic models, where g(I) is decreasing when I is large enough, i.e. g(I) = kI 1 + αI (1.1) This incidence rate seems more reasonable than the bilinear incidence rate because it includes the behavioral change and crowding effect of the infective individuals and prevents the unboundedness of the contact rate by choosing suitable parameters.
To incorporate the effect of the behavioral changes of the susceptible individuals, Liu et al. [3] proposed the general incidence rate g(I) = kI p 1 + αI q (1.2) where p and q are positive constants and α is a nonnegative constant. The special cases when p, q and k take different values have been used by many authors. For example, Ruan and Wang [4] studied the case when p = 2, q = 2 i.e. g(I) = kI 2 1+αI 2 and got some complicated dynamics phenomenons, such as saddle-node bifurcation, Hopf bifurcation, Bogdanov-Takens bifurcation and the existence of none, one and two limit cycles. Derrick and van den Driessche [5] , Hethcote [6], Alexander and Moghadas [7], etc. also used the general incidence rate.
If the function g(I) is nonmonotone, that is, g(I) is increasing when I is small and decreasing when I is large, it can be used to interpret the 'psychological' effect: for a very large number of infective individuals the infection force may decrease as the number of infective individuals increases, because in the presence of large number of infectives the population may tend to reduce the number of contacts per unit time [8]. To model this phenomenon, Xiao and Ruan [8] proposed a incidence rate g(I) = kI 1 + αI 2 (1.3) where kI measures the infection force of the disease and 1/(1 + αI 2 ) describes the psychological or inhibitory effect from the behavioral change of the susceptible individuals when the number of infective individuals is very large. This is important because the number of effective contacts between infective individuals and susceptible individuals decreases at high infective levels due to the crowding of infective individuals or due to the protection measures by the susceptible individuals. Notice that when α = 0, the nonmonotone incidence rate (1.3) becomes the bilinear incidence rate [8]. They used this incidence rate in an SIR model and got the result that this model does not exhibit complicated dynamics as other epidemic models with other types of incidence rates reported in. In this paper, we use this incidence rate in an SLIR model. In fact, epidemic models are inevitably affected by environmental white noise which is an important component in realism, because it can provide an additional degree of realism in compared to their deterministic counterparts. Therefore, many stochastic models for the epidemic populations have been developed. In addition, both from a biological and from a mathematical perspective, there are different possible approaches to include random effects in the model. Here, we mainly mention three approaches. The first one is through time Markov chain model to consider environment noise in HIV epidemic (see, e.g., [9][10][11][12]). The second is with parameters perturbation. There is an intensive literature on this area, such as [13][14][15][16][17][18][19]. The last important issue to model stochastic epidemic system is to robust the positive equilibria of deterministic model. In this situation, it is mainly to investigate whether the stochastic system preserves the asymptotic stability properties of the positive equilibria of deterministic models, see [20][21][22]. Recently, Yang et.al [19] discusses the stochastic SIR and SEIR epidemic models with saturated incidence rate βS /(1 + αI). In this paper, we introduce randomness in the SLIR model with the second approaches as [13] and [15].
The rest of the paper is organized as follows: in Section 2, the deterministic SLIR mathematical model is formulated, boundedness of solutions and existence of a positively invariant and attracting set are shown. In Section 3, the basic reproductive number, the conditions to the existence of possible equilibria of the system and their stability are established. In section 4, we obtain the analytic results of dynamics of the SDE model. Finally, a brief discussion is given in Section 5.

The deterministic model
In this section, we formulate an SLIR model with incidence rate (1.3). Figure 1 shows the model diagram. The total population at time t, denoted by N(t), is divided into four classes: susceptible (S ), latent (L), infectious (I) and treatment (R). All recruitment is into the susceptible class, and occurs at a constant rate Λ. We assume that an individual may be infected only through contacts with infectious individuals. The natural death rate is µ. The infectious class has an additional death rate due to the disease with rate constant d. Thus individuals leave class L for class I at rate k. Latent and Infectious individuals are treated with rate constant r 1 and r 2 , entering the treatment class,respectively. A fraction p of newly infected individuals moves to the latent class (L), and the remaining fraction 1 − p, moves to the active class (I). The incidence rate is g(I) = βI 1+αI 2 . It is assumed that individuals in the latent class do not transmit infection. Combining all the aforementioned assumptions, the model is given by the following system of differential equations: By adding all Eq (2.1), the dynamics of the total population N(t) is given by: Since dN/dt < 0 for N > Λ/µ, then, without loss of generality, we can consider only solutions of (2.1) in the following positively subset of R 4 : With respect to model system (2.1), we have the following result: Proposition 2.1. The compact set Ω ε is a positively invariant and absorbing set that attracts all solutions of Eq (2.1) in R 4 .
Proof. Define a Lyapunov function as W(t) = S (t) + L(t) + I(t) + R(t), then we have: Hence, that dW dt ≤ 0 for W > Λ µ . Ω ε is a positively invariant set. On the other hand, solving the differential inequality Eq (2.3) yields: W(0) is the initial condition of W(t). Thus, as t → ∞, one has that 0 ≤ W(t) ≤ Λ µ . To analysis the system (2.1), we notice that the variable R does not participate in the first three equations. Thus we can consider only the equations for S , L and I, i.e. the following system:

Dynamics of the deterministic model
In this section, the model is analyzed in order to obtain the basic reproduction number, conditions for the existence and uniqueness of non-trivial equilibria and asymptotic stability of equilibria.

Basic reproductive number
In this subsection, we define the basic reproduction number R 0 of system (2.4). R 0 is the average number of secondary infections that occur when one infective is introduced into a completely susceptible population [23]. For many deterministic epidemiology models, an infection can get started in a fully susceptible population if and only if R 0 > 1. Thus the basic reproductive number R 0 is often considered as the threshold quantity that determines when an infection can invade and persist in a new host population [1].
The disease-free equilibrium of system (2.4) is X 0 = (S 0 , 0, 0) with S 0 = Λ/µ. In order to compute the basic reproduction number, it is important to distinguish new infections from all other class transitions in the population. The infected classes are L and I. Following Van den Driessche and Watmough [24], we can rewrite system (2.4) aṡ where x = (L, I, S ), F is the rate of appearance of new infections in each class, V + is the rate of transfer into each class by all other means and V − is the rate of transfer out of each class. Hence, The jacobian matrices of F and V at the disease-free equilibrium X 0 = (0, 0, Λ/µ) can be partitioned as where F and V correspond to the derivatives of F and V with respect to the infected classes: The basic reproduction number is defined, following Van den Driessche and Watmough [24], as the spectral radius of the next generation matrix, FV −1 :

Stability of the disease-free equilibrium
We have the following result about the global stability of the disease free equilibrium: Theorem 3.1. When R 0 > 1, the disease free equilibrium X 0 is unstable.When R 0 ≤ 1, the disease free equilibrium X 0 is globally asymptotically stable in Ω ε ; this implies the global asymptotic stability of the disease free equilibrium on the nonnegative orthant R 3 . This means that the disease naturally dies out.
Proof. The Jacobian matrix of system (2.4) at X 0 is and the characteristic equation is From (3.1), we get the discriminant of f (λ) is Therefore, (3.1) has three real roots. If R 0 < 1, we have (3.1) has three negative real roots and hence X 0 is locally stable. If R 0 > 1, (3.1) has at least one positive real root and hence X 0 is unstable. When R 0 ≤ 1, we can define the following Lyapunov-LaSalle function Its time derivative along the trajectories of system (2.4) satisfieṡ Thus R 0 ≤ 1 implies thatV ≤ 0. By LaSalle's invariance principle, the largest invariant set in Ω ε contained in {(S , L, I, R) ∈ Ω ε ,V = 0} is reduced to the disease-free equilibrium X 0 . This proves the global asymptotic stability of the disease-free equilibrium on Ω ε [25]. Since Ω ε is absorbing, this proves the global asymptotic stability on the nonnegative octant for R 0 ≤ 1. It should be stressed that we need to consider a positively invariant compact set to establish the stability of X 0 since V is not positive definite. Generally, the LaSalle's invariance principle only proves the attractivity of the equilibrium. Considering Ω ε permits to conclude for the stability [25][26][27]. This achieves the proof.

Existence and uniqueness of endemic equilibrium
To find the positive equilibrium, let Then we have the following result: When R 0 > 1, there exist an unique endemic equilibrium X * = (S * , L * , I * ) for the system (2.4) where S * , L * , and I * are defined as in (3.2) which is in the nonnegative octant R 3 + .

Locally stability of endemic equilibrium
In this section, we denote µ + k + r 1 = A, µ + d + r 2 = B for writing convenience.
Proof. The Jacobian matrix of system (2.4) at endemic equilibrium X * = (S * , L * , I * ) is The characteristic equation of J(X * ) is According to direct calculation we have a, c > 0 and ab > c when R 0 > 1 . Therefore the endemic equilibrium X * is locally asymptotically stable in Ω ε by Routh-Hurwitz criterion.

Globally stability of endemic equilibrium
In this section, we briefly outline a general mathematical framework developed in [28] for proving global stability, which will be used to prove Theorem 3.10. Consider the autonomous dynamical system where f : D → R n open set and f ∈ C 1 (D). Letx be an equilibrium of (3.4), i.e. f (x) = 0. We recall thatx is said to be globally stable in D if it is locally stable and all trajectories in D converge tox. Assume that the following hypothesis hold: (H1) D is simply connected; (H2) There exists a compact absorbing set Γ ⊂ D; (H3) Eq (3.4) has a unique equilibriumx in D.
The Global Stability Problem arising in [28] is to find conditions under which the global stability of x with respect to D is implied by its local stability. In [28], the main idea of this geometric approach to the global stability problem is as follows: Assume that (3.4) satisfies a condition in D, which precludes the existence of periodic solutions and suppose that this condition is robust, in the sense that it is also satisfied by ordinary differential equations that are C 1 -close to (3.4); Then every nonwandering point of (3.4) is an equilibrium, as otherwise, by the C 1 closing lemma of Pugh [29,30], we can perturb (3.4) near such a nonequilibrium nonwandering point to get a periodic solution. As a special case, every omega limit point of (3.4) is an equilibrium. This implies that x attracts points in D. As a consequence, its global stability is implied by the local stability. For this purpose, we introduce the notion of robustness of a Bendixson criterion under local C 1 perturbations of f .
Thus, for example, any equilibrium, alpha limit point, or omega limit point is nonwandering.
and | · | denotes a vector norm on R n and the operator norm that it induces for linear mappings from R n to R n . Definition 3.7. A Bendixson criterion is said to be robust under C 1 local perturbations of f at x 0 if, for each sufficiently small > 0 and neighborhood U of x 0 , it is also satisfied by each C 1 localperturbations g such that supp( f − g) ⊂ U.
The following is the local version of the global-stability principle proved by [28].
Theorem 3.8. Suppose that assumptions (H2) and (H3) hold and that (3.4) satisfies a Bendixson criterion that is robust under C 1 local perturbations of f at all nonequilibrium nonwandering points for (3.4), thenx is globally asymptotically stable with respect to D, provided it is stable.
The following is to introduce the quantityq 2 given in [28]. Assume that (3.4) has a compact absorbing set K ⊂ D . Then every solution x(t, x 0 ) of (3.4) exists for any t > 0. The quantityq 2 is well defined:q and x → P(x) is a n 2 × n 2 matrix-valued function, which is C 1 in D and A f = (DP)( f ) or, equivalently, P f is the matrix obtained by replacing each entry a i j in P by its directional derivative in the direction of f , ∂ f ∂x [2] is the second additive compound matrix of the Jacobian matrix ∂ f ∂x of f . If | · | is a vector norm on R n 2 , then ρ(B) is the Lozinskǐ measure with respect to | · |, which is defined by Under assumptions (H1) and (H2), [28] proved thatq 2 < 0 is a Bendixson criterion for (3.4) and it is robust under C 1 local perturbations of f at all nonequilibrium nonwandering points for (3.4) by means of the local version of C 1 closing lemma of Pugh [29,30]. Then we have the following theorem Theorem 3.9. Under assumptions (H1), (H2), and (H3), the unique equilibriumx is globally asymptotically stable in D ifq 2 < 0 .
The criterionq 2 < 0 provides the flexibility of a choice of n 2 × n 2 arbitrary function in addition to the choice of vector norms | · | in deriving suitable conditions. It has been remarked in [28] that, under the assumptions of Theorem 3.9, the classical result of Lyapunov, the classical Bendixson-Dulac condition, and the criterion in [31] are obtained as special cases [28]. In addition, it has also been stated in [28] that, under the assumptions of Theorem 3.9 the conditionq 2 < 0 also implies the local stability of the equilibriumx, because, assuming the contrary,x is both the alpha and the omega limit point of a homoclinic orbit that is ruled out by the conditionq 2 < 0.
Global stability of endemic equilibrium.
We will focus on investigating the globally asymptotical stability of the unique endemic equilibrium X * (S * , L * , I * ). The main theorem of the method requires the use of Lozinskǐ Logarithmic norm.
Theorem 3.10. Assume that R 0 > 1. Then the unique endemic equilibrium X * is globally asymptotically stable in R + 3 . Proof. We set P as the following diagonal matrix: Then P is C 1 and nonsingular in • Ω ε . Let f denote the vector field of (2.4). Then the second compound matrix J [2] of the Jacobian matrix of system (2.4) can be calculated as follows Then the matrix B = P f P −1 + PJ [2] P −1 in (3.6) can be written in matrix form where Let (u, v, w) denote the vectors in R 3 , we select a norm in R 3 as |(u, v, w)| = max{|u|, |v| + |w|} and let µ denote the Lozinskǐ measure with respect to this norm. Following the method in [32], we have the estimateμ(B) ≤ sup{g 1 , g 2 }, where Therefore The uniform persistence constant c can be adjusted so that there exists T > 0 independent of (S (0), L(0), I(0)) ∈ K, the compact absorbing set, such that Substituting (3.9) into (3.7) and using (3.10), we obtain, for t > T , Thereforeμ(B) ≤ L L + S S for t > T by (3.8) and (3.11). Along each solution (S (t), L(t), I(t)) to (2.4) such that (S (0), L(0), I(0)) ∈ K and for t > T , we have which impliesq 2 ≤ 0 from (3.5), proving Theorem 3.10.
In this section we proposed a nonmonotone and nonlinear incidence rate of the form βIS /(1 + αI 2 ), which is increasing when I is small and decreasing when I is large. It can be used to interpret the 'psychological' effect: the number of effective contacts between infective individuals and susceptible individuals decreases at high infective levels due to the quarantine of infective individuals or the protection measures by the susceptible individuals.
Recall that the parameter α describes the psychological effect of the general public toward the infectives. Though the basic reproduction number R 0 does not depend on α explicitly, numerical simulations indicate that when the disease is endemic, the steady state value I * of the infectives decreases as α increases (see Figure 2). From the steady state expression (3.2) we can see that I * approaches zero as α tends to infinity.

Stochastic differential equation model
In this section we consider the stochastic model associated with the deterministic model in (2.1). We introduce randomness into the model (2.1) by replacing the parameters µ by µ → µ + σ i dB i (t)(i = 1, 2, 3, 4) with the second approaches as [13] and [15], where B i (t)(i = 1, 2, 3, 4) are independent standard Brownian motions with B i (0) = 0 and σ 2 i (i = 1, 2, 3, 4) denote the intensities of the white noise. Then equation (2.1) becomes (4.1) Throughout this section, unless otherwise specified, let (Ω, F , P) be a complete probability space with a filtration {F t } t≥0 . The Brown motions are defined on the complete probability space (Ω, F , P). Denote R n + = {x ∈ R n : x i > 0, for all 1 ≤ i ≤ n}. In general, consider n-dimensional stochastic differential equation [35]: with initial value x(t 0 ) = x 0 ∈ R n . B(t) denotes n-dimensional standard Brownian motions defined on the above probability space. Define the differential operator L associated with (4.2) by If L acts on a function V ∈ C 2,1 (S h ×R + ;R + ), then 4.1. Existence and uniqueness of the positive solution of (4.1) In this subsection we first show the solution of system (4.1) is global and positive. In order for a stochastic differential equation to have a unique global (i.e. no explosion in a finite time) solution for any given initial value, the coefficients of the equation are generally required to satisfy the linear growth condition and local Lipschitz condition [34,35]. However, the coefficients of system (4.1) do not satisfy the linear growth condition (as the incidence is nonlinear), so the solution of system (4.1) may explode at a finite time [34,35]. In this section, using Lyapunov analysis method (mentioned in [36]), we show the solution of system (4.1) is positive and global.

Asymptotic behavior around the disease-free equilibrium of the deterministic model (2.1)
Obviously, X 0 = ( Λ µ , 0, 0, 0) is the solution of system (2.1), which is called the disease-free equilibrium. If R 0 < 1, then X 0 is globally asymptotically stable, which means the disease will vanish after some period of time. Therefore, it is interesting to study the disease-free equilibrium for controlling infectious disease. But, there is no disease-free equilibrium in system (4.1). It is natural to ask how we can consider the disease will extinct. In this subsection we mainly estimate the average of oscillation around X 0 in time to exhibit whether the disease will die out. , Proof. Define C 2 -function V 1 , V 2 : R + → R + and V 3 , V 4 , V 5 : R 2 + → R + , respectively by Along the trajectories of system (4.1), we have

Asymptotic behavior around the endemic equilibrium of the deterministic model (2.4)
In this section, we will show there is a unique distribution for system (4.1) instead of asymptotically stable equilibrium(see [37]). We only consider the first three equation of system (4.1) because the variable R does not participate in the first three equations. Before giving the main theorem, we first give a lemma (see [38]).
Let X(t) be a regular temporally homogeneous Markov process in E l ⊂ R l described by the stochastic differential equation: σ r (X)dB r (t), and the diffusion matrix is defined as follows Lemma 4.3. (see [38]) We assume that there exists a bounded domain U ⊂ E l with regular boundary, having the following properties: • In the domain U and some neighborhood thereof, the smallest eigenvalue of the diffusion matrix A(x) is bounded away from zero.
• If x ∈ E l \U, the mean time τ at which a path issuing from x reaches the set U is finite,and sup x∈K E x τ < ∞ for every compact subsect K ⊂ E l .
then, the Markov process X(t) has a stationary distribution µ(·) with density in E l such that for any Borel set B ⊂ E l , lim t→∞ P(t, x, B) =μ(B), and P x lim for all x ∈ E l and f (x) being a function integrable with respect to the measureμ.

Exponential stability of system (4.1)
In this section, we investigate the exponential decay of the global solution of system (4.1) as the intensity of white noise is great. It can be shown below, even if the endemic equilibrium exists in the system (2.1), the stochastic effect may make washout more likely in system (4.1).
It has been shown that X(t) converges weakly to distribution ν, so does S (t) as t → ∞.

Discussion
Based on the model proposed in Xiao and Ruan [8], we proposed an SLIR model with a nonmonotone and nonlinear incidence rate of the form βIS /(1 + αI 2 ) which is increasing when I is small and decreasing when I is large. It can be used to interpret the 'psychological' effect: the number of effective contacts between infective individuals and susceptible individuals decreases at high infective levels due to the quarantine of infective individuals or the protection measures by the susceptible individuals.
We have provided a complete description of the asymptotic behaviour of the solutions of the deterministic model (2.1) and (2.4). Interestingly, this model does not exhibit complicated dynamics as other epidemic models with other types of incidence rates reported in [1,[3][4][5][6]8]. In terms of the basic reproduction number R 0 = βΛ µ k(1−p)+(µ+k+r 1 )p (µ+k+r 1 )(µ+d+r 2 ) , our main results indicate that when R 0 < 1, the disease-free equilibrium is globally attractive. When R 0 > 1, the endemic equilibrium exists and is globally stable.
A stochastic differential equation (SDE) is formulated for describing the disease in this case. We prove the existence and uniqueness of the solution of this SDE. We proved the positivity of the solutions. Then, we investigate the stability of the model. We illustrated the dynamical behavior of the SLIR model according to R 0 < 1. We proved that the infective tends asymptotically to zero exponentially almost surely as R 0 < 1. We also proved that the SLIR model has the ergodic property as the fluctuation is small, where the positive solution converges weakly to the unique stationary distribution.