- Email: [email protected]
Spread of disease in a patchy environment
Chaos, Solitons & Fractals 45 (2012) 942–949
Contents lists available at SciVerse ScienceDirect
Chaos, Solitons & Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos
Spread of disease in a patchy environment Junli Liu School of Science, Xi’an Polytechnic University, Xi’an 710048, PR China
a r t i c l e
i n f o
Article history: Received 12 February 2011 Accepted 19 March 2012 Available online 28 April 2012
a b s t r a c t For a two patches SIR model, it is shown that its dynamic behavior is determined by several quantities. We have shown that if R0 < 1, then the disease-free equilibrium is globally asymptotically stable, otherwise it is unstable. Some sufficient conditions for the local stability of boundary equilibria are obtained. Numerical simulations indicate that travel between patches can reduces oscillations in both patches; may enhances oscillations in both patches; or travel switches oscillations from one patch to another. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction An SIR model was proposed by Cooke [1] for epidemics which spread in a human population via a vector. Susceptible individuals receive the infection from infectious vectors, and susceptible vectors receive the infection from infectious individuals. It is assumed that when a susceptible vector is infected by a person, there is a time s > 0 during which the infectious agents develop in the vector and it is only after that time that the infected vector becomes itself infectious. It is also assumed that the vector population is very large and at any time t the infectious vector population is simply proportional to the infectious human population at time t s. Thus, if we denote S(t), I(t), and R(t) as the numbers of the population susceptible to the disease, of infective members, and of members who have been removed from the possibility of infection through full immunity, respectively. Assume that the recruitment into the susceptible class is A. Let b be the transmission coefficient with mass action incidence, l be the natural death rate, c be the recovery rate of infectives, be the disease-reduced death rate. Then we obtain
8 _ > < SðtÞ ¼ A bSðtÞIðt sÞ lSðtÞ; _IðtÞ ¼ bSðtÞIðt sÞ ðl þ c þ ÞIðtÞ; > :_ RðtÞ ¼ cIðtÞ lRðtÞ:
ð1:1Þ
E-mail address: [email protected] 0960-0779/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.chaos.2012.03.007
where S(t) + I(t) + R(t) = N(t) denotes the number of a population at time t. For the detailed biological meanings, we refer to [1,4–7]. System (1.1) has one disease-free equilibrium lA ; 0; 0 , since lþ1cþ is the average infection period, we can see that e 0 ¼ bA R lðlþcþÞ is a reproduction number of the disease. Stability of (1.1) was studied in [4,5,8–10]. It is shown that, e 0 < 1, the disease-free equilibrium P0 A ; 0; 0 is when R l
globally asymptotically stable for any time delay s. If e 0 > 1; P 0 becomes unstable and there is a positive equiR lS Þ librium P S ; AbSlS ; cðA with S ¼ lþbcþ, and P⁄ is locally lbS asymptotically stable for any time delay s. As for the global asymptotic stability on P⁄, some sufficient conditions were obtained in [5,9]. It is also shown that system (1.1) is pere 0 > 1 [10,11]. manent if R System (1.1) assumes that the mix of population is homogeneous. Because of the dispersal among populations, infection will be carried to many places, such as the SARS outbreak of 2002–2003. Much has been done in terms of modeling spatial spread of diseases. For example, Wang and Ruan [12] studied the global spread pattern of 2002– 03 SARS outbreak. Sattenspiel and Herring [13] used a model to the consideration of travel between populations in the Canadian subarctic. Sattenspiel and Dietz [14] studied a model for the spread of infectious diseases incorporating a mobility process among discrete geographic regions. In [15], Colizza and Vespignani used the basic reaction– diffusion equations to describe the metapopulation system
943
J. Liu / Chaos, Solitons & Fractals 45 (2012) 942–949
at the mechanistic level and derive an early stage dynamics approximation for the subpopulation invasion dynamics. In [16], Colizza et al. considered the problem of the predictability of epidemic scenarios, within the global stochastic modeling framework. By the level of sampling of the underlying travel pattern by infected and latent individuals, the overall predictability can be determined. Bianco and Shaw studied the stability of a multistrain disease model in two coupled populations in [19]. The metapopulation dynamics of infectious diseases has generated a lot of models and results considering population moved among regions with endemic disease [17,18,20–22]. To investigate the spatial effects, in the next section, we will construct a two patches SIR model, and assume the susceptible and recovered individuals will travel between patches. We discuss the existence of equilibria in Section 3, and show that R0 < 1 implies the global stability of the disease-free equilibrium. In Section 4, we consider the local stability of boundary equilibria and endemic equilibrium, and we show the permanence of the model. We conclude with some numerical simulations and discussions in Section 5 about the role of travel in a disease outbreak.
Proof. From system (2.1), the total population N(t) satis_ fies NðtÞ ¼ A1 þ A2 l N 1 ðtÞ l N 2 ðtÞ 1 I1 ðtÞ 2 I2 ðtÞ. 1
2
It follows form I1(t) < N1(t), I2(t) < N2(t) that A1 þ A2 _ 6 A1 þ A2 minfl1 ; l2 g maxfl1 þ 1 ; l2 þ 2 gNðtÞ 6 NðtÞ h n o 2 NðtÞ. Hence, NðtÞ 2 min Nð0Þ; maxflAþ1 þA ; max 1 1 ;l2 þ2 g n oi A1 þA2 Nð0Þ; minfl ;l g for t P 0. This lemma is complete. 1
2
3. Existence of equilibria and stability of the disease-free equilibrium It is clear that system (2.1) always has a disease-free equilibrium E0(S01, 0, 0, S02, 0, 0) for all s P 0, with S01 ¼ A1 l2 þðA1 þA2 Þa2 l1 l2 þl1 a2 þa1 l2 ;
1 þA2 Þa1 S02 ¼ lA2ll1þþðA l a2 þa1 l . As variables R1, R2 do 1
2
1
2
not appear in the first, second, fourth and fifth equations in (2.1). Thus, for local stability, we consider the system in Si and Ii for i = 1, 2. Linearizing these equations about E0 gives the characteristic equation
½ðz þ l1 þ a1 Þðz þ l2 þ a2 Þ a1 a2 gðzÞ ¼ 0;
ð3:1Þ
where 2. Model formulation
gðzÞ ¼ ðz þ l1 þ c1 þ 1 b1 S01 ezs Þðz þ l2 þ c2
In this section, we will construct a model with two patches, for each patch, consider an SIR model as (1.1). We also assume that the disease renders infectives too sick to travel. Let ai and ci denote the travel rates of susceptible individuals and recovered individuals from patch i to patch j, i, j = 1, 2, respectively. The dynamics is governed by the delay differential system
8 > S_ ðtÞ ¼ A1 b1 S1 ðtÞI1 ðt sÞ l1 S1 ðtÞ þ a2 S2 ðtÞ a1 S1 ðtÞ; > > 1 > > _ > > > I1 ðtÞ ¼ b1 S1 ðtÞI1 ðt sÞ ðl1 þ c1 þ 1 ÞI1 ðtÞ; > > < R_ ðtÞ ¼ c I ðtÞ l R ðtÞ þ c R ðtÞ c R ðtÞ; 1 2 2 1 1 1 1 1 1 _ 2 ðtÞ ¼ A2 b S2 ðtÞI2 ðt sÞ l S2 ðtÞ þ a1 S1 ðtÞ a2 S2 ðtÞ; > S > 2 2 > > > > I_ ðtÞ ¼ b S ðtÞI ðt sÞ ðl þ c þ ÞI ðtÞ; > > 2 2 2 2 2 2 2 2 > > :_ R2 ðtÞ ¼ c2 I2 ðtÞ l2 R2 ðtÞ þ c1 R1 ðtÞ c2 R2 ðtÞ: ð2:1Þ
Lemma 2.1. The solution of (2.1) with an initial condition (2.2) is nonnegative for all t P 0. Lemma 2.2. The total population N(t) of system (2.1) is h n o 1 þA2 ; bounded and is in the interval min Nð0Þ; maxflAþ 1 1 ;l2 þ2 g n o A1 þA2 max Nð0Þ; minfl ;l g . 2
Clearly, (z + l1 + a1)(z + l2 + a2) a1a2 = 0 admits zeros with negative real parts only. Therefore, the stability of E0 is determined by the distribution of roots of g(z) = 0. For s = 0, g(z) = 0 reduces to
ðz þ l1 þ c1 þ 1 b1 S01 Þðz þ l2 þ c2 þ 2 b2 S02 Þ ¼ 0: ð3:3Þ Then when s = 0, E0 is locally stable if and only if b1S01 < l1 + c1 + 1 and b2S02 < l2 + c2 + 2. By Theorem 2.1 of [2], we have that if b1S01 < l1 + c1 + 1 and b2S02 < l2 + c2 + 2, then E0 is locally asymptotically stable for all s P 0; if b1S01 > l1 + c1 + 1 or b2S02 > l2 + c2 + 2, then E0 is unstable for all s P 0. e 01 ¼ b1 S01 ; R e 02 ¼ b2 S02 ; R0 ¼ maxf R e 01 ; R e 02 g. Define R l þc þ1 l þc þ2 1
1
2
2
ð2:2Þ
where R6þ0 ¼ fðx1 ; x2 ; x3 ; x4 ; x5 ; x6 Þ 2 R6 : xi P 0; i ¼ 1; 2; 3; 4; 5; 6g. The total population in each patch is Ni(t) = Si(t) + Ii (t) + Ri(t), i = 1, 2, then the overall total population is N(t) = N1(t) + N2(t). It is straightforward to show that.
1
ð3:2Þ
Theorem 3.1. If R0 < 1, then E0 is globally asymptotically stable for all s P 0, whereas if R0 > 1, E0 is unstable for all s P 0.
The initial conditions for system (2.1) are
ð/1 ðhÞ; /2 ðhÞ; /3 ðhÞ; /4 ðhÞ; /5 ðhÞ; /6 ðhÞÞ 2 C ¼ C ½s; 0; R6þ0 ;
þ 2 b2 S02 ezs Þ:
Proof. From the above analysis, it is left to show that E0 is globally attractive if R0 < 1. From system (2.1) and Lemma 2.1, we get
(
S_ 1 ðtÞ 6 A1 l1 S1 ðtÞ þ a2 S2 ðtÞ a1 S1 ðtÞ; S_ 2 ðtÞ 6 A2 l S2 ðtÞ þ a1 S1 ðtÞ a2 S2 ðtÞ:
ð3:4Þ
2
Consider an auxiliary system
u_ 1 ðtÞ ¼ A1 l1 u1 ðtÞ þ a2 u2 ðtÞ a1 u1 ðtÞ; u_ 2 ðtÞ ¼ A2 l2 u2 ðtÞ þ a1 u1 ðtÞ a2 u2 ðtÞ:
ð3:5Þ
It is easy to show that (3.5) has a globally asymptotically stable equilibrium u1 ; u2 ¼ ðS01 ; S02 Þ, that is, limt?1 u1(t) = S01, limt?1u2(t) = S02. By standard comparison theorem, for sufficiently small e > 0, there is a t0 > 0 such that
944
J. Liu / Chaos, Solitons & Fractals 45 (2012) 942–949
S1(t) < S01 + e and S2(t) < S02 + e for all t P t0. For R0 < 1, we can choose e such that
b1 ðS01 þ eÞ < l1 þ c1 þ 1 ;
b2 ðS02 þ eÞ
< l2 þ c2 þ 2 :
e1; e e 2 Þ exists if and (i) Boundary equilibria E1 ðe S 1 ; eI 1 ; R S 2 ; 0; R e 01 > 1; and E2 ðb b1; b b 2 Þ exists if and only if R S 1 ; 0; R S 2 ; bI 2 ; R e 02 > 1. Here, only if R
l þ c1 þ 1 e l1 l2 þ l1 a2 þ a1 l2 e e ; I1 ¼ ð R 01 1Þ; S1 ¼ 1 b1 b1 ðl2 þ a2 Þ ðl2 þ c2 Þc1eI 1 b A2 þ a1 ðl1 þ c1 þ 1 Þ e1 ¼ ;e S ¼ 1 ; R ðl2 þ a2 Þb1 l1 l2 þ l1 c 2 þ c 1 l2 2 c1 c1eI 1 b A1 þ a2 ðl2 þ c2 þ 2 Þ e2 ¼ ;b S ¼ 2 ; R ðl1 þ a1 Þb2 l1 l2 þ l1 c 2 þ c 1 l2 1 e c 2 c2 I 2 l þ c2 þ 2 b1 ¼ ;b S ¼ 2 ; R l1 l2 þ l1 c 2 þ c 1 l2 2 b2 l l þ l a þ a l 1 2 e 1 2 bI 2 ¼ 1 2 ð R 02 1Þ; b2 ðl1 þ a1 Þ e ðl1 þ c1 Þc2 I 2 b2 ¼ : R l1 l2 þ l1 c 2 þ c 1 l2
ð3:6Þ
Then from (2.1), we get
(
I_1 ðtÞ 6 b1 ðS01 þ eÞI1 ðt sÞ ðl1 þ c1 þ 1 ÞI1 ðtÞ; I_2 ðtÞ 6 b ðS02 þ eÞI2 ðt sÞ ðl þ c þ 2 ÞI2 ðtÞ 2
2
ð3:7Þ
2
for t P t0. Consider another auxiliary equation
x_ ¼ b1 ðS01 þ eÞxðt sÞ ðl1 þ c1 þ 1 ÞxðtÞ:
ð3:8Þ
Then x(t) ? 0 as t ? 1 if b1(S01 + e) < l1 + c1 + 1, again by comparison theorem, we have I1(t) ? 0. Similarly, we can show that I2(t) ? 0 when b2(S02 + e) < l2 + c2 + 2. Applying the theory of asymptotically autonomous systems (see [3]) to system (2.1) yields R1(t) ? 0, R2(t) ? 0, S1(t) ? S01 and S2(t) ? S02 as t ? 1. Therefore, if R0 < 1, then E0 is globally asymptotically stable. This completes the proof. If we assume each patch is isolated, then the corresponding basic reproduction number in each patch is
R01 ¼
A1 b 1
l1 ðl1 þ c1 þ 1 Þ
;
R02 ¼
A2 b2
l2 ðl2 þ c2 þ 2 Þ
Corollary 3.2. If R0 < 1, then E0 is the only equilibrium of system (2.1).
l2 þ a2 þ AA2 a1 2 l1 þ a1 þ AA1 2a1 e R : l2 a1 R01 ; R 02 ¼ l2 þ a2 þ l1 l1 þ a1 þ ll1 a2 2 02
Proof. R0 < 1 implies that boundary equilibria E1 and E2 do not exist and
e 01 P R01 and R e 02 6 R02 . The reNote that if AA2 a1 2 P ll2 a1 , then R 1 e 02 P R02 and verse inequality gives AA2 1a2 6 ll2 a1 , then R 1 e 01 , if a1 < A2 a2 R01 e 01 6 R01 . Also from the definition of R R A1
e 01 > 1. Similar result holds for ð1 R01 Þðl2 þ a2 ÞÞ ll1 , then R 2 e R 02 . Therefore, we have the following relationship of R0 to R01 and R02. Corollary 3.1. If R01 > 1 and R02 > 1, then R0 > 1. If R01 < 1 and R02 < 1, then R0 > 1 provided one of the following holds
R01 ð1 R01 Þðl2 þ a2 Þ ll1 . (i) a1 < 2 (ii) a2 < AA1 2a1 R02 ð1 R02 Þðl1 þ a1 Þ ll2 . A2 a2 A1
1
Besides the disease-free equilibrium E0, the other equilibria are stated in the following theorem. First, we define some notations
R1;2 R2;1
l1 þ c1 þ 1
a1 þ l1 ; I1 ¼ ðR1;2 1Þ; b1 b1 ðl þ c2 Þc1 I1 þ c2 c2 I2 l þ c2 þ 2 ; S2 ¼ 2 ; R1 ¼ 2 l1 ðl2 þ c2 Þ þ c1 l2 b2 a2 þ l2 ðl þ c1 Þc2 I2 þ c1 c1 I1 ðR2;1 1Þ; R2 ¼ 1 : I2 ¼ b2 l1 ðl2 þ c2 Þ þ c1 l2
S1 ¼
:
e 01 and R e 02 can be written as Then R
e 01 ¼ R
(ii) System (2.1) has a unique endemic equilibrium E S1 ; I1 ; R1 ; S2 ; I2 ; R2 if and only if R1,2 > 1 and R2,1 > 1. Here,
A1 b1 b1 a2 ðl2 þ c2 þ 2 Þ ¼ þ ; ðl1 þ c1 þ 1 Þða1 þ l1 Þ b2 ðl1 þ c1 þ 1 Þða1 þ l1 Þ A2 b2 b2 a1 ðl1 þ c1 þ 1 Þ ¼ þ : ðl2 þ c2 þ 2 Þða2 þ l2 Þ b1 ðl2 þ c2 þ 2 Þða2 þ l2 Þ
S01 <
l1 þ c1 þ 1 b1
¼ S1 ;
S02 <
l2 þ c2 þ 2 b2
¼ S2 :
ð3:9Þ
From the first and the fourth equations of (2.1), we have
A1 þ A2 ¼ l1 S01 þ l2 S02 :
ð3:10Þ
If the endemic equilibrium E⁄ exists, it follows from adding the equations in (2.1) that
A1 þ A2 ¼ l1 S1 þ l2 S2 þ ðl1 þ 1 ÞI1 þ ðl2 þ 2 ÞI2 þ l1 R1 þ l2 R2 : This gives A1 þ A2 > l1 S1 þ l2 S2 . By (3.9) and (3.10), l1 S1 þ l2 S2 < l1 S01 þ l2 S02 < l1 S1 þ l2 S2 , which is a contradiction. Therefore, the endemic equilibrium E⁄ does not exist. Thus, if R0 < 1, the only equilibrium of (2.1) is E0.
4. Stability analysis and permanence In this section, we will discuss the stability of boundary equilibria and endemic equilibrium. 4.1. Local stability of E1 and E2
Theorem 3.2. The system (2.1) can have other three equilibria, and we have the following results:
Set d1 = l1 + c1 + 1, d2 = l2 + c2 + 2. Theorem 3.2 shows e 01 > 1, then E1 exists. Linearizing the Si, Ii, i = 1, that when R
945
J. Liu / Chaos, Solitons & Fractals 45 (2012) 942–949
2 equations in (2.1) about equilibrium E1 gives the characteristic equation
0. Based on the above discussion, we know that if e 01 > 1 þ l2 þa2 , (4.3) has no positive root, and so f(z) has R l þa1
S 2 ezs þ l2 þ c2 þ 2 Þf ðzÞ ¼ 0; ðz b2 e
no purely imaginary zero, thus E1 is locally asymptotically stable.
1
ð4:1Þ
where
l þa
S 1 ezs ½z2 þ ðl1 f ðzÞ ¼ z3 þ m1 z2 þ m2 z þ m3 b1 e þ a1 þ l2 þ a2 Þz þ l1 ðl2 þ a2 Þ þ a1 l2 ;
ð4:2Þ
with m1 ¼ l1 þ a1 þ l2 þ a2 þ d1 þ b1eI 1 , m2 ¼ ðl2 þ a2 Þ ðl1 þ b1eI 1 Þ þ a1 l2 þ d1 ðl1 þ a1 þ l2 þ a2 þ b1eI 1 Þ; m3 ¼ d1 ½ðl2 þ a2 Þðl1 þ b1eI 1 Þ þ a1 l2 . Thus, the roots of (4.1) consist of the roots of z b2 e S 2 ezs þ l2 þ c2 þ 2 ¼ 0 and the roots of f(z) = 0. Again by Theorem 2.1 of [2], we have that if b2 e S 2 < l2 þ c2 þ 2 (i.e., R2,1 < 1), then z b2 e S 2 ezs þ l2 þ c2 þ 2 ¼ 0 has roots with negative real parts only, if b2 e S 2 > l2 þ c2 þ 2 (i.e., R2,1 > 1), z b2 e S 2 ezs þ l2 þ c2 þ 2 ¼ 0 has root with positive real part. Then E1 is unstable if R2,1 > 1. Next we assume that R2,1 < 1 and consider the roots of f(z) = 0. If s = 0, then f ðzÞ ¼ z3 þ ðl1 þ a1 þ l2 þ a2 þ b1eI 1 Þ z2 þ ½ðl2 þ a2 Þðl1 þ b1eI 1 Þ þ a1 l2 þ d1 b1eI 1 z þ d1 ðl2 þ a2 Þb1eI 1 . By Routh-Hurwitz theorem, all roots of f(z) = 0 have negative real parts if s = 0. We look for purely imaginary roots z = ix of f(z). Substituting z = ix into f(z) and then setting X = x2 give the cubic equation
X 3 þ m21 2m2 X 2 þ ½m22 2m1 m3 ðb e S 1 Þ2 ðl þ a1 þ l þ a2 Þ2 X þ m2 1
1
3
2
S 1 Þ2 ðl1 ðl2 þ a2 Þ þ a1 l2 Þ2 ¼ 0: ðb1 e
ð4:3Þ
We can prove that m21 2m2 ¼ ðl1 þ a1 þ b1eI 1 Þ2 þ ðl2 þ a2 Þ2 þ 2a1 a2 þ d21 > 0; m23 ðb1 e S 1 Þ2 ðl1 ðl2 þ a2 Þ þ a1 l2 Þ2 ¼ ½m3 þ b1 e S 1 ðl1 ðl2 þ a2 Þ þ a1 l2 Þðl2 þ a2 Þb1eI 1 d1 > 0.
S 1 Þ2 ðl1 þ a1 þ l2 þ a2 Þ2 m22 2m1 m3 ðb1 e e 2 þ d2 b eI 1 ð2l þ 2a1 ¼ ðl1 l2 þ l1 a2 þ a1 l2 Þ2 R 1 1 01 1 þ b1eI 1 Þ 2d21 ðl1 l2 þ l1 a2 þ a1 l2 Þ:
ð4:4Þ
We know that
ð4:5Þ
b1eI 1 ð2l1 þ 2a1 þ b1eI 1 Þ 2ðl1 l2 þ l1 a2 þ a1 l2 Þ > 2 e l1 þa1 Þ e 01 > 1þ 1 , thus if R ðl1 l2 þ l1 a2 þ a1 l2 Þ ð R 01 1Þð l þa2
then
2
l2 þa2 l1 þa1 ,
we get
m22
In this case, the characteristic equation is given by
z þ p1 z3 þ p2 z2 þ p3 z þ p4 ezs ðq0 z3 þ q1 z2 þ q2 z þ q3 Þ 4
þ e2zs ðm0 z2 þ m1 z þ m2 Þ ¼ 0;
ð4:6Þ
b1 I1
b2 I2 ,
with p1 ¼ l1 þ a1 þ d1 þ þ l2 þ a2 þ d2 þ p2 ¼ d1 ðl1 þ a1 þ b1 I1 Þ þ d2 l2 þ a2 þ b2 I2 þ l1 þ d1 þ b1 I1 l2 þ a2 þ d2 þ b2 I2 þ a1 l2 þ d2 þ b2 I2 , p3 ¼ d1 l1 þ a1 þ b1 I1 Þ l2 þ d2 þ b2 I2 þ a2 d1 l1 þ b1 I1 þ d2 l2 þ a2 þ b2 I2 Þ l1 þ d1 þ b1 I1 þ a1 d2 l2 þ b2 I2 , p4 ¼ d1 d2 l1 þ a1 þ b1 I1 Þ l2 þ b2 I2 þ a2 l1 þ b1 I1 , q0 ¼ b1 S1 þ b2 S2 , q1 ¼ b1 S1 l1 þ a1 þ l2 þ a2 þ d2 þ b2 I2 þ b2 S2 l1 þ a1 þ l2 þ a2 þ d1 þ b1 I1 Þ, q2 ¼ b1 S1 d2 l2 þ a2 þ b2 I2 þ l1 l2 þ
a2 þ d2 þ b2 I2 Þ þ a1 l2 þ d2 þ b2 I2 þ b2 S2 d1 l1 þ a1 þ b1 I1
þl l þ a1 þ d1 þ b1 I1 þ a2 l1 þ d1 þ b1 I1 , q3 ¼ b1 S1 d2 2 1 l2 þ a2 þ b2 I2 l1 þ l2 þ b2 I2 a1 þ b2 S2 d1 l1 þ a1 þ b1 I1 Þl2 þ l1 þ b1 I1 a2 , m0 ¼ b1 b2 S1 S2 , m1 ¼ b1 b2 S1 S2 ðl1 þ a1 þ l2 þ a2 Þ, m2 ¼ b1 b2 S1 S2 ðl1 l2 þ l1 a2 þ a1 l2 Þ. When s = 0, (4.6) reduces to
ð4:7Þ b2 I2 ;
2
b1eI 1 l2 þ a2
¼ 2b1eI 1 ðl1 þ a1 Þ
1
4.2. Local stability of E⁄
where h ¼ l1 þ a1 þ þ l2 þ a2 þ h2 ¼ d1 b1 I1 þ 1 d2 b2 I2 þ l1 þ b1 I1 l2 þ a2 þ b2 I2 þ a1 l2 þ b2 I2 ; h3 ¼ d1 b1 I1 l2 þ a2 þ b2 I2 þ d2 b2 I2 l1 þ a1 þ b1 I1 ; h4 ¼ d1 d2 b1 I1 b2 I2 . By Routh–Hurwitz theorem, (4.7) has roots with negative real parts only requires that h1 > 0, h3 > 0, h4 > 0
þ a1 l2 þ 2a1 a2
2ðl1 l2 þ l1 a2 þ a1 l2 Þ l2 þ a2 e 01 1Þ; ðl þ a1 Þð R
e 02 > 1 þ l1 þa1 and R1,2 < 1, then E2 is locally Theorem 4.2. If R l2 þa2 asymptotically stable for all s P 0, whereas if R1,2 > 1, then E2 is unstable. Numerical simulations show that Hopf bifurcation giving periodic solutions is possible (see Fig. 4).
b1 I1
> 2ðl1 l2 þ l1 a2 þ a1 l2
¼
Similar results hold for boundary equilibrium E2.
z4 þ h1 z3 þ h2 z2 þ h3 z þ h4 ¼ 0;
b1eI 1 b1eI 1 ð2l1 þ 2a1 þ b1eI 1 Þ ¼ ½ðl l þ l1 a2 l2 þ a2 1 2 e 01 þ l l þ l a2 þ a1 l2 Þ R 1 2 1
þ a1 a2 Þ
e 01 > 1 þ 2 2 and R2,1 < 1, then E1 is Theorem 4.1. If R l1 þa1 locally asymptotically stable for all s P 0, whereas if R2,1 > 1, then E1 is unstable.
2m1 m3 ðb1 e S 1 Þ2 ðl1 þ a1 þ l2 þ a2 Þ2 >
and h3 ðh1 h2 h3 Þ h1 h4 > 0. Then we only need to show the last inequality.Direct computation gives h3 ðh1 h2 2 h3 Þ h1 h4 ¼ h1 h3 l1 þ b1 I1 l2 þ a2 þ b2 I2 þ a1 l2 þ b2 2 I2 Þ þ l1 þ a1 þ b1 I1 l2 þ a2 þ b2 I2 d1 b1 I1 d2 b2 I2 > 0, ⁄ then E is locally asymptotically stable for s = 0. Numerical simulations show that Hopf bifurcation giving periodic solutions is possible (see Fig. 3). 4.3. Permanence The following lemma is useful in the proof of persistence of system (2.1).
946
J. Liu / Chaos, Solitons & Fractals 45 (2012) 942–949
(a)
(b)
e 01 ¼ 17:3217; R e 02 ¼ 2:151; R1;2 ¼ 12:8788; R2;1 ¼ Fig. 1. The disease persists in patch one and dies out in patch two with a2 = 0.01, a1 = 0.005, then R 0:9498. (a) Time series of the infectives I1(t); (b) time series of the infectives I2(t).
(a)
(b)
e 01 ¼ 13:9906; R e 02 ¼ 2:9355; Fig. 2. The disease persists in patch one, the disease oscillates periodically in patch two with a2 = 0.01,a1 = 0.01, then R R1;2 ¼ 8:8542; R2;1 ¼ 1:0234. (a) Time series of the infectives I1(t); (b) time series of the infectives I2(t).
Lemma 4.1 (See [23]). Consider the following equation
_ xðtÞ ¼ axðt sÞ bxðtÞ;
ð4:8Þ
where a, b, s > 0; x(t) > 0 for s 6 t 6 0. We have (i) if a < b, then limt?+1x(t) = 0; (ii) if a > b, then limt?+1x(t) = +1. Using the same method as that in [24], we can give the persistence of system (2.1). e 01 ¼ b1 S01 > 1; R e 02 ¼ b2 S02 > 1, Theorem 4.3. Let R l1 þc1 þ1 l2 þc2 þ2 then there is a positive constant g such that every solution (S1(t), I1(t), R1(t), S2(t), I2(t), R2(t)) of system (2.1) with /2(h) > 0, /5(h) > 0 for h 2 [s, 0] satisfies
lim inf Si ðtÞ P g; t!1
lim inf Ii ðtÞ P g; t!1
permanence, that is, if the conditions in Theorem 4.3 hold, then (2.1) is permanent.
lim inf Ri ðtÞ P g; t!1
i ¼ 1; 2: Remember that Lemma 2.2 tells us system (2.1) is a dissipative system, then uniform persistence is equivalent to
5. Numerical simulations and discussions To investigate (2.1) numerically, we take parameter values as: A1 = 0.25, b1 = 0.3, l1 = 0.006, c1 = 0.7, 1 = 0.05, A2 = 0.15, b2 = 0.08, l2 = 0.006, c2 = 0.8, 2 = 0.05, s = 40. The initial value is (/1(h), /2(h), /3(h), /4(h), /5(h), /6(h)) = (1, 1, 1, 1, 1, 1) for h 2 [s, 0]. We fix the travel rate from patch 2 to patch 1 by assuming a2 = 0.01, we then change the travel rate from patch 1 to patch 2, i.e., a1 is varied from 0.005 to 0.2. In this case, R01 = 16.5344 > 1, R02 = 2.3364 > 1, this means that the disease persists in both patches without travel. As a1 varies, several possible outcomes appear. For a1 = 0.005, the infectives tend to endemic level in patch one (see Fig. 1(a)), while disease dies out in patch two (see Fig. 1(b)). For a1 = 0.01, oscillation appears in patch two (see Fig. 2(a)), disease still persists in patch one (see Fig. 2(b)). Increasing a1 to 0.07, we see oscillations appear in both patches, and
J. Liu / Chaos, Solitons & Fractals 45 (2012) 942–949
thus disease persists (see Fig. 3(a) and (b)). Further increases a1 to 0.2, then disease dies out in patch one (see Fig. 4(a)), and oscillations appears in patch two (see Fig. 4(b)). Then we take different values as: b1 = 0.01, b2 = 0.03, a1 = 0.02, the other parameter values are as before. Then we have R01 = 0.5511 < 1, R02 = 0.8762 < 1, then disease will die out in each patch without travel, Fig. 5(a) indicates that disease dies out in patch one, but the disease is endemic in patch two with travel (see Fig. 5(b)). Bifurcations of the asymptotic state of the infectives as a function of a1 or a2 are plotted in Figs. 6 and 7. The parameter values are the same as those in Fig. 1 except for a1 and a2. Fig. 6(a) and (b) give the bifurcation diagrams of I1(t) and I2(t) as a function of a1. In Fig. 6, we fix a2 = 0.01. We see that as a1 increases, the disease in patch one persists firstly, then oscillates periodically and the amplitude of I1(t) increases as a1 increases, and finally the disease becomes extinct (see Fig. 6(a)). While for patch two, the disease dies out when a1 < 0.005, and after that the disease oscillates and the amplitude of I2(t) becomes bigger and bigger (see Fig. 6(b)).
(a)
947
In Fig. 7, we fix a1 = 0.005. For patch one, as a2 increases, the disease firstly persists, then oscillates periodically, the amplitude of the periodic solution I1(t) does not vary much as a2 increases (see Fig. 7(a)). For patch two, the disease firstly oscillates periodically and the amplitude of the periodic solution I2(t) becomes smaller and smaller until a2 = 0.01, after that then the disease dies out (see Fig. 7(b)). Fig. 6 shows that when we fix the travel rate a2 of susceptible individuals in patch two, then as the travel rate a1 of susceptible individuals in patch one increases, the disease will decreases in patch one and finally die out, while for patch two, the decrease increases monotonically. This is very reasonable, when a1 increases, more susceptibles move to patch two from patch one, then in patch one the number of people to be infected is smaller, hence, the spread of disease will reduce. In contrast to patch one, the disease in patch two intensifies. Fig. 7 tells us the similar situation. In this paper, we discussed a two patches SIR model, and observed that travel between patches can cause substantial changes in the behavior of the system. The analysis to system (2.1) shows that its dynamic behavior is deter-
(b)
e 01 ¼ 4:2297; R e 02 ¼ 5:2344; R1;2 ¼ 1:864; R2;1 ¼ 1:9065. (a) Time series of Fig. 3. Periodic oscillations occur in both patches with a2 = 0.01, a1 = 0.07, then R the infectives I1(t); (b) time series of the infectives I2(t).
(a)
(b)
e 01 ¼ 1:6841; R e 02 ¼ 5:8339; R1;2 ¼ 0:6877; Fig. 4. The disease dies out in patch one, periodic oscillations occur in patch two with a2 = 0.01,a1 = 0.2, then R R2;1 ¼ 3:8201. (a) Time series of the infectives I1(t); (b) time series of the infectives I2(t).
948
J. Liu / Chaos, Solitons & Fractals 45 (2012) 942–949
(a)
(b)
e 01 ¼ 0:3368; R e 02 ¼ 1:4441; R1;2 ¼ Fig. 5. The disease dies out in patch one and persists in patch two with b1 = 0.01, b2 = 0.03, a2 = 0.01, a1 = 0.02, then R 0:2724; R2;1 ¼ 3:6405. (a) Time series of the infectives I1(t); (b) time series of the infectives I2(t).
(a)
(b)
Fig. 6. (a) Bifurcation diagram of I1(t) as a function of a1; (b) bifurcation diagram of I2(t) as a function of a1.
(a)
(b)
Fig. 7. (a) Bifurcation diagram of I1(t) as a function of a2; (b) bifurcation diagram of I2(t) as a function of a2.
mined by several different reproduction numbers: R01, R02 for the two patches in isolation, the travel modified reproe 01 ; R e 02 . The existence and stability of duction number R boundary equilibria and the endemic equilibrium also
depend on the quantities R1,2 and R2,1. We have show that e 01 ; R e 02 g < 1, then the disease-free equilibif R0 ¼ maxf R rium E0 is globally asymptotically stable, otherwise it is unstable. From the theoretical analysis, we know that
J. Liu / Chaos, Solitons & Fractals 45 (2012) 942–949
e 01 > 1; R e 02 > 1; R1;2 < 1; R2;1 < 1, then boundary when R equilibria E1, E2 exists, and E0 is unstable, endemic equilibrium E⁄ does not exists, thus it is possible to have both E1 and E2 locally stable, and thus have bistability. Some sufficient conditions for the local stability of boundary equilibria E1 and E2 are obtained. Numerical simulations indicate that travel between patches can reduce oscillations (Fig. 1); enhance oscillations (Fig. 3); travel switches oscillations (Fig. 4). Our results show that, for the parameter values considered, travel restriction (change the travel rate) does not necessarily always have a positive impact on the overall spread of disease. An increase of the dispersal rates from one patch to the other sometimes may reduce the disease spread in one patch while intensify the disease spread in the other. In our model, we assume that the infectives do not travel between patches, corresponding to either a very severe disease so that infectives are not able to travel or travel is forbidden in order to control a disease outbreak. If we consider the model in which infectives can travel between patches, then the analysis is more difficult. Acknowledgements This research was partially supported by the National Natural Science Foundation of China (No. 11001215 and 11101323), the Shaanxi Technology Department Foundation (No. 2010JQ1016) and the Doctoral Scientific Research Foundation of Xi’an Polytechnic University (No. BS1002). We would like to thank the editor and reviewers for their careful reading and valuable comments which led to an improvement of our original manuscript.
[4] [5]
[6] [7] [8] [9]
[10] [11] [12] [13]
[14] [15]
[16]
[17] [18] [19]
[20] [21]
References [1] Cooke KL. Stability analysis for a vector disease model. Rocky Mountain J Math 1979;9:31–42. [2] Kuang Y. Delay differential equations with applications in population dynamics. Boston: Academic Press; 1993. [3] Castillo-Chavez C, Thieme HR. Asymptotically autonomous epidemic models. In: Arino O, Axelrod D, Kimmel M, Langlais M, (Eds.),
[22] [23] [24]
949
Mathematical population dynamics: analysis of heterogeneity, I. Theory of Epidemics, Wuerz, Winnepeg, Canada, 1995. pp. 33–50. Beretta E, Takeuchi Y. Convergence results in SIR epidemic model with varying population sizes. Nonlinear Anal 1997;28:1909–21. Takeuchi Y, Ma W, Beretta E. Global asymptotic properties of a delay SIR epidemic model with finite incubation times. Nonlinear Anal 2000;42:931–47. Anderson RM, May RM. Population biology of infectious diseases: Part I. Nature 1979;280:361–7. Hethcote HW. Qualitative analyses of communicable disease models. Math Biosei 1976;28:335–56. Beretta E, Takeuchi Y. Global stability of an SIR epidemic model with time delays. J Math Biol 1995;33:250–60. Beretta E, Hara T, Ma W, Takeuchi Y. Global asymptotically stability of an SIR epidemic model with distributed time delay. Nonlinear Anal 2001;47:4107–15. Ma W, Takeuchi Y, Hara T, Beretta E. Permanence of an SIR epidemic model with distributed time delays. Tohoku Math J 2002;54:581–91. Ma W, Song M, Takeuchi Y. Global stability of an SIR epidemic model with time delay. Appl Math Lett 2004;17:1141–5. Wang W, Ruan S. Simulating the SARS outbreak in Beijing with limited data. J Theor Biol 2004;227:369–79. Sattenspiel L, Herring DA. Structrued epidemic models and the spread of influenza in the central Canada subarctic. Human Biol 1998;70:91–115. Sattenspiel L, Dietz K. A structured epidemic model incorporating geographic mobility among regions. Math Biosci 1995;128:71–91. Colizza V, Vespignani A. Epidemic modeling in metapopulation systems with heterogeneous coupling pattern: theory and simulations. J Theor Biol 2008;251:450–67. Colizza V, Barrat A, Barthelemy M, Vespignani A. Epidemic predicatability in meta-population models with heterogeneous couplings: the impact of disease parameter values. Int J Bif Chaos 2007;17:2491–500. Liebovitch LS, Schwartz IB. Migration induced epidemics: dynamics of flux-based multipatch models. Phys Lett A 2004;332:256–67. Keeling MJ, Rohani P. Estimating spatial coupling in epidemiological systems: a mechanistic approach. Ecol Lett 2002;5:20–9. Bianco S, Shaw LB. Asymmetry in the presence of migration stabilizes multistrain disease outbreaks. Bull Math Biol 2011;73:248–60. Wang W, Mulone G. Threshold of disease transmission in a patch environment. J Math Anal Appl 2003;285:321–35. Wang W, Zhao X-Q. An epidemic model in a patchy environment. Math Biosci 2004;190:97–112. Wang W, Zhao X-Q. An age-structured epidemic model in a patchy environment. SIAM J Appl Math 2005;65:1597–614. Xiao Y, Chen L. Modeling and analysis of a predator-prey model with disease in the prey. Math Biosci 2001;171:59–82. Liu J, Zhang T, Teng Z. Spread of disease involving time delay and population dispersal. Dynam Syst 2008;23:267–82.
Contents lists available at SciVerse ScienceDirect
Chaos, Solitons & Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos
Spread of disease in a patchy environment Junli Liu School of Science, Xi’an Polytechnic University, Xi’an 710048, PR China
a r t i c l e
i n f o
Article history: Received 12 February 2011 Accepted 19 March 2012 Available online 28 April 2012
a b s t r a c t For a two patches SIR model, it is shown that its dynamic behavior is determined by several quantities. We have shown that if R0 < 1, then the disease-free equilibrium is globally asymptotically stable, otherwise it is unstable. Some sufficient conditions for the local stability of boundary equilibria are obtained. Numerical simulations indicate that travel between patches can reduces oscillations in both patches; may enhances oscillations in both patches; or travel switches oscillations from one patch to another. Ó 2012 Elsevier Ltd. All rights reserved.
1. Introduction An SIR model was proposed by Cooke [1] for epidemics which spread in a human population via a vector. Susceptible individuals receive the infection from infectious vectors, and susceptible vectors receive the infection from infectious individuals. It is assumed that when a susceptible vector is infected by a person, there is a time s > 0 during which the infectious agents develop in the vector and it is only after that time that the infected vector becomes itself infectious. It is also assumed that the vector population is very large and at any time t the infectious vector population is simply proportional to the infectious human population at time t s. Thus, if we denote S(t), I(t), and R(t) as the numbers of the population susceptible to the disease, of infective members, and of members who have been removed from the possibility of infection through full immunity, respectively. Assume that the recruitment into the susceptible class is A. Let b be the transmission coefficient with mass action incidence, l be the natural death rate, c be the recovery rate of infectives, be the disease-reduced death rate. Then we obtain
8 _ > < SðtÞ ¼ A bSðtÞIðt sÞ lSðtÞ; _IðtÞ ¼ bSðtÞIðt sÞ ðl þ c þ ÞIðtÞ; > :_ RðtÞ ¼ cIðtÞ lRðtÞ:
ð1:1Þ
E-mail address: [email protected] 0960-0779/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.chaos.2012.03.007
where S(t) + I(t) + R(t) = N(t) denotes the number of a population at time t. For the detailed biological meanings, we refer to [1,4–7]. System (1.1) has one disease-free equilibrium lA ; 0; 0 , since lþ1cþ is the average infection period, we can see that e 0 ¼ bA R lðlþcþÞ is a reproduction number of the disease. Stability of (1.1) was studied in [4,5,8–10]. It is shown that, e 0 < 1, the disease-free equilibrium P0 A ; 0; 0 is when R l
globally asymptotically stable for any time delay s. If e 0 > 1; P 0 becomes unstable and there is a positive equiR lS Þ librium P S ; AbSlS ; cðA with S ¼ lþbcþ, and P⁄ is locally lbS asymptotically stable for any time delay s. As for the global asymptotic stability on P⁄, some sufficient conditions were obtained in [5,9]. It is also shown that system (1.1) is pere 0 > 1 [10,11]. manent if R System (1.1) assumes that the mix of population is homogeneous. Because of the dispersal among populations, infection will be carried to many places, such as the SARS outbreak of 2002–2003. Much has been done in terms of modeling spatial spread of diseases. For example, Wang and Ruan [12] studied the global spread pattern of 2002– 03 SARS outbreak. Sattenspiel and Herring [13] used a model to the consideration of travel between populations in the Canadian subarctic. Sattenspiel and Dietz [14] studied a model for the spread of infectious diseases incorporating a mobility process among discrete geographic regions. In [15], Colizza and Vespignani used the basic reaction– diffusion equations to describe the metapopulation system
943
J. Liu / Chaos, Solitons & Fractals 45 (2012) 942–949
at the mechanistic level and derive an early stage dynamics approximation for the subpopulation invasion dynamics. In [16], Colizza et al. considered the problem of the predictability of epidemic scenarios, within the global stochastic modeling framework. By the level of sampling of the underlying travel pattern by infected and latent individuals, the overall predictability can be determined. Bianco and Shaw studied the stability of a multistrain disease model in two coupled populations in [19]. The metapopulation dynamics of infectious diseases has generated a lot of models and results considering population moved among regions with endemic disease [17,18,20–22]. To investigate the spatial effects, in the next section, we will construct a two patches SIR model, and assume the susceptible and recovered individuals will travel between patches. We discuss the existence of equilibria in Section 3, and show that R0 < 1 implies the global stability of the disease-free equilibrium. In Section 4, we consider the local stability of boundary equilibria and endemic equilibrium, and we show the permanence of the model. We conclude with some numerical simulations and discussions in Section 5 about the role of travel in a disease outbreak.
Proof. From system (2.1), the total population N(t) satis_ fies NðtÞ ¼ A1 þ A2 l N 1 ðtÞ l N 2 ðtÞ 1 I1 ðtÞ 2 I2 ðtÞ. 1
2
It follows form I1(t) < N1(t), I2(t) < N2(t) that A1 þ A2 _ 6 A1 þ A2 minfl1 ; l2 g maxfl1 þ 1 ; l2 þ 2 gNðtÞ 6 NðtÞ h n o 2 NðtÞ. Hence, NðtÞ 2 min Nð0Þ; maxflAþ1 þA ; max 1 1 ;l2 þ2 g n oi A1 þA2 Nð0Þ; minfl ;l g for t P 0. This lemma is complete. 1
2
3. Existence of equilibria and stability of the disease-free equilibrium It is clear that system (2.1) always has a disease-free equilibrium E0(S01, 0, 0, S02, 0, 0) for all s P 0, with S01 ¼ A1 l2 þðA1 þA2 Þa2 l1 l2 þl1 a2 þa1 l2 ;
1 þA2 Þa1 S02 ¼ lA2ll1þþðA l a2 þa1 l . As variables R1, R2 do 1
2
1
2
not appear in the first, second, fourth and fifth equations in (2.1). Thus, for local stability, we consider the system in Si and Ii for i = 1, 2. Linearizing these equations about E0 gives the characteristic equation
½ðz þ l1 þ a1 Þðz þ l2 þ a2 Þ a1 a2 gðzÞ ¼ 0;
ð3:1Þ
where 2. Model formulation
gðzÞ ¼ ðz þ l1 þ c1 þ 1 b1 S01 ezs Þðz þ l2 þ c2
In this section, we will construct a model with two patches, for each patch, consider an SIR model as (1.1). We also assume that the disease renders infectives too sick to travel. Let ai and ci denote the travel rates of susceptible individuals and recovered individuals from patch i to patch j, i, j = 1, 2, respectively. The dynamics is governed by the delay differential system
8 > S_ ðtÞ ¼ A1 b1 S1 ðtÞI1 ðt sÞ l1 S1 ðtÞ þ a2 S2 ðtÞ a1 S1 ðtÞ; > > 1 > > _ > > > I1 ðtÞ ¼ b1 S1 ðtÞI1 ðt sÞ ðl1 þ c1 þ 1 ÞI1 ðtÞ; > > < R_ ðtÞ ¼ c I ðtÞ l R ðtÞ þ c R ðtÞ c R ðtÞ; 1 2 2 1 1 1 1 1 1 _ 2 ðtÞ ¼ A2 b S2 ðtÞI2 ðt sÞ l S2 ðtÞ þ a1 S1 ðtÞ a2 S2 ðtÞ; > S > 2 2 > > > > I_ ðtÞ ¼ b S ðtÞI ðt sÞ ðl þ c þ ÞI ðtÞ; > > 2 2 2 2 2 2 2 2 > > :_ R2 ðtÞ ¼ c2 I2 ðtÞ l2 R2 ðtÞ þ c1 R1 ðtÞ c2 R2 ðtÞ: ð2:1Þ
Lemma 2.1. The solution of (2.1) with an initial condition (2.2) is nonnegative for all t P 0. Lemma 2.2. The total population N(t) of system (2.1) is h n o 1 þA2 ; bounded and is in the interval min Nð0Þ; maxflAþ 1 1 ;l2 þ2 g n o A1 þA2 max Nð0Þ; minfl ;l g . 2
Clearly, (z + l1 + a1)(z + l2 + a2) a1a2 = 0 admits zeros with negative real parts only. Therefore, the stability of E0 is determined by the distribution of roots of g(z) = 0. For s = 0, g(z) = 0 reduces to
ðz þ l1 þ c1 þ 1 b1 S01 Þðz þ l2 þ c2 þ 2 b2 S02 Þ ¼ 0: ð3:3Þ Then when s = 0, E0 is locally stable if and only if b1S01 < l1 + c1 + 1 and b2S02 < l2 + c2 + 2. By Theorem 2.1 of [2], we have that if b1S01 < l1 + c1 + 1 and b2S02 < l2 + c2 + 2, then E0 is locally asymptotically stable for all s P 0; if b1S01 > l1 + c1 + 1 or b2S02 > l2 + c2 + 2, then E0 is unstable for all s P 0. e 01 ¼ b1 S01 ; R e 02 ¼ b2 S02 ; R0 ¼ maxf R e 01 ; R e 02 g. Define R l þc þ1 l þc þ2 1
1
2
2
ð2:2Þ
where R6þ0 ¼ fðx1 ; x2 ; x3 ; x4 ; x5 ; x6 Þ 2 R6 : xi P 0; i ¼ 1; 2; 3; 4; 5; 6g. The total population in each patch is Ni(t) = Si(t) + Ii (t) + Ri(t), i = 1, 2, then the overall total population is N(t) = N1(t) + N2(t). It is straightforward to show that.
1
ð3:2Þ
Theorem 3.1. If R0 < 1, then E0 is globally asymptotically stable for all s P 0, whereas if R0 > 1, E0 is unstable for all s P 0.
The initial conditions for system (2.1) are
ð/1 ðhÞ; /2 ðhÞ; /3 ðhÞ; /4 ðhÞ; /5 ðhÞ; /6 ðhÞÞ 2 C ¼ C ½s; 0; R6þ0 ;
þ 2 b2 S02 ezs Þ:
Proof. From the above analysis, it is left to show that E0 is globally attractive if R0 < 1. From system (2.1) and Lemma 2.1, we get
(
S_ 1 ðtÞ 6 A1 l1 S1 ðtÞ þ a2 S2 ðtÞ a1 S1 ðtÞ; S_ 2 ðtÞ 6 A2 l S2 ðtÞ þ a1 S1 ðtÞ a2 S2 ðtÞ:
ð3:4Þ
2
Consider an auxiliary system
u_ 1 ðtÞ ¼ A1 l1 u1 ðtÞ þ a2 u2 ðtÞ a1 u1 ðtÞ; u_ 2 ðtÞ ¼ A2 l2 u2 ðtÞ þ a1 u1 ðtÞ a2 u2 ðtÞ:
ð3:5Þ
It is easy to show that (3.5) has a globally asymptotically stable equilibrium u1 ; u2 ¼ ðS01 ; S02 Þ, that is, limt?1 u1(t) = S01, limt?1u2(t) = S02. By standard comparison theorem, for sufficiently small e > 0, there is a t0 > 0 such that
944
J. Liu / Chaos, Solitons & Fractals 45 (2012) 942–949
S1(t) < S01 + e and S2(t) < S02 + e for all t P t0. For R0 < 1, we can choose e such that
b1 ðS01 þ eÞ < l1 þ c1 þ 1 ;
b2 ðS02 þ eÞ
< l2 þ c2 þ 2 :
e1; e e 2 Þ exists if and (i) Boundary equilibria E1 ðe S 1 ; eI 1 ; R S 2 ; 0; R e 01 > 1; and E2 ðb b1; b b 2 Þ exists if and only if R S 1 ; 0; R S 2 ; bI 2 ; R e 02 > 1. Here, only if R
l þ c1 þ 1 e l1 l2 þ l1 a2 þ a1 l2 e e ; I1 ¼ ð R 01 1Þ; S1 ¼ 1 b1 b1 ðl2 þ a2 Þ ðl2 þ c2 Þc1eI 1 b A2 þ a1 ðl1 þ c1 þ 1 Þ e1 ¼ ;e S ¼ 1 ; R ðl2 þ a2 Þb1 l1 l2 þ l1 c 2 þ c 1 l2 2 c1 c1eI 1 b A1 þ a2 ðl2 þ c2 þ 2 Þ e2 ¼ ;b S ¼ 2 ; R ðl1 þ a1 Þb2 l1 l2 þ l1 c 2 þ c 1 l2 1 e c 2 c2 I 2 l þ c2 þ 2 b1 ¼ ;b S ¼ 2 ; R l1 l2 þ l1 c 2 þ c 1 l2 2 b2 l l þ l a þ a l 1 2 e 1 2 bI 2 ¼ 1 2 ð R 02 1Þ; b2 ðl1 þ a1 Þ e ðl1 þ c1 Þc2 I 2 b2 ¼ : R l1 l2 þ l1 c 2 þ c 1 l2
ð3:6Þ
Then from (2.1), we get
(
I_1 ðtÞ 6 b1 ðS01 þ eÞI1 ðt sÞ ðl1 þ c1 þ 1 ÞI1 ðtÞ; I_2 ðtÞ 6 b ðS02 þ eÞI2 ðt sÞ ðl þ c þ 2 ÞI2 ðtÞ 2
2
ð3:7Þ
2
for t P t0. Consider another auxiliary equation
x_ ¼ b1 ðS01 þ eÞxðt sÞ ðl1 þ c1 þ 1 ÞxðtÞ:
ð3:8Þ
Then x(t) ? 0 as t ? 1 if b1(S01 + e) < l1 + c1 + 1, again by comparison theorem, we have I1(t) ? 0. Similarly, we can show that I2(t) ? 0 when b2(S02 + e) < l2 + c2 + 2. Applying the theory of asymptotically autonomous systems (see [3]) to system (2.1) yields R1(t) ? 0, R2(t) ? 0, S1(t) ? S01 and S2(t) ? S02 as t ? 1. Therefore, if R0 < 1, then E0 is globally asymptotically stable. This completes the proof. If we assume each patch is isolated, then the corresponding basic reproduction number in each patch is
R01 ¼
A1 b 1
l1 ðl1 þ c1 þ 1 Þ
;
R02 ¼
A2 b2
l2 ðl2 þ c2 þ 2 Þ
Corollary 3.2. If R0 < 1, then E0 is the only equilibrium of system (2.1).
l2 þ a2 þ AA2 a1 2 l1 þ a1 þ AA1 2a1 e R : l2 a1 R01 ; R 02 ¼ l2 þ a2 þ l1 l1 þ a1 þ ll1 a2 2 02
Proof. R0 < 1 implies that boundary equilibria E1 and E2 do not exist and
e 01 P R01 and R e 02 6 R02 . The reNote that if AA2 a1 2 P ll2 a1 , then R 1 e 02 P R02 and verse inequality gives AA2 1a2 6 ll2 a1 , then R 1 e 01 , if a1 < A2 a2 R01 e 01 6 R01 . Also from the definition of R R A1
e 01 > 1. Similar result holds for ð1 R01 Þðl2 þ a2 ÞÞ ll1 , then R 2 e R 02 . Therefore, we have the following relationship of R0 to R01 and R02. Corollary 3.1. If R01 > 1 and R02 > 1, then R0 > 1. If R01 < 1 and R02 < 1, then R0 > 1 provided one of the following holds
R01 ð1 R01 Þðl2 þ a2 Þ ll1 . (i) a1 < 2 (ii) a2 < AA1 2a1 R02 ð1 R02 Þðl1 þ a1 Þ ll2 . A2 a2 A1
1
Besides the disease-free equilibrium E0, the other equilibria are stated in the following theorem. First, we define some notations
R1;2 R2;1
l1 þ c1 þ 1
a1 þ l1 ; I1 ¼ ðR1;2 1Þ; b1 b1 ðl þ c2 Þc1 I1 þ c2 c2 I2 l þ c2 þ 2 ; S2 ¼ 2 ; R1 ¼ 2 l1 ðl2 þ c2 Þ þ c1 l2 b2 a2 þ l2 ðl þ c1 Þc2 I2 þ c1 c1 I1 ðR2;1 1Þ; R2 ¼ 1 : I2 ¼ b2 l1 ðl2 þ c2 Þ þ c1 l2
S1 ¼
:
e 01 and R e 02 can be written as Then R
e 01 ¼ R
(ii) System (2.1) has a unique endemic equilibrium E S1 ; I1 ; R1 ; S2 ; I2 ; R2 if and only if R1,2 > 1 and R2,1 > 1. Here,
A1 b1 b1 a2 ðl2 þ c2 þ 2 Þ ¼ þ ; ðl1 þ c1 þ 1 Þða1 þ l1 Þ b2 ðl1 þ c1 þ 1 Þða1 þ l1 Þ A2 b2 b2 a1 ðl1 þ c1 þ 1 Þ ¼ þ : ðl2 þ c2 þ 2 Þða2 þ l2 Þ b1 ðl2 þ c2 þ 2 Þða2 þ l2 Þ
S01 <
l1 þ c1 þ 1 b1
¼ S1 ;
S02 <
l2 þ c2 þ 2 b2
¼ S2 :
ð3:9Þ
From the first and the fourth equations of (2.1), we have
A1 þ A2 ¼ l1 S01 þ l2 S02 :
ð3:10Þ
If the endemic equilibrium E⁄ exists, it follows from adding the equations in (2.1) that
A1 þ A2 ¼ l1 S1 þ l2 S2 þ ðl1 þ 1 ÞI1 þ ðl2 þ 2 ÞI2 þ l1 R1 þ l2 R2 : This gives A1 þ A2 > l1 S1 þ l2 S2 . By (3.9) and (3.10), l1 S1 þ l2 S2 < l1 S01 þ l2 S02 < l1 S1 þ l2 S2 , which is a contradiction. Therefore, the endemic equilibrium E⁄ does not exist. Thus, if R0 < 1, the only equilibrium of (2.1) is E0.
4. Stability analysis and permanence In this section, we will discuss the stability of boundary equilibria and endemic equilibrium. 4.1. Local stability of E1 and E2
Theorem 3.2. The system (2.1) can have other three equilibria, and we have the following results:
Set d1 = l1 + c1 + 1, d2 = l2 + c2 + 2. Theorem 3.2 shows e 01 > 1, then E1 exists. Linearizing the Si, Ii, i = 1, that when R
945
J. Liu / Chaos, Solitons & Fractals 45 (2012) 942–949
2 equations in (2.1) about equilibrium E1 gives the characteristic equation
0. Based on the above discussion, we know that if e 01 > 1 þ l2 þa2 , (4.3) has no positive root, and so f(z) has R l þa1
S 2 ezs þ l2 þ c2 þ 2 Þf ðzÞ ¼ 0; ðz b2 e
no purely imaginary zero, thus E1 is locally asymptotically stable.
1
ð4:1Þ
where
l þa
S 1 ezs ½z2 þ ðl1 f ðzÞ ¼ z3 þ m1 z2 þ m2 z þ m3 b1 e þ a1 þ l2 þ a2 Þz þ l1 ðl2 þ a2 Þ þ a1 l2 ;
ð4:2Þ
with m1 ¼ l1 þ a1 þ l2 þ a2 þ d1 þ b1eI 1 , m2 ¼ ðl2 þ a2 Þ ðl1 þ b1eI 1 Þ þ a1 l2 þ d1 ðl1 þ a1 þ l2 þ a2 þ b1eI 1 Þ; m3 ¼ d1 ½ðl2 þ a2 Þðl1 þ b1eI 1 Þ þ a1 l2 . Thus, the roots of (4.1) consist of the roots of z b2 e S 2 ezs þ l2 þ c2 þ 2 ¼ 0 and the roots of f(z) = 0. Again by Theorem 2.1 of [2], we have that if b2 e S 2 < l2 þ c2 þ 2 (i.e., R2,1 < 1), then z b2 e S 2 ezs þ l2 þ c2 þ 2 ¼ 0 has roots with negative real parts only, if b2 e S 2 > l2 þ c2 þ 2 (i.e., R2,1 > 1), z b2 e S 2 ezs þ l2 þ c2 þ 2 ¼ 0 has root with positive real part. Then E1 is unstable if R2,1 > 1. Next we assume that R2,1 < 1 and consider the roots of f(z) = 0. If s = 0, then f ðzÞ ¼ z3 þ ðl1 þ a1 þ l2 þ a2 þ b1eI 1 Þ z2 þ ½ðl2 þ a2 Þðl1 þ b1eI 1 Þ þ a1 l2 þ d1 b1eI 1 z þ d1 ðl2 þ a2 Þb1eI 1 . By Routh-Hurwitz theorem, all roots of f(z) = 0 have negative real parts if s = 0. We look for purely imaginary roots z = ix of f(z). Substituting z = ix into f(z) and then setting X = x2 give the cubic equation
X 3 þ m21 2m2 X 2 þ ½m22 2m1 m3 ðb e S 1 Þ2 ðl þ a1 þ l þ a2 Þ2 X þ m2 1
1
3
2
S 1 Þ2 ðl1 ðl2 þ a2 Þ þ a1 l2 Þ2 ¼ 0: ðb1 e
ð4:3Þ
We can prove that m21 2m2 ¼ ðl1 þ a1 þ b1eI 1 Þ2 þ ðl2 þ a2 Þ2 þ 2a1 a2 þ d21 > 0; m23 ðb1 e S 1 Þ2 ðl1 ðl2 þ a2 Þ þ a1 l2 Þ2 ¼ ½m3 þ b1 e S 1 ðl1 ðl2 þ a2 Þ þ a1 l2 Þðl2 þ a2 Þb1eI 1 d1 > 0.
S 1 Þ2 ðl1 þ a1 þ l2 þ a2 Þ2 m22 2m1 m3 ðb1 e e 2 þ d2 b eI 1 ð2l þ 2a1 ¼ ðl1 l2 þ l1 a2 þ a1 l2 Þ2 R 1 1 01 1 þ b1eI 1 Þ 2d21 ðl1 l2 þ l1 a2 þ a1 l2 Þ:
ð4:4Þ
We know that
ð4:5Þ
b1eI 1 ð2l1 þ 2a1 þ b1eI 1 Þ 2ðl1 l2 þ l1 a2 þ a1 l2 Þ > 2 e l1 þa1 Þ e 01 > 1þ 1 , thus if R ðl1 l2 þ l1 a2 þ a1 l2 Þ ð R 01 1Þð l þa2
then
2
l2 þa2 l1 þa1 ,
we get
m22
In this case, the characteristic equation is given by
z þ p1 z3 þ p2 z2 þ p3 z þ p4 ezs ðq0 z3 þ q1 z2 þ q2 z þ q3 Þ 4
þ e2zs ðm0 z2 þ m1 z þ m2 Þ ¼ 0;
ð4:6Þ
b1 I1
b2 I2 ,
with p1 ¼ l1 þ a1 þ d1 þ þ l2 þ a2 þ d2 þ p2 ¼ d1 ðl1 þ a1 þ b1 I1 Þ þ d2 l2 þ a2 þ b2 I2 þ l1 þ d1 þ b1 I1 l2 þ a2 þ d2 þ b2 I2 þ a1 l2 þ d2 þ b2 I2 , p3 ¼ d1 l1 þ a1 þ b1 I1 Þ l2 þ d2 þ b2 I2 þ a2 d1 l1 þ b1 I1 þ d2 l2 þ a2 þ b2 I2 Þ l1 þ d1 þ b1 I1 þ a1 d2 l2 þ b2 I2 , p4 ¼ d1 d2 l1 þ a1 þ b1 I1 Þ l2 þ b2 I2 þ a2 l1 þ b1 I1 , q0 ¼ b1 S1 þ b2 S2 , q1 ¼ b1 S1 l1 þ a1 þ l2 þ a2 þ d2 þ b2 I2 þ b2 S2 l1 þ a1 þ l2 þ a2 þ d1 þ b1 I1 Þ, q2 ¼ b1 S1 d2 l2 þ a2 þ b2 I2 þ l1 l2 þ
a2 þ d2 þ b2 I2 Þ þ a1 l2 þ d2 þ b2 I2 þ b2 S2 d1 l1 þ a1 þ b1 I1
þl l þ a1 þ d1 þ b1 I1 þ a2 l1 þ d1 þ b1 I1 , q3 ¼ b1 S1 d2 2 1 l2 þ a2 þ b2 I2 l1 þ l2 þ b2 I2 a1 þ b2 S2 d1 l1 þ a1 þ b1 I1 Þl2 þ l1 þ b1 I1 a2 , m0 ¼ b1 b2 S1 S2 , m1 ¼ b1 b2 S1 S2 ðl1 þ a1 þ l2 þ a2 Þ, m2 ¼ b1 b2 S1 S2 ðl1 l2 þ l1 a2 þ a1 l2 Þ. When s = 0, (4.6) reduces to
ð4:7Þ b2 I2 ;
2
b1eI 1 l2 þ a2
¼ 2b1eI 1 ðl1 þ a1 Þ
1
4.2. Local stability of E⁄
where h ¼ l1 þ a1 þ þ l2 þ a2 þ h2 ¼ d1 b1 I1 þ 1 d2 b2 I2 þ l1 þ b1 I1 l2 þ a2 þ b2 I2 þ a1 l2 þ b2 I2 ; h3 ¼ d1 b1 I1 l2 þ a2 þ b2 I2 þ d2 b2 I2 l1 þ a1 þ b1 I1 ; h4 ¼ d1 d2 b1 I1 b2 I2 . By Routh–Hurwitz theorem, (4.7) has roots with negative real parts only requires that h1 > 0, h3 > 0, h4 > 0
þ a1 l2 þ 2a1 a2
2ðl1 l2 þ l1 a2 þ a1 l2 Þ l2 þ a2 e 01 1Þ; ðl þ a1 Þð R
e 02 > 1 þ l1 þa1 and R1,2 < 1, then E2 is locally Theorem 4.2. If R l2 þa2 asymptotically stable for all s P 0, whereas if R1,2 > 1, then E2 is unstable. Numerical simulations show that Hopf bifurcation giving periodic solutions is possible (see Fig. 4).
b1 I1
> 2ðl1 l2 þ l1 a2 þ a1 l2
¼
Similar results hold for boundary equilibrium E2.
z4 þ h1 z3 þ h2 z2 þ h3 z þ h4 ¼ 0;
b1eI 1 b1eI 1 ð2l1 þ 2a1 þ b1eI 1 Þ ¼ ½ðl l þ l1 a2 l2 þ a2 1 2 e 01 þ l l þ l a2 þ a1 l2 Þ R 1 2 1
þ a1 a2 Þ
e 01 > 1 þ 2 2 and R2,1 < 1, then E1 is Theorem 4.1. If R l1 þa1 locally asymptotically stable for all s P 0, whereas if R2,1 > 1, then E1 is unstable.
2m1 m3 ðb1 e S 1 Þ2 ðl1 þ a1 þ l2 þ a2 Þ2 >
and h3 ðh1 h2 h3 Þ h1 h4 > 0. Then we only need to show the last inequality.Direct computation gives h3 ðh1 h2 2 h3 Þ h1 h4 ¼ h1 h3 l1 þ b1 I1 l2 þ a2 þ b2 I2 þ a1 l2 þ b2 2 I2 Þ þ l1 þ a1 þ b1 I1 l2 þ a2 þ b2 I2 d1 b1 I1 d2 b2 I2 > 0, ⁄ then E is locally asymptotically stable for s = 0. Numerical simulations show that Hopf bifurcation giving periodic solutions is possible (see Fig. 3). 4.3. Permanence The following lemma is useful in the proof of persistence of system (2.1).
946
J. Liu / Chaos, Solitons & Fractals 45 (2012) 942–949
(a)
(b)
e 01 ¼ 17:3217; R e 02 ¼ 2:151; R1;2 ¼ 12:8788; R2;1 ¼ Fig. 1. The disease persists in patch one and dies out in patch two with a2 = 0.01, a1 = 0.005, then R 0:9498. (a) Time series of the infectives I1(t); (b) time series of the infectives I2(t).
(a)
(b)
e 01 ¼ 13:9906; R e 02 ¼ 2:9355; Fig. 2. The disease persists in patch one, the disease oscillates periodically in patch two with a2 = 0.01,a1 = 0.01, then R R1;2 ¼ 8:8542; R2;1 ¼ 1:0234. (a) Time series of the infectives I1(t); (b) time series of the infectives I2(t).
Lemma 4.1 (See [23]). Consider the following equation
_ xðtÞ ¼ axðt sÞ bxðtÞ;
ð4:8Þ
where a, b, s > 0; x(t) > 0 for s 6 t 6 0. We have (i) if a < b, then limt?+1x(t) = 0; (ii) if a > b, then limt?+1x(t) = +1. Using the same method as that in [24], we can give the persistence of system (2.1). e 01 ¼ b1 S01 > 1; R e 02 ¼ b2 S02 > 1, Theorem 4.3. Let R l1 þc1 þ1 l2 þc2 þ2 then there is a positive constant g such that every solution (S1(t), I1(t), R1(t), S2(t), I2(t), R2(t)) of system (2.1) with /2(h) > 0, /5(h) > 0 for h 2 [s, 0] satisfies
lim inf Si ðtÞ P g; t!1
lim inf Ii ðtÞ P g; t!1
permanence, that is, if the conditions in Theorem 4.3 hold, then (2.1) is permanent.
lim inf Ri ðtÞ P g; t!1
i ¼ 1; 2: Remember that Lemma 2.2 tells us system (2.1) is a dissipative system, then uniform persistence is equivalent to
5. Numerical simulations and discussions To investigate (2.1) numerically, we take parameter values as: A1 = 0.25, b1 = 0.3, l1 = 0.006, c1 = 0.7, 1 = 0.05, A2 = 0.15, b2 = 0.08, l2 = 0.006, c2 = 0.8, 2 = 0.05, s = 40. The initial value is (/1(h), /2(h), /3(h), /4(h), /5(h), /6(h)) = (1, 1, 1, 1, 1, 1) for h 2 [s, 0]. We fix the travel rate from patch 2 to patch 1 by assuming a2 = 0.01, we then change the travel rate from patch 1 to patch 2, i.e., a1 is varied from 0.005 to 0.2. In this case, R01 = 16.5344 > 1, R02 = 2.3364 > 1, this means that the disease persists in both patches without travel. As a1 varies, several possible outcomes appear. For a1 = 0.005, the infectives tend to endemic level in patch one (see Fig. 1(a)), while disease dies out in patch two (see Fig. 1(b)). For a1 = 0.01, oscillation appears in patch two (see Fig. 2(a)), disease still persists in patch one (see Fig. 2(b)). Increasing a1 to 0.07, we see oscillations appear in both patches, and
J. Liu / Chaos, Solitons & Fractals 45 (2012) 942–949
thus disease persists (see Fig. 3(a) and (b)). Further increases a1 to 0.2, then disease dies out in patch one (see Fig. 4(a)), and oscillations appears in patch two (see Fig. 4(b)). Then we take different values as: b1 = 0.01, b2 = 0.03, a1 = 0.02, the other parameter values are as before. Then we have R01 = 0.5511 < 1, R02 = 0.8762 < 1, then disease will die out in each patch without travel, Fig. 5(a) indicates that disease dies out in patch one, but the disease is endemic in patch two with travel (see Fig. 5(b)). Bifurcations of the asymptotic state of the infectives as a function of a1 or a2 are plotted in Figs. 6 and 7. The parameter values are the same as those in Fig. 1 except for a1 and a2. Fig. 6(a) and (b) give the bifurcation diagrams of I1(t) and I2(t) as a function of a1. In Fig. 6, we fix a2 = 0.01. We see that as a1 increases, the disease in patch one persists firstly, then oscillates periodically and the amplitude of I1(t) increases as a1 increases, and finally the disease becomes extinct (see Fig. 6(a)). While for patch two, the disease dies out when a1 < 0.005, and after that the disease oscillates and the amplitude of I2(t) becomes bigger and bigger (see Fig. 6(b)).
(a)
947
In Fig. 7, we fix a1 = 0.005. For patch one, as a2 increases, the disease firstly persists, then oscillates periodically, the amplitude of the periodic solution I1(t) does not vary much as a2 increases (see Fig. 7(a)). For patch two, the disease firstly oscillates periodically and the amplitude of the periodic solution I2(t) becomes smaller and smaller until a2 = 0.01, after that then the disease dies out (see Fig. 7(b)). Fig. 6 shows that when we fix the travel rate a2 of susceptible individuals in patch two, then as the travel rate a1 of susceptible individuals in patch one increases, the disease will decreases in patch one and finally die out, while for patch two, the decrease increases monotonically. This is very reasonable, when a1 increases, more susceptibles move to patch two from patch one, then in patch one the number of people to be infected is smaller, hence, the spread of disease will reduce. In contrast to patch one, the disease in patch two intensifies. Fig. 7 tells us the similar situation. In this paper, we discussed a two patches SIR model, and observed that travel between patches can cause substantial changes in the behavior of the system. The analysis to system (2.1) shows that its dynamic behavior is deter-
(b)
e 01 ¼ 4:2297; R e 02 ¼ 5:2344; R1;2 ¼ 1:864; R2;1 ¼ 1:9065. (a) Time series of Fig. 3. Periodic oscillations occur in both patches with a2 = 0.01, a1 = 0.07, then R the infectives I1(t); (b) time series of the infectives I2(t).
(a)
(b)
e 01 ¼ 1:6841; R e 02 ¼ 5:8339; R1;2 ¼ 0:6877; Fig. 4. The disease dies out in patch one, periodic oscillations occur in patch two with a2 = 0.01,a1 = 0.2, then R R2;1 ¼ 3:8201. (a) Time series of the infectives I1(t); (b) time series of the infectives I2(t).
948
J. Liu / Chaos, Solitons & Fractals 45 (2012) 942–949
(a)
(b)
e 01 ¼ 0:3368; R e 02 ¼ 1:4441; R1;2 ¼ Fig. 5. The disease dies out in patch one and persists in patch two with b1 = 0.01, b2 = 0.03, a2 = 0.01, a1 = 0.02, then R 0:2724; R2;1 ¼ 3:6405. (a) Time series of the infectives I1(t); (b) time series of the infectives I2(t).
(a)
(b)
Fig. 6. (a) Bifurcation diagram of I1(t) as a function of a1; (b) bifurcation diagram of I2(t) as a function of a1.
(a)
(b)
Fig. 7. (a) Bifurcation diagram of I1(t) as a function of a2; (b) bifurcation diagram of I2(t) as a function of a2.
mined by several different reproduction numbers: R01, R02 for the two patches in isolation, the travel modified reproe 01 ; R e 02 . The existence and stability of duction number R boundary equilibria and the endemic equilibrium also
depend on the quantities R1,2 and R2,1. We have show that e 01 ; R e 02 g < 1, then the disease-free equilibif R0 ¼ maxf R rium E0 is globally asymptotically stable, otherwise it is unstable. From the theoretical analysis, we know that
J. Liu / Chaos, Solitons & Fractals 45 (2012) 942–949
e 01 > 1; R e 02 > 1; R1;2 < 1; R2;1 < 1, then boundary when R equilibria E1, E2 exists, and E0 is unstable, endemic equilibrium E⁄ does not exists, thus it is possible to have both E1 and E2 locally stable, and thus have bistability. Some sufficient conditions for the local stability of boundary equilibria E1 and E2 are obtained. Numerical simulations indicate that travel between patches can reduce oscillations (Fig. 1); enhance oscillations (Fig. 3); travel switches oscillations (Fig. 4). Our results show that, for the parameter values considered, travel restriction (change the travel rate) does not necessarily always have a positive impact on the overall spread of disease. An increase of the dispersal rates from one patch to the other sometimes may reduce the disease spread in one patch while intensify the disease spread in the other. In our model, we assume that the infectives do not travel between patches, corresponding to either a very severe disease so that infectives are not able to travel or travel is forbidden in order to control a disease outbreak. If we consider the model in which infectives can travel between patches, then the analysis is more difficult. Acknowledgements This research was partially supported by the National Natural Science Foundation of China (No. 11001215 and 11101323), the Shaanxi Technology Department Foundation (No. 2010JQ1016) and the Doctoral Scientific Research Foundation of Xi’an Polytechnic University (No. BS1002). We would like to thank the editor and reviewers for their careful reading and valuable comments which led to an improvement of our original manuscript.
[4] [5]
[6] [7] [8] [9]
[10] [11] [12] [13]
[14] [15]
[16]
[17] [18] [19]
[20] [21]
References [1] Cooke KL. Stability analysis for a vector disease model. Rocky Mountain J Math 1979;9:31–42. [2] Kuang Y. Delay differential equations with applications in population dynamics. Boston: Academic Press; 1993. [3] Castillo-Chavez C, Thieme HR. Asymptotically autonomous epidemic models. In: Arino O, Axelrod D, Kimmel M, Langlais M, (Eds.),
[22] [23] [24]
949
Mathematical population dynamics: analysis of heterogeneity, I. Theory of Epidemics, Wuerz, Winnepeg, Canada, 1995. pp. 33–50. Beretta E, Takeuchi Y. Convergence results in SIR epidemic model with varying population sizes. Nonlinear Anal 1997;28:1909–21. Takeuchi Y, Ma W, Beretta E. Global asymptotic properties of a delay SIR epidemic model with finite incubation times. Nonlinear Anal 2000;42:931–47. Anderson RM, May RM. Population biology of infectious diseases: Part I. Nature 1979;280:361–7. Hethcote HW. Qualitative analyses of communicable disease models. Math Biosei 1976;28:335–56. Beretta E, Takeuchi Y. Global stability of an SIR epidemic model with time delays. J Math Biol 1995;33:250–60. Beretta E, Hara T, Ma W, Takeuchi Y. Global asymptotically stability of an SIR epidemic model with distributed time delay. Nonlinear Anal 2001;47:4107–15. Ma W, Takeuchi Y, Hara T, Beretta E. Permanence of an SIR epidemic model with distributed time delays. Tohoku Math J 2002;54:581–91. Ma W, Song M, Takeuchi Y. Global stability of an SIR epidemic model with time delay. Appl Math Lett 2004;17:1141–5. Wang W, Ruan S. Simulating the SARS outbreak in Beijing with limited data. J Theor Biol 2004;227:369–79. Sattenspiel L, Herring DA. Structrued epidemic models and the spread of influenza in the central Canada subarctic. Human Biol 1998;70:91–115. Sattenspiel L, Dietz K. A structured epidemic model incorporating geographic mobility among regions. Math Biosci 1995;128:71–91. Colizza V, Vespignani A. Epidemic modeling in metapopulation systems with heterogeneous coupling pattern: theory and simulations. J Theor Biol 2008;251:450–67. Colizza V, Barrat A, Barthelemy M, Vespignani A. Epidemic predicatability in meta-population models with heterogeneous couplings: the impact of disease parameter values. Int J Bif Chaos 2007;17:2491–500. Liebovitch LS, Schwartz IB. Migration induced epidemics: dynamics of flux-based multipatch models. Phys Lett A 2004;332:256–67. Keeling MJ, Rohani P. Estimating spatial coupling in epidemiological systems: a mechanistic approach. Ecol Lett 2002;5:20–9. Bianco S, Shaw LB. Asymmetry in the presence of migration stabilizes multistrain disease outbreaks. Bull Math Biol 2011;73:248–60. Wang W, Mulone G. Threshold of disease transmission in a patch environment. J Math Anal Appl 2003;285:321–35. Wang W, Zhao X-Q. An epidemic model in a patchy environment. Math Biosci 2004;190:97–112. Wang W, Zhao X-Q. An age-structured epidemic model in a patchy environment. SIAM J Appl Math 2005;65:1597–614. Xiao Y, Chen L. Modeling and analysis of a predator-prey model with disease in the prey. Math Biosci 2001;171:59–82. Liu J, Zhang T, Teng Z. Spread of disease involving time delay and population dispersal. Dynam Syst 2008;23:267–82.