- Email: [email protected]
Analysis of a SEIV epidemic model with a nonlinear incidence rate
Applied Mathematical Modelling 33 (2009) 2919–2926
Contents lists available at ScienceDirect
Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm
Analysis of a SEIV epidemic model with a nonlinear incidence rate q Li-Ming Cai a,b,*, Xue-Zhi Li a a b
Department of Mathematics, Xinyang Normal University, Xinyang 464000, PR China Beijing Institute of Information Control, Beijing 100037, PR China
a r t i c l e
i n f o
Article history: Received 17 August 2007 Received in revised form 3 January 2008 Accepted 7 January 2008 Available online 15 January 2008
Keywords: Epidemic model Nonlinear incidence rate Uniformly persistence Global stability
a b s t r a c t In this paper, a SEIV epidemic model with a nonlinear incidence rate is investigated. The model exhibits two equilibria, namely, the disease-free equilibrium and the endemic equilibrium. It is shown that if the basic reproduction number R0 < 1, the disease-free equilibrium is globally asymptotically stable and in such a case the endemic equilibrium does not exist. Moreover, we show that if the basic reproduction number R0 > 1, the disease is uniformly persistent and the unique endemic equilibrium of the system with saturation incidence is globally asymptotically stable under certain conditions. Ó 2008 Elsevier Inc. All rights reserved.
1. Introduction Vaccination is a commonly used method for controlling disease, e.g., hepatitis B, measles, or influenza, etc. But it is practically impossible to vaccinate all susceptible individuals in a given community, especially in countries where such vaccines are not easily available or affordable. Some clinical results [1] show that vaccines can only offer a temporary immunity to the diseases. Thus once a vaccine wanes from the body of the vaccinated person, the person becomes susceptible to the disease again. Therefore, it is of vital importance to determine the optimal vaccine coverage needed to eradicate a disease. In modelling of communicable disease, the incidence rate (the rate of new infections) is considered to play a vital role in ensuring that the model can give a reasonable qualitative description of the disease dynamics [2]. Bilinear and standard incidence rate have been frequently used in classical epidemiological models [3]. However, there are several reasons for using nonlinear incidence rates such as saturating and nearly bilinear. For example, Yorke and London [4] showed that an incidence rate gðIÞS ¼ bð1 cIÞIS with positive c and time dependent b is consistent with the results of the simulations for measles outbreaks. To incorporate the effect of behavioral changes, Liu and coworkers [5,6] investigated a nonlinear incidence rate given by bIp S=ð1 þ aIq Þ, with b; a; p; q > 0. Mathematical models with nonlinear incidence rates are numerous in the literature (see [7–14]), and we refer the reader to [3] for a general reference. In this paper, incorporating a general nonlinear incidence rate and a waning preventive vaccines, we consider a fourdimensional model, which consists of the susceptible individuals ðSÞ, exposed individuals but not yet infectious ðEÞ, infectious individuals ðIÞ and vaccinated-treated individuals ðVÞ.
q This work is supported by the National Natural Science Foundation of China (10671166), the Natural Science Foundation of Henan Province (2007110028) and Cadreman youth teacher of Xinyang Normal University. * Corresponding author. Address: Department of Mathematics, Xinyang Normal University, Xinyang 464000, PR China. E-mail address: [email protected] (L.-M. Cai).
0307-904X/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2008.01.005
2920
L.-M. Cai, X.-Z. Li / Applied Mathematical Modelling 33 (2009) 2919–2926
The model equations are given as follows:
dSðtÞ bSðtÞIðtÞ ¼ ð1 pÞP lSðtÞ þ xVðtÞ; dt uðIÞ dEðtÞ bSðtÞIðtÞ ¼ ðl þ rÞEðtÞ; dt uðIÞ dIðtÞ ¼ rEðtÞ ðl þ cÞIðtÞ; dt dVðtÞ ¼ pP lVðtÞ þ cIðtÞ xVðtÞ; dt
ð1Þ
where P is the recruitment rate of individuals (including newborns and immigrants) into the susceptible population; p is the fraction of recruited individuals who are vaccinated; b is the rate at which susceptible individuals become infected by those who are infectious; l is the natural death rate; r is the rate at which exposed individuals become infectious; c is the rate at which infected individuals are treated or recovered; and x is the rate at which vaccine wanes. The nonlinear incidence is p q assumed to be of the form bSðtÞIðtÞ uðIÞ , which generalizes the incidence rate bI S=ð1 þ aI Þ. Obviously it is a bilinear incidence rate _ ðIÞ P 0, which implies uðIÞ P 1 for IðtÞ P 0. 1=uðIÞ may as uðIÞ ¼ 1. We assume that function uðIÞ satisfies uð0Þ ¼ 1 and u measure the psychological or inhibition effect from the behavior change of the susceptible individuals when the number of the infectious individuals increases. This is because the number of effective contacts between infective individuals and susceptible individuals decreases at high infective levels due to the quarantine of infective individuals or due to the protection measures by the susceptible individuals. Summing the equations in (1), we obtain that the total population N satisfies the differential equation
dN=dt ¼ P lN; whose solution is given by the formula
NðtÞ ¼ N0 elt þ
P
l
ð1 elt Þ:
Thus, we assume that the initial value is N 0 ¼ S0 þ E0 þ I0 þ V 0 ¼ P l in order to have a population of constant size (that is, SðtÞ þ EðtÞ þ IðtÞ þ VðtÞ ¼ N P l ). Obviously, the state variables ðSðtÞ; EðtÞ; IðtÞ; VðtÞÞ remain in the biologically meaningful 4 set X ¼ fðS; E; I; VÞ 2 R4þ j 0 6 S þ E þ I þ V 6 P l g for ðSð0Þ; Eð0Þ; Ið0Þ; Vð0ÞÞ 2 Rþ , which is a positively invariant region. Using VðtÞ ¼ P=l SðtÞ EðtÞ IðtÞ to eliminate VðtÞ from the equations in (1) leads to the following reduced threedimensional model:
dSðtÞ bSðtÞIðtÞ P SðtÞ EðtÞ IðtÞ ; ¼ ð1 pÞP lSðtÞ þ x dt uðIÞ l dEðtÞ bSðtÞIðtÞ ¼ ðl þ rÞEðtÞ; dt uðIÞ dIðtÞ ¼ rEðtÞ ðl þ cÞIðtÞ: dt
ð2Þ
Let N 1 ¼ S þ E þ I. From the reduced model (2), we have
dN1 ½lð1 pÞ þ xP ðl þ xÞN1 cI: ¼ dt l From the above equation, it can be seen that, in the absence of the disease ðI ¼ 0Þ, N 1 ! ððlð1 pÞ þ xÞPÞ=ðlðl þ xÞÞ. Since the spread of the disease in the population will reduce N 1 , it follows that N 1 2 ½0; ððlð1 pÞ þ xÞPÞ=ðlðl þ xÞÞ. Noting that X is a positively invariant region for the original model, it can be seen that
X1 ¼ ðS; E; IÞ : S P 0; E P 0; I P 0; S þ E þ I 6
½lð1 pÞ þ xP lðl þ xÞ
is also a positively invariant region for the model (2), and the model (2) is obviously well-posed in X1 . The aim of this paper is to use the model (2) to investigate the global dynamics and to predict the optimal vaccination convergence needed to ensure the disease does not spread. The paper is organized as follows. In the next section, the existence and stability of equilibria is investigated. In Section 3, the disease persistence is discussed. In Section 4, global asymptotic stability of the endemic equilibrium is also investigated. The paper ends with brief remarks. 2. Existence and stability of equilibria Now we investigate the model (2) by finding its equilibria and studying their stability. Steady states of the model (2) satisfy the following equations:
2921
L.-M. Cai, X.-Z. Li / Applied Mathematical Modelling 33 (2009) 2919–2926
ð1 pÞP
bSI P S E I ¼ 0; lS þ x uðIÞ l
bSI
ð3Þ
ðl þ rÞE ¼ 0; uðIÞ rE ðl þ cÞI ¼ 0: It is easy to check that the Eq. (3) always has the disease-free equilibrium P 0 equation of (3), we obtain
E¼
½xþlð1pÞP lðlþxÞ ; 0; 0
. From the second and third
lþc ðl þ rÞðl þ cÞuðIÞ I; S ¼ : br r
ð4Þ
After substituting (4) into the first equation of (3), we obtain the following equation for I:
HðIÞ ¼
½lð1 pÞ þ xP
l
ðl þ r þ xÞðl þ cÞ þ rx
r
I
ðl þ xÞðl þ rÞðl þ cÞ uðIÞ: br
ð5Þ
It can be easily seen that the function HðIÞ is negative for large positive I. Next, we determine the sign of its derivative:
ðl þ r þ xÞðl þ cÞ þ rx ðl þ xÞðl þ rÞðl þ cÞ _ u_ ðIÞ < 0 ðsince u_ ðIÞ P 0Þ: HðIÞ ¼ br r xP Let R0 ¼ lðrlb½þlrð1pÞþ . Since uð0Þ ¼ 1, it follows that ÞðlþcÞðlþxÞ
Hð0Þ ¼
½lð1 pÞ þ xP
l
ðl þ xÞðl þ rÞðl þ cÞ > 0 for R0 > 1: br
ð6Þ
Therefore, if R0 > 1, then the unique positive root of the equation HðIÞ ¼ 0 always exists. By (4), it follows that the model (2) has the unique positive equilibrium P ðS ; E ; I Þ for R0 > 1. In the following, firstly we investigate the stability of the disease-free equilibrium P0 . The linearization of the model (2) at the equilibrium P 0 gives the following characteristic equation:
½x þ lð1 pÞP ¼ 0: ðk þ l þ xÞ k2 þ ð2l þ r þ cÞk þ ðl þ rÞðl þ cÞ rb lðl þ xÞ
ð7Þ
It is obvious that the characteristic Eq. (7) has always a negative eigenvalue k ¼ l x. The other eigenvalues of Eq. (7) are determined by the equation
k2 þ ð2l þ r þ cÞk þ ðl þ rÞðl þ cÞ rb
½x þ lð1 pÞP ¼ 0: lðl þ xÞ
ð8Þ
P It is easy to see that all roots of Eq. (8) have negative real parts if and only if ðl þ rÞðl þ cÞ rb ½xþllðlð1pÞ > 0, i:e:; R0 < 1. If þxÞ R0 ¼ 1, one eigenvalue of Eq. (8) is 0 and it is simple. If R0 > 1, one of roots of Eq. (8) has positive real parts. Thus, we first establish the following lemma.
Lemma 1. If R0 < 1, the disease-free equilibrium P 0 is locally asymptotically stable; If R0 ¼ 1, P0 is stable; If R0 > 1, P 0 is unstable. To obtain the global attractivity of the disease-free equilibrium P0 , we need the following lemma. Let
f1 ¼ lim inf f ðhÞ; t!1 hPt
f 1 ¼ lim sup f ðhÞ: t!1 hPt
Lemma 2 [15]. Assume that a bounded real-valued function f : ½0; 1Þ ! R be twice differentiable with bounded second derivative. Let k ! 1, t k ! 1 and f ðt k Þ converges to f 1 or f1 . Then limk!1 f 0 ðtk Þ ¼ 0. Theorem 1. If R0 < 1, then the disease-free equilibrium P 0 is globally asymptotically stable. Proof. From the above discussion, we have obtained that the unique disease-free equilibrium P 0 of the model (2) is locally _ 6 ð1 pÞP þ xP=l ðl þ xÞS for asymptotically stable whenever R0 < 1. From the first equation of (2), we have SðtÞ uðIÞ P 1. A solution of the equation X_ ¼ ð1 pÞP þ xP=l ðl þ xÞX, is a supper solution of SðtÞ (That is, XðtÞ P SðtÞ for P all t P 0). Noting that XðtÞ ! ½xþllðlð1pÞ as t ! 1, it follows that for a given þxÞ
SðtÞ 6 XðtÞ 6
S1 6
½xþlð1pÞP lðlþxÞ
1
þ for t P t0 . Thus, S
½x þ lð1 pÞP : lðl þ xÞ
6
½xþlð1pÞP lðlþxÞ
þ . Letting
> 0,
there exists a t 0 such that
! 0, we have ð9Þ
2922
L.-M. Cai, X.-Z. Li / Applied Mathematical Modelling 33 (2009) 2919–2926
Consequently, the second equation of (2) can be expressed as
½x þ lð1 pÞP _ EðtÞ 6b þ IðtÞ ðl þ rÞEðtÞ for t P t 0 : l ð l þ xÞ
ð10Þ
Using (10) and the third equation of (2), we have the following system:
E_ I_
! 6Q
E I
with Q ¼
! P þ ðl þ rÞ b ½xþllðlð1pÞ þxÞ
r
ðl þ cÞ
:
ð11Þ
Let M 2 Rþ such that M > maxfl þ r; l þ cg. Thus, Q þ MI22 is a strictly positive matrix where I22 is the identity matrix. It is clear that if x1 and x2 are the eigenvalues of Q, then x1 þ M and x2 þ M are the eigenvalues of Q þ MI22 . It follows from the Perron–Frobenius theorem [16] that Q þ MI22 has a simple positive eigenvalue equal to spectral radius (dominant eigenvalue) and a corresponding eigenvector e > 0 (that is, all components of the eigenvector e are positive). This implies that x1 and x2 are both real. If x1 þ M is the dominant eigenvalue of Q þ MI22 , then x1 > x2 and eQ ¼ x1 e. Obviously, x1 and x2 are the roots of the equation
½x þ lð1 pÞP W 2 þ ð2l þ r þ cÞW þ ðl þ rÞðl þ cÞ rb þ ¼ 0: l ð l þ xÞ Since R0 < 1, for
ð12Þ
> 0 sufficiently small, we have
½x þ lð1 pÞP ðl þ rÞðl þ cÞ rb þ > 0: lðl þ xÞ Therefore, the coefficients of the quadratic equation (12) are positive. Thus, x1 and x2 are both negative. It follows from (11) that for t P t 0 ,
d ðe ½EðtÞ; IðtÞÞ 6 x1 e ½EðtÞ; IðtÞ; dt where ‘‘” denotes the dot product of the two vectors e and ½EðtÞ; IðtÞ. Integrating the above inequality gives
0 6 e ½EðtÞ; IðtÞ 6 e ½Eðt 1 Þ; Iðt1 Þex1 ðtt1 Þ for t P t 1 P t 0 : Since x1 < 0, it follows that e ½EðtÞ; IðtÞ ! 0; as t ! 1. Using e > 0, we conclude that
½EðtÞ; IðtÞ ! ð0; 0Þ as t ! 1:
ð13Þ
½xþlð1pÞP lðlþxÞ .
Now, we show limt!1 SðtÞ ¼ In fact, by Lemma 2, we choose a sequence t n ! 1; sn ! 1ðn ! 1Þ such that _ n Þ ! 0; Sðt _ n Þ ! 0. Noticing that EðtÞ; IðtÞ ! 0 for t ! 1, and from the first equation of system (2), Sðsn Þ ! S1 ; Sðtn Þ ! S1 ; Sðs we obtain
xP ðl þ xÞ lim sup SðtÞ; t!1 l xP ðl þ xÞ lim inf SðtÞ ¼ 0: ð1 pÞP þ t!1 l ð1 pÞP þ
ð14Þ
It follows from (14) that
lim SðtÞ ¼
t!1
½x þ lð1 pÞP : l ð l þ xÞ
Hence, incorporating Lemma 1, the disease-free equilibrium P 0 is globally asymptotically stable if R0 < 1. This completes the proof. h Now we investigate the local stability of the endemic equilibrium P ðS ; E ; I Þ. The linearization of the model (2) about the equilibrium P gives the following characteristic equation:
ðlþrÞE
k þ l þ x þ S
ðlþSr ÞE
0
x
0 x þ ubSðI Þ 1 I uuðIðI Þ Þ
kþlþr
0 ubSðI Þ 1 I uuðIðI Þ Þ
r
kþlþc
¼ 0:
ð15Þ
Thus, we have
k3 þ Q 1 k2 þ Q 2 k þ Q 3 ¼ 0;
ð16Þ
L.-M. Cai, X.-Z. Li / Applied Mathematical Modelling 33 (2009) 2919–2926
2923
where
ðl þ rÞE Q 1 ¼ 3l þ x þ r þ c þ > 0; S ðl þ rÞE rbS I u0 ðI Þ xðl þ rÞE > 0; þ Q2 ¼ l þ x þ ð2l þ r þ cÞ þ S u2 ðI Þ S bS I u0 ðI Þ ðl þ rÞE ðl þ rÞE brðl þ rÞE Q 3 ¼ ðl þ xÞ r ðl þ cÞ þ rx þ þx > 0: 2 u ðI Þ S S uðI Þ After a little algebraic calculation, we have
ðl þ rÞE ðl þ rÞE ðl þ rÞE Q 1 Q 2 Q 3 ¼ ð2l þ r þ cÞ l þ x þ 3l þ x þ r þ c þ þ ðl þ xÞ x S S S 0 0 0 ðl þ rÞE bS I u ðI Þ ðl þ rÞE bS I u ðI Þ bS I u ðI Þ þ lþ r þx þr ðl þ cÞ þ r2 S u2 ðI Þ S u2 ðI Þ uðI Þ uðI Þ bðl þ rÞE ðl þ rÞE ¼ ðl þ xÞð2l þ r þ cÞ 3l þ x þ r þ c þ r uðI Þ S ðl þ rÞE ðl þ rÞE ðl þ rÞE þ ð l þ r Þð2 l þ þ r þ x Þ þ ð l þ c Þ 3 l þ x þ r þ c þ S S S 0 ðl þ rÞE ðl þ rÞE bS I u ðI Þ ðl þ rÞE bS I u0 ðI Þ þ ð l þ xÞ x þ r2 þ lþ r þx 2 S S u ðI Þ S u2 ðI Þ 0 bS I u ðI Þ ðl þ cÞ > 0: þr u2 ðI Þ By the Routh–Hurwitz theorem, it follows that all the roots of the Eq. (16) have negative real parts. Hence, P ðS ; E ; I Þ is locally asymptotically stable. From the above discussion, we can summarize the following conclusion. Theorem 2. If R0 > 1, then system (2) has a unique equilibrium P ðS ; E ; I Þ, which is locally asymptotically stable. 3. Persistence of disease In the following, we shall apply Theorem 4.6 in [15] to study the persistence of the disease, which has the same significance as the global stability of the endemic equilibrium in the epidemiological sense. The persistence of the population has been investigated recently by some authors [17–22]. To keep the same sign as the conditions stated in Theorem 4.6 in [15], we define X 1 ¼ intðR3þ Þ and X 2 ¼ bdðR3þ Þ. Theorem 3. If R0 > 1, then there exists lim inf t!1 XðtÞ > e for X ¼ S; E; I.
e > 0, independent of initial conditions satisfying E0 þ I0 > 0; such that
Proof. In the previous section, we have shown X1 is bounded and positively invariant. Thus, there exists a compact set M0 , in which all solutions of system (2) initiated in R3þ ultimately enter and remain forever after. The compactness conditions (C 4:2 ) of Theorem 4.6 [15] is easily verified for this set M 0 . Let xðx0 Þ be the x-limit set of the solution xðt; x0 Þ of system (2) starting in x0 2 R3þ . We need to show the following set holds:
X2 ¼
[
xðyÞ; where Y ¼ fx0 2 X 2 j xðt; x0 Þ 2 X 2 8t > 0g:
y2Y
From system (2), it follows that all solutions start in bdðR3þ Þ but not on the S axis ultimately leave bdðR3 Þ. This implies that Y ¼ fðS; E; IÞ 2 bdðR3þ Þ j E ¼ I ¼ 0g. Furthermore, we observe that X2 ¼ fP0 g when all solutions initiated on the S axis converge to P0 . Then P 0 is a covering of X2 , which is isolated and acyclic. In the following, if it is verified that P 0 is a weak repeller for X 1 , then the proof of Theorem is complete. By definition, P0 is a weak repeller for X 1 if every solution starts in x0 2 X 1 ,
lim sup dðxðt; x0 Þ; P0 Þ > 0:
ð17Þ
t!1
To show that (17) holds, we verify the following:
W s ðP0 Þ s
\
intðR3þ Þ ¼ ;;
ð18Þ
where W ðP 0 Þ denotes the stable manifold of P0 . To do so, suppose that (17) does not hold for some solution xðt; x0 Þ starting in x0 2 X 1 . Since that the closed positive octant is positively invariant for system (2), it follows that lim inf t!1 dðxðt; x0 Þ; P0 Þ ¼ lim supt!1 dðxðt; x0 Þ; P 0 Þ ¼ 0 and limt!1 xðt; x0 Þ ¼ E0 . This is obviously impossible if (18) holds.
2924
L.-M. Cai, X.-Z. Li / Applied Mathematical Modelling 33 (2009) 2919–2926
We now verify that (18) holds. In fact, since the disease-free equilibrium P0 is unstable for R0 > 1. In particular, the Jacobian matrix of system (2) (JðP0 Þ) possesses one eigenvalue with positive real part, which denotes as kþ , and two eigenvalues with negative real part, which we respectively, denote as k and ðl þ xÞ (Note that k may be equal to ðl þ xÞ.) We shall proceed by determining the location of Es ðP0 Þ (the stable eigenspace of P0 ). Clearly, ð1; 0; 0ÞT is an eigenvector of JðP0 Þ associated to ðl þ xÞ. If k – ðl þ xÞ, then the eigenvector associated to k has the following: ð0; n1 ; n2 ÞT , where n1 and n2 satisfy the eigenvector equation P ðl þ rÞ b ½xþllðlð1pÞ þxÞ
r
!
n1
n2
ðl þ cÞ
¼ k
n1
n2
ð19Þ
:
If k ¼ ðl þ xÞ, then k is a repeated eigenvalue, and an associated generalized eigenvector will have the following structure: ð; n1 ; n2 ÞT , where the value of is irrelevant for what follows and n1 and n2 satisfy (19). In the rest of the proof, if we show that in both cases (i:e:; k – ðl þ xÞ and k ¼ ðl þ xÞ), the vector ðn1 ; n2 ÞT R R2þ , then the proof of Theorem 3 can complete. In fact, by the definition of an irreducible matrix, the matrix in (19) is an irreducible Metzler matrix. By adding a sufficiently large positive multiple of the identify matrix to the matrix in (19) results in a nonnegative irreducible matrix in (19). Thus the conditions of the Perron–Frobenius theorem [16] are satisfied. By the Perron–Frobenius theorem, the matrix in (19) possesses a simple real eigenvalue, which is large than the real part of any other eigenvalue (also called the dominant eigenvalue). Clearly, the dominant eigenvalue here is kþ . But the Perron– Frobenius theorem also implies that every eigenvector does not belong to the closed positive octant since it is not associated T T with the dominant eigenvalue. This means that ðn1 ; n2 ÞT R R2þ . Therefore, Es ðP0 Þ intðR3þ Þ ¼ ;, and thus, W s ðP0 Þ intðR3þ Þ ¼ ;. This completes the proof. h
4. Global stability of the endemic equilibrium In this section, we choose uðIÞ ¼ 1 þ aI; ða > 0Þ for the model (2). Thus, we obtain the following system with saturation incidence:
dSðtÞ bSðtÞIðtÞ P SðtÞ EðtÞ IðtÞ ; ¼ ð1 pÞP lSðtÞ þ x dt 1 þ aI l dEðtÞ bSðtÞIðtÞ ¼ ðl þ rÞEðtÞ; dt 1 þ aI dIðtÞ ¼ rEðtÞ ðl þ cÞIðtÞ: dt
ð20Þ
Obviously, if R0 > 1, the model (20) has a unique endemic equilibrium P . Now we use the geometrical approach of Li and Muldowney in [23] to investigate the global stability of the endemic equilibrium P in the feasible region X1 . Let x#f ðxÞ 2 Rn be a C 1 function for x in an open set X1 Rn : Consider the differential equation
x0 ¼ f ðxÞ:
ð21Þ
Let xðt; x0 Þ denote the solution of (21) satisfying xð0; x0 Þ ¼ x0 . We make the following two assumptions. ðH1 Þ There exists a compact absorbing set K X1 . x in X1 . ðH2 Þ Eq. (21) has a unique equilibrium Lemma 4 (see [23]). Suppose that assumptions ðH1 Þ and ðH2 Þ hold. Assume that (21) satisfies a Bendison criterion, which is x is globally stable in X1 robust under C 1 local perturbations of f at all non-equilibrium nonwandering points for (21). Then provided that it is stable. , then the endemic equilibrium P of ¼ minfr=2; ððbðlð1 pÞ þ xÞ þ xÞpÞ=ðlðl þ xÞÞg. If R0 > 1 and x 6 x Theorem 4. Let x system (20) is global stable in X1 . Proof. By Theorem 2, if R0 > 1, then P is the unique equilibrium in the interior of X1 . Hence the model (2) satisfies the assumption ðH1 Þ. We note that uniform persistence of (20), together with the boundedness of solutions, implies the existence of a compact absorbing set X1 (see [24,25]). This verifies the assumption ðH2 Þ. of of the model (20) is The Jacobian matrix J ¼ ox
0 B J¼B @
1þbIaI l x
x
ð1þbIaIÞ2 x
bI 1þaI
ðl þ rÞ
bI ð1þaIÞ2
0
r
ðl þ cÞ
1 C C A
2925
L.-M. Cai, X.-Z. Li / Applied Mathematical Modelling 33 (2009) 2919–2926
and its second additive compound matrix is
0 B J ½2 ¼ B @
bS ð1þaIÞ2
1
1þbIaI 2l x r
ð1þbSaIÞ2
r
1þbIaI 2l x
x
0
bI 1þaI
2l r c
þx
C C: A 0
0
0
0
Set the function PðxÞ ¼ PðS; E; IÞ ¼ diagf1; EI ; EIg. Then Pf P1 ¼ diagf0; EE II ; EE II g, and the matrix B ¼ Pf P1 þ PJ ½2 P 1 can be written in block form as
B¼
B11
B12
B21
B22
;
where
B11 B21
" ! # bI bS I bS I ; ¼ þx ; 2l x r; B12 ¼ 2 E 2 1 þ aI E ð1 þ aIÞ ð1 þ aIÞ " # 0 rE E0 I bI 2l x c x : ¼ I ; B22 ¼ E I 1þaI bI 0 0 E I 0 I 2l r c 1þaI E 3
Let ðu; v ; wÞ be the vectors in R3 ffi Rð2Þ . We select a norm in R3 as j ðu; v ; wÞ j¼ maxfj u j; j u þ w jg and let measure with respect to this norm. Following the method in [26], we have lðBÞ 6 supfg 1 ; g 2 g, where
g 1 ¼ l1 ðB11 Þþ j B12 j;
l be the Loziskii
g 2 ¼j B21 j þl1 ðB22 Þ;
where j B12 j; j B21 j are matrix norms with respect to the l1 vector norm, and l1 denotes the Loziskii measure with respect to this l1 norm. More specifically,
bI rE ; 2l x r; j B21 j¼ 1 þ aI I ( ! ) bS I bS I ¼ j¼ max þx ; E ð1 þ aIÞ2 E ð1 þ aIÞ2
l1 ðB11 Þ ¼ j B12
bS ð1 þ aIÞ
! þx 2
I : E
To calculate l1 ðB22 Þ, we add the absolute value of the off-diagonal elements to the diagonal one in each column of B22 , and then take the maximum of two sums. This leads to
l1 ðB22 Þ ¼ max
0 E I0 E0 I 0 E0 I 0 2l x c; 2l r c þ x ¼ 2l c þ þ maxfx; x rg: E I E I E I
From the second and the third equation of (20), we have
bSI E0 ¼ þ l þ r; Eð1 þ aIÞ E
rE I
¼
I0 þ l þ c: I
Therefore, we have
g1 ¼
bI bSI xI bI xI þ 6 l x : 2l x r þ þ 2 1 þ aI E 1 þ E a I Eð1 þ aIÞ
ð22Þ
It follows from x 6 r=2 that
g 2 ¼ 2l c þ
E0 I0 rE E0 6 l x þ : þ maxfx; x rg þ I E I E
ð23Þ
Note that it follows from x 6 ðbðlð1 pÞ þ xÞPÞ=ðlðl þ xÞÞ that
xI E
bI 6 0: 1 þ aI
Thus, from (22) and (23), we have
lðBÞ ¼
E0 xI bI E0 ¼ l x: l x þ max 0; E 1 þ aI E E
ð24Þ
Along each solution ðSðtÞ; EðtÞ; IðtÞÞ of system (20) with ðSð0Þ; Eð0Þ; Ið0ÞÞ 2 K, where K is the compact absorbing set. We have
1 t
Z 0
t
1 t
lðBÞds 6 log
EðtÞ ðl þ xÞ; Eð0Þ
2926
L.-M. Cai, X.-Z. Li / Applied Mathematical Modelling 33 (2009) 2919–2926
which implies that
2 ¼ lim sup sup q t!1
x2K
1 t
Z
t
lðBÞds 6
0
lþx 2
< 0:
, then the endemic equilibrium P of the model (20) is globally stable in X1 . This According to Lemma 4, if R0 > 1 and x 6 x completes the proof. h
5. Concluding remarks In this paper, we investigate the dynamical behavior of a SEIV epidemic model that incorporates a nonlinear incidence rate and a waning preventive vaccine. For the model (2), we obtain the basic reproduction number R0 . From the definition of R0 , it can be seen that if
pv ¼ ðl þ xÞ
rbP lðl þ cÞðl þ rÞ ; lrbP
then R0 ¼ 1. Since R0 is a decreasing function of pv , it follows that if p > pv , then R0 < 1; if p < pv , then R0 > 1. Thus, the condition for disease eradication from the community is satisfied if p > pv ; the disease will persist if p < pv . Note that if the vaccinated number p ¼ 0, then R0 ¼ ðrbPÞ=ðlðl þ cÞðl þ rÞÞ: Thus R0 can be rewritten as R0 ¼ ð1 lp=ðl þ xÞÞR0 . The expression for R0 can show vaccination coverage p with a vaccine that induces immunity, which wanes with average duration of protection ð1=xÞ in a population with average life expectancy ð1=lÞ. Therefore, we can obtain the optimal vaccine convergence pv ¼ lþlx ð1 R1 Þ. Clearly this expression also shows vaccination coverage (pv ) is an increasing function of x. 0
Thus, it is necessary and important for public health management to control an epidemic by increasing the duration of the loss of immunity induced by vaccination ð1=xÞ, which reduces the threshold vaccination convergence (pv ). Acknowledgements The authors would like to thank Prof. M. Cross and the anonymous referees for their useful feedback and helpful comments and suggestions which have improved the manuscript. The authors would also like to thank Prof. Michael Y. Li for his help. References [1] M.A. Teitelbaulm, M. Edmunds, Immunization and vaccine-preventable illness, Unites States, 1992–1997, Stat. Bull. Metrop. Insur. Co. 80 (2) (1999) 13–20. [2] S.A. Levin, T.G. Hallam, L.J. Gross, Applied Mathematical Ecology, Springer, New York, 1990. [3] Z. Ma, Y. Zhou, W. Wang, Z. Jin, Mathematical Models and Dynamics of Infectious Diseases, China Sciences Press, Beijing, 2004. [4] J.A. Yorke, W.P. London, Recurrent outbreaks of measles, chickenpox and mumps II, Am. J. Epidemiol. 98 (1973) 469–482. [5] W.M. Liu, H.W. Hethcote, S.A. Levin, Dynamical behavior of epidemiological models with nonlinear incidence rates, J. Math. Biol. 25 (1987) 359–380. [6] W.M. Liu, S.A. Levin, Y. Iwasa, Influence of nonlinear incidence rates upon the behavior of SIRS epidemiological models, J. Math. Biol. 23 (1986) 187– 204. [7] V. Capasso, G. Serio, A generalization of the Kermack–McKendrick deterministic epidemic model, Math. Biosci. 42 (1978) 43–61. [8] S. Ruan, W. Wang, Dynamical behavior of an epidemic model with a nonlinear incidence rate, J. Differen. Equat. 188 (2003) 135–163. [9] D. Xiao, S. Ruan, Global analysis of an epidemic model with nonmonotone incidence rate, Math. Biosci. 208 (2007) 419–429. [10] T. Zhang, Z. Teng, Pulse vaccination delayed SEIRS epidemic model with saturation incidence, Appl. Math. Model. 32 (2008) 1403–1416. [11] H.W. Hethcote, The mathematics of infectious diseases, SIAM Rev. 42 (2000) 599–653. [12] P. van den Driessche, J. Watmough, A simple SIS epidemic model with a backward bifurcation, J. Math. Biol. 40 (2000) 525–540. [13] M.E. Alexander, S.M. Moghadas, Periodicity in an epidemic model with a generalized nonlinear incidence, Math. Biosci. 189 (2004) 75–96. [14] M.E. Alexander, S.M. Moghadas, Bifurcation analysis of an SIRS Epidemic model with generalized incidence, SIAM J. Appl. Math. 65 (2005) 1794–1816. [15] R.H. Thieme, Persistence under relaxed point-dissipativity (with application to an endemic model), SIAM J. Math. Anal. 24 (1993) 407–435. [16] R.S. Varga, Matrix Iterative Analysis, Prentice-Hall, Inc., Englewood Cliffs, NJ, 1962. [17] A. Fonda, Uniformly persistent semidynamical systems, Proc. Amer. Math. Soc. 104 (1988) 111–116. [18] W. Ma, M. Song, Y. Takeuchi, Global stability of an SIR epidemic model with time delay, Appl. Math. Lett. 17 (2004) 1141–1145. [19] P. De Leenheer, H.L. Smith, Virus dynamics: a global analysis, SIAM J. Appl. Math. 63 (2003) 1313–1327. [20] S. Liu, E. Beretta, A stage-structured predator–prey model of Beddington–DeAngelis type, SIAM J. Appl. Math. 66 (2006) 1101–1129. [21] R. Xu, M.A.J. Chaplain, F.A. Davidson, Persistence and global stability of a ratio-dependent predator–prey model with stage structure, Appl. Math. Comput. 158 (2004) 729–744. [22] L. Cai, X. Song, Permanence and stability of a predator–prey system with stage-structure for predator, J. Comput. Appl. Math. 201 (2007) 356–366. [23] M.Y. Li, J.S. Muldowney, A geometric approach to the global-stability problems, SIAM J. Math. Anal. 27 (1996) 1070–1083. [24] L. Wang, M.Y. Li, Mathematical analysis of the global dynamics of a model for HIV infection of CD4þ T cells, Math. Biosci. 200 (2006) 44–57. [25] G.J. Butler, P. Waltman, Persistence in dynamical systems, Proc. Am. Math. Soc. 96 (1986) 425–430. [26] R.H. Martin Jr., Logarithmic norms and projections applied to linear differential systems, J. Math. Anal. Appl. 45 (1974) 432–454.
Contents lists available at ScienceDirect
Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm
Analysis of a SEIV epidemic model with a nonlinear incidence rate q Li-Ming Cai a,b,*, Xue-Zhi Li a a b
Department of Mathematics, Xinyang Normal University, Xinyang 464000, PR China Beijing Institute of Information Control, Beijing 100037, PR China
a r t i c l e
i n f o
Article history: Received 17 August 2007 Received in revised form 3 January 2008 Accepted 7 January 2008 Available online 15 January 2008
Keywords: Epidemic model Nonlinear incidence rate Uniformly persistence Global stability
a b s t r a c t In this paper, a SEIV epidemic model with a nonlinear incidence rate is investigated. The model exhibits two equilibria, namely, the disease-free equilibrium and the endemic equilibrium. It is shown that if the basic reproduction number R0 < 1, the disease-free equilibrium is globally asymptotically stable and in such a case the endemic equilibrium does not exist. Moreover, we show that if the basic reproduction number R0 > 1, the disease is uniformly persistent and the unique endemic equilibrium of the system with saturation incidence is globally asymptotically stable under certain conditions. Ó 2008 Elsevier Inc. All rights reserved.
1. Introduction Vaccination is a commonly used method for controlling disease, e.g., hepatitis B, measles, or influenza, etc. But it is practically impossible to vaccinate all susceptible individuals in a given community, especially in countries where such vaccines are not easily available or affordable. Some clinical results [1] show that vaccines can only offer a temporary immunity to the diseases. Thus once a vaccine wanes from the body of the vaccinated person, the person becomes susceptible to the disease again. Therefore, it is of vital importance to determine the optimal vaccine coverage needed to eradicate a disease. In modelling of communicable disease, the incidence rate (the rate of new infections) is considered to play a vital role in ensuring that the model can give a reasonable qualitative description of the disease dynamics [2]. Bilinear and standard incidence rate have been frequently used in classical epidemiological models [3]. However, there are several reasons for using nonlinear incidence rates such as saturating and nearly bilinear. For example, Yorke and London [4] showed that an incidence rate gðIÞS ¼ bð1 cIÞIS with positive c and time dependent b is consistent with the results of the simulations for measles outbreaks. To incorporate the effect of behavioral changes, Liu and coworkers [5,6] investigated a nonlinear incidence rate given by bIp S=ð1 þ aIq Þ, with b; a; p; q > 0. Mathematical models with nonlinear incidence rates are numerous in the literature (see [7–14]), and we refer the reader to [3] for a general reference. In this paper, incorporating a general nonlinear incidence rate and a waning preventive vaccines, we consider a fourdimensional model, which consists of the susceptible individuals ðSÞ, exposed individuals but not yet infectious ðEÞ, infectious individuals ðIÞ and vaccinated-treated individuals ðVÞ.
q This work is supported by the National Natural Science Foundation of China (10671166), the Natural Science Foundation of Henan Province (2007110028) and Cadreman youth teacher of Xinyang Normal University. * Corresponding author. Address: Department of Mathematics, Xinyang Normal University, Xinyang 464000, PR China. E-mail address: [email protected] (L.-M. Cai).
0307-904X/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.apm.2008.01.005
2920
L.-M. Cai, X.-Z. Li / Applied Mathematical Modelling 33 (2009) 2919–2926
The model equations are given as follows:
dSðtÞ bSðtÞIðtÞ ¼ ð1 pÞP lSðtÞ þ xVðtÞ; dt uðIÞ dEðtÞ bSðtÞIðtÞ ¼ ðl þ rÞEðtÞ; dt uðIÞ dIðtÞ ¼ rEðtÞ ðl þ cÞIðtÞ; dt dVðtÞ ¼ pP lVðtÞ þ cIðtÞ xVðtÞ; dt
ð1Þ
where P is the recruitment rate of individuals (including newborns and immigrants) into the susceptible population; p is the fraction of recruited individuals who are vaccinated; b is the rate at which susceptible individuals become infected by those who are infectious; l is the natural death rate; r is the rate at which exposed individuals become infectious; c is the rate at which infected individuals are treated or recovered; and x is the rate at which vaccine wanes. The nonlinear incidence is p q assumed to be of the form bSðtÞIðtÞ uðIÞ , which generalizes the incidence rate bI S=ð1 þ aI Þ. Obviously it is a bilinear incidence rate _ ðIÞ P 0, which implies uðIÞ P 1 for IðtÞ P 0. 1=uðIÞ may as uðIÞ ¼ 1. We assume that function uðIÞ satisfies uð0Þ ¼ 1 and u measure the psychological or inhibition effect from the behavior change of the susceptible individuals when the number of the infectious individuals increases. This is because the number of effective contacts between infective individuals and susceptible individuals decreases at high infective levels due to the quarantine of infective individuals or due to the protection measures by the susceptible individuals. Summing the equations in (1), we obtain that the total population N satisfies the differential equation
dN=dt ¼ P lN; whose solution is given by the formula
NðtÞ ¼ N0 elt þ
P
l
ð1 elt Þ:
Thus, we assume that the initial value is N 0 ¼ S0 þ E0 þ I0 þ V 0 ¼ P l in order to have a population of constant size (that is, SðtÞ þ EðtÞ þ IðtÞ þ VðtÞ ¼ N P l ). Obviously, the state variables ðSðtÞ; EðtÞ; IðtÞ; VðtÞÞ remain in the biologically meaningful 4 set X ¼ fðS; E; I; VÞ 2 R4þ j 0 6 S þ E þ I þ V 6 P l g for ðSð0Þ; Eð0Þ; Ið0Þ; Vð0ÞÞ 2 Rþ , which is a positively invariant region. Using VðtÞ ¼ P=l SðtÞ EðtÞ IðtÞ to eliminate VðtÞ from the equations in (1) leads to the following reduced threedimensional model:
dSðtÞ bSðtÞIðtÞ P SðtÞ EðtÞ IðtÞ ; ¼ ð1 pÞP lSðtÞ þ x dt uðIÞ l dEðtÞ bSðtÞIðtÞ ¼ ðl þ rÞEðtÞ; dt uðIÞ dIðtÞ ¼ rEðtÞ ðl þ cÞIðtÞ: dt
ð2Þ
Let N 1 ¼ S þ E þ I. From the reduced model (2), we have
dN1 ½lð1 pÞ þ xP ðl þ xÞN1 cI: ¼ dt l From the above equation, it can be seen that, in the absence of the disease ðI ¼ 0Þ, N 1 ! ððlð1 pÞ þ xÞPÞ=ðlðl þ xÞÞ. Since the spread of the disease in the population will reduce N 1 , it follows that N 1 2 ½0; ððlð1 pÞ þ xÞPÞ=ðlðl þ xÞÞ. Noting that X is a positively invariant region for the original model, it can be seen that
X1 ¼ ðS; E; IÞ : S P 0; E P 0; I P 0; S þ E þ I 6
½lð1 pÞ þ xP lðl þ xÞ
is also a positively invariant region for the model (2), and the model (2) is obviously well-posed in X1 . The aim of this paper is to use the model (2) to investigate the global dynamics and to predict the optimal vaccination convergence needed to ensure the disease does not spread. The paper is organized as follows. In the next section, the existence and stability of equilibria is investigated. In Section 3, the disease persistence is discussed. In Section 4, global asymptotic stability of the endemic equilibrium is also investigated. The paper ends with brief remarks. 2. Existence and stability of equilibria Now we investigate the model (2) by finding its equilibria and studying their stability. Steady states of the model (2) satisfy the following equations:
2921
L.-M. Cai, X.-Z. Li / Applied Mathematical Modelling 33 (2009) 2919–2926
ð1 pÞP
bSI P S E I ¼ 0; lS þ x uðIÞ l
bSI
ð3Þ
ðl þ rÞE ¼ 0; uðIÞ rE ðl þ cÞI ¼ 0: It is easy to check that the Eq. (3) always has the disease-free equilibrium P 0 equation of (3), we obtain
E¼
½xþlð1pÞP lðlþxÞ ; 0; 0
. From the second and third
lþc ðl þ rÞðl þ cÞuðIÞ I; S ¼ : br r
ð4Þ
After substituting (4) into the first equation of (3), we obtain the following equation for I:
HðIÞ ¼
½lð1 pÞ þ xP
l
ðl þ r þ xÞðl þ cÞ þ rx
r
I
ðl þ xÞðl þ rÞðl þ cÞ uðIÞ: br
ð5Þ
It can be easily seen that the function HðIÞ is negative for large positive I. Next, we determine the sign of its derivative:
ðl þ r þ xÞðl þ cÞ þ rx ðl þ xÞðl þ rÞðl þ cÞ _ u_ ðIÞ < 0 ðsince u_ ðIÞ P 0Þ: HðIÞ ¼ br r xP Let R0 ¼ lðrlb½þlrð1pÞþ . Since uð0Þ ¼ 1, it follows that ÞðlþcÞðlþxÞ
Hð0Þ ¼
½lð1 pÞ þ xP
l
ðl þ xÞðl þ rÞðl þ cÞ > 0 for R0 > 1: br
ð6Þ
Therefore, if R0 > 1, then the unique positive root of the equation HðIÞ ¼ 0 always exists. By (4), it follows that the model (2) has the unique positive equilibrium P ðS ; E ; I Þ for R0 > 1. In the following, firstly we investigate the stability of the disease-free equilibrium P0 . The linearization of the model (2) at the equilibrium P 0 gives the following characteristic equation:
½x þ lð1 pÞP ¼ 0: ðk þ l þ xÞ k2 þ ð2l þ r þ cÞk þ ðl þ rÞðl þ cÞ rb lðl þ xÞ
ð7Þ
It is obvious that the characteristic Eq. (7) has always a negative eigenvalue k ¼ l x. The other eigenvalues of Eq. (7) are determined by the equation
k2 þ ð2l þ r þ cÞk þ ðl þ rÞðl þ cÞ rb
½x þ lð1 pÞP ¼ 0: lðl þ xÞ
ð8Þ
P It is easy to see that all roots of Eq. (8) have negative real parts if and only if ðl þ rÞðl þ cÞ rb ½xþllðlð1pÞ > 0, i:e:; R0 < 1. If þxÞ R0 ¼ 1, one eigenvalue of Eq. (8) is 0 and it is simple. If R0 > 1, one of roots of Eq. (8) has positive real parts. Thus, we first establish the following lemma.
Lemma 1. If R0 < 1, the disease-free equilibrium P 0 is locally asymptotically stable; If R0 ¼ 1, P0 is stable; If R0 > 1, P 0 is unstable. To obtain the global attractivity of the disease-free equilibrium P0 , we need the following lemma. Let
f1 ¼ lim inf f ðhÞ; t!1 hPt
f 1 ¼ lim sup f ðhÞ: t!1 hPt
Lemma 2 [15]. Assume that a bounded real-valued function f : ½0; 1Þ ! R be twice differentiable with bounded second derivative. Let k ! 1, t k ! 1 and f ðt k Þ converges to f 1 or f1 . Then limk!1 f 0 ðtk Þ ¼ 0. Theorem 1. If R0 < 1, then the disease-free equilibrium P 0 is globally asymptotically stable. Proof. From the above discussion, we have obtained that the unique disease-free equilibrium P 0 of the model (2) is locally _ 6 ð1 pÞP þ xP=l ðl þ xÞS for asymptotically stable whenever R0 < 1. From the first equation of (2), we have SðtÞ uðIÞ P 1. A solution of the equation X_ ¼ ð1 pÞP þ xP=l ðl þ xÞX, is a supper solution of SðtÞ (That is, XðtÞ P SðtÞ for P all t P 0). Noting that XðtÞ ! ½xþllðlð1pÞ as t ! 1, it follows that for a given þxÞ
SðtÞ 6 XðtÞ 6
S1 6
½xþlð1pÞP lðlþxÞ
1
þ for t P t0 . Thus, S
½x þ lð1 pÞP : lðl þ xÞ
6
½xþlð1pÞP lðlþxÞ
þ . Letting
> 0,
there exists a t 0 such that
! 0, we have ð9Þ
2922
L.-M. Cai, X.-Z. Li / Applied Mathematical Modelling 33 (2009) 2919–2926
Consequently, the second equation of (2) can be expressed as
½x þ lð1 pÞP _ EðtÞ 6b þ IðtÞ ðl þ rÞEðtÞ for t P t 0 : l ð l þ xÞ
ð10Þ
Using (10) and the third equation of (2), we have the following system:
E_ I_
! 6Q
E I
with Q ¼
! P þ ðl þ rÞ b ½xþllðlð1pÞ þxÞ
r
ðl þ cÞ
:
ð11Þ
Let M 2 Rþ such that M > maxfl þ r; l þ cg. Thus, Q þ MI22 is a strictly positive matrix where I22 is the identity matrix. It is clear that if x1 and x2 are the eigenvalues of Q, then x1 þ M and x2 þ M are the eigenvalues of Q þ MI22 . It follows from the Perron–Frobenius theorem [16] that Q þ MI22 has a simple positive eigenvalue equal to spectral radius (dominant eigenvalue) and a corresponding eigenvector e > 0 (that is, all components of the eigenvector e are positive). This implies that x1 and x2 are both real. If x1 þ M is the dominant eigenvalue of Q þ MI22 , then x1 > x2 and eQ ¼ x1 e. Obviously, x1 and x2 are the roots of the equation
½x þ lð1 pÞP W 2 þ ð2l þ r þ cÞW þ ðl þ rÞðl þ cÞ rb þ ¼ 0: l ð l þ xÞ Since R0 < 1, for
ð12Þ
> 0 sufficiently small, we have
½x þ lð1 pÞP ðl þ rÞðl þ cÞ rb þ > 0: lðl þ xÞ Therefore, the coefficients of the quadratic equation (12) are positive. Thus, x1 and x2 are both negative. It follows from (11) that for t P t 0 ,
d ðe ½EðtÞ; IðtÞÞ 6 x1 e ½EðtÞ; IðtÞ; dt where ‘‘” denotes the dot product of the two vectors e and ½EðtÞ; IðtÞ. Integrating the above inequality gives
0 6 e ½EðtÞ; IðtÞ 6 e ½Eðt 1 Þ; Iðt1 Þex1 ðtt1 Þ for t P t 1 P t 0 : Since x1 < 0, it follows that e ½EðtÞ; IðtÞ ! 0; as t ! 1. Using e > 0, we conclude that
½EðtÞ; IðtÞ ! ð0; 0Þ as t ! 1:
ð13Þ
½xþlð1pÞP lðlþxÞ .
Now, we show limt!1 SðtÞ ¼ In fact, by Lemma 2, we choose a sequence t n ! 1; sn ! 1ðn ! 1Þ such that _ n Þ ! 0; Sðt _ n Þ ! 0. Noticing that EðtÞ; IðtÞ ! 0 for t ! 1, and from the first equation of system (2), Sðsn Þ ! S1 ; Sðtn Þ ! S1 ; Sðs we obtain
xP ðl þ xÞ lim sup SðtÞ; t!1 l xP ðl þ xÞ lim inf SðtÞ ¼ 0: ð1 pÞP þ t!1 l ð1 pÞP þ
ð14Þ
It follows from (14) that
lim SðtÞ ¼
t!1
½x þ lð1 pÞP : l ð l þ xÞ
Hence, incorporating Lemma 1, the disease-free equilibrium P 0 is globally asymptotically stable if R0 < 1. This completes the proof. h Now we investigate the local stability of the endemic equilibrium P ðS ; E ; I Þ. The linearization of the model (2) about the equilibrium P gives the following characteristic equation:
ðlþrÞE
k þ l þ x þ S
ðlþSr ÞE
0
x
0 x þ ubSðI Þ 1 I uuðIðI Þ Þ
kþlþr
0 ubSðI Þ 1 I uuðIðI Þ Þ
r
kþlþc
¼ 0:
ð15Þ
Thus, we have
k3 þ Q 1 k2 þ Q 2 k þ Q 3 ¼ 0;
ð16Þ
L.-M. Cai, X.-Z. Li / Applied Mathematical Modelling 33 (2009) 2919–2926
2923
where
ðl þ rÞE Q 1 ¼ 3l þ x þ r þ c þ > 0; S ðl þ rÞE rbS I u0 ðI Þ xðl þ rÞE > 0; þ Q2 ¼ l þ x þ ð2l þ r þ cÞ þ S u2 ðI Þ S bS I u0 ðI Þ ðl þ rÞE ðl þ rÞE brðl þ rÞE Q 3 ¼ ðl þ xÞ r ðl þ cÞ þ rx þ þx > 0: 2 u ðI Þ S S uðI Þ After a little algebraic calculation, we have
ðl þ rÞE ðl þ rÞE ðl þ rÞE Q 1 Q 2 Q 3 ¼ ð2l þ r þ cÞ l þ x þ 3l þ x þ r þ c þ þ ðl þ xÞ x S S S 0 0 0 ðl þ rÞE bS I u ðI Þ ðl þ rÞE bS I u ðI Þ bS I u ðI Þ þ lþ r þx þr ðl þ cÞ þ r2 S u2 ðI Þ S u2 ðI Þ uðI Þ uðI Þ bðl þ rÞE ðl þ rÞE ¼ ðl þ xÞð2l þ r þ cÞ 3l þ x þ r þ c þ r uðI Þ S ðl þ rÞE ðl þ rÞE ðl þ rÞE þ ð l þ r Þð2 l þ þ r þ x Þ þ ð l þ c Þ 3 l þ x þ r þ c þ S S S 0 ðl þ rÞE ðl þ rÞE bS I u ðI Þ ðl þ rÞE bS I u0 ðI Þ þ ð l þ xÞ x þ r2 þ lþ r þx 2 S S u ðI Þ S u2 ðI Þ 0 bS I u ðI Þ ðl þ cÞ > 0: þr u2 ðI Þ By the Routh–Hurwitz theorem, it follows that all the roots of the Eq. (16) have negative real parts. Hence, P ðS ; E ; I Þ is locally asymptotically stable. From the above discussion, we can summarize the following conclusion. Theorem 2. If R0 > 1, then system (2) has a unique equilibrium P ðS ; E ; I Þ, which is locally asymptotically stable. 3. Persistence of disease In the following, we shall apply Theorem 4.6 in [15] to study the persistence of the disease, which has the same significance as the global stability of the endemic equilibrium in the epidemiological sense. The persistence of the population has been investigated recently by some authors [17–22]. To keep the same sign as the conditions stated in Theorem 4.6 in [15], we define X 1 ¼ intðR3þ Þ and X 2 ¼ bdðR3þ Þ. Theorem 3. If R0 > 1, then there exists lim inf t!1 XðtÞ > e for X ¼ S; E; I.
e > 0, independent of initial conditions satisfying E0 þ I0 > 0; such that
Proof. In the previous section, we have shown X1 is bounded and positively invariant. Thus, there exists a compact set M0 , in which all solutions of system (2) initiated in R3þ ultimately enter and remain forever after. The compactness conditions (C 4:2 ) of Theorem 4.6 [15] is easily verified for this set M 0 . Let xðx0 Þ be the x-limit set of the solution xðt; x0 Þ of system (2) starting in x0 2 R3þ . We need to show the following set holds:
X2 ¼
[
xðyÞ; where Y ¼ fx0 2 X 2 j xðt; x0 Þ 2 X 2 8t > 0g:
y2Y
From system (2), it follows that all solutions start in bdðR3þ Þ but not on the S axis ultimately leave bdðR3 Þ. This implies that Y ¼ fðS; E; IÞ 2 bdðR3þ Þ j E ¼ I ¼ 0g. Furthermore, we observe that X2 ¼ fP0 g when all solutions initiated on the S axis converge to P0 . Then P 0 is a covering of X2 , which is isolated and acyclic. In the following, if it is verified that P 0 is a weak repeller for X 1 , then the proof of Theorem is complete. By definition, P0 is a weak repeller for X 1 if every solution starts in x0 2 X 1 ,
lim sup dðxðt; x0 Þ; P0 Þ > 0:
ð17Þ
t!1
To show that (17) holds, we verify the following:
W s ðP0 Þ s
\
intðR3þ Þ ¼ ;;
ð18Þ
where W ðP 0 Þ denotes the stable manifold of P0 . To do so, suppose that (17) does not hold for some solution xðt; x0 Þ starting in x0 2 X 1 . Since that the closed positive octant is positively invariant for system (2), it follows that lim inf t!1 dðxðt; x0 Þ; P0 Þ ¼ lim supt!1 dðxðt; x0 Þ; P 0 Þ ¼ 0 and limt!1 xðt; x0 Þ ¼ E0 . This is obviously impossible if (18) holds.
2924
L.-M. Cai, X.-Z. Li / Applied Mathematical Modelling 33 (2009) 2919–2926
We now verify that (18) holds. In fact, since the disease-free equilibrium P0 is unstable for R0 > 1. In particular, the Jacobian matrix of system (2) (JðP0 Þ) possesses one eigenvalue with positive real part, which denotes as kþ , and two eigenvalues with negative real part, which we respectively, denote as k and ðl þ xÞ (Note that k may be equal to ðl þ xÞ.) We shall proceed by determining the location of Es ðP0 Þ (the stable eigenspace of P0 ). Clearly, ð1; 0; 0ÞT is an eigenvector of JðP0 Þ associated to ðl þ xÞ. If k – ðl þ xÞ, then the eigenvector associated to k has the following: ð0; n1 ; n2 ÞT , where n1 and n2 satisfy the eigenvector equation P ðl þ rÞ b ½xþllðlð1pÞ þxÞ
r
!
n1
n2
ðl þ cÞ
¼ k
n1
n2
ð19Þ
:
If k ¼ ðl þ xÞ, then k is a repeated eigenvalue, and an associated generalized eigenvector will have the following structure: ð; n1 ; n2 ÞT , where the value of is irrelevant for what follows and n1 and n2 satisfy (19). In the rest of the proof, if we show that in both cases (i:e:; k – ðl þ xÞ and k ¼ ðl þ xÞ), the vector ðn1 ; n2 ÞT R R2þ , then the proof of Theorem 3 can complete. In fact, by the definition of an irreducible matrix, the matrix in (19) is an irreducible Metzler matrix. By adding a sufficiently large positive multiple of the identify matrix to the matrix in (19) results in a nonnegative irreducible matrix in (19). Thus the conditions of the Perron–Frobenius theorem [16] are satisfied. By the Perron–Frobenius theorem, the matrix in (19) possesses a simple real eigenvalue, which is large than the real part of any other eigenvalue (also called the dominant eigenvalue). Clearly, the dominant eigenvalue here is kþ . But the Perron– Frobenius theorem also implies that every eigenvector does not belong to the closed positive octant since it is not associated T T with the dominant eigenvalue. This means that ðn1 ; n2 ÞT R R2þ . Therefore, Es ðP0 Þ intðR3þ Þ ¼ ;, and thus, W s ðP0 Þ intðR3þ Þ ¼ ;. This completes the proof. h
4. Global stability of the endemic equilibrium In this section, we choose uðIÞ ¼ 1 þ aI; ða > 0Þ for the model (2). Thus, we obtain the following system with saturation incidence:
dSðtÞ bSðtÞIðtÞ P SðtÞ EðtÞ IðtÞ ; ¼ ð1 pÞP lSðtÞ þ x dt 1 þ aI l dEðtÞ bSðtÞIðtÞ ¼ ðl þ rÞEðtÞ; dt 1 þ aI dIðtÞ ¼ rEðtÞ ðl þ cÞIðtÞ: dt
ð20Þ
Obviously, if R0 > 1, the model (20) has a unique endemic equilibrium P . Now we use the geometrical approach of Li and Muldowney in [23] to investigate the global stability of the endemic equilibrium P in the feasible region X1 . Let x#f ðxÞ 2 Rn be a C 1 function for x in an open set X1 Rn : Consider the differential equation
x0 ¼ f ðxÞ:
ð21Þ
Let xðt; x0 Þ denote the solution of (21) satisfying xð0; x0 Þ ¼ x0 . We make the following two assumptions. ðH1 Þ There exists a compact absorbing set K X1 . x in X1 . ðH2 Þ Eq. (21) has a unique equilibrium Lemma 4 (see [23]). Suppose that assumptions ðH1 Þ and ðH2 Þ hold. Assume that (21) satisfies a Bendison criterion, which is x is globally stable in X1 robust under C 1 local perturbations of f at all non-equilibrium nonwandering points for (21). Then provided that it is stable. , then the endemic equilibrium P of ¼ minfr=2; ððbðlð1 pÞ þ xÞ þ xÞpÞ=ðlðl þ xÞÞg. If R0 > 1 and x 6 x Theorem 4. Let x system (20) is global stable in X1 . Proof. By Theorem 2, if R0 > 1, then P is the unique equilibrium in the interior of X1 . Hence the model (2) satisfies the assumption ðH1 Þ. We note that uniform persistence of (20), together with the boundedness of solutions, implies the existence of a compact absorbing set X1 (see [24,25]). This verifies the assumption ðH2 Þ. of of the model (20) is The Jacobian matrix J ¼ ox
0 B J¼B @
1þbIaI l x
x
ð1þbIaIÞ2 x
bI 1þaI
ðl þ rÞ
bI ð1þaIÞ2
0
r
ðl þ cÞ
1 C C A
2925
L.-M. Cai, X.-Z. Li / Applied Mathematical Modelling 33 (2009) 2919–2926
and its second additive compound matrix is
0 B J ½2 ¼ B @
bS ð1þaIÞ2
1
1þbIaI 2l x r
ð1þbSaIÞ2
r
1þbIaI 2l x
x
0
bI 1þaI
2l r c
þx
C C: A 0
0
0
0
Set the function PðxÞ ¼ PðS; E; IÞ ¼ diagf1; EI ; EIg. Then Pf P1 ¼ diagf0; EE II ; EE II g, and the matrix B ¼ Pf P1 þ PJ ½2 P 1 can be written in block form as
B¼
B11
B12
B21
B22
;
where
B11 B21
" ! # bI bS I bS I ; ¼ þx ; 2l x r; B12 ¼ 2 E 2 1 þ aI E ð1 þ aIÞ ð1 þ aIÞ " # 0 rE E0 I bI 2l x c x : ¼ I ; B22 ¼ E I 1þaI bI 0 0 E I 0 I 2l r c 1þaI E 3
Let ðu; v ; wÞ be the vectors in R3 ffi Rð2Þ . We select a norm in R3 as j ðu; v ; wÞ j¼ maxfj u j; j u þ w jg and let measure with respect to this norm. Following the method in [26], we have lðBÞ 6 supfg 1 ; g 2 g, where
g 1 ¼ l1 ðB11 Þþ j B12 j;
l be the Loziskii
g 2 ¼j B21 j þl1 ðB22 Þ;
where j B12 j; j B21 j are matrix norms with respect to the l1 vector norm, and l1 denotes the Loziskii measure with respect to this l1 norm. More specifically,
bI rE ; 2l x r; j B21 j¼ 1 þ aI I ( ! ) bS I bS I ¼ j¼ max þx ; E ð1 þ aIÞ2 E ð1 þ aIÞ2
l1 ðB11 Þ ¼ j B12
bS ð1 þ aIÞ
! þx 2
I : E
To calculate l1 ðB22 Þ, we add the absolute value of the off-diagonal elements to the diagonal one in each column of B22 , and then take the maximum of two sums. This leads to
l1 ðB22 Þ ¼ max
0 E I0 E0 I 0 E0 I 0 2l x c; 2l r c þ x ¼ 2l c þ þ maxfx; x rg: E I E I E I
From the second and the third equation of (20), we have
bSI E0 ¼ þ l þ r; Eð1 þ aIÞ E
rE I
¼
I0 þ l þ c: I
Therefore, we have
g1 ¼
bI bSI xI bI xI þ 6 l x : 2l x r þ þ 2 1 þ aI E 1 þ E a I Eð1 þ aIÞ
ð22Þ
It follows from x 6 r=2 that
g 2 ¼ 2l c þ
E0 I0 rE E0 6 l x þ : þ maxfx; x rg þ I E I E
ð23Þ
Note that it follows from x 6 ðbðlð1 pÞ þ xÞPÞ=ðlðl þ xÞÞ that
xI E
bI 6 0: 1 þ aI
Thus, from (22) and (23), we have
lðBÞ ¼
E0 xI bI E0 ¼ l x: l x þ max 0; E 1 þ aI E E
ð24Þ
Along each solution ðSðtÞ; EðtÞ; IðtÞÞ of system (20) with ðSð0Þ; Eð0Þ; Ið0ÞÞ 2 K, where K is the compact absorbing set. We have
1 t
Z 0
t
1 t
lðBÞds 6 log
EðtÞ ðl þ xÞ; Eð0Þ
2926
L.-M. Cai, X.-Z. Li / Applied Mathematical Modelling 33 (2009) 2919–2926
which implies that
2 ¼ lim sup sup q t!1
x2K
1 t
Z
t
lðBÞds 6
0
lþx 2
< 0:
, then the endemic equilibrium P of the model (20) is globally stable in X1 . This According to Lemma 4, if R0 > 1 and x 6 x completes the proof. h
5. Concluding remarks In this paper, we investigate the dynamical behavior of a SEIV epidemic model that incorporates a nonlinear incidence rate and a waning preventive vaccine. For the model (2), we obtain the basic reproduction number R0 . From the definition of R0 , it can be seen that if
pv ¼ ðl þ xÞ
rbP lðl þ cÞðl þ rÞ ; lrbP
then R0 ¼ 1. Since R0 is a decreasing function of pv , it follows that if p > pv , then R0 < 1; if p < pv , then R0 > 1. Thus, the condition for disease eradication from the community is satisfied if p > pv ; the disease will persist if p < pv . Note that if the vaccinated number p ¼ 0, then R0 ¼ ðrbPÞ=ðlðl þ cÞðl þ rÞÞ: Thus R0 can be rewritten as R0 ¼ ð1 lp=ðl þ xÞÞR0 . The expression for R0 can show vaccination coverage p with a vaccine that induces immunity, which wanes with average duration of protection ð1=xÞ in a population with average life expectancy ð1=lÞ. Therefore, we can obtain the optimal vaccine convergence pv ¼ lþlx ð1 R1 Þ. Clearly this expression also shows vaccination coverage (pv ) is an increasing function of x. 0
Thus, it is necessary and important for public health management to control an epidemic by increasing the duration of the loss of immunity induced by vaccination ð1=xÞ, which reduces the threshold vaccination convergence (pv ). Acknowledgements The authors would like to thank Prof. M. Cross and the anonymous referees for their useful feedback and helpful comments and suggestions which have improved the manuscript. The authors would also like to thank Prof. Michael Y. Li for his help. References [1] M.A. Teitelbaulm, M. Edmunds, Immunization and vaccine-preventable illness, Unites States, 1992–1997, Stat. Bull. Metrop. Insur. Co. 80 (2) (1999) 13–20. [2] S.A. Levin, T.G. Hallam, L.J. Gross, Applied Mathematical Ecology, Springer, New York, 1990. [3] Z. Ma, Y. Zhou, W. Wang, Z. Jin, Mathematical Models and Dynamics of Infectious Diseases, China Sciences Press, Beijing, 2004. [4] J.A. Yorke, W.P. London, Recurrent outbreaks of measles, chickenpox and mumps II, Am. J. Epidemiol. 98 (1973) 469–482. [5] W.M. Liu, H.W. Hethcote, S.A. Levin, Dynamical behavior of epidemiological models with nonlinear incidence rates, J. Math. Biol. 25 (1987) 359–380. [6] W.M. Liu, S.A. Levin, Y. Iwasa, Influence of nonlinear incidence rates upon the behavior of SIRS epidemiological models, J. Math. Biol. 23 (1986) 187– 204. [7] V. Capasso, G. Serio, A generalization of the Kermack–McKendrick deterministic epidemic model, Math. Biosci. 42 (1978) 43–61. [8] S. Ruan, W. Wang, Dynamical behavior of an epidemic model with a nonlinear incidence rate, J. Differen. Equat. 188 (2003) 135–163. [9] D. Xiao, S. Ruan, Global analysis of an epidemic model with nonmonotone incidence rate, Math. Biosci. 208 (2007) 419–429. [10] T. Zhang, Z. Teng, Pulse vaccination delayed SEIRS epidemic model with saturation incidence, Appl. Math. Model. 32 (2008) 1403–1416. [11] H.W. Hethcote, The mathematics of infectious diseases, SIAM Rev. 42 (2000) 599–653. [12] P. van den Driessche, J. Watmough, A simple SIS epidemic model with a backward bifurcation, J. Math. Biol. 40 (2000) 525–540. [13] M.E. Alexander, S.M. Moghadas, Periodicity in an epidemic model with a generalized nonlinear incidence, Math. Biosci. 189 (2004) 75–96. [14] M.E. Alexander, S.M. Moghadas, Bifurcation analysis of an SIRS Epidemic model with generalized incidence, SIAM J. Appl. Math. 65 (2005) 1794–1816. [15] R.H. Thieme, Persistence under relaxed point-dissipativity (with application to an endemic model), SIAM J. Math. Anal. 24 (1993) 407–435. [16] R.S. Varga, Matrix Iterative Analysis, Prentice-Hall, Inc., Englewood Cliffs, NJ, 1962. [17] A. Fonda, Uniformly persistent semidynamical systems, Proc. Amer. Math. Soc. 104 (1988) 111–116. [18] W. Ma, M. Song, Y. Takeuchi, Global stability of an SIR epidemic model with time delay, Appl. Math. Lett. 17 (2004) 1141–1145. [19] P. De Leenheer, H.L. Smith, Virus dynamics: a global analysis, SIAM J. Appl. Math. 63 (2003) 1313–1327. [20] S. Liu, E. Beretta, A stage-structured predator–prey model of Beddington–DeAngelis type, SIAM J. Appl. Math. 66 (2006) 1101–1129. [21] R. Xu, M.A.J. Chaplain, F.A. Davidson, Persistence and global stability of a ratio-dependent predator–prey model with stage structure, Appl. Math. Comput. 158 (2004) 729–744. [22] L. Cai, X. Song, Permanence and stability of a predator–prey system with stage-structure for predator, J. Comput. Appl. Math. 201 (2007) 356–366. [23] M.Y. Li, J.S. Muldowney, A geometric approach to the global-stability problems, SIAM J. Math. Anal. 27 (1996) 1070–1083. [24] L. Wang, M.Y. Li, Mathematical analysis of the global dynamics of a model for HIV infection of CD4þ T cells, Math. Biosci. 200 (2006) 44–57. [25] G.J. Butler, P. Waltman, Persistence in dynamical systems, Proc. Am. Math. Soc. 96 (1986) 425–430. [26] R.H. Martin Jr., Logarithmic norms and projections applied to linear differential systems, J. Math. Anal. Appl. 45 (1974) 432–454.