Effects of changes in temperature on Zika dynamics and control.

When a rare pathogen emerges to cause a pandemic, it is critical to understand its dynamics and the impact of mitigation measures. We use experimental data to parametrize a temperature-dependent model of Zika virus (ZIKV) transmission dynamics and analyse the effects of temperature variability and control-related parameters on the basic reproduction number (R0) and the final epidemic size of ZIKV. Sensitivity analyses show that these two metrics are largely driven by different parameters, with the exception of temperature, which is the dominant driver of epidemic dynamics in the models. Our R0 estimate has a single optimum temperature (≈30°C), comparable to other published results (≈29°C). However, the final epidemic size is maximized across a wider temperature range, from 24 to 36°C. The models indicate that ZIKV is highly sensitive to seasonal temperature variation. For example, although the model predicts that ZIKV transmission cannot occur at a constant temperature below 23°C (≈ average annual temperature of Rio de Janeiro, Brazil), the model predicts substantial epidemics for areas with a mean temperature of 20°C if there is seasonal variation of 10°C (≈ average annual temperature of Tampa, Florida). This suggests that the geographical range of ZIKV is wider than indicated from static R0 models, underscoring the importance of climate dynamics and variation in the context of broader climate change on emerging infectious diseases.

Introduction

Vector-borne viruses (arboviruses) are emerging threats to both human and animal health. The global expansions of dengue virus (DENV), West Nile virus (WNV), chikungunya (CHIKV) and most recently Zika virus (ZIKV) are prominent examples of how quickly mosquito-transmitted viruses can emerge and spread through naive host populations. Currently, 3.9 billion people living within 120 countries are at risk of mosquito-borne arboviral diseases [1] with effects on human well-being that can be devastating (e.g. death, illness, as well as social and human ramifications of Zika-induced microcephaly and other congenital disorders) [2]. Anticipating and preventing outbreaks of emerging mosquito-borne viruses across these host populations is a major challenge.

Despite growing research to develop new therapeutics and vaccines, mitigating arbovirus disease spread still depends on conventional mosquito control methods, often with mixed success. Developing tools that allow us to successfully predict outbreaks of these viruses and efficiently target current and future interventions to specific times and locations can aid effective mosquito and disease control. Such efforts are often limited by gaps in knowledge on the relationships among mosquito vectors, pathogens and the environment, especially for emerging arboviruses such as CHIKV and ZIKV. Even in well-researched disease systems (e.g. malaria and DENV), key transmission parameters are only estimated from a few studies [3-5].

Variation in environmental temperature has a strong impact on the environmental suitability for transmission risk across a diversity of vector-borne disease systems [6-9]. Mosquitoes are small ectothermic organisms, and their fitness [10,11], life history [12-17] and vectorial capacity [3-5,16,18-21] exhibit nonlinear, unimodal relationships with environmental temperature. Recent work by Tesla et al. [19] demonstrates such temperature–transmission relationships for ZIKV, a recently emerging pathogen. These temperature–transmission relationships have significant ramifications on how disease transmission varies seasonally, across geographical locations, and with future climate and land use change. Control tools being considered for use within integrated vector management strategies may also be affected by temperature, such as conventional chemical insecticides that target a diverse range of insect pests [22-27], including mosquitoes [28,29]. Furthermore, there is evidence that temperature could modify the efficacy of novel control interventions, such as mosquito lines transinfected with the intracellular bacteria Wolbachia [30-33].

Several modelling frameworks have been used to predict environmental suitability for vector-borne disease transmission, including, most recently, temperature-dependent R0 models [3-5,18,19] and compartmental models of vector-borne disease dynamics [8,34,35]. The parameter R0 is broadly considered to be the most important summary statistic in epidemiology and disease ecology. It is defined as the expected number of new human (respectively, mosquito) cases generated by a single infectious human (respectively, mosquito) introduced into a fully susceptible human (respectively, mosquito) population throughout the period within which that human (respectively, mosquito) is infectious [36]. As a simple metric, it can easily incorporate the nonlinear influence of multiple temperature-dependent mosquito and pathogen traits, and has been applied to define the thermal optimum and limits for malaria [4,5,37], DENV, CHIKV [3,38,39], ZIKV [3,19] and Ross River virus [18]. However, temperature-dependent R0 formulations only define the relative risk of disease emergence and do not predict the final epidemic size (or incidence). The derivation, interpretation, and validation of R0 models are thus problematic in highly variable systems [40]. Dynamical models of transmission that track densities of infectious individuals over time, on the other hand, can more readily capture the impact of varying environmental conditions.

