Temperature-dependent transmission

Data gathered from the literature [9,13–30] revealed that all mosquito traits relevant to transmission—biting rate, egg-to-adult survival and development rate, adult lifespan, and fecundity—respond strongly to temperature and peak between 23°C and 34°C for the two mosquito species (Ae. aegypti in Fig 1 and Ae. albopictus in Fig A in S1 Text). DENV extrinsic incubation and vector competence peak at 35°C [31–37] and 31–32°C [31,32,34,38], respectively, in both mosquitoes—temperatures at which mosquito survival is low, limiting transmission potential (Fig 1, Fig A in S1 Text). Appropriate thermal response data were not available for CHIKV and ZIKV extrinsic incubation and vector competence.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Thermal responses of Ae. aegypti and DENV traits that drive transmission (data sources listed in Table B in S2 Text).

Informative priors based on data from additional Aedes spp. and flavivirus studies helped to constrain uncertainty in the model fits (see Materials and Methods; Table C in S2 Text). Points and error bars indicate the data means and standard errors (for display only; models were fit from the raw data). Black solid lines are the mean model fits; red dashed lines are the 95% credible intervals. Thermal responses for Ae. albopictus are shown in Fig A in S1 Text.

https://doi.org/10.1371/journal.pntd.0005568.g001

We estimated the posterior distribution of R₀(T) and used it to calculate key temperature values that indicate suitability for transmission: the mean and 95% credible intervals (95% CI) on the critical thermal minimum, maximum, and optimum temperature for transmission by the two mosquito species. At constant temperatures, Ae. aegypti transmission peaked at 29.1°C (95% CI: 28.4–29.8°C), and declined to zero below 17.8°C (95% CI: 14.6–21.2°C) and above 34.6°C (95% CI: 34.1–35.6°C) (Fig 2). Ae. albopictus transmission peaked at 26.4°C (95% CI: 25.2–27.4°C) and declined to zero below 16.2°C (95% CI: 13.2–19.9°C) and above 31.6°C (95% CI: 29.4–33.7°C) (Fig 2). Overall, the thermal response curve for Ae. albopictus is shifted towards lower temperatures than Ae. aegypti, so Ae. albopictus transmission is better suited to cooler environments. For a more realistic scenario in which daily temperature ranged over 8°C, the transmission peak, minimum, and maximum were slightly lower for both Ae. aegypti (28.5°C, 13.5°C, 34.2°C, respectively) and Ae. albopictus (26.1°C, 11.9°C, and 28.3°C, respectively). The lower thermal maximum under fluctuating temperatures occurs because we incorporated empirically supported irreversible lethal effects of temperatures that exceed thermal maxima for survival (see Materials and Methods).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Relative R₀ across constant temperatures (°C; top) for Ae. albopictus (light blue) and Ae. aegypti (dark blue), and histograms of the posterior distributions of the critical thermal minimum (bottom left), temperature at peak transmission (bottom middle), and critical thermal maximum (bottom right; all in °C).

Solid lines: mean posterior estimates; dashed lines: 95% credible intervals. R₀ curves normalized to a 0–1 scale for ease of comparison and visualization.

https://doi.org/10.1371/journal.pntd.0005568.g002

The posterior distribution of R₀(T) allows us to evaluate uncertainty in key temperature values that define the transmission range, including critical thermal minimum, maximum, and optimum. Uncertainty was higher for the critical thermal minimum for transmission than for the maximum or optimum, and the two mosquito species overlapped most for this outcome (Fig 2, bottom panels). This result occurred because several trait thermal responses increase gradually from low to mid temperatures but decline more steeply at high temperatures (Fig 1), so uncertainty is greatest at low temperatures. Ae. aegypti has a substantially higher optimum and maximum temperature than Ae. albopictus (Fig 2) due to its greater rates of adult survival at high temperatures (see Supplementary Materials for sensitivity analyses).

Model validation

