Estimating Past, Present, and Future Trends in the Global Distribution and Abundance of the Arbovirus Vector Aedes aegypti Under Climate Change Scenarios

Background: Aedes aegypti is the principal vector for several important arbovirus diseases, including dengue, chikungunya, yellow fever, and Zika. While recent empirical research has attempted to identify the current global distribution of the vector, the seasonal, and longer-term dynamics of the mosquito in response to trends in climate, population, and economic development over the twentieth and the twenty-first century remains to be elucidated. Methods: In this study, we use a process-based mathematical model to estimate global vector distribution and abundance. The model is based on the lifecycle of the vector and its dependence on climate, and the model sensitivity to socio-economic development is tested. Model parameters were generally empirically based, and the model was calibrated to global databases and time series of occurrence and abundance records. Climate data on temperature and rainfall were taken from CRU TS3.25 (1901–2015) and five global circulation models (CMIP5; 2006–2099) forced by a high-end (RCP8.5) and a low-end (RCP2.6) emission scenario. Socio-economic data on global GDP and human population density were from ISIMIP (1950–2099). Findings: The change in the potential of global abundance in A. aegypti over the last century up to today is estimated to be an increase of 9.5% globally and a further increase of 20 or 30% by the end of this century under a low compared to a high carbon emission future, respectively. The largest increase has occurred in the last two decades, indicating a tipping point in climate-driven global abundance which will be stabilized at the earliest in the mid-twenty-first century. The realized abundance is estimated to be sensitive to socioeconomic development. Interpretation: Our data indicate that climate change mitigation, i.e., following the Paris Agreement, could considerably help in suppressing risks of increased abundance and emergence of A. aegypti globally in the second half of the twenty-first century.

. Model A framework for Aedes aegypti mosquito population dynamics. Colour codes are: shaded green is precipitation dependent or under-water stages, red coloured words temperature dependent. Three vector populations are in boxes. Model parameters are those boxes without borders, which describe all the vital rates, and environmental influences through the carrying capcity.
All the mosquitoes' vital rates (birth, death, and transition) depend on temperature and/or precipitation. In addition, the egg-to-larva hatching fraction depends on the larvae population and environmental carrying capacity due to competition for survival. The whole lifecycle from egg to adult takes one to a few weeks on average depending on temperature and precipitation.
Three ordinary differential equations are used to describe the rates of change of individuals for each of the life-history stages: Here L, P, and A are the total female larval, pupal, and adult population per area, with the area being each 0.5 x 0.5 degree grid from which we generate our global maps. All the variables and parameters used in these equations are explained in Table S1-3. From the top down, the three equations describe the rates of change of their respective densities. Female adults lay eggs at rate of which a fraction ( ) are viable, of which a fraction f become female. Larvae die at rate and develop into pupae at rate l . Similarly, pupae die at rate p and emerge as female adults at rate p . Adults die at rate a .
The egg-to-larva hatching fraction is expressed as the product of two factors: the fraction of eggs hatching q at low larval population density and ( ) = 1 − / the larval density-dependent reduction factor for the viable eggs to survive and develop into larvae if larval population density is under the environmental carrying capacity C. In principle, L can be greater than C in which case we consider the hatching rate to be zero, but it did not happen in any of the cases we consider.
Model B includes human contributions to both the larval breeding sites and blood meals for the mosquito's reproduction, the fecundity rate of adult females as the shaded yellow colour ( Figure S1B). GDP is also included in contributing to the mosquito's larval sites. Figure S1B. Model B framework for Aedes aegypti mosquito population dynamics. Shaded yellow indicates human population or activity related parameters. Three vector populations are in boxes. Additional model parameters (shaded yellow) are used to describe the environmental and human influences.
In Model B, only Equation (1A) is changed to account for human population and gross-domestic product (GDP) per capita: Here, one additional factor is added to the first term to describe the eggs laid by the female adults at rate ℎ. Larvae die at rate and develop into pupae at rate l . Similarly, pupae die at rate p and emerge as female adults at rate p . Adults die at rate a .
The adult female fecundity rate ℎ is expressed as the product of two factors. The first is the oviposition rate as determined in laboratory studies when sufficient blood meals were provided. Under natural conditions, the availability of blood meals is limited by human population density ρ, as described by the second factor h (human blood-meal factor), which describes the reduction in oviposition rate due to lack of available blood meals.
The colour codes are similar to the framework in Figure S1, with red representing temperaturedependent variables, green rainfall-dependent variables, yellow human-dependent parameters or variables, grey vector population density, and white constant or mixed variable types. We used a pseudo-unit [M] for mosquito numbers (2) with all other units conforming to the SI system.

