- Email: [email protected]
Global dynamics of an SEIR epidemic model with saturating contact rate
Mathematical Biosciences 185 (2003) 15–32 www.elsevier.com/locate/mbs
Global dynamics of an SEIR epidemic model with saturating contact rate Juan Zhang, Zhien Ma
*
Department of Mathematics, Xi’an Jiaotong University, Xi’an 710049, People’s Republic of China Received 11 June 2001; received in revised form 4 June 2002; accepted 10 June 2003
Abstract Heesterbeek and Metz [J. Math. Biol. 31 (1993) 529] derived an expression for the saturating contact rate of individual contacts in an epidemiological model. In this paper, the SEIR model with this saturating contact rate is studied. The basic reproduction number R0 is proved to be a sharp threshold which completely determines the global dynamics and the outcome of the disease. If R0 6 1, the disease-free equilibrium is globally stable and the disease always dies out. If R0 > 1, there exists a unique endemic equilibrium which is globally stable and the disease persists at an endemic equilibrium state if it initially exists. The contribution of the saturating contact rate to the basic reproduction number and the level of the endemic equilibrium is also analyzed. 2003 Elsevier Inc. All rights reserved. Keywords: SEIR model; Saturating contact rate; Asymptotically autonomous system; Competitive system; Orbital asymptotical stability
1. Introduction The incidence of a disease is the number of new cases per unit time and plays an important role in the study of mathematical epidemiology. Thieme and Castillo-Chavez [2] argued that the general form of a population size dependent incidence should be written as bCðN Þ NS I, where S and I are respectively the numbers of susceptibles and infectives at time t, b is the probability per unit time of transmitting the infection between two individuals taking part in a contact, and CðN Þ is the ÔunknownÕ probability for an individual to take part in a contact. Thus, CðN Þ is usually called the contact rate, and bCðN Þ, which is the average number of adequate contacts of an individual per *
Corresponding author. Tel.: +86-29 266 9227; fax: +86-29 323 7910. E-mail addresses: [email protected] (J. Zhang), [email protected] (Z. Ma).
0025-5564/$ - see front matter 2003 Elsevier Inc. All rights reserved. doi:10.1016/S0025-5564(03)00087-7
16
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
unit time, is said to be the adequate contact rate. An adequate contact is a contact which is sufficient for transmission of the infection from an infective to a susceptible. In most of the literature, the adequate contact rate frequently takes two forms. One is linearly proportional to the total population size N, or bN , so that the corresponding incidence is bilinear form bN NS I ¼ bSI, the other a constant k, the corresponding incidence k NS I is called standard form. When the total population size N is not too large, since the number of contacts made by an individual per unit time should increase as the total population size N increases, the linear adequate contact rate bN would be suitable. But when the total population size is quite large, since the number of contacts made by an infective per unit time should be limited, or should grow less rapidly as the total population size N increases, the linear adequate contact rate bN is not suitable and the constant adequate contact rate k may be more realistic. Hence the two adequate contact rates mentioned above are actually two extreme cases for the total population size N being very small and very large respectively. More generally, it may be reasonable to assume that the adequate contact rate is a function CðN Þ of the total population size N . Some reasonable demands on CðNÞ are that it should be a non-decreasing should be a non-increasing function of N . Furthermore, CðNÞ should function of N and that CðNÞ N behave linearly in N , for small N, and it should be independent of N , for large N . A problem that has been around for a long time in mathematical epidemiology is to give a mechanistic description of the saturation in the number of new contacts of infectives per unit time made by an individual. Frequently a Holling-type argument, borrowed from predator–prey systems, is thought to be a solution to the problem. However, as Heesterbeek and Metz explain in [1], on closer examination the application of the usual Holling argument to epidemic models cannot be justified. Of course, many other function forms can, and have been, suggested that have these properties, but a mechanistically derived form was lacking. Using a mechanistic argument, Heesterbeek and Metz [1] derived the expression for the saturating contact rate of individual contacts in a population that mixes randomly, that is CðNÞ ¼
bN pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; 1 þ bN þ 1 þ 2bN
ð1:1Þ
then we see that, for N small, CðNÞ bN , whereas for N large, CðNÞ 1. Furthermore, CðN Þ is Þ is non-increasing. non-decreasing and CðN N In this paper, we will consider a SEIR model with the saturating contact rate CðNÞ defined by (1.1). The demographic structure used in this SEIR model has recruitment and deaths. Let A be constant recruitment rate into the population. Let natural deaths occur at a rate proportional to the population size, so the death rate term is lN where l is the death rate constant. Thus, in the ¼ A lN . For this absence of disease, the differential equation for the population size N is dN dt demographic structure, the total population size N ðtÞ approaches A=l for any non-zero initial population size. Then the constant A=l is the carrying capacity. Although N varies in the finite interval ð0; A=lÞ, it can still be far greater (for example, A is very large and l is very small) than the value at which CðN Þ reaches its saturating state. Thus the saturation effect of the saturating contact rate CðNÞ can still take place, that is, it is reasonable that the saturating contact rate CðN Þ is used in the SEIR model with recruitment. Most of the research literature assumes that the disease incubation period is negligible so that once infected, each susceptible individual (in the class S) instantaneously becomes infectious (in the class I) and later recovers (in the class R) with a permanent or temporary acquired immunity.
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
17
A compartmental model based on these assumptions is customarily called an SIR or SIRS model respectively. The SIR and SIRS type models are widely studied [3–12]. In particular, Thieme [3] considered an SIRS epidemiological model with population size dependent contact rate and exponential demographics. Brauer and van den Driessche [4] have obtained the threshold and asymptotical behavior results of an SIR model with population size dependent contact rate and immigration of infectives. Many diseases, however, incubate inside the hosts for a period of time before the hosts become infectious. Models that are more general than the SIR or SIRS types need to be studied to investigate the role of an incubation period in disease transmission. Using a compartmental approach, one may assume that a susceptible individual first goes through a latent period (and is said to become exposed or in the class E) after infection, before becoming infectious. The resulting models are of SEIR or SEIRS types, respectively, depending on whether the acquired immunity is permanent or otherwise. Greenhalgh [13] considered SEIR models that incorporate density dependence in the death rate. Cooke and van den Driessche [14] introduced and studied SEIRS models with two delays. Recently, Greenhalgh [15] studied Hopf bifurcations in models of the SEIRS type with density dependent contact rate and death rate. Li and Muldowney [16] and Li et al. [17] studied the global dynamics of the SEIR models with a non-linear incidence rate and with a standard incidence, respectively. Li et al. [18] analyzed the global dynamics of a SEIR model with vertical transmission and a bilinear incidence. Research on epidemic models of SEIR or SEIRS types with the saturating contact rate CðN Þ defined by (1.1) are scarce in the literature. To the authorÕs knowledge, the present paper is the first one that gives a rigorous treatment of the global stability of the endemic equilibrium and establishes the complete global dynamics for the SEIR model with the saturating contact rate CðN Þ given by (1.1) and disease-related death rate. Epidemic models of this type are notorious for the fact that the method of Lyapunov functions has rarely worked for the proof of the global stability of the endemic equilibrium. In the present paper, to establish the global stability of the unique endemic equilibrium, we first introduce a change of variables by which our four-dimensional SEIR model can be reduced to a three-dimensional asymptotical autonomous differential system, and later use the theory of the asymptotical autonomous differential systems and the ideas of Li and Muldowney [16] who use the general theory of competitive systems and the relations, which is due to [19], between additive compound matrices and differential equations. This paper is organized as follows. In Section 2, the SEIR model with the saturating contact rate CðN Þ defined by (1.1) is formulated and the global stability of the disease-free equilibrium is obtained. In Section 3, the existence and uniqueness of the endemic equilibrium P is proved by means of geometric methods. At the same time the basic reproduction number is determined. Section 4 obtains the local stability of the endemic equilibrium P . Section 5 gives a rigorous treatment of the global stability of the endemic equilibrium. Thus, the global dynamics of the SEIR model considered here have been obtained completely.
2. Model formulation The total population size N ðtÞ is divided into four distinct epidemiological subclasses of individuals which are susceptible, exposed, infectious, and recovered (with permanent immunity), with sizes denoted by SðtÞ, EðtÞ, IðtÞ, and RðtÞ, respectively. Our assumptions on the dynamical transfer of the population are demonstrated in the following (Fig. 1).
18
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
Fig. 1. Transfer diagram for the SEIR model with saturating contact rate.
Here the parameters l, b, e0 , c0 are all positive constants, a0 is a non-negative constant and represents the death rate due to disease (the disease-related death rate). The constant l is the rate for natural death. b is the probability per unit time of transmitting the infection between two individuals taking part in a contact. c0 is the rate constant for recovery, so that c10 is the mean infective period. Also e0 is the rate constant at which the exposed individuals become infective, so that e10 is the mean latent period. The recovered individuals are assumed to acquire permanent immunity, so there is no transfer from the R class back to the S class. The positive constant A is the constant recruitment rate into the population, so that A=l represents a carrying capacity, or maximum pffiffiffiffiffiffiffiffiffiffi , bN derived in [1], where b is a positive constant, is possible population size. CðN Þ ¼ 1þbN þbN hðNÞ 1þ2bN the saturating contact rate of individual contacts in a population that mixes randomly. Then the SEIR model with the saturating contact rate CðNÞ is derived based on the basic assumptions and using the transfer diagram and described by the following system of differential equations 8 dS a0 SI ¼ A hðNÞ lS; > dt > > < dE a0 SI ¼ hðN Þ e0 E lE; dt ð2:1Þ dI > ¼ e0 E c0 I lI a0 I; > > dt : dR ¼ c0 I lR; dt pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where a0 ¼ bb, hðNÞ ¼ 1 þ bN þ 1 þ 2bN . Let ldt ¼ ds, we obtain the following system analogous to (2.1). 8 dS A aSI ¼ l hðN Þ S; > ds > > dE < aSI ¼ hðN ð1 þ eÞE; ds Þ ð2:2Þ dI > ¼ eE ð1 þ c þ aÞI; > ds > : dR ¼ cI R; ds where a ¼ al0 , e ¼ el0 , c ¼ cl0 , a ¼ al0 . The total population size N ðtÞ can be determined by ¼ Al N aI. We use N as a NðtÞ ¼ SðtÞ þ EðtÞ þ IðtÞ þ RðtÞ, or from the differential equation dN ds variable in place of the variable S to give the following system: 8 dE ¼ aðN EIRÞI ð1 þ eÞE; > > ds hðNÞ > < dI ¼ eE ð1 þ c þ aÞI; ds ð2:3Þ dR > ¼ cI R; > ds > : dN A ¼ l N aI: ds The system (2.3) is equivalent to the system (2.2). This allows us to attack (2.2) by studying the system (2.3). From biological considerations, we study (2.3) in the closed set
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
19
C ¼ fðE; I; R; NÞ 2 R4þ : 0 6 E þ I þ R 6 N 6 A=lg; where R4þ denotes the non-negative cone of R4 including its lower dimensional faces. It can be verified that C is positively invariant with respect to (2.3). oC and C denote the boundary and the interior of C, respectively. The point P0 ¼ ð0; 0; 0; A=lÞ 2 C is the disease-free equilibrium of (2.3), and it exists for all nonnegative values of the parameters. Any equilibrium in C corresponds to the disease being endemic and is called an endemic equilibrium. Equilibria are obtained by setting the right side of each of the four differential equations equal to zero, giving the following algebraic system 8 aðN E I RÞI ¼ ð1 þ eÞEhðN Þ; > > < eE ¼ ð1 þ c þ aÞI; ð2:4Þ cI ¼ R; > > : A N ¼ aI: l I and substitution of this into the first From the second of these equations, we have E ¼ 1þcþa e equation gives ð1 þ eÞð1 þ c þ aÞ hðNÞI: aðN E I RÞI ¼ e From this, we see that either I ¼ 0 (the disease-free equilibrium P0 ) or dx hðN Þ aðN E I RÞ ¼ e (the endemic equilibria), where d ¼ 1 þ c þ a, x ¼ 1 þ e. If I ¼ 0, it is easy to deduce that E ¼ R ¼ 0; N ¼ Al, the disease-free equilibrium P0 ¼ ð0; 0; 0; A=lÞ. Let R0 ¼
aeA e e0 ¼ bCðA=lÞ ¼ bCðA=lÞ ; ldxhðA=lÞ dxl ðl þ c0 þ a0 Þðl þ e0 Þ
ð2:5Þ
R0 is the number of secondary infections caused by a typical newly infected individual entering the population at the disease-free equilibrium during his or her entire infectious period, so R0 is called the basic reproduction number. Theorem 2.1. The disease-free equilibrium P0 ¼ ð0; 0; 0; A=lÞ of (2.3) is globally asymptotically stable in C if R0 6 1 and it is unstable if R0 > 1. Proof. Consider the function L ¼ eE þ xI. Its derivative along the solutions of (2.3) is dL aeðN E I RÞI ¼ dxI ds hðN Þ aN aeðE þ I þ RÞI ¼ eI dxI hðN Þ hðN Þ CðN Þ 6 dxI R0 1 : CðA=lÞ
ð2:6Þ
20
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
bN Since N 6 Al and CðN Þ ¼ hðNÞ is non-decreasing, thus dL 6 0 if R0 6 1. Furthermore, the Lyapunov– ds Lasalle Theorem (cf. [20], pp. 296–297) implies that all paths in C approach the largest positively ¼ 0. Here, M is the set where I ¼ 0. On the boundary of C invariant subset of the set M where dL ds ¼ R so R ¼ Rð0Þes ! 0 as s ! þ1, and dN ¼ Al N so where I ¼ 0, we have E ¼ 0, dR ds ds A A s A N ¼ l þ ½N ð0Þ le ! l as s ! þ1. Hence, all solution paths in C approach the disease-free equilibrium P0 . The Jacobian matrix of (2.3) at P0 is 1 0 aA 0 0 x lhðA=lÞ B e d 0 0 C C J0 ¼ B @ 0 c 1 0 A 0 a 0 1
and because there is only one non-zero element in the third column and in the fourth column, respectively, and they are both negative, we may reduce the question of whether the eigenvalues of this matrix have negative real part to the same question for the 2 · 2 matrix aA x lhðA=lÞ J01 ¼ : e d The trace of J01 is clearly negative, and the condition that the determinant of J01 is positive is R0 < 1. Thus the disease-free equilibrium is unstable if R0 > 1. Thus, these results establish the theorem. Theorem 2.1 completely determines the global dynamics of (2.3) in C when R0 6 1. Its epidemiological implication is that the sum of the exposed and the infectious subclasses in the population vanishes over time so the disease dies out. In Section 5, we show that the disease persists when R0 > 1.
3. Existence and uniqueness of endemic equilibria Global stability of P0 in C excludes the existence of equilibria other than P0 , thus the study of endemic equilibria is restricted to the case R0 > 1. C satisfy: From Section 2, the coordinates of an equilibrium P 2 8 aeðN E I RÞ ¼ dxhðNÞ; > > < eE ¼ dI; ð3:1Þ cI ¼ R; > > : A N ¼ aI; l
from the last three equations we have d A 1 A N ; I¼ N ; E¼ ae l a l
c R¼ a
A N : l
ð3:2Þ
Substituting (3.2) into the first equation in (3.1) gives a½dxN ðdx aeÞðA=lÞ ¼ adxhðNÞ:
ð3:3Þ
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
Noting that hðN Þ ¼ 1 þ bN þ
21
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ 2bN , (3.3) becomes
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi F ðN Þ,ða abÞdxN ½aðdx aeÞðA=lÞ þ adx adx 1 þ 2bN ¼ 0
ð3:4Þ
because R0 ¼
aeA aeA aeA a pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 ¼ ¼ ; ldxhðA=lÞ ð1 þ c þ aÞð1 þ eÞðl þ bA þ l 1 þ 2bðA=lÞÞ aebA ab
thus, a ab > 0 when R0 > 1. It is easy to see that F ð0Þ ¼ ½adx þ aðdx aeÞðA=lÞ adx < 0 and pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi F ðA=lÞ ¼ ða abÞdxðA=lÞ ½adx þ aðdx aeÞðA=lÞ adx 1 þ 2bðA=lÞ ¼ aaeðA=lÞ adxhðA=lÞ ¼ adxhðA=lÞðR0 1Þ: When R0 > 1, then F ðA=lÞ > 0. So Eq. (3.4) has at least one root N 2 ð0; A=lÞ. The derivative of F ðN Þ is dF ðN Þ ab ¼ dx a ab pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dN 1 þ 2bN að2abaÞ að2abaÞ . Together with which shows that F ðN Þ is decreasing if N < 2bðaabÞ 2 and increasing if N > 2bðaabÞ2
F ð0Þ < 0, we may claim that the root N 2 ð0; A=lÞ of Eq. (3.4) is unique. Therefore, we have the following theorem. Theorem 3.1. When R0 > 1, the system (2.3) has a unique endemic equilibrium P ¼ ðE ; I ; R ; N Þ. Here E ¼ aed ðAl N Þ, I ¼ 1a ðAl N Þ, and R ¼ ac ðAl N Þ.
4. Local stability of the endemic equilibrium P* Theorem 4.1. If R0 > 1, the endemic equilibrium P of the system (2.3) is locally asymptotically stable. Proof. The matrix of the linearization of the system (2.3) at the equilibrium P ¼ ðE ; I ; R ; N Þ is 0 1 aI dx aI aI x hðN hðN hðN q Þ Þ Þ e B e d 0 0 C C; J ¼ B @ 0 c 1 0 A 0 a 0 1
22
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
where
aI abI ðE þ I þ R Þ 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ q ¼ 1 þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h2 ðN Þ 1 þ 2bN hðN Þ 1 þ 2bN aI 2abI ðE þ I þ R Þ 6 þ h2 ðN Þ hðN Þ aI 2aI bN 6 þ hðN Þ hðN Þ hðN Þ 3aI : 6 hðN Þ The characteristic equation of J is
ðk þ 1Þðk3 þ a1 k2 þ a2 k þ a3 Þ ¼ 0; where
aI a1 ¼ ðd þ xÞ þ 1 þ > 0; hðN Þ aI a2 ¼ ðd þ xÞ 1 þ >0 hðN Þ
and a3 ¼ Then
aI ðd þ ec þ eÞ þ aeq > 0: hðN Þ
2 aI aI 3aI aI ae ðd þ ec þ eÞ a1 a2 a3 P ðd þ xÞ 1 þ þ ðd þ xÞ 1 þ hðN Þ hðN Þ hðN Þ hðN Þ 2 aI aI 2 ¼ ðd þ xÞ þ ðd þ xÞ 1 þ þ ½ðd þ xÞ2 3ae d e ec hðN Þ hðN Þ 2 aI aI 2 P ðd þ xÞ þ ðd þ xÞ 1 þ þ ð4dx 3ae d e ecÞ hðN Þ hðN Þ 2 aI aI 2 ð3dx 2aeÞ þ ¼ ðd þ xÞ þ ðd þ xÞ 1 þ hðN Þ hðN Þ 2
P 0: Hence the Routh–Hurwitz conditions are satisfied. Thus it follows that the endemic equilibrium P of (2.3), which exists if R0 > 1, is always locally asymptotically stable. The following results, which may be found in [21], are used in the next section for the study of the global stability of the endemic equilibrium P . Consider the following systems x0 ¼ f ðt; xÞ; ð4:1Þ y 0 ¼ gðyÞ;
ð4:2Þ
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
23
where f and g are continuous and locally Lipschitz in x 2 Rn and solutions exist for all positive time. System (4.1) is called asymptotically autonomous with limit system (4.2) if f ðt; xÞ ! gðxÞ as t ! 1 uniformly for x 2 Rn . Lemma 4.1. Let e be a locally asymptotically stable equilibrium of (4.2) and x be the x-limit set of a forward bounded solution xðtÞ of (4.1). If x contains a point y0 such that the solution of (4.2), with yð0Þ ¼ y0 converges to e as t ! 1, then x ¼ feg, i.e., xðtÞ ! e as t ! 1. Corollary 4.1. If solutions of the system (4.1) are bounded and the equilibrium e of the limit system (4.2) is globally asymptotically stable, then any solution xðtÞ of the system (4.1) satisfies xðtÞ ! e as t ! 1. 5. Global stability of the endemic equilibrium P* For the system (2.2) equal to the system (2.3), making the change of variable c aþc M ¼ ðS þ E þ IÞ þ R; a a then dM ¼ cA M so that the following system (5.2) is equivalent to (2.2). ds al 8 dS A aSI ¼ l hðNÞ S; > ds > > dE < aSI ¼ hðN xE; ds Þ dI > ¼ eE dI; > ds > : dM ¼ Ac M; al ds
ð5:1Þ
ð5:2Þ
a ðS þ E þ I þ MÞ, x ¼ 1 þ e, and d ¼ 1 þ a þ c. The last equation of the system where N ¼ aþc (5.2) implies that Ac Ac s þ Mð0Þ ð5:3Þ MðsÞ ¼ e : al al
Substituting (5.3) into the first three equations of the system (5.2) gives the following threedimensional non-autonomous differential system 8 dS A aSI ¼ l hðNÞ S; > ds > > > dE < ¼ aSI xE; ds hðN Þ ð5:4Þ dI ¼ eE > ds > dI; > > : N ¼ a S þ E þ I þ Ac þ Mð0Þ Ac es : aþc
al
al
a ðS þ E þ I þ Ac Þ,X as s ! þ1, the system (5.4) is a three-dimensional asympSince N ! aþc al totically autonomous differential system with limit system 8 dS A aSI > < ds ¼ l hðX Þ S; dE aSI ¼ hðX xE; ð5:5Þ ds Þ > : dI ¼ eE dI: ds
24
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
a Here X ¼ aþc ðS þ E þ I þ Ac Þ. It can be verified that the region al
T ¼ fðS; E; IÞ 2 Rþ 3 : 0 6 S þ E þ I 6 A=lg is positively invariant with respect to (5.5). Here R3þ denotes the non-negative cone of R3 including its lower dimensional faces. Thus, the system (5.5) is a bounded system. The main aim of this section is to prove the following Theorem 5.1. Once Theorem 5.1 is proved, together with Theorem 2.1, we can claim that the basic reproduction number R0 is a sharp threshold parameter and that the global dynamical behaviors of the system (2.3) and the outcome of the disease are completely determined. In other words, if R0 6 1, the disease-free equilibrium P0 is globally asymptotically stable so that the disease dies out, while if R0 > 1, the endemic equiC so that the disease remains endemic. librium P is globally asymptotically stable in Theorem 5.1. When R0 > 1, the unique endemic equilibrium P of the system (2.3) is globally asymptotically stable in the interior of C. Moreover, P attracts all trajectories in C except those starting with the invariant E ¼ I ¼ 0 which converge to P0 along this hyperplane. By inspecting the vector field given by (2.3), we see that all trajectories starting from the boundary oC of C enter the interior C of C except those on the N -axis which converge to P0 along this invariant axis. Thus we need only to prove that C is the attractive region of P . In accordance with Corollary 4.1, it is sufficient to prove that the equilibrium Q of (5.5), which is corresponding to the equilibrium P of (2.3), is globally asymptotically stable in the interior of the region T . For this purpose, we need to study the dynamics of the system (5.5). Recalled that R0 is defined by (2.5), first of all, it is easy to obtain the following result: Lemma 5.1. If R0 6 1, Q0 ¼ ðA=l; 0; 0Þ is the only equilibrium in T and it is globally asymptotically stable. If R0 > 1, Q0 becomes unstable while Q ¼ ðS ; E ; I Þ, which corresponds to the equilibrium P of the system (2.3), emerges as a unique equilibrium in the interior of T and it is locally asA N E I cA , E ¼ de I , I ¼ al 1a N , and N is the unique ymptotically stable, here S ¼ aþc a al root of (3.4). Lemma 5.2. (1) When R0 > 1, Q0 is the only omega limit point of the system (5.5) on the boundary of T and cannot be the omega limit point of any orbit which starts in the interior of T , of the system (5.5). (2) The system (5.5) is uniformly persistent when R0 > 1. Proof. It is easy to see that the vector field of (5.5) is transversal to the boundary of T on all its faces except the S-axis which is invariant with respect to (5.5). On the S-axis the equation for S is dS ¼ Al S, which implies that SðsÞ ! A=l as s ! þ1. Therefore, Q0 is the only x-limit point on ds the boundary of T . Let L ¼ eE þ xI. Its derivative along the solutions of (5.5) is dL aeSI S hðA=lÞ ¼ dxI ¼ dxI R0 1 : ds hðX Þ A=l hðX Þ > 0 for S sufficiently close to A=l except when I ¼ 0. Hence when R0 > 1, the If R0 > 1, then dL ds trajectories starting from the points which are close enough to Q0 must leave a neighborhood of Q0 except those on the invariant S-axis. This proves the first conclusion of the lemma.
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
25
System (5.5) is said to be uniformly persistent (see [3,22]) if there exists a constant 0 < c < 1 such that any solution ðSðsÞ; EðsÞ; IðsÞÞ with ðSð0Þ; Eð0Þ; Ið0ÞÞ 2 T satisfies lim inf jðSðtÞ; EðtÞ; IðtÞÞj P c: t!1
The second conclusion of the lemma follows from a uniform persistence result, Theorem 4.3, in [23]. Since the maximal invariant set on the boundary oT of T is the S-axis which is isolated, the hypothesis ðH Þ of [23] holds for (5.5). The condition for uniform persistence in Theorem 4.3 of [23] is equivalent to Q0 being unstable. Thus, when R0 > 1, the system (5.5) is uniformly persistent in T. The definitions and results on competitive systems, which will be used in the following lemma, can be found in Appendix B. Lemma 5.3. The system (5.5) is competitive in the convex region T . Also any non-empty compact omega limit set of (5.5) in the interior of T is either a closed orbit or the equilibrium Q . Proof. The Jacobian matrix of (5.5) at a point P ¼ ðS; E; IÞ 2 T is 0 1 0 ðX Þ 0 ðX Þ aSIh0 ðX Þ a aI a aS a 1 hðX þ aSIh hðX þ aSIh Þ Þ h2 ðX Þ aþc h2 ðX Þ aþc h2 ðX Þ aþc B C 0 0 0 aI aS J ¼@ aSIh2 ðX Þ a x aSIh2 ðX Þ a aSIh2 ðX Þ a A; hðX Þ
h ðX Þ
0
aþc
h ðX Þ
e
aþc
hðX Þ
h ðX Þ
aþc
ð5:6Þ
d
b ffi where h0 ðX Þ ¼ b þ pffiffiffiffiffiffiffiffiffi > 0. 1þ2bX pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a Note that X ¼ aþc ðS þ E þ I þ Ac Þ and hðX Þ ¼ 1 þ bX þ 1 þ 2bX , then al
aI aSIh0 ðX Þ a aI aIXh0 ðX Þ aI pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > 0; 2 P 2 ¼ hðX Þ h ðX Þ a þ c hðX Þ h ðX Þ hðX Þ 1 þ 2bX
ð5:7Þ
aS aSIh0 ðX Þ a aS aSXh0 ðX Þ aS pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > 0: 2 P 2 ¼ hðX Þ h ðX Þ a þ c hðX Þ h ðX Þ hðX Þ 1 þ 2bX
ð5:8Þ
Choosing the matrix H as H ¼ diagð1; 1; 1Þ, we can see that the system (5.5) is competitive in the convex region T , with respect to the partial ordering defined by the orthant fðS; E; IÞ 2 R3 : S 6 0; E P 0; I 6 0g. Suppose that X is an omega limit set of (5.5) in the interior of T . If X does not contain Q , since Q is the only interior equilibrium, then it contains no equilibria. Theorem B.1 implies that X is a closed orbit. Suppose that X contains Q , since Q is locally asymptotically stable whenever it exists in the interior of T , any orbit that gets arbitrarily close to Q must converge to Q . Thus X ¼ Q . Lemma 5.3 is proved. From Lemma 5.3, in the absence of closed orbits, all trajectories in T converge to Q when R0 > 1. Thus, the key to verifying the global stability of Q is to rule out the existence of periodic solutions of (5.5). This is achieved by showing that any periodic solution to (5.5) is orbitally asymptotically stable. The following lemma is the main tool to prove the global stability of the equilibrium Q . Its proof depends on the criterion for the asymptotical orbital stability of period
26
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
orbits given by Muldowney [19]. For the sake of convenience, some definitions and conclusions about the asymptotical orbital stability of period solutions and compound matrices which are relevant to this lemma give are shown in Appendix A. Lemma 5.4. The orbit of any non-constant periodic solution to the system (5.5), if it exists, is asymptotically orbitally stable with asymptotic phase. Proof. Suppose that pðsÞ ¼ ðSðsÞ; EðsÞ; IðsÞÞ is a periodic solution to (5.5) with period T > 0, whose orbit q is contained in the interior of T . Obviously, SðsÞ; EðsÞ, and IðsÞ are all bounded functions and have positive lower boundary for all s P 0 from Lemma 5.2. The Jacobian matrix of (5.5) concerned with this periodic orbit q is given by (5.6). The second compound system, whose definition can be found in Appendix A, with respect to the periodic solution pðsÞ ¼ ðSðsÞ; EðsÞ; IðsÞÞ of (5.5) is the following 8 aSIh0 ðX Þ a dU aI aS > ¼ 1 þ x þ U þ ðV þ W Þ; 2 > hðX Þ hðX Þ h ðX Þ aþc > < ds 0 0 ðX Þ a ðX Þ a dV aI ð5:9Þ ¼ eU 1 þ d þ hðX aSIh W; V þ aSIh ds Þ h2 ðX Þ aþc h2 ðX Þ aþc > > 0 0 > : dW ¼ aI aSIh ðX Þ a V d þ x þ aSIh ðX Þ a W : ds hðX Þ h2 ðX Þ aþc h2 ðX Þ aþc Using Theorem A.1 shown in Appendix A, Lemma 5.4 will follow if we can prove that the system (5.9) is asymptotically stable. For this purpose, we shall use the following Lyapunov function EðsÞ ðjV ðsÞj þ jW ðsÞjÞ : LðsÞ ¼ LðUðsÞ; V ðsÞ; W ðsÞ; SðsÞ; EðsÞ; IðsÞÞ ¼ sup jU ðsÞj; IðsÞ From Lemma 5.2, we obtain that the orbit q of pðsÞ remains at a positive distance from the boundary of T , therefore EðsÞ P k
and IðsÞ P k
with 0 < k < 1:
Hence the function L is well defined along pðsÞ and there exists a constant c > 0 such that LðU ; V ; W ; S; E; IÞ P c supfjU j; jV j þ jW jg:
ð5:10Þ
The right-hand derivative of LðsÞ exists. Together with (5.7) and (5.8), direct calculation yields aI aS aSIh0 ðX Þ a Dþ jU ðsÞj 6 1 þ x þ jU ðsÞj þ 2 ðjV ðsÞj þ jW ðsÞjÞ hðX Þ hðX Þ h ðX Þ a þ c aI aS aSIh0 ðX Þ a I E ðjV ðsÞj þ jW ðsÞjÞ 6 1þxþ jU ðsÞj þ 2 hðX Þ hðX Þ h ðX Þ a þ c E I and
aI aSIh0 ðX Þ a aSIh0 ðX Þ a Dþ jV ðsÞj 6 ejU ðsÞj 1 þ d þ 2 jV ðsÞj þ 2 jW ðsÞj; hðX Þ h ðX Þ a þ c h ðX Þ a þ c aI aSIh0 ðX Þ a aSIh0 ðX Þ a 2 jV ðsÞj d þ x þ 2 jW ðsÞj: Dþ jW ðsÞj 6 hðX Þ h ðX Þ a þ c h ðX Þ a þ c
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32 dI Let E0 ¼ dE and I 0 ¼ ds , then ds
27
E0 EI 0 E 2 ðjV ðsÞj þ jW ðsÞjÞ þ Dþ ðjV ðsÞj þ jW ðsÞjÞ I I I 0 E I0 E eE 6 ðjV ðsÞj þ jW ðsÞjÞ þ jU ðsÞj I E I I E E ð1 þ dÞ jV ðsÞj ðd þ xÞ jW ðsÞj I I 0 0 eE E I E jU ðsÞj þ ðjV ðsÞj þ jW ðsÞjÞ: 6 1d I I E I
E Dþ ðjV ðsÞj þ jW ðsÞjÞ ¼ I
Thus, we can obtain Dþ LðsÞ 6 supfg1 ðsÞ; g2 ðsÞgLðsÞ;
ð5:11Þ
where aI aSI aSI 2 h0 ðX Þ a þ hðX Þ EhðX Þ Eh2 ðX Þ a þ c aSI 6 1xþ EhðX Þ
g1 ðsÞ ¼ 1 x
and eE E0 I 0 þ 1 d: I E I We rewrite the last two equations of (5.5) as g2 ðsÞ ¼
aSI E0 ¼ þ x; EhðX Þ E eE I 0 ¼ þ d; I I we find E0 1: E Thus, from Eq. (5.11) and GronwallÕs inequality, we obtain supfg1 ðsÞ; g2 ðsÞg 6
LðsÞ 6
Lð0ÞEðsÞ s ALð0Þ s e 6 e ; Eð0Þ lEð0Þ
which implies that LðsÞ ! 0 as s ! 1. By (5.10) it turns out that LðU ðsÞ; V ðsÞ; W ðsÞÞ ! 0
as s ! 1:
This implies that the system (5.9) is asymptotically stable and therefore the periodic solution pðsÞ is asymptotically orbitally stable with asymptotical phase by Theorem A.1 in Appendix A. This proves Lemma 5.4.
28
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
As mentioned above, by Lemma 5.4, in the absence of closed orbits, all trajectories in T converge to Q when R0 > 1. Thus, we have the following theorem whose proof is similar to that of Theorem 2.1 in [16] and Theorem 6.5 in [17]. Theorem 5.2. If R0 > 1, the positive equilibrium Q of the system (5.5) is globally asymptotically stable in the interior of T . Thus, the main result of this section, Theorem 5.1, follows from Corollary 4.1. Therefore, the global dynamics of the SEIR model considered in this paper are completely determined. 6. Conclusions and discussion This paper has considered a SEIR model that incorporates recruitment and exponential natural death, as well as disease-related death, so that the total population size may vary in time. A distinguishing feature of the SEIR model considered here is that there is a saturating contact rate bN pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : CðNÞ ¼ 1 þ bN þ 1 þ 2bN The basic reproduction number R0 defined by (2.5) of this SEIR model (2.2) is a sharp threshold parameter which completely determines the global dynamics of the system (2.2) and the outcome of the disease. If R0 6 1, the disease-free equilibrium is globally stable so that the disease always dies out, and if R0 > 1, the disease-free equilibrium becomes unstable while the endemic equilibrium emerges as the unique positive equilibrium and it is globally stable in the interior of the feasible region C. In the special case when the mean latent period e10 ! 0 or the mean infective period c10 ! 1, the SEIR model here reduces to an SIR model or an SEI model with recruitment and saturating contact rate CðNÞ, respectively. Thus, the global dynamic results of these two models can be obtained as special cases when e10 ¼ 0 and c0 ¼ 0 respectively. Especially, to prove the global stability of the endemic equilibrium, we make a change of variable by which our fourdimensional model can be reduced to a three-dimensional asymptotical autonomous system with limit system. The similar change of variable M ¼ ac N þ R can reduce the SIR model with a general contact rate bðN Þ (the incidence is bðN ÞSI) in [4] by Brauer and van den Driessche to the following two-dimensional asymptotical autonomous differential system dI aþc cA cA dt ða þ c þ dÞI ð6:1Þ ¼ pA þ bðNÞI N I Mð0Þ e dt a ad ad and
dN dt
¼ A dN aI with the following limit system. dI aþc cA ¼ pA þ bðNÞI N I ða þ c þ dÞI dt a ad
ð6:2Þ
¼ A dN aI. The Dulac function BðI; N Þ ¼ 1I can exclude the existence of period soluand dN dt tions to (6.2). The global stability of the endemic equilibrium follows from Corollary 4.1 and the Poincare–Bendixson theorem. This change of variable can also reduce the proof of the global stability of the endemic equilibrium of the system (5) in [4]. Additionally, because of the complex pffiffiffiffiffiffiffiffiffiffi used in our SEIR model, many mathematical techniques structure of the incidence rate 1þbN þbbSI 1þ2bN
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
29
are used in the present paper to overcome mathematical difficulties that are not present in the SEIR model with bilinear incidence or standard incidence. The saturating contact rate CðN Þ defined by (1.1) is increasing, approximately a linear function bN for small N , and asymptotically approaches a constant value 1 as N ! 1. Using notation from Section 2, the SEIR model with recruitment and the standard incidence b NS I or the bilinear incidence bbN NS I ¼ bbSI can be respectively described by the following systems dS A bSI ¼ S; ds l N dE bSI ¼ xE; ds N dI ¼ eE dI; ds dR ¼ cI R ds
ð6:3Þ
and dS A ¼ bbSI S; ds l dE ¼ bbSI xE; ds ð6:4Þ dI ¼ eE dI; ds dR ¼ cI R: ds Thus both (6.3) and (6.4) are two extreme cases of the SEIR model in this paper. Noting that R0 ¼
aAe be0 bA=l pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; ¼ ldxhðA=lÞ ðl þ e0 Þðl þ a0 þ c0 Þ ð1 þ bA=l þ 1 þ 2bA=lÞ
be0 ,R01 P R0 and if N is very if letting N ! 1 (then, Al ! 1), since CðN Þ ! 1, R0 ! ðlþe0 Þðlþa 0 þc0 Þ small, since CðN Þ is approximately a linear function bN , R0 is approximately be0 bA ,R02 P R0 . It is easy to verify that R01 and R02 are the basic reproduction numbers ðlþe0 Þðlþa0 þc0 Þ l I ; d I ; I ; cI1 Þ and ðAl dx I ; d I ; I ; cI2 Þ are the of the system (6.3) and (6.4), respectively. ðAl dx e 1 e 1 1 e 2 e 2 2 endemic equilibrium of (6.3) when R01 > 1 and the endemic equilibrium (6.4) when R02 > 1, rea ae A , I2 ¼ 1a ðA=l N2 Þ, and N2 ¼ bb þ ð1 dx Þ l. spectively, where I1 ¼ 1a ðA=l N1 Þ, N1 ¼ bðdxaeÞA dxðbaÞl Substituting N1 and N2 into (3.4) gives pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi F ðN1 Þ ¼ adxð1 þ 1 þ 2bN1 Þ < 0
and
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi F ðN2 Þ ¼ abðdx aeÞA=l a2 dx=b adx 1 þ 2bN2 < 0:
Thus, N1 < N and N2 < N , then I1 > I and I2 > I . Thus, although the threshold behavior and dynamic behavior of the SEIR model (2.2) in this paper are similar to those of the SEIR model (6.3) and (6.4), the basic reproduction number R0 of (2.2) is less than R01 and R02 , the endemic
30
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
equilibrium population size N is grater than those of (6.3) and (6.4), and the infective level of the endemic equilibrium state is less than those of (6.3) and (6.4). The epidemiological and demographical phenomena as well as their interaction can be observed in our model. For example, when the disease persists, the disease-related deaths can reduce the equilibrium population size below its maximum carrying capacity. The model considered here cannot give rise to periodic solutions. Acknowledgements This research is supported by the National Natural Sciences Foundation of PeopleÕs Republic of China grant no. 19971066. We also wish to thank two anonymous referees for helpful criticism and suggestions. Appendix A Let D Rn be an open set, and gðyÞ 2 Rn be a C 1 function defined in D. We consider the autonomous system in Rn y 0 ¼ gðyÞ: ðA:1Þ Let yðt; y0 Þ denote the solution of (A.1) such that yð0; y0 Þ ¼ y0 . The linear variational equation of (A.1) with respect to yðt; y0 Þ is given by og ðA:2Þ Y 0 ðtÞ ¼ ðyðt; y0 ÞÞY ðtÞ; oy where og is the Jacobian matrix of g. Suppose the autonomous system (A.1) has a periodic solution oy y ¼ pðtÞ with least period T > 0 and orbit q ¼ fpðtÞ : 0 6 t 6 T g. This orbit is orbitally stable if for each > 0, there exists a r > 0 such that any solution yðtÞ, for which the distance of yð0Þ from q is less than r remains at a distance less than from q for all t P 0. It is asymptotically orbitally stable if the distance of yðtÞ from q also tends to zero as t ! 1. This orbit q is asymptotically orbitally stable with asymptotic phase if it is asymptotically orbitally stable and there is a h > 0 such that any solution yðtÞ for which the distance of yð0Þ from q is less than h, satisfies jyðtÞ pðt t Þj ! 0 as t ! 1 for some t which may depend on yð0Þ. The following is a criterion given in [19] for the asymptotical orbital stability of a periodic orbit to the general autonomous system (A.1). Theorem A.1. A sufficient condition for a period orbit q ¼ fpðtÞ : 0 6 t 6 T g of (A.1) to be asymptotically orbitally stable with asymptotic phase is that the linear system ! ½2 og Z 0 ðtÞ ¼ ðpðtÞÞ ZðtÞ ðA:3Þ oy be asymptotically stable. ½2
Eq. (A.3) is called the second compound equation of (A.2) and the matrix og ðpðtÞÞ is the second oy at y ¼ pðtÞ of g. Generally speaking, for a n n matrix compound matrix of the Jacobian matrix og oy
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
31
A and integer1 6k 6 n, the kth additive compound matrix of A is denoted by A½k . This is a N N n matrix, N ¼ , defined by k A½k ¼ Dþ ðI þ hAÞðkÞ jh¼0 ; where BðkÞ is the kth exterior power of a n n matrix B and Dþ denotes the right-hand derivative. A survey on the definition and properties of additive compound matrices together with their connections to differential equations may be found in [19,24]. Further applications may be found in [25,26]. The following are the examples of various compound matrices when n ¼ 3. 0 1 a11 a12 a13 A½1 ¼ A ¼ @ a21 a22 a23 A; a31 a32 a33 0 1 a11 þ a22 a23 a13 a11 þ a33 a12 A A½2 ¼ @ a32 a31 a21 a22 þ a33 and A½3 ¼ a11 þ a22 þ a33 ¼ tr A:
Appendix B All concepts and results in this appendix can be found in [16,17,27,28]. The autonomous system (A.1) is said to be competitive in D if, for some diagonal matrix H has non-positive off-diagonal elements H ¼ diagðr1 ; . . . ; rn Þ where each ri is either 1 or )1, H og oy for all y 2 D. If D is convex, the flow of such a competitive system preserves, for t < 0, the partial ordering in Rn defined by the orthant K ¼ fðy1 ; . . . ; yn Þ 2 Rn : ri yi P 0; i ¼ 1; 2; . . . ; ng. The concept of competitiveness used above is more general than the one in [28] in that the partial ordering is not necessarily defined by the standard orthant of Rn . However, by a linear change of variables z ¼ Hy, a competitive system as defined above can be transformed into a system that is ÔcompetitiveÕ in the sense of [28]. An important characteristic of a three-dimensional competitive system is the following Poincare–Bendixson property. Theorem B.1. Suppose that n ¼ 3, D is convex, (A.1) is competitive in D and X is a non-empty compact omega limit set of (A.1). If X contains no equilibria, then X must be a closed orbit (cf. [28] Theorem 1). References [1] J.A.P. Heesterbeek, J.A.J. Metz, The saturating contact rate in marriage and epidemic models, J. Math. Biol. 31 (1993) 529. [2] H.R. Thieme, C. Castillo-Chavez, On the role of variable infectivity in the dynamics of the human immunodeficiency virus epidemic, in: C. Castillo-Chavez (Ed.), Mathematical and Statistical Approaches to AIDS Epidemiology, Lecture Notes in Biomathematics, vol. 83, Springer, Berlin, 1989, p. 157.
32
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
[3] H.R. Thieme, Epidemic and demographic interaction in the spread of potentially fatal diseases in growing populations, Math. Biosci. 111 (1992) 99. [4] F. Brauer, P. van den Driessche, Models for transmission of disease with immigration of infectives, Math. Biosci. 171 (2001) 143. [5] J. Mena-Lorca, H.W. Hethcote, Dynamic models of infectious diseases as regulators of population sizes, J. Math. Biol. 30 (1992) 693. [6] F. Brauer, Epidemic models in populations of varying size, in: C. Castillo-Chavez, S.A. Levin, C. Shoemaker (Eds.), Mathematical Approaches to Problems in Resource Management and Epidemiology, Lecture Notes in Biomathematics, vol. 81, Springer, p. 109, 1989. [7] L.Q. Gao, H.W. Hethcote, Disease transmission models with density-dependent demographics, J. Math. Biol. 30 (1992) 717. [8] D. Greenhalgh, Some threshold and stability results for epidemic models with a density dependent death rate, Theor. Pop. Biol. 42 (1992) 130. [9] H.J. Bremermann, H.R. Thieme, A competitive exclusion principle for pathogen virulence, J. Math. Biol. 27 (1989) 179. [10] H.W. Hethcote, P. van den Driessche, Some epidemiological models with nonlinear incidence, J. Math. Biol. 29 (1991) 271. [11] W.R. Derrick, P. van den Driessche, A disease transmission model in a non-constant population, J. Math. Biol. 31 (1993) 495. [12] S. Busenberg, P. van den Driessche, Analysis of a disease transmission model with varying population size, J. Math. Biol. 29 (1990) 257. [13] D. Greenhalgh, Some results for a SEIR epidemic model with density dependence in the death rate, IMA J. Math. Appl. Med. Biol. 9 (1992) 67. [14] K. Cooke, P. van den Driessche, Analysis of an SEIRS epidemic model with two delays, J. Math. Biol. 35 (1996) 240. [15] D. Greenhalgh, Hopf bifurcation in epidemic models with a latent period and nonpermanent immunity, Math. Comput. Model. 25 (1997) 85. [16] M.Y. Li, J.S. Muldowney, Global stability for the SEIR model in epidemiology, Math. Biosci. 125 (1995) 155. [17] M.Y. Li, J.R. Graef, L. Wang, J. Karsai, Global dynamics of a SEIR model with varying total population size, Math. Biosci. 160 (1999) 191. [18] M.Y. Li, H.L. Smith, L. Wang, Global dynamics of an SEIR epidemic model with vertical transmission, SIAM J. Appl. Math. 62 (2001) 58. [19] J.S. Muldowney, Compound matrices and ordinary differential equations, Rocky Mount. J. Math. 20 (1990) 857. [20] J.K. Hale, Ordinary Differential Equations, John Wiley, New York, 1969. [21] H.R. Thieme, Convergence results and a Poincare–Bendixson trichotomy for asymptotically automous differential equations, J. Math. Biol. 30 (1992) 755. [22] G. Butler, P. Waltman, Persistence in dynamical systems, Proc. Am. Math. Soc. 96 (1986) 425. [23] H.I. Freedman, S.G. Ruan, M.X. Tang, Uniform persistence and flows near a closed positively invariant set, J. Dynam. Diff. Equat. 6 (1994) 583. [24] D. London, On derivations arising in differential equations, Linear Multilin. Algeb. 4 (1976) 179. [25] Y. Li, J.S. Muldowney, Evolution of surface functionals and differential equations, in: J. Wiener, J. Hale (Eds.), Ordinary and Delay Differential Equations, Pitman Research Notes in Mathematics, vol. 272, Longman, England, 1992, p. 144. [26] Y. Li, J.S. Muldowney, On BendixsonÕs criterion, J. Diff. Equat. 106 (1993) 27. [27] H.L. Smith, Systems of ordinary differential equations which generate an order preserving flow. A survey of results, SIAM Rev. 30 (1988) 87. [28] M.W. Hirsch, Systems of differential equations which are competitive or cooperative. IV: Structural stabilities in three dimensional systems, SIAM J. Math. Anal. 21 (1990) 1225.
Global dynamics of an SEIR epidemic model with saturating contact rate Juan Zhang, Zhien Ma
*
Department of Mathematics, Xi’an Jiaotong University, Xi’an 710049, People’s Republic of China Received 11 June 2001; received in revised form 4 June 2002; accepted 10 June 2003
Abstract Heesterbeek and Metz [J. Math. Biol. 31 (1993) 529] derived an expression for the saturating contact rate of individual contacts in an epidemiological model. In this paper, the SEIR model with this saturating contact rate is studied. The basic reproduction number R0 is proved to be a sharp threshold which completely determines the global dynamics and the outcome of the disease. If R0 6 1, the disease-free equilibrium is globally stable and the disease always dies out. If R0 > 1, there exists a unique endemic equilibrium which is globally stable and the disease persists at an endemic equilibrium state if it initially exists. The contribution of the saturating contact rate to the basic reproduction number and the level of the endemic equilibrium is also analyzed. 2003 Elsevier Inc. All rights reserved. Keywords: SEIR model; Saturating contact rate; Asymptotically autonomous system; Competitive system; Orbital asymptotical stability
1. Introduction The incidence of a disease is the number of new cases per unit time and plays an important role in the study of mathematical epidemiology. Thieme and Castillo-Chavez [2] argued that the general form of a population size dependent incidence should be written as bCðN Þ NS I, where S and I are respectively the numbers of susceptibles and infectives at time t, b is the probability per unit time of transmitting the infection between two individuals taking part in a contact, and CðN Þ is the ÔunknownÕ probability for an individual to take part in a contact. Thus, CðN Þ is usually called the contact rate, and bCðN Þ, which is the average number of adequate contacts of an individual per *
Corresponding author. Tel.: +86-29 266 9227; fax: +86-29 323 7910. E-mail addresses: [email protected] (J. Zhang), [email protected] (Z. Ma).
0025-5564/$ - see front matter 2003 Elsevier Inc. All rights reserved. doi:10.1016/S0025-5564(03)00087-7
16
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
unit time, is said to be the adequate contact rate. An adequate contact is a contact which is sufficient for transmission of the infection from an infective to a susceptible. In most of the literature, the adequate contact rate frequently takes two forms. One is linearly proportional to the total population size N, or bN , so that the corresponding incidence is bilinear form bN NS I ¼ bSI, the other a constant k, the corresponding incidence k NS I is called standard form. When the total population size N is not too large, since the number of contacts made by an individual per unit time should increase as the total population size N increases, the linear adequate contact rate bN would be suitable. But when the total population size is quite large, since the number of contacts made by an infective per unit time should be limited, or should grow less rapidly as the total population size N increases, the linear adequate contact rate bN is not suitable and the constant adequate contact rate k may be more realistic. Hence the two adequate contact rates mentioned above are actually two extreme cases for the total population size N being very small and very large respectively. More generally, it may be reasonable to assume that the adequate contact rate is a function CðN Þ of the total population size N . Some reasonable demands on CðNÞ are that it should be a non-decreasing should be a non-increasing function of N . Furthermore, CðNÞ should function of N and that CðNÞ N behave linearly in N , for small N, and it should be independent of N , for large N . A problem that has been around for a long time in mathematical epidemiology is to give a mechanistic description of the saturation in the number of new contacts of infectives per unit time made by an individual. Frequently a Holling-type argument, borrowed from predator–prey systems, is thought to be a solution to the problem. However, as Heesterbeek and Metz explain in [1], on closer examination the application of the usual Holling argument to epidemic models cannot be justified. Of course, many other function forms can, and have been, suggested that have these properties, but a mechanistically derived form was lacking. Using a mechanistic argument, Heesterbeek and Metz [1] derived the expression for the saturating contact rate of individual contacts in a population that mixes randomly, that is CðNÞ ¼
bN pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; 1 þ bN þ 1 þ 2bN
ð1:1Þ
then we see that, for N small, CðNÞ bN , whereas for N large, CðNÞ 1. Furthermore, CðN Þ is Þ is non-increasing. non-decreasing and CðN N In this paper, we will consider a SEIR model with the saturating contact rate CðNÞ defined by (1.1). The demographic structure used in this SEIR model has recruitment and deaths. Let A be constant recruitment rate into the population. Let natural deaths occur at a rate proportional to the population size, so the death rate term is lN where l is the death rate constant. Thus, in the ¼ A lN . For this absence of disease, the differential equation for the population size N is dN dt demographic structure, the total population size N ðtÞ approaches A=l for any non-zero initial population size. Then the constant A=l is the carrying capacity. Although N varies in the finite interval ð0; A=lÞ, it can still be far greater (for example, A is very large and l is very small) than the value at which CðN Þ reaches its saturating state. Thus the saturation effect of the saturating contact rate CðNÞ can still take place, that is, it is reasonable that the saturating contact rate CðN Þ is used in the SEIR model with recruitment. Most of the research literature assumes that the disease incubation period is negligible so that once infected, each susceptible individual (in the class S) instantaneously becomes infectious (in the class I) and later recovers (in the class R) with a permanent or temporary acquired immunity.
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
17
A compartmental model based on these assumptions is customarily called an SIR or SIRS model respectively. The SIR and SIRS type models are widely studied [3–12]. In particular, Thieme [3] considered an SIRS epidemiological model with population size dependent contact rate and exponential demographics. Brauer and van den Driessche [4] have obtained the threshold and asymptotical behavior results of an SIR model with population size dependent contact rate and immigration of infectives. Many diseases, however, incubate inside the hosts for a period of time before the hosts become infectious. Models that are more general than the SIR or SIRS types need to be studied to investigate the role of an incubation period in disease transmission. Using a compartmental approach, one may assume that a susceptible individual first goes through a latent period (and is said to become exposed or in the class E) after infection, before becoming infectious. The resulting models are of SEIR or SEIRS types, respectively, depending on whether the acquired immunity is permanent or otherwise. Greenhalgh [13] considered SEIR models that incorporate density dependence in the death rate. Cooke and van den Driessche [14] introduced and studied SEIRS models with two delays. Recently, Greenhalgh [15] studied Hopf bifurcations in models of the SEIRS type with density dependent contact rate and death rate. Li and Muldowney [16] and Li et al. [17] studied the global dynamics of the SEIR models with a non-linear incidence rate and with a standard incidence, respectively. Li et al. [18] analyzed the global dynamics of a SEIR model with vertical transmission and a bilinear incidence. Research on epidemic models of SEIR or SEIRS types with the saturating contact rate CðN Þ defined by (1.1) are scarce in the literature. To the authorÕs knowledge, the present paper is the first one that gives a rigorous treatment of the global stability of the endemic equilibrium and establishes the complete global dynamics for the SEIR model with the saturating contact rate CðN Þ given by (1.1) and disease-related death rate. Epidemic models of this type are notorious for the fact that the method of Lyapunov functions has rarely worked for the proof of the global stability of the endemic equilibrium. In the present paper, to establish the global stability of the unique endemic equilibrium, we first introduce a change of variables by which our four-dimensional SEIR model can be reduced to a three-dimensional asymptotical autonomous differential system, and later use the theory of the asymptotical autonomous differential systems and the ideas of Li and Muldowney [16] who use the general theory of competitive systems and the relations, which is due to [19], between additive compound matrices and differential equations. This paper is organized as follows. In Section 2, the SEIR model with the saturating contact rate CðN Þ defined by (1.1) is formulated and the global stability of the disease-free equilibrium is obtained. In Section 3, the existence and uniqueness of the endemic equilibrium P is proved by means of geometric methods. At the same time the basic reproduction number is determined. Section 4 obtains the local stability of the endemic equilibrium P . Section 5 gives a rigorous treatment of the global stability of the endemic equilibrium. Thus, the global dynamics of the SEIR model considered here have been obtained completely.
2. Model formulation The total population size N ðtÞ is divided into four distinct epidemiological subclasses of individuals which are susceptible, exposed, infectious, and recovered (with permanent immunity), with sizes denoted by SðtÞ, EðtÞ, IðtÞ, and RðtÞ, respectively. Our assumptions on the dynamical transfer of the population are demonstrated in the following (Fig. 1).
18
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
Fig. 1. Transfer diagram for the SEIR model with saturating contact rate.
Here the parameters l, b, e0 , c0 are all positive constants, a0 is a non-negative constant and represents the death rate due to disease (the disease-related death rate). The constant l is the rate for natural death. b is the probability per unit time of transmitting the infection between two individuals taking part in a contact. c0 is the rate constant for recovery, so that c10 is the mean infective period. Also e0 is the rate constant at which the exposed individuals become infective, so that e10 is the mean latent period. The recovered individuals are assumed to acquire permanent immunity, so there is no transfer from the R class back to the S class. The positive constant A is the constant recruitment rate into the population, so that A=l represents a carrying capacity, or maximum pffiffiffiffiffiffiffiffiffiffi , bN derived in [1], where b is a positive constant, is possible population size. CðN Þ ¼ 1þbN þbN hðNÞ 1þ2bN the saturating contact rate of individual contacts in a population that mixes randomly. Then the SEIR model with the saturating contact rate CðNÞ is derived based on the basic assumptions and using the transfer diagram and described by the following system of differential equations 8 dS a0 SI ¼ A hðNÞ lS; > dt > > < dE a0 SI ¼ hðN Þ e0 E lE; dt ð2:1Þ dI > ¼ e0 E c0 I lI a0 I; > > dt : dR ¼ c0 I lR; dt pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi where a0 ¼ bb, hðNÞ ¼ 1 þ bN þ 1 þ 2bN . Let ldt ¼ ds, we obtain the following system analogous to (2.1). 8 dS A aSI ¼ l hðN Þ S; > ds > > dE < aSI ¼ hðN ð1 þ eÞE; ds Þ ð2:2Þ dI > ¼ eE ð1 þ c þ aÞI; > ds > : dR ¼ cI R; ds where a ¼ al0 , e ¼ el0 , c ¼ cl0 , a ¼ al0 . The total population size N ðtÞ can be determined by ¼ Al N aI. We use N as a NðtÞ ¼ SðtÞ þ EðtÞ þ IðtÞ þ RðtÞ, or from the differential equation dN ds variable in place of the variable S to give the following system: 8 dE ¼ aðN EIRÞI ð1 þ eÞE; > > ds hðNÞ > < dI ¼ eE ð1 þ c þ aÞI; ds ð2:3Þ dR > ¼ cI R; > ds > : dN A ¼ l N aI: ds The system (2.3) is equivalent to the system (2.2). This allows us to attack (2.2) by studying the system (2.3). From biological considerations, we study (2.3) in the closed set
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
19
C ¼ fðE; I; R; NÞ 2 R4þ : 0 6 E þ I þ R 6 N 6 A=lg; where R4þ denotes the non-negative cone of R4 including its lower dimensional faces. It can be verified that C is positively invariant with respect to (2.3). oC and C denote the boundary and the interior of C, respectively. The point P0 ¼ ð0; 0; 0; A=lÞ 2 C is the disease-free equilibrium of (2.3), and it exists for all nonnegative values of the parameters. Any equilibrium in C corresponds to the disease being endemic and is called an endemic equilibrium. Equilibria are obtained by setting the right side of each of the four differential equations equal to zero, giving the following algebraic system 8 aðN E I RÞI ¼ ð1 þ eÞEhðN Þ; > > < eE ¼ ð1 þ c þ aÞI; ð2:4Þ cI ¼ R; > > : A N ¼ aI: l I and substitution of this into the first From the second of these equations, we have E ¼ 1þcþa e equation gives ð1 þ eÞð1 þ c þ aÞ hðNÞI: aðN E I RÞI ¼ e From this, we see that either I ¼ 0 (the disease-free equilibrium P0 ) or dx hðN Þ aðN E I RÞ ¼ e (the endemic equilibria), where d ¼ 1 þ c þ a, x ¼ 1 þ e. If I ¼ 0, it is easy to deduce that E ¼ R ¼ 0; N ¼ Al, the disease-free equilibrium P0 ¼ ð0; 0; 0; A=lÞ. Let R0 ¼
aeA e e0 ¼ bCðA=lÞ ¼ bCðA=lÞ ; ldxhðA=lÞ dxl ðl þ c0 þ a0 Þðl þ e0 Þ
ð2:5Þ
R0 is the number of secondary infections caused by a typical newly infected individual entering the population at the disease-free equilibrium during his or her entire infectious period, so R0 is called the basic reproduction number. Theorem 2.1. The disease-free equilibrium P0 ¼ ð0; 0; 0; A=lÞ of (2.3) is globally asymptotically stable in C if R0 6 1 and it is unstable if R0 > 1. Proof. Consider the function L ¼ eE þ xI. Its derivative along the solutions of (2.3) is dL aeðN E I RÞI ¼ dxI ds hðN Þ aN aeðE þ I þ RÞI ¼ eI dxI hðN Þ hðN Þ CðN Þ 6 dxI R0 1 : CðA=lÞ
ð2:6Þ
20
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
bN Since N 6 Al and CðN Þ ¼ hðNÞ is non-decreasing, thus dL 6 0 if R0 6 1. Furthermore, the Lyapunov– ds Lasalle Theorem (cf. [20], pp. 296–297) implies that all paths in C approach the largest positively ¼ 0. Here, M is the set where I ¼ 0. On the boundary of C invariant subset of the set M where dL ds ¼ R so R ¼ Rð0Þes ! 0 as s ! þ1, and dN ¼ Al N so where I ¼ 0, we have E ¼ 0, dR ds ds A A s A N ¼ l þ ½N ð0Þ le ! l as s ! þ1. Hence, all solution paths in C approach the disease-free equilibrium P0 . The Jacobian matrix of (2.3) at P0 is 1 0 aA 0 0 x lhðA=lÞ B e d 0 0 C C J0 ¼ B @ 0 c 1 0 A 0 a 0 1
and because there is only one non-zero element in the third column and in the fourth column, respectively, and they are both negative, we may reduce the question of whether the eigenvalues of this matrix have negative real part to the same question for the 2 · 2 matrix aA x lhðA=lÞ J01 ¼ : e d The trace of J01 is clearly negative, and the condition that the determinant of J01 is positive is R0 < 1. Thus the disease-free equilibrium is unstable if R0 > 1. Thus, these results establish the theorem. Theorem 2.1 completely determines the global dynamics of (2.3) in C when R0 6 1. Its epidemiological implication is that the sum of the exposed and the infectious subclasses in the population vanishes over time so the disease dies out. In Section 5, we show that the disease persists when R0 > 1.
3. Existence and uniqueness of endemic equilibria Global stability of P0 in C excludes the existence of equilibria other than P0 , thus the study of endemic equilibria is restricted to the case R0 > 1. C satisfy: From Section 2, the coordinates of an equilibrium P 2 8 aeðN E I RÞ ¼ dxhðNÞ; > > < eE ¼ dI; ð3:1Þ cI ¼ R; > > : A N ¼ aI; l
from the last three equations we have d A 1 A N ; I¼ N ; E¼ ae l a l
c R¼ a
A N : l
ð3:2Þ
Substituting (3.2) into the first equation in (3.1) gives a½dxN ðdx aeÞðA=lÞ ¼ adxhðNÞ:
ð3:3Þ
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
Noting that hðN Þ ¼ 1 þ bN þ
21
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ 2bN , (3.3) becomes
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi F ðN Þ,ða abÞdxN ½aðdx aeÞðA=lÞ þ adx adx 1 þ 2bN ¼ 0
ð3:4Þ
because R0 ¼
aeA aeA aeA a pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6 ¼ ¼ ; ldxhðA=lÞ ð1 þ c þ aÞð1 þ eÞðl þ bA þ l 1 þ 2bðA=lÞÞ aebA ab
thus, a ab > 0 when R0 > 1. It is easy to see that F ð0Þ ¼ ½adx þ aðdx aeÞðA=lÞ adx < 0 and pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi F ðA=lÞ ¼ ða abÞdxðA=lÞ ½adx þ aðdx aeÞðA=lÞ adx 1 þ 2bðA=lÞ ¼ aaeðA=lÞ adxhðA=lÞ ¼ adxhðA=lÞðR0 1Þ: When R0 > 1, then F ðA=lÞ > 0. So Eq. (3.4) has at least one root N 2 ð0; A=lÞ. The derivative of F ðN Þ is dF ðN Þ ab ¼ dx a ab pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi dN 1 þ 2bN að2abaÞ að2abaÞ . Together with which shows that F ðN Þ is decreasing if N < 2bðaabÞ 2 and increasing if N > 2bðaabÞ2
F ð0Þ < 0, we may claim that the root N 2 ð0; A=lÞ of Eq. (3.4) is unique. Therefore, we have the following theorem. Theorem 3.1. When R0 > 1, the system (2.3) has a unique endemic equilibrium P ¼ ðE ; I ; R ; N Þ. Here E ¼ aed ðAl N Þ, I ¼ 1a ðAl N Þ, and R ¼ ac ðAl N Þ.
4. Local stability of the endemic equilibrium P* Theorem 4.1. If R0 > 1, the endemic equilibrium P of the system (2.3) is locally asymptotically stable. Proof. The matrix of the linearization of the system (2.3) at the equilibrium P ¼ ðE ; I ; R ; N Þ is 0 1 aI dx aI aI x hðN hðN hðN q Þ Þ Þ e B e d 0 0 C C; J ¼ B @ 0 c 1 0 A 0 a 0 1
22
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
where
aI abI ðE þ I þ R Þ 1 pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ q ¼ 1 þ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h2 ðN Þ 1 þ 2bN hðN Þ 1 þ 2bN aI 2abI ðE þ I þ R Þ 6 þ h2 ðN Þ hðN Þ aI 2aI bN 6 þ hðN Þ hðN Þ hðN Þ 3aI : 6 hðN Þ The characteristic equation of J is
ðk þ 1Þðk3 þ a1 k2 þ a2 k þ a3 Þ ¼ 0; where
aI a1 ¼ ðd þ xÞ þ 1 þ > 0; hðN Þ aI a2 ¼ ðd þ xÞ 1 þ >0 hðN Þ
and a3 ¼ Then
aI ðd þ ec þ eÞ þ aeq > 0: hðN Þ
2 aI aI 3aI aI ae ðd þ ec þ eÞ a1 a2 a3 P ðd þ xÞ 1 þ þ ðd þ xÞ 1 þ hðN Þ hðN Þ hðN Þ hðN Þ 2 aI aI 2 ¼ ðd þ xÞ þ ðd þ xÞ 1 þ þ ½ðd þ xÞ2 3ae d e ec hðN Þ hðN Þ 2 aI aI 2 P ðd þ xÞ þ ðd þ xÞ 1 þ þ ð4dx 3ae d e ecÞ hðN Þ hðN Þ 2 aI aI 2 ð3dx 2aeÞ þ ¼ ðd þ xÞ þ ðd þ xÞ 1 þ hðN Þ hðN Þ 2
P 0: Hence the Routh–Hurwitz conditions are satisfied. Thus it follows that the endemic equilibrium P of (2.3), which exists if R0 > 1, is always locally asymptotically stable. The following results, which may be found in [21], are used in the next section for the study of the global stability of the endemic equilibrium P . Consider the following systems x0 ¼ f ðt; xÞ; ð4:1Þ y 0 ¼ gðyÞ;
ð4:2Þ
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
23
where f and g are continuous and locally Lipschitz in x 2 Rn and solutions exist for all positive time. System (4.1) is called asymptotically autonomous with limit system (4.2) if f ðt; xÞ ! gðxÞ as t ! 1 uniformly for x 2 Rn . Lemma 4.1. Let e be a locally asymptotically stable equilibrium of (4.2) and x be the x-limit set of a forward bounded solution xðtÞ of (4.1). If x contains a point y0 such that the solution of (4.2), with yð0Þ ¼ y0 converges to e as t ! 1, then x ¼ feg, i.e., xðtÞ ! e as t ! 1. Corollary 4.1. If solutions of the system (4.1) are bounded and the equilibrium e of the limit system (4.2) is globally asymptotically stable, then any solution xðtÞ of the system (4.1) satisfies xðtÞ ! e as t ! 1. 5. Global stability of the endemic equilibrium P* For the system (2.2) equal to the system (2.3), making the change of variable c aþc M ¼ ðS þ E þ IÞ þ R; a a then dM ¼ cA M so that the following system (5.2) is equivalent to (2.2). ds al 8 dS A aSI ¼ l hðNÞ S; > ds > > dE < aSI ¼ hðN xE; ds Þ dI > ¼ eE dI; > ds > : dM ¼ Ac M; al ds
ð5:1Þ
ð5:2Þ
a ðS þ E þ I þ MÞ, x ¼ 1 þ e, and d ¼ 1 þ a þ c. The last equation of the system where N ¼ aþc (5.2) implies that Ac Ac s þ Mð0Þ ð5:3Þ MðsÞ ¼ e : al al
Substituting (5.3) into the first three equations of the system (5.2) gives the following threedimensional non-autonomous differential system 8 dS A aSI ¼ l hðNÞ S; > ds > > > dE < ¼ aSI xE; ds hðN Þ ð5:4Þ dI ¼ eE > ds > dI; > > : N ¼ a S þ E þ I þ Ac þ Mð0Þ Ac es : aþc
al
al
a ðS þ E þ I þ Ac Þ,X as s ! þ1, the system (5.4) is a three-dimensional asympSince N ! aþc al totically autonomous differential system with limit system 8 dS A aSI > < ds ¼ l hðX Þ S; dE aSI ¼ hðX xE; ð5:5Þ ds Þ > : dI ¼ eE dI: ds
24
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
a Here X ¼ aþc ðS þ E þ I þ Ac Þ. It can be verified that the region al
T ¼ fðS; E; IÞ 2 Rþ 3 : 0 6 S þ E þ I 6 A=lg is positively invariant with respect to (5.5). Here R3þ denotes the non-negative cone of R3 including its lower dimensional faces. Thus, the system (5.5) is a bounded system. The main aim of this section is to prove the following Theorem 5.1. Once Theorem 5.1 is proved, together with Theorem 2.1, we can claim that the basic reproduction number R0 is a sharp threshold parameter and that the global dynamical behaviors of the system (2.3) and the outcome of the disease are completely determined. In other words, if R0 6 1, the disease-free equilibrium P0 is globally asymptotically stable so that the disease dies out, while if R0 > 1, the endemic equiC so that the disease remains endemic. librium P is globally asymptotically stable in Theorem 5.1. When R0 > 1, the unique endemic equilibrium P of the system (2.3) is globally asymptotically stable in the interior of C. Moreover, P attracts all trajectories in C except those starting with the invariant E ¼ I ¼ 0 which converge to P0 along this hyperplane. By inspecting the vector field given by (2.3), we see that all trajectories starting from the boundary oC of C enter the interior C of C except those on the N -axis which converge to P0 along this invariant axis. Thus we need only to prove that C is the attractive region of P . In accordance with Corollary 4.1, it is sufficient to prove that the equilibrium Q of (5.5), which is corresponding to the equilibrium P of (2.3), is globally asymptotically stable in the interior of the region T . For this purpose, we need to study the dynamics of the system (5.5). Recalled that R0 is defined by (2.5), first of all, it is easy to obtain the following result: Lemma 5.1. If R0 6 1, Q0 ¼ ðA=l; 0; 0Þ is the only equilibrium in T and it is globally asymptotically stable. If R0 > 1, Q0 becomes unstable while Q ¼ ðS ; E ; I Þ, which corresponds to the equilibrium P of the system (2.3), emerges as a unique equilibrium in the interior of T and it is locally asA N E I cA , E ¼ de I , I ¼ al 1a N , and N is the unique ymptotically stable, here S ¼ aþc a al root of (3.4). Lemma 5.2. (1) When R0 > 1, Q0 is the only omega limit point of the system (5.5) on the boundary of T and cannot be the omega limit point of any orbit which starts in the interior of T , of the system (5.5). (2) The system (5.5) is uniformly persistent when R0 > 1. Proof. It is easy to see that the vector field of (5.5) is transversal to the boundary of T on all its faces except the S-axis which is invariant with respect to (5.5). On the S-axis the equation for S is dS ¼ Al S, which implies that SðsÞ ! A=l as s ! þ1. Therefore, Q0 is the only x-limit point on ds the boundary of T . Let L ¼ eE þ xI. Its derivative along the solutions of (5.5) is dL aeSI S hðA=lÞ ¼ dxI ¼ dxI R0 1 : ds hðX Þ A=l hðX Þ > 0 for S sufficiently close to A=l except when I ¼ 0. Hence when R0 > 1, the If R0 > 1, then dL ds trajectories starting from the points which are close enough to Q0 must leave a neighborhood of Q0 except those on the invariant S-axis. This proves the first conclusion of the lemma.
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
25
System (5.5) is said to be uniformly persistent (see [3,22]) if there exists a constant 0 < c < 1 such that any solution ðSðsÞ; EðsÞ; IðsÞÞ with ðSð0Þ; Eð0Þ; Ið0ÞÞ 2 T satisfies lim inf jðSðtÞ; EðtÞ; IðtÞÞj P c: t!1
The second conclusion of the lemma follows from a uniform persistence result, Theorem 4.3, in [23]. Since the maximal invariant set on the boundary oT of T is the S-axis which is isolated, the hypothesis ðH Þ of [23] holds for (5.5). The condition for uniform persistence in Theorem 4.3 of [23] is equivalent to Q0 being unstable. Thus, when R0 > 1, the system (5.5) is uniformly persistent in T. The definitions and results on competitive systems, which will be used in the following lemma, can be found in Appendix B. Lemma 5.3. The system (5.5) is competitive in the convex region T . Also any non-empty compact omega limit set of (5.5) in the interior of T is either a closed orbit or the equilibrium Q . Proof. The Jacobian matrix of (5.5) at a point P ¼ ðS; E; IÞ 2 T is 0 1 0 ðX Þ 0 ðX Þ aSIh0 ðX Þ a aI a aS a 1 hðX þ aSIh hðX þ aSIh Þ Þ h2 ðX Þ aþc h2 ðX Þ aþc h2 ðX Þ aþc B C 0 0 0 aI aS J ¼@ aSIh2 ðX Þ a x aSIh2 ðX Þ a aSIh2 ðX Þ a A; hðX Þ
h ðX Þ
0
aþc
h ðX Þ
e
aþc
hðX Þ
h ðX Þ
aþc
ð5:6Þ
d
b ffi where h0 ðX Þ ¼ b þ pffiffiffiffiffiffiffiffiffi > 0. 1þ2bX pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi a Note that X ¼ aþc ðS þ E þ I þ Ac Þ and hðX Þ ¼ 1 þ bX þ 1 þ 2bX , then al
aI aSIh0 ðX Þ a aI aIXh0 ðX Þ aI pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > 0; 2 P 2 ¼ hðX Þ h ðX Þ a þ c hðX Þ h ðX Þ hðX Þ 1 þ 2bX
ð5:7Þ
aS aSIh0 ðX Þ a aS aSXh0 ðX Þ aS pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi > 0: 2 P 2 ¼ hðX Þ h ðX Þ a þ c hðX Þ h ðX Þ hðX Þ 1 þ 2bX
ð5:8Þ
Choosing the matrix H as H ¼ diagð1; 1; 1Þ, we can see that the system (5.5) is competitive in the convex region T , with respect to the partial ordering defined by the orthant fðS; E; IÞ 2 R3 : S 6 0; E P 0; I 6 0g. Suppose that X is an omega limit set of (5.5) in the interior of T . If X does not contain Q , since Q is the only interior equilibrium, then it contains no equilibria. Theorem B.1 implies that X is a closed orbit. Suppose that X contains Q , since Q is locally asymptotically stable whenever it exists in the interior of T , any orbit that gets arbitrarily close to Q must converge to Q . Thus X ¼ Q . Lemma 5.3 is proved. From Lemma 5.3, in the absence of closed orbits, all trajectories in T converge to Q when R0 > 1. Thus, the key to verifying the global stability of Q is to rule out the existence of periodic solutions of (5.5). This is achieved by showing that any periodic solution to (5.5) is orbitally asymptotically stable. The following lemma is the main tool to prove the global stability of the equilibrium Q . Its proof depends on the criterion for the asymptotical orbital stability of period
26
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
orbits given by Muldowney [19]. For the sake of convenience, some definitions and conclusions about the asymptotical orbital stability of period solutions and compound matrices which are relevant to this lemma give are shown in Appendix A. Lemma 5.4. The orbit of any non-constant periodic solution to the system (5.5), if it exists, is asymptotically orbitally stable with asymptotic phase. Proof. Suppose that pðsÞ ¼ ðSðsÞ; EðsÞ; IðsÞÞ is a periodic solution to (5.5) with period T > 0, whose orbit q is contained in the interior of T . Obviously, SðsÞ; EðsÞ, and IðsÞ are all bounded functions and have positive lower boundary for all s P 0 from Lemma 5.2. The Jacobian matrix of (5.5) concerned with this periodic orbit q is given by (5.6). The second compound system, whose definition can be found in Appendix A, with respect to the periodic solution pðsÞ ¼ ðSðsÞ; EðsÞ; IðsÞÞ of (5.5) is the following 8 aSIh0 ðX Þ a dU aI aS > ¼ 1 þ x þ U þ ðV þ W Þ; 2 > hðX Þ hðX Þ h ðX Þ aþc > < ds 0 0 ðX Þ a ðX Þ a dV aI ð5:9Þ ¼ eU 1 þ d þ hðX aSIh W; V þ aSIh ds Þ h2 ðX Þ aþc h2 ðX Þ aþc > > 0 0 > : dW ¼ aI aSIh ðX Þ a V d þ x þ aSIh ðX Þ a W : ds hðX Þ h2 ðX Þ aþc h2 ðX Þ aþc Using Theorem A.1 shown in Appendix A, Lemma 5.4 will follow if we can prove that the system (5.9) is asymptotically stable. For this purpose, we shall use the following Lyapunov function EðsÞ ðjV ðsÞj þ jW ðsÞjÞ : LðsÞ ¼ LðUðsÞ; V ðsÞ; W ðsÞ; SðsÞ; EðsÞ; IðsÞÞ ¼ sup jU ðsÞj; IðsÞ From Lemma 5.2, we obtain that the orbit q of pðsÞ remains at a positive distance from the boundary of T , therefore EðsÞ P k
and IðsÞ P k
with 0 < k < 1:
Hence the function L is well defined along pðsÞ and there exists a constant c > 0 such that LðU ; V ; W ; S; E; IÞ P c supfjU j; jV j þ jW jg:
ð5:10Þ
The right-hand derivative of LðsÞ exists. Together with (5.7) and (5.8), direct calculation yields aI aS aSIh0 ðX Þ a Dþ jU ðsÞj 6 1 þ x þ jU ðsÞj þ 2 ðjV ðsÞj þ jW ðsÞjÞ hðX Þ hðX Þ h ðX Þ a þ c aI aS aSIh0 ðX Þ a I E ðjV ðsÞj þ jW ðsÞjÞ 6 1þxþ jU ðsÞj þ 2 hðX Þ hðX Þ h ðX Þ a þ c E I and
aI aSIh0 ðX Þ a aSIh0 ðX Þ a Dþ jV ðsÞj 6 ejU ðsÞj 1 þ d þ 2 jV ðsÞj þ 2 jW ðsÞj; hðX Þ h ðX Þ a þ c h ðX Þ a þ c aI aSIh0 ðX Þ a aSIh0 ðX Þ a 2 jV ðsÞj d þ x þ 2 jW ðsÞj: Dþ jW ðsÞj 6 hðX Þ h ðX Þ a þ c h ðX Þ a þ c
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32 dI Let E0 ¼ dE and I 0 ¼ ds , then ds
27
E0 EI 0 E 2 ðjV ðsÞj þ jW ðsÞjÞ þ Dþ ðjV ðsÞj þ jW ðsÞjÞ I I I 0 E I0 E eE 6 ðjV ðsÞj þ jW ðsÞjÞ þ jU ðsÞj I E I I E E ð1 þ dÞ jV ðsÞj ðd þ xÞ jW ðsÞj I I 0 0 eE E I E jU ðsÞj þ ðjV ðsÞj þ jW ðsÞjÞ: 6 1d I I E I
E Dþ ðjV ðsÞj þ jW ðsÞjÞ ¼ I
Thus, we can obtain Dþ LðsÞ 6 supfg1 ðsÞ; g2 ðsÞgLðsÞ;
ð5:11Þ
where aI aSI aSI 2 h0 ðX Þ a þ hðX Þ EhðX Þ Eh2 ðX Þ a þ c aSI 6 1xþ EhðX Þ
g1 ðsÞ ¼ 1 x
and eE E0 I 0 þ 1 d: I E I We rewrite the last two equations of (5.5) as g2 ðsÞ ¼
aSI E0 ¼ þ x; EhðX Þ E eE I 0 ¼ þ d; I I we find E0 1: E Thus, from Eq. (5.11) and GronwallÕs inequality, we obtain supfg1 ðsÞ; g2 ðsÞg 6
LðsÞ 6
Lð0ÞEðsÞ s ALð0Þ s e 6 e ; Eð0Þ lEð0Þ
which implies that LðsÞ ! 0 as s ! 1. By (5.10) it turns out that LðU ðsÞ; V ðsÞ; W ðsÞÞ ! 0
as s ! 1:
This implies that the system (5.9) is asymptotically stable and therefore the periodic solution pðsÞ is asymptotically orbitally stable with asymptotical phase by Theorem A.1 in Appendix A. This proves Lemma 5.4.
28
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
As mentioned above, by Lemma 5.4, in the absence of closed orbits, all trajectories in T converge to Q when R0 > 1. Thus, we have the following theorem whose proof is similar to that of Theorem 2.1 in [16] and Theorem 6.5 in [17]. Theorem 5.2. If R0 > 1, the positive equilibrium Q of the system (5.5) is globally asymptotically stable in the interior of T . Thus, the main result of this section, Theorem 5.1, follows from Corollary 4.1. Therefore, the global dynamics of the SEIR model considered in this paper are completely determined. 6. Conclusions and discussion This paper has considered a SEIR model that incorporates recruitment and exponential natural death, as well as disease-related death, so that the total population size may vary in time. A distinguishing feature of the SEIR model considered here is that there is a saturating contact rate bN pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi : CðNÞ ¼ 1 þ bN þ 1 þ 2bN The basic reproduction number R0 defined by (2.5) of this SEIR model (2.2) is a sharp threshold parameter which completely determines the global dynamics of the system (2.2) and the outcome of the disease. If R0 6 1, the disease-free equilibrium is globally stable so that the disease always dies out, and if R0 > 1, the disease-free equilibrium becomes unstable while the endemic equilibrium emerges as the unique positive equilibrium and it is globally stable in the interior of the feasible region C. In the special case when the mean latent period e10 ! 0 or the mean infective period c10 ! 1, the SEIR model here reduces to an SIR model or an SEI model with recruitment and saturating contact rate CðNÞ, respectively. Thus, the global dynamic results of these two models can be obtained as special cases when e10 ¼ 0 and c0 ¼ 0 respectively. Especially, to prove the global stability of the endemic equilibrium, we make a change of variable by which our fourdimensional model can be reduced to a three-dimensional asymptotical autonomous system with limit system. The similar change of variable M ¼ ac N þ R can reduce the SIR model with a general contact rate bðN Þ (the incidence is bðN ÞSI) in [4] by Brauer and van den Driessche to the following two-dimensional asymptotical autonomous differential system dI aþc cA cA dt ða þ c þ dÞI ð6:1Þ ¼ pA þ bðNÞI N I Mð0Þ e dt a ad ad and
dN dt
¼ A dN aI with the following limit system. dI aþc cA ¼ pA þ bðNÞI N I ða þ c þ dÞI dt a ad
ð6:2Þ
¼ A dN aI. The Dulac function BðI; N Þ ¼ 1I can exclude the existence of period soluand dN dt tions to (6.2). The global stability of the endemic equilibrium follows from Corollary 4.1 and the Poincare–Bendixson theorem. This change of variable can also reduce the proof of the global stability of the endemic equilibrium of the system (5) in [4]. Additionally, because of the complex pffiffiffiffiffiffiffiffiffiffi used in our SEIR model, many mathematical techniques structure of the incidence rate 1þbN þbbSI 1þ2bN
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
29
are used in the present paper to overcome mathematical difficulties that are not present in the SEIR model with bilinear incidence or standard incidence. The saturating contact rate CðN Þ defined by (1.1) is increasing, approximately a linear function bN for small N , and asymptotically approaches a constant value 1 as N ! 1. Using notation from Section 2, the SEIR model with recruitment and the standard incidence b NS I or the bilinear incidence bbN NS I ¼ bbSI can be respectively described by the following systems dS A bSI ¼ S; ds l N dE bSI ¼ xE; ds N dI ¼ eE dI; ds dR ¼ cI R ds
ð6:3Þ
and dS A ¼ bbSI S; ds l dE ¼ bbSI xE; ds ð6:4Þ dI ¼ eE dI; ds dR ¼ cI R: ds Thus both (6.3) and (6.4) are two extreme cases of the SEIR model in this paper. Noting that R0 ¼
aAe be0 bA=l pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ; ¼ ldxhðA=lÞ ðl þ e0 Þðl þ a0 þ c0 Þ ð1 þ bA=l þ 1 þ 2bA=lÞ
be0 ,R01 P R0 and if N is very if letting N ! 1 (then, Al ! 1), since CðN Þ ! 1, R0 ! ðlþe0 Þðlþa 0 þc0 Þ small, since CðN Þ is approximately a linear function bN , R0 is approximately be0 bA ,R02 P R0 . It is easy to verify that R01 and R02 are the basic reproduction numbers ðlþe0 Þðlþa0 þc0 Þ l I ; d I ; I ; cI1 Þ and ðAl dx I ; d I ; I ; cI2 Þ are the of the system (6.3) and (6.4), respectively. ðAl dx e 1 e 1 1 e 2 e 2 2 endemic equilibrium of (6.3) when R01 > 1 and the endemic equilibrium (6.4) when R02 > 1, rea ae A , I2 ¼ 1a ðA=l N2 Þ, and N2 ¼ bb þ ð1 dx Þ l. spectively, where I1 ¼ 1a ðA=l N1 Þ, N1 ¼ bðdxaeÞA dxðbaÞl Substituting N1 and N2 into (3.4) gives pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi F ðN1 Þ ¼ adxð1 þ 1 þ 2bN1 Þ < 0
and
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi F ðN2 Þ ¼ abðdx aeÞA=l a2 dx=b adx 1 þ 2bN2 < 0:
Thus, N1 < N and N2 < N , then I1 > I and I2 > I . Thus, although the threshold behavior and dynamic behavior of the SEIR model (2.2) in this paper are similar to those of the SEIR model (6.3) and (6.4), the basic reproduction number R0 of (2.2) is less than R01 and R02 , the endemic
30
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
equilibrium population size N is grater than those of (6.3) and (6.4), and the infective level of the endemic equilibrium state is less than those of (6.3) and (6.4). The epidemiological and demographical phenomena as well as their interaction can be observed in our model. For example, when the disease persists, the disease-related deaths can reduce the equilibrium population size below its maximum carrying capacity. The model considered here cannot give rise to periodic solutions. Acknowledgements This research is supported by the National Natural Sciences Foundation of PeopleÕs Republic of China grant no. 19971066. We also wish to thank two anonymous referees for helpful criticism and suggestions. Appendix A Let D Rn be an open set, and gðyÞ 2 Rn be a C 1 function defined in D. We consider the autonomous system in Rn y 0 ¼ gðyÞ: ðA:1Þ Let yðt; y0 Þ denote the solution of (A.1) such that yð0; y0 Þ ¼ y0 . The linear variational equation of (A.1) with respect to yðt; y0 Þ is given by og ðA:2Þ Y 0 ðtÞ ¼ ðyðt; y0 ÞÞY ðtÞ; oy where og is the Jacobian matrix of g. Suppose the autonomous system (A.1) has a periodic solution oy y ¼ pðtÞ with least period T > 0 and orbit q ¼ fpðtÞ : 0 6 t 6 T g. This orbit is orbitally stable if for each > 0, there exists a r > 0 such that any solution yðtÞ, for which the distance of yð0Þ from q is less than r remains at a distance less than from q for all t P 0. It is asymptotically orbitally stable if the distance of yðtÞ from q also tends to zero as t ! 1. This orbit q is asymptotically orbitally stable with asymptotic phase if it is asymptotically orbitally stable and there is a h > 0 such that any solution yðtÞ for which the distance of yð0Þ from q is less than h, satisfies jyðtÞ pðt t Þj ! 0 as t ! 1 for some t which may depend on yð0Þ. The following is a criterion given in [19] for the asymptotical orbital stability of a periodic orbit to the general autonomous system (A.1). Theorem A.1. A sufficient condition for a period orbit q ¼ fpðtÞ : 0 6 t 6 T g of (A.1) to be asymptotically orbitally stable with asymptotic phase is that the linear system ! ½2 og Z 0 ðtÞ ¼ ðpðtÞÞ ZðtÞ ðA:3Þ oy be asymptotically stable. ½2
Eq. (A.3) is called the second compound equation of (A.2) and the matrix og ðpðtÞÞ is the second oy at y ¼ pðtÞ of g. Generally speaking, for a n n matrix compound matrix of the Jacobian matrix og oy
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
31
A and integer1 6k 6 n, the kth additive compound matrix of A is denoted by A½k . This is a N N n matrix, N ¼ , defined by k A½k ¼ Dþ ðI þ hAÞðkÞ jh¼0 ; where BðkÞ is the kth exterior power of a n n matrix B and Dþ denotes the right-hand derivative. A survey on the definition and properties of additive compound matrices together with their connections to differential equations may be found in [19,24]. Further applications may be found in [25,26]. The following are the examples of various compound matrices when n ¼ 3. 0 1 a11 a12 a13 A½1 ¼ A ¼ @ a21 a22 a23 A; a31 a32 a33 0 1 a11 þ a22 a23 a13 a11 þ a33 a12 A A½2 ¼ @ a32 a31 a21 a22 þ a33 and A½3 ¼ a11 þ a22 þ a33 ¼ tr A:
Appendix B All concepts and results in this appendix can be found in [16,17,27,28]. The autonomous system (A.1) is said to be competitive in D if, for some diagonal matrix H has non-positive off-diagonal elements H ¼ diagðr1 ; . . . ; rn Þ where each ri is either 1 or )1, H og oy for all y 2 D. If D is convex, the flow of such a competitive system preserves, for t < 0, the partial ordering in Rn defined by the orthant K ¼ fðy1 ; . . . ; yn Þ 2 Rn : ri yi P 0; i ¼ 1; 2; . . . ; ng. The concept of competitiveness used above is more general than the one in [28] in that the partial ordering is not necessarily defined by the standard orthant of Rn . However, by a linear change of variables z ¼ Hy, a competitive system as defined above can be transformed into a system that is ÔcompetitiveÕ in the sense of [28]. An important characteristic of a three-dimensional competitive system is the following Poincare–Bendixson property. Theorem B.1. Suppose that n ¼ 3, D is convex, (A.1) is competitive in D and X is a non-empty compact omega limit set of (A.1). If X contains no equilibria, then X must be a closed orbit (cf. [28] Theorem 1). References [1] J.A.P. Heesterbeek, J.A.J. Metz, The saturating contact rate in marriage and epidemic models, J. Math. Biol. 31 (1993) 529. [2] H.R. Thieme, C. Castillo-Chavez, On the role of variable infectivity in the dynamics of the human immunodeficiency virus epidemic, in: C. Castillo-Chavez (Ed.), Mathematical and Statistical Approaches to AIDS Epidemiology, Lecture Notes in Biomathematics, vol. 83, Springer, Berlin, 1989, p. 157.
32
J. Zhang, Z. Ma / Mathematical Biosciences 185 (2003) 15–32
[3] H.R. Thieme, Epidemic and demographic interaction in the spread of potentially fatal diseases in growing populations, Math. Biosci. 111 (1992) 99. [4] F. Brauer, P. van den Driessche, Models for transmission of disease with immigration of infectives, Math. Biosci. 171 (2001) 143. [5] J. Mena-Lorca, H.W. Hethcote, Dynamic models of infectious diseases as regulators of population sizes, J. Math. Biol. 30 (1992) 693. [6] F. Brauer, Epidemic models in populations of varying size, in: C. Castillo-Chavez, S.A. Levin, C. Shoemaker (Eds.), Mathematical Approaches to Problems in Resource Management and Epidemiology, Lecture Notes in Biomathematics, vol. 81, Springer, p. 109, 1989. [7] L.Q. Gao, H.W. Hethcote, Disease transmission models with density-dependent demographics, J. Math. Biol. 30 (1992) 717. [8] D. Greenhalgh, Some threshold and stability results for epidemic models with a density dependent death rate, Theor. Pop. Biol. 42 (1992) 130. [9] H.J. Bremermann, H.R. Thieme, A competitive exclusion principle for pathogen virulence, J. Math. Biol. 27 (1989) 179. [10] H.W. Hethcote, P. van den Driessche, Some epidemiological models with nonlinear incidence, J. Math. Biol. 29 (1991) 271. [11] W.R. Derrick, P. van den Driessche, A disease transmission model in a non-constant population, J. Math. Biol. 31 (1993) 495. [12] S. Busenberg, P. van den Driessche, Analysis of a disease transmission model with varying population size, J. Math. Biol. 29 (1990) 257. [13] D. Greenhalgh, Some results for a SEIR epidemic model with density dependence in the death rate, IMA J. Math. Appl. Med. Biol. 9 (1992) 67. [14] K. Cooke, P. van den Driessche, Analysis of an SEIRS epidemic model with two delays, J. Math. Biol. 35 (1996) 240. [15] D. Greenhalgh, Hopf bifurcation in epidemic models with a latent period and nonpermanent immunity, Math. Comput. Model. 25 (1997) 85. [16] M.Y. Li, J.S. Muldowney, Global stability for the SEIR model in epidemiology, Math. Biosci. 125 (1995) 155. [17] M.Y. Li, J.R. Graef, L. Wang, J. Karsai, Global dynamics of a SEIR model with varying total population size, Math. Biosci. 160 (1999) 191. [18] M.Y. Li, H.L. Smith, L. Wang, Global dynamics of an SEIR epidemic model with vertical transmission, SIAM J. Appl. Math. 62 (2001) 58. [19] J.S. Muldowney, Compound matrices and ordinary differential equations, Rocky Mount. J. Math. 20 (1990) 857. [20] J.K. Hale, Ordinary Differential Equations, John Wiley, New York, 1969. [21] H.R. Thieme, Convergence results and a Poincare–Bendixson trichotomy for asymptotically automous differential equations, J. Math. Biol. 30 (1992) 755. [22] G. Butler, P. Waltman, Persistence in dynamical systems, Proc. Am. Math. Soc. 96 (1986) 425. [23] H.I. Freedman, S.G. Ruan, M.X. Tang, Uniform persistence and flows near a closed positively invariant set, J. Dynam. Diff. Equat. 6 (1994) 583. [24] D. London, On derivations arising in differential equations, Linear Multilin. Algeb. 4 (1976) 179. [25] Y. Li, J.S. Muldowney, Evolution of surface functionals and differential equations, in: J. Wiener, J. Hale (Eds.), Ordinary and Delay Differential Equations, Pitman Research Notes in Mathematics, vol. 272, Longman, England, 1992, p. 144. [26] Y. Li, J.S. Muldowney, On BendixsonÕs criterion, J. Diff. Equat. 106 (1993) 27. [27] H.L. Smith, Systems of ordinary differential equations which generate an order preserving flow. A survey of results, SIAM Rev. 30 (1988) 87. [28] M.W. Hirsch, Systems of differential equations which are competitive or cooperative. IV: Structural stabilities in three dimensional systems, SIAM J. Math. Anal. 21 (1990) 1225.