We used generalized linear models (GLM) to ask whether the predicted relationship between temperature and transmission, R₀(T), was consistent with observed human cases of DENV, CHIKV, and ZIKV. Specifically, we assessed whether R₀(T) was an important predictor of the probability of autochthonous transmission occurring and of the incidence given that transmission occurred. We also controlled for human population size, virus species, and two socioeconomic factors. (Note that we focused on testing the R₀(T) model, rather than on constructing the best possible statistical model of human case data.) To do this, we used the version of the Ae. aegypti R₀(T) model that includes 8°C daily temperature range, along with country-scale weekly case reports of DENV, CHIKV, and ZIKV in the Americas and the Caribbean between 2014–2016. We first addressed the fact that countries with larger populations have greater opportunities for (large) epidemics by creating two predictors that incorporate scaled R₀(T) and population size. In the models of the probability of autochthonous transmission occurring we used the product of the posterior probability that R₀(T) > 0 (which we notate as GR₀) and the log of population size (p) to give log(p)*GR₀. (Here, and throughout, log denotes the natural logarithm.) In the models of incidence, given that transmission does occur, we used the log of the product of the posterior mean of R₀(T) and population size, log(p*R₀(T)). To control for several socioeconomic factors that might obscure the impact of temperature, we also included log of gross domestic product (GDP) and log of percent of GDP in tourism (using logs because the predictors were highly skewed, to stabilize variance). These are potential indicators of investment in and/or success of vector control and infrastructure improvements that prevent transmission. By comparing models that included the R₀(T) metric alone, socioeconomic factors alone, or both, we tested whether R₀(T) was an important predictor of observed transmission occurrence and incidence (see Table D in S2 Text). Note that R₀(T) is out of sample for all validation analyses because it is derived and calculated strictly from laboratory data on mosquitoes, and we perform a validation analyses for R₀(T) using independent case incidence reports. For this validation step we assessed model adequacy for the transmission data in two ways. First we used the full dataset for case incidence reports to select the best model (Table D in S2 Text) and to determine whether or not our predicted value of relative R₀(T) based on laboratory data was included in the model (“within sample” analysis). Second we used a bootstrapping approach where models were fit on subsets of the case incidence data that were randomly sampled and then predictive accuracy of the competing models (Table D in S2 Text) was assessed on left-out data (“out of sample” analysis).

For the probability of autochthonous transmission occurring, the model that included both the R₀(T) predictor and socioeconomic predictors had overwhelming support based on Bayesian Information Criterion (BIC; model PA5 relative probability = 1, Table D in S2 Text). Based on deviance explained, the models that included R₀(T), with or without the socioeconomic predictors out-performed the model that did not include R₀(T) (Table D in S2 Text; Fig 3A, Fig B in S1 Text). In analyses of out-of-sample accuracy, models that included the R₀(T) metric (with or without the socioeconomic factors) were surprisingly accurate. They predicted the probability of transmission with 86–91% out-of-sample accuracy for DENV (Table D in S2 Text). For CHIKV and ZIKV, models that included the R₀(T) metric or population alone had 66–69% out-of-sample accuracy (Table D in S2 Text). There were no significant differences in out-of-sample accuracy between the top four models, but for both DENV and CHIKV/ZIKV the best model was significantly better than the worst model [see supplementary code in 39 for full results]. The lower out-of-sample accuracy for CHIKV and ZIKV likely reflects the much lower frequency of positive values and the lower total sample size of this dataset. All results were similar for a set of models that separated GR₀ from population size, so for simplicity we show the model predictors that combines GR₀ and population size here (see Table D in S2 Text and [39] for results of other models). Further, from a biological perspective, the combined model better describes what we know about disease systems: if either the probability of R₀(T) being greater than zero is small or population size is very small, transmission is unlikely to occur. Together, these analyses suggest that R₀(T) is an important predictor of transmission occurrence, but that CHIKV and ZIKV need further data to better explain the probability of transmission occurrence (Fig 3A, Fig B in S1 Text).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Ae. aegypti R₀(T) and population size predict the probability and magnitude of transmission of DENV, CHIKV, and ZIKV across the Americas.

A, log(p)*GR₀ (the posterior probability that R₀(T) > 0 times the log of population size) versus the probability of local transmission in the data. B, log(p*R₀(T)) (log of R₀(T) times the population size) versus the log of incidence, given that it exceeds the threshold for local transmission. Tick-marks and points: human transmission occurrence and incidence data, respectively, by country-week in the Americas and Caribbean. Lines and shaded areas: mean and 95% CI from GLM fits for DENV (blue) and CHIKV and ZIKV (red). For simplicity, we show the models that only include the covariates log(p)*GR₀ or log(p*R₀(T)), respectively, and do not include the socioeconomic covariates (models PA6 and IM4 in Table D in S2 Text). For each case report data point, log(p)*GR₀ and log(p*R₀(T)) were calculated at the mean temperature 10 weeks prior to the reporting week [40].

