pdf file - National Center for Theoretical Sciences
Impact of Travel Between Patches for Spatial Spread of Disease Ying-Hen Hsieh Department of Applied Mathematics National Chung Hsing University Taichung, Taiwan P. van den Driessche Department of Mathematics and Statistics University of Victoria Victoria, BC, Canada Lin Wang Department of Mathematics and Statistics University of Victoria Victoria, BC, Canada April 21, 2006
1
Corresponding author: P. van den Driessche Department of Mathematics and Statistics University of Victoria Victoria, BC, Canada Fax: +1-250 721 8962 Email: [email protected]
2
Abstract A multi-patch model is proposed to study the impact of travel on the spatial spread of disease between patches with different level of disease prevalence. (i)
The basic reproduction number for the ith patch in isolation, R0 , is obtained along with the basic reproduction number of the system of patches, R0 . Inequalities describing the relationship between these numbers are also given. For a two-patch model with one high prevalence patch and one low prevalence patch, results pertaining to the dependence of R0 on the travel rates between the two patches are obtained. These results show that, while banning travel from the low prevalence patch to the high prevalence patch always contributes to disease control, banning travel of symptomatic travelers only from the high prevalence patch to the low prevalence patch could adversely affect the containment of the outbreak under certain ranges of parameter values. Moreover, banning all travel from the high prevalence region to the low prevalence region could result in the low prevalence patch becoming disease-free, while the high prevalence patch becomes even more disease-prevalent, with the resulting number of infectives in patch 1 alone exceeding the combined number of infectives in both regions if border control had not been in place. Our results demonstrate that border control, if properly implemented, could be useful to stop the spatial spread of disease. Keywords: multi-patch model, spatial spread, basic reproduction number, travel rate, border control, influenza.
1
Introduction
Circulating in Asia since late 2003 and causing more than one hundred human deaths around the world by April 2006 (WHO, 2006a), the H5N1 avian influenza virus jumped to Africa and the Middle East, spread through Eu3
rope, affecting wild bird populations in 13 additional countries within a single month in early 2006 (WOAH, 2006). Seasonal migration alone probably cannot explain the westward spread in Europe, so imports of poultry and pet birds must also be considered as factors that might lead to the spread of H5N1 through countries and continents (Butler, 2006). Furthermore, its initial appearance on the African continent marks a huge leap in its geographical range, and opens up a whole new front where the vast bird reservoir could potentially spark a pandemic of human-to-human infections. The 2003 severe acute respiratory syndrome (SARS) global epidemic demonstrated the ability of infectious disease spreading due to modern globalization to countries in several continents within a matter of days. In its aftermath and with the potential threat of a flu pandemic, several models describing spatial spread of infectious diseases have been proposed. These include the use of multi-patch compartment models to explore the dynamics of spatial spread of disease(see, for example, Arino and van den Driessche, 2003, 2006; Arino et al. 2005; Hyman and LaForce, 2003; Ruan et. al., 2006; Salmani and van den Driessche, 2006; Sattenspiel and Herring, 2003; Wang and Mulone, 2003; Wang and Zhao, 2004). Moreover, previous modeling work has shown that compartmental models incorporating multiple patches (cities, countries, etc.) describe epidemic spread more accurately than non-spatial models, even at early epidemic phases. Examples of such include foot-and-mouth disease (Chowell et. al., 2006), tuberculosis in possums (Fulford et. al., 2002), in addition to models of infectious diseases of humans. In this work we propose a multi-patch model to study the spread of influenza among the patches. We add a subpopulation of partially immune individuals to account for this important feature of influenza. Although the model is intended to be used for theoretical studies of the spread of human influenza, it can also be used as a model for studying the spread of enzootic diseases such as avian flu (H5N1) among birds. Partial immunity might 4
not be as important for studies of infectious disease involving poultry bird populations where the birds are slaughtered for food after a fixed time, it nonetheless could be an important consideration for studies involving wild bird populations. The paper is organized as follows. We formulate the model in Section 2, and in Section 3, we give results regarding the basic reproduction number of the model. In Section 4, some global analysis is obtained for the model with two patches. Numerical simulations, with parameter values relevant for influenza, to complement our analytical results are given in Section 5. Finally in the last section, from our results, we discuss the impact of travel on the spatial spread of disease.
2
Model Formulation
In this section we formulate a model describing the spread of a disease in a population with n patches taking into consideration travel among patches (Arino and van den Driessche, 2006) and also including partially immune individuals (Hyman and LaForce, 2003). The population in patch i is divided into compartments of susceptible, incubating (infected but not yet showing symptoms), infective (infected with symptoms), recovered and partially immune individuals. Let Si , Ei , Ii , Ri , Pi denote, respectively, the associated population size. Then the total population size in patch i is Ni = Si + Ei + Ii + Ri + Pi for i = 1, 2, . . . , n. A partially immune compartment between the recovered and the susceptible compartments is introduced in an influenza model (Hyman and LaForce, 2003) to account for people who have partial immunity to the current strain of influenza from a previous infection by an earlier strain. We include this partially immune compartment, and assume that individuals in this compartment may become infected again (but with a reduced rate). Also the partial immunity wanes, with these 5
individuals returning to the susceptible compartment. For patch i, let Ai be the recruitment, αi , γi , δi and ηi be the progression rates of the incubating, infective, recovered and partially immune individuals, respectively, di > 0 be the natural death rate and ǫi be the disease-related death rate. We assume that individuals do not change their disease state during travel, and mK ij for K = S, E, I, R, P are the constant travel rates from patch j to patch i for i 6= j of susceptible, incubating, infective, recovered and partially immune individuals, respectively, with mK ii = 0. All parameters are assumed to be positive except that ǫi , δi , ηi can be zero, and the travel rate matrices M K = [mK ij ] for K = S, E, I, R, P are irreducible. For more detailed discussion on travel rates, see, e.g., Arino and van den Driessche (2006). The number of new individuals infected by infectives per unit time in patch i is given by βi (Ni )Si Ii . The term βi (Ni )Si is the product of βi (Ni )Ni , the average number of contacts made by each individual in patch i per unit time, and
Si , Ni
the proportion of susceptibles. It is assumed that infectivity βi (Ni ) is
a continuously differentiable, non-increasing function of Ni with βi (0) < ∞. Note that the above assumptions encompass the widely used mass action and standard incidence disease transmission terms defined for Ni > 0 as well as many other forms of saturating incidence. Reduced infectivity of incubating and partially immune individuals is given by σi βi (Ni ) and νi βi (Ni ), respectively, where σi , νi ∈ [0, 1) and it is assumed that infective individuals can reinfect those who are partially immune. The model flow chart for patch i omitting travel and natural death is given in Figure 1. The above assumptions lead to the following SEIRP model for i = 1, 2, ..., n:
6
dSi dt
= Ai − βi (Ni )Si (Ii + σi Ei ) − di Si + ηi Pi +
n X
mSij Sj
−
j=1
dEi dt
= βi (Ni )[Si (Ii + σi Ei ) + νi Pi Ii ] − (di + αi )Ei +
n X j=1
n X
mE ij Ej
−
dIi dt
= αi Ei − (γi + εi + di )Ii +
mIij Ij −
dRi dt
= γi Ii − (di + δi )Ri +
mR ij Rj −
mIji Ii
n X
mR ji Ri
j=1
j=1
dPi dt
mE ji Ei
j=1
j=1
n X
n X
n X j=1
j=1
n X
mSji Si
= δi Ri − (di + ηi )Pi − νi βi (Ni )Pi Ii +
n X j=1
mPij Pj
−
n X
mPji Pi
j=1
(2.1)
Figure 1: Flow diagram of the model. Initially each variable is assumed to be non-negative with Si (0) > 0 and Pn i=1 Ei (0) + Ii (0) > 0. It follows that for a given set of non-negative initial
conditions, there is a unique solution to system (2.1). The total population Pn ¯ size in all n patches is N (t) = i=1 Ni (t). Let d = min{d1 , d2 , . . . , dn } Pn and A = i=1 Ai . Then the following result, which can be proved in a similar way to that of Theorem 1.1 in Salmani and van den Driessche (2006),
indicates that the model is well posed and all variables remain non-negative and bounded. 7
Theorem 2.1. Consider system (2.1) with non-negative initial conditions. Then for each i = 1, 2, . . . , n, Ei (t), Ii (t), Ri (t), Pi (t) remain non-negative, Si (t) and Ni (t) remain positive and the total population N (t) is in the interval ¯ N (0)}]. (0, max{A/d, Notice that our patch model is an extension of those in Arino and van den Driessche (2006), Salmani and van den Driessche (2006). It includes a new compartment, namely, partially immune individuals in each patch, and also the probability of disease transmission by incubating individuals. These additional features may be important for modeling diseases such as influenza and SARS. In Hyman and LaForce (2003), which does not include an incubating compartment in a model for the spread of influenza, travel is assumed to be independent of disease status and symmetric. Moreover, the disease does not cause death, and the population of each patch remains constant. Notice also that our patch model keeps track of individuals present in patch i at time t, but does not keep track of where an individual resides. Models that include this are developed in Arino and van den Driessche (2003, 2006) and Sattenspiel and Herring (2003), and Ruan et. al. (2006) use such a model to study the effect of travel on the spread of SARS.
3
The basic reproduction number
In this section, we derive a formula to compute the basic reproduction number R0 for the general model (2.1) and then give a lower bound for R0 . Both lower and upper bounds of R0 are given for a special case. A disease free equilibrium (DFE) is a steady state solution of system (2.1) with Si = Ni∗ > 0 and all other variables Ei , Ii , Ri , Pi equal to 0 for i = 1, 2, . . . , n. Let S ∗ = (N1∗ , N2∗ , . . . , Nn∗ )T . Then from the first equation in (2.1), there is a DFE if and only if S ∗ satisfies CS ∗ = A with A = P (A1 , A2 , . . . , An )T and C = D − M S , where D = diag( nj=1 mSji + di ). Note 8
that C is irreducible, has positive column sums d1 , d2 , . . . , dn and negative off-diagonal entries. Thus C is a non-singular M -matrix (page 141, Berman and Plemmons, 1979), and therefore C −1 > 0. Hence S ∗ = C −1 A > 0 is the unique solution of CS ∗ = A. This shows that the DFE always exists and is unique. In the absence of disease system (2.1) reduces to just the first equation and S ∗ is stable. Next we consider the local stability of the DFE for system (2.1). To this end, we order the infected variables by E1 , E2 , . . . , En , I1 , I2 , . . . , In and make use of the result in van den Driessche and Watmough (2002) to obtain # " # " diag(σi βi (Ni∗ )Ni∗ ) diag(βi (Ni∗ )Ni∗ ) F11 F12 = F = 0 0 0 0 with V =
"
V11
0
−V21 V22
#
,
where V11 = diag(di + αi +
n X
E mE ji ) − M ,
j=1
V22 = diag(γi + ǫi + di +
n X
mIji ) − M I ,
j=1
V21 = diag(αi ), i = 1, 2, . . . , n. Note that V11 and V22 are both irreducible non-singular M −matrices with positive column sums and hence V11−1 > 0, V22−1 > 0.
(3.1)
The basic reproduction number for the system, denoted by R0 is then the spectral radius of F V −1 , i.e., R0 = ρ{F V −1 } where " # #" −1 F F V 0 11 12 11 F V −1 = . −1 0 0 V22 V21 V11−1 V22−1 9
Therefore R0 = ρ(F11 V11−1 + F12 V22−1 V21 V11−1 ).
(3.2)
The first term accounts for infection from incubating individuals, while the second term accounts for infection from infective individuals who survive the incubating compartment. Travel rates influence the average time spent in the incubating and infective compartments. Note that R0 does not depend on νi , δi , ηi for i = 1, 2, . . . , n. The basic reproduction number gives an important threshold for the disease, as shown in the following result. Theorem 3.1. Consider model (2.1). If R0 < 1, then the DFE is locally asymptotically stable and if R0 > 1, the DFE is unstable. Moreover, if the disease transmission is standard incidence, then the DFE is globally asymptotically stable provided that R0 < 1. Proof. It follows from Theorem 2 of van den Driessche and Watmough (2002) that the DFE is locally asymptotically stable if R0 < 1 and is unstable if R0 > 1. If the disease transmission is standard incidence, then βi (Ni )Ni = βi for i = 1, 2, . . . , n. Note that Si ≤ Ni and Si + νi Pi ≤ Ni . This gives the inequality n
n
X X dEi mE mE E − ≤ βi Ii + σi βi Ei − (di + αi )Ei + ji Ei . ij j dt j=1 j=1 Consider the linear system dEi dt
= βi Ii + σi βi Ei − (di + αi )Ei +
n X
mE ij Ej −
dIi dt
= αi Ei − (γi + εi + di )Ii +
j=1
mIij Ij
mE ji Ei
j=1
j=1
n X
n X
−
n X
(3.3)
mIji Ii .
j=1
The right hand side of the above system has F − V as its coefficient matrix. Again, by proof of Theorem 2 of van den Driessche and Watmough (2002), 10
each eigenvalue of F − V has negative real part if R0 < 1. Thus any solution of (3.3) satisfies lim Ei = 0 and lim Ii = 0 for i = 1, 2, . . . , n. Using a t→∞
t→∞
comparison theorem (Thorem B.1, Smith and Waltman, 1995), each solution of system (2.1) satisfies lim Ei = 0 and lim Ii = 0 for i = 1, 2, . . . , n. Using t→∞
t→∞
a similar argument as in the proof of Theorem 2.2, Salmani and van den Driessche (2006), lim Ri (t) = 0, and similarly lim Pi (t) = 0. Thus the DFE t→∞
t→∞
is globally asymptotically stable provided that R0 < 1. Let ai = γi + ǫi + di , bi = βi (Ni∗ )Ni∗ , ci = di + αi for i = 1, 2, . . . , n. In the case that there is no travel between patch i and all other patches, the basic reproduction number in patch i in isolation is given by (i)
R0 =
σi bi bi αi + . ci ai c i
(3.4)
Define a modified reproduction number (modified by travel) in patch i, e (i) R 0 =
ci +
σb Pini
E j=1 mji
+
(ai +
bi αi I j=1 mji )(ci
Pn
+
Pn
j=1
mE ji )
.
(3.5)
The following result gives bounds for the basic reproduction number for system (2.1) in terms of the numbers defined in (3.4) and (3.5) for each patch. Theorem 3.2. For system (2.1), e (i) R0 ≥ max R 0 .
(3.6)
1≤i≤n
Furthermore, if ai = a, αi = α, σi = σ and di = d for i = 1, 2, . . . , n, then (i) (i) e0 , min R0 ≤ R0 ≤ max R(i) (3.7) max max R 0 . 1≤i≤n
1≤i≤n
1≤i≤n
Proof For j = 1, 2, let Vjj [1′ ] denote the matrix Vjj with row and column
1 deleted, Y = [yij ] and Z = [zij ] denote V11−1 and V22−1 , respectively. Let W = [wij ] = G + H, where G = [gij ] = F11 Y , H = [hij ] = F12 ZV21 Y. It follows from (3.1) that yij > 0, zij > 0 for i, j = 1, 2, . . . , n. Then, by Corollary 8.1.20, Horn and Johnson (1985), R0 = ρ(W ) ≥ wii = gii + hii for i = 1, 2, . . . , n. 11
′
11 [1 ] Note that g11 = σ1 b1 y11 = σ1 b1 detV . By virtue of Fischer’s inequality detV11
(page 117, Horn and Johnson, 1991), it follows that ! n X ′ mE detV11 ≤ c1 + j1 detV11 [1 ]. j=1
Therefore, g11 ≥ Similarly,
c1 +
σb P1n1
h11 = b1 α1 z11 y11 +
j=1
n X
mE j1
.
b1 αk z1k yk1 ,
k=2
thus,
detV22 [1′ ] detV11 [1′ ] . detV22 detV11 Again, by Fischer’s inequality, it follows that h11 ≥ b1 α1 z11 y11 = b1 α1
b1 α1 . h11 ≥ P Pn a1 + j=1 mIj1 c1 + nj=1 mE j1
Adding these shows that
(1)
e0 . R0 ≥ R
Similarly, it can be shown that
(i)
and this gives (3.6).
e0 R0 ≥ R
for i = 2, 3, . . . , n
If ai = a, αi = α, σi = σ and di = d (thus ci = c) for i = 1, 2, . . . , n, P then wij = σbi yij + αbi nk=1 zik ykj for i, j = 1, 2, . . . , n. Without loss of
generality, assume that 0 < b1 ≤ b2 ≤ · · · ≤ bn . From the fact that the
matrix V11 has each column sum equal to c > 0 and the matrix V22 has P Pn each column sum equal to a > 0, it follows that ni=1 yij = 1c , i=1 zij = 12
1 a
for j = 1, 2, . . . , n. Therefore, for the matrix W , the sum of column j is
given by n X
wij =
n X
σbi yij +
≤ σbn = σbn
αbi
i=1
i=1
i=1
n X
n X
i=1 n X
yij + αbn yij + αbn
i=1
n X
zik ykj
i=1 k=1 n n X X k=1
σbn αbn + c ac (n) = R0 .
zik ykj
k=1 n n XX
zik
i=1
!
ykj
=
Similarly,
Pn
i=1
wij ≥
σb1 c
+
αb1 ac
(1)
= R0 . From Theorem 8.1.22, Horn and
Johnson (1985), ρ(W ) lies between the minimum and maximum column sums of W , thus (i)
(i)
min R0 ≤ R0 = ρ(W ) ≤ max R0 .
1≤i≤n
1≤i≤n
Combining with (3.6) immediately gives the desired (3.7) The analytical results give the basic reproduction numbers in isolation of the individual patches as upper and lower bounds for the basic reproduction number of the system. However, its usefulness is rather limited in practice since the range for the bounds could be too large. In the next section, we consider the simple case of only two patches, and obtain more explicit results giving insight about the impact of travel on the spatial spread of disease between patches.
13
4
The Model with 2 Patches
From now on we consider the special case of (2.1) with only 2 patches, i.e., we consider the system dS1 dt dE1 dt dI1 dt dR1 dt dP1 dt dS2 dt dE2 dt dI2 dt dR2 dt dP2 dt
= A1 − β1 (N1 )S1 (I1 + σ1 E1 ) − d1 S1 + η1 P1 + mS12 S2 − mS21 S1 E = β1 (N1 )[S1 (I1 + σ1 E1 ) + ν1 P1 I1 ] − (d1 + α1 )E1 + mE 12 E2 − m21 E1
= α1 E1 − (γ1 + ε1 + d1 )I1 + mI12 I2 − mI21 I1 R = γ1 I1 − (d1 + δ1 )R1 + mR 12 R2 − m21 R1
= δ1 R1 − (d1 + η1 )P1 − ν1 β1 (N1 )P1 I1 + mP12 P2 − mP21 P1 = A2 − β2 (N2 )S2 (I2 + σ2 E2 ) − d2 S2 + η2 P2 + mS21 S1 − mS12 S2 E = β2 (N2 )[S2 (I2 + σ2 E2 ) + ν2 P2 I2 ] − (d2 + α2 )E2 + mE 21 E1 − m12 E2
= α2 E2 − (γ2 + ε2 + d2 )I2 + mI21 I1 − mI12 I2 R = γ2 I2 − (d2 + δ2 )R2 + mR 21 R1 − m12 R2
= δ2 R2 − (d2 + η2 )P2 − ν2 β2 (N2 )P2 I2 + mP21 P1 − mP12 P2 (4.1)
The DFE of (4.1) is given by S∗ =
"
N1∗ N2∗
#
=
"
d1 + mS21
−mS12
−mS21
d2 + mS12
#−1 "
A1 A2
#
.
Also F11 = diag(σ1 b1 , σ2 b2 ), F12 = diag(b1 , b2 ), V21 = diag(α1 , α2 ), and V11 =
"
c1 + m E 21
−mE 12
−mE 21
c2 + m E 12
#
, V22 =
"
a1 + mI21
−mI12
−mI21
a2 + mI12
#
.
Hence R0 = ρ(W ), where W = F11 V11−1 + F12 V22−1 V21 V11−1 = 14
"
w11 w12 w21 w22
#
(4.2)
with w11 =
I I E σ1 b1 (c2 + mE α1 b1 (c2 + mE 12 ) 12 )(a2 + m12 ) + α2 b1 m12 m21 + E E I I c1 c2 + c1 mE (c1 c2 + c1 mE 12 + c2 m21 12 + c2 m21 )(a1 a2 + a1 m12 + a2 m21 )
w12 =
I I E σ1 b1 mE α1 b1 mE 12 12 (a2 + m12 ) + α2 b1 m12 (c1 + m21 ) + E E I I c1 c2 + c1 mE (c1 c2 + c1 mE 12 + c2 m21 12 + c2 m21 )(a1 a2 + a1 m12 + a2 m21 )
w21 =
E I σ2 b2 mE α1 b2 mI21 (c2 + mE 21 12 ) + α2 b2 m21 (a1 + m21 ) + E E I I c1 c2 + c1 mE (c1 c2 + c1 mE 12 + c2 m21 12 + c2 m21 )(a1 a2 + a1 m12 + a2 m21 )
w22
I E I α1 b2 mE σ2 b2 (c1 + mE 12 m21 + α2 b2 (c1 + m21 )(a1 + m21 ) 21 ) . + = I I E E c1 c2 + c1 mE (c1 c2 + c1 mE 12 + c2 m21 12 + c2 m21 )(a1 a2 + a1 m12 + a2 m21 )
It follows from (3.2) that R0 =
p 1 w11 + w22 + (w11 − w22 )2 + 4w12 w21 . 2
(4.3)
From Theorem 3.1, we have the following result.
Theorem 4.1. For the two-patch model, assume that the disease transmission is standard incidence. Then the DFE is globally asymptotically stable if R0 < 1 and is unstable if R0 > 1, where R0 is given by (4.3). As expected, an increase in σi or bi increases R0 . To investigate how R0 changes with the other parameters, we first assume that all travel rates for I incubating and infective individuals are equal, namely, mE ij = mij = m for
i, j = 1, 2. The proof of the following result is given in the Appendix. I Theorem 4.2. If mE ij = mij = m for i, j = 1, 2 and σ1 = σ2 = σ, α1 = α2 =
α, a1 = a2 = a, c1 = c2 = c, b1 6= b2 , then an increase in a or m decreases R0 . Remark 4.1. With the conditions of the above theorem, as the travel rate becomes large, the basic reproduction number R0 in (4.3) approaches the (1)
(2)
(1)
(2)
mean value of R0 and R0 , i.e. R0 → 12 (R0 + R0 ) as m → ∞. This can be seen from the form of W in the proof of Theorem 4.2. In this limit case, the two patches merge into one. 15
Using a similar technique as in proof of Theorem 4.2, we can prove the following result in the case that infectives of both patches are too sick to travel. I Theorem 4.3. If mE ij = m, mij = 0 for i, j = 1, 2 and σ1 = σ2 , α1 = α2 , a1 =
a2 , c1 = c2 , b1 6= b2 , then
∂R0 ∂m
< 0, i.e., the basic reproduction number R0
decreases as m increases.
5
Numerical Simulations
Assuming standard incidence disease transmission, we present some numerical simulations for 2 patches to illustrate how R0 changes with travel rates, with the choice of parameters relevant for human influenza (Ferguson et. al., 2005). In the next section, we present simulations for restricted travel. The model parameters with time unit as one day are taken as: α1 = α2 = 1 1.48
≈ 0.68, γ1 = γ2 =
1 2.6
≈ 0.38 (the average incubating time is 1.48 days
and the average infective time is 2.6 days), ǫ1 = ǫ2 = 4.0 × 10−3 , d1 = d2 =
1 70×365
≈ 3.91 × 10−5 , b1 = 0.6, b2 = 0.1, σ1 = σ2 = 0.04. These
parameters yields the respective basic reproduction numbers in isolation of (1)
(2)
R0 ≈ 1.58 > 1 and R0 ≈ 0.26 < 1. Hence we consider the hypothetical scenario of disease spread between a high-prevalence endemic region (patch 1) and a low-prevalence region where a minor outbreak could be eradicated I (patch 2). We first keep mE 21 = m21 fixed with values at 0.6, 0.8 and 1. We I let mE 12 = m12 vary from 0 to 1. The curves of R0 are given in Figure 2. A I E I 3-d plot of R0 versus mE 12 = m12 and m21 = m21 is given in Figure 3.
Figure 2 shows that by keeping the same travel rate from patch 1 to patch 2, an increase in the travel rate from patch 2 (the patch with lower disease transmission rate) to patch 1 (the patch with higher disease transmission rate) results in an increase in R0 . Consequently, we can conclude that travel between the two patches may cause the disease to become endemic in both 16
1.2 1.1
E
m21=0.6
1
mE =0.8 21
0.9
R
0
0.8
E
m21=1.0
0.7 0.6 0.5 0.4 0.3 0.2 0
0.2
0.4
0.6 E 12
0.8
1
I 12
m =m
I E I Figure 2: R0 vs mE 12 = m12 for fixed m21 = m21 . The three curves from top I to bottom correspond to : mE 21 = m21 = 0.6, 0.8, 1.0
2
R
0
1.5 1 0.5 0 1 1 0.5 E 12
I 12
m =m
0.5 0
0
E 21
I 21
m =m
I E I Figure 3: R0 vs mE 12 = m12 and m21 = m21 .
17
patches if the travel rate from patch 2 to patch 1 is large enough. The numerical simulation in Figure 3 shows that, for the chosen parameters and I E I keeping mE 12 = m12 fixed, an increase in m21 = m21 leads to a decrease in
R0 . In particular, R0 could conceivably decrease to a value less than one if the travel rate from the high-prevalence patch to the low-prevalence patch, I namely mE 21 = m21 , is high enough.
If all travel rates from one patch to the other are the same, i.e. mE 12 = I mI12 = mE 21 = m21 = m, then an increase in m leads to a decrease in R0
(Figure 4), as predicted by Theorem 4.2. Note that, for these parameter values, R0 < 1 if m > 0.4. 1.6 1.5 1.4
R0
1.3 1.2 1.1 1 0.9 0
0.2
0.4
0.6
0.8
1
m I E I Figure 4: R0 vs m, where m = mE 12 = m12 = m21 = m21 .
Taking m = 0.5 and subsequently R0 < 1, the disease will die out eventually. To investigate how the infective population sizes in the two patches change over time, we numerically simulate the system (4.1) with i = 1, 2, and plot I1 and I2 against time in Figure 5. In addition to the parameter values given at the beginning of this section, we take A1 = A2 = 100, η1 = η2 = δ1 = δ2 = 0.2, ν1 = ν2 = 0.1. The infective population sizes in both 18
patches increase initially then slowly decrease to zero since the DFE is globally asymptotically stable if R0 < 1. 70
60
I1 and I2
50
40
30
I
1
I2
20
10
0
0
20
40
60
80
100
t
Figure 5: Numerical solution of system (4.1) with parameters as in the text and m = 0.5 showing I1 and I2 vs time t. Initial conditions are: S1 (0) = 400, E1 (0) = 50, I1 (0) = 20, R1 (0) = P1 (0) = 0; S2 (0) = 800, E2 (0) = 30, I2 (0) = 10, R2 (0) = P2 (0) = 0. Now we use the same parameter values except that m = 0.2, and subsequently R0 > 1. Numerical simulation of this case given in Figure 6 shows that in both patches there is an initial increase in the number of infectives before the infective populations decrease to their endemic levels. However, the existence and stability of such an endemic equilibrium remain unproved analytically. Some interesting observations can be made from our results regarding the role that travel plays in the spatial spread of a disease. Figures 2 and 3 demonstrate the possibility that, for a lower-prevalence patch with a minor (2)
disease outbreak (basic reproduction numbers in isolation R0 less than one) open travel with a high-prevalence patch could lead to the disease becoming 19
45 40
I1 and I2
35
I
1
30 25
I
2
20 15 10 0
20
40
60
80
100
t
Figure 6: Numerical solution of system (4.1) showing I1 and I2 vs time t. Parameter values and initial conditions are the same as in Figure 5 except that m = 0.2. endemic. On the other hand, for a high-prevalence patch with endemic dis(1)
ease in isolation (basic reproduction numbers in isolation R0 greater than one), open travel with a high-prevalence patch could eradicate the disease. Essentially, under appropriate parameter ranges, travel between patches dilutes the overall prevalence to the point that it could either lessen the severity of an endemic patch or worsen a minor outbreak region. Further evidence can be found in Figures 4 and 5 in which assuming all travel rates are equal, the disease can be eventually eradicated if the travel rates are sufficiently large.
6
Discussion of Travel Restrictions
To consider hypothetical intervention scenarios, we let mI21 = 0 and keep all other parameters the same as in Figure 2. This models a situation in which the authority bans all travel of symptomatic travelers from patch 1, the high 20
prevalence region, to patch 2, the low prevalence region. In comparison with Figure 2, Figure 7 shows this gives an overall increase in the value of the basic reproduction number R0 . More significantly, there is a range of parameter values for mI12 = mE 12 in which stopping all travel of symptomatic travelers from patch 1 to patch 2 could adversely impact the epidemic by driving the basic reproduction number above one, thus prolonging the epidemic that otherwise would be eradicated. 1.5 mE =0.6 21
1.4 1.3 1.2
mE =0.8 21
mE =1.0 21
R
0
1.1 1 0.9 0.8 0.7 0
0.2
0.4
E 12
I 0.6 12
m =m
0.8
1
I I Figure 7: R0 vs mE 12 = m12 for fixed m21 = 0. The three curves from top to
bottom correspond to : mE 21 = 0.6, 0.8, 1.0 Conversely, Figure 8 shows that stopping all travel of symptomatic travelers from patch 2 to patch 1 (letting mI12 = 0) would alleviate the epidemic, lowering the basic reproduction number to below one, thus successfully preventing further spread of the outbreak. Therefore, the policy of border control to ban travel of symptomatic travelers only from the high prevalence patch to the low prevalence patch could affect the containment of the outbreak adversely (Figures 2 and 7). However, banning travel of symptomatic travelers only from the low prevalence 21
0.7 0.65 0.6
R0
0.55 0.5 0.45 0.4 0.35
0
0.2
0.4
mE
0.6
0.8
1
12
I E I Figure 8: R0 vs mE 12 with m12 = 0 for fixed m21 = m21 = 0.6
. patch to the high prevalence patch always has a positive impact (Figures 2 and 8). We note here that by banning all travels (both exposed but asymptomatic and symptomatic individuals) would result in an outbreak in the high prevalence patch. We now consider the case in which the disease is endemic, as in Figure 6. We again let mI21 = 0 and keep all other parameters the same as in Figure 6. Figure 9 shows that the disease is still endemic in both patches, but with significantly larger numbers of infectives in both patches. Thus once again, banning all travel of symptomatic travelers from the high prevalence region to the low prevalence region is detrimental to the intervention and control of the outbreak. If we let mI21 = mE 21 = 0, i.e., banning all travel from the high prevalence region to the low prevalence region (Figure 10), patch 2 becomes disease-free, while patch 1 is even more prevalent with the number of infectives in patch 1 alone exceeding the total number of infectives in the two regions if border control had not been in place. Note that in this case the infected equations for patch 2 uncouple. Matrix W given by (4.2) has 22
100
80 I
I1 and I2
1
60
40 I
2
20
0 0
20
40
60
80
100
t
Figure 9: Numerical solution of system (4.1) with two patches showing I1 and I2 vs time t. Parameters as in Figure 6 except that mI21 = 0. w21 = 0 and so is reducible. Thus the reproduction number in patch 1 is (1) e (2) R0 > 1 and the reproduction number in patch 2 is R 0 < 1.
Now we let mI12 = 0 and keep all other parameters the same as in Figure 6,
i.e., banning travel of symptomatic travelers from the low prevalence patch
to the high prevalence patch only, the disease is eradicated in both regions (Figure 11). This indicates the importance of border control out of a low prevalence patch and into a high prevalence patch. We could also consider a slightly different scenario in which both patches have the same infectivity, but initially one patch is disease free while the other patch has infectives, i.e., E2 (0) = I2 (0) = 0 and E1 (0) + I1 (0) > 0. However, since the analytical results do not depend on the initial conditions, the resulting simulations are asymptotically similar to Figures 9-11. During the 2003 SARS outbreak, travel warnings to all affected areas were issued by WHO to prevent travelers from entering and becoming infected. Dozens of countries also issued border control either banning travelers from
23
140 120
I1 and I2
100
I
1
80 60 40 20 0 0
I
2
20
40
60
80
100
t
Figure 10: Numerical solution of system (4.1) with two patches showing I1 and I2 vs time t. Parameters are the same as in Figure 6 except that mI21 = mE 21 = 0.
35 30 I
I1 and I2
25
1
20 I
15
2
10 5 0 0
20
40
60
80
100
t
Figure 11: Numerical solution of system (4.1) with two patches showing I1 and I2 vs time t. Parameters as in Figure 6 except that mI12 = 0.
24
entering or placing them under quarantine. From a simple public health point of view, it is imperative for the low-prevalence region to stop travel of sick (both exposed and asymptomatic) individuals from the high-prevalence region for intervention purposes. See (Ruan et. al., 2006) for a study of a patch model for the spread of SARS. Our results show that border control does not necessarily always have a positive impact on the overall spread of disease and it is more important to ban travel from a low prevalence patch to a high prevalence patch. Moreover, an increase of travel rates in the opposite direction (from a high prevalence patch to a low prevalence patch), while theoretically alleviating the spatial spread of the disease, is not likely an implementable policy. Furthermore, letting mI12 = mE 12 = 0, i.e., banning all travel from the low prevalence patch to the high prevalence patch, numerical simulations in Figure 12 show R0 decreasing drastically as mI21 = mE 21 increase from 0 to 1. Note that in this case the infected equations in patch 1 uncouple, and the matrix W in (4.2) is reducible, thus (2) e (1) R0 = max{R 0 , R0 }. On the other hand, banning all travel from the high
prevalence patch to the low prevalence patch (mI21 = mE 21 = 0) results in (1)
R0 = R0 ≈ 1.58 > 1 and the outbreak becomes endemic. On the basis of our model assumptions and the parameters considered, we could therefore conclude that, retrospectively, WHO acted properly with travel warnings for travelers to avoid all but essential travel to affected areas (WHO, 2006b). Moreover, screening at border for travelers in and out is important mainly for the purpose of quick identification of sick individuals. As a final remark, we note again that, while the model is proposed with the spread of human influenza in mind, it also can be used, with some appropriate modifications and parameter changes, as the basis of a theoretical model to study the spread of enzootic diseases such as avian flu (H5N1) among birds, both that of importation of poultry and pet birds, as well as wild migratory bird populations. 25
1.6 1.4 1.2
R
0
1 0.8 0.6 0.4 0.2 0
0.2
0.4
0.6 I
0.8
1
E
m21=m21
E I Figure 12: R0 vs mI21 = mE 21 for m12 = m12 = 0.
Appendix: Proof of Theorem 4.2 Proof. It is clear that an increase in a decreases each entry in the matrix W , and so decreases the basic reproduction number R0 . Under the parameter assumptions W = where u=
"
b1 u b 1 v b2 v b 2 u
#
,
α(2m2 + (a + c)m + ac) σ(c + m) + c(c + 2m) ac(c + 2m)(a + 2m)
v=
αm(a + c + 2m) σm + c(c + 2m) ac(c + 2m)(a + 2m)
By (4.3), R0 is the larger root of the quadratic equation λ2 − (b1 + b2 )uλ + b1 b2 (u2 − v 2 ) = 0. It follows from the above expressions that u+v =
α + σa α + σ(a + 2m) , u−v = , ac (c + 2m)(a + 2m) 26
(6.1)
giving u2 − v 2 = (u + v)(v − v) = It follows from u + v =
α+σa ac
that
∂u ∂m
α + σa α + σ(a + 2m) . ac (c + 2m)(a + 2m)
∂v = − ∂m . Taking partial derivative with
respect to m for (6.1) gives ∂(u2 − v 2 ) ∂u ∂λ − b1 b2 (2λ − (b1 + b2 )u) ∂m = λ(b1 + b2 ) ∂m ∂m α + σa ∂u = λ(b1 + b2 ) − 2b1 b2 ac ∂m From the definition of u, α(4m + a + c) + σ(4m2 + 4am + a2 ) ∂u =− < 0. ∂m [(c + 2m)(a + 2m)]2 Note that 2λ|λ=R0 > (b1 + b2 )u. Therefore, to show that show
∂λ | ∂m λ=R0
∂R0 ∂m
< 0, i.e., to
< 0, it suffices to show that for λ = R0 , (b1 + b2 )λ − 2b1 b2
α + σa > 0. ac
We first claim that λ > 12 (b1 + b2 )(u + v). Since λ=
i p 1h (b1 + b2 )u + (b1 + b2 )2 u2 − 4b1 b2 (u2 − v 2 ) , 2
it is equivalent to show that
(b1 + b2 )2 u2 − 4b1 b2 (u2 − v 2 ) > (b1 + b2 )2 v 2 , or (b1 − b2 )2 (u2 − v 2 ) > 0. This is automatically true since u > v > 0. It follows from λ > 12 (b1 + b2 )(u + v) that 1 α + σa 1 (b1 + b2 )2 . (b1 + b2 )λ > (b1 + b2 )2 (u + v) = 2 2 ac 27
Thus as required (b1 + b2 )λ − 2b1 b2
α + σa 1 α + σa (b1 − b2 )2 > 0. ac 2 ac
Acknowledgements. YHH is supported by grant (NSC 94-2115-M005-006) from the National Science Council of Taiwan. YHH is also grateful to the Canadian government for their generous financial support to fund YHH’s visit to University of Victoria under a Faculty Research Award (623-2-FRP2005-04). PvdD is partially supported by NSERC of Canada, MITACS, and LW is supported by PIMS and MITACS PDF fellowships.
References Arino, J. and van den Driessche, P., 2006. Disease spread in metapopulations. Fields Institute Communications, in: Nonlinear Dynamics and Evolution Equations, 48. Arino, J. and van den Driessche, P., 2003. A multi-city epidemic model. Mathematical Population Studies, 10, 175-193. Arino, J., Jordan,R. and van den Driessche, P., 2005. Quarantine in a multispecies epidemic model with spatial dynamics. Math Biosci. Dec 9, [Epub ahead of print]. Berman, A. and Plemmons, R. J., 1979. Non-negative Matrices in the Mathematical Sciences. Academic Press, New York. Butler, D., 2006. Doubts hang over source of bird flu spread, Nature, 439, 772.
28
Chowell, G., Rivas, A. L., Hengartner, N. W., Hyman, J. M. and CastilloChavez, C., 2006. The role of spatial mixing in the spread of foot-and-mouth disease. Prev. Vet. Med. 73, 297-314. Ferguson, N. M., Cummings, D. A., Cauchemez, S., Fraser, C., Riley, S., Meeyai, A., Iamsirithaworn, S. and Burke, D. S., 2005. Strategies for containing an emerging influenza pandemic in Southeast Asia, Nature, 437, 209214. Fulford, G. R., Roberts, M. G. and Heesterbeek, J. A. P., 2002. The metapopulation dynamics of an infectious disease:tuberculosis in possums. Theor. Popul. Biol. 61, 15-29. Horn, R. A. and Johnson, C. R., 1991. Topics in Matrix Analysis, Cambridge University Press, New York. Horn, R. A. and Johnson, C. R., 1985. Matrix Analysis, Cambridge University Press, New York. Hyman, J. M. and LaForce, T., 2003. Modeling the spread of influenza among cities, 211-236, Bioterrorism (Edited by H.T. Banks and C. Castillo-Chavez), SIAM, New York. Ruan, S., Wang, W. and Levin, S. A., 2006. The effect of global travel on the spread of SARS, Math. Bios. and Eng. 3, 205-218. Salmani, M. and van den Driessche, P., 2006. A model for disease transmission in a patchy environment. Disc. and Cont. Dyna. Syst. -Series B, 6, 185-202. Sattenspiel, L. and Herring, D. A., 2003. Simulating the effect of quarantine on the spread of the 1918-19 flu in central Canada. Bull Math Biol. 65, 1-26. Smith, H. L. and Waltman, P., 1995. The Theory of the Chemostat, Cam29
bridge University Press, New York. van den Driessche, P. and Watmough, J., 2002. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math. Biosci. 180, 29-48. Wang, W. and Mulone, G., 2003. Threshold of disease transmission in a patchy environment, J. Math. Anal. Appl. 285, 321-335. Wang, W. and Zhao, X.-Q., 2004. An epidemic model in a patchy environment, Math. Biosci. 190, 97-112. World Health Organization (WHO), 2006a, Cumulative Number of Confirmed Human Cases of Avian Influenza A/(H5N1) Reported to WHO. Available at: http:// www. who. int/ csr/disease/avian influenza/country/ cases table 2006 04 19/en/index.html. Accessed April 19. World Health Organization (WHO), 2006b, Summary of WHO measures related to international travel. Available at: http:// www. who. int/csr/sars/ travelupdate/en/. Accessed April 4. World Organization for Animal Health (WOAH), 13 April 2006. UPDATE ON AVIAN INFLUENZA IN ANIMALS (TYPE H5). Available at: http://www.oie.int/downld/avian%20INFLUENZA/A AI-Asia.htm. Accessed April 19.
30
1
Corresponding author: P. van den Driessche Department of Mathematics and Statistics University of Victoria Victoria, BC, Canada Fax: +1-250 721 8962 Email: [email protected]
2
Abstract A multi-patch model is proposed to study the impact of travel on the spatial spread of disease between patches with different level of disease prevalence. (i)
The basic reproduction number for the ith patch in isolation, R0 , is obtained along with the basic reproduction number of the system of patches, R0 . Inequalities describing the relationship between these numbers are also given. For a two-patch model with one high prevalence patch and one low prevalence patch, results pertaining to the dependence of R0 on the travel rates between the two patches are obtained. These results show that, while banning travel from the low prevalence patch to the high prevalence patch always contributes to disease control, banning travel of symptomatic travelers only from the high prevalence patch to the low prevalence patch could adversely affect the containment of the outbreak under certain ranges of parameter values. Moreover, banning all travel from the high prevalence region to the low prevalence region could result in the low prevalence patch becoming disease-free, while the high prevalence patch becomes even more disease-prevalent, with the resulting number of infectives in patch 1 alone exceeding the combined number of infectives in both regions if border control had not been in place. Our results demonstrate that border control, if properly implemented, could be useful to stop the spatial spread of disease. Keywords: multi-patch model, spatial spread, basic reproduction number, travel rate, border control, influenza.
1
Introduction
Circulating in Asia since late 2003 and causing more than one hundred human deaths around the world by April 2006 (WHO, 2006a), the H5N1 avian influenza virus jumped to Africa and the Middle East, spread through Eu3
rope, affecting wild bird populations in 13 additional countries within a single month in early 2006 (WOAH, 2006). Seasonal migration alone probably cannot explain the westward spread in Europe, so imports of poultry and pet birds must also be considered as factors that might lead to the spread of H5N1 through countries and continents (Butler, 2006). Furthermore, its initial appearance on the African continent marks a huge leap in its geographical range, and opens up a whole new front where the vast bird reservoir could potentially spark a pandemic of human-to-human infections. The 2003 severe acute respiratory syndrome (SARS) global epidemic demonstrated the ability of infectious disease spreading due to modern globalization to countries in several continents within a matter of days. In its aftermath and with the potential threat of a flu pandemic, several models describing spatial spread of infectious diseases have been proposed. These include the use of multi-patch compartment models to explore the dynamics of spatial spread of disease(see, for example, Arino and van den Driessche, 2003, 2006; Arino et al. 2005; Hyman and LaForce, 2003; Ruan et. al., 2006; Salmani and van den Driessche, 2006; Sattenspiel and Herring, 2003; Wang and Mulone, 2003; Wang and Zhao, 2004). Moreover, previous modeling work has shown that compartmental models incorporating multiple patches (cities, countries, etc.) describe epidemic spread more accurately than non-spatial models, even at early epidemic phases. Examples of such include foot-and-mouth disease (Chowell et. al., 2006), tuberculosis in possums (Fulford et. al., 2002), in addition to models of infectious diseases of humans. In this work we propose a multi-patch model to study the spread of influenza among the patches. We add a subpopulation of partially immune individuals to account for this important feature of influenza. Although the model is intended to be used for theoretical studies of the spread of human influenza, it can also be used as a model for studying the spread of enzootic diseases such as avian flu (H5N1) among birds. Partial immunity might 4
not be as important for studies of infectious disease involving poultry bird populations where the birds are slaughtered for food after a fixed time, it nonetheless could be an important consideration for studies involving wild bird populations. The paper is organized as follows. We formulate the model in Section 2, and in Section 3, we give results regarding the basic reproduction number of the model. In Section 4, some global analysis is obtained for the model with two patches. Numerical simulations, with parameter values relevant for influenza, to complement our analytical results are given in Section 5. Finally in the last section, from our results, we discuss the impact of travel on the spatial spread of disease.
2
Model Formulation
In this section we formulate a model describing the spread of a disease in a population with n patches taking into consideration travel among patches (Arino and van den Driessche, 2006) and also including partially immune individuals (Hyman and LaForce, 2003). The population in patch i is divided into compartments of susceptible, incubating (infected but not yet showing symptoms), infective (infected with symptoms), recovered and partially immune individuals. Let Si , Ei , Ii , Ri , Pi denote, respectively, the associated population size. Then the total population size in patch i is Ni = Si + Ei + Ii + Ri + Pi for i = 1, 2, . . . , n. A partially immune compartment between the recovered and the susceptible compartments is introduced in an influenza model (Hyman and LaForce, 2003) to account for people who have partial immunity to the current strain of influenza from a previous infection by an earlier strain. We include this partially immune compartment, and assume that individuals in this compartment may become infected again (but with a reduced rate). Also the partial immunity wanes, with these 5
individuals returning to the susceptible compartment. For patch i, let Ai be the recruitment, αi , γi , δi and ηi be the progression rates of the incubating, infective, recovered and partially immune individuals, respectively, di > 0 be the natural death rate and ǫi be the disease-related death rate. We assume that individuals do not change their disease state during travel, and mK ij for K = S, E, I, R, P are the constant travel rates from patch j to patch i for i 6= j of susceptible, incubating, infective, recovered and partially immune individuals, respectively, with mK ii = 0. All parameters are assumed to be positive except that ǫi , δi , ηi can be zero, and the travel rate matrices M K = [mK ij ] for K = S, E, I, R, P are irreducible. For more detailed discussion on travel rates, see, e.g., Arino and van den Driessche (2006). The number of new individuals infected by infectives per unit time in patch i is given by βi (Ni )Si Ii . The term βi (Ni )Si is the product of βi (Ni )Ni , the average number of contacts made by each individual in patch i per unit time, and
Si , Ni
the proportion of susceptibles. It is assumed that infectivity βi (Ni ) is
a continuously differentiable, non-increasing function of Ni with βi (0) < ∞. Note that the above assumptions encompass the widely used mass action and standard incidence disease transmission terms defined for Ni > 0 as well as many other forms of saturating incidence. Reduced infectivity of incubating and partially immune individuals is given by σi βi (Ni ) and νi βi (Ni ), respectively, where σi , νi ∈ [0, 1) and it is assumed that infective individuals can reinfect those who are partially immune. The model flow chart for patch i omitting travel and natural death is given in Figure 1. The above assumptions lead to the following SEIRP model for i = 1, 2, ..., n:
6
dSi dt
= Ai − βi (Ni )Si (Ii + σi Ei ) − di Si + ηi Pi +
n X
mSij Sj
−
j=1
dEi dt
= βi (Ni )[Si (Ii + σi Ei ) + νi Pi Ii ] − (di + αi )Ei +
n X j=1
n X
mE ij Ej
−
dIi dt
= αi Ei − (γi + εi + di )Ii +
mIij Ij −
dRi dt
= γi Ii − (di + δi )Ri +
mR ij Rj −
mIji Ii
n X
mR ji Ri
j=1
j=1
dPi dt
mE ji Ei
j=1
j=1
n X
n X
n X j=1
j=1
n X
mSji Si
= δi Ri − (di + ηi )Pi − νi βi (Ni )Pi Ii +
n X j=1
mPij Pj
−
n X
mPji Pi
j=1
(2.1)
Figure 1: Flow diagram of the model. Initially each variable is assumed to be non-negative with Si (0) > 0 and Pn i=1 Ei (0) + Ii (0) > 0. It follows that for a given set of non-negative initial
conditions, there is a unique solution to system (2.1). The total population Pn ¯ size in all n patches is N (t) = i=1 Ni (t). Let d = min{d1 , d2 , . . . , dn } Pn and A = i=1 Ai . Then the following result, which can be proved in a similar way to that of Theorem 1.1 in Salmani and van den Driessche (2006),
indicates that the model is well posed and all variables remain non-negative and bounded. 7
Theorem 2.1. Consider system (2.1) with non-negative initial conditions. Then for each i = 1, 2, . . . , n, Ei (t), Ii (t), Ri (t), Pi (t) remain non-negative, Si (t) and Ni (t) remain positive and the total population N (t) is in the interval ¯ N (0)}]. (0, max{A/d, Notice that our patch model is an extension of those in Arino and van den Driessche (2006), Salmani and van den Driessche (2006). It includes a new compartment, namely, partially immune individuals in each patch, and also the probability of disease transmission by incubating individuals. These additional features may be important for modeling diseases such as influenza and SARS. In Hyman and LaForce (2003), which does not include an incubating compartment in a model for the spread of influenza, travel is assumed to be independent of disease status and symmetric. Moreover, the disease does not cause death, and the population of each patch remains constant. Notice also that our patch model keeps track of individuals present in patch i at time t, but does not keep track of where an individual resides. Models that include this are developed in Arino and van den Driessche (2003, 2006) and Sattenspiel and Herring (2003), and Ruan et. al. (2006) use such a model to study the effect of travel on the spread of SARS.
3
The basic reproduction number
In this section, we derive a formula to compute the basic reproduction number R0 for the general model (2.1) and then give a lower bound for R0 . Both lower and upper bounds of R0 are given for a special case. A disease free equilibrium (DFE) is a steady state solution of system (2.1) with Si = Ni∗ > 0 and all other variables Ei , Ii , Ri , Pi equal to 0 for i = 1, 2, . . . , n. Let S ∗ = (N1∗ , N2∗ , . . . , Nn∗ )T . Then from the first equation in (2.1), there is a DFE if and only if S ∗ satisfies CS ∗ = A with A = P (A1 , A2 , . . . , An )T and C = D − M S , where D = diag( nj=1 mSji + di ). Note 8
that C is irreducible, has positive column sums d1 , d2 , . . . , dn and negative off-diagonal entries. Thus C is a non-singular M -matrix (page 141, Berman and Plemmons, 1979), and therefore C −1 > 0. Hence S ∗ = C −1 A > 0 is the unique solution of CS ∗ = A. This shows that the DFE always exists and is unique. In the absence of disease system (2.1) reduces to just the first equation and S ∗ is stable. Next we consider the local stability of the DFE for system (2.1). To this end, we order the infected variables by E1 , E2 , . . . , En , I1 , I2 , . . . , In and make use of the result in van den Driessche and Watmough (2002) to obtain # " # " diag(σi βi (Ni∗ )Ni∗ ) diag(βi (Ni∗ )Ni∗ ) F11 F12 = F = 0 0 0 0 with V =
"
V11
0
−V21 V22
#
,
where V11 = diag(di + αi +
n X
E mE ji ) − M ,
j=1
V22 = diag(γi + ǫi + di +
n X
mIji ) − M I ,
j=1
V21 = diag(αi ), i = 1, 2, . . . , n. Note that V11 and V22 are both irreducible non-singular M −matrices with positive column sums and hence V11−1 > 0, V22−1 > 0.
(3.1)
The basic reproduction number for the system, denoted by R0 is then the spectral radius of F V −1 , i.e., R0 = ρ{F V −1 } where " # #" −1 F F V 0 11 12 11 F V −1 = . −1 0 0 V22 V21 V11−1 V22−1 9
Therefore R0 = ρ(F11 V11−1 + F12 V22−1 V21 V11−1 ).
(3.2)
The first term accounts for infection from incubating individuals, while the second term accounts for infection from infective individuals who survive the incubating compartment. Travel rates influence the average time spent in the incubating and infective compartments. Note that R0 does not depend on νi , δi , ηi for i = 1, 2, . . . , n. The basic reproduction number gives an important threshold for the disease, as shown in the following result. Theorem 3.1. Consider model (2.1). If R0 < 1, then the DFE is locally asymptotically stable and if R0 > 1, the DFE is unstable. Moreover, if the disease transmission is standard incidence, then the DFE is globally asymptotically stable provided that R0 < 1. Proof. It follows from Theorem 2 of van den Driessche and Watmough (2002) that the DFE is locally asymptotically stable if R0 < 1 and is unstable if R0 > 1. If the disease transmission is standard incidence, then βi (Ni )Ni = βi for i = 1, 2, . . . , n. Note that Si ≤ Ni and Si + νi Pi ≤ Ni . This gives the inequality n
n
X X dEi mE mE E − ≤ βi Ii + σi βi Ei − (di + αi )Ei + ji Ei . ij j dt j=1 j=1 Consider the linear system dEi dt
= βi Ii + σi βi Ei − (di + αi )Ei +
n X
mE ij Ej −
dIi dt
= αi Ei − (γi + εi + di )Ii +
j=1
mIij Ij
mE ji Ei
j=1
j=1
n X
n X
−
n X
(3.3)
mIji Ii .
j=1
The right hand side of the above system has F − V as its coefficient matrix. Again, by proof of Theorem 2 of van den Driessche and Watmough (2002), 10
each eigenvalue of F − V has negative real part if R0 < 1. Thus any solution of (3.3) satisfies lim Ei = 0 and lim Ii = 0 for i = 1, 2, . . . , n. Using a t→∞
t→∞
comparison theorem (Thorem B.1, Smith and Waltman, 1995), each solution of system (2.1) satisfies lim Ei = 0 and lim Ii = 0 for i = 1, 2, . . . , n. Using t→∞
t→∞
a similar argument as in the proof of Theorem 2.2, Salmani and van den Driessche (2006), lim Ri (t) = 0, and similarly lim Pi (t) = 0. Thus the DFE t→∞
t→∞
is globally asymptotically stable provided that R0 < 1. Let ai = γi + ǫi + di , bi = βi (Ni∗ )Ni∗ , ci = di + αi for i = 1, 2, . . . , n. In the case that there is no travel between patch i and all other patches, the basic reproduction number in patch i in isolation is given by (i)
R0 =
σi bi bi αi + . ci ai c i
(3.4)
Define a modified reproduction number (modified by travel) in patch i, e (i) R 0 =
ci +
σb Pini
E j=1 mji
+
(ai +
bi αi I j=1 mji )(ci
Pn
+
Pn
j=1
mE ji )
.
(3.5)
The following result gives bounds for the basic reproduction number for system (2.1) in terms of the numbers defined in (3.4) and (3.5) for each patch. Theorem 3.2. For system (2.1), e (i) R0 ≥ max R 0 .
(3.6)
1≤i≤n
Furthermore, if ai = a, αi = α, σi = σ and di = d for i = 1, 2, . . . , n, then (i) (i) e0 , min R0 ≤ R0 ≤ max R(i) (3.7) max max R 0 . 1≤i≤n
1≤i≤n
1≤i≤n
Proof For j = 1, 2, let Vjj [1′ ] denote the matrix Vjj with row and column
1 deleted, Y = [yij ] and Z = [zij ] denote V11−1 and V22−1 , respectively. Let W = [wij ] = G + H, where G = [gij ] = F11 Y , H = [hij ] = F12 ZV21 Y. It follows from (3.1) that yij > 0, zij > 0 for i, j = 1, 2, . . . , n. Then, by Corollary 8.1.20, Horn and Johnson (1985), R0 = ρ(W ) ≥ wii = gii + hii for i = 1, 2, . . . , n. 11
′
11 [1 ] Note that g11 = σ1 b1 y11 = σ1 b1 detV . By virtue of Fischer’s inequality detV11
(page 117, Horn and Johnson, 1991), it follows that ! n X ′ mE detV11 ≤ c1 + j1 detV11 [1 ]. j=1
Therefore, g11 ≥ Similarly,
c1 +
σb P1n1
h11 = b1 α1 z11 y11 +
j=1
n X
mE j1
.
b1 αk z1k yk1 ,
k=2
thus,
detV22 [1′ ] detV11 [1′ ] . detV22 detV11 Again, by Fischer’s inequality, it follows that h11 ≥ b1 α1 z11 y11 = b1 α1
b1 α1 . h11 ≥ P Pn a1 + j=1 mIj1 c1 + nj=1 mE j1
Adding these shows that
(1)
e0 . R0 ≥ R
Similarly, it can be shown that
(i)
and this gives (3.6).
e0 R0 ≥ R
for i = 2, 3, . . . , n
If ai = a, αi = α, σi = σ and di = d (thus ci = c) for i = 1, 2, . . . , n, P then wij = σbi yij + αbi nk=1 zik ykj for i, j = 1, 2, . . . , n. Without loss of
generality, assume that 0 < b1 ≤ b2 ≤ · · · ≤ bn . From the fact that the
matrix V11 has each column sum equal to c > 0 and the matrix V22 has P Pn each column sum equal to a > 0, it follows that ni=1 yij = 1c , i=1 zij = 12
1 a
for j = 1, 2, . . . , n. Therefore, for the matrix W , the sum of column j is
given by n X
wij =
n X
σbi yij +
≤ σbn = σbn
αbi
i=1
i=1
i=1
n X
n X
i=1 n X
yij + αbn yij + αbn
i=1
n X
zik ykj
i=1 k=1 n n X X k=1
σbn αbn + c ac (n) = R0 .
zik ykj
k=1 n n XX
zik
i=1
!
ykj
=
Similarly,
Pn
i=1
wij ≥
σb1 c
+
αb1 ac
(1)
= R0 . From Theorem 8.1.22, Horn and
Johnson (1985), ρ(W ) lies between the minimum and maximum column sums of W , thus (i)
(i)
min R0 ≤ R0 = ρ(W ) ≤ max R0 .
1≤i≤n
1≤i≤n
Combining with (3.6) immediately gives the desired (3.7) The analytical results give the basic reproduction numbers in isolation of the individual patches as upper and lower bounds for the basic reproduction number of the system. However, its usefulness is rather limited in practice since the range for the bounds could be too large. In the next section, we consider the simple case of only two patches, and obtain more explicit results giving insight about the impact of travel on the spatial spread of disease between patches.
13
4
The Model with 2 Patches
From now on we consider the special case of (2.1) with only 2 patches, i.e., we consider the system dS1 dt dE1 dt dI1 dt dR1 dt dP1 dt dS2 dt dE2 dt dI2 dt dR2 dt dP2 dt
= A1 − β1 (N1 )S1 (I1 + σ1 E1 ) − d1 S1 + η1 P1 + mS12 S2 − mS21 S1 E = β1 (N1 )[S1 (I1 + σ1 E1 ) + ν1 P1 I1 ] − (d1 + α1 )E1 + mE 12 E2 − m21 E1
= α1 E1 − (γ1 + ε1 + d1 )I1 + mI12 I2 − mI21 I1 R = γ1 I1 − (d1 + δ1 )R1 + mR 12 R2 − m21 R1
= δ1 R1 − (d1 + η1 )P1 − ν1 β1 (N1 )P1 I1 + mP12 P2 − mP21 P1 = A2 − β2 (N2 )S2 (I2 + σ2 E2 ) − d2 S2 + η2 P2 + mS21 S1 − mS12 S2 E = β2 (N2 )[S2 (I2 + σ2 E2 ) + ν2 P2 I2 ] − (d2 + α2 )E2 + mE 21 E1 − m12 E2
= α2 E2 − (γ2 + ε2 + d2 )I2 + mI21 I1 − mI12 I2 R = γ2 I2 − (d2 + δ2 )R2 + mR 21 R1 − m12 R2
= δ2 R2 − (d2 + η2 )P2 − ν2 β2 (N2 )P2 I2 + mP21 P1 − mP12 P2 (4.1)
The DFE of (4.1) is given by S∗ =
"
N1∗ N2∗
#
=
"
d1 + mS21
−mS12
−mS21
d2 + mS12
#−1 "
A1 A2
#
.
Also F11 = diag(σ1 b1 , σ2 b2 ), F12 = diag(b1 , b2 ), V21 = diag(α1 , α2 ), and V11 =
"
c1 + m E 21
−mE 12
−mE 21
c2 + m E 12
#
, V22 =
"
a1 + mI21
−mI12
−mI21
a2 + mI12
#
.
Hence R0 = ρ(W ), where W = F11 V11−1 + F12 V22−1 V21 V11−1 = 14
"
w11 w12 w21 w22
#
(4.2)
with w11 =
I I E σ1 b1 (c2 + mE α1 b1 (c2 + mE 12 ) 12 )(a2 + m12 ) + α2 b1 m12 m21 + E E I I c1 c2 + c1 mE (c1 c2 + c1 mE 12 + c2 m21 12 + c2 m21 )(a1 a2 + a1 m12 + a2 m21 )
w12 =
I I E σ1 b1 mE α1 b1 mE 12 12 (a2 + m12 ) + α2 b1 m12 (c1 + m21 ) + E E I I c1 c2 + c1 mE (c1 c2 + c1 mE 12 + c2 m21 12 + c2 m21 )(a1 a2 + a1 m12 + a2 m21 )
w21 =
E I σ2 b2 mE α1 b2 mI21 (c2 + mE 21 12 ) + α2 b2 m21 (a1 + m21 ) + E E I I c1 c2 + c1 mE (c1 c2 + c1 mE 12 + c2 m21 12 + c2 m21 )(a1 a2 + a1 m12 + a2 m21 )
w22
I E I α1 b2 mE σ2 b2 (c1 + mE 12 m21 + α2 b2 (c1 + m21 )(a1 + m21 ) 21 ) . + = I I E E c1 c2 + c1 mE (c1 c2 + c1 mE 12 + c2 m21 12 + c2 m21 )(a1 a2 + a1 m12 + a2 m21 )
It follows from (3.2) that R0 =
p 1 w11 + w22 + (w11 − w22 )2 + 4w12 w21 . 2
(4.3)
From Theorem 3.1, we have the following result.
Theorem 4.1. For the two-patch model, assume that the disease transmission is standard incidence. Then the DFE is globally asymptotically stable if R0 < 1 and is unstable if R0 > 1, where R0 is given by (4.3). As expected, an increase in σi or bi increases R0 . To investigate how R0 changes with the other parameters, we first assume that all travel rates for I incubating and infective individuals are equal, namely, mE ij = mij = m for
i, j = 1, 2. The proof of the following result is given in the Appendix. I Theorem 4.2. If mE ij = mij = m for i, j = 1, 2 and σ1 = σ2 = σ, α1 = α2 =
α, a1 = a2 = a, c1 = c2 = c, b1 6= b2 , then an increase in a or m decreases R0 . Remark 4.1. With the conditions of the above theorem, as the travel rate becomes large, the basic reproduction number R0 in (4.3) approaches the (1)
(2)
(1)
(2)
mean value of R0 and R0 , i.e. R0 → 12 (R0 + R0 ) as m → ∞. This can be seen from the form of W in the proof of Theorem 4.2. In this limit case, the two patches merge into one. 15
Using a similar technique as in proof of Theorem 4.2, we can prove the following result in the case that infectives of both patches are too sick to travel. I Theorem 4.3. If mE ij = m, mij = 0 for i, j = 1, 2 and σ1 = σ2 , α1 = α2 , a1 =
a2 , c1 = c2 , b1 6= b2 , then
∂R0 ∂m
< 0, i.e., the basic reproduction number R0
decreases as m increases.
5
Numerical Simulations
Assuming standard incidence disease transmission, we present some numerical simulations for 2 patches to illustrate how R0 changes with travel rates, with the choice of parameters relevant for human influenza (Ferguson et. al., 2005). In the next section, we present simulations for restricted travel. The model parameters with time unit as one day are taken as: α1 = α2 = 1 1.48
≈ 0.68, γ1 = γ2 =
1 2.6
≈ 0.38 (the average incubating time is 1.48 days
and the average infective time is 2.6 days), ǫ1 = ǫ2 = 4.0 × 10−3 , d1 = d2 =
1 70×365
≈ 3.91 × 10−5 , b1 = 0.6, b2 = 0.1, σ1 = σ2 = 0.04. These
parameters yields the respective basic reproduction numbers in isolation of (1)
(2)
R0 ≈ 1.58 > 1 and R0 ≈ 0.26 < 1. Hence we consider the hypothetical scenario of disease spread between a high-prevalence endemic region (patch 1) and a low-prevalence region where a minor outbreak could be eradicated I (patch 2). We first keep mE 21 = m21 fixed with values at 0.6, 0.8 and 1. We I let mE 12 = m12 vary from 0 to 1. The curves of R0 are given in Figure 2. A I E I 3-d plot of R0 versus mE 12 = m12 and m21 = m21 is given in Figure 3.
Figure 2 shows that by keeping the same travel rate from patch 1 to patch 2, an increase in the travel rate from patch 2 (the patch with lower disease transmission rate) to patch 1 (the patch with higher disease transmission rate) results in an increase in R0 . Consequently, we can conclude that travel between the two patches may cause the disease to become endemic in both 16
1.2 1.1
E
m21=0.6
1
mE =0.8 21
0.9
R
0
0.8
E
m21=1.0
0.7 0.6 0.5 0.4 0.3 0.2 0
0.2
0.4
0.6 E 12
0.8
1
I 12
m =m
I E I Figure 2: R0 vs mE 12 = m12 for fixed m21 = m21 . The three curves from top I to bottom correspond to : mE 21 = m21 = 0.6, 0.8, 1.0
2
R
0
1.5 1 0.5 0 1 1 0.5 E 12
I 12
m =m
0.5 0
0
E 21
I 21
m =m
I E I Figure 3: R0 vs mE 12 = m12 and m21 = m21 .
17
patches if the travel rate from patch 2 to patch 1 is large enough. The numerical simulation in Figure 3 shows that, for the chosen parameters and I E I keeping mE 12 = m12 fixed, an increase in m21 = m21 leads to a decrease in
R0 . In particular, R0 could conceivably decrease to a value less than one if the travel rate from the high-prevalence patch to the low-prevalence patch, I namely mE 21 = m21 , is high enough.
If all travel rates from one patch to the other are the same, i.e. mE 12 = I mI12 = mE 21 = m21 = m, then an increase in m leads to a decrease in R0
(Figure 4), as predicted by Theorem 4.2. Note that, for these parameter values, R0 < 1 if m > 0.4. 1.6 1.5 1.4
R0
1.3 1.2 1.1 1 0.9 0
0.2
0.4
0.6
0.8
1
m I E I Figure 4: R0 vs m, where m = mE 12 = m12 = m21 = m21 .
Taking m = 0.5 and subsequently R0 < 1, the disease will die out eventually. To investigate how the infective population sizes in the two patches change over time, we numerically simulate the system (4.1) with i = 1, 2, and plot I1 and I2 against time in Figure 5. In addition to the parameter values given at the beginning of this section, we take A1 = A2 = 100, η1 = η2 = δ1 = δ2 = 0.2, ν1 = ν2 = 0.1. The infective population sizes in both 18
patches increase initially then slowly decrease to zero since the DFE is globally asymptotically stable if R0 < 1. 70
60
I1 and I2
50
40
30
I
1
I2
20
10
0
0
20
40
60
80
100
t
Figure 5: Numerical solution of system (4.1) with parameters as in the text and m = 0.5 showing I1 and I2 vs time t. Initial conditions are: S1 (0) = 400, E1 (0) = 50, I1 (0) = 20, R1 (0) = P1 (0) = 0; S2 (0) = 800, E2 (0) = 30, I2 (0) = 10, R2 (0) = P2 (0) = 0. Now we use the same parameter values except that m = 0.2, and subsequently R0 > 1. Numerical simulation of this case given in Figure 6 shows that in both patches there is an initial increase in the number of infectives before the infective populations decrease to their endemic levels. However, the existence and stability of such an endemic equilibrium remain unproved analytically. Some interesting observations can be made from our results regarding the role that travel plays in the spatial spread of a disease. Figures 2 and 3 demonstrate the possibility that, for a lower-prevalence patch with a minor (2)
disease outbreak (basic reproduction numbers in isolation R0 less than one) open travel with a high-prevalence patch could lead to the disease becoming 19
45 40
I1 and I2
35
I
1
30 25
I
2
20 15 10 0
20
40
60
80
100
t
Figure 6: Numerical solution of system (4.1) showing I1 and I2 vs time t. Parameter values and initial conditions are the same as in Figure 5 except that m = 0.2. endemic. On the other hand, for a high-prevalence patch with endemic dis(1)
ease in isolation (basic reproduction numbers in isolation R0 greater than one), open travel with a high-prevalence patch could eradicate the disease. Essentially, under appropriate parameter ranges, travel between patches dilutes the overall prevalence to the point that it could either lessen the severity of an endemic patch or worsen a minor outbreak region. Further evidence can be found in Figures 4 and 5 in which assuming all travel rates are equal, the disease can be eventually eradicated if the travel rates are sufficiently large.
6
Discussion of Travel Restrictions
To consider hypothetical intervention scenarios, we let mI21 = 0 and keep all other parameters the same as in Figure 2. This models a situation in which the authority bans all travel of symptomatic travelers from patch 1, the high 20
prevalence region, to patch 2, the low prevalence region. In comparison with Figure 2, Figure 7 shows this gives an overall increase in the value of the basic reproduction number R0 . More significantly, there is a range of parameter values for mI12 = mE 12 in which stopping all travel of symptomatic travelers from patch 1 to patch 2 could adversely impact the epidemic by driving the basic reproduction number above one, thus prolonging the epidemic that otherwise would be eradicated. 1.5 mE =0.6 21
1.4 1.3 1.2
mE =0.8 21
mE =1.0 21
R
0
1.1 1 0.9 0.8 0.7 0
0.2
0.4
E 12
I 0.6 12
m =m
0.8
1
I I Figure 7: R0 vs mE 12 = m12 for fixed m21 = 0. The three curves from top to
bottom correspond to : mE 21 = 0.6, 0.8, 1.0 Conversely, Figure 8 shows that stopping all travel of symptomatic travelers from patch 2 to patch 1 (letting mI12 = 0) would alleviate the epidemic, lowering the basic reproduction number to below one, thus successfully preventing further spread of the outbreak. Therefore, the policy of border control to ban travel of symptomatic travelers only from the high prevalence patch to the low prevalence patch could affect the containment of the outbreak adversely (Figures 2 and 7). However, banning travel of symptomatic travelers only from the low prevalence 21
0.7 0.65 0.6
R0
0.55 0.5 0.45 0.4 0.35
0
0.2
0.4
mE
0.6
0.8
1
12
I E I Figure 8: R0 vs mE 12 with m12 = 0 for fixed m21 = m21 = 0.6
. patch to the high prevalence patch always has a positive impact (Figures 2 and 8). We note here that by banning all travels (both exposed but asymptomatic and symptomatic individuals) would result in an outbreak in the high prevalence patch. We now consider the case in which the disease is endemic, as in Figure 6. We again let mI21 = 0 and keep all other parameters the same as in Figure 6. Figure 9 shows that the disease is still endemic in both patches, but with significantly larger numbers of infectives in both patches. Thus once again, banning all travel of symptomatic travelers from the high prevalence region to the low prevalence region is detrimental to the intervention and control of the outbreak. If we let mI21 = mE 21 = 0, i.e., banning all travel from the high prevalence region to the low prevalence region (Figure 10), patch 2 becomes disease-free, while patch 1 is even more prevalent with the number of infectives in patch 1 alone exceeding the total number of infectives in the two regions if border control had not been in place. Note that in this case the infected equations for patch 2 uncouple. Matrix W given by (4.2) has 22
100
80 I
I1 and I2
1
60
40 I
2
20
0 0
20
40
60
80
100
t
Figure 9: Numerical solution of system (4.1) with two patches showing I1 and I2 vs time t. Parameters as in Figure 6 except that mI21 = 0. w21 = 0 and so is reducible. Thus the reproduction number in patch 1 is (1) e (2) R0 > 1 and the reproduction number in patch 2 is R 0 < 1.
Now we let mI12 = 0 and keep all other parameters the same as in Figure 6,
i.e., banning travel of symptomatic travelers from the low prevalence patch
to the high prevalence patch only, the disease is eradicated in both regions (Figure 11). This indicates the importance of border control out of a low prevalence patch and into a high prevalence patch. We could also consider a slightly different scenario in which both patches have the same infectivity, but initially one patch is disease free while the other patch has infectives, i.e., E2 (0) = I2 (0) = 0 and E1 (0) + I1 (0) > 0. However, since the analytical results do not depend on the initial conditions, the resulting simulations are asymptotically similar to Figures 9-11. During the 2003 SARS outbreak, travel warnings to all affected areas were issued by WHO to prevent travelers from entering and becoming infected. Dozens of countries also issued border control either banning travelers from
23
140 120
I1 and I2
100
I
1
80 60 40 20 0 0
I
2
20
40
60
80
100
t
Figure 10: Numerical solution of system (4.1) with two patches showing I1 and I2 vs time t. Parameters are the same as in Figure 6 except that mI21 = mE 21 = 0.
35 30 I
I1 and I2
25
1
20 I
15
2
10 5 0 0
20
40
60
80
100
t
Figure 11: Numerical solution of system (4.1) with two patches showing I1 and I2 vs time t. Parameters as in Figure 6 except that mI12 = 0.
24
entering or placing them under quarantine. From a simple public health point of view, it is imperative for the low-prevalence region to stop travel of sick (both exposed and asymptomatic) individuals from the high-prevalence region for intervention purposes. See (Ruan et. al., 2006) for a study of a patch model for the spread of SARS. Our results show that border control does not necessarily always have a positive impact on the overall spread of disease and it is more important to ban travel from a low prevalence patch to a high prevalence patch. Moreover, an increase of travel rates in the opposite direction (from a high prevalence patch to a low prevalence patch), while theoretically alleviating the spatial spread of the disease, is not likely an implementable policy. Furthermore, letting mI12 = mE 12 = 0, i.e., banning all travel from the low prevalence patch to the high prevalence patch, numerical simulations in Figure 12 show R0 decreasing drastically as mI21 = mE 21 increase from 0 to 1. Note that in this case the infected equations in patch 1 uncouple, and the matrix W in (4.2) is reducible, thus (2) e (1) R0 = max{R 0 , R0 }. On the other hand, banning all travel from the high
prevalence patch to the low prevalence patch (mI21 = mE 21 = 0) results in (1)
R0 = R0 ≈ 1.58 > 1 and the outbreak becomes endemic. On the basis of our model assumptions and the parameters considered, we could therefore conclude that, retrospectively, WHO acted properly with travel warnings for travelers to avoid all but essential travel to affected areas (WHO, 2006b). Moreover, screening at border for travelers in and out is important mainly for the purpose of quick identification of sick individuals. As a final remark, we note again that, while the model is proposed with the spread of human influenza in mind, it also can be used, with some appropriate modifications and parameter changes, as the basis of a theoretical model to study the spread of enzootic diseases such as avian flu (H5N1) among birds, both that of importation of poultry and pet birds, as well as wild migratory bird populations. 25
1.6 1.4 1.2
R
0
1 0.8 0.6 0.4 0.2 0
0.2
0.4
0.6 I
0.8
1
E
m21=m21
E I Figure 12: R0 vs mI21 = mE 21 for m12 = m12 = 0.
Appendix: Proof of Theorem 4.2 Proof. It is clear that an increase in a decreases each entry in the matrix W , and so decreases the basic reproduction number R0 . Under the parameter assumptions W = where u=
"
b1 u b 1 v b2 v b 2 u
#
,
α(2m2 + (a + c)m + ac) σ(c + m) + c(c + 2m) ac(c + 2m)(a + 2m)
v=
αm(a + c + 2m) σm + c(c + 2m) ac(c + 2m)(a + 2m)
By (4.3), R0 is the larger root of the quadratic equation λ2 − (b1 + b2 )uλ + b1 b2 (u2 − v 2 ) = 0. It follows from the above expressions that u+v =
α + σa α + σ(a + 2m) , u−v = , ac (c + 2m)(a + 2m) 26
(6.1)
giving u2 − v 2 = (u + v)(v − v) = It follows from u + v =
α+σa ac
that
∂u ∂m
α + σa α + σ(a + 2m) . ac (c + 2m)(a + 2m)
∂v = − ∂m . Taking partial derivative with
respect to m for (6.1) gives ∂(u2 − v 2 ) ∂u ∂λ − b1 b2 (2λ − (b1 + b2 )u) ∂m = λ(b1 + b2 ) ∂m ∂m α + σa ∂u = λ(b1 + b2 ) − 2b1 b2 ac ∂m From the definition of u, α(4m + a + c) + σ(4m2 + 4am + a2 ) ∂u =− < 0. ∂m [(c + 2m)(a + 2m)]2 Note that 2λ|λ=R0 > (b1 + b2 )u. Therefore, to show that show
∂λ | ∂m λ=R0
∂R0 ∂m
< 0, i.e., to
< 0, it suffices to show that for λ = R0 , (b1 + b2 )λ − 2b1 b2
α + σa > 0. ac
We first claim that λ > 12 (b1 + b2 )(u + v). Since λ=
i p 1h (b1 + b2 )u + (b1 + b2 )2 u2 − 4b1 b2 (u2 − v 2 ) , 2
it is equivalent to show that
(b1 + b2 )2 u2 − 4b1 b2 (u2 − v 2 ) > (b1 + b2 )2 v 2 , or (b1 − b2 )2 (u2 − v 2 ) > 0. This is automatically true since u > v > 0. It follows from λ > 12 (b1 + b2 )(u + v) that 1 α + σa 1 (b1 + b2 )2 . (b1 + b2 )λ > (b1 + b2 )2 (u + v) = 2 2 ac 27
Thus as required (b1 + b2 )λ − 2b1 b2
α + σa 1 α + σa (b1 − b2 )2 > 0. ac 2 ac
Acknowledgements. YHH is supported by grant (NSC 94-2115-M005-006) from the National Science Council of Taiwan. YHH is also grateful to the Canadian government for their generous financial support to fund YHH’s visit to University of Victoria under a Faculty Research Award (623-2-FRP2005-04). PvdD is partially supported by NSERC of Canada, MITACS, and LW is supported by PIMS and MITACS PDF fellowships.
References Arino, J. and van den Driessche, P., 2006. Disease spread in metapopulations. Fields Institute Communications, in: Nonlinear Dynamics and Evolution Equations, 48. Arino, J. and van den Driessche, P., 2003. A multi-city epidemic model. Mathematical Population Studies, 10, 175-193. Arino, J., Jordan,R. and van den Driessche, P., 2005. Quarantine in a multispecies epidemic model with spatial dynamics. Math Biosci. Dec 9, [Epub ahead of print]. Berman, A. and Plemmons, R. J., 1979. Non-negative Matrices in the Mathematical Sciences. Academic Press, New York. Butler, D., 2006. Doubts hang over source of bird flu spread, Nature, 439, 772.
28
Chowell, G., Rivas, A. L., Hengartner, N. W., Hyman, J. M. and CastilloChavez, C., 2006. The role of spatial mixing in the spread of foot-and-mouth disease. Prev. Vet. Med. 73, 297-314. Ferguson, N. M., Cummings, D. A., Cauchemez, S., Fraser, C., Riley, S., Meeyai, A., Iamsirithaworn, S. and Burke, D. S., 2005. Strategies for containing an emerging influenza pandemic in Southeast Asia, Nature, 437, 209214. Fulford, G. R., Roberts, M. G. and Heesterbeek, J. A. P., 2002. The metapopulation dynamics of an infectious disease:tuberculosis in possums. Theor. Popul. Biol. 61, 15-29. Horn, R. A. and Johnson, C. R., 1991. Topics in Matrix Analysis, Cambridge University Press, New York. Horn, R. A. and Johnson, C. R., 1985. Matrix Analysis, Cambridge University Press, New York. Hyman, J. M. and LaForce, T., 2003. Modeling the spread of influenza among cities, 211-236, Bioterrorism (Edited by H.T. Banks and C. Castillo-Chavez), SIAM, New York. Ruan, S., Wang, W. and Levin, S. A., 2006. The effect of global travel on the spread of SARS, Math. Bios. and Eng. 3, 205-218. Salmani, M. and van den Driessche, P., 2006. A model for disease transmission in a patchy environment. Disc. and Cont. Dyna. Syst. -Series B, 6, 185-202. Sattenspiel, L. and Herring, D. A., 2003. Simulating the effect of quarantine on the spread of the 1918-19 flu in central Canada. Bull Math Biol. 65, 1-26. Smith, H. L. and Waltman, P., 1995. The Theory of the Chemostat, Cam29
bridge University Press, New York. van den Driessche, P. and Watmough, J., 2002. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math. Biosci. 180, 29-48. Wang, W. and Mulone, G., 2003. Threshold of disease transmission in a patchy environment, J. Math. Anal. Appl. 285, 321-335. Wang, W. and Zhao, X.-Q., 2004. An epidemic model in a patchy environment, Math. Biosci. 190, 97-112. World Health Organization (WHO), 2006a, Cumulative Number of Confirmed Human Cases of Avian Influenza A/(H5N1) Reported to WHO. Available at: http:// www. who. int/ csr/disease/avian influenza/country/ cases table 2006 04 19/en/index.html. Accessed April 19. World Health Organization (WHO), 2006b, Summary of WHO measures related to international travel. Available at: http:// www. who. int/csr/sars/ travelupdate/en/. Accessed April 4. World Organization for Animal Health (WOAH), 13 April 2006. UPDATE ON AVIAN INFLUENZA IN ANIMALS (TYPE H5). Available at: http://www.oie.int/downld/avian%20INFLUENZA/A AI-Asia.htm. Accessed April 19.
30
