The network level reproduction number for infectious diseases with both vertical and horizontal transmission

Mathematical Biosciences 243 (2013) 67–80

Contents lists available at SciVerse ScienceDirect

Mathematical Biosciences journal homepage: www.elsevier.com/locate/mbs

The network level reproduction number for infectious diseases with both vertical and horizontal transmission Ling Xue ⇑, Caterina Scoglio Department of Electrical & Computer Engineering, Kansas State University, KS 66506, United States

a r t i c l e

i n f o

Article history: Received 26 July 2012 Received in revised form 30 January 2013 Accepted 6 February 2013 Available online 28 February 2013 Keywords: Reproduction number Vertical and horizontal transmission Heterogeneous networks Rift Valley fever Multiple species

a b s t r a c t A wide range of infectious diseases are both vertically and horizontally transmitted. Such diseases are spatially transmitted via multiple species in heterogeneous environments, typically described by complex meta-population models. The reproduction number, R0 , is a critical metric predicting whether the disease can invade the meta-population system. This paper presents the reproduction number for a generic disease vertically and horizontally transmitted among multiple species in heterogeneous networks, where nodes are locations, and links reflect outgoing or incoming movement flows. The metapopulation model for vertically and horizontally transmitted diseases is gradually formulated from two species, twonode network models. We derived an explicit expression of R0 , which is the spectral radius of a matrix reduced in size with respect to the original next generation matrix. The reproduction number is shown to be a function of vertical and horizontal transmission parameters, and the lower bound is the reproduction number for horizontal transmission. As an application, the reproduction number and its bounds for the Rift Valley fever zoonosis, where livestock, mosquitoes, and humans are the involved species are derived. By computing the reproduction number for different scenarios through numerical simulations, we found the reproduction number is affected by livestock movement rates only when parameters are heterogeneous across nodes. To summarize, our study contributes the reproduction number for vertically and horizontally transmitted diseases in heterogeneous networks. This explicit expression is easily adaptable to specific infectious diseases, affording insights into disease evolution. Ó 2013 Elsevier Inc. All rights reserved.

1. Introduction Communicable diseases are readily transmitted from one region to another [1,2]. Population travel continues to influence the temporal and spatial spread of infectious diseases [1,3]. Observation of the introduction of infectious agents resulting in spatial spreading of effective infections in different locations at different times [3], revealed great economic losses, many animal and human cases, and deaths. Noteworthy examples include the fourteenth century plague in Europe [1,4] and the sixteenth century smallpox epidemic in the New World [1]. More recent epidemics, including HIV/AIDS and West Nile virus in North America [5], and SARS in Asia [6], show infections spreading over vast regions and even jumping continents [7]. Many communicable diseases are propagated by two distinct mechanisms: vertical and horizontal transmission [8]. Vertical transmission occurs when infection is passed from mother to a portion of offspring [8,9], often transmitted by insect eggs and/or plant seeds [10]. A variety of diseases are transmitted vertically and horizontally, including the human diseases rubella, Hepatitis ⇑ Corresponding author. Tel.: +1 7855324646. E-mail addresses: [email protected] (L. Xue), [email protected] (C. Scoglio). 0025-5564/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.mbs.2013.02.004

B, Chagas disease, and AIDS [10,11]. Vertical transmission is a proven factor in the size and persistence of Rift Valley fever (RVF) epidemic [12]. The prevalence of vertical transmission establishes it as a crucial biological mechanism [11], potentially affecting infectious spreading in elaborate ways [13]. Therefore, vertical transmission acts to maintain the spread of infection [13,14]. The logical complement of vertical transmission is horizontal transmission. For animal and human diseases [10], horizontal transmission is often through direct or indirect contact with infectious hosts or infectious vectors, such as biting insects [10]. Spatially structured models, such as meta-population models or multiple-patch models are widely used in epidemiological modeling to capture the effect of space [15]. Meta-population models describe systems containing spatially discrete sub-populations connected by the movement of individuals between a set of patches or nodes [16,17]. Modeling the dynamics of large metapopulations is complex, presenting challenges during analysis [18]. One approach considers the mobility of individuals between discrete regions [18], creating a directed network where nodes represent locations and links are movements between locale [18]. The importance of tracking mobility rates and movement patterns is highlighted in the foot-and-mouth outbreak of 2001 in the United Kingdom [7]. There, infected cattle were widely

68

L. Xue, C. Scoglio / Mathematical Biosciences 243 (2013) 67–80

distributed before the movement ban was announced [19], prompting the necessary development of a transportation network capturing the spatial spread of foot-and-mouth disease [7]. Numerical tools are widely used to obtain quantitative results and analytic tools are used to understand model behaviors [3]. The reproduction number, R0 , defined as the average number of new infected individuals produced by one infectious individual, in a population with only susceptibles [20], is arguably the most important quantity in communicable disease modeling [20]. Theoretically, R0 plays an important role in analyzing the dynamics of an epidemic [20]. It is a quantity commonly used to estimate the dynamics of emerging infectious diseases at the beginning of an outbreak, aiding in the design of control strategies for established infections [20]. The next generation method developed by [21, 22, Chapter 5] and popularized by van den Driessche and Watmough [23] is one of many methods applied to compute the reproduction number for compartmental models. This method manages matrix size by including only infected and asymptomatically infected states [24]. The next generation matrix relates the number of new cases in various generations and provides the basis of defining and computing the reproduction number [20]. The very little work on the reproduction number for metapopulations with vertical transmission we encountered included the modeling of horizontal and vertical transmission dynamics of a parasite with two differential equations [25]. In this special case, the reproduction number is the sum of the reproduction numbers for both types of transmission, and does not hold for a more complicated situation, such as in the model [26], where the next generation matrices for the two types of transmission are not both scalars. As far as we know, an insightful explicit expression of R0 for multiple species meta-population model with complex transmission has not yet been presented. This paper presents the computation of the reproduction number and its bounds for compartmental models considering diseases with complex transmission. We consider meta-populations consisting of discrete, well-mixed subpopulations. We assume that individuals move between different nodes and the disease can be transmitted within a node. An nnode compartmental model incorporates h species, of which g species transmit a disease both vertically and horizontally and other species only transmit horizontally. All sojourn times are taken to be exponentially distributed, and vertical transmission is restricted to the egg stage with exponential duration. Presented here is a general network-level model applicable when studying the temporal-spatial propagation of an infectious disease with multi-species, vertical and horizontal transmission, where the reproduction number is derived as a function of the two types of transmission parameters. Finally, the exact value and bounds of the reproduction number for the RVF metapopulation model are computed and factors affecting the reproduction number are analyzed. We found the upper bound depends on both horizontal and vertical transmission, while the lower bound is determined solely by horizontal transmission. The contribution of our work is summarized as follows: 1. An explicit expression of the reproduction number considering vertical and horizontal transmission in a general multi-species, meta-population model is derived. 2. This formula for the reproduction number is applied to an RVF meta-population model to compute R0 and its bounds. 3. Numerical simulations show that livestock movement rates only affect R0 for heterogeneous networks relative to disease parameters. Our work facilitates computation of the exact reproduction number in a meta-population model with complex disease transmission.

The paper is organized as follows. Section 2 describes the next generation matrix approach used to derive an explicit expression of the reproduction number, and presents the general meta-population model beginning with two species, two-node network models, as well as computing the reproduction number. In Section 3, we apply our R0 formula to the RVF meta-population model, computing R0 and its bounds. The effects of livestock movement, heterogeneities of parameters, and the size of a network on the reproduction number are also studied through simulations. Section 4 provides a summary and discussion of mathematical derivations and simulation results. 2. The reproduction number for diseases with both vertical and horizontal transmission One frequently used method computes the reproduction number as the spectral radius of the next generation matrix [22, Chapter 5],[27,20]. For the ease of computation, only the compartments corresponding to infected and asymptomatically infected compartments are considered [20]. First, the original nonlinear ODE system is decomposed into two column vectors F ¼ ðF i Þ and V=V i Þ, where th F i is the i row of F representing the rate at which new infections appear in compartment i, and V i is the ith row of V. Moreover, þ  Vi ¼ V i  V i , where V i represents the rate at which individuals transfer out of compartment i, and V þ i is the rate at which individuals transfer into compartment i [23]. Assume that the number of infected and asymptomatically infected compartments is m. The Jacobian matrices F denoting transmission, and V denoting transition [20] are defined as:

F¼

  @F i ðx0 Þ ; @xj

V¼

  @V i ðx0 Þ ; @xj

ð1Þ

where x0 represents the disease free equilibrium (DFE), and xj is the number or proportion of infected individuals in compartment j; j ¼ 1; 2;    ; m. The spectral radius of a matrix A is denoted by qðAÞ. The reproduction number, R0 , is defined as qðFV 1 Þ [21]. To understand entries of FV 1 , called the next generation matrix, consider the consequence of an infected individual introduced into compartment k in a population at DFE [23]. The ði; jÞ entry of F represents the rate at which new infected individuals in compartment i are produced by infected individuals in compartment j [23]. The ðj; kÞ entry of V 1 represents the average time that an infected individual stays in compartment j [23]. Hence, the ði; kÞ entry of FV 1 represents the expected number of new infections in compartment i resulting from the infected individual originally introduced into compartment k [23], where i; k ¼ 1; 2;    ; m. Note that matrix F is nonnegative and V is proved to be a nonsingular M-matrix [23]. Recall that an n  n matrix A is an M-matrix if it can be expressed in the form A ¼ sI  B, such that matrix B is non-negative, and s P qðBÞ [28]. Next, we illustrate computational procedures for finding R0 using the next generation matrix method for susceptibleexposed-infectious-recovered (SEIR) compartmental models, assuming a disease is transmittable within a species and between different species, and movement rates for all species are independent of disease status. Daily time steps are used in all models. 2.1. Models for two species in two nodes We present two applications of a simplified system for a disease involving two species in a two-node network with movement between the two nodes. In the first example, R0 is computed while assuming only horizontal transmission is taking place. In the second example, the first model is extended by introducing vertical transmission into one species. The reproduction number is then computed.