To better understand potential climate effects on control strategies for ZIKV, we developed a temperature-dependent dynamical model based on recent experimental work characterizing temperature–trait relationships between ZIKV vector competence, extrinsic incubation rate, and the per capita daily mosquito mortality rate [41]. Unlike other published results in the literature (e.g. R0 model in Tesla et al. [19]), we model ZIKV transmission dynamics between humans and vectors and the flow of humans and vectors between various classes explicitly through a compartmental SEIR model for the human population and SEI model for the vector population. Because the model is dynamic, we are also able to account for seasonal temperature variation. The model and analysis differ from Huber et al. [8] in that, through numerical and sensitivity analyses, we explicitly analyse the simultaneous effects of parameters that are influenced by control measures (including vaccination) on both the basic reproduction number R0 and the final epidemic size. Our analysis thus addresses the following questions: (1) How do the thermal optima and ranges for R0 compare to those for the human final epidemic size? (2) How does seasonal temperature variation affect the final epidemic size relative to a constant temperature environment? (3) Which parameters have the greatest impact on R0 and the final epidemic size that can inform control efforts? (4) Are different thermal environments more or less suitable for specific control strategies?

Our results show that R0 and the final epidemic size were largely driven by different sets of parameters, with the exception of temperature being the dominant driver of both. Furthermore, the human final epidemic size was maximized across a wider range of temperatures than what would have been predicted from the temperature-dependent R0 model. The human final epidemic size was highly sensitive to seasonal temperature variation, suggesting the potential invasion map of ZIKV may be wider than previously reported. Furthermore, the effectiveness of potential control strategies (e.g. vaccines, drug treatment, and insecticides, assessed through model parameters that are influenced by those strategies) is predicted to be sensitive to such differences in seasonal temperature variation.

Methods

We construct a temperature-dependent compartmentalized model of ZIKV dynamics with and without seasonal temperature change. Where possible, model parameters are estimated from the most recent laboratory experiments on temperature effects on the life cycle of the virus [3,19]. We first compare how temperature dependence affects R0, the human ‘final epidemic size’ (total number of infected individuals over the course of the epidemic), and key ecological characteristics of the system, such as extrinsic incubation period, the probability of transmission from the mosquito to the human, the probability of transmission from the human to the mosquito, and daily rates of mosquito and egg to adult survival. To understand how temperature change can influence the effects of control measures, we then analyse the combined effects of temperature and parameters that correspond to disease control measures on R0 and the final epidemic size. These ‘control parameters’ include vaccination, recovery, vector biting rates, vector-to-human transmission, vector carrying capacity, egg survival, and adult mosquito survival. Through a Latin hypercube sampling-based sensitivity analysis [42], we identify key parameters that most drive the epidemiological outcomes (R0 and the final epidemic size). All simulations were carried out using Matlab R2019b.

The basic dynamic model

The model apportions humans into four groups based on ZIKV infection status that changes over time, t: susceptible S (not infected), latent E (contracted the virus, but not yet infectious), infectious I (contracted the virus and can transmit it), and recovered R with lifelong immunity. The mosquito population is divided into similar classes, where the state variables have subscript v, but without an immune class since it is assumed that infectious mosquitoes do not clear the virus once it is in the salivary glands. The total human and mosquito populations are N = S + E + I + R and N = S + E + I.

