http://orcid.org/0000-0003-4496-8809
Figures
Abstract
West Nile virus (WNV)—a mosquito-borne arbovirus—entered the USA through New York City in 1999 and spread to the contiguous USA within three years while transitioning from epidemic outbreaks to endemic transmission. The virus is transmitted by vector competent mosquitoes and maintained in the avian populations. WNV spatial distribution is mainly determined by the movement of residential and migratory avian populations. We developed an individual-level heterogeneous network framework across the USA with the goal of understanding the long-range spatial distribution of WNV. To this end, we proposed three distance dispersal kernels model: 1) exponential—short-range dispersal, 2) power-law—long-range dispersal in all directions, and 3) power-law biased by flyway direction —long-range dispersal only along established migratory routes. To select the appropriate dispersal kernel we used the human case data and adopted a model selection framework based on approximate Bayesian computation with sequential Monte Carlo sampling (ABC-SMC). From estimated parameters, we find that the power-law biased by flyway direction kernel is the best kernel to fit WNV human case data, supporting the hypothesis of long-range WNV transmission is mainly along the migratory bird flyways. Through extensive simulation from 2014 to 2016, we proposed and tested hypothetical mitigation strategies and found that mosquito population reduction in the infected states and neighboring states is potentially cost-effective.
Author summary
The underlying pattern of West Nile virus (WNV) geographic spread across the United States is not completely clear, which is a necessary step for continental or state level mitigation strategies to reduce WNV transmission. We report a network model that explains the geographic spread of WNV in the United States. West Nile virus is a mosquito-borne pathogen that infects many avian species with different movement ranges. From our research, we found that migration patterns and routes play an essential role in the WNV spatial distribution. The virus spreads in all directions at short distances because of local birds and short-distance migratory birds. However, the virus also disperses long distances along the avian migratory routes. Our model is designed to be flexible and therefore can be used to explore spreading patterns of other infectious diseases in other geographic locations.
Citation: Moon SA, Cohnstaedt LW, McVey DS, Scoglio CM (2019) A spatio-temporal individual-based network framework for West Nile virus in the USA: Spreading pattern of West Nile virus. PLoS Comput Biol 15(3): e1006875. https://doi.org/10.1371/journal.pcbi.1006875
Editor: Carrie A. Manore, Los Alamos National Laboratory, UNITED STATES
Received: May 1, 2018; Accepted: February 17, 2019; Published: March 13, 2019
This is an open access article, free of all copyright, and may be freely reproduced, distributed, transmitted, modified, built upon, or otherwise used by anyone for any lawful purpose. The work is made available under the Creative Commons CC0 public domain dedication.
Data Availability: The data underlying the results presented in the study are available from the Centers for Disease Control and Prevention (CDC)(https://www.cdc.gov/westnile/statsmaps/index.html), the National Centers for Environmental Information (https://www.ncdc.noaa.gov/cag/), and the eBird database (https://ebird.org/barchart?r=US-KS&bmo=1&emo=12&byr=1900&eyr=2017&spp=amerob).
Funding: The work has been financially supported by the NSF/NIH/USDA/BBSRC Ecology and Evolution of Infectious Diseases (EEID) Program through USDA-NIFA Award 2015-67013-23818, the United States Department of Agriculture award ABADRU-427647, and the National Science Foundation award CIF-1423411. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
West Nile disease (WND) is a vector-borne zoonosis which may result from infection by West Nile virus (WNV), a member of the family Flaviviridae, genus Flavivirus. This virus is the most common cause of arboviral disease in the United States [1]. From 1999 to 2017, more than 48 thousands WNV disease cases were reported to the Centers for Disease Control and Prevention (CDC) and more than two thousands of these reported cases resulted in death [2]. WNV is maintained in an enzootic transmission cycle between competent mosquitoes and birds. Birds are the reservoir and amplifying host for this virus. The US Centers for Diseases Control and Prevention (CDC) has identified WNV infection in more than three hundred species of birds. Infected bird movement is likely a key factor that affects the geographic spread of WNV, especially given the different habitats and routes of various species. Although many bird species may be infected with WNV, the American robin is considered an important amplifier of WNV and maybe a driver geographic spread because WNV-infected American robins have low mortality and high viremia [3, 4]. Members of the Culex genus of mosquito are the principal vectors of this virus in the United States [5]. Humans, horses, and other mammals can be infected with WNV. However, these infections result in relatively low virus titers (viremia) therefore the infected animals and people are considered dead-end hosts (not capable of infecting feeding mosquitoes). Therefore, they do not have any epidemiological impact on WNV transmission or geographic spread [6].
To understand the transmission dynamics of WNV, several mathematical models have been developed [3, 7–10]. These models predict the threshold conditions for WNV spreading in different scenarios. However, most of these models do not consider the spatial dynamics of WNV. Space or geographic spread has a significant role in WNV disease dynamics and modeling of WNV spatial spreading is complex because of the interactions of multiple potential mosquito vectors, avian amplifiers, and mammalian hosts. Liu et al. [9] developed a patchy model to analyze the spatial spreading of WNV, where patches are geographical space. They assumed patches are identical, spatial dispersal of birds and mosquitoes are symmetric within patches, and movement of birds and mosquitoes are only one-dimensional. According to this investigation, long-range dispersal of infected bird populations determines the spatial spread of WNV, not the dispersal of infected mosquito populations. Other investigators proposed a reaction-diffusion model [10], where they have spatially extended the non-spatial model of Wonham et al. [8] to mathematically estimate the spread of WNV. Here, diffusion terms in the reaction-diffusion partial differential equations represent vector mosquito and host bird population movements. They identified traveling wave solutions in their model and calculated the rate of spatial spread of infection. Durand et al. [11] developed a discrete time deterministic meta-population model in order to analyze the circulation of WNV between Southern Europe and West Africa. Another spatial model proposed by Maidana and Yang [12] used a system of partial differential reaction-diffusion equations. They also calculated the speed of disease dissemination by investigating the traveling wave solution of their model. They concluded, mosquito movements do not play an important role in disease dissemination. In addition, they included vertical transmission in their model and determined that vertical transmission is not an important factor for the spatial spread of WNV.
Most WNV spread models are mathematical deterministic compartmental models. However WNV spread is highly stochastic because of the demography and movement of hosts and vectors varies between different locations. The major weaknesses of these models are the number and complexity of the compartments required to account for the many host and vector populations. In turn, the number of compartments increases the number of unknown parameters. Approximation of these parameters in any biological system is very challenging and prone to estimation errors which can create inaccuracies in the model outputs.
We developed an individual-based heterogeneous network framework to understand WNV geographic spread. To build the network framework, we used the American Robin population density across the contiguous United States. The demographic characteristics of avian host populations and vector populations are not homogenous geographically, so we used a heterogeneous network framework. The transmission intensity of WNV depends on the abundance of WNV-infected vector mosquitoes in a given location. Mosquito population numbers fluctuate with local weather and season throughout the year, therefore we used a temperature dependent transmission rate. Although dead-end hosts cannot spread WNV to mosquitoes, we have quantified WNV case data only for humans, which we used to estimate unknown parameters.
To understand the WNV spatial distribution, we proposed distance dispersal kernels, which describes the probability of dispersal with respect to distances. In this framework, we proposed three types of distance dispersal kernels: 1) exponential, 2) power-law, and 3) power-law biased by flyway. Then we compared the three distance kernels using approximate Bayesian computation based on sequential Monte Carlo sampling (ABC-SMC) method [13–18]. After conducting an extensive simulation for 2014-2016, we observed that an adapted fat-tailed or power-law kernel, which has long-distance links in specified directions can best describe the WNV human case data. We tested this network framework for the best kernel with the human case data and found that simulated results for more than 41 states of 49 states are consistent with the reported WNV cases. Our results support previous work on WNV spreading [3], which also modeled WNV spreading with migratory birds. We validate our work computationally from human incidence data. We proposed several theoretical mitigation strategies to control WNV and calculated their estimated costs. From the analysis of mitigation strategies, we suggest that potentially effective mitigation policies would include the application of mitigation control in areas with active transmission and in immediate neighboring states.
Materials and methods
In this section, we present data sources, an epidemic model for WNV, then develop a network framework for WNV geographic spread across the United States. At the end of this section, we present a statistical tool, approximate Bayesian computation with sequential Monte Carlo sampling (ABC-SMC) for parameter estimation and model selection.
Data
The study area of this research was the contiguous United States where WNV is considered endemic. We modeled WNV case distributions for 2014-2016. We used three data sets each year to develop our model. The first dataset contained the average monthly temperatures. Mosquito vector abundance correlated with temperature. Temperature data was from the National Centers for Environmental Information [19]. The second dataset contains American Robin population data from eBird [20]. This is a database for bird abundance and distribution, which is formed by the Cornell Lab of Ornithology and National Audubon Society. We used total observation of American Robin in each state of the USA for each month. The robin data set was used to train the network model. The American Robin is abundant throughout the United States and is a preferred food source for many WNV-competent mosquito species [21]. Based on host feeding patterns of the Culex genus of mosquitoes, robins are the most common WNV amplifying host [22–24]. Other important susceptible birds, such as American crow were not used because although they are an indicator species (high crow mortality), they are unlikely to spread virus geographically as they are mostly a residential species. In addition, as an indication of epidemic start point, we used WNV human incidence data. Many species of birds have long-distance migration during the spring and fall. Therefore the network does not focus on one long-distance migrating bird species but aggregates all species along the known flyways. To estimate model parameters we used human case data for WNV from CDC [2], which is the third dataset.
WNV epidemic model
To explore WNV long-distance spatial distribution in the USA, we used an individual-based heterogeneous network framework. In this framework, birds are on the individual level, a node represents an individual bird and connection between nodes is the possibility of virus dispersal from one infected bird to another susceptible bird by mosquito vectors. Links or connections are formed by movement of birds or movement of vectors. If there is no link between nodes then infected birds and insects are not moving virus between nodes. All virus transmission occurs by local competent vector mosquitoes. There is some evidence of bird-to-bird transmission, but it likely does not contribute to or maintain outbreaks. We split the bird population into four compartments; susceptible, exposed, infected, and recovered. Although, in the literature most mathematical models do not consider the exposed avian class when modeling WNV [8, 12, 25, 26]. Birds transmit virus to mosquitoes when a susceptible mosquito vector takes an infected blood meal, then the mosquito becomes infectious after the extrinsic incubation period (EIP), or the time needed for the virus to spreads from the mosquito mid gut to the salivary glands; usually this process takes 7 to 14 days [3, 27]. In addition, an infected bird can infect many mosquitoes simultaneously and also an infected mosquito can bite many susceptible or infected birds. Therefore, there is some delay in the system, to represent this delay we added the exposed class. We estimated exposed period from data by using the approximate Bayesian computation with sequential Monte Carlo sampling (ABC-SMC) method. After the exposed period, birds entered the infected compartment and an infected bird transitions to recovered after 4-5 days. To simulate this model, we used generalized epidemic mean-field (GEMF) framework developed by the Network Science and Engineering (NetSE) group at Kansas State University [28]. In GEMF, each node stays in a different state and the joint state of all nodes follows a Markov process [28–30]. The node level description of this Markov process is:
(1)
(2)
(3)
Here, X(t) is the joint state of all individual nodes at time t. xi(t) is a node state, xi(t) = C means node i is in C compartment at time t, C = 0, 1, 2, 3 corresponds to susceptible, exposed, infected, and recovered compartment. Yi is the number of infected neighbors of node i, β(T) is the transmission rate from one infected bird to one susceptible bird, which is a function of temperature, λ is the rate for exposed to infectious state, and finally, a node recovers from infectious state at a rate δ.
Zoonotic spillover transmission.
To model disease transmission from the bird population to human population, we added a zoonotic spillover transmission compartment. We modeled occurrence of human cases as a Poisson process [26, 31]. This part of the framework can be expressed as the following equation:
(4)
In this equation,
is number of infected human cases at n sub-network in s time steps, where s = 1, 2, 3…‥ are the discrete time steps,
is infected bird population in sub-network n, and η is a scaler quantity, accounts for the contact rate and probability of pathogen transmission from bird to human. We calculated WNV spilling over to humans by using a Poisson random number generator.
Temporal transmission rate and environmental conditions.
The transmission rate for WNV is sensitive to weather data as mosquito abundance depends on the environmental conditions. Temperature, precipitation, landscape features, daylight conditions etc. are environmental conditions, which has an impact on the transmission dynamics of WNV [32]. In this research, we considered average monthly temperature data, optimal mosquito season [33], and suitable temperature range for co-occurrence of WNV and competent mosquito species. Temperature plays a very important role in the transmission dynamics of WNV because mosquito longevity and EIP are sensitive to temperature. Mosquito longevity and EIP decrease with the increase of temperature. However, there is no straightforward relationship of vectorial capacity for WNV with temperature. If incubation period decreases more than longevity, then mosquitos will be infective longer. However if longevity decreases more than incubation period, then mosquitos will not be able to transmit the virus. We used information about rainfall in this research implicitly through optical mosquito season. Optimal mosquito season of any location is estimated from monthly average temperature and rainfall data for that location [33]. In this model, we used a simple linear relation of transmission rate with temperature in a temperature window from 12°C to 32°C in the optimal mosquito season. Outside this window, transmission rate is very low. Suitable temperature for co-occurrence of WNV and Culex pipiens is around 12° to 27°C and for Culex quinquefasciatus is 20°C to 32°C [33]. Survival rate to adult stage for Culex quinquefasciatus is significantly high when temperature is in 20°C to 30°C [34]. For Culex tarsalis favorable temperature for WNV development start after 14°C [35], however larval survival reduced after 30°C temperature [36]. To compute the transmission rate of any link from node a to node b, we used temperature of the location of node b. Transmission rate for a location l is, βl(T) = β∘(Tlm − T∘); here, β∘ is the proportional constant, what we estimated by using ABC-SMC method, Tlm is the average temperature for month m in location l and T∘ is the threshold temperature. Threshold temperature for this model is 12°C. As the temperature is space dependent, our transmission rate also differs across the network. This individual-level heterogeneous network model gives us this flexibility to use different transmission rate at a time for different parts of the network.
Network framework
For the spatial dynamic characteristics of WNV transmission, we built a network framework, which has 49 sub-networks one for each adjoining states of the contiguous United States plus the District of Columbia. The number of nodes in each sub-network is proportional to the size of the avian population in that state [20]. We considered the mosquito season June-October for the simulation period. Although the mosquito season is not the same for all states, mosquitoes are active from June to September in all of the states at these times [33].
The network for the avian population is (V, E). Here, V is the set of nodes, which is the union of nodes of all sub-network, V = SN1 ∪ SN2 ∪ SN3 ∪ ………‥ ∪ SN49, here SNi is a set of nodes in the sub-network i and E is the set of links among individual nodes. To build sub-networks, we used the total number of observations of American Robin for states per month in the simulation time period. , here,
is the total number of observations of American Robins in state i in month mj, N0 is the error term and N0 ∼ N(5, 2) for this model. m1 is the first month after May and m2 is the last month before October when the average monthly temperature is greater than T0. Sc is the scaling constant.
In each sub-network, we assumed that nodes are connected through Erdos-Renyi (n,p) random network topology [37]. In this network topology, we created links randomly among nodes with a probability p. Here, n is the number of nodes in a sub-network and p is the probability to form an edge. We set the probability p = R * log(n)/n, here R is a constant (R ≥ 2), as this value is more than the threshold value for the connectedness of an Erdos-Renyi graph [38], so nodes of a sub-network are locally connected. We will refer these networks as a local network in the subsequent sections of this paper. To build connections among sub-networks, we considered long-distance dispersal kernels [7, 39], which describe the probability of dispersal with respect to distances. Dispersal kernels provide a simple model of dispersal to model dispersal events. For long-distance events, we used three types of kernel models; 1) Exponential, 2) power-law, and 3) power-law-flyway, which is a power-law kernel biased by flyway. The dispersal phenomenon in this work is not conserved because of long-distance movement of migratory birds and seasonality within bird populations. Some long-distance migratory birds can disperse outside the contiguous United States or outside the network nodes, which are discrete points. The connection probability between two nodes does not represent the probability that a single dispersal event happens rather it represents the probability of contact and subsequent pathogen transmission between them. A simple caricature of the network is shown in Fig 1. There are three sub-networks, A, B, and C. The links, which formed local networks are shown by solid lines. These links are introduced by Erdos-Renyi (n,p) network topology. Dashed lines are inter-links among sub-networks. These links established by using long-distance dispersal kernels.
Here, A, B, C are three sub-networks. Solid lines represent intra-links in a sub-network and dashed lines represent inter-sub-network links.
Exponential distance kernel.
In this distance kernel, connection probability among sub-networks will decrease exponentially with distance. Probability to form a link is:
(5)
Here, dij is the distance between sub-network i and j, Ke is the shape parameter of exponential distribution kernel. For distance between two states, we took the distance between their centroids. The network with the exponential dispersal kernel was created as follows:
- Step 1. Calculate the distance among sub-networks. dij is the distance between sub-network i and j.
- Step 2. Calculate P(dij), this is the probability to form a link between sub-network i and j.
- Step 3. Generate a random number rand for each pair of nodes (a,b), where a ∈ i and b ∈ j.
- Step 4. If rand < P(dij) then an undirected link will form between node a and b.
Inter-links among sub-networks, generated by exponential distance kernel are shown in Fig 2a.
a) for exponential distance kernel, b) for power-law distance kernel, and c) for power-law distance kernel biased by flyway. Gray links represent undirected links and orange links represent directed links (for spring migration –northbound; for late summer/fall migration –southbound). Intra-links are not visible here. These are one realization of the stochastic networks, which are rescaled by 0.1 for better visualization.
Power-law distance kernel.
Power-Law, heavy-tailed, or fat-tailed distribution allows occasional long-range transmissions of infection with frequent short-range transmissions. In this fat-tailed distance kernel, there is a greater chance of creating links over the same long-distances compared to the exponential kernel. Power-law transmission kernel was used previously to model spatial dynamics of several infectious diseases, for example, in plant epidemiology [40], in 2001 foot-and-mouth disease epidemic [39], and also, in human diseases [41]. In power-law connections [42], the probability of connectivity among sub-networks will decrease with distance according to the following equation:
(6)
Here dmin is minimum distance among sub-networks and Kpl is the power-law parameter. The process to build this network is similar to a network for exponential kernel with the only difference being the calculation of P(dij). Inter-links among sub-networks for power-law distance kernel are shown in Fig 2b.
Power-law distance kernel biased by flyway.
To form this distance kernel, we included the migratory behavior of birds. Migratory birds can spread pathogens during the migration periods [43, 44]. According to the United States Fish and Wildlife Services and Flyway Councils, there are four flyways in the United States; the Atlantic flyway (AF), the Mississippi flyway (MF), the Central flyway (CF), and the Pacific flyway (PF) [45]. Although flyways overlap and the migratory patterns are very complex, these migratory routes play a vital role in the long-distance spreading of WNV [46]. To build this distance kernel, we considered two types of links among sub-networks; 1) links which are formed for residential or short-distance migratory bird movements and 2) links which are formed for long-distance migratory bird movements. For the first type of links, we used an estimated movement range of 500 km [47], these connections are unrelated to flyways. For the second type of connections, we considered two migration periods; spring migration (April—June) and late summer/fall migration (July—September) [30]; during the spring migration, we established long links from south to north and in late summer/fall migration, the reverse. To establish any long link, we picked two sub-network and establish a link if they were in the same flyway with probability P(dij) (Eq 6), these links were directional and direction was imposed with respect to migratory period. Inter-links among sub-networks for this kernel were shown in Fig 2c. The algorithm to create this network was:
- Step 1. Calculate the distance among sub-networks. dij is the distance between sub-network i and j.
- Step 2. Calculate P(dij) using Eq 6, this is the probability to form a link between states i and j.
- Step 3. Generate a random number rand for each pair of nodes (a, b), where a ∈ i and b ∈ j.
- Step 4. If rand < P(dij) and dij < 500km then an undirected link will form between node a and b.
- Step 5. If rand < P(dij) and dij > 500km and states i and j are in the same flyway then an directed link will form between node a and b according to the migration period.
Temporal network behavior.
Bird populations are not constant in any region, they change with time because of bird movement. To consider this fact, this study adds a node property, namely, Activity. This property can hold two values: 1 = Active and 0 = Inactive. In the entire network, only Active node can contribute to the spreading of the WNV. By controlling this property, we varied the size of the active node population in any sub-network with respect to the variation of the avian population in that region. The length of the simulation each year was five months (June—October). Then, each month nodes are activated randomly according to the total number of birds observed in that region in that month.
ABC-SMC for parameter estimation and model comparison
In this framework, we adopted approximate Bayesian computation based on a sequential Monte Carlo sampling (ABC-SMC) method for parameter estimation and model selection [13–18].
Parameter estimation.
ABC-SMC is a computational method of Bayesian statistics that combines a particle filtering method with summary statistics. This method is ideal for a stochastic complex model where likelihood function is intractable or computationally expensive to evaluate. ABC estimates the posterior distribution of parameters from data. Let, θ is a parameter vector to be estimated. The goal of the ABC is to approximate the posterior distribution, Π(θ|d) ∝ f(d|θ)Π(θ), where prior distribution of parameters Π(θ) are given and f(d|θ) is the likelihood of θ given the data d. This method samples parameter values from their prior distribution through subsequent SMC rounds. Intermediate distribution of the parameter is Π(θ|dist(x, d) ≤ ϵi);i = 1, 2, ….P. The target posterior distribution is Π(θ|dist(x, d) ≤ ϵP). Here, x is the simulated data set, dist is the distance function, ϵ is the tolerance and P is the number of SMC rounds or the number of populations, where ϵP < …‥ < ϵ2 < ϵ1 [48]. This is an adapted sequential importance sampling. In each SMC round, it uses perturbation kernel to sample a parameter set. After each simulation of the model, the model output and data are compared using some goodness-of-fit metrics. A parameter set is accepted if the distance between the model output and data is less than the tolerance level. The accepted parameter set is a particle and accepted particles form a population for that SMC round. We used two goodness-of-fit metric or distance function in this research. The first goodness-of-fit metric is squared root of the sum of squared error between observed incidence data and simulated incidence data for any proposed parameter set. The first goodness-of-fit metric for this model is:
(7)
Here, x(i,j) is simulated incidence model data for i week and for j location. The second goodness-of-fit metric is the absolute difference between the number of infected states from observed data and simulated data, infected state defined as a state where at least one infected individual has reported. The ABC-SMC algorithm, we adopted for this model from Toni et al. [13], which has given in S2 Text. We used this algorithm separately for estimating parameters for this three distance dispersal kernel network models. As our models are an event based stochastic simulation, we simulated them 30 times with GEMF for each particle to get 30 realizations of the system. Then we take the average of these realizations. As the average over the multiple runs of a stochastic system holds more information than a single stochastic run.
Model comparison.
In many areas, researchers deal with model selection. Bayesian theory is a comprehensive method to make inference about models from data. Approximate Bayesian computation was used in many research areas for model selection [49]. To compare among three distance kernels, this investigation used ABC-SMC model selection framework [13, 50, 51]. For given data d, the marginal posterior probability of model m is:
(8)
Here, Pr(d|m) is the marginal likelihood and Pr(m) is the prior probability of the model. We used a uniform distribution for prior distribution of unknown parameters. For each model, we have four unknown parameters; network parameter K (Ke is the network parameter for the exponential kernel and Kpl is the network parameter for the both power-law kernels), constant for transmission rate β0, transition rate from exposed to infectious state λ, and zoonotic transmission spillover rate η. In each population, we took 1000 particles. We used Bayes factor to compare a model with another model. For model mi and mj, Bayes factor [52] is,
(9)
Here, Pr(mi) is the prior and Pr(mi|d) is the marginal posterior distribution of model mi. The Bayes factor is a summary evidence in favor of one model over another supported by the data. If Bij is in range 1-3, we can conclude that summary of the evidence against mj in favor of mi is very weak. If Bij is in range 3-20, we can conclude that summary of the evidence against mj in favor of mi is positive [52]. The ABC-SMC model selection algorithm is very similar to the algorithm for parameter estimation. Here, m is the model indicator, m ∈ 1, 2,…‥, M, M is the number of model. In this research, we had three network models (M = 3) to compare.
m = 1: exponential kernel network model,
m = 2: power-law kernel network model, and
m = 3: power-law kernel influenced by flyway network model.
In each population, the model selection algorithm starts by sampling the model parameter m from the prior distribution Π(m). Then the algorithm proposes a new set of parameters (particle) from the sets of parameters of the model m from the previous population. The Bayes factor was calculated from the final population of m. The algorithm for model selection has given in S2 Text. Although ABC-SMC is an accurate statistical tool for parameter estimation and model selection, however, the results of this method are sensitive to summary statistics [53]. For our case, no summary statistics were required because we used the entire set of data and we compared the simulated and observed dataset directly by using goodness-of-fit or distance metric. A full dataset is sufficient to get the consistent result from approximate Bayesian Computation [54].
Mitigation strategies
The role of mosquito populations in WNV transmission is expressed by disease transmission rate β. This framework used different transmission rates in different parts of the network corresponding to the local mosquito abundance. Using this heterogeneous feature in the framework, we evaluated theoretical mosquito population management measures to reduce the outbreak size or transmission rates in the state level. Some states such as Kansas, do not have statewide mosquito surveillance or management, but in these theoretical scenarios, it is assumed they can develop or benefit from effective statewide mosquito management programs. The framework will simply estimate how much the mosquito abundance is reduced or maintained based on the theoretical outcomes of coordinated control. Furthermore, we realize mosquito control is generally conducted on a county or municipal level, but the human case data is only available on a state level. Therefore the recommendations are for the lowest resolution of the data, which is state level but applies to counties and municipalities as well. If vector management is increased in a sub-network, then transmission rates will be changed by, , here βr is the reduced transmission rate and RF is the reduction factor. Then management costs will be Cost = RF * NSc, here NSc is the number of states where control measures were applied. We considered supplemental management measures with the existing management measures. We used two types of mitigation strategies across the United States, 1) dynamic infected place tracing strategy and 2) static ranked based strategy.
In the infected place tracing, we traced the infected states, then plan the mitigation strategies according to them. For this type of mitigation strategies, we considered three cases; 1) case-1: only infected: applied control only in the infected states; 2) case-2: infected & first neighbors: applied control in the infected states with its first neighboring states (whose distance is less than 500km), and 3) case-3: infected & first neighbors & second neighbors: applied control in the infected states with its first neighboring states, and also with its second neighboring states (whose distance is in 500 − 1000km). For infected tracing control measure, we kept track of infected places monthly. If SNi sub-network is infected for month t, then control measures were applied for the month t + 1 based on these three cases.
In the static ranked based mitigation strategy, we ranked the states by different variables (for example, temperature, size of the avian population etc.). For this strategy, we considered three cases; 1) temp.: states ranked by temperature, 2) pop.: states ranked by avian population size, and 3) temp. & pop.: states ranked by temperature and avian population size both, then we applied management measures in the top 30% of the states.
Results
We developed a novel flexible individual based heterogeneous network framework to test three WNV dispersal kernels across the contiguous United States based on human case data distributions. We used this framework for the year 2014, 2015, and 1016. The results for network formulation, parameter estimation, and dispersal kernels selection using Bayesian inference are given below for the year 2015 and the results for other two years are given in the S1 Text.
Network framework
In this spatial-temporal individual-based heterogeneous network framework, we used three distance kernel models. The fundamental basic WNV epidemic model is the same for all the three network kernels. In the entire network, there are 49 sub-networks representing the 48 adjoining contiguous states plus the District of Columbia. All sub-network nodes are locally connected. The topology of the local network is Erdos-Renyi. The total nodes for the year 2015 was |V| = 7657 and the scaling constant is Sc = 0.02. Here, E = El ∪ Edd; |El| is the number of total intra-links for all local networks, which is around 167000-170000 and |Edd| is the number of total inter-links among sub-networks. The description of sub-networks is provided in Table B in the S3 Text. We started the epidemic from states with the highest human incidence prior to June. We started the epidemic for the year 2015 by adding two infected nodes, one in sub-network SN4 (California) and another in sub-network SN42 (Texas). Connections among sub-networks are developed by distance dispersal kernels. Parameters for these kernels are estimated from the ABC-SMC method.
ABC-SMC for parameter estimation and model comparison
Parameter estimation.
ABC-SMC parameter estimation was applied to three dispersal kernel network models separately. For each set of prior distributions, convergence to the posterior distribution was achieved after 13-15 SMC rounds. Convergence of the posterior distributions was monitored by visual inspection of the outputs from consecutive SMC rounds. The prior distribution for exponential network parameter was, Ke ∼ U(0.1, 0.3), for power-law Kpl ∼ U(2, 4), for power-law biased by flyway was Kpl ∼ U(2, 4). Prior distribution for constant of transmission rate β0, transition rate from exposed to infectious λ, and human spillover rate η is same for three kernel models; β0 ∼ U(0, 15), λ ∼ U(0.025, 10) and η ∼ U(0, 50). Perturbation kernels were also uniform, PK = αU(−1, 1), with α = 0.5(maxθp−1−minθp−1), here θp−1 is the set of a parameter values in the previous population. We used weekly human case data for 49 locations, as observed data. The estimated parameters for this three dispersal kernel network models for 2015 are presented in Table 1.
Model comparison.
ABC-SMC for model selection allows us to estimate posterior model distributions. We used this algorithm to compare the three distance kernels. Prior distributions and perturbation kernels are the same for both the model selection and the parameter estimation algorithm. Here we used one more prior distribution for discrete model parameter; m ∼ U(1, 3). The tolerance vector for ABC-SMC model selection algorithm is, ϵ = {2200, 2000, 1800, 1600, 1400, 1200, 1100, 1000}. The target and intermediate distributions of model parameters are shown in Fig 3.
Model-1 represents exponential kernel, model-2 represents power-law kernel, and model-3 represents power-law influenced by flyway kernel. Here, Population-8 is the approximation of the final marginal posterior distribution of model parameter m and population 1-7 are intermediate distributions. Population-0 is the discrete uniform prior distribution, which is not shown here.
We calculated the Bayes factor from the marginal posterior distribution of m, which we took from the final or last population. In the final population for 2015, exponential distance kernel model (m = 1) was selected for 64 times, power-law distance kernel (m = 2) was selected for 95 times and power-law influenced by flyway distance kernel model (m = 3) was selected for 841 times. Bayes factor B3,1 = 841/64 = 13.1406, B3,2 = 841/95 = 8.8526. In the marginal posterior distribution of three models, there is positive evidence in favor of power-law influenced by flyway distance kernel when compared with other two models [13]. The distribution of parameters for power-law influenced by flyway for 2015 are presented in Fig 4. Calculation of the Bayes factor for 2014 and 2016 are provided in the S1 Text.
a) Network Parameter K; b) constant for transmission rate β0; c) transition rate from exposed to infectious node λ, and d) human spillover η.
Performance of the power-law-flyway network model
To test the performance of this framework, we used estimated parameters from Table 1 for power-law kernel influenced by flyway. We set the parameters value; Kpl = 2.3147, β0 = 0.0059day-1, λ = 0.0721day-1, and δ = 0.2031day-1. The simulation period for the avian population model is from week-23 to week-44. The output of avian population was used as the input of zoonotic spillover compartment. Then we compared the output of zoonotic spillover compartment with human case data for week 24 to week 45. We considered a one-week lag between WNV incidence in birds and WNV incidence in humans. In humans, WNV-infected individuals (approximately 20%) develop a mild febrile illness after 3–6 days [56]. Peak of reporting of dead birds is one week prior than the reporting peak of human incidence [57]. In Fig 5, the mean simulated human case from the 49 sub-networks is compared with the weekly human case data for 2015 for the contiguous USA. The absolute errors between them are shown here. From this whisker plot, we can see that the median of the absolute error for the states is close to zero. In Fig 5, the largest outlier is California (marked by black circles). These outliers result from a mismatch between the simulated peak human incidence time and the observed human incidence peak time possibly because the very long state (north to south) has weather which is very different in southern California (warmer and drier) than northern California (cooler and wetter) causing a difference between peak mosquito seasons in the southern and northern parts.
Mean of 1000 realizations has used as the simulated data. On the blue boxes, the red horizontal lines show the median and the bottom and top edges of the boxes indicate 25th and 75th percentile respectively. The whiskers show the ranges of data points not considered outliers and outliers are showing by red + symbol. Californian outliers are marked by black circles.
We compared the total yearly incidence of human WNV from this model with the state level reported case data. The results are shown in Fig 6. For 2015, we found that the case data for 42 of 49 locations were within the simulation results. The states where human cases were different from the simulation results were over-reported states (Nevada) and under-reported states (Louisiana, Mississippi, Nebraska, North Dakota, South Dakota, and Washington). The possible reason for this mismatch are reporting error or overwintering of virus in birds or mosquitoes or another bird species (not robins) is the key reservoir species for that state
States name are given in the short form. Simulated results are represented with a box plot in which the red horizontal lines show the median and the bottom and top edges of the boxes indicate 25th and 75th percentile respectively, The whiskers show the ranges of data points not considered outliers and outliers are showing by red + symbol. Broken scale is used for sake of visualization.
To build a disease prevalence map, we grouped the states in four categories; 1) higher prevalence —incidence is more than 100, 2) intermediate prevalence—incidence is in between 50-99, 3) moderate prevalence—incidence is in between 25-49 and 4) low prevalence—incidence is less than 25. To group the states, we used the median of the simulation results. The disease prevalence map from the model are presented in Fig 7a and from observed data are presented in Fig 7b. Among 49 locations, 40 locations are in the same prevalence group in both maps.
The darker regions imposed greater prevalence. States are divided into four groups by incidence number; group-1: more than 99, group-2: 50-99. group-3: 25-49, and group-4: less than 25 incidences. a) States are divided by the median of the output of 1000 simulations, b) states are divided by observed data.
Mitigation strategies
We applied mitigation strategies on the power-law-flyway kernel network model to find the optimal mitigation plan. Fig 8a shows the number of infected states or epidemic size for dynamic infected places tracing. Epidemic size decreased faster with increased reduction factor for case-2(infected & first neighbors) and case-3(infected & first neighbors & second neighbors) than case-1(only infected). The number of states where control measures were applied is displayed in Fig 9, which is proportional to cost. Therefore, the cost was minimal for case-2 than other two cases for RF > 2. From the cost analysis, we concluded that, although the cost for case-1 is less at the beginning of the yearly outbreak, we need to apply management only in the infected places, however by the end of the year the total cost for case-2 will smaller because of the smaller epidemic size.
a) Dynamic infected places tracing; case-1: control measures are applied only in the infected states, case-2: control measures are applied in the infected states plus in their first neighboring states, case-3: control measures are applied in the infected places plus in their first and second neighboring states, and b) static ranked based strategy –states are ranked by; 1) temperature (Temp.), 2) avian population size (Pop.), 3) both(Temp & Pop.), then control measures are applied in the top 30% states. Log scale has used in x-axis for better visualization.
Log scale has used in x-axis for better visualization.
The results of the static ranked based mitigation strategy measure are presented in Fig 8b. We observed that, before RF = 4.5, number of infected states for temp. & pop. dropped earlier than others. Number of infected states or epidemic size was smaller for temp. than pop. after RF > 3, infected population of a sub-network are more positively correlated with temperature. The NSc is always the same for these three cases. For all mitigation strategies, minimum epidemic size could be 2, as we started the epidemic from two states.
Discussion
We proposed an individual-based heterogeneous network framework and tested three dispersal kernels to understand the spatial spread patterns of WNV human case data across the contiguous United States.
This framework requires fewer parameters and has more flexibility to represent the spatial-temporal dynamics of WNV. Adding parameters can make the framework more realistic, for example, more competent bird species, landscape features for habitat preferences of host and vector species, daylight conditions [32], pathogen invasion from outside of USA, variable susceptibility among different hosts and vectors, WNV strain variability, mosquito and virus overwintering, vertical transmission, human movement characteristics etc. However, inclusion of too many factors increases model complexity which makes model optimization difficult given the availability of limited observational data. On the other hand, a simple model may insufficient to represent WNV spatial dynamics. Computational models need to be developed and parameters calculated with sufficient detail to be biologically accurate if they are used to evaluate epidemic management measures. However, for most biological systems, reliable parameter information is unknown. Unknown parameters or inaccurate assumptions add uncertainty to the model. Our framework has only four parameters to estimate (network Parameter K, transmission rate β, transition rate from exposed to infectious state, λ, and human spillover, η). This framework has compartments only for the avian population (susceptible, exposed, infected, and recovered), and is not species specific. We reduced the compartments for vector population by implementing them implicitly through transmission rate between infected nodes and susceptible nodes. The presented framework and dispersal kernel network model has an intermediate complexity that approximate Bayesian computation based on sequential Monte Carlo sampling (ABC-SMC) method successfully calibrated and estimated the parameters with the available data. If more data becomes available, it is possible to add them in this model for improved performance of the model.
Furthermore, this framework is flexible and therefore can represent various hosts and vectors including with population seasonality, which plays an important role in WNV dynamics. For host population seasonality, we added a node property Activity, this property allows us to control active host populations in the network in a specific time period. We added vector seasonality in this framework with a temperature dependent transmission rate. This framework proposed one exponential and two fat-tailed distance kernel models for long-distance transmission of WNV with each model having increasing complexity and similarities to natural avian movement. WNV spatial distribution is very complex because WNV can infect more than 300 bird species, some of which are residential birds and short-distance migrators which disperse less than 500 km distances (short connections) whereas some species are long-distance migratory birds creating long connections. The long-distance migratory birds are the long-distance dispersal (LDD) agents for WNV. Previous studies tried to analyze spreading of WNV using a traveling wave with constant velocity, however, WNV spread more rapidly across the North America than would be expected from the assumption of constant velocity traveling wave [58]. Likely this is because traveling wave models unlike distance dispersal kernel models for WNV spreading do not capture the long-distance migrating birds which can have various migratory ranges and distances. Distance dispersal kernels have more flexibility to represent the different bird migration distances and can account for accelerating invasions. However, exponential kernels produce short-connections and therefore like traveling waves are limited to constant expansion, unlike fat-tailed power-law kernels which can generate accelerating invasions by creating the long-distance connections from migratory birds [59]. However, a general fat-tailed power-law kernel makes long-distance links in every direction which does not follow the incidence of WNV. Instead, a power-law-flyway kernel can be used to produce the long connections in the direction of flyways and short links in other directions. Bayesian inference was used to test which of the three kernel models best described WNV distribution on the network for three most recent years (2014- 2016). The power-law-flyway kernel best described the distribution of WNV cases because the long-range WNV transmission was concentrated mainly along the migratory bird flyways. The general power-law kernel overestimated the incidence data in some states because it was creating long-distance links in all directions.
The performance for the power-law-flyway dispersal kernel model was evaluated for the three most recent years (2014-2016) when WNV was endemic in the USA. The observed case data for the 49 locations were within the range of the simulated results for 41 states for 2014 (Fig B in S1 Text), 42 states for 2015 (Fig 6), and 45 states for 2016 (Fig D in S1 Text). For all three years, the simulated results were similar to the observed data, except in Colorado, Louisiana, Mississippi, Nevada, Nebraska, North Dakota, and Washington. Nevada was over-reported for 2015 and all others were under-reported. The power law flyway dispersal kernel network model reported more WNV human incidence in Nevada than reported cases, one possible reason for over-reporting cases in Nevada has rural areas, which tend to under report human cases, whereas mosquito control districts and health departments, focused in urban areas, must test birds and mosquitoes, which explains why CDC reported WNV infected mosquitoes in 25% of counties in Nevada. The under-reported states had more human cases than predicted by the model. Under-reporting by the power-law-flyway kernel network model is likely because overwintering of the virus in some states (for example, Louisiana, Mississippi etc.), which was not considered. The overwintering infected Culex mosquitoes can stay in hibernacula such as sewers, houses, caves, and other warm areas in urban, suburban, and rural areas and initiate the outbreak in the spring. Furthermore, there may be under-reporting of cases by the model if robins are not the main reservoir species in a state, which would be predicted between gulf coast states (Louisiana and Mississippi) and northern states such as North and South Dakota and Washington.
Stochastic simulations are useful tools to select the optimal future mitigation strategy after outbreaks of invasive species and pathogens. The foot-and-mouth disease (FMD) epidemics in 2001 in the United Kingdom developed by Keeling et al. [60], and mitigation strategies for pandemic influenza in the United States [61] are two well developed models with similarities to the current model. These models explore possible control measures such as culling, vaccination etc. for FMD [60], and vaccination, quarantine etc for influenza [61]. Most of these strategies can be examined with the network framework however, avian culling or vaccination for WNV control is not feasible. Vector control (or mosquito control) is a viable mitigation strategy for WNV, which is not considered by the other two models (FMD and influenza). To be applicable to any pathogens and inclusive of new mitigation methods, the mitigation strategies are non-specific and the predicted effectiveness of the mitigation methods can be adjusted to other methods. In the planning of the mitigation strategies, there is a trade-off between control measures effectiveness and their cost both monetary and loss of life. A stochastic simulation tool can decide the optimal mitigation strategy by dealing with this trade-off.
Mitigation strategies for WNV were tested using the power-law-flyway dispersal kernel network model. The mosquito management measures are not specific to larvae or adults, rather simply generally accepted best practices to reduce mosquito abundance for the purpose of reducing pathogen transmission. The mitigation strategy analysis proposes supplemental measures in addition to the existing mosquito management in each state because the states had yearly reported WNV cases despite the existing management methods. To reduce WNV spread, a theoretical policy would be management in neighboring regions and not exclusively in the infected places. Although this approach can cost more at the beginning of the epidemic season however at the end, it can reduce total cost by decreasing the size of the epidemic. If management measures are applied only in the infected states, it is not possible to control the epidemic because of long-distance migratory birds. This is statewide management in a unified effort. We acknowledge that states do not conduct mosquito management in this way, but to test the spillover it was necessary to do the simulation in this way because only state-level data was available. The infected place tracing mitigation technique has been used to control other diseases (for example, FMD, influenza etc.), although their host population and control measure means are different, however, the main concept behind the mitigation techniques are similar. The findings from this research to control WNV epidemic can be useful to select optimal mitigation strategies for other pathogens.
This research showed that the inclusion of directional long-distance dispersal of migratory birds improves model representations of the spatial patterns of WNV spread in the United States. The simulation of our framework in the context of long-distance directional dispersal suggested that cooperation and communication can facilitate early treatment and reduced outbreak sizes because of reduced WNV dispersal by American robins.
Supporting information
S1 Text. Simulation results for 2014 and 2016.
https://doi.org/10.1371/journal.pcbi.1006875.s001
(PDF)
S2 Text. Approximate Bayesian computation based on sequential Monte Carlo sampling (ABC-SMC) method.
https://doi.org/10.1371/journal.pcbi.1006875.s002
(PDF)
S3 Text. The sub-networks description for 2014-2016.
https://doi.org/10.1371/journal.pcbi.1006875.s003
(PDF)
Acknowledgments
We would like to thank Dr. Christopher Mundt for his excellent guidance in our EEID (Ecology and Evolution of Infectious Diseases) project. We also would like to thank the developers of the GEMF software [28].
References
- 1. Burakoff A. West Nile Virus and Other Nationally Notifiable Arboviral Diseases—United States, 2016. MMWR Morbidity and Mortality Weekly Report. 2018;67.
- 2.
Centers for Disease Control and Prevention;. Available from: https://www.cdc.gov/westnile/index.html.
- 3. Bergsman LD, Hyman JM, Manore CA. A mathematical model for the spread of West Nile virus in migratory and resident birds. Math Biosci Eng. 2016;13(2):401–24 pmid:27105987
- 4. Komar N, Langevin S, Hinten S, Nemeth N, Edwards E, Hettler D, et al. Experimental infection of North American birds with the New York 1999 strain of West Nile virus. Emerging infectious diseases. 2003;9(3):311. pmid:12643825
- 5. Turell MJ, Dohm DJ, Sardelis MR, O’guinn ML, Andreadis TG, Blow JA. An update on the potential of North American mosquitoes (Diptera: Culicidae) to transmit West Nile virus. Journal of medical entomology. 2005;42(1):57–62. pmid:15691009
- 6. Pauli G, Bauerfeind U, Blümel J, Burger R, Drosten C, Gröner A, et al. West nile virus. Transfusion medicine and hemotherapy. 2013;40(4):265. pmid:24179475
- 7.
Keeling MJ, Rohani P. Modeling infectious diseases in humans and animals. Princeton University Press; 2011.
- 8. Wonham MJ, de Camino-Beck T, Lewis MA. An epidemiological model for West Nile virus: invasion analysis and control applications. Proceedings of the Royal Society of London B: Biological Sciences. 2004;271(1538):501–507.
- 9. Liu R, Shuai J, Wu J, Zhu H. Modeling spatial spread of West Nile virus and impact of directional dispersal of birds. Mathematical Biosciences and Engineering. 2006;3(1):145. pmid:20361815
- 10. Lewis M, Rencławowicz J, Van den Driessche P. Traveling waves and spread rates for a West Nile virus model. Bulletin of mathematical biology. 2006;68(1):3–23. pmid:16794919
- 11. Durand B, Balança G, Baldet T, Chevalier V. A metapopulation model to simulate West Nile virus circulation in Western Africa, Southern Europe and the Mediterranean basin. Veterinary research. 2010;41(3):32. pmid:20167194
- 12. Maidana NA, Yang HM. Spatial spreading of West Nile Virus described by traveling waves. Journal of theoretical biology. 2009;258(3):403–417. pmid:19167405
- 13. Toni T, Welch D, Strelkowa N, Ipsen A, Stumpf MP. Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems. Journal of the Royal Society Interface. 2009;6(31):187–202.
- 14. Csilléry K, Blum MG, Gaggiotti OE, François O. Approximate Bayesian computation (ABC) in practice. Trends in ecology & evolution. 2010;25(7):410–418.
- 15. Beaumont MA. Approximate Bayesian computation in evolution and ecology. Annual review of ecology, evolution, and systematics. 2010;41:379–406.
- 16. Sunnåker M, Busetto AG, Numminen E, Corander J, Foll M, Dessimoz C. Approximate bayesian computation. PLoS computational biology. 2013;9(1):e1002803. pmid:23341757
- 17. Del Moral P, Doucet A, Jasra A. Sequential monte carlo samplers. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 2006;68(3):411–436.
- 18. Sisson SA, Fan Y, Tanaka MM. Sequential monte carlo without likelihoods. Proceedings of the National Academy of Sciences. 2007;104(6):1760–1765.
- 19.
National Centers for Environmental Information;. Available from: https://www.ncdc.noaa.gov/.
- 20.
Clements J, Schulenberg T, Iliff M, Roberson D, Fredericks T, Sullivan B, et al. The eBird/Clements checklist of birds of the world: v2015. URL: http://www.birds.cornell.edu/clementschecklist/download/IOC. 2015.
- 21. Apperson CS, Hassan HK, Harrison BA, Savage HM, Aspen SE, Farajollahi A, et al. Host feeding patterns of established and potential mosquito vectors of West Nile virus in the eastern United States. Vector-Borne and Zoonotic Diseases. 2004;4(1):71–82. pmid:15018775
- 22. Kilpatrick AM, Daszak P, Jones MJ, Marra PP, Kramer LD. Host heterogeneity dominates West Nile virus transmission. Proceedings of the Royal Society of London B: Biological Sciences. 2006;273(1599):2327–2333.
- 23. Savage HM, Aggarwal D, Apperson CS, Katholi CR, Gordon E, Hassan HK, et al. Host choice and West Nile virus infection rates in blood-fed mosquitoes, including members of the Culex pipiens complex, from Memphis and Shelby County, Tennessee, 2002–2003. Vector-Borne and Zoonotic Diseases. 2007;7(3):365–386. pmid:17767413
- 24. Hamer GL, Kitron UD, Goldberg TL, Brawn JD, Loss SR, Ruiz MO, et al. Host selection by Culex pipiens mosquitoes and West Nile virus amplification. The American journal of tropical medicine and hygiene. 2009;80(2):268–278. pmid:19190226
- 25. Cruz-Pacheco G, Esteva L, Montaõ-Hirose JA, Vargas C. Modelling the dynamics of West Nile virus. Bulletin of mathematical biology. 2005;67(6):1157. pmid:16125762
- 26. DeFelice NB, Little E, Campbell SR, Shaman J. Ensemble forecast of human West Nile virus cases and mosquito infection rates. Nature Communications. 2017;8. pmid:28233783
- 27. Sardelis MR, Turell MJ, Dohm DJ, O’Guinn ML. Vector competence of selected North American Culex and Coquillettidia mosquitoes for West Nile virus. Emerging infectious diseases. 2001;7(6):1018 pmid:11747732
- 28. Sahneh FD, Scoglio C, Van Mieghem P. Generalized epidemic mean-field model for spreading processes over multilayer complex networks. IEEE/ACM Transactions on Networking. 2013;21(5):1609–1620.
- 29. Scoglio CM, Bosca C, Riad MH, Sahneh FD, Britch SC, Cohnstaedt LW, et al. Biologically informed individual-based network model for Rift Valley fever in the US and evaluation of mitigation strategies. PloS one. 2016;11(9):e0162759. pmid:27662585
- 30. Riad MH, Scoglio CM, McVey DS, Cohnstaedt LW. An individual-level network model for a hypothetical outbreak of Japanese encephalitis in the USA. Stochastic environmental research and risk assessment. 2017;31(2):353–367.
- 31. Rabinowitz PM, Galusha D, Vegso S, Michalove J, Rinne S, Scotch M, et al. Comparison of human and animal surveillance data for H5N1 influenza A in Egypt 2006–2011. PloS one. 2012;7(9):e43851. pmid:23028474
- 32. Nasrinpour HR, Reimer AA, Friesen MR, McLeod RD. Data preparation for West Nile virus agent-based modelling: Protocol for processing bird population estimates and incorporating ArcMap in AnyLogic. JMIR research protocols. 2017;6(7). pmid:28716770
- 33. Blagrove MS, Caminade C, Waldmann E, Sutton ER, Wardeh M, Baylis M. Co-occurrence of viruses and mosquitoes at the vectors’ optimal climate range: An underestimated risk to temperate regions? PLoS neglected tropical diseases. 2017;11(6):e0005604. pmid:28617853
- 34. Rueda L, Patel K, Axtell R, Stinner R. Temperature-dependent development and survival rates of Culex quinquefasciatus and Aedes aegypti (Diptera: Culicidae). Journal of medical entomology. 1990;27(5):892–898. pmid:2231624
- 35. Reisen WK, Fang Y, Martinez VM. Effects of temperature on the transmission of West Nile virus by Culex tarsalis (Diptera: Culicidae). Journal of medical entomology. 2006;43(2):309–317 pmid:16619616
- 36. Reisen WK. Effect of temperature on Culex tarsalis (Diptera: Culicidae) from the Coachella and San Joaquin valleys of California. Journal of medical entomology. 1995;32(5):636–645. pmid:7473618
- 37. Erdos P, Rényi A. On the evolution of random graphs. Publ Math Inst Hung Acad Sci. 1960;5(1):17–60.
- 38. Barabási AL. Network science book. Network Science. 2014;625.
- 39. Ster IC, Ferguson NM. Transmission parameters of the 2001 foot and mouth epidemic in Great Britain. PLoS One. 2007;2(6):e502.
- 40. Soubeyrand S, Held L, Höhle M, Sache I. Modelling the spread in space and time of an airborne plant disease. Journal of the Royal Statistical Society: Series C (Applied Statistics). 2008;57(3):253–272.
- 41. Meyer S, Held L, et al. Power-law models for infectious disease spread. The Annals of Applied Statistics. 2014;8(3):1612–1639.
- 42. Newman ME. Power laws, Pareto distributions and Zipf’s law. Contemporary physics. 2005;46(5):323–351.
- 43. Lam TTY, Ip HS, Ghedin E, Wentworth DE, Halpin RA, Stockwell TB, et al. Migratory flyway and geographical distance are barriers to the gene flow of influenza virus among North American birds. Ecology letters. 2012;15(1):24–33. pmid:22008513
- 44. Fourment M, Darling AE, Holmes EC. The impact of migratory flyways on the spread of avian influenza virus in North America. BMC evolutionary biology. 2017;17(1):118. pmid:28545432
- 45.
Lincoln FC. Migration of birds. 16. Government Printing Office; 1999.
- 46. Tsiodras S, Kelesidis T, Kelesidis I, Bauchinger U, Falagas ME. Human infections associated with wild birds. Journal of Infection. 2008;56(2):83–98. pmid:18096237
- 47. Rappole JH, Hubalek Z. Migratory birds and West Nile virus. Journal of applied microbiology. 2003;94(s1):47–58.
- 48. Brooks-Pollock E, Roberts GO, Keeling MJ. A dynamic model of bovine tuberculosis spread and control in Great Britain. Nature. 2014;511(7508):228. pmid:25008532
- 49. Barnes CP, Silk D, Stumpf MP. Bayesian design strategies for synthetic biology. Interface focus. 2011;1(6):895–908. pmid:23226588
- 50. Toni T, Stumpf MP. Simulation-based model selection for dynamical systems in systems and population biology. Bioinformatics. 2009;26(1):104–110. pmid:19880371
- 51. Abdessalem AB, Dervilis N, Wagg D, Worden K. Model selection and parameter estimation in structural dynamics using approximate Bayesian computation. Mechanical Systems and Signal Processing. 2018;99:306–325.
- 52. Kass RE, Raftery AE. Bayes factors. Journal of the american statistical association. 1995;90(430):773–795.
- 53. Didelot X, Everitt RG, Johansen AM, Lawson DJ, et al. Likelihood-free estimation of model evidence. Bayesian analysis. 2011;6(1):49–76.
- 54. Marin JM, Pillai NS, Robert CP, Rousseau J. Relevant statistics for Bayesian model choice. Journal of the Royal Statistical Society: Series B (Statistical Methodology). 2014;76(5):833–859.
- 55. Komar N. West Nile virus: epidemiology and ecology in North America. Advances in virus research. 2003;61:185–234. pmid:14714433
- 56. Dauphin G, Zientara S. West Nile virus: recent trends in diagnosis and vaccine development. Vaccine. 2007;25(30):5563–5576. pmid:17292514
- 57. Johnson GD, Eidson M, Schmit K, Ellis A, Kulldorff M. Geographic prediction of human onset of West Nile virus using dead crow clusters: an evaluation of year 2002 data in New York State. American Journal of Epidemiology. 2005;163(2):171–180. pmid:16306307
- 58. Mundt CC, Sackett KE, Wallace LD, Cowger C, Dudley JP. Long-distance dispersal and accelerating waves of disease: empirical relationships. The American Naturalist. 2009;173(4):456–466. pmid:19249979
- 59. Kot M, Lewis MA, van den Driessche P. Dispersal data and the spread of invading organisms. Ecology. 1996;77(7):2027–2042.
- 60. Keeling MJ, Woolhouse ME, Shaw DJ, Matthews L, Chase-Topping M, Haydon DT, et al. Dynamics of the 2001 UK foot and mouth epidemic: stochastic dispersal in a heterogeneous landscape. Science. 2001;294(5543):813–817. pmid:11679661
- 61. Germann TC, Kadau K, Longini IM, Macken CA. Mitigation strategies for pandemic influenza in the United States. Proceedings of the National Academy of Sciences. 2006;103(15):5935–5940.