A Network-Based Meta-Population Approach to Model Rift Valley Fever Epidemics
Ling Xue, H. Morgan Scott, Lee. Cohnstaedt, Caterina Scoglio
AA Network-Based Meta-Population Approach to Model Rift Valley fever Epidemics
Ling Xue a , H. M. Scott b , Lee W. Cohnstaedt c , Caterina Scoglio a, ∗ a Department of Electrical & Computer Engineering,Kansas State University, U.S. 66506 b Department of Diagnostic Medicine / Pathobiology,Kansas State University, U.S. 66506 c Center for Grain and Animal Health Research,United States Department of Agriculture, U.S. 66502
Abstract
Rift Valley fever virus (RVFV) has been expanding its geographical distribution with important implications for both human andanimal health. The emergence of Rift Valley fever (RVF) in the Middle East, and its continuing presence in many areas of Africa,has negatively impacted both medical and veterinary infrastructures and human morbidity, mortality, and economic endpoints.Furthermore, worldwide attention should be directed towards the broader infection dynamics of RVFV, because suitable host,vector and environmental conditions for additional epidemics likely exist on other continents; including Asia, Europe and theAmericas. We propose a new compartmentalized model of RVF and the related ordinary di ff erential equations to assess diseasespread in both time and space; with the latter driven as a function of contact networks. Humans and livestock hosts and twospecies of vector mosquitoes are included in the model. The model is based on weighted contact networks, where nodes of thenetworks represent geographical regions and the weights represent the level of contact between regional pairings for each setof species. The inclusion of human, animal, and vector movements among regions is new to RVF modeling. The movementof the infected individuals is not only treated as a possibility, but also an actuality that can be incorporated into the model. Wehave tested, calibrated, and evaluated the model using data from the recent 2010 RVF outbreak in South Africa as a case study;mapping the epidemic spread within and among three South African provinces. An extensive set of simulation results shows thepotential of the proposed approach for accurately modeling the RVF spreading process in additional regions of the world. Thebenefits of the proposed model are twofold: not only can the model di ff erentiate the maximum number of infected individualsamong di ff erent provinces, but also it can reproduce the di ff erent starting times of the outbreak in multiple locations. Finally, theexact value of the reproduction number is numerically computed and upper and lower bounds for the reproduction number areanalytically derived in the case of homogeneous populations. Keywords: networks, meta-population, deterministic model, Rift Valley fever (RVF), mitigation,
Aedes mosquitoes,
Culex mosquitoes
1. Introduction
Rift Valley fever (RVF) is a viral zoonosis with enormous health and economic impacts on domestic animals and humans[26], in countries where the disease is endemic and in others where sporadic epidemics and epizootics have occurred. Anoutbreak in South Africa in 1951 was estimated to have infected 20 ,
000 people and killed 100 ,
000 sheep and cattle [12, 35].In Egypt in 1977, there were 18 ,
000 human cases with 698 deaths resulting from the disease [12, 35]. While RVF is endemicin Africa, it also represents a threat to Europe and Western hemispheres [7, 18]. In 1997 − ,
000 human cases and 478 death. The first recorded outbreak outside of Africa occurred in the Arabianpeninsula in 2000 − − . − ∗ Corresponding Author
Email addresses: [email protected] (Ling Xue), [email protected] (H. M. Scott),
[email protected] (Lee W. Cohnstaedt), [email protected] (Caterina Scoglio) a r X i v : . [ q - b i o . O T ] S e p edes and Culex are believed to be the main arthropod vectors [7] . Rift Valley fever virus can be transferred vertically fromfemales to their eggs in some species of the
Aedes mosquitoes [18, 27]. The disease has been shown to be endemic in semi-aridzones, such as northern Senegal [45, 7, 28], and RVF epidemics often appears at 5 −
15 year cycles [28]. As noted earlier, RVFVhas already spread outside Africa, to Yemen and Saudi Arabia [7, 6]. The species of vectors that are capable of transmittingRVFV have a wide global distribution [20] and there is therefore a distinct possibility for the virus to spread out of its currentlyexpanding geographic range [9]. A pathways analysis [21] has shown that the RVF virus might be introduced into the UnitedStates in several di ff erent ways [21, 24] and that analysis identified several regions of the United States that are most susceptibleto RVFV introduction. It is therefore desirable to develop e ff ective models to better understand the potential dynamics of RVF inheretofore una ff ected regions and then develop e ffi cient mitigation strategies in case this virus appears in the Western hemisphere[18]. Such preparedness can help avoid a rapid spread of the virus throughout North America, as happened with the West Nilevirus during the last decade [7, 18].A RVF disease risk mapping model was developed by [1]. The authors observed sea surface temperature (SST) patterns,cloud cover, rainfall, and ecological indicators (primarily vegetation) via satellite data to evaluate di ff erent aspects of climatevariability and their relationships to disease outbreaks in Africa and the Middle East [3, 2]. The researchers successfully predictedareas where outbreaks of RVF in humans and animals were expected using climate data for the Horn of Africa from December2006 to May 2007. An ordinary di ff erential equation (ODE) mathematical model was developed by [18]. The model is bothan individual-based and deterministic model. The authors analyzed the stability of the model and tested the importance of themodel parameters. However, neither human population parameters nor spatial (or, network) aspects are explicitly incorporatedin the model. Another theoretical mathematical model on RVFV dynamic transmission was proposed [29]. This model is also anindividual based model. The most important parameters to initial disease transmission and the endemic equilibrium have beencarried out.In this paper, we present a novel model incorporating Aedes and
Culex mosquito vector, and livestock and human hostpopulations. Our model is based on weighted contact networks, where nodes of the networks represent geographical regionsand weights represent the level of contact between regional pairs for each vector or host species. Environmental factors such asrainfall, temperature, wind and evaporation are incorporated into the model. For each subpopulation, a set of ordinary di ff erentialequations describes the dynamics of the population in a specific geographical location, and the transitions among the di ff erentcompartments, after contracting the virus. We compute the lower and upper bounds of the reproduction number for homogeneouspopulations, explain their biological meaning, and numerically compare the bounds with exact values.We test, calibrate, and evaluate the model using the recent 2010 RVF outbreak in South Africa as a case study, mapping theepidemic spread in three South African provinces: Free State, Northern Cape, and Eastern Cape. An extensive set of simulationresults shows the potential of the proposed approach to accurately describe the spatialtemporal evolution of RVF epidemics.The paper is organized as follows: 1) in Section 2, we describe our compartmentalized mathematical model, present thelower bound and upper bound of the reproduction number for homogeneous populations. 2) in Section 3, we introduce the casestudy using outbreak data from South Africa, 2010, 3) in Section 4, we conclude our work. In the Appendix, we show how wederive the bounds for the reproduction number for homogeneous populations.
2. Compartmentalized Mathematical Model
We have constructed Compartmentalized Mathematical Models based on the principle of RVFV transmission. The parame-ters used in the model are shown in Table 1.
The principle of RVFV transmission between di ff erent species is shown in Figure 1. All the Aedes ssp. and Culex spp. we aregoing to discuss only include the mosquitoes that are competent vectors of Rift Valley fever. In this paper, the Culex parametersare based on Culex Tarsalis mosquitoes and Aedes parameters are based on Ae. dorsalis mosquitoes. The main vectors, Aedes and
Culex mosquitoes and the main hosts, livestock and humans are considered in the model. We use an SEI compartmentalmodel in which individuals are either in a susceptible (S) state, an exposed (E) state, or an infected state (I) for both
Aedes and
Culex mosquitoes, and an SEIR compartmental model in which individuals are either in a susceptible (S) state, an exposed (E)state, an infected state (I), or a recovered (R) state for both livestock and human populations. Infectious
Aedes mosquitoes cannot only transmit RVFV to susceptible livestock and humans but also to their own eggs [18, 27].
Culex mosquitoes acquire thevirus during blood meals on an infected animal and then amplify the transmission of RVFV through blood meals on livestockand humans [44]. Direct ruminant-to-human contact is the major (though not only) way for humans to acquire the infection[1, 11]. Accidental RVFV infections have been recorded in laboratory sta ff handling blood and tissue from infected animals[1]. Usually, humans are thought of as dead end hosts that do not contribute significantly to propagation of the epidemic [7].There has been no direct human-to-human transmission of RVFV in field conditions recorded thus far [21]. The mosquitoeswill not spontaneously recover once they become infectious [18]. Livestock and humans either perish from the infection orrecover [18]. All four species have a specified incubation period [44]. The model is based on a daily time step. Aedes and
Culex S a , exposed E a and infected I a compartments. The subscript a = Aedes and a = Culex . The size of each adult mosquito population is N = S + E + I for adult Aedes mosquitoes and N = S + E + I for adult Culex mosquitoes. The livestock and human hosts contain susceptible S b , exposed E b , infected I b and recovered R b individuals. The subscript b = b = N b = S b + E b + I b + R b . The four populations are modeled with a specified carrying capacity K , K , K , K respectively. Figure 1: Flow diagram of RVFV transmission with each species, namely,
Aedes mosquitoes,
Culex mosquitoes, livestock, and humans homogeneously mixed(the solid lines represent transition between compartments and the dash lines represent the transmission between di ff erent species) Mosquito Population Model d P d t = b ( N − q I ) − θ P (1)d Q d t = b q I − θ Q (2)d S d t = θ P − β S I / N − d S N / K (3)d E d t = β S I / N − ε E − d E N / K (4)d I d t = θ Q + ε E − d I N / K (5)3arameter Description Value Units Source β contact rate: Aedes to livestock 0 .
002 1 / day [18] β contact rate: livestock to Aedes . / day [18] β contact rate: livestock to Culex . / day [18] β contact rate: Culex to livestock 0 . / day [18] β contact rate: Aedes to humans 0 . / day Assume β contact rate: livestock to humans 0 . / day [29] β contact rate: Culex to humans 0 . / day Assume γ recover rate in livestock 0 .
14 1 / day [18] γ recover rate in humans 0 .
14 1 / day [39, 38, 37] d death rate of Aedes mosquitoes 0 .
025 1 / day [18] d death rate of livestock 1 / / day [18] d death rate of Culex mosquitoes 0 .
025 1 / day [18] d death rate of humans 1 / / day [39, 38, 37] b number of Aedes eggs laid per day 0 .
05 1 / day [18] b daily birthrate of livestock 0 . / day [18] b number of Culex eggs laid per day weather dependent 1 / day [19] b daily birthrate of humans 1 / / day [39, 38, 37]1 /(cid:15) incubation period in Aedes mosquitoes 6 days [18]1 /(cid:15) incubation period in livestock 4 days [18]1 /(cid:15) incubation period in Culex mosquitoes 6 days [18]1 /(cid:15) incubation period in humans 4 days [44] µ mortality rate in livestock 0 . / day [18] µ mortality rate in humans 0 . / day [39, 38, 37] q transovarial transmission rate in Aedes .
05 - [18]1 /θ development time of Aedes
15 days [18] θ development rate of Culex weather dependent 1 / day [19] K carrying capacity of Aedes mosquitoes 1000000000 - [32] K carrying capacity of livestock 10000000 - Assume K carrying capacity of Culex mosquitoes 1000000000 - [32] K carrying capacity of humans 10000000 - Assume f fraction of those working with animals 0 .
82 - [31] τ return rate 3 times / day [5] p reduction in ω i j due to infection - Assume Table 1: Parameters of the compartmentalized mathematical model N d t = θ ( P + Q ) − d N N / K (6)where: P = the number of uninfected Aedes mosquito eggs Q = the number of infected Aedes mosquito eggs S = the number of susceptible Aedes mosquitoes E = the number of exposed Aedes mosquitoes I = the number of infected Aedes mosquitoes N = the total number of Aedes mosquitoesThe above model is a modified SEI model with compartments P and Q. Compartments P and Q represent uninfected eggsand infected eggs respectively. The total number of eggs laid each day is b N with b q I infected eggs and b N − b q I unin-fected eggs [18]. After development period, θ P develop into susceptible adult mosquitoes and θ Q develop into infected adultmosquitoes [18]. There are d X N / K mosquitoes removed from compartment X due to natural death. Compartment X can beP, Q, S, E, and I here. The number of Aedes mosquitoes infected by livestock is denoted by β S I / N which is proportional tothe density of infected Aedes mosquitoes [18]. After incubation period, ε E Aedes mosquitoes transfer to infected compartment[18].
Mosquito Population Model d P d t = b N − θ P (7)d S d t = θ P − β S I / N − d S N / K (8)d E d t = β S I / N − ε E − d E N / K (9)d I d t = ε E − d I N / K (10)d N d t = θ P − d N N / K (11)where: P = the number of uninfected Culex mosquito eggs S = the number of susceptible Culex mosquitoes E = the number of exposed Culex mosquitoes I = the number of infected Culex mosquitoes N = the total number of Culex mosquitoesBesides compartment S , E , I , compartment P is added to represent uninfected eggs. Only uninfected eggs are included be-cause the female Culex mosquitoes do not transmit RVFV vertically [18]. The total number of eggs laid each day is b N . Thereare d X N / K Culex mosquitoes removed due to natural death. Compartment X can be P, S, E, and I here. After developmentperiod, θ P eggs develop into susceptible adult Culex mosquitoes and become secondary vectors [18]. The number of infected
Culex mosquitoes from contact with livestock is denoted by β S I / N which is proportional to the percentage of infected live-stock [18]. After incubation period, ε E Culex mosquitoes transfer from exposed compartment to infected compartment [18]. d S d t = b N − d S N / K − β S I / N − β S I / N (12)d E d t = β S I / N + β S I / N − ε E − d E N / K (13)d I d t = ε E − γ I − µ I − d I N / K (14)d R d t = γ I − d R N / K (15)5 N d t = b N − d N N / K − µ I (16)where: S = the number of susceptible livestock E = the number of exposed livestock I = the number of infected livestock N = the total number of livestockThere are b N livestock born, d X N / K livestock removed due to natural death [18], and µ I livestock dying of the infec-tion each day [18]. Compartment X can be S, E, I, and R here. Following incubation period, ε E livestock transfer from exposedcompartment to infected compartment [18]. The number of livestock infected by Aedes mosquitoes and
Culex mosquitoes aredenoted as β S I / N and β S I / N respectively [18]. Following infection period, γ I livestock recover from RVFV infec-tion [18]. d S d t = b N − β S I / N − f β S I / N − β S I / N − d S N / K (17)d E d t = β S I / N + f β S I / N + β S I / N − d E N / K − ε E (18)d I d t = ε E − γ I − µ I − d I N / K (19)d R d t = γ I − d R N / K (20)d N d t = b N − d N N / K − µ I (21)where: S = the number of susceptible humans E = the number of exposed humans I = the number of infected humans N = the total number of humansThere are b N humans born, d X N / K humans removed from compartment X due to natural death, and µ I humansdying of RVFV infection each day. Compartment X can be S, E, I, and R here. The number of humans that acquire the infectionfrom Aedes mosquitoes,
Culex mosquitoes, and livestock is β S I / N , β S I / N , and f β S I / N respectively. We assumeonly those who work with animals can be infected by animals. Therefore, a coe ffi cient f (0 < f <
1) which represents thefraction of humans working with animals is multiplied by β S I / N . After incubation period, ε E humans transfer to infectedcompartment and γ I humans transfer to recovered compartment after infection period. The equation (22) is used to model the development rate of
Culex mosquitoes [19]. The daily egg laying rate expressed inequation (23) is a function of moisture [19]. Moisture in equation (24) is obtained by summing the di ff erence of precipitation[30] and evaporation (mm) [25] over the proceeding 7 days [19]. In the equations (22) to (25), A , HA , HH , K , T H , E max , E var , E mean , b are parameters [19] which are described in Table 2. This model is specific for West Nile virus model in 2010 in thenorthern US. θ ( T emp , t ) = A ∗ ( T emp ( t ) + K )298 . ∗ exp [ HA . ∗ ( . − Temp ( t ) + K )]1 + exp [ HH . ∗ ( TH − Temp ( t ) + K )] (22) b ( T emp , precipitation , t ) = b + Emax + exp [ − Moisture ( t ) − EmeanEvar ] (23)6arameter Description Value Source A parameter of the development rate 0 .
25 [19] HA parameter of the development rate 28094 [19] HH parameter of the development rate 35692 [19] T H parameter of the development rate 298 . b minimum constant fecundity rate 3 [19] E max maximum daily egg laying rate 3 [19] E mean value at which moisture index = . E max E var the variance of the daily egg laying rate 12 [19] Table 2: Parameters of the model for
Culex
Moisture ( t ) = t (cid:88) D = t − participation ( D ) − evaporation ( D ) (24) Evaporation ( t ) = T emp ( t ) + . h ) / (100 − latitude )80 − T emp ( t ) + T emp ( t ) − T d ( t ))80 − T emp ( t ) mm / day (25)Where: T emp ( t ) = air temperature in units of o C [25] latitude = the latitude (degrees) [25] T d ( t ) = the mean dew-point in units of o C [25] h = the elevation (meters) [25] K = Kelvin parameter [25]
The reproduction number R is defined as: “ The average number of secondary cases arising from an average primary casein an entirely susceptible population” [14]. The reproduction number is used to predict whether the epidemic will spread or dieout. There are several methods used to compute R . One of these methods computes the reproduction number as the spectralradius [14 ? ] of the next generation matrix [14 ? ].The next generation matrix is defined as FV − , and the matrices F and V are determined as: F = [ ∂ F i ( x ) ∂ x j ] , V = [ ∂ V i ( x ) ∂ x j ]where x j is the number or proportion of infected individuals in compartment j , j = , , , ... , m , m being the total number ofinfected compartments, x is the disease free equilibrium vector, F i is the rate of appearance of new infections in compartment i , and V i = V − i − V + i with V − i denoting the transfer of individuals out of compartment i and V + i denoting the rate of transferof individuals into compartment i [16]. The ( i , j ) entry of F represents the rate at which infected individuals in compartment j produce infected individuals in compartment i [16]. The ( j , k ) entry of V − represents the average time that an individual spendsin compartment j , where i , j , k = , , , ... , m [16]. Finally, the ( i , k ) entry of FV − represents the expected number of infectedindividuals in compartment i produced by the infected individuals in compartment k [16].For our homogeneous population model, we found that R H (cid:54) R (cid:54) R H + q (26)where R H = (cid:115) ε ( b + ε )( b + γ + µ ) (cid:104) ε β β b ( b + ε ) + ε β β b ( b + ε ) (cid:105) (27)See the Appendix for the derivation details, the biological interpretation, and the comparison among exact values and boundsfor the reproduction number. 7 .2. Meta-Population Model A meta-population model is a model with several subpopulations. It assumes homogeneity within each subpopulation andheterogeneity among di ff erent subpopulations. The Aedes and
Culex mosquitoes in location i ( i = , , , · · · , n ), are distributedamong susceptible S ai , exposed E ai and infected I ai compartments. The subscript a = Aedes and a = Culex .The size of each adult mosquito population in location i is N i = S i + E i + I i for adult Aedes mosquitoes and N i = S i + E i + I i for adult Culex mosquitoes. The livestock and human hosts contain susceptible S bi , exposed E bi , infected I bi and recovered R bi individuals. The subscript b = b = i is N i = S i + E i + I i + R i for livestock hosts and N i = S i + E i + I i + R i for human hosts. The four populationsare modeled with a specified carrying capacity K , K , K , K respectively. We used weighted networks for each compartment of the four species as is shown in Figure 2. The superscripts of ω on theleft hand side of equations (28), (29), and (31) represent the movement of di ff erent species. The number ’1’ in the superscriptmeans the movement of Aedes or Culex population, ’2’ means the livestock movement, and ’3’ means the human movement inthe superscript. The subscript i j of ω i j means that the direction of the movement is from location i to location j . The di ff erencein the thickness of the lines represent the di ff erence in weight. Thicker lines represent the larger weight. The weight for eachpopulation is between 0 and 1. RVFV has been documented to be spread by wind [35]. Wind dispersal of mosquitoes haschanged geographic distribution and accelerated the spread of RVFV to new geographic areas [21]. Some locations can becomesecondary epidemic sites after the virus has been introduced (especially in irrigated areas, e.g. Gazeera in Sudan or rice valleys inthe center of Madagascar) [28]. Livestock trade and transport also can a ff ect the geographic distribution of RVF [7]. One criticalobjective in developing e ff ective models is to determine the major factors involved in the disease spreading process. Therefore,we parameterize the weight due to mosquito movement with wind [21, 8], livestock movement due to transportation to feedlotsor trade centers [40], and human mobility due to commuting [5] as shown in equations (28), (29), and (31), respectively. Themovement rate of infected livestock is reduced due to infection [44]. We use the wind data [41] in Bloemfontein, which is thecapital of Free State, as the wind of Free State Province, that of Kimberley, which is the capital of Northern Cape, as the wind ofNorthern Cape Province and that of Grahams town, which is the center of Eastern Cape Province, as the wind of Eastern CapeProvince. The distance vector is calculated with longitude and latitude in the center of each location. The number of animals sold[37] and the number of livestock in the feedlots [34] are factors of weight for livestock movement. Distance, human population,commuting rate, and return rate [39] a ff ect the weight for human movement. Weight for mosquito movement is decided bydistance and the projection of wind in the direction of distance vector [8]. ω i j = c (cid:126) W i · (cid:126) D i j | (cid:126) D i j | | (cid:126) D i j | (28) ω i j = c F M j F M i | (cid:126) D i j | (29) σ i j = c N α i N γ j e β | (cid:126) D ij | (30) ω i j = σ i j N i (31) ω i = n (cid:88) j = , j (cid:44) i ω i j (32)Here: (cid:126) W i = the wind vector in location i [8] (cid:126) D i j = the distance vector from location i to location j ω i j ( t ) = the weight for mosquitoes moving from location i to location j ω i j ( t ) = the weight for livestock moving from location i to location j σ i j ( t ) = the number of commuters between location i and location jF M i = the number of animals in markets and feedlots in location i2.2.2. Aedes Movement between Nodesd P i d t = b ( N i − q I i ) − θ P i (33)8 a) Mosquito movement network (Mosquitoes canmove from node i to node j , j , and j and vice versadue to wind. We assume mosquitoes do not return tothe node they are from.) (b) Livestock movement network (Livestock can movefrom node i to node j , j , and j and vice versa dueto trade. We assume livestock do not return to the nodethey are from.)(c) Human movement network (Humans can commutefrom node i to node j , j , and j and vice versa. Weassume humans return to the node they are from.)Figure 2: Network graphs with node i which has three neighbors as an example Q i d t = b q I i − θ Q i (34)d S i d t = θ P i − β S i I i / N i − d S i N i / K + n (cid:88) j = , j (cid:44) i ω ji S j − n (cid:88) j = , j (cid:44) i ω i j S i (35)d E i d t = β S i I i / N i − ε E i − d E i N i / K + n (cid:88) j = , j (cid:44) i ω ji E j − n (cid:88) j = , j (cid:44) i ω i j E i (36)d I i d t = θ Q i + ε E i − d I i N i / K + n (cid:88) j = , j (cid:44) i ω ji I j − n (cid:88) j = , j (cid:44) i ω i j I i (37)d N i d t = θ ( P i + Q i ) − d N i N i / K + n (cid:88) j = , j (cid:44) i ω ji S j − n (cid:88) j = , j (cid:44) i ω i j S i + n (cid:88) j = , j (cid:44) i ω ji E j − n (cid:88) j = , j (cid:44) i ω i j E i + n (cid:88) j = , j (cid:44) i ω ji I j − n (cid:88) j = , j (cid:44) i ω i j I i (38)The change in the number of Aedes mosquitoes due to mobility in compartment X is given as (cid:80) nj = , j (cid:44) i ω ji X j − (cid:80) nj = , j (cid:44) i ω i j X i [22]. Movement between Nodesd P i d t = b N i − θ P i (39)d S i d t = θ P i − β S i I i / N i − d S i N i / K + n (cid:88) j = , j (cid:44) i ω ji S j − n (cid:88) j = , j (cid:44) i ω i j S i (40)d E i d t = β S i I i / N i − ε E i − d E i N i / K + n (cid:88) j = , j (cid:44) i ω ji E j − n (cid:88) j = , j (cid:44) i ω i j E i (41)d I i d t = ε E i − d I i N i / K + n (cid:88) j = , j (cid:44) i ω ji I j − n (cid:88) j = , j (cid:44) i ω i j I i (42)d N i d t = θ P i − d N i N i / K + n (cid:88) j = , j (cid:44) i ω ji S j − n (cid:88) j = , j (cid:44) i ω i j S i + n (cid:88) j = , j (cid:44) i ω ji E j − n (cid:88) j = , j (cid:44) i ω i j E i + n (cid:88) j = , j (cid:44) i ω ji I j − n (cid:88) j = , j (cid:44) i ω i j I i (43)The change in the number of Culex mosquitoes in compartment X due to movement is given as (cid:80) nj = , j (cid:44) i ω ji X j − (cid:80) nj = , j (cid:44) i ω i j X i [22]. d S i d t = b N i − β S i I i / N i − β S i I i / N i − d S i N i / K + n (cid:88) j = , j (cid:44) i ω ji S j − n (cid:88) j = , j (cid:44) i ω i j S i (44)d E i d t = β S i I i / N i + β S i I i / N i − ε E i − d E i N i / K + n (cid:88) j = , j (cid:44) i ω ji E j − n (cid:88) j = , j (cid:44) i ω i j E i (45)d I i d t = p n (cid:88) j = , j (cid:44) i ω ji I j − p n (cid:88) j = , j (cid:44) i ω i j I i − d I i N i / K + ε E i − γ I i − µ I i (46)d R i d t = n (cid:88) j = , j (cid:44) i ω ji R j − n (cid:88) j = , j (cid:44) i ω i j R i + γ I i − d R i N i / K (47)10 N i d t = b N i − d N i N i / K − µ I i + n (cid:88) j = , j (cid:44) i ω ji S j − n (cid:88) j = , j (cid:44) i ω i j S i + n (cid:88) j = , j (cid:44) i ω ji E j − n (cid:88) j = , j (cid:44) i ω i j E i + p n (cid:88) j = , j (cid:44) i ω ji I j − p n (cid:88) j = , j (cid:44) i ω i j I i + n (cid:88) j = , j (cid:44) i ω ji R j − n (cid:88) j = , j (cid:44) i ω i j R i (48)The change in the number of animals due to movement in susceptible, exposed, and recovered compartment is (cid:80) nj = , j (cid:44) i ω ji X j − (cid:80) nj = , j (cid:44) i ω i j X i [22] for livestock. Concerning the animals in the infected compartment, we assume that the movement rate of theinfected livestock is p (0 < p <
1) of livestock in other compartments. This value of the movement rate has been selected in theabsence of further information. d S i d t = b N i − d S i N i / K − β S i I i / N i + σ i /τ − β f S i I i / N i + σ i /τ − β S i I i / N i + σ i /τ − n (cid:88) j = , j (cid:44) i β S i I j / N j σ i j /τ + σ i /τ − n (cid:88) j = , j (cid:44) i β f S i I j / N j σ i j /τ + σ i /τ − n (cid:88) j = , j (cid:44) i β S i I j / N j σ i j /τ + σ i /τ (49)d E i d t = β S i I i / N i + σ i /τ + β f S i I j / N i + σ i /τ + β S i I i / N i + σ i /τ + n (cid:88) j = , j (cid:44) i β S i I j / N j σ i j /τ + σ i /τ + n (cid:88) j = , j (cid:44) i β f S i I j / N j σ i j /τ + σ i j /τ + n (cid:88) j = , j (cid:44) i β S i I j / N j σ i j /τ + σ i /τ − d E i N i / K − ε E i (50)d I i d t = ε E i − γ I i − µ I i − d I i N i / K (51)d R i d t = γ I i − d R i N i / K (52)d N i d t = b N i − d N i N i / K − µ I i (53)The humans from location i can stay in location i or move to location j at time t [5]. The number of humans infected by Aedes mosquitoes,
Culex mosquitoes and livestock is β ( S ii I i N i + (cid:80) nj = , j (cid:44) i S i j I j N j ) [5], β ( S ii I i N i + (cid:80) nj = , j (cid:44) i β S i j I j N j ) [5], and f β ( S ii I i N i + (cid:80) nj = , j (cid:44) i S i j I j N j ) [5] respectively.where: S ii = the number of humans that are from location i and stay in location i at time t [5]. S i j = the number of humans that are from location i and stay in location j at time t [5]. ω i j = the commuting rate between subpopulation i and each of its neighbor j [5] ω i = daily total rate of commuting for population i [5]The change in the number of susceptible humans that are from location i and stay in location i is given [5] by the followingexpression. ∂ S ii ∂ t = n (cid:88) j = , j (cid:44) i τ S i j − n (cid:88) j = , j (cid:44) i ω i j S ii The change in the number of susceptible humans that are from location i and stay in location j is given [5] by the followingexpression. ∂ S i j ∂ t = ω i j S ii − τ S i j We can get the solution of S ii and S i j through the above two equations at the equilibrium. S ii = S i + ω i /τ (54) S i j = S i + ω i /τ ω i j /τ (55)11 . Case Study: South Africa 2010 We have used data from the South African RVF outbreak in 2010 as a case study.
Outbreak data for animals are obtained from [15, 43], while outbreak data for human subpopulations are collected from[13, 31]. As far as animal data is concerned, we chose to analyze RVF incidence in the sheep population. Because the granularityof human incidence data is provided at Province level, each node in the network represents a province. We selected threeprovinces: Free State (location 1), Northern Cape (location 2) and Eastern Cape (location 3), because they had the highest levelsof RVF incidence for humans. The curves of the incidence data are shown in Figure 3 using green histograms, while the redcurves represent simulations obtained with our model. From the data in Figure 3, it is possible to observe that the epidemicstarted first in the Free State Province and later in Northern Cape Province. The sustained heavy rainfall likely triggered theoutbreak, causing infected eggs to hatch in the Free State Province. Additionally, the number of animal and human cases inEastern Cape Province is smaller than the other two provinces.
The three parameters c , c and c are estimated using the least square approach, which is based on minimization of errorsbetween the incidence data of humans and the percentage of humans calculated by the mathematical model. At first, we establishan objective function. At each sample time, we calculate the di ff erence between the number of humans calculated by di ff erentialequations and that reported [36] during outbreaks in three provinces of South Africa from January, 2010. We calculate the squareof each di ff erence. Then, we add all the squares for each location in each day together to obtain the objective function as is shownbelow. Minimization of the objective function is initiated by providing initial values c , c and c for each parameter. Thedi ff erential equations are solved with each set of the parameters and the square errors between the number of infected humansobtained from the objective function and those from incidence data are calculated. The parameters c , c and c we used in themodel are c = . c = .
05 and c = . F = t f (cid:88) t = t n (cid:88) i = [( I i ( t ) − PR i ( t )) ] (56)In the equations above, n = the number of nodes t = starting time t f = end time I i ( t ) = human prevalence calculated by the model PR i ( t ) = human prevalence reportedTo conduct a sensitivity analysis of the parameters c , c , and c in equations (28), (29), and (31), we have changed eachparameter within ±
10% of the values c = . c = .
05 and c = . c = . c = .
05 and c = .
005 is represented as I Oi ( t ). The percentage of infected human obtainedfrom simulation with the parameters within ±
10% bound is represented as I i ( t ). The relative errors between the fractions ofinfected humans are calculated for each set of parameters, in each location, at time t as | I i ( t ) − I Oi ( t ) I i ( t ) | . The relative errors and thelower bound and upper bound of the parameters are shown in Figure 4.All the values of relative errors shown in Figure 4 are smaller than 10%, proving the model robustness with respect to limiteduncertainties in the parameter estimation. The rest of the parameters such as contact rate β , β , β , β , death rate d , d and recovery rate γ are the most significant parameters in [18]. Similarly, β , β , β and γ are also the most significantparameters in this model. To explore the behavior of RVFV, we conducted numerical simulations of an open system considering movement of thefour species among di ff erent locations. To test the validity of the model, we changed some parameters in the weights to seethe impact of each variation. If the number of infected eggs Q = Q =
0, and Q =
0, at the beginning infected eggsonly exist in location 1, Free State. However, our model considers movement of mosquitoes to other locations with wind. As aconsequence, infected animals and humans appear in all three locations as is shown in Figure 5. Therefore, the infection spreadsdue to movement of the four populations. 12
Time (days) N u m be r o f i n f e c t ed an i m a l s (a) Simulation result and incidence data for sheep in FreeState Province
50 100 150 20005101520 N u m be r o f i n f e c t ed hu m an s Time (days) (b) Simulation result and incidence data for humans inFree State Province
50 100 150 200050010001500200025003000 N u m be r o f i n f e c t ed an i m a l s Time (days) (c) Simulation result and incidence data for sheep inNorthern Cape Province
50 100 150 20005101520 N u m be r o f i n f e c t ed hu m an s Time (days) (d) Simulation result and incidence data for humans inNorthern Cape Province
50 100 150 200050010001500200025003000 N u m be r o f i n f e c t ed an i m a l s Time (days) (e) Simulation result and incidence data for sheep in East-ern Cape Province
50 100 150 20005101520 N u m be r o f i n f e c t ed hu m an s Time (days) (f) Simulation result and incidence data for humans inEastern Cape ProvinceFigure 3: Simulation results and incidence data from January, 2010 in South Africa (bars are data and lines are simulation results) −3 The value of c R e l a t i v e e rr o r (a) Relative error of the number ofinfected humans with di ff erent valueof c The value of c R e l a t i v e e rr o r (b) Relative error of the number ofinfected humans with di ff erent valueof c −3 −16 The value of c R e l a t i v e e rr o r (c) Relative error of the number ofinfected humans with di ff erent valueof c Figure 4: The relative error of the number of infected humans with changing one of the parameters c , c , and c Time(days) N u m be r o f c a s e s (a) Free State province N u m be r o f c a s e s Time(days) (b) Northern Cape province N u m be r o f c a s e s Time(days) (c) Eastern Cape provinceFigure 5: Simulation results with nonzero movement weights (the solid line represents livestock with y-axis on the left and the dash line represents humans withy-axis on the right)
If we also assume that at the beginning infected eggs only exist in location 1, Q = Q = Q =
0, and movementsof the four species from one location to another are not allowed, c = c =
0, and c =
0, then infected animals and humanswill not appear in location 2 and location 3 as is shown in Figure 6. We can test the mitigation strategy of movement ban withthis model. We performed the simulations to reproduce the RVF outbreak in the three South African Provinces. The simulationresults and the incidence data are shown in Figure 3. The model can di ff erentiate the maximum number of infected individualsamong the three di ff erent provinces, it also reproduces the di ff erent starting time of the outbreak in the three locations. With ahomogeneous population model, such as the one in [18], the spatial di ff erentiation is not possible.The animal incidence curves provided by the model were always an overestimation of the data, since underreporting isvery common during outbreaks. Finally, our approach in which the fractions of each subpopulation in each compartment areexpressed as continuous variables, requires a large number of cases to be accurate. For this reason, the incidence data for location3, Eastern Cape, are better approximated by a stochastic model. The model has shown the ability of fitting the data. The startingtime and trend of outbreak dynamics have been reproduced by the model.
4. Conclusions
A meta-population, network-based, deterministic RVF model is presented here. The animal, human and mosquito movementand their spatial distribution are considered by the model. The model successfully describes a real outbreak dynamics of RVFV,taking into account space and movement. When considering n locations or nodes ( n > = ∗ n di ff erential equationsand 21 ∗ n variables in our model, while there are only 14 equations with 14 variables in the model presented in [29] and 16 equa-tions with 16 variables presented in [18]. Greater accuracy of our model is obtained at the cost of an increased complexity. Thenovelty of our model is that it considers a weighted contact network to represent the movement of four species. Subpopulations atthe node level are also incorporated in our model. Additionally, parameters representing mosquito propagation and developmentare not constant but are the functions of climate factors. The model has been evaluated using data from the recent outbreak inSouth Africa. We reproduced not only the starting time but also the trend of RVFV transmission with time in di ff erent locations.14
50 100 150 200010002000
Time(days) N u m be r o f c a s e s (a) Free State province N u m be r o f c a s e s Time(days) (b) Northern Cape province
Time(days) N u m be r o f c a s e s (c) Eastern Cape provinceFigure 6: Simulation results with c = c = c = The model has shown to be very promising notwithstanding the limitation of the data. Due to the flexibility and accuracy ofthe proposed model, we can test and design multiple and di ff erent mitigation strategies in di ff erent locations at di ff erent times.The lower bound and upper bound of the reproduction number for homogeneous populations are shown to be very close to theexact value, and they provide insights on the biology of the spreading process. Future work in follow-up mathematical modelsincludes the development of a stochastic model, the study of the impact of climate changes on the epidemiology and control ofRVF, and the improvement of the mosquito movement model considering di ff usion equations. Moreover, the carrying capacitiesof mosquitoes will be considered dependent on climate factors in the future. Acknowledgments
This work has been supported by the DHS Center of Excellence for Emerging and Zoonotic Animal Diseases (CEEZAD),and by the National Agricultural Biosecurity Center (NABC) at Kansas State University. We would like to give special thanksto the anonymous editor and reviewer for their comments. We are grateful to Jason Coleman, Regina M. Beard, Livia Olsen,and Kelebogile Olifant for their help on bibliography research. We would like to give thanks to Duygu Balcan, Phillip Schumm,Faryad Darabi Sahneh, Anton Lyubinin, and Getahun Agga for the help in producing the paper, and Kristine Bennett for answer-ing questions on entomology.
Appendix
Exact Computation of R We compute R as the spectral radius of the next generation matrix of the entire system [14, 29], R = ρ ( FV − ). Beforeapplying the method we need to verify that the five assumptions in [16] are satisfied [23]. First, the equations in the system arereordered so that the first m ( m =
9) compartments correspond to infected individuals.d Q d t = b q I − θ Q (57)d E d t = β S I / N − ε E − d E N / K (58)d I d t = θ Q + ε E − d I N / K (59)d E d t = β S I / N + β S I / N − ε E − d E N / K (60)d I d t = ε E − γ I − µ I − d I N / K (61)d E d t = β S I / N − ε E − d E N / K (62)d I d t = ε E − d I N / K (63)d E d t = β S I / N + f β S I / N + β S I / N − d E N / K − ε E (64)15 I d t = ε E − γ I − µ I − d I N / K (65)d P d t = b ( N − q I ) − θ P (66)d P d t = b N − θ P (67)d S d t = θ P − β S I / N − d S N / K (68)d S d t = b N − d S N / K − β S I / N − β S I / N (69)d S d t = θ P − β S I / N − d S N / K (70)d S d t = b N − β S I / N − f β S I / N − β S I / N − d S N / K (71)d R d t = γ I − d R N / K (72)d R d t = γ I − d R N / K (73)The above system can be written as f k ( x ) = F k ( x ) − V k ( x ), k = x = [ Q E I E I E I E I P P S S S S R R ] T , is the number of individuals in each compartment.and X S = [ Q E I E I E I E I P P S S S S R R ] T = [0 0 0 0 0 0 0 0 0 b K d θ b K d θ b K d b K d b K d b K d T , is the set of disease free states. F ( x ), V − ( x ), and V + ( x ) are given in the following. F = b q I , V − = θ Q F = β S I / N , V − = ε E + d E N / K F = , V − = d I N / K , V + = θ Q + ε E F = β S I / N + β S I / N , V − = ε E + d E N / K F = , V − = γ I + µ I + d I N / K , V + = ε E F = β S I / N , V − = ε E + d E N / K F = , V − = d I N / K V + = ε E , F = β S I / N + f β S I / N , V − = d E N / K + ε E + β S I / N , F = , V − = γ I + µ I + d I N / K , V + = ε E F = , V − = b q I + θ P , V + = b N F = , V − = θ P , V + = b N F = , V − = β S I / N + d S N / K , V + = θ P F = , V − = d S N / K + β S I / N + β S I / N , V + = b N F = , V − = β S I / N + d S N / K , V + = θ P F = , V − = = β S I / N + f β S I / N , V + = b N + β S I / N + d S N / K F = , V − = d R N / K V + = γ I F = , V − = d R N / K , V + = γ I As it can been easily seen, the following five assumptions [16] are satisfied.(A1) if x (cid:62)
0, then F i , V + i , V − i (cid:62) i = , ..., x i =
0, then V − i =
0; in particular, if x ∈ X s , then V − i = i = , ..., F i = i >
9; there are no new infections in uninfected compartments.16A4) if x ∈ X s , then F i ( x ) = V + i ( x ) = i = , ..., F ( x ) is set to 0, then all eigenvalues of D f ( x ) have negative real parts.To construct the next generation matrix, we only consider infected and exposed compartments. The equations are trans-formed as follows. ddt Q E I E I E I E I = F − V = b q I β S I / N β S I / N + β S I / N β S I / N β S I / N + f β S I / N + β S I / N − θ Q d E N / K + ε E − θ Q + d I N / K − ε E d E N / K + ε E − ε E + d I N / K + γ I + µ I d E N / K + ε E d I N / K − ε E d E N / K + ε E d I N / K − ε E + γ I + µ I , The equation system is nonlinear; we linearize it, deriving the two Jacobian matrices. First, the partial derivative of F withrespect to each variable at the disease free equilibrium is as follows [16]. F = b q β S N β S N β S N β S N β S N f β S N β S N (74)Second, the partial derivative of V with respect to each variable at disease free equilibrium is as follows. V = θ b + ε − θ − ε b b + ε − ε b + γ + µ b + ε − ε b b + ε
00 0 0 0 0 0 0 − ε b + γ + µ (75)The inverse of matrix V is computed as follows. V − = θ b + ε b ε b ( b + ε ) 1 b b + ε ε ( b + ε )( b + γ + µ ) 1 b + γ + µ b + ε ε b ( b + ε ) 1 b b + ε
00 0 0 0 0 0 0 ε ( b + ε )( b + γ + µ ) 1 b + γ + µ Finally, the next generation matrix, which is the product FV − , is as follows.17 V − = A B A C D E F G H I J L M N P Q R S T where: A = q B = q ε b + ε C = ε β b K d ( b + ε )( b + γ + µ ) d b K D = β b K d ( b + γ + µ ) d b K E = β b K d b d K F = ε β b K d ( b + ε ) b d K G = β b K d b d K H = ε β b K d ( b + ε ) b d K I = β b K d b d K J = ε b K d β ( b + ε )( b + γ + µ ) d b K L = b K d β ( b + γ + µ ) d b K M = b K d β b d K N = ε β b K d ( b + ε ) b d K P = b d K β b d K Q = f ε b K d β ( b + ε )( b + γ + µ ) d b K R = f b K d β ( b + γ + µ ) d b K S = ε b K d β ( b + ε ) b d K T = b K d β b K d Recall that the reproduction number is the spectral radius of FV − , we compute the nine eigenvalues λ i of FV − to select the onewith maximum magnitude. 18 FV − − λ I | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) A − λ B A − λ C D − λ E F G − λ H I − λ J L − λ − λ M N P Q R S T − λ
00 0 0 0 0 0 0 0 − λ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = − λ [ λ − A λ − ( C F + J H ) λ + ( H A J + A C F − B C E )] (76) B C E = q ε ε β β b ( b + ε )( b + ε )( b + µ + γ ) (77) A C F = q ε ε β β b ( b + ε )( b + ε )( b + µ + γ ) (78)Because B C E = A C F as shown in Equation (78) and (77) , equation (76) can be rewritten as follows. − λ [ λ − A λ − ( C F + J H ) λ + H A J ] = | λ i | ( i = , , λ − A λ − ( C F + J H ) λ + H A J = λ − q λ − [ ε ε β β b ( b + ε )( b + ε )( b + µ + γ ) + ε ε β β b ( b + ε )( b + ε )( b + µ + γ ) ] λ + q ε ε β β b ( b + ε )( b + ε )( b + µ + γ ) = ff erent sets of parameters uniformly distributed within the rangein [18]. The histogram of the reproduction number is shown in Figure 7. From the histogram, we can see that R can be greateror smaller than 1. In particular, the mean is 1 .
17 and the maximum is 3 .
68, respectively. H i s t og r a m o f R Figure 7: Histogram of the reproduction number, the mean is 1 .
17, the maximum is 3 . pper and Lower Bound for R Although we are able to only obtain the exact expression of R numerically, we determine the lower bound and the upperbound of R in the following. | FV − − λ I | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) A − λ B A − λ C D − λ E F G − λ H I − λ J L − λ − λ M N P Q R S T − λ
00 0 0 0 0 0 0 0 − λ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = − λ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) A − λ B − λ C E F − λ H J − λ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = − λ |A − λ I | Where: A = A B C E F H J = q q ε b + ε ε β b K d ( b + ε )( b + γ + µ ) d b K β b K d b d K ε β b K d ( b + ε ) b d K ε β b K d ( b + ε ) b d K ε b K d β ( b + ε )( b + γ + µ ) d b K Matrix FV − has five zero roots. To find out the spectral radius of FV − , we only need to find out the spectral radius of matrix A which can be rewritten as follows. A = ε β b K d ( b + ε )( b + γ + µ ) d b K β b K d b d K ε β b K d ( b + ε ) b d K ε β b K d ( b + ε ) b d K ε b K d β ( b + ε )( b + γ + µ ) d b K + q q ε b + ε = B + CB is the first matrix and C is the second matrix. The matrix C has three multiple eigenvalues λ = λ = λ =
0. Since
Rank ( λ i I − C ) | λ = λ = λ = =
1, there are three linear independent eigenvectors corresponding to zero eigenvalue. Overall, matrix C has four linear independent eigenvectors. Therefore, matrix C can be diagonalized as follows. C = P q P − , where: P = − q ε b + ε q , P − = ε b + ε q , Substitute C in the expression of A = B + C . A = P ( P − B P + q ) P − = P ( X + Y ) P − where X = P − B P = ε ε β b K d ( b + ε )( b + ε )( b + γ + µ ) d b K
00 0 ε β b K d q ( b + ε )( b + γ + µ ) d b K β b K d b d K ε β b K d ( b + ε ) d b K ε β b K d ( b + ε )( b + γ + µ ) d b K , Y = q . Matrix Y is a 4 × Y i j = i (cid:44) j ) and 0 (cid:54) Y ii < q = max j v j j < ∞ , ( i = , , X is nonnegative. Therefore, ρ ( X ) (cid:54) ρ ( X + Y ) (cid:54) ρ ( X ) + q according to Theorem 1 . in [10]. The eigenvalues of matrix X + Y A because the two matrices are similar. Similarly, the eigenvalues of matrix X are the same as those ofmatrix B . ρ ( X + Y ) = ρ ( FV − ) = R ρ ( X ) = (cid:115) ε ( b + ε )( b + γ + µ ) (cid:104) ε β β b ( b + ε ) + ε β β b ( b + ε ) (cid:105) If we only count the horizontal transmission and denote the new F (resp. V ) by F H (resp. V H ) , F H and V H are as follows. F H = β S N β S N β S N β S N β S N β S N β S N , V H = d N K + ε − ε d N K d N K + ε − ε d N K + γ + µ d N K + ε − ε d N K d N K + ε
00 0 0 0 0 0 − ε d N K + γ + µ By calculation, R H = ρ ( F H V − H ) = (cid:115) ε ( b + ε )( b + γ + µ ) (cid:104) ε β β b ( b + ε ) + ε β β b ( b + ε ) (cid:105) = ρ ( X )Therefore, R H (cid:54) R (cid:54) R H + q (82)We denote R H and R H + q as R L and R U , respectively. Numerical Comparison among R , R L , and R U R H and R H + q are the lower bound and the upper bound of R , respectively. Therefore, R H + q < ⇒ R <
1, and R H > ⇒ R >
1. To verify that the derived bounds are tight, we perform extensive simulations using 5000 sets of parametersuniformly distributed within the range defined in [18]. In Figure 8(a) and 8(b), R U = R H + q vs. R and R L = R H vs. R areplotted, respectively. First, the di ff erence between the exact values and each bound is very small. In fact, the red and green pointslay very close to the line y = x . Additionally, in the case of the upper bound, the red points are just slightly above the line y = x ,while in the case of the lower bound, the green points are just slightly below the line y = x .21 R R U upper boundy=x (a) The reproduction number and its upper bound. R R lower boundy=x (b) The reproduction number and its lower bound.Figure 8: The reproduction number and its upper and lower bound iological Interpretation of Bounds for R The bounds for R , as given in inequalities (26), can be interpreted biologically as follows. The lower bound, R H , is thereproduction number for horizontal transmission because R H = ρ ( F H V − H ), where ρ ( F H V − H ) represents the spectral radius ofthe next generation matrix for horizontal transmission F H V − H . The upper bound is given by the sum of R H and a second termthat is only related to vertical transmission, i.e. from mothers to their o ff spring in the Aedes mosquito population. R H includes Aedes -livestock interaction and
Culex -livestock interaction. More specifically, we define the reproduction num-ber due to the interaction between
Aedes and livestock represented by R H ( A − L )0 as R H ( A − L )0 = (cid:115) ε ( b + ε )( b + γ + µ ) ε β β b ( b + ε ) R H ( A − L )0 can be rewritten as follows. R H ( A − L )0 = (cid:115)(cid:104) β d N ∗ K ε ( d N ∗ K + ε ) (cid:105)(cid:104) β ( d N ∗ K + γ + µ ) ε ( d N ∗ K + ε ) (cid:105) (83)because b = d N ∗ K b = d N ∗ K where: N ∗ = the total number of Aedes mosquitoes at disease free equilibrium. N ∗ = the total number of livestock at disease free equilibrium. R H ( A − L )0 consists of the product of four terms. Each infected Aedes mosquito can infect β d N ∗ K susceptible livestock throughoutits lifetime. Similarly, each infected livestock can infect β d N ∗ K + γ + µ susceptible Aedes mosquitoes during its lifetime. The prob-ability of
Aedes mosquitoes and livestock surviving through the incubation period to the point where they become infectiousis ε d N ∗ K + ε and ε d N ∗ K + ε , respectively. Therefore, R H ( A − L )0 is the geometric mean of the average number of secondary livestockinfections produced by one Aedes mosquito vector in the first square bracket in (83), and the average number of secondary
Aedes mosquito vector infections produced by one livestock host in the second square bracket in (83).Similarly, we define the reproduction number due to the interaction between
Culex and livestock represented by R H ( C − L )0 as R H ( C − L )0 = (cid:115) ε ( b + ε )( b + γ + µ ) ε β β b ( b + ε )We can rewrite R H ( C − L )0 as follows. R H ( C − L )0 = (cid:115)(cid:104) β d N ∗ K ε ( d N ∗ K + ε ) (cid:105)(cid:104) β ( d N ∗ K + γ + µ ) ε ( d N ∗ K + ε ) (cid:105) (84)because b = d N ∗ K where: N ∗ = the total number of Culex mosquitoes at disease free equilibrium.23 H ( C − L )0 also consists of the product of four terms. Each infected Culex mosquito can infect β d N ∗ K susceptible livestockthroughout its lifetime. Similarly, each infected livestock can infect β d N ∗ K + γ + µ susceptible Culex mosquitoes. The probability of
Culex mosquitoes surviving through the incubation period to the point where they become infectious is ε d N ∗ K + ε . Similarly, theprobability of livestock surviving through the incubation period to the point where they become infectious is ε d N ∗ K + ε . Therefore, R H ( C − L )0 is the geometric mean of the average number of the secondary livestock infections produced by one Culex mosquitovector in the first square bracket in (84), and the average number of secondary
Culex mosquito vector infections produced byone livestock in the second square bracket in (84).The expression (27) for R H , can be rewritten as R H = (cid:113) ( R H ( A − L )0 ) + ( R H ( C − L )0 ) , where the dependence of R H on R H ( A − L )0 and R H ( C − L )0 is shown in Figure 9. The square root is due to the vector-host-vector viral transmission path [18 ? , 33]. Obviously,the horizontal reproduction number increases with the increase of each of the four terms in R H ( A − L )0 and R H ( C − L )0 . H R O )(0 LAH R - )(0 LCH R - Figure 9: The interpretation of R H References [1] Anyamba, A., Chretien, J.P., Small, J., Tucker, C.J., Formenty, P.B., Richardson, J.H., Britch, S.C., Schnabel, D.C., Erickson, R.L., Linthicum, K.J., 2009.Prediction of a Rift Valley fever outbreak. Proceedings of the National Academy of Sciences 106, 955–959.[2] Anyamba, A., Chretien, J.P., Small, J., Tucker, C.J., Linthicum, K.J., 2006. Developing global climate anomalies suggest potential disease risks for2006-2007. International Journal of Health Geographics 5, 60.[3] Anyamba, A., Linthicum, K.J., Mahoney, R., Tucker, C.J., 2002. Mapping potential risk of Rift Valley fever outbreaks in african savannas using vegetationindex time series data. Photogrammetric Engineering and Remote Sensing 68, 137–145.[4] Anyamba, A., Linthicum, K.J., Tucker, C.J., 2001. Climate-disease connections: Rift Valley fever in kenya. Cadernos de saude publica / Ministerio daSaude, Fundacao Oswaldo Cruz, Escola Nacional de Saude Publica 17 Suppl, 133–140.[5] Balcan, D., Colizza, V., Goncalves, B., Hu, H., Ramasco, J.J., Vespignani, A., 2009. Multiscale mobility networks and the spatial spreading of infectiousdiseases. Proceedings of the National Academy of Sciences of the United States of America 106, 21484–21489.[6] Centers for Disease Control and Prevention (CDC), 2007. Rift Valley fever outbreak–kenya, november 2006-january 2007. Morbidity and MortalityWeekly Report 56, 73–76.[7] Chevalier, V., Lancelot, R., Thiongane, Y., Sall, B., Diaite, A., Mondet, B., 2005. Rift Valley fever in small ruminants, senegal, 2003. Emerging IinfectiousDiseases 11, 1693–1700.[8] Chowdhury, S.R., Scoglio, C., Hsu, W., 2010. Simulative modeling to control the foot and mouth disease epidemic. Procedia Computer Science 1,2261–2270.[9] Clements, A.C., Pfei ff er, D.U., Martin, V., Pittliglio, C., Best, N., Thiongane, Y., 2007. Spatial risk assessment of rift Valley fever in senegal. Vector Borneand Zoonotic Diseases (Larchmont, N.Y.) 7, 203–216.[10] Cohen, J., 1979. Random evolutions and the spectral radius of a non-negative matrix. Math. Proc. Camb. Phil. Soc. 86, 345–350.[11] Davies, F.G., Martin, V., 2006. Recognizing Rift Valley fever. Veterinaria Italiana 42, 31–53.[12] Department for Environment Food and Rural A ff airs, 2010. Rift Valley fever. . Accessed Nov 10, 2010.[13] Department of Agriculture Forestry and Fisheries of Republic of South Africa, 2010. Livestock number 96 to date. . Accessed Sep 26, 2010.[14] Diekmann, O., Heesterbeek, J.A.P., 2000. Mathematical Epidemiology of Infectious Diseases: Model Building, Analysis and Interpretation (Wiley Seriesin Mathematical & Computational Biology). Chapter 5, Wiley.[15] Disease BioPortal, 2010. http://fmdbioportal.ucdavis.edu . Accessed Nov 23, 2010.[16] van den Driessche, P., Watmough, J., 2002. Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission.Mathematical biosciences 180, 29–48.
17] Florida Department of Health, 2010. Rift Valley fever. . Accessed Nov 30, 2010.[18] Ga ff , H.D., Hartley, D.M., Leahy, N.P., 2007. An epidemiological model of Rift Valley fever. Electron. J. Di ff . Eqns 2007, 1–12.[19] Gong, H., Degaetano, A.T., Harrington, L.C., 2010. Climate-based models for west nile culex mosquito vectors in the northeastern us. InternationalJournal of Biometeorology 55, 435–446.[20] Gubler, D.J., 2002. The global emergence / resurgence of arboviral diseases as public health problems. Archives of Medical Research 33, 330–342.[21] Kasari, T.R., Carr, D.A., Lynn, T.V., Weaver, J.T., 2008. Evaluation of pathways for release of Rift Valley fever virus into domestic ruminant livestock,ruminant wildlife, and human populations in the continental united states. Journal of the American Veterinary Medical Association 232, 514–529.[22] Keeling, J., Rohani, P., 2008. Modeling infectious disease in humans and animals. Princeton University Press.[23] Kim, B.N., Gordillo, L.F., Kim, Y., 2010. A model for the transmission dynamics of orientia tsutsugamushi among its natural reservoirs. Journal oftheoretical biology 266, 154–161.[24] Konrad, S.K., s. N. Miller, Reeves, W.K., 2010. A spatially explicit degree-day model of Rift Valley fever transmission risk in the continental united states.GeoJournal .[25] Linacre, E.T., 1977. A simple formula for estimating evaporation rates in various climates using temperature data alone. Agricultural Meteorology 18,409–424.[26] Linthicum, K.J., Anyamba, A., Britch, S.C., Chretien, J.P., Erickson, R.L., Small, J., Tucker, C.J., Bennett, K.E., Mayer, R.T., Schmidtmann, E.T.,Andreadis, T.G., Anderson, J.F., Wilson, W.C., Freier, J.E., James, A.M., Miller, R.S., Drolet, B.S., Miller, S.N., Tedrow, C.A., Bailey, C.L., Strickman,D.A., Barnard, D.R., Clark, G.G., Zou, L., 2007. A Rift Valley fever risk surveillance system for Africa using remotely sensed data: potential for use onother continents. Veterinaria Italiana 43, 663–674.[27] Linthicum, K.J., Davies, F.G., Kairo, A., Bailey, C.L., 1985. Rift Valley fever virus (family bunyaviridae, genus phlebovirus). isolations from dipteracollected during an inter-epizootic period in kenya. The Journal of Hygiene 95, 197–209.[28] Martin, V., Chevalier, V., Ceccato, P., Anyamba, A., Simone, L.D., Lubroth, J., de La Rocque, S., Domenech, J., 2008. The impact of climate change onthe epidemiology and control of Rift Valley fever. Revue Scientifique et Technique (International O ffi ce of Epizootics) 27, 413–426.[29] Mpeshe, S.C., Haario, H., Tchuenche, J.M., 2011. A mathematical model of rift valley fever with human host. Acta Biotheoretica 59, 231–250.[30] National Climatic Data center, 2010. NOAA Satellite and Information Service. . Accessed Nov 22, 2010.[31] National Institute for Communicable Diseases, 2010. Interim report on the Rift Valley fever (RVF) outbreak in south africa. . Accessed Nov 23, 2010.[32] Newton, E., Reiter, P., 1992. A model of the transmission of dengue fever with an evaluation of the impact of ultra-low volume (ulv) insecticide applicationson dengue epidemics. The American Journal of Tropical Medicine and Hygiene 47, 709–720.[33] O. Diekmann, H.H., Metz, H., 1995. The legacy of Kermack and McKendrick. In D. Mollison, editor, Epidemic Models: Their Structure and Relation toData . pages 95–115. Cambridge University Press, Cambridge, UK.[34] Olivier.C.G, 2004. An Analysis of the South African Beef Supply Chain: from farm to folk. Master’s thesis. University of Johannesburg. South Africa.[35] Sellers, R.F., Pedgley, D.E., Tucker, M.R., 1982. Rift Valley fever, egypt 1977: disease spread by windborne insect vectors? The Veterinary record 110,73–77.[36] South Africa Department of Health, 2010. Press releases 2010. . Accessed Nov 30, 2010.[37] Statistics South Africa, 2010a. Agricultural Census (Census of Commercial Agriculture), 2007. . Accessed Nov 22, 2010.[38] Statistics South Africa, 2010b. Domestic tourism survey 2009. . Accessed Nov 23, 2010.[39] Statistics South Africa, 2010c. Mid-year population estimates. . Ac-cessed Nov 21, 2010.[40] Swanson, J.C., Morrow-Tesch, J., 2001. Cattle transport: Historical, research, and future perspectives. Journal of Animal Science 79, E102–E109.[41] Weather Underground, 2010. . Accessed Nov 20, 2010.[42] Woods, C.W., Karpati, A.M., Grein, T., McCarthy, N., Gaturuku, P., Muchiri, E., Dunster, L., Henderson, A., Khan, A.S., Swanepoel, R., Bonmarin, I.,Martin, L., Mann, P., Smoak, B.L., Ryan, M., Ksiazek, T.G., Arthur, R.R., Ndikuyeze, A., Agata, N.N., Peters, C.J., Force, W.H.O.H.F.T., 2002. Anoutbreak of Rift Valley fever in northeastern kenya, 1997-98. Emerging Infectious Diseases 8, 138–144.[43] World Animal Health Information Database, 2010. Summary of immediate notifications and follow-ups - 2010. . Accessed Oct 14, 2010.[44] World Health Organization, 2010. Rift Valley fever. . Accessed Nov 22, 2010.[45] Zeller, H.G., Fontenille, D., Traore-Lamizana, M., Thiongane, Y., Digoutte, J.P., 1997. Enzootic activity of Rift Valley fever virus in senegal. TheAmerican Journal of Tropical Medicine and Hygiene 56, 265–272.. Accessed Nov 22, 2010.[45] Zeller, H.G., Fontenille, D., Traore-Lamizana, M., Thiongane, Y., Digoutte, J.P., 1997. Enzootic activity of Rift Valley fever virus in senegal. TheAmerican Journal of Tropical Medicine and Hygiene 56, 265–272.