The model assumes a constant human population during the epidemic. Susceptible humans acquire the virus at rate (force of infection) , while susceptible mosquitoes acquire the virus at rate λ(I, N) = bβI/N, where b is the number of human bites per mosquito per unit time, β is the probability that an infectious mosquito successfully transmits the virus while taking a blood meal from a susceptible human (i.e. the transmission rate), and β is the probability that an infectious human successfully transmits the virus to a biting, susceptible mosquito (i.e. the infection rate). The respective average residence times of infected humans and mosquitoes in the latent classes are 1/σ and 1/σ, while the respective rates at which humans and mosquitoes become infectious are σ and σ. Humans are infectious for approximately 1/γ days before recovering with permanent immunity (γ is the per capita human recovery rate), while infectious mosquitoes remain infectious until they die. Mosquito recruitment occurs at a per capita rate f(I) = α(1 − (N/κ)), where κ is the carrying capacity (maximum number of mosquitoes a breeding site can support). Furthermore, α = θνϕ/μ consists of θ, or the number of eggs a female mosquito produces per day; ν, the probability of surviving from egg to adult; and ϕ, the rate at which an egg develops into an adult mosquito. Mosquitoes die naturally at per capita rate μ, where 1/μ is the average lifespan of mosquitoes. See figure 1 for a schematic of the model and table 1 for descriptions and baseline values of the parameters. The dynamic model for the Zika virus is described by the equationsDots denote differentiation with respect to time, t (in days). The dynamics of the total human population and mosquito populations are described, respectively, by the equations

Figure 1

Compartmental model of Zika virus transmission. Compartments are divided into humans (blue), and vectors (red), representing disease status, with transitions between compartments (rates) in solid lines. The transmission of Zika virus from humans to vectors is denoted by dashed lines, and from vectors to humans by the double dashed (short and long dashed) line. Rates of demographic change (births and deaths) in the vector population are denoted by dotted lines.

Table 1

Parameter definitions, baseline values, and ranges of values for the system (2.1) without temperature dependence.

parameter	description	value	range	source
γh−1	human infectious period	5 days	4–7 days	[43]
σh−1	intrinsic incubation period	5.9 days	3–14 days	[44,45]
σv−1	extrinsic incubation period	10 days	8–12 days	[46,47]
βvh	probability of a mosquito infecting a human	0.33	0.1–0.75	[48,49]
βhv	probability of a human infecting a mosquito	0.33	0.1–0.75	[48,49]
μv−1	mosquito lifespan	14 days	7–30 days	[50]
κv	mosquito carrying capacity	2 × 104	(1 − 5) × 104	assumed

Without Zika virus, the mosquito population grows according to equation (2.2), or , where is the initial mosquito population and for α > μ is the positive equilibrium obtained by setting the right-hand side of the equation to zero. Observe that , and that when α > μ, the total mosquito population relaxes on the equilibrium population () in the long run. Therefore, the equilibrium point is stable when α > μ and vanishes when α < μ. The case for which α < μ results in a trivial mosquito equilibrium represents a situation in which the mosquito population becomes extinct. Since is constant. To illustrate the dynamics of the system, we set N = 1000.

The next generation operator approach [36,51,52] is used to compute the basic reproduction number of the model (2.1). This involves re-writing equations (2.1) as two sub-systems—a disease-free sub-system (consisting of the equations for the susceptible human and vector classes and the recovered human class) and a disease sub-system (consisting of equations for the exposed and infectious human and mosquito classes—i.e. the second, third, sixth and seventh equations of the model (2.1)). The disease system is expressed as the difference of two vectors—a vector of new infections F = ((bβI/N)S, 0, (bβI/N)S, 0) and a vector of transitions V = (σE, − σE + γI, (σ + μ)E, − σE + μI). The corresponding matrices of new infections and transitions (given by the Jacobians of the vectors F and V), as well as the inverse of the matrix arewhere , α > μ. The next generation matrix is The spectrum (set of eigenvalues) of the next generation matrix is

The basic reproduction number of the model (2.1) is the spectral radius, i.e. the largest eigenvalue of the matrix . Hence, in the presence of the Zika virus, the basic reproduction number of system (2.1) is

The main difference between this R0 calculation and that from the Ross–MacDonald model [53,54] is in the probability that the mosquito survives the latent period. The Zika virus can spread when R0 > 1 and can be contained when R0 < 1.