https://doi.org/10.1371/journal.pntd.0005568.g003

R₀(T) was also an important predictor of incidence, given that autochthonous transmission did occur. Within-sample, incidence was best predicted by the model that included both R₀(T) and the socioeconomic predictors (model IM5 in Table D in S2 Text) based on BIC (relative probability = 1). The models that included R₀(T) out-performed those that did not based on deviance explained (Table D in S2 Text). In out-of-sample validation, the models that included R₀(T) explained the magnitude of incidence based on mean absolute percentage error (85–86% accuracy versus 83% accuracy for models that did not include R₀(T); Table D in S2 Text), but this difference was not statistically significant. For illustration, we show the simpler model that only contains the R₀(T) predictor in the main text (Fig 3B; model IM1 in Table D in S2 Text). Notably, the models that contained R₀(T) predicted incidence well for all three viruses, despite the lower incidence of CHIKV and ZIKV.

Although predicted R₀(T) correlated with the observed occurrence and magnitude of human incidence for all three viruses, these observed incidence metrics were higher for DENV than for CHIKV and ZIKV. While the reason for this difference is unclear, the most likely explanation is that DENV is much more established in the Americas, so it is more likely to be detected, diagnosed, and reported. Because ZIKV and CHIKV are newly emerging, they may not have fully saturated the region at this early stage.

The ability of the model to explain the probability and magnitude of transmission is notable given the coarse scale of the human incidence versus mean temperature data (i.e., country-scale means), the lack of CHIKV- and ZIKV-specific trait thermal response data to inform the model, the nonlinear relationship between transmission and incidence, and all the well-documented factors other than temperature that influence transmission. Together, these analyses show simple mechanistic models parameterized with laboratory data on mosquitoes and dengue virus are consistent with observed temperature suitability for transmission. Moreover, the similar responses of human incidence of ZIKV, CHIKV, and DENV to temperature suggest that the thermal ecology of their shared mosquito vectors is a key determinant of outbreak location, timing, and intensity.

Mapping climate suitability for transmission

The validated model can be used to predict where transmission is not excluded (posterior probability that R₀(T) > 0, a conservative estimate of transmission risk). Considering the number of months per year at which mean temperatures do not prevent transmission, large areas of tropical and subtropical regions, including Puerto Rico and parts of Florida and Texas, are currently suitable year-round or seasonally (Fig 4). These regions are fundamentally at risk for DENV, CHIKV, ZIKV, and other Aedes arbovirus transmission during a substantial part of the year (Fig 4). Indeed, DENV, CHIKV, and/or ZIKV local transmission has occurred in Texas, Florida, Hawaii, and Puerto Rico (www.cdc.gov). On the other hand, many temperate regions experience temperatures suitable for transmission three months or less per year (Fig 4). Temperature thus limits the potential for the viruses to generate extensive epidemics in temperate areas even where the vectors are present. Moreover, many temperate regions with seasonally suitable temperatures currently lack Ae. aegypti and Ae. albopictus mosquitoes, making vector transmission impossible (Fig 4, black line). The posterior distribution of R₀(T) also allows us to map months of risk with different degrees of uncertainty (e.g., 97.5%, 50%, and 2.5% posterior probability that that R₀ > 0), ranging from the most to least conservative (Fig D in S1 Text).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Map of predicted temperature suitability for virus transmission by Ae. albopictus and Ae. aegypti.

Color indicates the consecutive months in which temperature is permissive for transmission (predicted R₀ > 0) for Aedes spp. transmission based on the minimum likely range (> 97.5% posterior probability that R₀ > 0). Black lines indicate the CDC estimated range for the two Aedes spp. in the United States. Model suitability predictions combine temperature mean and 8°C daily variation and are informed by laboratory data (Fig 1, Fig A in S1 Text) and validated against field data (Fig 3).

https://doi.org/10.1371/journal.pntd.0005568.g004

Temperature-sensitive R₀ models

