Abstract

While estimates of the impact of climate change on health are necessary for health care planners and climate change policy makers, models to produce quantitative estimates remain scarce. This study describes a freely available dynamic simulation model parameterized for three West Nile virus vectors, which provides an effective tool for studying vectorborne disease risk due to climate change. The Dynamic Mosquito Simulation Model is parameterized with species-specific temperature-dependent development and mortality rates. Using downscaled daily weather data, this study estimates mosquito population dynamics under current and projected future climate scenarios for multiple locations across the country. Trends in mosquito abundance were variable by location; however, an extension of the vector activity periods, and by extension disease risk, was almost uniformly observed. Importantly, midsummer decreases in abundance may be offset by shorter extrinsic incubation periods, resulting in a greater proportion of infective mosquitoes. Quantitative descriptions of the effect of temperature on the virus and mosquito are critical to developing models of future disease risk.

1. Introduction

Robust quantitative models predicting vectorborne disease changes in response to the changing climate are lacking (Hosking and Campbell-Lendrum 2012; Rodo et al. 2013). Mosquitoborne diseases are particularly vulnerable to a changing climate because of the effect of weather conditions on pathogen and mosquito development, survival, and reproduction and also on mosquito distribution and behavior [see review by Paz (2015)]. Predicting how mosquitoes will respond is challenging because of the complexity arising from interactions between mosquitoes and hosts, which may be mediated by climate, host, and mosquitoes susceptibility to infection, which varies among populations, and pathogen and vector development rates, which are influenced by weather (Hartvigsen et al. 1998; Paz and Semenza 2013).

West Nile virus disease (WNV) is now the most common domestic arthropodborne viral disease (Reimann et al. 2008). Presumably introduced from the Middle East, the first WNV case in the United States occurred at the Bronx Zoo, New York, in 1999 (Jia et al. 1999). It subsequently spread across the country [see review by Reisen (2013)]. WNV is particularly challenging as it involves multiple hosts and multiple mosquito vectors [see, e.g., the review by Reisen (2013)]. While the virus has been detected in many mammal species, the disease is primarily one in which the virus is maintained in birds through a mosquito–avian transmission cycle, and humans and horses are dead-end hosts (they can suffer from the disease but do not transmit virus to mosquitoes). The disease exhibits seasonality with most of the human cases occurring in late summer and fall. Because of the primary host of concern being humans and also because of the preferred habitat of two primary vectors (Culex pipiens and Cx. quinquefasciatus), West Nile virus disease is associated with urban areas (Reisen 2013). The other common vector, Cx. tarsalis, is more commonly associated with agriculture but, being a stronger flyer, urban areas near agriculture are generally more affected (Reisen and Reeves 1990).

Species’ distributions models, which identify abiotic and biotic factors that support vector survival, have been used to predict a climate change–driven northerly range expansion of the primary West Nile virus vectors (Chen et al. 2013; Harrigan et al. 2014; Hongoh et al. 2012). Temperature and precipitation are important underpinnings of the temporal component as seen with seasonal abundance for Cx. tarsalis (Chen et al. 2012; Reisen et al. 2008) and Cx. pipiens (Ruiz et al. 2010). By integrating landscape (spatial) and climatic (temporal) information, researchers can identify entomologic risk: when and where human exposure to potentially infected vectors is most likely (Liu and Weng 2012).

Dynamic simulation models provide an alternative for studying vectorborne diseases because they utilize a mechanistic approach (Pearce and Merletti 2006). Rather than identifying suitable ecological niches, dynamic models simulate vector abundance over time utilizing stage-dependent temperature-driven development and mortality with location-specific daily weather data. They successfully simulated observed WNV vector abundance (Gong et al. 2011; Morin and Comrie 2010).

We focus on two major WNV U.S. vectors, Cx. pipiens and Cx. tarsalis, expanding the Dynamic Mosquito Simulation Model (DyMSiM) from the original STELLA-based model for a southern U.S. WNV vector, Cx. quinquefasciatus (Morin and Comrie 2010). DyMSiM is an agent-based compartmental model that estimates daily mosquito abundance using temperature-dependent development and mortality rates as well as availability of water for the immature stages of mosquito development: egg, larval, and pupal stages. DyMSiM can be modified for most mosquito species for which temperature-dependent development data are available and for locations where appropriate weather (daily minimum and maximum temperature and precipitation) information is known. We translated DyMSiM from the original STELLA program into MATLAB to broaden the accessibility of this model and enhance our ability to manipulate and add additional modules. The findings are interpreted with respect to climate-related changes in WNV transmission.

2. Methods

2.1. Model parameters

DyMSiM simulates mosquito development using observed and simulated weather data. Full model details are provided in Morin and Comrie (2010). Here, we explain how daily mosquito development rates for each life stage were estimated using daily mean temperatures.

2.1.1. Aquatic development

For the immature aquatic (all four larval instar stages, pupal, and total immature) development, temperature-dependent rates were calculated using an updated Sharpe and DeMichele biophysical model formula provided by Wagner et al. (1984):

 
formula

where r(T) is the mean development rate at a given temperature T (kelvins); R is the universal gas constant (1.987 cal deg−1 mole−1); RHO25 is the development rate at 25°C assuming no enzymatic activity; HA, HL, and HH are the enthalpies of activation; and TL and TH are the temperatures (K) at which the rate-controlling enzyme reaches 50% activity. Details of these parameters and Statistical Analysis System (SAS) code can be found in Wagner et al. (1984). We used the Wagner et al. SAS program that fits the best form of the model (2-, 4-, or 6-parameter models) using temperature-dependent development times reported in the literature for Cx. tarsalis (Bailey and Gieke 1968; Henn et al. 2008; Mead and Conner 1987; Milby and Meyer 1986; Miura and Takahashi 1988; Rosay 1972) and Cx. pipiens (Farid 1949; Gong et al. 2011; Horsfall 1955; Kramer 1915; Lang 1963; Mead and Conner 1987; Rosay 1972; Tekle 1960). The best-fit parameter estimates are provided in Table 1. Minimum and maximum temperatures where development ceased were set at 5° and 35°C for Cx. tarsalis (Miura and Takahashi 1988) and 2.25° and 42°C for Cx. pipiens (Farid 1949).