69

L. Xue, C. Scoglio / Mathematical Biosciences 243 (2013) 67–80

2.1.1. R0 for two species with only horizontal transmission Below, a compartmental model for an infectious disease incorporating four compartments ðJ ¼ S; E; I; RÞ, two species ðk ¼ 1; 2Þ, two nodes ði ¼ 1; 2Þ, and only horizontal transmission is presented. The differential equations representing the dynamic behavior are:

By (1), the Jacobian matrices for this model are:

dSki ¼ rki  b1ki Ski I1i =N1i  b2ki Ski I2i =N2i  dki Ski dt 2 2 X X þ xkji Skj  xkij Ski

where the symbol a represents the direct sum of matrices, i.e.,   A 0 AaB = for matrices A and B. The subscript of the zero blocks, 0 B 4  4, indicates the size of the block. Matrices A; Mk and X k are:

j¼1;j–i

ð2Þ

j¼1;j–i

ð3Þ

j¼1;j–i

2 2 X X dIki xkji Ikj  xkij Iki ¼ eki Eki  cki Iki  dki Iki þ dt j¼1;j–i j¼1;j–i 2 X

2 X

dRki xkji Rkj  xkij Rki : ¼ cki Iki  dki Rki þ dt j¼1;j–i j¼1;j–i

ð4Þ

ð5Þ

node i at DFE are denoted by J 0ki and N 0ki , respectively. To compute R0 using the next generation matrix method, we need to prove the existence and uniqueness of DFE. At DFE, S01i ¼ N 01i , and S02i ¼ N 02i , as E01i ¼ I01i ¼ R01i ¼ E02i ¼ I02i ¼ R02i ¼ 0. This is a special case of the model for Theorem 5 (see appendix), which determines the  T existence of a unique solution N 01i N 02i . The equations related to exposed and infected compartments are ordered:

E12

E21

E22

I11

I12

I21

I22 T ¼ F H  V H ; where

3 b211 S11 I21 =N21 þ b111 S11 I11 =N11 6 b212 S12 I22 =N22 þ b112 S12 I12 =N12 7 7 6 6 b S I =N þ b S I =N 7 6 121 21 11 11 21 7 221 21 21 7 6 6 b122 S22 I12 =N12 þ b222 S22 I22 =N22 7 7; FH ¼ 6 7 6 0 7 6 7 6 0 7 6 7 6 5 4 0 0 3 2 d11 E11 þ e11 E11 þ x112 E11  x121 E12 7 6 d12 E12 þ e12 E12 þ x121 E12  x112 E11 7 6 7 6 d21 E21 þ e21 E21 þ x212 E21  x221 E22 7 6 7 6 7 6 d E þ e E þ x E  x E 22 22 22 22 221 22 212 21 7 VH ¼ 6 6 e E þ d I þ c I þ x I  x I 7: 11 11 112 11 121 12 7 6 11 11 11 11 7 6 6 e12 E12 þ d12 I12 þ c12 I12 þ x121 I12  x112 I11 7 7 6 4 e21 E21 þ d21 I21 þ c21 I21 þ x212 I21  x221 I22 5 2

e22 E22 þ d22 I22 þ c22 I22 þ x221 I22  x212 I21

044

A

0

044

S0

2

7 5;

6 VH ¼ 4

2

b211 N110

21

S0

b112 N120

0

0

b221 N210

12

S0

21

S0

0

b122 N220

6 X 2 ¼4

0

2i¼1 ki

e



2k¼1 X k

3 7 5;

ð6Þ

3

0

7 7 7 7 0 7 S12 b212 N0 7 7 22 7 7; 7 7 0 7 7 7 7 7 5 0 S22 b222 N0 22

x112

d12 þ e12 þ x121

d21 þ e21 þ x212

x221

x212

d22 þ e22 þ x221

d11 þ c11 þ x112

x121

x112

d12 þ c12 þ x121

d21 þ c21 þ x212

x221

x212

d22 þ c22 þ x221

6 M2 ¼4

2



x121

2

6 X 1 ¼4

2k¼1

d11 þ e11 þ x112

6 M1 ¼4

2

2k¼1 M k

S0

0

b111 N110

6 11 6 6 6 6 6 0 6 6 A¼6 6 6 0 6 b S21 6 121 N011 6 6 6 4 0

3

12

The number of newborn individuals of species k in node i per day is denoted by r ki . The number of species k individuals in node i of compartment J is denoted by J ki , and the total number of species k individuals in node i is N ki ¼ Ski þ Eki þ Iki þ Rki . Total individuals of species k infected daily in node i by species 1 and species 2 are b1ki Ski I1i =N 1i and b2ki Ski I2i =N 2i , respectively. The number of deaths from each compartment J per day is dki J ki . After the incubation period, eki Eki individuals transfer to infected compartment daily. Following the infection period, cki Iki recover from the infection each day. Movement rates for species k individuals in compartment J P P in and out of node i are 2j¼1;j–i xkji J kj and 2j¼1;j–i xkij J ki , respectively. Species k quantity in compartment J and the total number in

d ½ E11 dt

6 FH ¼ 4

2

dEki ¼ b1ki Ski I1i =N1i þ b2ki Ski I2i =N2i  eki Eki  dki Eki dt 2 2 X X þ xkji Ekj  xkij Eki j¼1;j–i

2

3 7 5; 3 7 5;

ð7Þ

3 7 5; 3 7 5:

ð8Þ

Because the matrices M1 ; M2 ; X 1 , and X 2 are all invertible, we can readily check:

2 6 V 1 H ¼ 4

3

2k¼1 M 1 k

0

2k¼1 Z k

2k¼1 X 1 k

7 5;

1 2 where Z k ¼ X 1 k ði¼1 eki ÞM k . The spectral radius of the next gener-

ation matrix F H V 1 H is:

02 B6 qðF H V 1 H Þ ¼ q@4

044

A

0

044

32 76 54

2k¼1 M 1 k

0

2k¼1 Z k 2k¼1 X 1 k

31   7C 5A ¼ q Að2k¼1 Z k Þ :

Therefore,

   RH0 :¼ q F H V 1 ¼ q Að2k¼1 Z k Þ ; H

ð9Þ

where RH0 is the reproduction number for horizontal transmission. 2.1.2. R0 for two species with vertical transmission in one species We keep the model for species 2 (Eqs. (2)–(5)) with k ¼ 2, while extending the model for species 1 by incorporating vertical transmission. The model for species 1 is:

70

L. Xue, C. Scoglio / Mathematical Biosciences 243 (2013) 67–80

dP1i ¼r 1i  b1 q1i I1i  h1i P1i dt dQ 1i ¼b1i q1i I1i  h1i Q 1i dt dS1i ¼h1i P1i  b11i S1i I1i =N1i  b21i S1i I2i =N2i  d1i S1i dt 2 2 X X þ x1ji S1j  x1ij S1i j¼1;j–i

U¼

ð12Þ

ð13Þ

2 2 X X dR1i x1ji R1j  x1ij R1i : ¼c1i I1i  d1i R1i þ dt j¼1;j–i j¼1;j–i

ð15Þ

The number of eggs laid by species 1 per day is denoted as r1i , including b1i q1i I1i infected eggs, and r 1i  b1i q1i I1i uninfected eggs. After the development period, h1i P 1i eggs develop into susceptible adults, and h1i Q 1i eggs develop into infected adults daily. The interpretations of other terms are the same as corresponding terms described in Section 2.1.1. At DFE, Q 01i ¼ E01i ¼ I01i ¼ R01i ¼ E02i ¼ I02i ¼ R02i ¼ 0; S01i ¼ N 01i , and S02i ¼ N 02i . Since this is another special case of the model for Theo T rem 5, a unique solution N 01i N 02i exists. In our second model, the equations related to exposed and infected compartments are ordered:

2

E22

I11

I12

I21

I22 T

3

b11 q11 I11

7 6 b12 q12 I12 7 6 7 6 6 b211 S11 I21 =N21 þ b111 S11 I11 =N11 7 7 6 6 b S I =N þ b S I =N 7 12 7 6 212 12 22 22 112 12 12 7 6 6 b121 S21 I11 =N11 þ b221 S21 I21 =N21 7 7 F ¼6 6 b S I =N þ b S I =N 7; 22 7 6 122 22 12 12 222 22 22 7 6 7 6 0 7 6 7 6 0 7 6 7 6 5 4 0 3

e22 E22 þ d22 I22 þ c22 I22 þ x221 I22  x212 I21 By (1), the Jacobian matrices for this model are:

F¼

022

U 28

082

FH

" ;

V¼

2i¼1 h1i

028

W 82

VH

3

a h1 1i i¼1 2

!

1 V 1 H W a h1i

0 7 7 7; 7 5 V 1 H

2

0 Since M ðFV ÞM ¼ 0 F H V 1 H  " # 0 I  22 M ¼ W a2 h1 I88 , we have i¼1 1i



R0 ¼ q FV

1



¼q

3

!

i¼1

1

1

2

1 6 UV 1 UV 1 H W a h1i H 7 7 6 i¼1 7 6 1 ! FV ¼ 6 7: 7 6 2 4 F V 1 W a h1 F V 1 5 H H H H 1i

"

F H V 1 H

2

W

UV 1 H 2 1 Wðai¼1 h1 1i ÞUV H

a h1 1i i¼1

!

# ,

where

!

UV 1 H

:

ð16Þ

R0 is a function of vertical and horizontal transmission parameters.  2 1 Since F H V 1 UV 1 H and W ai¼1 h1i H are both nonnegative matrices, by Theorem 4 in appendix,

R0 P qðF H V 1 H Þ:

ð17Þ

2.2. R0 for multiple species in a general network The model presented in Section 2.1.2 is generalized to model diseases transmitted among all h species in node i ði ¼ 1; 2; . . . ; nÞ. Suppose a disease is transmitted by species k ðk ¼ 1; 2;    ; hÞ vertically and horizontally if 1 6 k 6 g and only horizontally otherwise. The dynamical behavior is given by the system with 4hn þ 2gn differential equations:

dPki ¼ ½r ki  bki qki Iki  hki P ki dðkÞ dt

ð18Þ

dQ ki ¼ ½bki qki Iki  hki Q ki dðkÞ dt

ð19Þ

n X

xkji Skj 

j¼1;j–i

h11 Q 11



6 6 V 1 ¼ 6 6 4

and the next generation matrix FV 1 are:

2

þ

7 6 h12 Q 12 7 6 7 6 7 6 d11 E11 þ e11 E11 þ x112 E11  x121 E12 7 6 7 6 d12 E12 þ e12 E12 þ x121 E12  x112 E11 7 6 7 6 7 6 d21 E21 þ e21 E21 þ x212 E21  x221 E22 7: V¼6 7 6 d E þ e E þ x E  x E 22 22 22 22 221 22 212 21 7 6 7 6 6 h11 Q 11  e11 E11 þ d11 I11 þ c11 I11 þ x112 I11  x121 I12 7 7 6 6 h Q  e E þ d I þ c I þ x I  x I 7 12 12 12 12 121 12 112 11 7 6 12 12 12 12 7 6 5 4 e21 E21 þ d21 I21 þ c21 I21 þ x212 I21  x221 I22



6 7 2 6 7 W ¼ 6  a h1i 7: 4 i¼1 5

h X dSki bmki Ski Imi =Nmi  dki Ski ¼ hki Pki dðkÞ þ rki ð1  dðkÞÞ  dt m¼1

0

2

2

1

i¼1

j¼1;j–i

E21

022 ;

i¼1

The matrix V

2 2 X X dI1i x1ji I1j  x1ij I1i ð14Þ ¼h1i Q 1i þ e1i E1i  c1i I1i  d1i I1i þ dt j¼1;j–i j¼1;j–i

E12

2

3

042

022

dE1i ¼b11i S1i I1i =N1i þ b21i S1i I2i =N2i  e1i E1i  d1i E1i dt 2 2 X X þ x1ji E1j  x1ij E1i

d ½Q Q 12 E11 dt 11 ¼ F  V; where

2

#

a b1i q1i

024

ð11Þ

j¼1;j–i

j¼1;j–i

"

ð10Þ

# :

Here F H and V H are the matrices in (6) and

n X

xkij Ski

ð20Þ

j¼1;j–i

h n X dEki X bmki Ski Imi =N mi  eki Eki  dki Eki þ xkji Ekj ¼ dt m¼1 j¼1;j–i



n X

xkij Eki

ð21Þ

j¼1;j–i n X dIki xkji Ikj ¼ hki Q ki dðkÞ þ eki Eki  cki Iki  dki Iki þ dt j¼1;j–i



n X

xkij Iki

ð22Þ

j¼1;j–i n n X X dRki xkji Rkj  xkij Rki : ¼ cki Iki  dki Rki þ dt j¼1;j–i j¼1;j–i

ð23Þ

The daily number of species k individuals infected by species m is bmki Ski Imi =Nmi . The daily numbers of species k individuals in com-

71

L. Xue, C. Scoglio / Mathematical Biosciences 243 (2013) 67–80

Pn partment J moving in and out of node i are j¼1;j–i xkji J kj and Pn x J , respectively. Other terms in the above equations have j¼1;j–i kij ki the same meanings as the corresponding ones in Section 2.1.1 (Eqs. (2)–(5)) and Section 2.1.2 (Eq. (10)) except dðkÞ defined below, which is used to differentiate the horizontally-transmitting species and the species exhibiting both types of transmission.

dðkÞ ¼

1 for 1 6 k 6 g; 0

for g þ 1 6 k 6 h:

To compute R0 using the next generation matrix method, we need to find matrices F and V, omitted here due to large size. In determining Jacobian matrices F and V, the infected variables are ordered by compartment, species, and node index, i.e.,

Q 11 ; Q 12 ; . . . ; Q 1n ; Q 21 ; Q 22 ; . . . Q 2n ; . . . ; Q g1 ; Q g2 ; . . . ; Q gn ; E11 ; E12 ; . . . ; E1n ; E21 ; E22 ; . . . ; E2n ; . . . ; Eh1 ; Eh2 ; . . . ; Ehn ; I11 ; I12 ; . . . ; I1n ; I21 ; I22 ; . . . ; I2n ; . . . ; Ih1 ; Ih2 ; . . . ; Ihn : At DFE, Q ki ¼ Eki ¼ Iki ¼ Rki ¼ 0, and Ski ¼ N ki . By Theorem 5 in  T appendix, a unique solution N0k1 N0k2    N 0kn exists. Since incorporating multiple species in multiple nodes leads to matrices F and V growing very large, the computation of R0 is simplified by decomposing the matrices into blocks, deriving block upper or lower triangular matrices as follows:

F¼



0gngn

U gn2hn

02hngn

FH

2



g

!

n

W 2hngn

VH

where

2 FH ¼



0hnhn

Ahnhn

0hnhn

0hnhn

6 6 VH ¼ 6 6 4

 ;

3 0hnhn 7 7 k¼1 7; ! 7 h n h 5 a Xk  a a eki h

a Mk

k¼1

i¼1

k¼1

2

" U¼

g

0gnhn

!

n

a a bki qki

# 0gnðhgÞn ;

i¼1

k¼1

3 ! 6 7 g n 6 7 7:  a a h W ¼6 ki 6 7 4 k¼1 i¼1 5 0ðhgÞngn 0hngn

The block matrix A in F H is written into an h  h block matrix A ¼ ðAkm Þ and its ðk; mÞ entry is an n  n diagonal matrix  n S0 Akm ¼ ai¼1 bmki N0ki . The matrices Mk and X k are: mi

3 fk1 xk21  xkn1 6 x n fk2  xkn2 7 k12 7 6 Mk ¼ 6 7; and X k ¼ M k þ aðcki  eki Þ; 4     5 i¼1 2

xk1n





fkn ð24Þ

Pn

where fki ¼ dki þ eki þ j¼1;j–i xkij . Since matrices Mk and X k are invertible, according to Theorem 7, V H and V are invertible. It is easy to check: 2 " V 1 H ¼

hk¼1 M 1 k

#

0

hk¼1 Z k hk¼1 X 1 k

g

n

!

3

6 a ah1 0gn2hn 7 ki 7 6 k¼1 i¼1 7 6 1 !! ; V ¼6 7; g 7 6 n 1 5 4 V 1 W a ah1 V ki H H k¼1

i¼1

ð25Þ

where Z k ¼ 2.1.2, R0 is:

X 1 k



ni¼1 ki

e



M1 k .

n

k¼1

i¼1

! UV 1 : H

ð26Þ

 > 1, we can Moreover, (17) still holds. If the lower bound q F H V 1 H conclude that a network may be invaded without computing the upper bound or the exact value of R0 . The term F H V 1 is related to horizontal transmission, and the  g H n  1 UV H is related to vertical transmission, term W ak¼1 ai¼1 h1 ki making R0 a function of vertical and horizontal transmission parameters. Generally, R0 depends on demographic, disease and movement factors, proving too complicated to compute or analyze [7]. The complexity of computing R0 using Eq. (26) depends on a specific model for a certain disease. For the general model, we can only provide the formula of R0 in Eq. (26) and its lower bound in Inequality (17). In the following section, Eq. (26) is applied to an RVF virus transmission meta-population model. Then, based on the assumptions for the RVF model, we compute R0 using Eq. (26) and further derive lower bound and upper bound, providing insights into the role of model parameters on R0 . 3. The application of proposed method to RVF meta-population model

3 0gn2hn 7 5;

6 a a hki V ¼ 4 k¼1 i¼1

;

g

1 R0 ¼ qðFV 1 Þ ¼ q F H V 1 H  W a a hki

!!

Similar to the derivation in Section

Rift Valley fever is an emerging mosquito-borne disease mainly affecting and colonizing domestic ruminants and humans [29,30]. Main vectors of RVF include Aedes and Culex mosquitoes [30]. Humans and ruminants are main hosts [30]. Aedes mosquitoes are believed to be initial source of RVF outbreaks [31], since RVF virus-carrying eggs can survive in drought area soil for many years, later breeding infected mosquitoes in flooded habitats [32,33]. Ruminants infected by mosquito bites [29] can transmit RVF virus to Aedes feeding on them as blood meals [30]. Culex mosquitoes also amplify RVF virus transmission by ingesting blood from infected ruminants [29]. Most humans acquire RVF virus infection when bitten by infected mosquitoes or during contact Table 1 with body fluid of infected ruminants [34]. Next, we derive R0 for an RVF meta-population model to study the role of parameters and networks on the reproduction number. 3.1. The network-based RVF meta-population model In this section, the general model in Eqs. (18)–(23) of Section 2.2 is applied to study the dynamics of RVF virus transmission with h ¼ 4; g ¼ 1. Aedes and Culex mosquito vectors are considered in the model, as are livestock and human hosts. The RVF model is less complex than the general model presented in Eqs. (18)–(23). Here, we assume only livestock can move in and out of nodes, and all mosquitoes do not recover. We consider disease-induced mortality for livestock and humans, and carrying capacity for mosquitoes and humans. Due to lack of transmission by humans or direct intra-species transmission, this RVF model contains fewer infection terms than those in the general model. See appendix for the full model (Eqs. (47)–(67)) and relative parameters (Table 2). The number of species k individuals ðk ¼ 1; 2; 3; 4Þ from node i ði ¼ 1; 2; . . . ; nÞ in compartment J is represented by J ki , where k ¼ 1 (resp. 2, 3, 4) represents Aedes mosquitoes (resp. livestock, Culex mosquitoes, and humans). The parameter r 2i is the number of livestock born daily in node i ði ¼ 1; 2; . . . ; nÞ. The daily numbers of new born Aedes mosquitoes, Culex mosquitoes, and humans are bki N ki . A node index is added at the end of the subscript of a parameter only when referring to a parameter for a specific node. For example, b12i represents the contact rate from Aedes mosquitoes (k ¼ 1) to livestock (k ¼ 2) in node i.