We constructed temperature-dependent models of transmission using a previously developed R₀ framework. We modeled transmission rate as the basic reproduction rate, R₀—the number of secondary infections that would originate from a single infected individual introduced to a fully susceptible population. In previous work on malaria, we adapted a commonly used expression for R₀ for vector transmission to include the temperature-sensitive traits that drive mosquito population density [12]: (1)

Here, (T) indicates that the trait is a function of temperature, T; a is the per-mosquito biting rate, b is the proportion of infectious bites that infect susceptible humans, c is the proportion of bites on infected humans that infect previously uninfected mosquitoes (i.e., b*c = vector competence), μ is the adult mosquito mortality rate (lifespan, lf = 1/μ), PDR is the parasite development rate (i.e., the inverse of the extrinsic incubation period, the time required between a mosquito biting an infected host and becoming infectious), EFD is the number of eggs produced per female mosquito per day, p_EA is the mosquito egg-to-adult survival probability, MDR is the mosquito immature development rate (i.e., the inverse of the egg-to-adult development time), N is the density of humans, and r is the human recovery rate. For each temperature-sensitive trait in each mosquito species, we fit either symmetric (Quadratic, -c(T–T₀)(T–T_m)) or asymmetric (Brière, cT(T–T₀)(T_m−T)^1/2) unimodal thermal response models to the available empirical data [48]. In both functions, T₀ and T_m are respectively the minimum and maximum temperature for transmission, and c is a positive rate constant.

We consider a normalized version of the R₀ equation such that it is rescaled to range from zero to one with the value of one occurring at the unimodal peak. Although absolute values of R₀ that are used to determine when transmission is stable depend on additional factors not captured in our model, the minimum and maximum temperatures for which R₀ > 0 map exactly onto our normalized equations, allowing us to accurately calculate whether or not transmission should be possible at all. Empirical estimates of absolute values of R₀ are difficult to obtain in any case, but it is much easier to determine whether transmission is occurring and for how long. While different model formulations for predicting R₀ versus temperature can produce results with different magnitudes and potentially different overall shapes [49], the temperatures for which R₀ is above or below zero (or one) are mostly model independent. For instance, two competing models differ only by whether or not the formula in Eq (1) is squared, but the square of a number (e.g., an absolute R₀ value) greater than one is always greater than one, and the square of a number less than one is always less than one. Therefore, the threshold temperatures at which absolute R₀ > 0 or absolute R₀ > 1 will be exactly the same for either choice of formula (Fig F in S1 Text). Similarly, because different expressions for R₀, including the square of Eq (1), map monotonically onto our function, they will produce identical estimates for the temperatures at which transmission declines to zero and peaks (Fig F in S1 Text). Consequently, our use of relative R₀ adequately describes the nonlinear relationship between mosquito and virus traits and transmission.

We fit the trait thermal responses in Eq (1) based on an exhaustive search of published laboratory studies that fulfilled the criterion of measuring a trait at three or more constant temperatures, ideally capturing both the rise and the fall of each unimodal curve (Tables S1-S2). Constant-temperature laboratory conditions are required to isolate the direct effect of temperature from confounding factors in the field and to provide a baseline for estimating the effects of temperature variation through rate summation [50]. We attempted to obtain raw data from each study, but if they were not available we collected data by hand from tables or digitized data from figures using WebPlotDigitizer [51]. We obtained raw data from Delatte [19] and Alto [21] for the Ae. albopictus egg-to-adult survival probability (pEA), mosquito development rate (MDR), gonotrophic cycle duration (GCD, which we assumed was equal to the inverse of the biting rate) and total fecundity (TFD) (Table D in S2 Text). Data did not meet the inclusion criterion for CHIKV or ZIKV vector competence (b, c) or extrinsic incubation period (EIP) in either Ae. albopictus or Ae. aegypti. Instead, we used DENV EIP and vector competence data, combined with sensitivity analyses.