Table 1.

Parameters for the Sharpe and DeMichele Biophysical Model amended by Wagner et al. (1984) and including the R2 measure of model fit.

Parameters for the Sharpe and DeMichele Biophysical Model amended by Wagner et al. (1984) and including the R2 measure of model fit.
Parameters for the Sharpe and DeMichele Biophysical Model amended by Wagner et al. (1984) and including the R2 measure of model fit.

2.1.2. Egg laying and development rate

The number of eggs laid per female was selected at random from a normal distribution with a mean of 160 and standard deviation of 21 for Cx. tarsalis (Moon 1976) and a mean of 255 and standard deviation of 25 for Cx. pipiens (Horsfall 1955; Vinogradova 2000). For Cx. tarsalis, the Sharpe and DeMichele model was used to estimate egg development (see Table 1; Miura and Takahashi 1988). For Cx. pipiens, a regression model was used because of the small number of data points precluding the biophysical model from converging (n = 5; R2 = 0.97); the development rate was set at

 
formula

where T is temperature in degrees Celsius. The same values for minimum and maximum temperatures where development ceased were used for egg development as for immature life stages.

2.1.3. Gonotrophic progression rate, extrinsic incubation period, and vector competence

Temperature-dependent gonotrophic development was defined by the following formulas:

 
formula

and

 
formula

where T is temperature (°C). The rate for Cx. tarsalis was reported in Hartley et al. (2012), while the rate for Cx. pipiens was derived from reported development data (Tekle 1960; n = 4 data points; R2 = 0.98). The minimum temperature to lay eggs was set to 3.67°C for Cx. tarsalis and 10°C for Cx. pipiens (Horsfall 1955). After a lag of approximately 2 days depending on the daily calculated gonotrophic progression rate, female mosquitoes were assumed to successfully find blood meals and lay eggs after completion of a temperature-dependent gonotrophic cycle.

The extrinsic incubation period (EIP) was set to occur above 14.35°C (when the EIP becomes positive) at a rate of per day for both species (estimate derived using Cx. tarsalis data; Reisen et al. 2006). Using temperature-dependent EIP and the gonotrophic cycle, we calculated the number of days to infectiousness for each cohort of mosquitoes. We then calculated the proportion of infected vectors.

Vector competence, the proportion of mosquitoes fed an infectious blood meal that become infected (detectable West Nile virus in the saliva), was set to 0.79 for both species (Goddard et al. 2002; Turell et al. 2001; Turell et al. 2002).

2.1.4. Mortality rates

During each life stage, mortality occurs, removing individuals from the population. For eggs and adults, the rates were temperature independent: 5% day−1 for eggs (Eisenberg et al. 1995) with no thresholds for Cx. tarsalis found in the literature and above 2° and below 42°C (Farid 1949) for Cx. pipiens; between 0° and 40°C, adults died at a rate of 16.6% day−1 for Cx. tarsalis (Mahmood et al. 2004; Moon 1976; Reisen et al. 1995) and 22% day−1 for Cx. pipiens (Gong et al. 2011; Tekle 1960).

The amount of permanent water was set by the user. Daily precipitation data are used to increase the amount of water available in the system, which then evaporates at a fixed rate. Ideally, the model can be amended to account for permanent water sources (accounting for human behavior like excess watering, unmaintained pools or fountains, and other sources like ponds and sewers, which are important habitats for Cx. pipiens), semipermanent sources (such as irrigation, which is important habitat for Cx. tarsalis), and containers that fill with precipitation and evaporate based on temperature and hours of daylight. In this version of the model, water levels are important for the calculation of carrying capacity [see Morin and Comrie (2010) for detail] and do not overflow and washout immature life stages with extreme events (Koenraadt and Harrington 2008).

2.1.5. Weather data

Long-term daily temperature and precipitation data from weather stations were downloaded from the U.S. Historical Climatology Network. Stations were selected where daily temperature and precipitation data were available for at least 95% of each year for the period 1970–2000 and include some diversity with respect to the range of each mosquito’s distribution [based on Darsie and Ward (2005); recreated as Figure 1]. This high threshold on the availability of meteorological data for urban areas within the current ranges of Cx. pipiens and Cx. tarsalis yielded four (Allentown, Pennsylvania; Groton, Connecticut; Fort Collins, Colorado; and Urbana, Illinois) and five (Chandler Heights, Arizona; El Paso, Texas; Fort Collins, Colorado; Redlands, California; and Urbana, Illinois) sites, respectively (starred locations in Figure 1).

Figure 1.

Comparison of simulated (a) Culex pipiens and (b) Culex tarsalis mosquito populations. Shaded areas show the U.S. distribution of the two species [redrawn from Darsie and Ward (2005)], and the black stars indicate locations from which weather data were acquired and the abundance simulated (for Cx. pipiens: Groton, Connecticut; Allentown, Pennsylvania; Urbana, Illinois; and Fort Collins, Colorado; for Cx. tarsalis: Urbana, Illinois; Fort Collins, Colorado; El Paso, Texas; Chandler, Arizona; Redlands, California). The left column of graphs for each location shows MC (2045–65) and the right column shows EOC (2080–99) estimates. The upper graph of each couplet shows a scaled estimate of the number of mosquitoes per trap night and the lower graph shows the percentage of infective mosquitoes. In all images, black dashed lines provide current estimates, the solid blue lines are SRA2 scenarios and dotted red lines are SRB1 scenarios.

Figure 1.

Comparison of simulated (a) Culex pipiens and (b) Culex tarsalis mosquito populations. Shaded areas show the U.S. distribution of the two species [redrawn from Darsie and Ward (2005)], and the black stars indicate locations from which weather data were acquired and the abundance simulated (for Cx. pipiens: Groton, Connecticut; Allentown, Pennsylvania; Urbana, Illinois; and Fort Collins, Colorado; for Cx. tarsalis: Urbana, Illinois; Fort Collins, Colorado; El Paso, Texas; Chandler, Arizona; Redlands, California). The left column of graphs for each location shows MC (2045–65) and the right column shows EOC (2080–99) estimates. The upper graph of each couplet shows a scaled estimate of the number of mosquitoes per trap night and the lower graph shows the percentage of infective mosquitoes. In all images, black dashed lines provide current estimates, the solid blue lines are SRA2 scenarios and dotted red lines are SRB1 scenarios.

