Vector-borne diseases constitute a significant source of the global disease burden (WHO, 2019). Vectors and pathogens are ectotherms whose life history traits and contact rates are directly affected by temperature. There is increasing evidence that climate warming will not only increase global mean temperatures but also hot extremes (IPCC, 2018; Sherwood et al., 2020). This means that diseases are likely to go extinct from areas that become too hot and spread to cooler climates at higher latitudes. Knowing when and where this will happen is crucial for developing cost-effective public health policies to combat vector-borne diseases.

For well over a century, vector-borne disease spread has been predicted using the Basic Reproductive Number (R₀) (Ross, 1911; MacDonald, 1955, 1957; Aron, 1988; May and Anderson, 1991), the expected number of secondary cases produced by a single case during its infectious period (Diekmann et al., 1990; Dietz, 1993). This metric was initially derived using the classical Ross-Macdonald model for Malaria with time delays (Ross, 1911; MacDonald, 1955; Diekmann et al., 1990; Dietz, 1993), under the strict assumption of a constant environment and equilibrium conditions. Specifically, R₀ is the dominant eigenvalue of the next generation matrix evaluated at the disease-free equilibrium (Diekmann et al., 1990; Dietz, 1993). Despite this history, R₀ computed as the dominant eigenvalue is increasingly being used to predict the effects of climate warming on vector-borne disease spread [e.g., Mordecai et al. (2013, 2017), Johnson et al. (2015), Shocket et al. (2018, 2020), Tesla et al. (2018); Cator et al. (2020)]. The problem is that climate warming violates both the constant environment and equilibrium assumptions.

It is not uncommon in theoretical ecology to apply models to real-world situations that violate model assumptions. However, doing this in the case of disease spread is risky, given the direct impacts of model predictions on public health. Climate warming is a phenomenon that is not only non-constant but also non-stationary, meaning that, unlike typical seasonal variation, it does not settle into a predictable pattern of temperature variation. In such an environment, there is no disease-free equilibrium, and the eigenvalue approach no longer holds. The R₀ formula used to predict climate warming, retains the eigenvalue approach but assumes the vector population to be growing exponentially in the absence of the disease (Mordecai et al., 2013, 2017; Johnson et al., 2015; Shocket et al., 2018, 2020; Tesla et al., 2018; Cator et al., 2020). This itself is a contradiction, since the eigenvalue approach is predicated on equilibrium assumptions. The studies using the R₀ formula attempt to circumvent this problem by assuming the vector population is in climatic equilibrium under exponential growth, i.e., for each climate there is a population size at which the vector's intrinsic growth rate is zero (Parham and Michael, 2010). But, climate warming is a dynamical process, which means that a climatic equilibrium is never reached.

These complications suggest the need for alternative approaches that can accommodate the non-stationary thermal environment generated by climate warming. Here, we propose a trait-based framework that provides such an alternative. This work goes above and beyond previous work in three important ways. First, we use mechanistic response functions, derived from first principles of thermodynamics, to characterize vector and pathogen trait responses to temperature. Second, we incorporate vector trait response functions into stage-structured models that realistically depict the developmental delays characteristic of vector life cycles. Third, we couple the vector dynamics model with trait-based epidemiological models of host-vector interactions to predict disease prevalence under the most recent IPCC predictions of hotter-than-average summers and warmer-than-average winters (IPCC, 2018). Our mechanistic approach allows us to make predictions about disease prevalence based solely on trait response data and completely independently of population-level information, a potentially powerful approach for formulating public health policy in the age of warming.

We defer the mathematical details to Appendices and focus instead on the biological implications. We start with a summary of past work, with emphasis on the use of R₀ to predict disease spread. We then explain why R₀ is not suitable for predicting disease spread under warming. We present an alternative approach, based on a dynamical model, for investigating climate warming effects on vector-borne diseases.