72

L. Xue, C. Scoglio / Mathematical Biosciences 243 (2013) 67–80

Table 1 Different scenarios for numerical simulations in four-node networks. Other parameters are kept the same and homogeneous across all nodes during all realizations. The superscripts i; j ¼ 1; 2; 3; 4 and i > j. No.

Parameter

Livestock movement rates

R0

1

b12i > b12j , b21i > b21j , b23i > b23j , b32i > b32j

2

d2i > d2j

increases decreases decreases

3

c2i > c2j

4

l2i > l2j

x2ji increases x2ij increases x2ji increases x2ij increases x2ji increases x2ij increases x2ji increases x2ij increases

increases decreases increases decreases increases

Table 2 Parameters in the model omitting the node index. Parameter

Description

Range value or

Units

Source

b12 b21 b23 b32 b14 b24 b34 1=c2 1=c4 1=d1 1=d2 1=d3 1=d4 b1 b3 b4 1=1 1=2 1=3 1=4

contact rate: Aedes to livestock contact rate: livestock to Aedes contact rate: livestock to Culex contact rate: Culex to livestock contact rate: Aedes to humans contact rate: livestock to humans contact rate: Culex to humans recover rate in livestock recover period in humans longevity of Aedes mosquitoes longevity of livestock longevity of Culex mosquitoes longevity of humans egg laying rate of Aedes mosquitoes egg laying rate of Culex mosquitoes birth rate of humans incubation period in Aedes mosquitoes incubation period in livestock incubation period in Culex mosquitoes incubation period in humans mortality rate in livestock transovarial transmission rate in Aedes mosquitoes development period of Aedes mosquitoes development period of Culex mosquitoes carrying capacity of Aedes mosquitoes carrying capacity of Culex mosquitoes carrying capacity of humans livestock recruitment rate livestock movement rate from node i to node j

ð0:0021; 0:2762Þ ð0:0021; 0:2429Þ ð0:0000; 0:3200Þ ð0:0000; 0:096Þ

1=day 1=day 1=day 1=day 1=day 1=day 1=day days days days days days days 1=day 1=day 1=day days days days days 1=day 1=day days days

[42–48] [42–46,49] [43–46,49,50] [43–46,50]

1=day 1=day

[52]

l2 q1 1=h1 1=h3 K1 K3 K4 r2i

x2ij

ð2; 5Þ ð4; 7Þ ð3; 60Þ ð360; 3600Þ ð3; 60Þ

ð4; 8Þ ð2; 6Þ ð4; 8Þ ð2; 6Þ ð0:025; 0:1Þ ð0; 0:1Þ ð5; 15Þ ð5; 15Þ 10000 10000 100000 1  1 0; n

[51] [52] [53,54,46] [55] [53,54,46] [53,54,46] [53,54,46] [47] [56] [47] [52] [51,56] [57] [46] [46]

2

3.2. The computation of R0 for RVF

2

The explicit expression of R0 in Eq. (26) is applied to the RVF meta-population model. The above assumptions allow us to obtain the lower and upper bounds of R0 .

6 FH ¼ 4

04n4n

A4n4n

04n4n

04n4n

3 7 5;

6 6 6 VH ¼ 6 6 6 4

3

4

aM k k¼1

ð4k¼1 ðni¼1 eki ÞÞ4n4n

04n4n 7 7 7 7: 7 7 4 5 aX k k¼1

3.2.1. Explicit expression of R0 for RVF First, we check if a unique solution N 0ki exists. At DFE,

ð27Þ 2

E0ki ¼ I0ki ¼ R0ki ¼ 0. By computation, S0ki ¼ N 0ki ¼ bkid K k for k ¼ 1; 3; 4, ki

where K k is the carrying capacity of species k. This is a special case of the model for Theorem 5, which generates a unique nonnegative solution for the total number of livestock in node i at DFE denoted  T by: N 021 N 022    N 02n . By (1), the Jacobian matrices for the RVF model are:

2 6 F¼4

0nn 08nn

U n8n FH

3

2

7 5;

6 V ¼4

ni¼1 h1i W 8nn

0n8n VH

3 7 5:

Each component of the R0 formula is computed as follows:

 U ¼ 0n4n

 ni¼1 ðb1i q1i Þ 0n3n ;

04nn

3

7 6 7 6 7 6 W ¼ 6 ðni¼1 h1i Þ 7: 7 6 5 4

ð28Þ

03nn 2

0

6 6 6 A21 6 A¼6 6 6 0 4 A41

A12

0

0

A23

A32

0

A42

A43

0

3

7 7 07 7 7; 7 07 5 0

ð29Þ

73

L. Xue, C. Scoglio / Mathematical Biosciences 243 (2013) 67–80

A12 ¼ ni¼1 b21i A32 ¼ A43 ¼

S01i

N02i S0 ni¼1 b23i 3i0 N2i S0 ni¼1 b34i 4i0 N3i

A21 ¼ ni¼1 b12i

;

; A41 ¼

ni¼1 b14i

S02i N01i

S04i N01i

;

A23 ¼ ni¼1 b32i

;

A42 ¼

ni¼1 b24i

S02i N03i

S04i N02i

where

;

;

:

þ e1i ; X 1 ¼ M 1  ni¼1 e1i ;  0 d N M 3 ¼ ni¼1 3iK 3 3i þ e3i ; X 3 ¼ M 3  ni¼1 e3i ;  0 d N M 4 ¼ ni¼1 4iK 4 4i þ e4i ; X 4 ¼ M 4  ni¼1 e4i ; 2



d1i N 01i K1

f21

6 6 6 x212 6 M2 ¼ 6 6 6  6 4 x21n

x221 f22  x22n

   x2n1

b1i ¼b1 ;

3

b23i ¼ b23 ;

b32i ¼ b32 :

ð33Þ

ð34Þ

i

7 7    x2n2 7 7 7; 7   7 7 5

where



Theorem 3. Under the condition of Theorem 2, R0 can be estimated by the following inequality:

v¼

f2n

Theorem 1. Consider the model presented in Section 3.1 (Eqs. (47)–(67)), we obtain

q

e1i ¼ e1 ; e2i ¼ e2 ; e3i ¼ e3

b21i ¼ b21 ;

Then

3.2.2. Deriving lower bound and upper bound for R0 Bounds of R0 are derived in many articles, among which are some following examples. Gao and Ruan present bounds of R0 for an SIS patch model [36] investigating effects of media coverage and human movement on the spread of infectious diseases, as well as a malaria model [37]. Hsieh et al. [38] derive bounds of R0 , describing the relationship between the reproduction numbers for the isolated ith patch and for the system. Salmani and Driessche [1] derive bounds for an SEIRS patch model. Arino [35] presents bounds of R0 for patch models considering multiple species. The reproduction number for an averaging process of mixed individuals or groups is estimated to be smaller than or equal to the reproduction number before mixing [39]. We derive bounds of R0 for RVF meta-population model in this section. In the following, we shall state main results and prove them in appendix.



b3i ¼ b3 ;

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi  ffi 1 1 1 vq X 1 M 6 v q X M 6 R þ maxðq1i Þ: 0 2 2 2 2

The reproduction number, R0 can be computed by plugging the above terms into Eq. (26). Typically R0 for a meta-population model is complicated [35]. Deriving some bounds on the value of R0 can be helpful [35]. In the following section, we derive lower and upper bounds for R0 .

F H V 1 H

Corollary 1. Suppose for all i, birth and incubation rates in mosquitoes and livestock, contact rates between livestock and mosquitoes are homogeneous for different nodes, i.e.,

b12i ¼b12 ;

X 2 ¼ M2 þ ni¼1 ðc2i þ l2i  e2i Þ:



ð32Þ

The difference between the lower bound and the upper bound in a network with heterogeneous corresponding parameters across nodes is larger than that in Inequality (30).

1 The matrices V 1 are in Eq. (25) with g = 1 and h = 4, H and V respectively. Below, matrices M k and X k relate to Aedes mosquitoes, livestock, Culex mosquitoes, and humans with k ¼ 1; 2; 3; 4, respectively.

M 1 ¼ ni¼1

e1i e2 b12i b21i e2 e3i b32i b23i þ : b1i ðb1i þ e1i Þ b3i ðb3i þ e3i Þ

vi ¼

 6 R0 6 q F H V 1 þ maxðq1i Þ: H

ð30Þ

i

e1 e2 b12 b21 e2 e3 b32 b23 þ : b1 ðb1 þ e1 Þ b3 ðb3 þ e3 Þ

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u minðvi Þ u i t 6 R0 maxðd2i þ e2 Þmaxðd2i þ c2i þ l2i Þ i

Theorem 2. For the model in Section 3.1 (Eqs. (47)–(67)), assume

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 minðvi ÞqðX 1 2 M 2 Þ 6 R0

i

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 6 maxðvi ÞqðX 1 2 M 2 Þ þ maxðq1i Þ; i

i

ð31Þ

ð36Þ

i

If the differences between mini ðvi Þand maxi ðvi Þ; mini ðd2i þ e2 Þ and maxi ðd2i þ e2 Þ; mini ðd2i þ c2i þ l2i Þ and maxi ðd2i þ c2i þ l2i Þ are large, then the difference between the lower bound and the upper bound may be large. However, the scalar lower bound and upper bound are easily computed. Moreover, if the lower bound is greater than 1, we can conclude that the network may be invaded without computing R0 or its upper bound. Corollary 2. Based on the condition of Corollary 1, we further assume that for all i, the death rate, mortality rate, and recovery rate in livestock, and transovarial transmission rate in Aedes mosquitoes are homogeneous for all nodes, i.e.,

d2i ¼ d2 ;

l2i ¼ l2 ; c2i ¼ c2 ; q1i ¼ q1 :

ð37Þ

Then

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

v

ðd2 þ e2 Þðd2 þ c2 þ l2 Þ

6 R0 rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 6

v

ðd2 þ e2 Þðd2 þ c2 þ l2 Þ

þ q1 :