Future weather data were generated using the Long Ashton Research Station Weather Generator (LARS-WG), version 5, a stochastic weather generator that utilizes current daily weather data to generate a time series of daily weather data with the same statistical properties as the recorded weather station data (Semenov and Stratonovitch 2010). Of the 14 available global climate models used in the International Panel on Climate Change (IPCC) Fourth Assessment Report1 incorporated into LARS-WG, we chose the National Center for Atmospheric Research’s Community Climate System Model, version 3 (CCSM), to generate a synthetic current baseline, midcentury (MC: 2045–65) and end of century (EOC: 2080–99) estimates under two projected climate scenarios—the IPCC’s SRA2 and SRB1 scenarios. We selected CCSM because this high-resolution (1.4° × 1.4°) model incorporates radiation, boundary physics, and precipitation physics with their atmospheric component and a land surface model incorporating soil and vegetation.

While there is some question regarding whether to use single models, multiple models for comparison, or means of multiple models, Winter and Nychka (2010) showed that model averages perform less well under situations where individual forecasts generally agree. To test whether the differences among the GCMs would influence the interpretation of our simulated mosquito abundances, we compared a subset of the results based on the CCSM to the United Kingdom’s HadCM3 GCMs for the midcentury and end of century using the original STELLA version of the model.

The two future climate scenarios (SRA2 and SRB1) represent the range of future climate assumptions (Bernstein et al. 2007). SRA2 assumes high human population growth with slow carbon dioxide reductions. Under this scenario, mean global temperatures are assumed to rise by 3.4°C from 1980 to 1999 to 2090 to 2099. SRB1 assumes the world population peaks in the midcentury with rapid changes in economic structure that reduce carbon dioxide emissions. In this scenario, the average global temperature increase over the next century is 1.8°C.

At each of the four Cx. pipiens and five Cx. tarsalis locations listed above for each time period (midcentury and end of century), 30 simulations of the 30 years of simulated data were run for each species and climate change scenario combination. A double-pass filter was applied to smooth the data, and a weighted scale factor was applied to scale the estimate to the observed data.

2.2. Model validation

Model accuracy was assessed using Colorado Mosquito Control, Inc.’s (CMC), 2003–09 surveillance data (Cx. tarsalis and Cx. pipiens: 34.8 trap nights per week, standard deviation 12.3, range 1–50) and the Connecticut Agricultural Experiment Station’s (CAES) 2000–10 surveillance data (Cx. pipiens: 1.7 traps per week, range 1–4). CMC’s surveillance data were from Fort Collins (40°33′33″N, 105°4′41″W), an urban area in north-central Colorado with the Rocky Mountain foothills to the west and irrigated agriculture to the east. CMC’s municipal surveillance program was established in response to the incursion of WNV into the region and serves a dual purpose of monitoring infection rates and mosquito abundance on a temporal scale. The CAES has been collecting mosquitoes for arbovirus surveillance in 91 locations across the state since the 1990s. The CAES data used here were from Bridgeport (41°11′11″N, 73°11′44″W), an urban area (Cx. pipiens habitat) that also had a complete set of meteorological data available. Ideally, additional locations could have been validated, but mosquito surveillance data consistently collected at a specific location (rather than in response to complaints) over long periods of time are rare.

Data analysis

Multiple comparisons of the abundance estimations between the observed and simulated data were used. Because the magnitude of Pearson’s correlation is not well defined with respect to the accuracy of a model prediction, we also reported Willmott’s d and root-mean-square error (Willmott et al. 1985). These latter indices consider the difference between observed and predicted values, where RMSE is low and Willmott’s d approaches 1 when the predicted values match the observed data.

To compare an epidemiologically important level of abundance, we dichotomized weeks at greater than 10 mosquitoes per trap night (TA from experience). We then visually compared the number of weeks across the different species–scenario–time period combinations.

Vectorial capacity C provides an overview of the components necessary for mosquitoborne disease transmission:

 
formula

where m is the density of mosquito to hosts, a is the probability a vector feeds on the host, P is the daily survival, and n is the EIP—the number of days it takes for the mosquito to become infectious (Macdonald 1957).

For a mosquito to transmit disease, it must take an infectious blood meal, survive the period of time it takes for the pathogen to replicate and infect the mosquito (EIP), and the mosquito must successfully find a second, susceptible host. Thus, mosquitoborne disease transmission is related to the number of mosquitoes m surviving the infectious period Pn and feeding on susceptible hosts a2. These are mediated by the number of infected hosts, host feeding preferences, susceptibility of mosquitoes and host to infection, and host behavior (e.g., birds eating mosquitoes or humans using repellant). These, in turn, are mediated by environmental factors including the influence of weather on host (e.g., birds taking flight or humans remaining indoors) and mosquito behavior, the availability of habitat for mosquitoes, and weather-related effects on mosquito survival and development. To address multiple components of transmission while restricting ourselves to those parameters for which reliable data were available, we focused on changes in 1) mosquito abundance (related to m), 2) the period of time (weeks of the year) when the mosquitoes were active (related to m), and 3) the fraction of infectious mosquitoes (related to Pn).

Notably missing in this model are hosts. This limits our interpretation to changes in entomologic risk rather than changes in disease incidence. However, at this time, good estimates for host abundance (e.g., species-specific bird populations), susceptibility across communities (both avian and human), and behavior (e.g., human mosquito repellant use or avian host aggregation) are unavailable.

3. Results

3.1. Model validation

Similar to a mosquito model for the northeastern United States by Gong et al. (2011), our model was able to simulate the observed interannual variability in vector populations and captures the peaks and troughs in annual mosquito abundance. The Pearson’s correlation coefficient for Cx. pipiens was low but significant for the single Connecticut site (Pearson’s r = 0.34, P < 0.001), and the model comparisons indicated a good fit (RMSE and Willmott’s d) between the simulation population and the observed data (Table 2).

