doi:10.1155/2008/636153 Research Article Analysis of a Delayed SIR Model with Nonlinear Incidence Rate

An SIR epidemic model with incubation time and saturated incidence rate is formulated, where the susceptibles are assumed to satisfy the logistic equation and the incidence term is of saturated form with the susceptible. The threshold value determining whether the disease dies out is found. The results obtained show that the global dynamics are completely determined by the values of the threshold value and time delay (i.e., incubation time length). If is less than one, the disease-free equilibrium is globally asymptotically stable and the disease always dies out, while if it exceeds one there will be an endemic. By using the time delay as a bifurcation parameter, the local stability for the endemic equilibrium is investigated, and the conditions with respect to the system to be absolutely stable and conditionally stable are derived. Numerical results demonstrate that the system with time delay exhibits rich complex dynamics, such as quasiperiodic and chaotic patterns.


Introduction
Epidemic models with nonlinear incidence have been studied by many authors, and related literature of SIR disease transmission model is quite large, where S denotes the number of individuals that are susceptible to infection, I denotes the number of individuals that are infectious, and R denotes the number of individuals that have been removed with immunity.For example, a detailed dynamical analysis of the nonlinear incidence rate βI p S q where β is the average number of contacts per infective per unit time is given by Liu et al. 1, 2 , Hethcote and van den Driessche 3 , Moghadas and Alexander 4 , Korobeinikov and Maini 5 , and others.After a study of the cholera epidemic spread in Bari in 1973, Capasso and Serio 6 introduced a saturated incidence rate g I S into epidemic models, the incidence rate seems more reasonable than βI p S q because the number of effective contacts between infective individuals and susceptible individuals may saturate at high infective levels due to crowding of infective individuals or due to the protection measures by the susceptible individuals 7 , such as the incidence rate β I p S/ 1 αI q see 1, 8 and references therein .Zhang and Chen 9 investigated a class of SIR epidemiological models under assumption that the susceptible satisfies the logistic equation and the incidence rate is of the form βIS q .More recently, Zhang et al. 10 , serve as an extended version to 9 , have carried out a long-term qualitative analysis incorporating incubation time delay into incidence rate in the case of q 1, that is, with the force of infection βS t I t−τ , which was proposed by Cooke 11 .The incubation period τ τ > 0 is a time, during which the infectious agents develop in the vector, and only after that time the infected vector becomes itself infectious.The detailed biological meanings and transmission mechanisms were given in 11 .The results obtained in 10 represent that, if the epidemic is persistent, introducing time delay changes the dynamical behaviors of the epidemic state.If the threshold value determining whether the disease dies out is larger than one and less than three, the endemic equilibrium is absolutely stable in the sense that it is asymptotically stable for all values of the delays 12 ; when it exceeds three, the endemic equilibrium is conditionally stable i.e., it is asymptotically stable for the delays in some intervals , and limit cycles arise by Hopf-type bifurcation with increasing time delay.
In 1978, May and Anderson 13 proposed the saturated incidence rate of the form β SI/ 1 αS , and used by some authors 14-17 , recently.The effect of saturation factor refer to α stems from epidemical control.In the absence of effective therapeutic treatment and vaccine, the epidemical control strategies are based on taking appropriate preventive measures.For example, if transmission vector is mosquito, these measures include mosquito reduction mechanisms and personal protection against exposure to mosquitos.Mosquito reduction mechanisms entail the elimination of mosquito breeding sites such as clearing culverts, roadside ditches, eliminating standing water, etc. , larvaciding killing of larvae before they become adults and adulticiding killing of adult mosquitoes by spraying .On the other hand, personal protection is based on preventing vector mosquitoes from biting humans by using mosquito repellents, avoiding locations where mosquitoes are biting, and using barrier methods such as window screens and long-sleeved clothing 18-20 .
From a practical point of view, instead of the bilinear incidence rate in 10 , we consider saturation incidence rate in this paper and assume the force of infection is in this version β SI t − τ / 1 αS which is saturated with the susceptible.The susceptible host population is also assumed to have the logistic growth with carrying capacity K, with a specific growth rate constant r.We can get a generalized SIR epidemiological model as follows: where K, r, α, γ, μ 1 , and μ 2 are positive constants.α is the parameter that measures the inhibitory effect, γ is the natural recovery rate of the infective individuals, μ 1 and μ 2 represent the per capita death rates of infectives and recovered, respectively.Notice that when α 0, the system 1.1 becomes the system of bilinear incidence rate in 10 , throughout this paper, we assume α / 0. By mathematical analysis, we derive a threshold value R 0 and prove that the values of R 0 and incubation time length completely determine the global dynamics of system 1.1 , that is, this two factors determine whether the disease approaches an endemic value or whether solutions oscillate.If R 0 ≤ 1, the disease-free equilibrium is globally asymptotically stable and the disease always dies out, whereas if R 0 > 1, the disease persists if it is initially present.By taking the incubation time delay as a bifurcation parameter, the local stability for the endemic equilibrium is investigated, and the conditions with respect to the system to be absolutely stable and conditionally stable are derived.Numerical simulations show that the system with time delay admits rich complex dynamic, and a sequence of periodic solutions will emanate with increasing time delay, which exhibits quite complex periodic and chaotic patterns.
We arrange our paper as follows.In Section 2, results on positivity and boundedness of solutions are presented.In addition, we also consider the equilibria of system 1.1 and give the threshold for the existence of endemic equilibrium.In Section 3, we consider the global stability of the disease-free equilibrium and obtain the necessary and sufficient conditions for the permanence of endemic equilibrium.The local stability analysis of system 1.1 is considered in Section 4. Some numerical results will be given as applications in Section 5.
It can be verified that the positive cone R 3 is positively invariant with respect to 1.1 from 21, Lemma 2.1 .Lemma 2.1.All feasible solutions of the system 1.1 are bounded and enter the region where μ m min 1, μ 1 , μ 2 , lim sup t → ∞ S t ≤ M : max{S 0 , K}.
Proof.From the first equation of 1.