ð38Þ

In this case, the lower and upper bounds of R0 correspond to the bounds for homogeneous populations presented in [26] and are tight [26]. Clearly, R0 for horizontal transmission,

RH0 ¼

i

i

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u maxðvi Þ u i þ maxðq1i Þ: 6t i minðd2i þ e2 Þminðd2i þ c2i þ l2i Þ

The difference between the lower and upper bounds is maxi ðq1i Þ with lower bound qðF H V 1 H Þ computed by Eq. (43).

e2i ¼ e2 for all i, then

ð35Þ

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   e2 e1 b12 b21 e3 b32 b23 ; þ ðb2 þ e2 Þðb2 þ c2 þ l2 Þ b1 ðb1 þ e1 Þ b3 ðb3 þ e3 Þ

ð39Þ

does not depend on livestock movement rates. Only bounds for R0 can theoretically be obtained. Based on numerical simulation

74

L. Xue, C. Scoglio / Mathematical Biosciences 243 (2013) 67–80 5

5

lower bound

upper bound

y=x

y=x

4

3

3

R0

R L0

U

4

2

2

1

1

0

0

1

2

3

4

0

5

0

1

2

R0

3

4

5

R0

Fig. 1. The reproduction number and its lower and upper bounds computed using Theorem 1 for one hundred simulation runs in one hundred-node heterogeneous networks.

upper bound y=x

5

5

lower bound y=x 4

3

3

R0

U

R L0

4

2

2

1

1

0

0

1

2

3

4

5

R0

0

0

1

2

R0

3

4

5

Fig. 2. The reproduction number and its lower and upper bounds computed using Theorem 2 for one hundred simulation runs in one hundred-node heterogeneous networks.

results, we conjecture that, given the conditions for Corollary 2, R0 does not depend on livestock movement rates.

3.2.3. Tightness of bounds for R0 A one hundred-node network with heterogeneous corresponding parameters among nodes is built to study the tightness of bounds. We uniformly distribute disease parameters for each node during one hundred runs within their respective ranges, given in Table 2. Then, R0 is numerically computed according to Eq. (26). Lower and upper bounds of R0 are computed according to Inequality (30) in Theorem 1. The reproduction number for horizontal transmission is computed according to Eq. (43). The lower bound of R0 (denoted by RL0 ) versus R0 in each run is shown in Fig. 1(a), and the upper bound of R0 (denoted by RU0 ) versus R0 in each run is shown in Fig. 1(b). In each run, the upper bound is slightly greater and the lower bound is slightly smaller than R0 . With the same network and the same set of parameters, the lower and upper bounds of R0 are computed using Inequality (31). The lower bound versus exact R0 is shown in Fig. 2(a), and the upper bound versus exact R0 is shown in Fig. 2(b). The bounds obtained by Inequality (31) in Theorem 2 are less tight than those obtained by Inequality (30) in Theorem 1, as qðF H V 1 H Þ is estimated by computing the spectral radius of a smaller size matrix. The bounds obtained by Inequality (36) in Theorem 3 can be even looser because 1 qðX 1 2 M 2 Þ is simply estimated by scalars.

The above bounds are for heterogeneous networks. The bounds in Corollary 2 (see Inequality (38)) apply to homogeneous networks, where the difference between the lower bound and the upper bound is the largest transovarial transmission rate of Aedes mosquitoes across nodes.

3.3. Assessing the role of parameters on R0 As an example, a two-node network demonstrates how bounds of R0 alter with livestock movement rates, if parameters d2i ; c2i , and l2i are heterogeneous, i.e., at least one of inequalities d2i – d2j ; c2i – c2j ; l2i – l2j holds for different i and j. In this example, M2 corresponds to the one in Eq. (7) and X 2 ¼ M2 þ 2i¼1 ðc2i þ l2i  e2i Þ. Since X 2 ; M 2 are both diagonal dom1 inant matrices, by Theorem 7, M 1 2 and X 2 are both nonnegative matrices. 1 According to Proposition 4.3 in [37], qðX 1 2 M 2 Þ is decreasing in x212 if

x212 ða2  a1 Þ > ða1 c1  a2 c2 Þ  ða2  a1 Þx221 and increasing otherwise, where a1 ¼ d21 þ e21 ; a2 ¼ d22 þ e22 ; c1 ¼ d21 þ c21 þ l21 and c2 ¼ d22 þ c22 þ l22 . In the case that 1 a1 ¼ a2 ; qðX 1 2 M 2 Þ is decreasing in x212 if c2 > c1 and increasing 2 c2  x221 is a critical point of otherwise. If a1 – a2 ; x212 :¼ a1ac12 a a1

L. Xue, C. Scoglio / Mathematical Biosciences 243 (2013) 67–80

75

Fig. 3. The reproduction number for four-node networks with different contact rates during one hundred runs.

Fig. 4. The reproduction number for four-node networks with different livestock death rates during one hundred runs.

Fig. 5. The reproduction number for four-node networks with different livestock recovery rates during one hundred runs.

1 1 1 qðX 1 2 M 2 Þ. Moreover, qðX 2 M 2 Þ reaches the maximum value at x212 if a2 > a1 and the minimum value at x212 otherwise.

To evaluate the impact of networks with corresponding homogeneous parameters across all nodes on the value of R0 computed using Eq. (26), we construct three networks with three, four, and

one hundred nodes, respectively. Simulation runs with varying livestock movement rates, and parameters in (33) and (37) held constant and homogeneous across nodes showed R0 is not affected by livestock movement rates during one hundred runs per network. Moreover, the values and bounds of R0 obtained through

76

L. Xue, C. Scoglio / Mathematical Biosciences 243 (2013) 67–80

Fig. 6. The reproduction number for four-node networks with different livestock mortality rates during one hundred runs.

numerical simulations are the same for networks with three, four, and one hundred nodes. Through extensive numerical simulations, we have observed that R0 does not depend on livestock movement rates or the number of nodes in a network when (33) and (37) hold. We run scenarios (see Table 1) one hundred times for each fournode network to study the impact of livestock movement rates on R0 . During one hundred realizations for each scenario, we increase livestock movement rates while keeping remaining parameters constant and homogeneous across all nodes. In Scenario 1, we set contact rates b12 ; b21 ; b23 , and b32 for node i larger than respective parameters for node j (i > j; i; j ¼ 1; 2; 3; 4). During each run, R0 increases while increasing livestock movement rates from node j to node i; x2ji , and decreases while increasing livestock movement rates from node i to node j; x2ij (see Fig. 3(a) and (b), respectively). In Scenario 2, under setting d2i > d2j ; R0 decreases when x2ji increases, and increases when x2ij increases (see Fig. 4(a) and (b), respectively). With livestock recovery rates c2i > c2j in Scenario 3; R0 decreases when x2ji increases, and increases when x2ij increases (see Fig. 5(a) and (b), respectively). Similarly, when livestock mortality rates l2i > l2j in Scenario 4; R0 decreases when x2ji increases, and increases with larger x2ij (see Fig. 6(a) and (b), respectively). Tuning the parameters in above scenarios yields R0 from below 1 to above 1. As a consequence, livestock movement rates are important in either leading to an epidemic outbreak or epidemic burnout. 4. Results and discussions We propose an explicit expression of R0 , which is formulated as a function of vertical and horizontal transmission parameters shown in Eq. (26). This formula facilitates computing R0 for many diseases that involve both vertical and horizontal transmission by replacing the spectral radius of the original next generation matrix with that of a smaller size matrix. The lower bound of R0 equals the reproduction number for horizontal transmission. We applied Eq. (26) to the RVF model, deriving R0 and its lower and upper bounds. We compared the tightness of different bounds, and analyzed the role of livestock movement rates and disease parameters on R0 through numerical simulations. The reproduction number for RVF meta-population model relates to the reproduction number for horizontal transmission, involving Aedes-livestock interaction and Culex-livestock interaction, and vertical transmission parameters. Different bounds of R0 for heterogeneous networks are given by Theorems 1–3 with decreasing tightness and increasing easiness. For homogeneous networks, the reproduction number for horizontal transmission in Eq. (39) and bounds of R0 given by Corollary 2 are proved inde-

pendent of livestock movement rates, and equal to corresponding terms for homogeneous populations presented in [26]. The lower bound is the reproduction number for horizontal transmission and upper bound is the sum of the reproduction number for horizontal transmission and the largest transovarial transmission rate of Aedes mosquitoes among nodes. Typically networks in the real world are heterogeneous. Rates of livestock death, incubation, mortality, recovery, and contact with mosquitoes can vary in different nodes due to climate, public health facilities, environment, and/or type of nodes (e.g., death rates of livestock in feedlots are higher than those in livestock premises). Variations in weather may affect values of some mosquito parameters, e.g., rainfall affects mosquito birth rates, and temperature affects mosquito incubation rates. Even if weather conditions are homogeneous across all nodes, different genera and/or species of mosquitoes can exhibit different rates of incubation, contact, death, birth, and/or birth. Numerical simulations show livestock movement rates between different nodes only affect R0 when the network is spatially heterogeneous regarding parameters. Changing livestock movement rates on heterogeneous networks results in R0 varying between values below and above the critical value 1. When other parameters remain homogeneous and constant, increasing livestock movement rates from nodes with smaller contact rates to those with larger contact rates increases R0 . If livestock movement rates are increased from nodes with smaller livestock death rates (or recovery rates, or mortality rates) to nodes with larger livestock death rates (or recovery rates, or mortality rates), R0 decreases. This observation helps us better envision effective mitigation strategies executing movement bans between some nodes and in some directions. Whatever heterogeneity exists between nodes, our same mathematical model in Eq. (18) through (23), and the explicit expression of R0 in (26), are applicable. Our formula for R0 presented in this paper can be used for numerous diseases models aside from RVF. Our work on RVF contributes computing R0 accurately by taking into account vertical transmission, which is important but ignored by modelers. We simplified the derivation of R0 by computing the spectral radius of a smaller size matrix than the original next generation matrix. Bounds of R0 facilitate estimating R0 of RVF metapopulation model. The simulation results on livestock movement rates and parameters are helpful in developing efficient mitigation strategies for RVF. Acknowledgement This material is based upon work supported by the U. S. Department of Homeland Security under Grant Award Number