Table 2.

Comparison between median numbers of mosquitoes per trap night averaged across the weeks of the available data. Observed and double pass–smoothed estimated weekly data [measures of agreement performed on log10(+1) transformed data] excluding missing data (i.e. no winter).

Comparison between median numbers of mosquitoes per trap night averaged across the weeks of the available data. Observed and double pass–smoothed estimated weekly data [measures of agreement performed on log10(+1) transformed data] excluding missing data (i.e. no winter).
Comparison between median numbers of mosquitoes per trap night averaged across the weeks of the available data. Observed and double pass–smoothed estimated weekly data [measures of agreement performed on log10(+1) transformed data] excluding missing data (i.e. no winter).

While light traps are the standard for Culex surveillance, they are limited in the information they provide. Discrepancies between observed and predicted abundance are common in models like these (Gong et al. 2011; Morin and Comrie 2010), but scaling the number of mosquitoes used to seed the model helped to narrow the difference. Moreover, while the light traps used for mosquito surveillance are a good estimate of relative vector abundance for the Culex species, they are limited in the physical area they sample (Brown et al. 2008). Thus, while the magnitude may be scaled to match the simulated to the observed, similarity in the seasonal patterns is more critical.

3.2. Projected changes in vector abundance

Changes in mosquito abundance were localized, with abundance increasing (e.g., Fort Collins, Colorado, for both species; Figure 1, top row of graphic) or decreasing (e.g., Chandler Heights, Arizona; Figure 1b) depending on the location. Generally, there was a latitudinal gradient where higher latitudes experienced increased abundance, while lower latitudes experienced a decrease in midsummer abundance. Warming in cool regions may be exhibiting a positive effect on mosquito development but, at the upper edges of thermal tolerance, heating had a negative effect. The most dramatic decrease is predicted for Cx. tarsalis in Chandler Heights, Arizona, where for all future scenarios midsummer abundance was reduced to near zero.

3.3. Projected changes in vector phenology

Comparing the 30 simulated years for each scenario (SRA2 and SRB1) and time period (midcentury and end of century), we showed a nearly uniform increase in the number of weeks with mosquito activity. The extension of the season is evident in the annual average abundance (top row of the Figure 1 graphs) as well as in the epidemiologically important increase in number of weeks with more than 10 mosquitoes per trap nights (Figure 2).

Figure 2.

Number of weeks with increased vector activity, defined as weeks where the scaled estimates of the number of mosquitoes per trap night are greater than 10.

Figure 2.

Number of weeks with increased vector activity, defined as weeks where the scaled estimates of the number of mosquitoes per trap night are greater than 10.

3.4. Projected changes in infectious vectors

Using temperature-dependent EIP and the gonotrophic cycle, we calculated the number of days to infectiousness for each cohort. We then calculated the proportion of infectious mosquitoes (Figure 1, lower panel of each graphic couplet). This provided nuance to the estimated changes in mosquito abundance; while populations may decrease with higher midsummer temperatures, because the number of bites to transmission would be reduced, the proportion of infectious mosquitoes increased. In Chandler Heights, Arizona, estimates of the infectious proportion reached near 40% under the end of century estimates compared with just 5% and ~3% for Urbana, Illinois (peak SRA2 and SRB1 end of century estimates, respectively).

3.5. Global climate model comparison

Using the STELLA version of the model, we compared the Cx. tarsalis abundance for Urbana, Illinois, based on the CCSM and the United Kingdom’s HadCM3 GCMs for the midcentury and end of century (results not shown). We selected Cx. trasalis in Urbana, Illinois, because this combination showed an increase in midcentury predictions but a decrease in midsummer abundance for the end of century simulations. There was no statistically significant difference in either the total annual abundance or the late summer (critical mosquito period) estimated by either model. A similar midsummer reduction in the number of mosquitoes is observed, though this decrease is more pronounced in the HadCM model. Both models predict a similar period of mosquito activity.

4. Discussion

We present a model parameterized using previously published field and laboratory data, which estimates mosquito population dynamics in response to daily weather data at selected sites across each mosquito species’ distribution. Like others, while the estimated number of vectors had to be scaled to numbers observed in light traps, we found a significant association between observed Culex mosquito population dynamics and ground level weather data (Gong et al. 2011; Morin and Comrie 2010). We add three important findings with respect to future climate and vectorborne disease to this body of literature: 1) climate-related changes in vector population dynamics will vary by location; 2) changes will occur in both vector abundance and length of the vector activity period; and 3) beyond changes in vector abundance, survival to transmission is important to predicting risk.

4.1. Changes in vector population dynamics will vary by location

Our model shows some areas (generally the higher latitudes: Colorado, Connecticut, and Illinois) will warm to support increased vector abundance, while others (generally currently hot areas: El Paso, Texas, and Chandler Heights, Arizona) may become too hot, thus negatively affecting mosquito development and survival. Based on the association between the probability of vector presence and current WNV cases at the county level, Harrigan et al. (2014) projected WNV expansion due to expansion of suitable climate. Their findings are synergistic to our city level modeling of seasonal changes, where the extension of the season and the increased abundance at more northerly latitudes likely translates to greater risk of West Nile virus disease. The midsummer reductions in abundance we observed in lower latitudes may translate to decreased mosquitoes activity and thus reduced disease risk.

The WNV vectors studied here are limited both in their presences and in their seasonal activity by elevation (Barker et al. 2010; Schurich et al. 2014). While the analysis presented here did not allow us to specifically evaluate vector distribution across elevations, we can expect that, just as cooler northerly latitudes may become more permissible, higher elevations may also become more permissible and conversely some of the lower elevations may become less permissible.

4.2. Changes will occur in both vector abundance and length of the vector activity period