The prevailing approach of describing trait responses to temperature is to fit a range of phenomenological functions (e.g., linear, quadratic) to the data, and to choose the best-fit function based on statistical criteria (e.g., AIC; Burnham and Anderson, 2002). We take the approach of deriving mechanistic response functions based on how temperature affects the underlying biochemical processes (Appendix B). This has the advantage that the functions are derived from first principles of thermodynamics, and, hence, independently of the data. We then fit the functions to the data and estimate parameters (Appendix C). An additional advantage of our approach is that both the response functions and their parameters have clear biological meaning (e.g., physiologically optimal temperature, thermal sensitivity, thermal tolerance limits). Because the response functions are mechanistic and derived independently of the data, a poor statistical fit can highlight data limitations that could lead to mis-specification of trait responses (e.g., characterizing a left-skewed function as exponential due to lack of data at high temperatures). Importantly, its mechanistic underpinnings make our approach less dependent on species-specific data. Several large-scale data analyses show that the qualitative nature (e.g., left-skewed, Gaussian) of trait responses to temperature is conserved across ectotherm taxa, and that their parameter values are thermodynamically constrained to fall within a narrow range (Gillooly et al., 2001, 2002; Brown et al., 2004; Savage et al., 2004; Dell et al., 2011). This holds the promise that our mechanistically-derived parameter estimates may yield model predictions that apply across different vectors and pathogens inhabiting different geographic areas and latitudes.

Figure 1 shows that these mechanistic functions provide a very good fit to temperature response data for the Malaria pathogen (Plasmodium) and its mosquito vectors (Anopheles species). Table 1 provides the resulting parameter estimates.

Incorporating trait response parameters for the Malaria vector and pathogen into our derivation of R₀ based on a self-regulated vector population [Equation (1)] and the previous derivation assuming an exponentially growing vector population (Equation (2)), we see that the latter is greater in magnitude, attains a maximum at a lower temperature, and has a lower high temperature threshold for disease spread compared to the former (Figure 2). This means that the R₀ formulation used in previous studies (Parham and Michael, 2010; Mordecai et al., 2013, 2017; Johnson et al., 2015; Shocket et al., 2018, 2020; Tesla et al., 2018) over-predicts disease spread in cooler environments and under-predicts it in warmer environments. This does not inspire confidence that the widely used formulation for R₀ accurately depicts disease spread, even in constant thermal environments. Given that R₀, even when correctly derived, cannot be used to predict disease spread under non-constant, non-equilibrium conditions (Aron, 1988; Diekmann et al., 1990; Dietz, 1993), we need to look for metrics for predicting disease spread under climate warming that do not depend on these restrictive assumptions.