Following Johnson et al. [52], we fit a thermal response for each trait using Bayesian models. We first fit Bayesian models for each trait thermal response using uninformative priors (T₀ ~ Uniform (0, 24), T_m ~ Uniform (25, 45), c ~ Gamma (1, 10) for Brière and c ~ Gamma (1, 1) for Quadratic fits) chosen to restrict each parameter to its biologically realistic range (i.e., T₀ < T_m and we assumed that temperatures below 0°C and above 45°C were lethal). Any negative values for all thermal response functions were truncated at zero, and thermal responses for probabilities (p_EA, b, and c) were also truncated at one. We modeled the observed data as arising from a normal distribution with the mean predicted by the thermal response function calculated at the observed temperature, and the precision τ, (τ = 1/σ), distributed as τ ~ Gamma (0.0001, 00001). We fit the models using Markov Chain Monte Carlo (MCMC) sampling in JAGS, using the R [53] package rjags [54]. For each thermal response, we ran five MCMC chains with a 5000-step burn-in and saved the subsequent 5000 steps. We thinned the posterior samples by saving every fifth sample and used the samples to calculate R₀ from 15–40°C, producing a posterior distribution of R₀ versus temperature. We summarized the relationship between temperature and each trait or overall R₀ by calculating the mean and 95% highest posterior density interval (HPD interval; a type of credible interval that includes the smallest continuous range containing 95% of the probability, as implemented in the coda package [55]) for each curve across temperatures.

We fit a second set of models for each mosquito species that used informative priors to reduce uncertainty in R₀ versus temperature and in the trait thermal responses. In these models, we used Gamma-distributed priors for each parameter T₀, T_m, c, and τ fit from an additional ‘prior’ dataset of Aedes spp. trait data that did not meet the inclusion criteria for the primary dataset (Table C in S2 Text). We found that these initial informative priors could have an overly strong influence on the posteriors, in some cases drawing the posterior distributions well away from the primary dataset, which was better controlled and met the inclusion criteria. We accounted for our lower confidence in this data set by increasing the variance in the informative priors, by multiplying all hyperparameters (i.e., the parameters of the Gamma distributions of priors for T₀, T_m, and c) by a constant k to produce a distribution with the same mean but 1/k times larger variance. We chose the value of k based on our relative confidence in the prior versus main data. Thus we chose k = 0.5 for b, c, and PDR and k = 0.01 for lf. This is the main model presented in the text (Fig 2). It is comparable to some but not all previous mechanistic models for Ae. aegypti and Ae. albopictus transmission (Fig E in S1 Text). Results of our main model, fit with informative priors, did not vary substantially from the model fit with uninformative priors (Figs G-H in S1 Text).

Incorporating daily temperature variation in transmission models

Because organisms do not typically experience constant temperature environments in nature, we incorporated the effects of temperature variation on transmission by calculating a daily average R₀ assuming a daily temperature range of 8°C, across the range of mean temperatures. This range is consistent with daily temperature variation in tropical and subtropical environments but lower than in most temperate environments. At each mean temperature, we used a Parton-Logan model to generate hourly temperatures and calculate each temperature-sensitive trait on an hourly basis [56]. We assumed an irreversible high-temperature threshold above which mosquitoes die and transmission is impossible [57,58]. We set this threshold based on hourly temperatures exceeding the critical thermal maximum (T_m in Tables A-B in S1 Text) for egg-to-adult survival or adult longevity by any amount for five hours or by 3°C for one hour. We averaged each trait over 24 hours to obtain a daily average trait value, which we used to calculate relative R₀ across a range of mean temperatures. We used this model in the validation against human cases (Fig 3) and the risk map (Fig 4).

Model validation with DENV, CHIKV, and ZIKV incidence data