Variants of the model (2.1) have been used to assess the impact of various control measures including insecticide-treated bed nets on vector-borne diseases such as malaria [55-60]. For the purposes of exploring control strategies, we also consider a variant of this model that includes vaccination, where vaccinated susceptible humans are assumed to enter the immune class directly. For this special case, the first and fourth equations of (2.1) are replaced by and , respectively, where δ represents the per capita human vaccination rate.

Introducing temperature

The majority of the parameters associated with the mosquito vector (θ, ν, ϕ, b, μ), as well as ZIKV transmission (β) and replication (σ), are known to be influenced by environmental and climate conditions [3,19]. We investigate the effects of temperature variation on the dynamics of the mosquito population and ZIKV transmission over time. We follow the approach in [8] and model temperature-dependent parameters with the functional forms presented in equation (2.4). We rely on values and ranges of temperature-dependent parameters from recent laboratory-generated analyses for Zika virus [19] and Ae. aegypti life-history parameters (i.e. the biting rate of mosquitoes, the number of eggs a female mosquito lays per day, the probability of an egg surviving to an adult mosquito, and the rate at which an egg develops into an adult mosquito) [3] as specified in table 2. As in [8], the temperature dependent parameters are based on the quadratic and Briere [61] functional forms, c(T − T)(T − T0) and , respectively, where T, c, T0 and T are the temperature, rate (or scaling) constant, critical thermal minimum temperature and critical thermal maximum temperature, respectively:The subscripts θ, ν, ϕ, b, β, σ and ls on c, T0, and T indicate that the temperature-related parameter (c, T0 or T) is associated with the corresponding original model parameter (θ, ν, ϕ, b, β or σ,). These relationships between some model parameters, model outcomes, and temperature are illustrated in figure 2.

Table 2

Parameter values for temperature-dependent functional forms. The temperature in degrees Celsius is denoted by T, the rate scaling factor is denoted by c (main Column 2), the minimum critical thermal temperature is denoted by T0 (main Column 3), and the maximum critical thermal temperature is denoted by T (main Column 4). The subscripts denote the corresponding traits, e.g. , and are for the number of eggs laid by a female mosquito per day, θ. See equation (2.4) for brief descriptions of the temperature-dependent parameters.

	c		T0		Tm
trait	mean	range	mean	range	mean	range	source
θv(T)	8.56 × 10−3	[3.78, 14.1] × 10−3	14.58	[8, 08, 20.60]	34.61	[34, 35.77]	[3]
νv(T)	−5.99 × 10−3	[−6.82, −5.13] × 10−3	13.56	[12.56, 14.51]	38.29	[38.29, 39.02]	[3]
ϕv(T)	7.86 × 10−5	[5.75, 9.93] × 10−5	11.36	[7.19, 15.03]	39.17	[39.17, 39.54]	[3]
bv(T)	2.02 × 10−4	[1.2, 2.8] × 10−4	13.35	[5.84, 14.82]	40.08	[36.60, 40.51]	[3]
βhv(T)	−3.54 × 10−3	[−5.6, −1.8] × 10−3	22.72	[21.09, 24]	38.38	[36.46, 40.25]	[19]
σv(T)	1.74 × 10−4	[5.4, 30.4] × 10−5	18.27	[7.68, 24]	42.31	[39.26, 45]	[19]
1/μv(T)	−3.02 × 10−1	[−4.68, −1.34] × 10−1	11.25	[6.3, 15.06]	37.22	[34.79, 39.57]	[19]

Figure 2

Effect of temperature T on mosquito and pathogen parameters: (a) mosquito biting rate, b; (b) the pathogen transmission probability from humans to mosquitoes, β; (c) the inverse of the pathogen extrinsic incubation period, σ; (d) the average mosquito lifespan, 1/μ; (e) the vector competence, vc = ββ; (f) the number of eggs laid by a female mosquito per day, θ; (g) the egg-to-adult mosquito development rate, ϕ; (h) the egg-to-adult mosquito survival probability, ν; (i) the basic reproduction number, R0; (j and k) the total infectious human population, I in thousands, and total vector population, I in hundreds of thousands; and (l) the infectious human population. Time t in (l) is in hundreds of days. The number of infectious individuals rises with temperature up to an optimal temperature between 29°C and 32°C. As temperatures are increased beyond the optimal, the number of infectious individuals falls.