77

L. Xue, C. Scoglio / Mathematical Biosciences 243 (2013) 67–80

2010-ST061-AG0001 and a grant through the Kansas Biosciences Authority. The views and conclusions contained in this publication are those of the authors and should not be interpreted as necessarily representing the official policies, either explicit or implicit, of the U. S. Department of Homeland Security and Kansas Biosciences Authority. We are grateful to the effort made by the anonymous reviewers. We would like to give thanks to Alice Trussell and Andrea Engelken for help on English proofreading.

Proof of Theorem 1. The left inequality is the same as (17). We now show that the right inequality holds. By (28) and (25),

W

 Y¼

ni¼1 h1 1i

X 1 1



  1 04n4n UV H ¼ Y

04n4n

of

ðni¼1 ðb1i q1i ÞÞX 1 1 .

1 are diagonal entries of Wðni¼1 h1 1i ÞUV H  n 1  1 Hence, W i¼1 h1i UV H ¼ PDP 1 for some P.

Here

04n4n

04n4n

04n4n

Q

 ;





 L¼

H4n4n J 4n4n 0n3n I3n3n

Q¼

03nn

0n3n

ðni¼1 ðb1i q1i ÞÞX 1 1

:

0 L4n4n

 ;

where H ¼



Inn ; 03nn

 J ¼



 0n3n ; I3n3n

ðni¼1 ðb1i q1i ÞÞX 1 1 03nn

1 ðni¼1 ðb1i q1i 1i ÞÞX 1 1 M1

e

03nn





0n3n : 03n3n



1 1 n qðFV 1 Þ ¼ qðF H V 1 H  W i¼1 h1i UV H Þ

¼ qðP 1 F H V 1 H P þ DÞ: We claim that P

L1

1

F H V 1 H P

is a nonnegative matrix. By calculation,

 0 ; L1  I3n3n : 0n3n

H1 1 L JH1

 03nn ¼ Inn

ð40Þ "

H1 ¼

ðni¼1 b1i1q1i ÞX 1

0n3n

03nn

I3n3n

# ;

It is clear that H1 ; L1 , and L1 JH1 are all nonnegative matrices. Hence, P 1 is a nonnegative matrix. We now show that F H V 1 H P is a nonnegative matrix.

F H V 1 H P ¼

d N0 ni¼1 1iK 1 1i

ð41Þ

and N 01i ¼ b1id K 1 , we further have 1i

    1 qðDÞ ¼ qðW ni¼1 h1 ni¼1 ðb1i q1i Þ X 1 1i UV H Þ ¼ q 1   ¼ q ni¼1 q1i 6 maxðq1i Þ:

"   A 4k¼1 Z k H þ Að4k¼1 X 1 k ÞJ 0

 # L A 4k¼1 X 1 k

1 1 n n A21 X 1 1 ði¼1 e1i ÞM 1 ði¼1 ðb1i q1i ÞÞX 1 1 1 n þ A21 X 1 1 ði¼1 ðb1i q1i e1i ÞÞX 1 M 1 ¼ 0:

i

This finishes the proof. Proof of Theorem 2. By Eqs. (25) and (27),

" F H V 1 H

¼

Að4k¼1 Z k Þ Að4k¼1 X 1 k Þ 0

#

0

:

Then

ð42Þ

3 2 3 0 B1 0 0 0 0 0 A12 Z 2 7 6 6   6 A21 Z 1 0 A23 Z 3 0 7 6 B2 0 B3 0 7 7 A 4k¼1 Z k ¼ 6 7 ¼: 6 7: 4 0 A32 Z 2 0 0 5 4 0 B4 0 0 5 2

A41 Z 1 A42 Z 2 A43 Z 3 0

;

0

B5 B6 B7 0

To compute the eigenvalues of Að4k¼1  Z k Þ, we first calculate the characteristic polynomial of A 4k¼1 Z k as follows.

  kIn        B2 jkI4n  A 4k¼1 Z k j ¼   0     B 5

 0      kIn B1 0       kIn B3 0     ¼ kn  B2 kIn B3        B4 kIn 0      0 B kI  4 n   B6 B7 kIn 0

B1

2 2 32 3 3  In kB1 0  0 B1 þ k2 B1 kB1 B3  kIn B1 0  2 2 2    6 6 76 7 7 6 6 76 7 7 n 6 n  6 6 7 ¼ k 6 0 I n 0 7 kIn B3 7 76 B2 kIn B3 7 ¼ k 6 B2 7 4 4 54 5 5       0 0 I  0 0 B4 kIn  B4 kIn n   2 32 3   B1 þ k2 B1 kB1 B3   B1 þ k2 B1 kB1 B3 In kB1 2 2 2 2 4     n    4 5 4 5 ¼ k jB2 j  ¼ k jB2 j     0 In  B4 kIn B4 kIn n

   B1 þ k2 B1 kðB1 B1  k2 B1 B1 þ B1 B3 Þ  2 4 2 4 2    ¼ k jB2 j    B4 0 n

 L is a nonnegative matrix and where A 4k¼1 X 1 k  n  1 1 Z k ¼ X k i¼1 eki Mk . Furthermore, the only possible negative en   tries of A 4k¼1 Z k H þ A 4k¼1 X 1 J are in its ð2; 1Þ and ð4; 1Þ k blocks. But the block in ð2; 1Þ-entry is

1 qðF H V 1 Þ 6 qðF H V 1 H Þ 6 R0 ¼ qðFV H Þ þ maxðq1i Þ:

By Eq. (29),

 n 1  1 1 1 Since F H V 1 F H V 1 , we have H  W i¼1 h1i UV H ¼ PðP H P þ DÞP

 P 1 ¼

1 qðFV 1 Þ 6 qðP 1 F H V 1 H PÞ þ qðDÞ ¼ qðF H V H Þ þ qðDÞ:

  4  RH0 ¼ qðF H V 1 H Þ ¼ q A k¼1 Z k :



03n3n

From linear algebra, each column of P can be chosen as an 1 eigenvector of Wðni¼1 h1 1i ÞUV H . By direct calculation,

P¼

Hence, F H V 1 H P is a nonnegative matrix. This proves the claim. By Theorem 2 in [40], we have

Therefore,

; where

Z

  n    i¼1 ðb1i q1i Þ X 1 ni¼1 ðb1i q1i e1i Þ M 1 0n3n 0n3n 1 1 ; Z¼ : 03nn 03n3n 03nn 03n3n

eigenvalues

D¼

1 1 n þ A41 X 1 1 ði¼1 ðb1i q1i e1i ÞÞX 1 M 1 ¼ 0:

i



Note that X 1 and X 1 1 are diagonal matrices. Moreover, the nonzero



1 1 n n A41 X 1 1 ði¼1 e1i ÞM 1 ði¼1 ðb1i q1i ÞÞX 1

Since X 1 ¼

Appendix



By assumption, X 1 and M 1 are both diagonal matrices. The last 1 1 1 equality follows X 1 1 M 1 ¼ M 1 X 1 . Similarly, the block in ð4; 1Þ-entry is

   kðB1 B1  k2 B1 B4 þ B1 B3 Þ B1  k2 B1  4 2 2 2    ¼ kn jB2 j    0 B4   n 2 1 1 1 1   ¼ k jB2 jjB4 j kðB1 B  k B B þ B B3 Þ 4

2

4

2

  1 1 1  ¼ k jB2 jjB4 j k2 B1 2 B4  ðB1 B 4 þ B 2 B 3 Þ 2n

  1  2 1  ¼ k2n jB2 jjB4 jjB1 2 jjB 4 j k I n  ðB4 B 2 B1 B 4 þ B 4 B3 Þ   : ¼ k2n  k2 In  ðB4 B2 B1 B1 4 þ B4 B 3 Þ

78

L. Xue, C. Scoglio / Mathematical Biosciences 243 (2013) 67–80

Matrix Að4k¼1 Z k Þ has 2n zero eigenvalues. The spectral radius of Að4k¼1 Z k Þ is the square root of the spectral radius of

for any matrix norm k  k. If A; B are both non-negative, then

B4 B2 B1 B1 4

of matrix norm, 0 6 kAk k 6 kðA þ BÞk k. Thus,

þ B4 B3 . By Eq. (42), we obtain

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qðF H V 1 qðB4 B2 B1 B1 H Þ ¼ 4 þ B4 B3 Þ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ qðB4 ðB2 B1 þ B3 B4 ÞB1 qðB2 B1 þ B3 B4 Þ: 4 Þ ¼

1

k!1

ð43Þ

 1 1 1 n e2 e1i b12i b21i X 1 M 1 ; B2 B1 ¼ ðni¼1 e1i e2 ÞA21 X 1 M A X M ¼  12 1 1 2 2 i¼1 b1i ðb1 þ e1i Þ 2 2 1 1 1 n e2 e3i b32i b23i B3 B4 ¼ ðni¼1 e2 e3i ÞA23 X 1 ÞX 1 M 1 : 3 M 3 A32 X 2 M 2 ¼ ði¼1 b3i ðb3i þ e3i Þ 2 2

1 ðX 1 2 M2 Þ

minðvi Þq i

The theorem follows from the Gelfand’s formula. Theorem 5. For the model presented in Section 2.2 (Eqs. (18) through (23)), a unique nonnegative solution for the total number of species k individuals in node i at DFE exists.

W ½ N k1

1 6 qðB2 B1 þ B3 B4 Þ 6 maxðvi ÞqðX 1 2 M 2 Þ:

Nk2

2

i

According to Theorem 1,

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1 minðvi ÞqðX 1 maxðvi ÞqðX 1 2 M 2 Þ 6 R0 6 2 M 2 Þ þ maxðq1i Þ: i

i

This finishes the proof. Proof of Corollary 1. By the assumptions in (33), we have mini ðvi Þ ¼ v ¼ maxi ðvi Þ. Corollary follows from Theorem 2. Proof of Theorem 3. According to Theorem 6,

