Data-Driven Modeling to Assess Receptivity for Rift Valley Fever Virus

Rift Valley Fever virus (RVFV) is an enzootic virus that causes extensive morbidity and mortality in domestic ruminants in Africa, and it has shown the potential to invade other areas such as the Arabian Peninsula. Here, we develop methods for linking mathematical models to real-world data that could be used for continent-scale risk assessment given adequate data on local host and vector populations. We have applied the methods to a well-studied agricultural region of California with 1 million dairy cattle, abundant and competent mosquito vectors, and a permissive climate that has enabled consistent transmission of West Nile virus and historically other arboviruses. Our results suggest that RVFV outbreaks could occur from February–November, but would progress slowly during winter–early spring or early fall and be limited spatially to areas with early increases in vector abundance. Risk was greatest in summer, when the areas at risk broadened to include most of the dairy farms in the study region, indicating the potential for considerable economic losses if an introduction were to occur. To assess the threat that RVFV poses to North America, including what-if scenarios for introduction and control strategies, models such as this one should be an integral part of the process; however, modeling must be paralleled by efforts to address the numerous remaining gaps in data and knowledge for this system.


Introduction
Rift Valley fever virus (RVFV; viral family Bunyaviridae, genus Phlebovirus) is a pathogen that causes febrile illness in domestic ruminants (sheep, cattle, and goats) and humans throughout Africa and parts of the Arabian Peninsula [1][2][3] that may be transmitted by several genera of mosquitoes [3][4][5][6]. Outbreaks often result in heavy economic costs through loss of livestock, especially when associated with an incursion into a new area [7,8]. Although never detected in the Western Hemisphere, RVFV is a threat to human and livestock health in North America and is included on select agent lists of the U.S. Department of Health and Human Services and the U.S. Department of Agriculture [9].
Mosquito species found to be vectors of RVFV with varying degrees of efficiency in laboratory settings [6,10] are known to be present throughout much of the U.S. [11], but other aspects of potential transmission cycles remain inadequately studied. To properly assess and mitigate the risk posed by a RVFV invasion, methods are needed to identify areas that are most likely to support transmission, the time periods when transmission is expected to pose a risk, and whether an introduced virus could become established. To date, such questions have been addressed by only a few analytic methods, including expert elicitation [12], basic GIS overlays of humans and vectors with a hypothetical host [13], and pathways analysis [14][15][16].
Process-based mathematical models provide a useful platform to coalesce disparate data, make logical assumptions concerning data gaps, and evaluate a range of potential scenarios. Gaff et al. developed a dynamical model for RVFV [17] that included livestock hosts and two genera of mosquitoes, Aedes and Culex, that respectively were or were not capable vertically transmitting RVFV. This model's structure has been extended in several important ways to 1) accommodate spatial structure through host or vector movements [18,19], 2) assess potential control methods [20,21], or 3) include humans [22] or asymptomatic livestock hosts [23]. These models have resulted in important advances in modeling RVFV, but their application is limited by the lack of appropriate data to inform parameters, many of which have been recycled between models, defined arbitrarily, or borrowed from literature on other arboviruses that may not apply for RVFV.
In the current study, we apply a unique and generalizable approach that links real-world data with the mathematical models, utilizing broadly available national-scale data where possible. To illustrate the methodology, we consider results for the southern Central Valley in California, an area with large, well-documented host and vector populations. We present two model-derived transmission metrics that quantify expectations for typical (R 0 , [24]) and maximal (E 0 , [25]) transmission from an initial diseasefree state, and we map these metrics according to the spatial pattern of vector abundance associated with various land uses. We also show how these metrics change in time as a function of temperature, thereby enabling the assessment of seasonal transmission risk. We also highlight several critical data gaps that must be addressed.