The obvious alternative to R₀ to predicting disease spread is prevalence, the proportion of infected hosts in the population (Diekmann et al., 1990; May and Anderson, 1991; Dietz, 1993). Prevalence has the advantage that it can be calculated in both stationary and nonstationary environments. Quantifying prevalence, however, requires a dynamical model of host-vector interactions. The host-vector model we have developed (Appendix A, Equation A.1) goes above and beyond previous models in incorporating all relevant time delays (latent period in the host, juvenile development delay in the vector, and the pathogen's extrinsic incubation period) as well as vector self-regulation at different life stages. It has the advantage that it can predict disease prevalence under any variable thermal regime based solely on how temperature affects the vector's and pathogen's life history traits, and thus completely independently of population-level (e.g., time series) data. However, the biological realism of the model comes at the cost of analytically intractability, and the numerical methods used for solving delay differential equations are out of reach for most ecologists and epidemiologists.

Here we adopt an approach that combines the biological realism of delay differential equation (DDE) models with the tractability of ordinary differential equation (ODE) models. A number of previous studies have used class-structured epidemiological models (e.g., Susceptible-Exposed-Infected-Removed), constructed using ordinary differential equations, to model such diseases as Zika, Malaria, Chikungunya, West Nile, and Leishmaniasis [e.g., Wonham et al. (2004), Bacaër and Guernaoui (2006), Bacaër (2007), Hartemink et al. (2009, 2011), Alonso et al. (2010), Gao et al. (2016)]. None of these studies, however, considered the temperature dependence of vector and pathogen traits. We develop a stage- and class-structured ODE model that explicitly incorporates mechanistic descriptions of vector and pathogen trait responses to temperature. The model depicts time delays in terms of instantaneous rates, and is therefore more tractable both analytically and numerically. We use the model to investigate climate warming effects on disease prevalence using the latest IPCC predictions (IPCC, 2018). We defer the mathematical details (Appendix D) and focus instead on parameterizing the model with trait response data to investigate disease prevalence under climate warming.

We use the ODE model (Appendix D, Equation D.1) to predict warming effects on disease prevalence, using seasonal variation as the basis for comparison. We depict seasonal variation using the sinusoidal function T(t) = M_T+A_TS(t) where t is time in days, M_T is the mean habitat temperature in K, A_T is the amplitude of seasonal fluctuations (AT=Tmax-Tmin2), and S(t)=sin2πtyr (or -cos2πtτ) with yr = 365 days. We incorporate latitudinal variation into our analyses by using climate data to derive a statistical relationship between mean habitat temperature and the amplitude of seasonal fluctuations. Based on 47 different latitudes ranging from the tropics to the high temperate (and including many locations in west Africa where Malaria is endemic), we obtained the following relationship: A_T = 151.46−0.49M_T (see online Figure E1 in Appendix E).

Pathogens (e.g., virus, bacteria, protists) and disease vectors (e.g., mosquitoes, flies, ticks) are ectotherms whose life history traits (e.g., birth, maturation, mortality) and contacts with hosts (e.g., biting rate, transmission probability) depend directly on the environmental temperature. Climate warming is therefore likely to have strong effects on disease spread and prevalence. Predicting warming effects on disease burden requires that we understand how temperature effects on vector and pathogen traits translate into population-level patterns of vector abundance and dynamics, and how these in turn influence vector species' interactions with their (mostly) endotherm hosts. Making accurate predictions about disease spread and prevalence is a high research priority, given the pressing need to develop cost-effective public health policies in an increasingly warmer world.

The typical approach to predicting effects of climate warming on vector-borne diseases is to use the basic reproductive number (R₀) to predict disease spread (Mordecai et al., 2013, 2017; Johnson et al., 2015; Shocket et al., 2018, 2020; Tesla et al., 2018; Cator et al., 2020). This metric, however, is derived under assumptions of constant thermal environments and equilibrium conditions (MacDonald, 1955, 1957; Aron, 1988; Diekmann et al., 1990; Dietz, 1993; Parham and Michael, 2010). It is not applicable to the nonstationary (i.e., non-constant, non-equilibrium) thermal environment generated by climate warming. Using R₀ to predict disease spread under warming runs the risk of making erroneous predictions about warming that could misguide policy decisions. The challenge is to find an alternative metric that both realistically depicts temperature effects on vector and pathogen traits, and yields reliable metrics of disease spread and prevalence.

Here, we present a theoretical approach that combines a stage-structured model of vector population dynamics with a class-structured epidemiological model of vector-host dynamics. It has the advantage of being able to predict disease spread and prevalence under any type of non-equilibrial, non-stationary environment, including seasonal variation and warming. The novelty of our approach is that we incorporate mechanistic descriptions of vector and pathogen trait responses to temperature into stage-structured models that realistically capture the complex life cycles of arthropod disease vectors. We attempt to strike a balance between biological realism, mathematical rigor, and tractability of numerical analyses. We report three key findings.

First, we find that fitting mechanistic temperature response functions do not require more data than do the current phenomenological approaches. Use of mechanistic response functions is important because it provides a common framework for characterizing vector and pathogen trait responses to temperature that can apply broadly across vector taxa, pathogen types, and geographic location. The functions we use, which are derived from first principles of thermodynamics, depict how temperature effects on underlying biochemical processes (e.g., reaction kinetics, enzyme inactivation) translate into phenotypic-level traits such as the vector's birth, maturation and biting rates, and the pathogen's incubation period within the vector. Several large-scale data analyses (Dell et al., 2011; Englund et al., 2011) have shown that the qualitative nature of these trait responses (e.g., Gaussian, left-skewed; Figure 1) are conserved across ectotherm taxa, and that their parameters are thermodynamically constrained to take only a narrow range of values (Gillooly et al., 2001, 2002; Savage et al., 2004). This means that one can use parameter values from related species when data are unavailable for focal species, a practice already employed in previous studies that use phenomenological trait response functions (Mordecai et al., 2013, 2017; Johnson et al., 2015).

Second, by deriving the Basic Reproductive Number (R₀) from an underlying model of vector-host dynamics that incorporates all relevant time delays (latent period in the host, juvenile development delay in the vector, and the pathogen's extrinsic incubation period), we show that the R₀ formulation commonly used to predict warming effects on disease spread (Mordecai et al., 2013, 2017; Johnson et al., 2015; Shocket et al., 2018, 2020; Tesla et al., 2018; Cator et al., 2020) is unsuitable for making such predictions. First, it excludes disease-induced mortality in the host, despite strong empirical evidence to the contrary. Second, it assumes an exponentially growing vector population, which violates the fundamental assumption in deriving R₀ as the dominant eigenvalue of the next generation matrix evaluated at the disease-free equilibrium. Third, it assumes that the exponentially growing vector population is in climatic equilibrium with the environment, making it inapplicable to the non-equilibrium conditions generated by climate warming. A comparison of the previous formulation of R₀ with one we have derived for a self-regulated vector population shows that the former overestimates disease spread in cooler environments and underestimates it in warmer environments, proving to be unreliable even in a constant thermal environment.

Our third finding is on the climate warming effects on vector-borne disease prevalence. The most recent IPCC predictions (IPCC, 2018) suggest an increase in mean temperatures as well as hot extremes. We find that the hotter summers scenario, which incorporates such extremes substantially lowers the upper thermal limit for disease persistence and reduces disease prevalence within the viable region. Because the lower thermal limit, which is set by infected vector population viability, is relatively inflexible, warming has the overall consequence of narrowing the thermal range for disease viability.

A closer inspection of vector-host dynamics reveals an interesting outcome: uninfected vector populations are largely unaffected by warming, while infected vector and host populations are greatly reduced even at low levels of warming. This finding was confirmed by a separate analysis of Malaria prevalence in several west African locations. The reason for this is as follows. Data show that density-independent mortality increases with increasing temperature above a critical threshold, and increases with decreasing temperature below it [Savage et al. (2004) and references in Gillooly et al. (2001, 2002)] (Figure 1). Because such mortality affects all vector classes regardless of their infection status, its effects are compounded as the uninfected vector class moves through the exposed class to the infected class (see Appendix D, Equation D.1). As a result, the infected vector population will necessarily have a narrower thermal limit for viability, and exhibit lower abundances at low and high temperature extremes, compared to its uninfected counterpart.

While higher mortality of the infected vector population at lower temperatures explains why the disease cannot spread to higher latitudes despite warming, higher vector mortality at higher temperatures cannot, in itself, explain the warming-induced decline in the infected host population. We need to also understand how warming affects the transition (maturation) and transmission rates. It is well-known, based on data from a large number of multicellular ectotherms, that the maturation rate exhibits a left-skewed response to temperature, with a rapid decline at high temperatures (Johnson and Lewin, 1946; Sharpe and DeMichele, 1977; Schoolfield et al., 1981; Van der Have and de Jong, 1996; Van der Have, 2002; Ratkowsky et al., 2005; Dell et al., 2011). The sporozoite maturation rate of Protozoan parasites (e.g., Plasmodium, Leishmania) within their vectors exhibits such a left-skewed response. Since it is directly dependent on the sporozoite maturation rate, the per-bite transmission probability of the pathogen from vector to host also exhibits a left-skewed temperature response (Fig. 1). Moreover, the difficulty of obtaining biting rate data has led to approximating the biting rate as the inverse of the gonadotropic cycle (egg maturation time of mosquitoes; Alonso et al., 2010; Hartemink et al., 2011; Mordecai et al., 2013; Johnson et al., 2015). But, the mosquito egg maturation rate itself exhibits a left-skewed response (Fig. 1). The key point is that the transition from exposed to infected vectors, and disease transmission from infected vector to susceptible host are both driven by four temperature response functions, all of which are left-skewed and decline rapidly at high temperature extremes. Not only does the high-temperature decline in the pathogen's maturation rate further reduce the number of infected vector individuals at any given time, but the concomitant decline in the biting rate and per-bite transmission probability means that fewer host individuals will be infected. The key point is that it is the interactive effect of the temperature response of mortality with those of transition (maturation) and transmission traits that drives the differential effect of warming on the uninfected vector class, and on the infected vector and host classes.

Given the inapplicability R₀ to the nonstationary environment generated by warming, a more informative approach is to use a trait-based model of vector-host interactions to calculate disease prevalence under the appropriate warming scenario. This requires incorporating temperature response functions of vector and pathogen traits into the dynamical model. Doing so, however, does not constitute an additional burden on the investigator since one needs to perform the same exercise in order to obtain R₀ predictions under warming (Mordecai et al., 2013, 2017; Johnson et al., 2015; Shocket et al., 2018, 2020; Tesla et al., 2018; Cator et al., 2020). As we have shown here, fitting mechanistic temperature response functions that accurately capture temperature effects on the biochemical processes underlying vector and pathogen traits does not require more data than the typical approach of fitting phenomenological functions. Moreover, it affords the advantage that parameters have clear biological meaning, and can be compared across different diseases as well as different studies of the same disease.

Previous studies attribute warming-induced declines in disease spread and prevalence to vector mortality. For instance, (Yamana et al., 2016) combined a decade of field observations with climate models to show that warming reduced mosquito survival and lowered the Malaria burden in west Africa. Erickson et al. (2012) found that warming facilitated the spread of dengue to higher latitudes (Illinois) by increasing season length and vector (Aedes albopitucs) period, but deterred it at lower latitudes (Texas, Georgia) due to hotter summers increasing larval mortality. The net outcome was a reduction in overall mosquito abundance, and a change in the abundance pattern from a single summer peak to a bimodal distribution with a spring and autumn peak. Morin and Comrie (2013) predicted a steeper summer decline of the west Nile virus vector (Culex sp.) for the central and Gulf states due to the longer and more intense summer warming projected for these regions, and a shallower decline for northerly locations and high elevations experiencing fewer high temperature extremes. Our trait-based analysis, which differs from these previous studies in using experimental data to parameterize mechanistic trait response functions, shows that when warming reduces infected host and vector populations it is through a complex interaction between adult mortality—the effects of which become compounded as the vector transitions from uninfected to infected—and the warming-induced decreases in the pathogen maturation rate, the vector's biting rate, and the per bite vector-to-host transmission probability. These findings highlight the crucial role that the pathogen's maturation rate and the vector's transmission traits play in driving warming effects on vector-borne diseases.

Our finding that uninfected vector populations can both persist at lower latitudes despite warming, and spread to higher latitudes because of it, has important evolutionary implications for predicting future disease burden. While the currently observed high vector mortality, combined with the reduction in maturation and transmission, may indicate a lowering of disease burden from lower latitudes and a curtailment of its spread to higher latitudes, this expectation does not accommodate vector adaptation to warming. Provided uninfected vector populations are not too small and the requisite genetic variation exists, it is entirely possible that the vector species' thermal reaction norms may evolve to withstand warmer climates and to spread to cooler climates. Evidence for rapid evolution of insecticide resistance in mosquitos and other arthropod vectors (Roberts and Andre, 1994; Hemingway and Ranson, 2000; Zaim and Guillet, 2002; Liu, 2015; Labbé et al., 2017) combined with evidence of genetic diversity within and between species (Tabachnick, 2013; Gloria-Soria et al., 2016; Powell, 2018) caution against being overly optimistic about warming-induced decreases in disease burden in lower latitudes. The same is true of the lower thermal limit for disease persistence. There is ample empirical evidence of rapid acclimation and adaptation of the lower critical temperature in insects (Storey and Storey, 1992; Addo-Bediako et al., 2000; Somero, 2010; Overgaard et al., 2011; Sunday et al., 2012; Araujo et al., 2013; Olson et al., 2013; Pintor et al., 2016), raising the possibility of greater disease spread to cooler climates than is expected based on ecological models alone. There is a crucial need to develop eco-evolutionary models of vector-host dynamics that can accommodate both thermal acclimation and adaptation.

Lastly, our analyses clearly show that using a DDE model with temperature- and time-dependent developmental delays provides the most biologically realistic depiction of vector and pathogen biology, and is therefore the best candidate for predicting disease spread and prevalence under climate warming. The problem, though, is that numerical implementation of such models is difficult and error-prone, unless the appropriate history conditions are correctly specified and the integral functions that calculate juvenile mortality during development are correctly implemented in a computer program (Murdoch et al., 2003). These challenges put DDE-based models out of reach for most ecologists and epidemiologists. The use of ODE equivalents presents an alternative that is more easily numerically implemented, but the representation of delays in terms of instantaneous maturation rates can lead potential under-prediction of disease spread and prevalence at high temperature extremes. The fact that the ODE models allow for calculations of both spread and prevalence in variable thermal environments make them more reliable than the sole use of R₀ formulations derived under constant-temperature assumptions. However, the reliability of ODE models needs to be verified by comparing predictions from time-varying DDE models with their ODE equivalents. While this is beyond the scope of our current study, it is an important future direction to pursue.

In conclusion, our work constitutes one of the first attempts at developing a mechanistic theoretical framework that integrates life history theory, thermal biology, population dynamics and epidemiology to investigate climate warming effects on vector-borne diseases. Its novelty lies in elucidating mechanisms from first principles, rather than inferring mechanisms from the data. We incorporate mechanistically-derived temperature response functions, parameterized with temperature response data, into dynamical models to make predictions about disease spread and prevalence that are completely independent of population-level (e.g., time series) data. This provides for greater degrees of freedom and statistical power in comparing time series of prevalence predicted from the model with equivalent data for any location. Our mechanistic approach has the potential to provide reliable information that can be used in formulating cost-effective public health policies.

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

PA designed the project and wrote this article with input from AE. AE and PA conducted the analyses. Both authors contributed to the article and approved the submitted version.

This research was supported by a National Science Foundation grant DEB-1949796 and a Guggenheim Fellowship to PA.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.