2.4
Thus, we have 0 Therefore, all feasible solutions of the system 1.1 are bounded and enter the region Ω ε .This completes the proof of Lemma 2.1.
Lemma 2.1 shows that the solutions of system 1.1 are bounded and, hence, lie in a compact set and are continuable for all positive time.

Permanence
Before starting our theorem, we give the following lemma.
From the first equation of 1.1 , then there exists a ε > 0 such that S t < K ε for some T 1 > 0 when t ≥ T 1 .Since S/ 1 αS is increasing function with respect to S, then from the second equation of 1.1 , we have Next, we shall consider the case when R 0 1.Noticing that R 0 1 is equal to β K/ 1 αK μ 1 γ.Since S t ≤ r K − S t S t , S t is always decreasing when above K.If S t should ever get below K then S t must stay strictly below K for all subsequent time.This implies there are two possible cases, either i S t → K from above as t → ∞, or ii there exists T such that S t < K for all t > T.
In the first of these cases, we have only to show that I t → 0. Integrating the first equation for S from τ to t τ in 1.1 , we get
In the second of these cases, consider the functional Then, for all t > T τ,

3.7
A direct application of Liapunov-LaSalle type theorem 22 shows that lim t → ∞ I t 0. By the third equation of 1.1 , we get that lim t → ∞ I t 0 implies lim t → ∞ R t 0. This proves R 0 ≤ 1 is the sufficient conditionfor lim t → ∞ S t , I t , R t K, 0, 0 .
According to 1.1 and the definitions for permanence in 23 , we have the following lemma.
Lemma 3.3.Permanence of S t , I t in system 1.1 implies that of R t .
Next, we represent our main results in this section.
In order to prove Theorem 3.4, we present uniform persistence theory for infinite dimensional systems from 24 .Let X be a complete metric space.Suppose that X 0 is open and dense in X and X Let T b t T t | X 0 and let A b be the global attractor for T b t .
Lemma 3.5 see 24 .Suppose that T t satisfies 3.8 and one has the following: iii A b x∈A b ω x is isolated and thus has an acyclic covering M, where Then X 0 is a uniform repeller with respect to X 0 , that is, there is an ε > 0 such that for any Proof of Theorem 3.4.We first prove that R 0 > 1 leads to the permanence of system 1.1 .
By Lemma 3.3, we only need to consider the following subsystem of 1.1 and prove that S t , I t in system 3.10 are permanent if and only if
We begin by showing that the boundary planes of R 2 { x, y : x ≥ 0, y ≥ 0} repel the positive solutions to system 3.10 uniformly.Let C −τ, 0 , R 2 denote the space of continuous functions mapping −τ, 0 into R 2 .We choose Next, we verify that the conditions of Lemma 3.5 are satisfied.By the definition of C 0 and C 0 and system 3.10 , it is easy to see that C 0 and C 0 are positively invariant.Moreover, conditions i and ii of Lemma 3.5 are clearly satisfied.Thus, we only need to verify conditions iii and iv .Since system 3.10 possesses two constant solutions in and we have Ṡ t | ϕ 1 ,ϕ 2 ∈C 1 ≡ 0, then we get S t | ϕ 1 ,ϕ 2 ∈C 1 ≡ 0 for all t ≥ 0, according to the second equation of 3.10 , we have This shows that invariant sets E 0 and E 1 are isolated invariant, then { E 0 , E 1 } is isolated and is an acyclic covering, satisfying condition iii of Lemma 3.5.Now, we show that W s E i ∩ C 0 ∅, i 0, 1.We only need to prove W s E 1 ∩ C 0 ∅, since the proof for W s E 0 ∩ C 0 ∅ is simple.
Assume the contrary, that is, W s E 1 ∩ C 0 / ∅, then there exists a positive solution S t , I t to system 3.10 with lim t → ∞ S t , I t K, 0 .Since R 0 > 1, then for a sufficiently small ε > 0 with μ 1 γ < β K − ε / 1 α K − ε , there exists a positive constant T T ε such that By the second equation of 3.10 and noting that β S/ 1 αS is an increasing function with respect to S, then we have