Comparing future climate scenarios with the current vector season, we show mosquitoes emerging earlier in the year and remaining active longer into the fall (Figures 1, 2). This is important in vectorborne diseases like WNV that have a period of amplification in their avian reservoir populations that may be followed by increases in human biting rates (Kilpatrick et al. 2006). The late summer is a critical period in WNV transmission when either through host switching (Kilpatrick et al. 2006) or through changes in the prevalence of infectious mosquitoes (Hamer et al. 2009), human cases are reported. Warmer temperatures have been associated with earlier onset of host-seeking and feeding activity among Cx. tarsalis; however, in 1 year of a study, but greater WNV activity was not observed (Reisen et al. 2010).

4.3. Survival to transmission is important to predicting risk

Our study highlights the importance of considering not just changes in mosquito abundance but changes in viral replication as well in order to get closer to a measure of human risk. We do not directly calculate vectorial capacity C because of insufficient estimates of virus transmission efficiencies and true avian and human biting rates. Nonetheless, the data presented address two components of vectorial capacity—the host biting rate ma and the fraction of infectious mosquitoes (which requires survival through the extrinsic incubation period). Our measures of mosquito abundance can be considered a surrogate for ma, the host biting rate (Dye 1986). Only mosquitoes surviving the EIP are able to transmit pathogens (Gubler et al. 2001), thus mosquito survival beyond the extrinsic incubation period Pn is an important component in determining whether disease will persist in a host population (Novoseltsev et al. 2012). It should be noted that while studies have shown a reduction in fecundity due to infection with WNV, mosquito survival was not affected (Styer et al. 2007).

Changes in abundance must be considered in conjunction with expected changes in the EIP. Hartley et al. (2012) found that, though survival decreased with warming, the effect of warming on the EIP and the gonotrophic cycle yielded a greater proportion of females surviving to transmission. This effect of shortening the EIP with increasing temperature was shown in vivo for WNV in both Cx. pipiens (Dohm et al. 2002; Kilpatrick et al. 2008) and Cx. tarsalis (Reisen et al. 2006). Similarly, our simulation suggests that the shortening of the EIP and gonotrophic cycle had a positive effect on the proportion of mosquitoes surviving to transmission. For example, in Chandler Heights, Arizona, while there may be a drop in mosquito abundance with the increased temperatures, the decrease in time to infectiousness results in a greater proportion of mosquitoes surviving the incubation period.

4.4. Limitations of results and recommended uses

Two important limitations of this model with respect to investigating the effect of climate change on WNV vector distribution include 1) the automatic recovery feature, where if the population drops below one, a single mosquito is returned to the pool of adults thus maintaining the system; and 2) the model has no evolutionary adaptation to new temperatures.

Recovery is incorporated into the model to maintain the population over the winter when mosquitoes in cooler climes either diapause because of the shortened day length or seek warmer locations and cease feeding and host-seeking activity. During this period, surveillance trapping yields few to no mosquitoes (Bolling et al. 2007; Spielman 2001). Recovery prevents the simulated populations from reaching zero despite inhospitable conditions and thus mosquito populations rebound when the conditions are again permissive (Morin and Comrie 2010). This mathematically necessary recovery is similar to biological diapause (effectively allowing mosquitoes to survive the winter and reappear in the spring) and has little effect on the calculated averages as the time at which mosquitoes are reintroduced varies between years. It could also reflect reinitialization of populations due to immigration. The current version of the model lets populations survive through the winter in permissible regions; simulations in regions that become prohibitive due to reduced precipitation would appear as successions of recovery and crashes and can therefore be detected.

An organism’s adaptation to the changing climate will influence its future abundance and distribution. Empirical data on the adaptive capacity of mosquitoes are limited. We know of one example where researchers investigated climate-related evolutionary changes in egg desiccation, a limiting factor for Aedes aegypti distribution. Relying on observed evolutionary changes in egg desiccation resistance in a model system (Drosophila), researchers showed adaptation in desiccation resistance of Ae. aegypti eggs may have led to the establishment of this important disease vector in heavily populated areas of Australia (Kearney et al. 2009). Bradshaw and Holzapfel (2006) report mosquito populations exhibiting rapid evolutionary shifts in the timing of diapause to take advantage of temperature-induced changes in season length. We know of no other models of disease vectors that incorporate evolutionary adaptation to a changing climate. Models will need to include adaptation to temperature, changes in seasonality, desiccation, and even snowmelt runoff for some species.

5. Conclusions

We expanded DyMSiM to incorporate two additional WNV vectors as well as mosquito competence and temperature-dependent EIP, yielding a quantitative tool for comparing population dynamics that vary based on meteorological data. This analysis provided three important epidemiological findings: 1) changes in mosquito population dynamics will vary by location, 2) mosquito activity periods are generally expected to increase across sites in the northern latitudes of the United States, and 3) the interaction between temperature on both the virus and the mosquito will be critical to understanding future disease risk.

While this model provides important insight to WNV transmission in the United States, it does not provide a future WNV risk map. Many factors beyond mosquito biology influence whether a disease occurs (Randolph and Rogers 2010). This model focuses on disease vectors, ignoring human behavior, which modifies disease risk. Engaging in positive public health measures (e.g., mosquito control, removal of potential breeding sites around the home, protecting water storage sources, utilization of repellants and protective clothing, being aware of host-seeking vector habits) has an influence on susceptibility and disease occurrence (Reiter et al. 2003). When extrapolating this model or building upon it to include disease risk, human behavior will need to be parameterized and integrated. Addressing the gap in quantitative assessments of the effects climate change will have on mosquitoborne disease requires integration of mosquito, viral, and host parameters including how they interact and how they will adapt to future climates (LaDeau et al. 2011).

Acknowledgments

The authors thank the individuals involved in the mosquito data collection at The CT Agricultural Experiment Station and at Colorado Mosquito Control, Inc. We thank Cory Morin for the development of the STELLA version of DyMSiM. Support for ACC and HEB came in part from a National Oceanic and Atmospheric Administration (NOAA) Regional Integrated Sciences and Assessments (RISA) program via the Climate Assessment for the Southwest (CLIMAS) program at the University of Arizona. Support for AY was partly provided by The University of Arizona Center for Insect Science and by the National Science Foundation under Award 0841234. The project described was supported by Grant Number K01AI101224 from the National Institute of Allergy and Infectious Diseases. Its contents are solely the responsibility of the authors and do not necessarily represent the official views of the National Institute of Allergy and Infectious Diseases.