1 1 6 qðX 1 2 M2 Þ maxðd2i þ e2i Þmaxðd2i þ c2i þ l2i Þ i

i

1 : minðd2i þ e2i Þminðd2i þ c2i þ l2i Þ i

r k2

   r kn T ;

ð45Þ 3

n X

xk1j 6 dk1 þ 6 j¼2 6 6 6 6 xk12 W¼6 6 6  6 6 6 4 xk1n

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1 1 maxðvi ÞqðX 1 minðvi ÞqðX 1 2 M 2 Þ 6 qðF H V H Þ 6 2 M 2 Þ:

i

T

   Nkn  ¼ ½ r k1

where

i

Therefore,

6

k!1

Proof. To solve the total number of species k individuals in each node at DFE, we need to solve the following system of equations.

vi in (32), we have

i

1

0 6 lim kAk kk 6 lim kðA þ BÞk kk :

Recall that A21 ; A12 ; X 1 ; M1 ; A23 ; A32 ; M 3 ; X 3 are all diagonal matrices. By the assumption that e2i ¼ e2 , for all i, we obtain

By the definition of

A 6 A þ B. Hence, 0 6 Ak 6 ðA þ BÞk for any k 2 N. By the property

ð44Þ

i

xk21 n X

dk2 þ



xk2j

j¼1;j–2

 xk2n

xkn1

7 7 7 7 7 7  xkn2 7: 7 7   7 7 n1 X 7    dkn þ xknj 5 j¼1

½ Nk1

Nk2

T    Nkn 

The variable vector is to be solved. We note that W is a diagonal dominant matrix of its column entries [41], i.e., P W ii P ni¼1;i–j W ij , for all i, where W ij is the ði; jÞ entry of W. By Theorem 1 in page 654 of [41], W is invertible. Moreover, by Theorem 7 in appendix, W 1 is nonnegative. Thus, there exists a unique nonnegative solution for the system of equations (45). The unique nonnegative solution is

 T T ½ Nk1 Nk2   Nkn  ¼ N0k1 N0k2   N0kn ¼ W 1 ½ r k1 r k2   r kn T :

By Theorem 2, we have

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u u minðvi Þ u i t 6 R0 maxðd2i þ e2i Þmaxðd2i þ c2i þ l2i Þ i

Theorem 6. Let Ak ðk ¼ 1; 2;    ; mÞ be an n  n diagonal dominant matrix with A1 P 0. Denote the ði; jÞentry of Ak by akij . Let P k P aLk ¼ minj ð i akij Þ > 0; aHk ¼ maxj ð i akij Þ, where j ¼ 1; 2;    ; n, then

i

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u maxðvi Þ u u i þ maxðq1i Þ: 6t i minðd2i þ e2i Þmind2i þ c2i þ l2i Þ i

ð

This finishes the proof. Proof of Corollary 2. By the conditions in (33) and (37), we have

minðvi Þ i

maxðd2i þ e2i Þmaxðd2i þ c2i þ l2i Þ i

¼

m m m Y Y Y 1 1 6 qð A1 : k Þ 6 H a aL k¼1 k k¼1 k¼1 k

¼

i

v

Proof. Clearly, 0 6 aLk C 6 CAk 6 aHk C, where C ¼ ½1 1    11n . Since A1 k P 0, we obtain

06

ðd2 þ e2 Þðd2 þ c2 þ l2 Þ

Similarly,

:

06

maxðvi Þ i

minðd2i þ e2i Þmind2i þ c2i þ l2i Þ i

ð

Corollary follows from Theorem 3. Theorem 4. If both A and B are non-negative square matrices, then qðAÞ 6 qðA þ BÞ. Proof. Recall that the Gelfand’s formula is that k

1 k

qðAÞ ¼ lim kA k ; k!1

C C 6 CA1 k 6 L: aHk ak

C CA1 CA1 C 1 1 k1 k1 6 6 CA A 6 6 L L : k k1 aHk aLk aHk aHk1 ak ak1

Following the same argument, m m m Y Y Y 1 1 C 6 C A1 C: k 6 H a aL k¼1 k k¼1 k¼1 k

By Corollary 1 in [40], any n  n nonnegative matrix A satisfies: n X min aij j

i¼1

!

! n X 6 qðAÞ 6 max aij : j

i¼1

ð46Þ

79

L. Xue, C. Scoglio / Mathematical Biosciences 243 (2013) 67–80

Qm

Because the entries of C k¼1 A1 k is the sum of each column of matrix Qm 1 k¼1 Ak , by Inequality (46), m m Y Y 1 6q A1 k H a k k¼1 k¼1

!

m Y 1 6 : L a k¼1 k

Proof. Note that M k is a diagonal dominant matrix of its column entries. By Theorem 1 in page 654 of [41], M k and X k are invertible. We now prove that M1 k is nonnegative. Matrix M k can be rewritten as follows. 2

2 3 fk1 xk21  xkn1 0 6 6 7 6 xk12 fk2  xkn2 7 6 xk12 6 6 7 Mk ¼ 6 7 ¼ ni¼1 fki  6 6  6     7 4 4 5 xk1n xk2n  fkn xk1n

3

xk21  xkn1 7 0  xkn2 7 7 



xk2n 

7 ¼: G  H:  7 5 0

Consequently,

2 G

¼

ni¼1 f1 ki

0

6 6 x f1 6 k12 k2 and G H ¼ 6 6  4

xk21 f1    xkn1 f1 k1 k1 0

1

1 k1n fkn

x Moreover, 0 <

Pn

1 j¼1 ðG HÞij

n n X X dS2i x2ji S2j  x2ij S2i ¼ r 2i  b12i S2i I1i =N 1i  b32i S2i I3i =N 3i  d2i S2i þ dt j¼1;j–i j¼1;j–i n X

Theorem 7. Matrices M k and X k in (24) are invertible, and M k and X k are nonnegative. Moreover, Matrices M 1 and X 1 are nonnegative k k matrices.

1

Livestock population model



3

7 7    xkn2 f1 k2 7 7:   7 5

xk2n f1  kn

0

< 1, for all i. Hence, qðG1 HÞ < 1, i.e.,

G1 H is convergent (see [28]). Obviously, G1 P 0, and G1 H P 0. By Theorem 1 in [28], M k is an M-matrix and M 1 k P 0. By the same argument, X k is an M-matrix and X 1 k P 0. This finishes the proof. Network-based RVF meta- population model Aedes population model

dP1i ¼ b1i ðN1i  q1i I1i Þ  h1i P1i dt

ð47Þ

dQ 1i ¼ b1i q1i I1i  h1i Q 1i dt

ð48Þ

dS1i ¼ h1i P1i  d1i S1i N1i =K 1  b21i S1i I2i =N2i dt

ð49Þ

dE1i ¼ b21i S1i I2i =N2i  e1i E1i  d1i E1i N1i =K 1 dt

ð50Þ

dI1i ¼ h1i Q 1i þ e1i E1i  d1i I1i N1i =K 1 dt

ð51Þ

dN1i ¼ h1i ðP1i þ Q 1i Þ  d1i N1i N1i =K 1 dt

ð52Þ

Culex population model

dP3i ¼ b3i N3i  h3i P 3i dt

ð53Þ

dS3i ¼ h3i P3i  b23i S3i I2i =N2i  d3i S3i N3i =K 3 dt

ð54Þ

dE3i ¼ b23i S3i I2i =N2i  e3i E3i  d3i E3i N3i =K 3 dt

ð55Þ

dI3i ¼ e3i E3i  d3i I3i N3i =K 3 dt

ð56Þ

dN3i ¼ h3i P3i  d3i N3i N 3i =K 3 dt

ð57Þ

n X

ð58Þ

dE2i ¼ b12i S2i I1i =N 1i þ b32i S2i I3i =N 3i  e2i E2i  d2i E2i þ x2ji E2j  x2ij E2i dt j¼1;j–i j¼1;j–i

ð59Þ

n n X X dI2i ¼ e2i E2i  c2i I2i  l2i I2i  d2i I2i þ x2ji I2j  x2ij I2i dt j¼1;j–i j¼1;j–i

ð60Þ

n n X X dR2i x2ji R2j  x2ij R2i ¼ c2i I2i  d2i R2i þ dt j¼1;j–i j¼1;j–i

ð61Þ

n n X X dN 2i x2ji N2j  x2ij N2i ¼ r 2i  l2i I2i  d2i N 2i þ dt j¼1;j–i j¼1;j–i

ð62Þ

Human population model dS4i ¼ b4i N 4i  b14i S4i I1i =N 1i  b24i S4i I2i =N 2i  b34i S4i I3i =N 3i  d4i S4i N 4i =K 4 dt dE4i ¼ b14i S4i I1i =N 1i þ b24i S4i I2i =N 2i þ b34i S4i I3i =N 3i  e4i E4i  d4i E4i N 4i =K 4 dt dI4i ¼ e4i E4i  c4i I4i  l4i I4i  d4i I4i N 4i =K 4 dt dR4i ¼ c4i I4i  d4i R4i N 4i =K 4 dt dN 4i ¼ b4i N 4i  l4i I4i  d4i N 4i N 4i =K 4 dt

ð63Þ ð64Þ ð65Þ ð66Þ ð67Þ