We further introduce seasonal variation in the system by modelling temperature through the sinusoidal functional form:where T is the mean annual temperature and T is the amplitude (divergence from mean temperature or mid-point between the lowest and highest annual temperatures, i.e. T = (T − T0)/2, where T and T0 are the respective average maximum and minimum temperatures across the year), and the period is 365 days. Here, we are assuming that at the onset of an epidemic outbreak, the phase is set to zero.

Control strategies

To analyse the relationship between temperature and Zika control, we identified the following parameters of the system that correspond to potential control strategies: vaccination (δ) decreases susceptibility and is directly incorporated into the models as described above; recovery rates (γ) can, for example, be increased through treatment with antiviral medication; vector biting rates (b) can be reduced through decreasing exposure to mosquitoes with personal protection or household improvements; vector-to-human transmission probability (β) can decrease with transmission-blocking Wolbachia; the vector carrying capacity (κ) can be reduced by eliminating vector breeding grounds near human habitats; egg–adult survival probability (ν) can be reduced through larvicides; and adult mosquito survival rate (μ) can be decreased through indoor spraying and the use of adulticides. We investigate the effects of the interaction between temperature and these parameters that are influenced by common disease control methods on R0 and the final epidemic size (total I).

Sensitivity analysis

Two types of sensitivity analyses—local and global—were used to explore the impact of temperature and selected parameters that are affected by disease control measures on the basic reproduction number (R0) and the human final epidemic size (total I). The local sensitivity analysis was conducted by varying only one parameter while holding all other parameters fixed, or varying both temperature and a control-sensitive parameter while holding the other parameters fixed. Each parameter that was varied was divided into 50, 100 and 250 equally spaced points within biologically feasible bounds. See figures 2–6 for results. As the human vaccination rate (δ) does not appear explicitly in the expression of R0, we could not assess how temperature modified the effect of this parameter on the basic reproduction number. However, we explore the impact of temperature on the human vaccination rate and associated implications for I, the human final epidemic size (figure 3).

Figure 6

The human final epidemic size is sensitive to the annual mean (oscillation) temperature, T, and several important control-related parameters: (a) the seasonal divergence of the annual temperature from the mean, T; (b) the human recovery rate, γ; (c) the human vaccination rate, δ; (d) the scaling factor of vector biting rate, ; (e) the probability of transmission from the mosquito to the human, β; (f ) the scaling factor of the probability of transmission from the human to the mosquito, ; (g) the scaling factor of the vector mortality rate, c; (h) the vector carrying capacity, κ; and (i) the egg survival probability scaling factor, . The annual mean (oscillation) temperature varies along the x-axis, the seasonal divergence of the annual temperature from the mean and the other parameters that are sensitive to control measures vary along the y-axes, and the colour scale indicates the total infectious humans. Apart from (a), where the temperature amplitude (T) is varying, the amplitude is set at 10°C for the other plots. Plot (a) shows that there can be large epidemics even when mean temperatures are low if the seasonal variation (the amplitude) is high enough, as would be found in subtropical and temperate regions.

Figure 3

Effect of temperature and parameters that are sensitive to control measures on the final epidemic size (total infectious human population) for the model without seasonal temperature variability: (a) vector biting rate, b; (b) human recovery rate, γ; (c) the probability of transmission from the mosquito to the human, β; (d) human vaccination rate, δ; (e) the probability of transmission from the human to the mosquito, β; (f) vector carrying capacity, κ; (g) vector mortality rate, μ; and (h) the vector egg–adult survival probability, ν. With the exception of the human vaccination and recovery rates, temperature does not substantially alter the effects of most parameters that are sensitive to control measures on the human final epidemic size (the colour bands have little gradient in the parameter space).