Study area
California's Great Central Valley extends for 700 km north to south through the center of the state and is home to extensive and varied agricultural lands that include irrigated crops, livestock operations, natural or restored wetlands, and urbanized areas. In this study, we considered the southern half of the Central Valley ( Figure 1), which contains very high densities of livestock (primarily dairy cows, with w1 million cattle within the study area) interspersed with managed wetlands and multicrop agriculture that can produce large populations of competent vectors (e.g., Culex tarsalis). Deer were rare on the valley floor and generally were restricted to surrounding higher-elevation foothills and mountains, and sheep typically were moved into the valley for grazing only during the cooler months of the year when transmission was expected to be minimal. The area is likely to be climatologically permissive for RVFV transmission as this is the warmest part of the valley and supports consistently high transmission of West Nile virus [26] and, previously, other arboviruses [27]. For the model, appropriate spatial dimensions were needed for patches that would represent the heterogeneity in land cover and host and vector densities at a fine enough scale that populations could be assumed to be well-mixed, given described ranges of vector movement [28]. To achieve this, we defined a uniform grid of 5|5 km squares (25 km 2 ) that covered the study area, and all model input variables were scaled to this grid. All model outputs were calculated by grid cell for each day of the year.

Data
Temperature. Daily maximum and minimum temperatures were acquired from the National Aeronautics and Space Administrations Terrestrial Observation and Prediction System (TOPS; http:// ecocast.arc.nasa.gov) [29], which uses weather and ecosystem models to combine ground-based and remotely sensed inputs to generate multiple measures of environmental conditions. For this study, the TOPS surfaces (1 km 2 resolution) for daily mean temperatures were averaged for each day of the year for the most recent 10-year period (2002-2011) and spatially within each 25 km 2 grid cell so that temperatures represented the typical pattern for each spatial location, summarized in Figure 2.
Vectors. Mosquitoes in the genera Aedes and Culex are important vectors of RVFV in enzootic areas, and we focused on two species that are both abundant within our study area [26,30] and likely to be capable of transmitting RVFV, albeit to differing degrees [6]. Aedes mosquitoes can be infected with RVFV either vertically from infected females [31] or horizontally via a blood meal from an infectious host. We focused on Aedes melanimon because it has been an important vector of vertically maintained California encephalitis virus [32], frequently feeds on mammals [33], and is likely to be a low-competent horizontal vector of RVFV based on results for the closely related species, Aedes dorsalis [6]. Culex vectors are able to transmit RVFV horizontally, but not vertically to their offspring, and here we considered Cx. tarsalis, which is a principal vector of several encephalitis viruses [34,35], feeds opportunistically on both birds and mammals [36,37], and is the most competent laboratory vector of RVFV studied in North America [6].
Annual patterns for the relative abundance of these vectors were assigned to each of several broad land use classes representing a generalization of the narrowly defined single-crop classes in the most recent USDA Cropland Data Layer for 2011 ( Figure  1; http://www.nass.usda.gov/research/Cropland/ SARS1a.htm). Our goal was to represent typical variation in daily abundance for each of the two mosquito genera in each land use class based on data from dry ice-baited CDC-style traps [38] operated by vector control agencies within the study area. Averaging trap counts spatially or over several years would not achieve this goal, especially for Aedes populations that hatch synchronously in response to flooding of eggs, resulting in an abundance spike when the cohort emerges. Such sharp peaks occur at similar but slightly different dates each year, with the result that a multi-year average would smooth the annual spikes into a rounded peak that is not representative of any year. To avoid this problem, a small number (3 to 5) of representative CO 2 -baited trap sites were identified within each land cover class, and their time series of trap counts were used to define a ''consensus'' time series for each class that captured the key features of its annual pattern (namely, seasonally varying rates of increase and decrease, and abundance minima and maxima; Figure 3). To apply these patterns spatially, each spatial grid cell was characterized by its dominant (i.e., most common) land cover class, and the relevant abundance patterns for Aedes and Culex were applied, multiplied by the gonotrophic period to scale the fraction of vectors that would be host-seeking on a given night to the total population size and by each grid cell's respective total host abundance (cows+birds) as defined below, based on the assumption that a trap represents one host (also see Mathematical Appendix Table S2 in Text S1). The scale factor for the gonotrophic period was used because traps represent only the fraction of all females that seeks hosts on a given night, equal to 1=(gonotrophic period) on average, but the model requires an

Author Summary
Rift Valley fever virus is a pathogen enzootic to sub-Saharan Africa, with epidemic transmission occurring sporadically between mosquitoes and mammals, notably livestock. The virus is regarded as a global threat to agriculture and human health because it has proven capable of expanding its range into western and northern Africa, Madagascar, and the Arabian Peninsula, and a recent study has shown that mosquitoes in North America are capable of transmitting the virus. Here, we used a set of mathematical equations to formulate a logical representation of potential transmission mechanisms, and we informed the model with real-world data and generalizable methods to define spatial and temporal variation in mosquito and host abundance. We applied these methods in California's warm, agricultural Central Valley, an area with a history of mosquito-borne virus transmission and a hub of California's dairy industry. Model-derived transmission estimates indicated broad potential for transient epidemics that could result in economic losses in livestock in all but the coldest winter months, but the greatest risk for intense, sustained transmission occurred during the summer when both vector abundance and temperatures were highest. We also highlight critical gaps in the data available to inform models for Rift Valley fever virus.  Livestock. RVFV antibodies and some viral isolates have been observed in a variety of African mammals, but viremias sufficient to infect biting vectors have been reported only in cattle, sheep, and goats [1], although high viremias may also occur in humans [39]. Our model includes two host populations, with the first represented by dairy cows considered to be competent for transmitting RVFV, and the second represented by birds considered to be incompetent ''sink'' hosts. Within our study area, dairy cows were the most abundant ungulates, and we calculated the total number of cows within each grid cell using data on dairy sizes obtained from the California Environmental Protection Agency's Central Valley Regional Water Quality Control Board.
Birds. Cx. tarsalis also feeds opportunistically upon birds [36,37], and we expected birds to be likely alternate hosts to represent the incompetent class. We acknowledge that Ae. melanimon also feeds frequently on hares and rabbits (Order Lagomorpha [33]), but data on their abundance within our study area were not available. To quantify the abundance of birds, we considered seven species of birds that are abundant near dairies and could divert bites from cattle: Red-Winged Blackbirds (Agelaius phoeniceus), Brown-Headed Cowbirds (Molothrus ater), Brewer's Blackbirds (Euphagus cyanocephalus), House Finches (Haemorhous mexicanus), House Sparrows (Passer domesticus), European Starlings (Sturnus vulgaris), and Rock Pigeons (Columba livia). Gridded data from the USGS Breeding Bird Survey (BBS; 21.5 km resolution) were downscaled to our 5-km grid by calculating a weighted average of the per-BBS route bird abundance for the species above within a 10-km buffer around each 5-km grid cell. Next, we rescaled per-route bird abundance to an estimate of the total birds within each grid cell by multiplying the per-route values by the ratio (&4) of estimated areal coverage of BBS routes to the area of the grid cell. BBS routes consist of 50 stops, spaced 800 m apart, with observers attempting to detect all birds within a 0.4-km radius at each stop; assuming they effectively observed half that radius (see discussion of detection in [40]), the area observed would be 6.28 km 2 .

Mathematical modeling
Process-based, dynamical mathematical models of virus transmission are built from knowledge of the interactions among virus, hosts, and vectors. In the case of RVFV in North America, such issues are uncertain. In the current study we extend previous work [17,41] to construct a mathematical model of RVFV (Figure 4), with the following assumptions regarding the anticipated epidemiology of RVFV in California.

1.
Hosts. We consider one prototypical competent host, which corresponds to domestic ungulates. Incubation periods and peak viremia in cattle, goat, and sheep species are similar, so the generic host could be any of these animals. Competent hosts become infected when fed upon by infectious vectors. Such hosts may then die from RVFV infection or recover, whereupon they have lifelong immunity. Incompetent or dead end hosts may be either mammals or birds, and are fed upon by vectors but do not become infectious and therefore serve as a sink for the virus and dampen transmission. 2. Vectors. We included two types of vectors, with the first (typified by Ae. melanimon) capable of vertical transmission to offspring and low-level competence for horizontal transmission, and the second (typified by Cx. tarsalis) being a highly competent horizontal vector that does not transmit the virus vertically. Vertical transmission of RVFV has been observed in fieldcollected larvae [31], but laboratory data from related viruses show a complex picture of the efficiency of vertical transmission. The literature suggests that vertical transmission is a lowfrequency event in nature on average. Therefore, we included vertical transmission in the model at low rates of occurrence. 3. Vector host selection. Vectors bite either RVFV competent or incompetent hosts in proportion to total host population size; i.e., we assume that there is no mosquito feeding preference.
Particularly for Ae. melanimon that feeds primarily on mammals, this assumption may seem questionable, but the very high densities of cattle in our study area ensured that they were the predominant bloodmeal source even with opportunistic feeding (cattle represented the majority of hosts in 94% of grid cells and w80% of hosts in 76% of grid cells).  Temperature dependence of the vectors' extrinsic incubation rate (i.e., the inverse of the extrinsic incubation period, EIP) and the gonotrophic cycle length (gonotrophic period, GP) were modeled based on published data as follows. We digitally extracted data points for temperatures at which a median EIP could be estimated from Figure 3 {1 . The relationship between EIP and temperature also has been studied for Culex mosquitoes [42,44], but heterogeneity among experiments and a paucity of comparable data points precluded construction of an analogous EIP model for this genus. Therefore, we modeled extrinsic incubation in Culex and Aedes by the same function. The gonotrophic period (i.e., the number of days between bloodmeals) was modeled as GP~2z({0:066z 0:018|temperature) {1 , using a published linear regression equation for the ovarian maturation rate [45] plus 2 days for oviposition and locating a bloodmeal host. The environmental carrying capacity could not be explicitly measured, and for both vectors, it was approximated daily using the vector abundance for the following day based on the typical abundance time series described above. This resulted in the desired inflation and deflation of the density-dependent birth rates in proportion to the rate of population growth or shrinkage, respectively, with a corresponding inverse impact on death rates.
We implemented the full model using a set of ordinary differential equations (mathematical details appear in the appendix). Using the methods described in van den Driessche and Watmough [46], we derived an expression for the basic reproduction ratio, R 0 , which represents the average number of secondary infections that arise from a single infectious individual (vector or host) introduced into a completely susceptible population [24,47], so that when R 0 v1, there are insufficient new cases per case for propagation and the pathogen cannot persist in the population. When R 0 §1, the pathogen is efficiently transmitted and becomes enzootic; elevated R 0 values indicate that transmission is more intense and that stochastic fadeout of the pathogen is less likely. For complex models of vectorborne infections, it has been demonstrated that outbreaks are possible for R 0 v1 under certain circumstances [48,49]. Because the model incorporates both vertical and horizontal transmission, R 0 was written as the sum of the R 0 values for each mode of transmission determined separately, R 0~R0(V ) zR 0(H) [17,18,50]. Details of the In addition to R 0 , we computed a recently described, metric, E 0 , that quantifies the reactivity, or epidemicity, of the system [25]. E 0 represents the maximum number of new infections produced by an infective individual at a disease free equilibrium, and like R 0 , epidemicity is a threshold quantity; when E 0 w1, transient epidemics (i.e., outbreaks that may eventually fade out) are possible regardless of the average system behavior predicted by R 0 . When E 0 v1, transmission, even for brief time periods, is not expected. Evaluating E 0 from our model allows us to investigate the potential for RVFV outbreaks in areas and times when R 0 suggests that efficient transmission is not possible. R 0 and E 0 are both functions of the model parameters shown in Mathematical Appendix Table S1 in Text S1.
Stochastic sampling from biologically relevant ranges of parameters was used to assess the sensitivity of R 0 and E 0 to the model parameters. The ranges for each parameter are presented in Mathematical Appendix Table S3 in Text S1. K A and K C , the vector carrying capacities, were computed from observed data as described above. Likewise, the EIP and vector GP values were functions of temperature. We assumed a uniform distribution for each parameter across ranges shown in Mathematical Appendix Table S3 in Text S1. The ranges of all the other parameters are from the references shown in Mathematical Appendix Table S1 in Text S1. Our model includes V~17 uncertain variables, so N~200 sets of sampled parameter values were generated by Latin hypercube sampling following the suggestion of Matala [51] that an N such that N=V w10 should suffice for the number of stochastic samples of complete parameter sets. Partial rank correlation coefficients (PRCC) were computed across ranges of parameters described in Mathematical Appendix Table S3 in Text S1 to assess the significance of each parameter with respect to R 0 and E 0 . Spatial analysis of temperatures, land cover, and host and vector abundance was carried out using R version 2.15 [52], ArcGIS 10.0 (ESRI, Redlands, CA, USA), and PostgreSQL 9.0 (http://www. postgresql.org) databases with added spatial capabilities of PostGIS (http://postgis.refractions.net). All code for mathematical modeling was written in R version 2.15 [52].

Potential for establishment
The seasonal patterns for the basic reproductive ratio for RVFV (R 0 ; Figure 5) indicated that the risk for sustained transmission increased rapidly by May, with R 0 exceeding 1 in the areas with both cattle and field crops ( Figure 5). Initially, these areas at risk consisted of a small number of grid cells in the center of the study area and a single cell near the southernmost end of the valley that could serve as early foci for transmission, and these areas remained at higher risk than other areas through the summer. Risk was greatest overall from late June-September, when a much broader area was at risk for sustained transmission that included grains and field crops and covered Tulare County, the core of California's dairy industry. In all areas, R 0 values v1 from October-April indicate that introductions from late fall-early spring would be unlikely to become established and that persistence of RVFV through winter may depend on mechanisms for long-term maintenance between epidemics (e.g., vertical transmission in vectors).

Potential for outbreaks
Our estimates of epidemicity for RVFV, E 0 , were much higher than the average expectations of R 0 (max = 95.3 and 3.2, respectively) and indicated that transient outbreaks could occur over a broader spatio-temporal window than that circumscribed by R 0 alone (Figure 6), although the relative seasonal patterns for the two metrics were strongly correlated (r~0:85,pv0:0001). E 0 values w1 indicate that transmission was possible from February-November in agricultural areas, although the number of cases probably would remain small for introductions during February-April or November. The highest transmission potential ( Figure 6) occurred in places and times where cattle were present and vectors reached high densities. Specifically, E 0 values were greatest in areas dominated by field crops and grains, which generally had the highest combined concentrations of Culex vectors and cattle in the study area, although the latter had little impact on epidemicity. Risk of RVFV transmission in urban areas was somewhat lower, with the highest level expected during spring, associated with the typical peak in Culex abundance at that time, followed by a slow decline in both E 0 and vector abundance through the end of the summer. Orchards and vineyards (trees/grapes) and grasslands had low risk for epidemics due to their low vector abundance, and fallow or barren habitats had no risk for outbreaks, even when cattle were present. In all areas, E 0 was closely linked to the abundance and carrying capacity of vectors, and E 0 was greatest during spring and summer when the abundance of vectors, primarily Culex, was highest. For the scenarios under study, other parameters had little impact on E 0 , resulting in negligible variation around the average seasonal pattern for each land use class ( Figure 6).

Wetlands
Seasonally-flooded wetlands were not included in our analysis because the relatively small fraction of the total study area that they occupied (Figure 1) did not include dairy cattle. However, these areas deserve special consideration because they were the only areas where both Aedes and Culex reached high abundance (Figure 3), making them a reasonable analogue to the dambo habitats where RVFV is enzootic in east Africa, although the timing of their flooding is linked to human water management rather than rainfall. If dairy cattle or other competent hosts were present, both horizontal and vertical transmission would be possible, presenting strong potential for both transient epidemics and establishment. To evaluate this possibility, we calculated E 0 and R 0 for the wetlands in the study area (Figure 1), with the addition of 1,000 cattle to each grid cell to simulate the effect of a moderate-sized dairy adjacent to the wetlands (median dairy size for the study area~1,140 cattle). Transmission potential was higher than for other land uses and reflected the abundance of Culex vectors as before. There were two annual peaks of 186 on June 14 and 489 on September 30 for E 0 and 2.84 on June 27 and 4.37 on September 27 for R 0 .

Sensitivity analysis
Sensitivity analysis revealed that both R 0 and E 0 were particularly responsive to the abundance (N A ,N C ), carrying capacity (K A ,K C ), and vector competence (r LA ,r LC ) of mosquito vectors. All of these relationships were positive, with the exception that R 0 had a negative correlation with mosquito abundance due to the density dependence of birth (and correspondingly, death) rates such that higher numbers of vectors resulted in smaller values for the K=N ratio, which also limited R 0 . R 0 was much more sensitive than E 0 to variation in temperature that affected RVFV extrinsic incubation and vector biting rates, which explained its broader range of values within each land use class ( Figure 5, lower panel). The abundance of hosts had little effect on either transmission metric. Complete results for the sensitivity analysis are presented in Mathematical Appendix Table S3 in Text S1.

Discussion
This study extends previous models of RVFV dynamics [17][18][19][20][21][22][23]53,54] by developing methods for linking our model and potentially others to real landscapes. We have demonstrated these methods in a well-studied and potentially RVFV-receptive region of California with high densities of dairy cattle and mosquitoes and a history of arbovirus transmission. Where possible, we have utilized readily available continental-scale data sources, and as a result, our methods could be employed to broadly assess risk for RVFV elsewhere in North America, given adequate knowledge of local patterns in vector and host population dynamics, although realistically addressing this last point is not a trivial hurdle for any broad-scale risk assessment.
Our model modifies the parameterization of several earlier RVFV models [17][18][19][20]23] that had the potential to downplay the role of vertical transmission in Aedes. In these models, the infected subadult stage (eggs or simply ''aquatic stage'') arose by applying the probability of vertical transmission (denoted q A herein) to infectious adult females (I A ). This would be appropriate, except that it implies that salivary transmission of RVFV by Aedes females is a necessary condition for vertical transmission to offspring. These processes occur by parallel infection pathways, although both require initial passage through the midgut barrier and a subsequent incubation period. The model we have presented continues to define infected eggs (Q A ) as originating from infectious females (I A ), which allows for the reasonable assumption that the incubation period for horizontal and vertical transmission are equal. However, we rescaled the parameter q A by the vector competence parameter r LA to match our understanding of the natural definition of q A , equal to the per capita probability of infection in offspring given that a female mosquito feeds on an infectious host (See Mathematical Appendix, Equation 2). Despite the conceptual aspects, previous models have assumed larger values for this vertical transmission parameter (most commonly 0.05) compared to the 0.001 we have chosen here based on the rarity of isolations from field-collected immature mosquitoes [31] and null findings in the laboratory [5,[55][56][57], which would negate some of the potential downward bias in model outcomes.
Previous dynamical models for RVFV have evaluated transmission outcomes for single parameter sets, generally combined with sensitivity analysis, but have not focused on realistic spatiotemporal variation in transmission risk for RVFV. One recent exception included two sets of parameters for ''dry'' vs. ''wet'' seasons [23], resulting in R 0 values of 0.80 and 2.28, respectively, which were within the range of our estimates (0. 23-3.24). The range for E 0 was much narrower (1.84-10.57 compared to our 0. 18-95.27), which is probably due to our broader range of parameters, particularly the data-driven seasonal variation and inequality of the abundance and carrying capacity of vectors in our model. Unlike in RVFV-enzootic areas of sub-Saharan Africa [53,58], California's rainfall occurs during the coolest months of the year and does not directly drive the population dynamics of most mosquitoes, so the deliberate distribution of surface water by humans (e.g., in irrigated agriculture and wetlands) creates and maintains aquatic mosquito habitats, and therefore transmission is expected to be greatest during the warm, ''dry'' season when both vector abundance and temperatures are high. R 0 values from single-patch [17] and metapopulation models for RVFV [19] (0:04{3:74 and &0{3:68, respectively) also agreed well with our estimates, although the ranges in these studies were calculated from parameter sets randomly drawn from uniform distributions rather than ''real-world'' scenarios.
We focused on two important transmission metrics, R 0 and E 0 , to convey the average and maximal expectations for transmission in initially disease-free, immunologically nave populations. Such a scenario is appropriate for the Western Hemisphere, which has no history of RVFV transmission, and these metrics were useful for defining the spatio-temporal range over which an introduced virus could cause outbreaks and the relative variation in risk over a typical year, which are driven primarily by seasonal patterns in vector abundance and temperature.
Our model-based estimates suggest that the risk for an RVFV outbreak -if it were to be introduced -should be a concern for California's dairy industry during all but the coldest winter months of December and January and is highest during the warm spring and summer seasons in the Central Valley. Epidemic risk increases rapidly once the abundance of mosquitoes begins to increase, typically in March-April, and E 0 has shown that the potential for RVFV circulation begins earlier and lasts longer than would be suggested by an analysis informed by R 0 alone. The theoretical threshold of 1.0 for R 0 should be viewed with caution as an absolute limit for transmission, especially for vectorborne diseases [48], and its values depend on the models and methods chosen for computation [59]. However, we believe the relative patterns are more important than the exact numerical values for determining when and where the risk for RVFV outbreaks would be highest.
In parameterizing the model, we were confronted with many gaps in knowledge and data. This contrasts with our recent study [20] that illustrated a similar methodology for WNV, a system where the competent vectors are well-known and previous studies could be invoked to estimate the EIP and biting rates as functions of environmental temperature. For RVFV, the competent North American vectors are only partially known [6], and the data available for quantifying extrinsic incubation rates at different temperatures are scanty and based on colonized mosquitoes [42,43]. Similarly, there is little data available for North American Aedes species on the variation in biting rates with temperature; here we have used a model based on egg maturation rates in Culex mosquitoes [45] and applied it to both genera of RVFV vectors. This is complicated by fact that Aedes may take multiple partial blood meals during a gonotrophic cycle, so a rate based on ovarian development periods may understate the frequency of vector-host contact.
In our model, we assumed that the vectors feed on the two types of hosts represented by cattle and birds at rates proportional to their availability. Most studies on Cx. tarsalis have found a broad host range, including mammals, with a greater number of blood meals taken on avian hosts compared to mammalian hosts [33,36,37,[60][61][62][63][64]. These patterns appear to be a function of host availability and possibly other factors such as body size or defensive behaviors. Studies utilizing forage ratios to relate feeding to proportionate availability have found that passerines were often underrepresented in relation to their densities [36,63,65], while cattle were frequently the most common -and sometimes overrepresented -mammalian host [36,37,60,61,65,66]. Studies on Ae. melanimon have indicated a consistent preference for feeding on mammals over birds (w90% of total blood meals; [33,60,62,66,67]). In our study area, cows generally constituted a high proportion of the total hosts in each area (i.e., cows represented w50% and w90% of all hosts in 94.5% and 59.7% of the areas modeled, respectively), which meant that the opportunistic feeding assumption resulted in consistently high rates of biting on cattle for both species. The impacts of host preferences on the potential for RVFV transmission is an interesting topic for consideration in future modeling studies.
To understand the potential role of Aedes, Culex, and possibly other genera in the transmission of RVFV or other pathogens to cattle, we need a better understanding of heterogeneities in biting pressure from the various mosquito species. Other studies have shown promise in this regard, in the use of antibodies to salivary antigens as biomarkers for vector exposure (e.g., [68][69][70][71][72]) and specifically in characterization of the sialotranscriptome for Cx. tarsalis [73].
There are also substantial uncertainties in the viral transmission cycle should RVFV find its way into North America. Which wildlife species would serve as competent hosts is a key unknown. White-tailed deer are very abundant in the eastern US and have been hypothesized as potential wildlife carriers in North America [13], but the ability of these animals to become infected and develop viremias high enough to infect biting vectors remains undocumented. The potential role of lagomorphs (hares and rabbits) is also unknown, although these are important in the transmission cycles of several bunyaviruses in the family Bunyaviridae in California [30,34]. Our study therefore did not include wildlife hosts.
Despite the potential for RVFV outbreaks, it is unclear whether the virus could persist between years to become enzootic in the U.S. This uncertainty is due in part to the numerous gaps in the data available to inform model parameters described above, but also to the seasonality of temperatures and vector populations that ensure that conditions are never constant. Additional study is needed within a stochastic modeling framework to understand the range of potential invasion mechanisms (e.g., mosquitoes vs. human or animal hosts), as well as the potential mechanisms for RVFV persistence and whether they could permit the virus to avoid the possibility of stochastic fadeout during North America's winter. Host and vector movements will be important for shortterm spread of RVFV [18,19], but their necessity for interannual persistence will depend in part on whether the African paradigm of inter-epizootic maintenance in vertically infected mosquitoes holds true. If vertical transmision turns out to be an inadequate persistence mechanism in North American Aedes or Culex, movement could increase the likelihood of RVFV persistence by bringing infectious vectors and hosts into contact with new susceptible subpopulations [54]. Data on livestock movement are limited, but new Bayesian methods are giving hope that more can be done with incomplete data [74].

Conclusion
We have developed novel, generalizable methods to link mathematical models for RVFV with broad-scale spatio-temporal data for realistic landscapes. These methods could be useful for prioritizing when and where to focus control strategies (e.g., vector control or cattle vaccination) during an invasion. Many gaps in both data and knowledge remain, but this is an important step toward understanding the potential seasonal transmission cycles of RVFV and other vectorborne pathogens that may invade temperate North America.

Supporting Information
Text S1 Mathematical appendix describing the model, including definitions for state variables and parameters and results for sensitivity analysis. (PDF)