2
Variables and their relations with weather, human population and GDP

Parameterization
The variables and functions used in Equations (1-3) are summarized as equations (4)(5)(6)(7)(8)(9)(10)(11), see Table  S1A for variables used in Model A and Table S1B for Model B. The parameters and values used in these equations are summarized in Table S2 and S3. Each function is motivated and explained here.
The human blood meal factor h, Equation (4) is used based on the fact that in nature, except for rainforests, humans are the main food that Aedes aegypti require for laying eggs (3). ρ1 corresponds to the human population density at which half blood meal is obtained. The value was taken as the minimum human density within the flying distance of Aedes aegypti, one person per 50 x 50 meter (Table S1B) (2). Equation (6) is the fraction of eggs that hatch to larvae, q, in the situation where there are sufficient larval sites for hatching. This under-water process takes place from both the rainfalldependent larval site (first term) and the rainfall-independent larval sites (q2see Table S2) (2).
Larval population is regulated by the environmental carrying capacity (C). C is expressed as the product of two factors, C=cd. Here, c as shown in Equation (7) is the carrying capacity per larval site and is mainly contributed by rainfall (W) with a small fraction from sources independent of rainfall (c2) (2). d as shown in Equation (5) is the number of larval sites that comes from human contribution, directly or indirectly through agriculture, polyculture and urbanization (4). Directly, humans create the untreated water containers that catch rain water, i.e., for water storages used for drinking, farming, and husbandry, or used tires and waste containers (4,5). In general, we assume that the higher the human population in a given area (ρ), the more the larval sites (d). The lower the economic level, the more the water containers for larval site. We assume that d is proportional to the human population and inversely proportional to the economic level, which we approximate by the square root of the ratio of GDP/capita ( ) to the world average value ( ). In Equation (5), the proportional coefficient, b, is used to represent the number of larval sites per person if the GDP/capita is at the world average value. b is determined by comparing the model output with the field data for Juaweiro, Brazilsee Table S2.
Equations (9-10) describe the mortality rates of larva and pupa. They depend on both temperature and precipitation since they are under-water stages. Heavy rainfall washes out the larval site through overflow and creates extra mortality if it exceeds the critical value Wc -30 mm/day (2). In the equations, θ(x) is the Heaviside function which is 1 if x≥0 and 0 otherwise. The temperature dependent part of the mortality μ(T) as well as the development/transition rate (σ) between compartments and the intrinsic oviposition rate ( ) are estimated from laboratory studies of Aedes aegypti (6). They are fitted using n-th degree polynomial as shown in Equation (8) with coefficients shown in Table S3, where y(T) stands for either μj(T), σj (T) or (T), and j stands for either l, p, or a, denoting the three stages of the mosquito's life.   (2015) (7). Table S3. Coefficients (pi) in temperature-dependent vector parameters obtained by fitting polynomial (2). The unit of coefficients pi is days −1 × (°C) −i . Figure S2 illustrates the method used in this study as a flow diagram from data input to model output.