Global sensitivity analysis is presented in figure 7. The analysis is carried out using the Latin hypercube sampling and partial rank correlation coefficient (PRCC) technique [42]. The process involves identifying a biologically feasible mean, minimum and maximum value for each of the parameters (e.g. [3,19]) and subdividing the range of each parameter into 1000 equal sub-intervals, assuming a uniform distribution between the minimum and maximum values of each parameter. We then sample at random and without replacement from the parameter distributions to generate an m × n Latin hypercube sampling matrix, whose m rows (i.e. 1000 rows) consist of different values for each of the model parameters and the n columns (corresponding to the number of parameters in the system) consist of different values for the same parameter. Thus, each row of the Latin hypercube sampling matrix provides a parameter regime that is used for computing the basic reproduction number, solving the dynamic system, and computing the human final epidemic size. The parameters, basic reproduction number, and the human final epidemic size are then ranked with partial correlation coefficients estimated for each parameter along with corresponding p-values. PRCCs range from −1 to 1 and are used to examine the correlation between model parameters and model outputs (R0 and the final epidemic size). This method thus identifies parameters with the most significant influence on model outputs; it does not quantify the effect of a change in a parameter on the output.

Figure 7

Global sensitivity analysis indicates the sensitivity of the basic reproduction number, R0 (a and c), and the final epidemic size, total I (b, d and e), to all model parameters. Bars indicate partial rank correlation coefficients (PRCC), illustrating the contribution of parameters to variability or uncertainty in the model outputs (R0 and I). (a–b) No temperature dependence; (b–c) temperature dependence but no temperature variation; (e) temperature dependence and temperature variation. Without accounting for temperature, models show different sets of drivers for R0 (a) than burden (b). When temperature is included, it is the dominant contributor to model sensitivity for all models (c–e).

Mapping seasonal control

We mapped the R0 as a function of monthly mean temperature (figure 5). Globally gridded monthly mean current temperatures were downloaded from WorldClim.org [62], at a 5 arc-minute resolution (approx. 10 km2 at the equator), and predicted rates as a function of temperature at 0.2°C were mapped to the global grids. All raster calculations and graphics were conducted in R, using package raster [63].

Figure 5

Estimated basic reproduction number R0, for different mean monthly temperatures, globally. The colour bar represents increasing R0 values with the lowest values denoted by dark blue and the highest values denoted by light yellow. A value of R0 below 1 means the disease cannot take off. Note the month-to-month variation in locations such as Australia and the southern USA, demonstrating the impact of seasonal temperature variation on R0, and therefore on the control level needed to reduce or eradicate disease.

Results

Impact of temperature on model parameters and key outputs

The models show unimodal relationships between temperature and the temperature-dependent parameters, resulting in an optimal temperature that maximizes parameter values and a critical minimum and maximum temperature at which parameter values go to zero. Figure 2 presents the effects of temperature on mosquito and pathogen parameters, the final epidemic size in humans (total number of infected individuals over the course of the epidemic) and mosquitoes, and the basic reproduction number, R0 (via temperature effects on mosquito and ZIKV parameters).

The response of R0 to temperature is strongly peaked (with an optimum around 30°C) as has been demonstrated for Zika and other systems (e.g. dengue, malaria, Ross River virus) [3,5,18,19]. By contrast, the relationship between the final epidemic size and temperature is flat; i.e. within the thermal range of disease transmission, the final epidemic size does not change (figure 2i versus 2j). At temperatures associated with lower epidemic peaks, there are longer epidemic periods, resulting in the same total number of infected individuals over the course of the epidemic (figure 2l).

Zika virus control

The effects of static temperature (i.e. not including seasonality) and control-related parameters on final epidemic size are presented in figure 3. The most striking feature of these plots is that the difference between very small and very large epidemics (represented by blue and red areas, respectively) is discrete for most parameters. Crossing these thresholds, the most significant changes occur not from incremental changes in control parameters at a given temperature, but when temperatures move into the suitable thermal range, in which case the models predict widespread transmission (this concept is more rigorously explored in the sensitivity analysis below). However, vaccination (δ) and recovery (γ) have more incremental effects (within their feasible ranges) on final epidemic size even within the suitable thermal range; i.e. proportion of the population that needs to be vaccinated is higher at optimal temperatures than at sub-optimal temperatures to achieve a given reduction in the overall final epidemic size (figure 3d). Thus, warming temperatures (for most countries) will require greater vaccination coverage and treatment rates in order to maintain control of Zika. The effects on the basic reproduction number (R0) of most parameters that are sensitive to control measures were more dependent on temperature (figure 4) than was the case for the final epidemic size (total I), with the greatest effect involving the clearance rate of infection (γ), the probability of transmission from an infectious human to a susceptible mosquito (β), mosquito carrying capacity (κ) and the mosquito mortality rate (μ).