3.14
According to the comparison principle, lim t → ∞ I t ∞ when R 0 > 1, contradicting I t < ε.Then we have W s E 1 ∩ C 0 ∅.At this time, we are able to conclude from Lemma 3.5 that C 0 repels the positive solutions of 3.10 uniformly.Incorporating the above results into Lemmas 3.3 and 3.5, we know that system 1.1 is permanent.
Next, we verify that permanence of system 1.1 indicates R 0 > 1. Assume that the contrary holds, that is, R 0 ≤ 1, then by Theorem 3.2, S t , I t , R t → K, 0, 0 , as t → ∞, contradicting permanence of system 1.1 .This proves Theorem 3.4.

Linearized analysis
The dynamics of model 1.1 are determined by the first two equations.Therefore, throughout the remainder of this paper, we consider the subsystem 3.10 , and rewrite it as follows:

4.1
Let E S, I be any equilibrium of 4.1 , linearized system of 4.1 at E S, I , we get

4.2
Then the characteristic equation of 4.1 at E is given by det At the equilibrium E 0 0, 0 , characteristic equation 4.3 reduces to λ − r λ μ 1 γ 0.

4.4
Obviously, 4.4 has a positive root λ r independent of any parameters.Hence, E 0 is always a unstable saddle point.
Theorem 4.1.For the system 4.1 , the equilibrium Proof.The characteristic equation at E 1 is Equation 4.5 has a negative real part characteristic root λ −r and roots of When R 0 < 1, there are no positive real roots ω.This shows that all roots of F λ 0 must have negative real parts, therefore, E 1 is an asymptotically stable equilibrium.
ii Assume that R 0 1, then λ 0 is a root of 4.6 .It is easy to verify that λ 0 is a simple characteristic root.If the other roots are λ α iω, then they must satisfy and we must have α ≤ 0. Therefore E 1 is linearly neutrally stable.
iii Assume that R 0 > 1, then F 0 < 0, and F ∞ ∞.Hence, F λ has at least one positive root and E 1 is unstable.
By the arguments to Theorems 3.2 and 4.1, we directly have the following corollary.
In the following, we will study the linear stability of the positive equilibrium E S * , I * of 4.1 .We can see that the characteristic roots of 4.3 at positive equilibrium E are the roots of det Since β I * / 1 αS * r 1 − 1/R 0 and β S * / 1 αS * μ 1 γ at S * , I * , we have P λ, τ Q λ, τ e −λτ 0, 4.10 where

4.11
When τ 0, the DDE 4.1 becomes ODE which has the same equilibria E as follows: 4.12 and 4.10 becomes /αS * , the system 4.12 is locally asymptotically stable.If R 0 > R cc , the unique positive equilibrium of 4.12 is unstable, system 4.12 becomes oscillatory in a stable limit cycle 25 , and this limit cycle is unique 26, 27 .
If λ iω ω > 0 is a root of 4.10 , then by separating the real and imaginary parts, we get

4.14
Squaring and adding both equations, then we have