Implementation
Results are presented as global maps and trend. The top row shows the process. Here the vector model was developed in two steps: Model A, climate driving the vector population ( Fig 6A); Model B, climate, human population and GDP driving together (Fig 6B).  Table 2 in Ref.
(2)* c 1 amount of rain to produce 1/2 of the max breeding sites mm/day 30 Table 2 in Ref.
(2)* q 0 amount of rain to allow 50% of eggs hatching mm/day 0.2 Table 2 in Ref.
(2)* q 2 hatching fraction from rain-independent sites 0.037 µ la additional mortality for larvae due to heavy rain day/mm 0.001 ** µ pa additional mortality for pupae due to heavy rain day/mm 0.001 ** W c critical rain volume to cause additional mortality mm/day 30 Table 2 in Ref. In Model A, only climate data are used as the input with two variables: mean temperature and precipitation. These two monthly datasets are then interpolated to generate continuous climate data.
Weather relations of eight model parameters are obtained from the literature (Equations (6-7, 8-10)) with some improvement to be consistent with the assumptions made. Next, the climate data are linked to these model parameters before entering into the model. Six of the eight model parameters depend on temperaturesee Figure S1. Precipitation (rainfall) contributes to the under-water stage of vector development, from mortality rate, egg-hatching fraction to environmental carrying capacity. Here the environmental carrying capacity for larvae is affected by the larval sites created by both rainfall, nature, and human activities (see Table S1B).
In Model B, human population and GDP are two extra datasets used as input into the model. They describe the socioeconomic factor that contributes to the vector's development through two model parameters (yellow boxes). First, they both contribute to the number of larval sites, assuming a direct proportion of the number of larval sites to human population and an inverse proportion to GDP per capita (Equation (5)). Second, human population provides blood meals to the urban vector; this is described as the human blood meal factor (Equation (4)). Like the climate data, the socioeconomic datasets are interpolated and linked to these two model parameters before being entered into the Model.
From the Climate Research Unit (CRU) online database, time series (CRU TS3.25) of gridded (0.5 x 0.5 degrees) monthly mean temperature and precipitation were downloaded from January 1, 1901 to December 31, 2016 (8). From ISMIP (9), we downloaded the future climate data from five global circulation models (CMIP5) (10) for RCP2.6 and RCP8. 5 (2006-2099), the yearly global human population  and GDP (1950GDP ( -2099 data. Both are gridded data at 0.5° and 5' resolutions. We use linear interpolation function to obtain continuous data before solving the equations. Programs Wolfram Mathematica 11 and R (package deSolve 1.20) (11) are used to obtain time series vector population and abundance values and to generate global maps.
To obtain the vector populations, we solve differential equations (Eq. (1-3)) numerically using the two computer programs mentioned above. From this, we summarize the total female adult vector population over a defined period to get vector abundance potential from Model A output or abundance density from Model B output. Finally, based on averaged seasonal or annual vector abundance over a decade, we map the vector abundance globally. The results from the two models A and B are compared to see the additional contribution to vector abundance from socioeconomic factors in addition to climate. In addition, we regressed the percentage change in global vector abundance potential against changes in global mean temperatures per decade over two centuries (1901-2099) using linear regression. The future two climate scenarios were used, RCP 2.6 and RCP 8.5. The resulting regression coefficient gave the average percentage change in vector abundance potential in relation to unit increase in temperature while the R 2 provided the correlation coefficient.

Optimization of model parameters
Two new model parameters were used in Model B to generate vector population density besides those used for Model A: b -the number of larval sites for larva at a given human population density if GDP/capita is the same as the world average and ρ1the half-saturation constant to allow gravid female mosquitos to feed human blood for laying eggs. Based on the published field study on Aedes aegypti population (7), we have calibrated these two parameters by comparing the model output of female adult with the field data on the same capture days (Table S2). See Figure S3 for the comparison of field data and our model 2 output.

Sensitivity analysis of model parameters
We have used nine time-varying parameters in the model B. Besides the two parameters related to human population and socioeconomics as described above, the rest seven parameters were climate related. They are listed in Table S1B as Equations 4, 6, 7B, 8-10. Some of these parameters were measured in laboratories (6) and some were assumed and validated to local dengue outbreak cases (1,2). Here we perform sensitivity analysis to see the effect of each parameter to the output of Model Bthe abundance of the adult vector population annual average over the recent 10 years for six locations.
Each parameter is varied 1000 times sampling around its original value following a normal distribution with sigma of 0.05. We used R-package (pse v.0.4.7) to run the model (14). The results were generated as empirical cumulative density functions and partial rank correlation coefficients (PRCC) between a set of model input parameters and the output of Aedes aegypti adult abundance. Figure S4 shows PRCC for six cities across different climate and locations. The top panel shows three areas with Aedes aegypti endemicity; the low panel corresponds to relatively mild climate areas where Aedes aegypti has been found in Madeira since 2005 (15) and Phoenix recently (16).
The result shows that adult mortality rate, μm, is consistently the most sensitive parameters in affecting the abundance of adult population. Environmental carrying capacity, C, for larvae and the development rate from larva to pupa, σl, are two important parameters for all areas except Madeira and Malaga where vector population is low due to low temperatures when using CRU TS3.25 dataset without coastal temperature corrections applied. In addition, three more parameters are sensitive for the low panel but not the upper one: the oviposition rate θm, the egg hatching fraction q, and h, the human blood meal factor.
Therefore, we have found that sensitivity of parameters to the adult abundance depends on the level of vector abundance. The high abundance areas is sensitive to the limiting factor from environment carrying capacity, while low abundance areas is sensitive to the oviposition rate and hatching fraction. Overall, the adult vector mortality rate is the sensitive parameter to the adult population, as expected. θm -oviposition rate, σl -development rate from larva to pupa, σp -development rate from pupa to adult, μl, μp, μmlarval, pupal and adult mosquito mortality rate, qeggs hatching fraction, C -environmental carrying capacity and hthe human blood meal factor.

Validationcomparing global map of invasion vs. occurrence data
Validation of our model has been carried out through comparing the model output -Aedes aegypti decadal abundance predictions -to the reported data on global occurrence of Aedes aegypti (18). As shown in Figure S5 To a large degree, we find a good overlap between our Aedes aegypti abundance predictions and empirical observations of its occurrence in Asia, America and Australia. The overlap is better for the vector density output from Model B than the abundance potential from Model A. The discrepancy is more in Africa and southeast Asia, as well as Northwest of Australia. Our density model output matches Brazil very well but not in Africa and Southeast Asia. This is probably due to incomplete monitoring and data collection system for the Aedes aegypti data. The monitoring system is much better in Brazil than Asia and Africa. Therefore, we believe that our prediction is reasonable. In Australia, the discrepancy can be due to other reasons, such as lacking introduction of vectors and their spontaneous die out. Those factors are not included in our deterministic mechanical models.
The differences between Fig. S5A and S5B indicates the human contribution to the Aedes aegypti population growth. Fig. S5A is from Model A which is climate driven. Fig. S5A showing the vector population potential captures largely the occurrence data globally except some small areas such as part of Australia and middle of USA. This means that vector population development is largely controlled by climate. Figure S5. Comparison of model outputs on global map of Aedes aegypti abundance potential (A) and density (B) (continuous colors) and occurrence data (black dots) (18). Different drivers were used in the models: using climate only (Fig. S5A -Model A) and climate, human population and GDPpc (Fig. B -Model B). Parameter inputs were based on interpolated CRU TS3.25 monthly mean temperature and precipitation data. Female Aedes aegypti population potential abundance per larval site ( Figure S5A, linear scale) and density ( Figure S5B, log scale) were averaged over the period 2006-2015. Figure S6 shows the percentage of observed occurrence areas at each category of potential abundance predicted by the model A. We found that within the occurrence data area, 97.5% of the global occurrences are within areas with potential abundance predicted above 300 and about 55% of occurrence data lies within abundance potential between 900-1200. Therefore, our model A predicted potential abundance agrees well with global observations of vector occurrence data.
The occurrence data need to be interpreted with caution as they are not systematically sampled. Thus, more data points in a specific region is likely to reflect the frequency of measuring and reporting rather than the actual occurrence. See Kraemer et al for more information and discussion on data collection and validity (18). Figure S6. Distribution of the occurrence data 18 observations by categories of global potential abundance. Parameter inputs were based on Interpolated CRU TS3.25 monthly mean temperature and precipitation data and potential abundance was averaged for the period 2006-2015.