Figure 4

Effect on the basic reproduction number, R0, of temperature and parameters that are sensitive to control measures for the model without seasonal temperature variability: (a) the mosquito biting rate, b; (b) the human recovery rate, γ; (c) the probability of transmission from the mosquito to the human, β; (d) the probability of transmission from the human to the mosquito, β; (e) the mosquito carrying capacity, κ; and (f) the mosquito mortality rate, μ. (Vaccination δ and the egg-to-adult mosquito survival probability, ν, do not appear explicitly in the R0 model.) The effect of these parameters on R0 is more dependent on temperature than their effects on the final epidemic size (e.g. the colour bands are less vertical and more diagonal in parts of the parameter space).

Seasonal variation

Seasonal temperature variation affects outcomes by providing transient temperatures (variation from the annual mean; equation (2.5)) where the basic reproduction number can rise above (or fall below) 1, allowing for transmission to occur (or cease). At constant temperatures, epidemics only occur in humans between ≈23−37°C (figure 3). Such average annual temperatures are only found in tropical countries. However, the model shows the potential for epidemics in areas with mean temperatures below 23°C if there is adequate seasonal variation. This would be the case, for example, for a subtropical area, such as Tampa, FL, with a mean annual temperature of ≈20°C and an amplitude of ≈10°C (figure 6a). In contrast to the models without seasonal variation (figure 3), the models with seasonal variation (figure 6) indicate that the effectiveness of parameters on the final epidemic size is generally sensitive to changes in temperature (e.g. the colour bands in the subplots of figure 6 are diagonal in more of the parameter space than they are in figure 3). Figure 5 shows how the thermal conditions that are suitable for Zika (where R0 > 1) change with seasonal temperature variation across the globe.

Global uncertainty and sensitivity analysis

A global sensitivity analysis using Latin hypercube sampling showed that R0 and the final epidemic size are largely sensitive to different parameters (as indicated by differences in figures 3 and 4). However, temperature is a dominant driver of variation in both the basic reproduction number (R0) and the final epidemic size (total I) when it is included in the model (figure 7c–e). The human recovery rate, γ, was a consistently influential driver of the final epidemic size. By contrast, the basic reproduction number was not sensitive to recovery in the models with and without variable temperature. While R0 was also sensitive to vector competence (β and β), biting rate (b), and mosquito lifespan (1/μ), total infection burden was far less sensitive to these parameters and was mainly sensitive to human recovery rate (γ) (figure 7a,b).

Discussion

We are interested in what drives arbovirus epidemics, with Zika as a model, and how to reduce the burden of these diseases, focusing on temperature and key parameters that correspond to existing or potential control methods (e.g. pesticides, reduced breeding habitats, vaccines or treatment). We investigated temperature-dependent dynamic transmission models that incorporated recent empirical estimates of the relationships between temperature and Zika infection, transmission, and mosquito lifespan [19]. These dynamical models that can measure final epidemic size and account for temperature variation generate qualitatively different results from static R0 models. Temperature had an overwhelmingly strong impact on both R0 and the final epidemic size (total infectious individuals, equivalent to area under the I epidemic curve), but the response was much more gradual and had a clear optimum for R0, while the final epidemic size responded as a threshold function (figure 2i,j). This is because, while epidemics have a higher peak at the maximum R0 (at optimal temperatures), the epidemics are longer at sub-optimal temperatures (lower R0). Thus, Zika virus is capable of spreading efficiently through the host population (high I) across a broad range of temperatures for which R0 > 1, spanning from 17 to 37°C in constant environments (figure 6a) [8]. This is broadly consistent with the high seroprevalence of Zika found in a number of countries [64,65]. This suitable temperature region expanded and shifted toward cooler mean temperatures under seasonally varying environments (figure 6a).