References [1] M. Salmani, P. van den Driessche, A model for disease transmission in a patchy environment, Discrete and Continuous Dynamical Systems – Series B 6 (2006) 185. [2] W. Wang, X.Q. Zhao, An epidemic model in a patchy environment, Mathematical Biosciences 190 (2004) 97. [3] J. Arino, J.R. Davis, D. Hartley, R. Jordan, J.M. Miller, P. van den Driessche, A multi-species epidemic model with spatial dynamics, Mathematical Medicine and Biology 22 (2005) 129. [4] R.B. Stothers, Climatic and demographic consequences of the massive volcanic eruption of 1258, Climatic Change 45 (2000) 361. [5] L.R. Petersen, J.T. Roehrig, West Nile virus: a reemerging global pathogen, Emerging Infectious Diseases 7 (2001) 611. [6] World Health Organization Multicentre Collaborative Network for Severe and Acute Respiratory Syndrome Diagnosis, A multicentre collaboration to investigate the cause of severe acute respiratory syndrome, Lancet 361 (2003) 1730–1733. [7] J. Arino, R. Jordan, P. van den Driessche, Quarantine in a multi-species epidemic model with spatial dynamics, Mathematical Biosciences 206 (2007) 46. [8] S. Busenberg, K.L. Cooke, Models of Vertically Transmitted Diseases with Sequential–Continuous Dynamics, Academic Press, 1982. [9] M. El-Doma, Analysis of an age-dependent SIS epidemic model with vertical transmission and proportionate mixing assumption, Mathematical and Computer Modelling 29 (1999) 31. [10] M.Y. Li, H.L. Smith, L.C. Wang, Global dynamics of an SEIR epidemic model with vertical transmission, SIAM Journal on Applied Mathematics 62 (2001) 58. [11] S. Busenberg, K.L. Cooke, Vertically Transmitted Diseases: Models and Dynamics, Springer, Berlin New York, 1993. [12] N. Chitnis, J.M. Hyman, C.A. Manore, Modelling vertical transmission in vectorborne diseases with applications to Rift Valley fever, Journal of Biological Dynamics 7 (2013) 11. [13] R.M. Anderson, R.M. May, The population-dynamics of micro-parasites and their invertebrate hosts, Philosophical Transactions of the Royal Society of London Series B-Biological Sciences 291 (1981) 451. [14] S. Busenberg, K.L. Cooke, The population-dynamics of two vertically transmitted infections, Theoretical Population Biology 33 (1988) 181. [15] P. Rohani, D.J. Earn, B.T. Grenfell, Opposite patterns of synchrony in sympatric disease metapopulations, Science 286 (1999) 968. [16] N.J. Gotelli, C.M. Taylor, Testing metapopulation models with stream-fish assemblages, Evolutionary Ecology Research 1 (1999) 835. [17] C.M. Taylor, R.J. Hall, Metapopulation models for seasonally migratory animals, Biology Letters 8 (2012) 477. [18] J. Arino, P. van den Driessche, The basic reproduction number in a multi-city compartmental epidemic model, Positive Systems, Proceedings 294 (2003) 135. [19] M.J. Keeling, M.E.J. Woolhouse, D.J. Shaw, L. Matthews, M. Chase-Topping, D.T. Haydon, S.J. Cornell, J. Kappey, J. Wilesmith, B.T. Grenfell, Dynamics of the 2001 UK foot and mouth epidemic: stochastic dispersal in a heterogeneous landscape, Science 294 (2001) 813. [20] O. Diekmann, J.A.P. Heesterbeek, M.G. Roberts, The construction of nextgeneration matrices for compartmental epidemic models, Journal of the Royal Society Interface 7 (2010) 873.

80

L. Xue, C. Scoglio / Mathematical Biosciences 243 (2013) 67–80

[21] O. Diekmann, J. Heesterbeek, J. Metz, On the definition and the computation of the basic reproduction ratio R0 in models for infectious diseases in heterogeneous populations, Journal of Mathematical Biology 28 (1990) 365. [22] O. Diekmann, J.A.P. Heesterbeek, Mathematical Epidemiology of Infectious Diseases, Wiley, Chichester, 2000. [23] P. van den Driessche, J. Watmough, Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission, Mathematical Biosciences 180 (2002) 29. [24] J. Li, D. Blakeley, R. Smith, The failure of R0, Computational and Mathematical Methods in Medicine (2011) 17, http://dx.doi.org/10.1155/2011/527610. [25] M. Lipsitch, M.A. Nowak, D. Ebert, R.M. May, The population dynamics of vertically and horizontally transmitted parasites, Proceedings: Biological sciences/The Royal Society 260 (1995) 321. [26] L. Xue, H.M. Scott, L.W. Cohnstaedt, C. Scoglio, A network-based metapopulation approach to model Rift Valley fever epidemics, Journal of Theoretical Biology 306 (2012) 129. [27] J. Heffernan, R. Smith, L. Wahl, Perspectives on the basic reproductive ratio, Journal of the Royal Society Interface 2 (2005) 281. [28] R.J. Plemmons, M-matrix characterizations 1: nonsingular M-matrices, Linear Algebra and its Applications 18 (1977) 175. [29] R. Flick, M. Bouloy, et al., Rift Valley fever virus, Current Molecular Medicine 5 (2005) 827. [30] V. Chevalier, M. Pépin, L. Plée, R. Lancelot, Rift Valley fever – a threat for Europe?, Euro Surveillance 15 (2010) 19506 [31] M.B. Crabtree, R.J.K. Crockett, B.H. Bird, S.T. Nichol, B.R. Erickson, B.J. Biggerstaff, K. Horiuchi, B.R. Miller, Infection and transmission of Rift Valley fever viruses lacking the NSs and/or NSm genes in mosquitoes: potential role for NSm in mosquito infection, Plos Neglected Tropical Diseases 6 (2012) e1639. [32] G. Gerdes, Rift Valley fever, Veterinary Clinics of North America – Food Animal Practice 18 (2002) 549. [33] R.F. Breiman, B. Minjauw, S.K. Sharif, P. Ithondeka, M.K. Njenga, Rift Valley fever: scientific pathways toward public health prevention and response, The American Journal of Tropical Medicine and Hygiene 83 (2010) 1. [34] K.J. Linthicum, A. Anyamba, C.J. Tucker, P.W. Kelley, M.F. Myers, C.J. Peters, Climate and satellite indicators to forecast Rift Valley fever epidemics in Kenya, Science 285 (1999) 397. [35] J. Arino, Diseases in metapopulations, Modeling and Dynamics of Infectious Diseases 11 (2009) 64. [36] D. Gao, S. Ruan, An SIS patch model with variable transmission coefficients, Mathematical Biosciences 232 (2011) 110. [37] D. Gao, S. Ruan, A multipatch malaria model with logistic growth populations, SIAM Journal on Applied Mathematics 72 (2012) 819. [38] Y.H. Hsieh, P. van den Driessche, L. Wang, Impact of travel between patches for spatial spread of disease, Bulletin of Mathematical Biology 69 (2007) 1355. [39] F.R. Adler, The effects of averaging on the basic reproduction ratio, Mathematical Biosciences 111 (1992) 89. [40] J. Cohen, Random evolutions and the spectral radius of a non-negative matrix, Mathematical Proceedings of the Cambridge Philosophical Society 86 (1979) 345.

[41] G. Cheng, X. Cheng, T. Huang, T. Tam, Some bounds for the spectral radius of the Hadamard product of matrices, Applied Mathematics E-Notes 5 (2005) 202. [42] D.V. Canyon, J.L.K. Hii, R. Muller, The frequency of host biting and its effect on oviposition and survival in Aedes aegypti, Bulletin of Entomological Research 89 (1999) 35 (Diptera: Culicidae). [43] R.O. Hayes, C.H. Tempelis, A.D. Hess, W.C. Reeves, Mosquito host preference studies in Hale County, Texas, American Journal of Tropical Medicine and Hygiene 22 (1973) 270. [44] C.J. Jones, J.E. Lloyd, Mosquitos feeding on sheep in southeastern Wyoming, Journal of the American Mosquito Control Association 1 (1985) 530. [45] L.A. Magnarelli, Host feeding patterns of Connecticut mosquitos, American Journal of Tropical Medicine and Hygiene 26 (1977) 547 (Diptera: Culicidae). [46] H.D. Pratt, C.G. Moore, Vector-borne disease control: mosquitoes of public health importance and their control, U.S. Department of Health and Human Services, Atlanta, GA, 1993. [47] M.J. Turell, C.L. Bailey, J.R. Beaman, Vector competence of a Houston Texas strain of Aedes Albopictus for Rift Valley fever virus, Journal of the American Mosquito Control Association 4 (1988) 94. [48] M.J. Turell, M.E. Faran, M. Cornet, C.L. Bailey, Vector competence of senegalese Aedes fowleri (Diptera: Culicidae) for Rift Valley fever virus, Journal of Medical Entomology 25 (1988) 262. [49] M.J. Turell, C.L. Bailey, Transmission studies in mosquitoes (Diptera: Culicidae) with disseminated Rift Valley fever virus infections, Journal of Medical Entomology 24 (1987) 11. [50] J.W. Wekesa, B. Yuval, R.K. Washino, Multiple blood feeding by Anopheles freeborni and Culex tarsalis (Diptera: Culicidae): spatial and temporal variation, Journal of Medical Entomology 34 (1997) 219. [51] B.J. Erasmus, J.A.W. Coetzer, The symptomatology and pathology of Rift Valley fever in domestic animals, Contributions to Epidemiolgy and Biostatistics 3 (1981) 77. [52] S.C. Mpeshe, H. Haario, J.M. Tchuenche, A mathematical model of Rift Valley fever with human host, Acta Biotheoretica 59 (2011) 231. [53] M. Bates, The natural history of mosquitoes, American Journal of Public Health 39 (1949) 1592. [54] C.G. Moore, R.G. McLean, C.J. Mitchell, R.S. Nasci, T.F. Tsai, C.H. Caslisher, A.A. Marfin, P.S. Moorse, D.J. Gubler, Guidelines for Arbovirus surveillance programs in the United Sates, Centers for Disease Control and Prevention, April 1993. [55] O.M. Radostits, third ed., Herd Healthy: Food Animal Production Medicine, Saunders, 2001. [56] C.J. Peters, K.J. Linthicum, Rift Valley fever, second ed., in: G.B. Beran (Ed.), Handbook of Zoonoses, Section B: viral, CRC Press, Inc., Boca Raton, Fl, 1994, p. 125. [57] J.E. Freier, L. Rosen, Vertical transmission of dengue viruses by mosquitoes of the aedes scutellaris group, The American Journal of Tropical Medicine and Hygiene 37 (1987) 640.