To validate the model, we used data on human cases of DENV, CHIKV, and ZIKV at the country scale and mean temperature during the transmission window. Using statistical models (as described below), we estimated the effects of predicted R₀(T) on the probability of local transmission and the magnitude of incidence, controlling for population size and several socioeconomic factors. We downloaded and manually entered Pan American Health Organization (PAHO) weekly case reports for DENV and CHIKV for all countries in the Americas (North, Central, and South America and the Caribbean Islands) from week 1 of 2014 to week 8 of 2015 for CHIKV and from week 52 of 2013 to week 47 of 2015 for DENV (www.paho.org). ZIKV weekly case reports for reporting districts (e.g., provinces) within Colombia, Mexico, El Salvador, and the US Virgin Islands were available from the CDC Epidemic Prediction Initiative (https://github.com/cdcepi/) from November 28, 2015 to April 2, 2016. We aggregated the ZIKV data into country-level weekly case reports to match the spatial resolution of the DENV, CHIKV, and covariate data.

Temperature data collection

We matched the DENV, CHIKV, and ZIKV incidence data with temperature using daily temperature data from METAR stations in each country, averaged at the country level by epidemic week. A previous study found a six-week lagged relationship between temperature and oviposition for Aedes aegypti in Ecuador [40]. Assuming that the subsequent transmission, disease development, medical care-seeking, and case reporting in humans takes an additional four weeks, we assumed a priori a ten-week lag between temperature and incidence (i.e., mean temperature for the week that is ten weeks prior to each case report). METAR stations are internationally standardized weather reporting stations that report hourly temperature and precipitation measures. Outlier weather stations were excluded if they reported a daily maximum temperature below 5°C or a daily minimum temperature above 40°C during the study period, extremes that would certainly eliminate the potential for transmission in a local area. Because case data are reported at the country level, we needed a collection of weather stations in each country that accurately represent weather conditions in the areas where transmission occurs, excluding extreme areas where transmission is unlikely. For the study period of October 1, 2013 through April 30, 2016, we downloaded daily temperature data for each station from Weather Underground using the weatherData package in R [59]. We removed all data from Chile because it spans so much latitude and the terrain is so diverse that its country-level mean is unlikely to be very representative of the temperature where an outbreak occurred.

Socioeconomic covariate data

We accessed available data on projected 2016 gross domestic product (GDP) for countries of interest via the International Monetary Fund’s World Economic Outlook Database (http://www.imf.org/external/ns/cs.aspx?id=28). The direct and total contributions of tourism to GDP in 2016 were compiled from World Travel and Tourism Council economic impact reports (http://www.wttc.org/research/economic-research/economic-impact-analysis/country-reports/#undefined). We retrieved population size data for 2013–2015 from the United Nations Population Division (https://esa.un.org/unpd/wpp/Download/Standard/Population/) and averaged them across the three years for each country. Throughout the analyses below, unless otherwise specified, we used the natural log of the population size and of GDP as our predictors. We have two reasons for this choice. The first is that, intuitively, the relative order of magnitude of the population/GDP is more important in determining observed outbreak sizes or probabilities than their absolute sizes. Second, population sizes and GDPs across countries tend to exhibit clumped patterns with a few outliers that are much larger than the others. From a statistical perspective, using the un-transformed populations (or GDPs) results in those few large/rich countries having very high leverage in the analysis, and thus potentially skewing the results. Taking a log of the population better balances these predictors and is the standard accepted approach when using these kinds of predictors in regression models.

Validation analyses with human incidence versus temperature datasets

To validate the R₀(T) model while controlling for population and socio-economic factors, we used generalized linear models (GLMs) on the weekly case count data. Importantly, we focused on testing whether the case counts were consistent with the transmission–temperature relationship predicted from our model, rather than on maximizing the variation explained in the statistical model. We are more specifically interested in understanding autochthonous transmission (i.e., locally acquired, not just imported cases). We set country-level thresholds for the number of cases defining autochthonous transmission for our three diseases separately, based on current transmission understanding: seven cases of CHIKV, 70 cases of DENV, and three cases of ZIKV. We derived these thresholds in the following way. First, we looked for data on outbreaks of travel related cases in countries that are not expected to experience any local transmission. For instance, in 2014 Canada experienced 320 confirmed, travel-related cases of chikungunya (http://www.phac-aspc.gc.ca/publicat/ccdr-rmtc/15vol41/dr-rm41-01/rapid-eng.php), equivalent to an average of more than six cases per week. Thus, to be conservative in our estimates, we set the threshold of transmission as seven cases/week for CHIKV. The reported weekly cases of DENV transmission in our study sample are considerably higher than for CHIKV (mean DENV incidence was nearly 100 times higher mean CHIKV incidence). We chose a moderately high threshold of 70 cases in a week (i.e., 10 times higher than the CHIKV threshold based on Canadian cases) to reflect higher overall incidence and increased potential for travel related cases. We examined the sensitivity of the results to choice of threshold by varying it from 25 to 100, and we found qualitatively similar results for all thresholds that we tested. As ZIKV is not as well established as either CHIKV or DENV at this time, smaller numbers of cases may indicate autochthonous transmission. Consequently, we chose a threshold of three cases for ZIKV (approximately half the CHIKV threshold). Further, the results were fairly sensitive to the ZIKV threshold as many locations have small numbers of cases. Since higher thresholds exclude a very large proportion of available case data making analysis impossible, we used the slightly less conservative threshold of three cases for autochthonous transmission of ZIKV. The resulting data consisted of zeros for no transmission and positive case counts when transmission is presumed to be occurring. To model these data, we used a hurdle model that first uses logistic regression on the presence/absence of local transmission data to understand the factors correlated with local transmission occurring or not (PA analysis). Then we modeled the log of incidence (number of new cases per reporting week) for positive values with a gamma generalized linear models (incidence analysis).

We were interested in understanding whether R₀(T) was an important predictor of human transmission occurrence and incidence, after controlling for potentially confounding factors like population size and socioeconomic conditions. To do this, we fit a series of models with different subsets of predictors that included R₀(T) and population size, the socioeconomic variables, or both (see Table D in S2 Text for full models). To control for human population size, we created new metrics based on R₀(T) and population size to use for validation against the PAHO incidence data. We define GR₀, which is the posterior probability that R₀(T) > 0. We use log(p)*GR₀, where p is the population size, as the relevant R₀-based predictor for the PA analysis. For the incidence analysis, we instead use log(p*R₀(T)) as the predictor. In all cases log refers to the natural logarithm. For simplicity, we refer to these as the R₀(T) metrics hereafter and in the Results.

In both the PA and incidence analyses, we first used the full data sets to examine which of the candidate models best described the data. Randomized quantile residuals indicated that the logistic and gamma GLM models were performing adequately. We compared the approximate model probabilities, calculated from the BIC scores, as well as the proportion of deviance explained (D²) from each model. Next we examined the performance of the models in predicting out of sample, for both PA and incidence analyses. To do this we created 1000 random partitions, where 90% of the data were used to train the model and 10% were used for testing. In the PA analyses we classified each partition based on presence/absence, with separate classification thresholds for DENV versus CHIKV/ZIKV as these grouping had much different probabilities of occurrence. We assessed the performance of the model for the PA analysis based on the mean misclassification rate. In the incidence analyses we assessed the model performance based on the predictive mean absolute percentage error (MAPE). Since differences in prediction success between the models in both the PA and incidence analyses were not statistically significant, we present the simpler models that only include the R₀(T) metrics in the main text (Fig 3) and the models that additionally include socioeconomic covariates in the Supplementary Information (Figs B-C in S1 Text). We plotted the model predictions as a function of the R₀(T) metrics together with the observed data for the PA and incidence analyses using the R package visreg [60].

The residuals of the incidence model exhibit “inverse trumpeting,” in which residual variation is larger at low than high predicted incidence (Fig I in S1 Text). This occurs in part because we forced the model to go through the origin, i.e., no transmission when R₀(T) or the population size is equal to zero. However, the data did sometimes show transmission where we did not expect it, potentially because of imported cases, errors in reporting, or small pockets of transmission suitability in countries or times that are otherwise unsuitable on average. More local-scale case reporting that separates autochthonous from travel-associated cases would be needed to tease apart the source of this error.

Mapping temperature suitability for transmission

Using the validated model, we were interested in where the temperature was suitable for Ae. aegypti and/or Ae. albopictus transmission for some or all of the year to predict the potential geographic range of outbreaks in the Americas. We visualized the minimum, median, and maximum extent of transmission based on probability of occurrence thresholds from the R₀ models for both mosquitoes. We calculated the number of consecutive months in which the posterior probability of R₀ > 0 exceeds a threshold of 0.025, 0.5, or 0.975 for both mosquito species, representing the maximum, median, and minimum likely ranges, respectively. The minimum range is shown in Fig 4 and all three ranges are overlaid in Fig D in S1 Text. This analysis indicates the predicted seasonality of temperature suitability for transmission geographically, but does not indicate its magnitude. To generate the maps, we cropped monthly mean temperature rasters from 1950–2000 for all twelve months (Worldclim; www.worldclim.org/) to the Americas (R, raster package, crop function) and assigned cells values of one or zero depending on whether the probability that R₀ > 0 exceeded the threshold at the temperatures in those cells. We then synthesized the monthly grids into a single raster that reflected the maximum number of consecutive months where cell values equaled one. The resulting rasters were plotted in ArcGIS 10.3, overlaying the three cutoffs (Fig D in S1 Text). We employed this process for both mosquito species.