These results have two key implications. First, large epidemics can occur under realistic, seasonally varying, temperature environments even in regions where the mean temperature alone would be expected to suppress transmission, for example in a location with a mean of 20°C and a seasonal amplitude of 10°C (e.g. Tampa, FL). Second, temperature determines both upper and lower thresholds for whether or not epidemics are possible [8]. However, within the predicted suitable temperature range defined by R0, the final epidemic size is largely limited by the density of susceptible hosts (figures 2 and 6a) [8]. More broadly, the results highlight the important principle that metrics of transmission (e.g. R0) have a nonlinear relationship with the human final epidemic size (total I) and contribute distinct implications for our understanding of the transmission process.

Whether or not temperature affects the potential for disease control is an important applied question for designing public health campaigns (either via vector control, reduction in host biting rate, vaccination, or drug administration). Temperature did not strongly affect the impact of most control-related parameters on the final epidemic size when the models did not include seasonal variation. In the more realistic situation where models included seasonal variation, the effectiveness of most parameters that are sensitive to disease control measures depended on temperature. In all models, human vaccination rate required to control epidemics varied strongly with mean temperature (figures 3–6). Achieving herd immunity and thereby suppressing transmission via vaccination is more difficult when temperatures are highly suitable (20−35°C under constant temperatures or 15−32°C under varying temperatures; figures 3–6). By contrast, the effects of the human recovery rate (γ) and the vector mortality rate (μ) on R0 were sensitive to temperature, but their effects on the final epidemic size were not sensitive to temperature.

Similar to previous work on dengue [8], our results show that Zika can invade and cause large outbreaks during the summer in seasonally varying environments with lower average temperatures, such as temperate regions of the USA, Europe, and Asia. This implies that differences in the size of epidemics in tropical versus temperate locations occur not just because of differences in temperature (and its impacts on R0) but also because of differences in vector breeding habitat availability, humidity, human mosquito exposure, and other socio-ecological factors. Much of the globe—including regions in temperate, subtropical, and tropical climates—is already suitable for Zika transmission for all or part of the year, and climate change is likely to expand this suitability geographically and seasonally [66]. However, processes that increase the density of susceptible human populations and their exposure to mosquitoes, including urbanization and urban poverty, human population growth, and the growth and geographical expansion of vector populations, are likely to expand the burden of Zika even more dramatically in the future.

There are several limitations of such a modelling study. First, the parameters are determined by a combination of laboratory-based estimates as well as from the literature on dengue, instead of being fitted to empirical data on the spatio-temporal dynamics of Zika from the field. Such divergent approaches can generate different parameter estimates. Furthermore, the projections of the model on potential geographical distribution of Zika are based on the average of temperatures by country and season with constant parameters. In reality, there is substantial heterogeneity of temperature and parameters over time and space, which have important implications for disease dynamics. For this reason, further investigation of Zika models that are calibrated from field based dynamics will be valuable for a fuller understanding of effects of temperature variation on Zika control.

Conclusion

The unexpected emergence and global expansion of Zika in 2015–2017 and its association with Zika congenital syndrome and Guillain–Barre syndrome revealed once again how poorly prepared the global community is for the looming and expanding threat of vector-borne diseases. Given the recent history of the rapid global expansion of Aedes aegypti-transmitted viruses (e.g. DENV, CHIKV and ZIKV), as well as the challenges of controlling these epidemics, understanding the ecological drivers of transmission and their effects on potential disease control tools is crucial for improving preparedness for future vector-borne disease emergence. If a Zika vaccine becomes available, then the precisely defined temperature thresholds for large epidemics predicted in our model imply that vaccination targets should be set based on climate. By contrast, because other potential interventions that would reduce vector population sizes, biting rates, and human recovery rates act more independently of temperature, targets could be set based on other socio-ecological factors in a given epidemic setting. This dynamic temperature-dependent modelling framework, which depends most strongly on vector and host parameters that are virus-independent, may be a useful first step for responding to future Aedes-borne disease epidemics.