4.15
, there is no positive real ω satisfying 4.15 , thus eigenvalues of 4.10 do not approach the imaginary axis for any τ > 0. This shows that E is absolutely stable when , there is a unique positive ω 0 satisfying 4.15 .That is, 4.10 has a unique pair of purely imaginary roots ±iω 0 .
Summarizing the discussion above, we have the following conclusion.Theorem 4.3.For system 4.1 , one has , then E is conditionally stable, that is, there is a critical delay value τ 0 such that E is asymptotically stable when τ ∈ 0, τ 0 and unstable when τ > τ 0 .Furthermore, system 4.1 undergoes Hopf bifurcation at E when τ τ n , n 0, 1, 2, . ..; iii if R 0 > R cc , then there is also a critical delay value τ 0 such that the periodic solution is still stable when τ ∈ 0, τ 0 , however, there are a sequence of periodic solutions emanate when τ τ n , n 0, 1, 2, . . . .These results are summarized in Table 1.
Remark 4.4.In fact, in the case of R 0 > R cc , we have known that ODE 4.12 has a stable periodic solution 25 .If we consider the impact of the incubation time on ODE 4.12 , that is, DDE 4.1 , from above discussion, we can see that there are a sequence of periodic solutions bifurcate from the positive equilibrium E when the time delay takes the critical delay τ n such that previous stable periodic solution losses stability, which will lead to complex dynamic phenomena.This can be seen from the simulation results Section 5.
In addition, we want to mention that Theorem 4.3 ii cannot determine the direction and stability of bifurcation periodic solutions, this can be done by analyzing the high-order terms in terms of 28 .The method is very complex and trivial, here we omit it.
From the point of biology, in comparison with the results of 10 , we see that if there is the inhibition effect from the behavioral change of the susceptible individuals when the infective increases i.e., we take use of the saturation incidence rate , then the threshold R c decline and less than 3.It is crucial for the government to take the corresponding control measures and policies against the disease when the epidemic outbreaks.

Numerical results
The main goal of the previous section was to qualitatively characterize the dynamic behaviors of system 1.1 at long term.In this section, we confirm our previous theoretical analysis in Section 4 and demonstrate that the local behaviors in the regions of the parameters space correspond to complex population dynamics to system 1.1 .The objective is to explore the possibility of chaotic behavior in system 1.1 .It is difficult to test whether there exists chaos in a time-delayed system, but numerical simulation analysis is a valid method for such a system.Extensive numerical simulations are carried out for different values of saturation parameter α and recover rate γ.The quality results are as follows.
First, we consider the property of system 1.1 in the regions of the parameter space corresponds to complex population dynamics in the case of R 0 > R cc .To illustrate the transition from the periodic pattern to chaotic pattern, we concentrate on the regions of small and large time delay as an example.We consider the set of parameter values as r 0.1, K 80, β 0.1, μ 1 0.5, μ 2 0.1, α 0.05, and γ 0.5.By calculating, there exists the relation  of R 0 > R cc R 0 4.00, R cc 3.00 for system 1.1 .Then, from Theorem 4.3 iii , the system 1.1 has a stable period solution if the time delay is less than the critical delay τ 0 .0.23 see Figure 1 , and the periodic solution will lost stability when the time delay is greater than the critical delay τ 0 .0.23, and then a typical chaos was observed with increasing the time delay see Figure 2 .This phenomenon has been verified by the bifurcation diagram via delay τ, as shown in Figure 3.We increase the recovery rate γ, let γ 0.75, and fix the other parameters as above.By calculating, the system 1.1 satisfies R c < R 0 < R cc R c 2.23, R 0 2.40, R cc 2.6 .In this context, by Theorem 4.3 ii , system 1.1 is conditionally stable at unique positive equilibrium E 33.33, 1.56, 11.67 , that is, there exists a critical delay τ 0 .0.31 such that E is stable if τ < τ 0 see Figure 4 , E will lose stability by an Hopf bifurcation if τ > τ 0 , as shown in Figure 5.We find that the periodic solution, quasiperiod, and chaos patterns emerge with increasing time delay.This may be clear from the bifurcation diagram see Figure 7 .Now, we fix the parameter γ 0.75, change α, and let α 0.055 i.e., we take some measures to protect on susceptibles , and the other parameters are also as above.By calculating, system 1.1 always exists a relationship 1 < R 0 < R c R 0 2.00, R c 2.19 .By Theorem 4.3 i , we know that system 1.1 is absolutely stable at unique positive equilibrium E 40, 1.60, 12 for any value of the time delay, as shown in Figure 8.Our numerical simulations have demonstrated the validity of our theoretical analysis, that is, the values of threshold value R 0 and incubation time τ length completely determine the dynamics of system 1.1 .
It is necessary to indicate that the system 1.1 is realistic at the initial phase of the disease emergence since the number of the infected is rare.However, when the number of the infected is large, it is more reasonable that one should replace the term r 1 − S/K S in the first equation of 1.1 by r 1 − S I R /K S. Hence, a profound understanding for this case is still desirable and could motivate further investigations.