References

References
Bailey
,
S. F.
, and
P. A.
Gieke
,
1968
: A study of the effect of water temperatures on rice field mosquito development. Proc. and Papers of the 36th Annual Conf. of the California Mosquito Control Association, Fresno, CA, California Mosquito Control Association,
53
61
.
Barker
,
C. M.
,
B. F.
Eldridge
, and
W. K.
Reisen
,
2010
:
Seasonal abundance of Culex tarsalis and Culex pipiens complex mosquitoes (Diptera: Culicidae) in California
.
J. Med. Entomol.
,
47
,
759
768
, doi:.
Bernstein
,
L.
, and Coauthors
,
2007
: Climate Change 2007: Synthesis Report. Cambridge University Press, 73 pp.
Bolling
,
B. G.
,
C. G.
Moore
,
S. L.
Anderson
,
C. D.
Blair
, and
B. J.
Beaty
,
2007
:
Entomological studies along the Colorado Front Range during a period of intense West Nile virus activity
.
J. Amer. Mosq. Control Assoc.
,
23
,
37
46
, doi:.
Bradshaw
,
W. E.
, and
C. M.
Holzapfel
,
2006
:
Evolutionary response to rapid climate change
.
Science
,
312
,
1477
1478
, doi:.
Brown
,
H. E.
,
M.
Paladini
,
R. A.
Cook
,
D.
Kline
,
D.
Barnard
, and
D.
Fish
,
2008
:
Effectiveness of mosquito traps in measuring species abundance and composition
.
J. Med. Entomol.
,
45
,
517
521
, doi:.
Chen
,
C. C.
,
T.
Epp
,
E.
Jenkins
,
C.
Waldner
,
P. S.
Curry
, and
C.
Soos
,
2012
:
Predicting weekly variation of Culex tarsalis (Diptera: Culicidae) West Nile virus infection in a newly endemic region, the Canadian prairies
.
J. Med. Entomol.
,
49
,
1144
1153
, doi:.
Chen
,
C. C.
,
E.
Jenkins
,
T.
Epp
,
C.
Waldner
,
P. S.
Curry
, and
C.
Soos
,
2013
:
Climate change and West Nile virus in a highly endemic region of North America
.
Int. J. Environ. Res. Public Health
,
10
,
3052
3071
, doi:.
Darsie
,
R. F.
, and
R. A.
Ward
,
2005
: Identification and Geographical Distribution of the Mosquitoes of North America, North of Mexico. University of Florida Press, 383 pp.
Dohm
,
D. J.
,
M. L.
O’Guinn
, and
M. J.
Turell
,
2002
:
Effect of environmental temperature on the ability of Culex pipiens (Diptera: Culicidae) to transmit West Nile virus
.
J. Med. Entomol.
,
39
,
221
225
, doi:.
Dye
,
C.
,
1986
:
Vectorial capacity: Must we measure all its components?
Parasitol. Today
,
2
,
203
209
, doi:.
Eisenberg
,
J. N.
,
W. K.
Reisen
, and
R. C.
Spear
,
1995
:
Dynamic model comparing the bionomics of two isolated Culex tarsalis (Diptera: Culicidae) populations: Model development
.
J. Med. Entomol.
,
32
,
83
97
, doi:.
Farid
,
M. A.
,
1949
:
Relationships between certain populations of Culex pipiens Linnaeus and Culex quinquefasciatus say in the United States
.
Amer. J. Epidemiol.
,
49
,
83
100
.
Goddard
,
L. B.
,
A. E.
Roth
,
W. K.
Reisen
, and
T. W.
Scott
,
2002
:
Vector competence of California mosquitoes for West Nile virus
.
Emerging Infect. Dis.
,
8
,
1385
1391
, doi:.
Gong
,
H.
,
A. T.
DeGaetano
, and
L. C.
Harrington
,
2011
:
Climate-based models for West Nile Culex mosquito vectors in the northeastern US
.
Int. J. Biometeor.
,
55
,
435
446
, doi:.
Gubler
,
D. J.
,
P.
Reiter
,
K. L.
Ebi
,
W.
Yap
,
R.
Nasci
, and
J. A.
Patz
,
2001
:
Climate variability and change in the United States: Potential impacts on vector- and rodent-borne diseases
.
Environ. Health Perspect.
,
109
,
223
233
, doi:.
Hamer
,
G. L.
,
U. D.
Kitron
,
T. L.
Goldberg
,
J. D.
Brawn
,
S. R.
Loss
,
M. O.
Ruiz
,
D. B.
Hayes
, and
E. D.
Walker
,
2009
:
Host selection by Culex pipiens mosquitoes and West Nile virus amplification
.
Amer. J. Trop. Med. Hyg.
,
80
,
268
278
.
Harrigan
,
R. J.
,
H. A.
Thomassen
,
W.
Buermann
, and
T. B.
Smith
,
2014
:
A continental risk assessment of West Nile virus under climate change
.
Global Change Biol.
,
20
,
2417
2425
, doi:.
Hartley
,
D. M.
,
C. M.
Barker
,
A.
Le Menach
,
T.
Niu
,
H. D.
Gaff
, and
W. K.
Reisen
,
2012
:
Effects of temperature on emergence and seasonality of West Nile virus in California
.
Amer. J. Trop. Med. Hyg.
,
86
,
884
894
, doi:.
Hartvigsen
,
G.
,
A.
Kinzig
, and
G.
Peterson
,
1998
:
Use and analysis of complex adaptive systems in ecosystem science: Overview of special section
.
Ecosystems
,
1
,
427
430
, doi:.
Henn
,
J. B.
,
M. E.
Metzger
,
J. A.
Kwan
,
J. E.
Harbison
,
C. L.
Fritz
,
J.
Riggs-Nagy
,
M.
Shindelbower
, and
V. L.
Kramer
,
2008
:
Development time of Culex mosquitoes in stormwater management structures in California
.
J. Amer. Mosq. Control Assoc.
,
24
,
90
97
, doi:.
Hongoh
,
V.
,
L.
Berrang-Ford
,
M. E.
Scott
, and
L. R.
Lindsay
,
2012
:
Expanding geographical distribution of the mosquito, Culex pipiens, in Canada under climate change
.
Appl. Geogr.
,
33
,
53
62
, doi:.
Horsfall
,
W. R.
,
1955
: Mosquitoes: Their Bionomics and Relation to Disease. Ronald Press, 723 pp.
Hosking
,
J.
, and
D.
Campbell-Lendrum
,
2012
:
How well does climate change and human health research match the demands of policymakers? A scoping review
.
Environ. Health Perspect.
,
120
,
1076
1082
, doi:.
Jia
,
X. Y.
, and Coauthors
,
1999
:
Genetic analysis of West Nile New York 1999 encephalitis virus
.
Lancet
,
354
,
1971
1972
.
Kearney
,
M.
,
W. P.
Porter
,
C.
Williams
,
S.
Ritchie
, and
A. A.
Hoffmann
,
2009
:
Integrating biophysical models and evolutionary theory to predict climatic impacts on species’ ranges: The dengue mosquito Aedes aegypti in Australia
.
Funct. Ecol.
,
23
,
528
538
, doi:.
Kilpatrick
,
A. M.
,
L. D.
Kramer
,
M. J.
Jones
,
P. P.
Marra
, and
P.
Daszak
,
2006
:
West Nile virus epidemics in North America are driven by shifts in mosquito feeding behavior
.
PLoS Biol.
,
4
,
606
610
, doi:.
Kilpatrick
,
A. M.
,
M. A.
Meola
,
R. M.
Moudy
, and
L. D.
Kramer
,
2008
:
Temperature, viral genetics, and the transmission of West Nile virus by Culex pipiens mosquitoes
.
PLoS Pathog.
,
4
,
e1000092
, doi:.
Koenraadt
,
C. J.
, and
L. C.
Harrington
,
2008
:
Flushing effect of rain on container-inhabiting mosquitoes Aedes aegypti and Culex pipiens (Diptera: Culicidae)
.
J. Med. Entomol.
,
45
,
28
35
, doi:.
Kramer
,
S. D.
,
1915
:
The effect of temperature on the life cycle of Musca domestica and Culex pipiens
.
Science
,
41
,
874
877
, doi:.
LaDeau
,
S. L.
,
G. E.
Glass
,
N. T.
Hobbs
,
A.
Latimer
, and
R. S.
Ostfeld
,
2011
:
Data-model fusion to better understand emerging pathogens and improve infectious disease forecasting
.
Ecol. Appl.
,
21
,
1443
1460
, doi:.
Lang
,
C. A.
,
1963
:
The effect of temperature on the growth and chemical composition of the mosquito
.
J. Insect Physiol.
,
9
,
279
286
, doi:.
Liu
,
H.
, and
Q.
Weng
,
2012
:
Environmental factors and risk areas of West Nile virus in southern California, 2007-2009
.
Environ. Model. Assess.
,
17
,
441
452
, doi:.
Macdonald
,
G.
,
1957
: The Epidemiology and Control of Malaria. Oxford University Press, 201 pp.
Mahmood
,
F.
,
W. K.
Reisen
,
R. E.
Chiles
, and
Y.
Fang
,
2004
:
Western equine encephalomyelitis virus infection affects the life table characteristics of Culex tarsalis (Diptera: Culicidae)
.
J. Med. Entomol.
,
41
,
982
986
, doi:.
Mead
,
S. S.
, and
G. E.
Conner
,
1987
: Temperature-related growth and mortality rates. Proc. and Papers of the 55th Annual Conf. of the California Mosquito Control Association, Fresno, CA, California Mosquito Control Association,
133
137
.
Milby
,
M. M.
, and
R. P.
Meyer
,
1986
:
The influence of constant versus fluctuating water temperatures on the preimaginal development of Culex tarsalis
.
J. Amer. Mosq. Control Assoc.
,
2
,
7
10
.
Miura
,
T.
, and
R. M.
Takahashi
,
1988
: Development and survival rates of immature stages of Culex tarsalis (Diptera: Culicidae) in central California rice fields. Proc. and Papers of the 56th Annual Conf. of the California Mosquito Control Association, San Mateo, CA, California Mosquito Control Association,
168
179
.
Moon
,
T. E.
,
1976
:
A statistical model of the dynamics of a mosquito vector (Culex tarsalis) population
.
Biometrics
,
32
,
355
368
, doi:.
Morin
,
C. W.
, and
A. C.
Comrie
,
2010
:
Modeled response of the West Nile virus vector Culex quinquefasciatus to changing climate using the dynamic mosquito simulation model
.
Int. J. Biometeor.
,
54
,
517
529
, doi:.
Novoseltsev
,
V. N.
,
A. I.
Michalski
,
J. A.
Novoseltseva
,
A. I.
Yashin
,
J. R.
Carey
, and
A. M.
Ellis
,
2012
:
An age-structured extension to the vectorial capacity model
.
PLoS One
,
7
,
e39479
, doi:.
Paz
,
S.
,
2015
:
Climate change impacts on West Nile virus transmission in a global context
.
Philos. Trans. Roy. Soc. London
,
B370
,
20130561
, doi:.
Paz
,
S.
, and
J. C.
Semenza
,
2013
:
Environmental drivers of West Nile fever epidemiology in Europe and western Asia—A review
.
Int. J. Environ. Res. Public Health
,
10
,
3543
3562
, doi:.
Pearce
,
N.
, and
F.
Merletti
,
2006
:
Complexity, simplicity, and epidemiology
.
Int. J. Epidemiol.
,
35
,
515
519
, doi:.
Randolph
,
S. E.
, and
D. J.
Rogers
,
2010
:
The arrival, establishment and spread of exotic diseases: Patterns and predictions
.
Nat. Rev. Microbiol.
,
8
,
361
371
, doi:.
Reimann
,
C. A.
,
E. B.
Hayes
,
C.
DiGuiseppi
,
R.
Hoffman
,
J. A.
Lehman
,
N. P.
Lindsey
,
G. L.
Campbell
, and
M.
Fischer
,
2008
:
Epidemiology of neuroinvasive arboviral disease in the United States, 1999–2007
.
Amer. J. Trop. Med. Hyg.
,
79
,
974
979
.
Reisen
,
W. K.
,
2013
:
Ecology of West Nile virus in North America
.
Viruses
,
5
,
2079
2105
, doi:.
Reisen
,
W. K.
, and
W. C.
Reeves
,
1990
: Bionomics and ecology of Culex tarsalis and other potential mosquito vector species. Epidemiology and Control of Mosquito Borne Arboviruses in California, 1943-1987, W. C. Reeves et al., Ed., California Mosquito and Vector Control Association, 254–329.
Reisen
,
W. K.
,
H. D.
Lothrop
, and
J. L.
Hardy
,
1995
:
Bionomics of Culex tarsalis (Diptera: Culicidae) in relation to arbovirus transmission in southeastern California
.
J. Med. Entomol.
,
32
,
316
327
, doi:.
Reisen
,
W. K.
,
Y.
Fang
, and
V. M.
Martinez
,
2006
:
Effects of temperature on the transmission of West Nile virus by Culex tarsalis (Diptera: Culicidae)
.
J. Med. Entomol.
,
43
,
309
317
, doi:.
Reisen
,
W. K.
,
D.
Cayan
,
M.
Tyree
,
C. A.
Barker
,
B.
Eldridge
, and
M.
Dettinger
,
2008
:
Impact of climate variation on mosquito abundance in California
.
J. Vector Ecol.
,
33
,
89
98
, doi:.
Reisen
,
W. K.
,
T.
Thiemann
,
C. M.
Barker
,
H.
Lu
,
B.
Carroll
,
Y.
Fang
, and
H. D.
Lothrop
,
2010
:
Effects of warm winter temperature on the abundance and gonotrophic activity of Culex (Diptera: Culicidae) in California
.
J. Med. Entomol.
,
47
,
230
237
, doi:.
Reiter
,
P.
, and Coauthors
,
2003
:
Texas lifestyle limits transmission of dengue virus
.
Emerg. Infect. Dis.
,
9
,
86
89
, doi:.
Rodo
,
X.
, and Coauthors
,
2013
:
Climate change and infectious diseases: Can we meet the needs for better prediction?
Climatic Change
,
118
,
625
640
, doi:.
Rosay
,
B.
,
1972
: Comparative growth rates of aquatic stages of ten mosquito species (Diptera: Culicidae) at two constant temperatures. Proc. and Papers of the 25th Annual Meeting of the Utah Mosquito Abatement Association, Salt Lake City, UT, Utah Mosquito Abatement Association, 31–44.
Ruiz
,
M. O.
, and Coauthors
,
2010
:
Local impact of temperature and precipitation on West Nile virus infection in Culex species mosquitoes in northeast Illinois, USA
.
Parasites Vectors
,
3
,
19
, doi:.
Schurich
,
J. A.
,
S.
Kumar
,
L.
Eisen
, and
C. G.
Moore
,
2014
:
Modeling Culex tarsalis abundance on the northern Colorado Front Range using a landscape-level approach
.
J. Amer. Mosq. Control Assoc.
,
30
,
7
20
, doi:.
Semenov
,
M. A.
, and
P.
Stratonovitch
,
2010
:
Use of multi-model ensembles from global climate models for assessment of climate change impacts
.
Climate Res.
,
41
,
1
14
, doi:.
Spielman
,
A.
,
2001
:
Structure and seasonality of Nearctic Culex pipiens populations
.
Ann. N. Y. Acad. Sci.
,
951
,
220
234
, doi:.
Styer
,
L. M.
,
M. A.
Meola
, and
L. D.
Kramer
,
2007
:
West Nile virus infection decreases fecundity of Culex tarsalis females
.
J. Med. Entomol.
,
44
,
1074
1085
, doi:.
Tekle
,
A.
,
1960
:
The physiology of hibernation and its role in the geographical distribution of populations of the Culex pipiens complex
.
Amer. J. Trop. Med. Hyg.
,
9
,
321
330
.
Turell
,
M. J.
,
M. L.
O’Guinn
,
D. J.
Dohm
, and
J. W.
Jones
,
2001
:
Vector competence of North American mosquitoes (Diptera: Culicidae) for West Nile virus
.
J. Med. Entomol.
,
38
,
130
134
, doi:.
Turell
,
M. J.
,
M. L.
O’Guinn
,
D. J.
Dohm
,
J. P.
Webb
Jr
., and
M. R.
Sardelis
,
2002
:
Vector competence of Culex tarsalis from Orange County, California, for West Nile virus
.
Vector-Borne Zoonotic Dis.
,
2
,
193
196
, doi:.
Vinogradova
,
E. B.
,
2000
: Chapter 2: Physiology and ecology of Culex p. pipiens. Culex Pipiens Pipiens Mosquitoes: Taxonomy, Distribution, Ecology, Physiology, Genetics, Applied Importance and Control, Pensoft, 46–115.
Wagner
,
T. L.
,
H.
Wu
,
P. J. H.
Sharpe
,
R. M.
Schoolfield
, and
R. N.
Coulson
,
1984
:
Modeling insect development rates: A literature review and application of a biophysical model
.
Ann. Entomol. Soc. Amer.
,
77
,
208
220
, doi:.
Willmott
,
C. J.
,
S. G.
Ackleson
,
R. E.
Davis
,
J. J.
Feddema
,
K. M.
Klink
,
D. R.
Legates
,
J.
O’Donnell
, and
C. M.
Rowe
,
1985
:
Statistics for the evaluation and comparison of models
.
J. Geophys. Res.
,
90
,
8995
9005
, doi:.
Winter
,
C. L.
, and
D.
Nychka
,
2010
:
Forecasting skill of model averages
.
Stochastic Environ. Res. Risk Assess.
,
24
,
633
638
, doi:.

Footnotes

1

At the time of this analysis, the IPCC Fifth Assessment Report scenarios were not yet incorporated into LARS-WG.