Effects of Mosquito Biology on Modeled Chikungunya Virus Invasion Potential in Florida

Arboviruses transmitted by Aedes aegypti and Aedes albopictus have been introduced to Florida on many occasions. Infrequently, these introductions lead to sporadic local transmission and, more rarely, sustained local transmission. Both mosquito species are present in Florida, with spatio-temporal variation in population composition. We developed a two-vector compartmental, deterministic model to investigate factors influencing Chikungunya virus (CHIKV) establishment. The model includes a nonlinear, temperature-dependent mosquito mortality function based on minimum mortality in a central temperature region. Latin Hypercube sampling was used to generate parameter sets used to simulate transmission dynamics, following the introduction of one infected human. The analysis was repeated for three values of the mortality function central temperature. Mean annual temperature was consistently important in the likelihood of epidemics, and epidemics increased as the central temperature increased. Ae. albopictus recruitment was influential at the lowest central temperature while Ae. aegypti recruitment was influential at higher central temperatures. Our results indicate that the likelihood of CHIKV establishment may vary, but overall Florida is permissive for introductions. Model outcomes were sensitive to the specifics of mosquito mortality. Mosquito biology parameters are variable, and improved understanding of this variation will improve our ability to predict the outcome of introductions.


Introduction
Chikungunya virus (CHIKV) is an Alphavirus in the family Togiviridae, native to Africa, which causes human disease and is transmitted by Aedes mosquitoes, primarily invasive Stegomyia. Zoonotic cycles between sylvatic Aedes species and non-human primates occur in Africa [1], with spillover to the human population, but elsewhere the transmission cycle is primarily between Ae. aegypti or Ae. albopictus (in its invasive range) and humans. Aedes aegypti is considered the most common vector, but outbreaks have occurred where Aedes albopictus is the dominant vector, such as in La Reúnion in 2005, India (2006, and Italy (2007) [1,2]. Three major lineages are recognized (East, Central, and South African (ECSA), Asian, and West African (WA)), with an additional sublineage of ECSA (Indian Ocean lineage, IOL). The genotypes of these lineages have diverged, and virulence in humans varies among lineages [3]. In addition, a mutation in the envelope protein (E1-A226V) that enhances transmission and reduces the extrinsic incubation period (EIP) in Ae. albopictus was found in the ECSA lineage [4], but the lineage is heterogeneous for the presence of the mutation. Globally, CHIKV transmission is ongoing [5] but can increase to outbreak levels and invade new areas. In 2014, there were high levels of CHIKV transmission in the Americas and the Caribbean, leading to concern about introduction of the virus to new areas. The risk of introductions and outbreaks continues, as illustrated impacts the probability of survival through the extrinsic incubation period, considering other complex, biologically relevant relationships is important in addressing the risk of invasion in diverse locations.
The distribution of Ae. aegypti and Ae. albopictus across Florida is highly variable in both space and time [28][29][30][31][32]. The two species share many characteristics, including larval development in a variety of small water-holding containers such as discarded tires or plant pots. However, they also differ in aspects important to arbovirus transmission, such as host preference [33][34][35]. The domesticated form of Ae. aegypti prefers to feed on humans (e.g., [33,34]) while Ae. albopictus is more opportunistic, but will readily feed on humans (e.g., [35][36][37]). Both are competent vectors for DENV, Zika virus (ZIKV), and CHIKV, although there are population and strain differences in virus-mosquito interactions. With this background, the risk of establishment and local transmission in Florida is expected to be highly variable, and a better understanding of the factors affecting this risk is critically needed. Process-based epidemiological modeling is an effective way to consider many combinations of factors to assess their contribution to outcomes such as local transmission and number of infected people. A two-vector WNV model for Florida demonstrated that the presence of two vector species can influence West Nile virus dynamics [19]. However, that study focused on species with different seasonal patterns but similar biology and transmission, while Aedes-transmitted arboviruses in Florida are likely to be more similar in seasonal patterns but have differences in biology and transmission. Given the variation in the major CHIKV vectors in Florida, understanding potential invasions requires inclusion of both vectors and varying relative abundance in models to assess how these factors influence arboviral dynamics.
To address questions about the likelihood of successful establishment of CHIKV into Florida against a background of fluctuating populations of two mosquito species, we adapted a previous two-vector arbovirus model [19]. This model allows different population sizes, both absolute and relative, of Ae. aegypti and Ae. albopictus, CHIKV strain variation, and includes a nonlinear relationship between adult mosquito mortality and temperature. We used the model to investigate factors influencing risk of invasion of CHIKV into different locations in Florida, and to consider the impact of fixed parameters on model outcomes.

Model Structure
The basic model structure follows an earlier deterministic two-vector, multiple host arbovirus model [19], which was developed from one-vector, two-host arbovirus models [17,18]. Two vectors were retained but the host structure was simplified to consider only one host type, as there are no known alternate hosts for CHIKV in Florida ( Figure 1). The model is SIR (Susceptible-Infectious-Recovered) for hosts (humans) and SEI (Susceptible-Exposed or latent-Infectious) for vectors. Virus genotype differences are included by altering relevant parameters. Below, we detail the changes made to the structure, equations, and parameters of the model to consider CHIKV introductions into Florida. Other model structure and equations can be found in the references above. The model has a time step of 1 day and is run for 730 days following virus introduction. A combination of stratified random sampling and systematic variation was used for sensitivity analysis, detailed below. Model parameters and ranges used in the sensitivity analysis are given in Table 1.  The previous models considered only one location, while the intent of this model was to consider risk of invasion across Florida. Temperature is modeled with an annual cycle of daily mean temperatures. We analyzed historical yearly means and amplitudes and determined that, for this set of weather station data in Florida, means and ranges are correlated ( Figure S1, r = −0.96). Therefore, we used the annual mean temperature (T mean ) as a parameter (Table 1) in sensitivity analyses, then calculated T range using a linear equation fit to the data, daily temperature was then calculated as Day 0 of the simulation is 1 January; the offset in the cosine function (26.5) aligns the temperature function maximum with historical annual peak temperatures. T mean was varied to reflect spatial and temporal variation in Florida (Table 1).

Human Population
To focus attention on the dynamics of the mosquito-virus interaction, the human population was simplified. Human populations are relatively stable over short time periods, so the human population size was fixed at the initial value (H tot ) and there were no births, immigration, mortality, or emigration. Although neuroinfection and deaths have been reported [38], CHIKV frequently has a low fatality rate [39] so disease-induced mortality was not included. Human population size does vary in communities across Florida, and the size of the human population interacting with a mosquito population is often unknown and dynamic. To focus attention on the differences in the mosquito populations, human population size was fixed at 50,000 in this analysis. The mosquito population size (see below) is described as mosquitoes per human host per year. No latent period in humans was included, as this was thought to be relatively short [39] and, with no human mortality, would act only to slightly slow or delay the dynamics. The duration of infectiousness in humans may extend Viruses 2020, 12, 830 6 of 21 beyond the symptomatic period (approximately 1 week [40]), and more recently has been suggested to last longer than in other arboviruses [41]). In addition, recent work suggests variability between people and PCR detection of virus for 4-16 days [42]). To focus attention on the mosquito elements of the model, for this analysis duration of infectiousness was set at 8 days.

Mosquito Population Dynamics and Seasonality
Mosquito dynamics are highly variable, with both environmental and biotic influences. To focus on the mosquito-virus interactions, we simplified mosquito dynamics. We assumed that water availability would be the primary driver of recruitment to the adult population, and that both species would respond similarly. This model forces seasonal recruitment using repeated year-round pulses (iv days apart) and two large peaks (~N(q i , σ i ); i = [1,2] for spring, summer peaks) in adult recruitment ( [19], Figure S2). Abundance is specified as ρ max,j , the number of adult female mosquitoes of species j per human host over the entire year. The total number of mosquitoes of species j is then given by R tot,j = H tot *ρ max,j where H tot is the number of humans, and is partitioned between the pulses and peaks by the proportion in the baseline pulses (p base,j ) and in the first peak (p δ1,j ); the proportion in the second peak is thus Recruitment to the adult population was thus (for each mosquito species j, j = [aegypti, albopictus]):

Mosquito Mortality
The previous models used a simple linear increase in mortality with temperature over a relatively narrow temperature range. To include a wider temperature range, increases in mortality at colder temperatures needed to be considered. One study showed relatively little change in mortality over a wide temperature range for Ae. aegypti [27], increasing outside the optimal range. Laboratory studies may underestimate mortality, and this study did not extend to temperatures likely to be immediately lethal [27]. The authors fit higher order polynomials to their data, with the best fit with a 4th order polynomial. However, the parameters of this polynomial equation are not easily interpreted biologically. We modeled temperature dependent mortality for both species with a piecewise linear function, resulting in low mortality in an optimal survival temperature range and linear increases on each side. The minimum mortality (µ min, j ) and slope with respect to temperature (µ sl,j ) varied separately for each species (j = [1, 2]) but the center of the optimum survival temperature range (Temp c ) did not. The optimum survival range was narrow (Temp c ± W; W = 2 • C) to allow more variation in the slope while still reaching high levels of mortality ( Figure 2). The mortality rate was bounded (µ min , 1) although the maximum of 1 may not be reached under all parameter combinations. Thus, Viruses 2020, 12, 830 7 of 21 Viruses 2020, 12, x 7 of 21

Virus Development in Mosquitoes and Vector Competence
The mutations in CHIKV that allow increased vector competence in Ae. albopictus affect the EIP (rate of viral development) and the proportion of insects that can acquire and transmit the virus [4]. Virus development in mosquitoes, or the duration of the EIP, is well known to be temperature dependent (e.g., DENV, [16,43,44]; West Nile virus, [45]; CHIKV, [46,47]), and there was evidence that virus development is faster for CHIKV than for some other Aedes-transmitted arboviruses [46,48]. We retained the previous structure of a linear increase with temperature and added a lower bound of 10 °C, below which no development occurs. Parameters for this equation vary between mosquito species and are given by γ22,j and γsl,j, virus development at 22 °C, and the slope of the line, respectively. The EIP is then 1/γj and thus is temperature dependent and can vary between mosquito species. The ability of the mosquitoes to acquire (bj) and transmit the virus (βj) is species-specific and varied independently. In this analysis, the vector competence parameters were not correlated, to assess their independent effects on dynamics.

Model Equations
The model structure consists of nine differential equations following changes in susceptible, infected, and recovered humans and susceptible, exposed (latently infected) and infectious mosquitoes of each species (j): As detailed above, mosquito recruitment, mortality and virus development are temperature and thus time dependent but not formally noted in Equation (8) as such for clarity.
In a deterministic model, a virus can persist at extremely low levels where, in a stochastic situation, it would have a high probability of going extinct. To address this, if the combined virus prevalence went below 0.5 (i.e., + ∑ ( + ) < 0.5 ) the virus was assumed to die out and all state variables including infection reset to 0.

Virus Development in Mosquitoes and Vector Competence
The mutations in CHIKV that allow increased vector competence in Ae. albopictus affect the EIP (rate of viral development) and the proportion of insects that can acquire and transmit the virus [4]. Virus development in mosquitoes, or the duration of the EIP, is well known to be temperature dependent (e.g., DENV, [16,43,44]; West Nile virus, [45]; CHIKV, [46,47]), and there was evidence that virus development is faster for CHIKV than for some other Aedes-transmitted arboviruses [46,48]. We retained the previous structure of a linear increase with temperature and added a lower bound of 10 • C, below which no development occurs. Parameters for this equation vary between mosquito species and are given by γ 22,j and γ sl,j , virus development at 22 • C, and the slope of the line, respectively. The EIP is then 1/γ j and thus is temperature dependent and can vary between mosquito species. The ability of the mosquitoes to acquire (b j ) and transmit the virus (β j ) is species-specific and varied independently. In this analysis, the vector competence parameters were not correlated, to assess their independent effects on dynamics.

Model Equations
The model structure consists of nine differential equations following changes in susceptible, infected, and recovered humans and susceptible, exposed (latently infected) and infectious mosquitoes of each species (j): As detailed above, mosquito recruitment, mortality and virus development are temperature and thus time dependent but not formally noted in Equation (8) as such for clarity.
In a deterministic model, a virus can persist at extremely low levels where, in a stochastic situation, it would have a high probability of going extinct. To address this, if the combined virus prevalence went below 0.5 (i.e., H in + j (L j + Y j ) < 0.5) the virus was assumed to die out and all state variables including infection reset to 0.

Parameter Values and Sensitivity Analysis
Latin hypercube sampling (LHC), a stratified random sampling method (e.g., [17,18,49]) was used for sensitivity analysis. Briefly, for each parameter of interest a probability density function (pdf) is chosen. The probability axis is divided into Z equi-probable intervals and one value chosen from each interval. These probabilities are then back-transformed into parameter values. Values for each parameter are randomly resorted into Z parameter sets used to run the model; here Z = 250. We used uniform pdfs where little information was available or the intent was to extend over biologically or environmentally feasible ranges, defining the maximum and minimum values. Triangular distributions were used when data were available to identify a most probable value (peak of the distribution) along with maximums and minimums. A preliminary LHC analysis was conducted and showed very high virus activity; parameters were adjusted to reduce outbreak activity to allow analysis and reflect observations of infrequent establishment (few observed local transmission events relative to observed travel-associated imported cases). Distributions and parameter values are given in Table 1 and discussed below, along with fixed values for parameters not varied in this analysis. In addition to this analysis, Temp c (center of optimal survival range) was varied systematically, (10,16,22) repeating the same LHC parameter sets with each Temp c value (for a total of 750 simulation runs). For each LHC set (each value of Temp c ), epidemic presence or absence was analyzed using multiple logistic regression and the size of the epidemic (MaxH i ) analyzed using multiple regression. Parameters were ranked by the p-value of the coefficients and the top five parameters identified as factors of interest and compared between Temp c values.
LHC and other variations on randomized sampling are unlikely to include sets of parameter values of interest, for example comparing outcomes when a virus strain adapted to one mosquito species is introduced into a mosquito population dominated by that species or another. To consider these events, eight sets of parameter values (see Figure 6 for parameter values) were developed to compare introductions of specific virus strains into defined mosquito populations. For this analysis, human population size was reduced to 5000. Adaptation between viral strain and mosquito species was defined as higher vector competence (b j and β j ). To limit the parameters considered, EIP was not varied in this comparison. Mosquito populations with varied proportions of the two species and selected mosquito biology parameters, either identical or species-specific (a j , µ min,j , µ sl,j ), were modeled and the outcome of introductions assessed. The size of the epidemic was compared between parameter sets. Baseline model runs had mosquito biology parameters set identically between the species. The baseline outcomes were compared to outcomes with species-specific population sizes and vector competence.

Mosquito Population Dynamics and Seasonality
We assume the two species respond similarly to water (rainfall, human generated water) and temperature for hatching and larval development and thus the seasonality of recruitment to the adult population will be similar. However, population abundance of both species varies across space and time in Florida (e.g., [29,31]) and one goal was to consider how this spatio-temporal variation in relative population sizes affects the likelihood of local transmission. Therefore, parameters describing seasonality (iv, q i, σ i (i = [1, 2]) were set equal for the two species, while parameters affecting relative abundance were allowed to vary between species. Diapause was not considered, although diapause in Ae. albopictus varies with latitude [50]. The two seasonal peaks occurred on days 165 and 245 (q 1 , q 2 , respectively) with σ i = [7,15]. To consider variation in wetting events in the dry season, iv was included in the LHS analysis and ranged from 10 to 50 days. Reflecting the rainy season in Florida (approximately June through October) but also considering human generated water in containers (e.g., irrigation, watering plants) and differential response of the two species to desiccation [51], the proportion of mosquito recruitment in the year-round baseline pulses was fixed at p base,j = [0.13, 0.15] [j = aegypti, albopictus] and in the spring peak at p δ1,j = [0.25, 0.2]. Recruitment to each population, ρ max,j was sampled independently with a triangular distribution, over a wide range (100-20,000 mosquitoes/species/host/year). This resulted in both the total number of mosquitoes and the relative Viruses 2020, 12, 830 9 of 21 population size of the two species varying between runs in the sensitivity analysis. Note that although seasonality in recruitment was the same for the two species, temperature-dependent adult mortality was not (see above, [52]). Therefore, adult populations could show differences between species in seasonal abundance. For these and other mosquito biology parameters, the goal was to consider effects of variation expected across different environments and seasons and how differences in the two species may affect transmission.

Mosquito Mortality
Although there are numerous studies of mortality rates in both species (reviewed in [52]), the relationship between laboratory measurements and field mortality is poorly understood. We assume field mortality rates will be higher, and therefore extend the range of observed mortality rates to include higher values. Temperature dependence is well documented in the lab, although less so in the field. For both species, the minimum mortality was fixed (µ min,j = [0.1, 0.06]) with higher mortality rates in Ae. aegypti. The slope of mortality with respect to temperature (µ sl"j ) was included in the sensitivity analysis, with the same range (0.05-0.15) used for both species but values chosen separately, so in each of the LHS runs the species would differ in their sensitivity to temperature. We then systematically varied the central point of the temperature range with minimum mortality (constant between species) over three temperatures (Temp c = 10, 16,22), repeating the LHC parameter sets.
For the virus strain analysis, Temp c = 10 and 22 only. Minimum mortality was set higher in Ae. aegypti than in Ae. albopictus to reflect shorter expected life span, but temperature sensitivity was lower in Ae. aegypti as it is expected to be less sensitive to high temperatures (µ min,j = [0.08,0.06] and µ sl,j = [0.006, 0.08]).

Biting Rate
The biting rate, a j , (bites on humans per day) is the inverse of days between blood meals on humans, α j . Non-human hosts are not explicitly modeled but included by adjusting a j . Aedes aegypti can feed more frequently (e.g., [53,54]) and is more likely to feed on humans (e.g., [33,34]), while Aedes albopictus can be more opportunistic (e.g., [35]) and so may feed on humans less frequently. However, high variability between populations, particularly for Ae. albopictus, has been observed. Thus, α j is higher and has a wider distribution for Ae. albopictus than for Ae. aegypti, with triangular pdfs used for both to concentrate sampling near the most probable value (Ae. aegypti range 1-5, center 3; Ae. albopictus range 2-10, center 5). The model was constructed to facilitate later inclusion of explicit host preference [17,18] or seasonal shifts in host feeding behavior, but was constrained to a j = 1/α j for this analysis. LHC sampling was done on α j then back-transformed for simulation.
In the virus strain analysis, the shorter interval between blood meals for Ae. aegypti was retained, with α = 3d and for Ae. albopictus α = 4d. The difference between the species was decreased from the most probable values (above) to focus attention on the virus transmission parameters (below).

Vector Competence and Virus Development in Mosquitoes
Studies of the vector competence and virus development of these two species and CHIKV strains show variable results and a range of competency (10-80% transmitting virus; e.g., [4,46,55]). Some CHIKV strains have acquired mutations that increase the competence of Ae. albopictus (although Ae. aegypti may still be the more competent vector). The ability of the mosquitoes to acquire and transmit the virus (b j ) is therefore species-specific. There is less information available on the rate of human infection following infectious bites (β j ), so the range and most likely values for β j were set the same as b j but varied independently in the LHC sensitivity analysis. Preliminary analysis with ranges including higher transmission rates reported in the literature generated high levels of outbreaks, restricting analysis and not reflective of the low rate of local transmission following introductions. Accordingly, rates were adjusted down for this analysis (Table 1), while retaining wide ranges to explore the effect on transmission dynamics. Later work supports the wide range of transmission values, with geographic, species, and virus strain variation (e.g., [56]). Notable for CHIKV is the short extrinsic incubation period relative to other arboviruses (e.g., [4,46,48,55]). Virus development in mosquitoes is temperature dependent in other systems (e.g., DENV, [16,43]; West Nile virus, [45]). Limited information available during parameter estimation suggested CHIKV development would also be temperature dependent (e.g., [46]). Therefore, we assumed the CHIKV development rate was a linear function of temperature. Parameters for the linear function are γ 22 (virus development at 22 • C, fixed at 0.25 for both species) and γ sl (slope with respect to temperature, ranging from 0.004 to 0.02 with a most probable value of 0.015 for each species, chosen independently).
In the strain analysis, our goal was to explore the potential effects of virus-mosquito adaptations leading to higher transmission. Therefore, b, β parameters were varied to reflect different CHIKV strains The virus was introduced by one infectious human, on the day specified by t crjt , as this is considered a common route of introduction for CHIKV [57]. The day of introduction spanned most of the year with a uniform distribution, to consider the outcome of introduction at any point in the temperature cycle. This is plausible for introductions, as outbreaks elsewhere and travel occur throughout the year. Day of introduction was not varied in the strain analysis and was fixed at day 200.
The model was programmed in Matlab (The Mathworks, Natick, MA, USA) using the differential equation solver function ode45. Model analysis was conducted in Matlab using multiple regression functions to generate output needed. Matlab code for the ODE model is provided in the Supplemental Materials.

Outcome Variables
The peak number of human cases (MaxH i ) was determined for each simulation run and those runs where MaxH i > 2 were scored as epidemic (epi = 1) or not (epi = 0). Further outcome variables were analyzed for epidemic runs only. The time of the peak in humans (tMaxH i ) and the maximum number of humans infected (MaxH i ) were determined for each simulation. The lag between virus introduction and the time of peak virus activity was calculated as lag = (tMaxH i -t crit ).

Analysis
Logistic regression was used to determine the effect of input parameters on the presence/absence of epidemics in the LHC parameter sets, separately for each value of Temp c . Only main effects were considered, with the full regression model. Parameters were ranked using p-values and the five most influential parameters compared between Temp c values. Given the deterministic model format, all parameters will influence outcomes and so conventional statistical significance is not appropriate; p-values are a convenient way to rank the input parameters in importance. Multiple regression was used to assess input parameter effects on peak case numbers (MaxH i ) and lag for epidemic runs only, with the most influential parameters identified as for epi.

Epidemic Behavior
Most simulation runs demonstrated one of three behaviors: no outbreak and no local transmission occurred, a single outbreak occurred that infected all susceptibles in a relatively short time period ( Figure 3A), or two outbreaks occurred (usually relatively close together in time) that also ultimately infected all hosts ( Figure 3B). Timing of the outbreaks varied within these behaviors (Table 2). A few parameter sets showed intermediate or alternate behaviors. Changes in Temp c alone could change this behavior. For example, in one parameter set the outcome was a single epidemic fairly soon after introduction with Temp c = 22 • C ( Figure 3A). Decreasing Temp c to 16 shifted the outcome behavior to two outbreaks ( Figure 4A), while decreasing Temp c further to 10 recovered a single epidemic peak but with a delay in the time of the epidemic ( Figure 4B). The number of parameter sets exhibiting outbreaks increased with Temp c , from 75 (30%) to 178 (72%) ( Table 2).

Epidemic Behavior
Most simulation runs demonstrated one of three behaviors: no outbreak and no local transmission occurred, a single outbreak occurred that infected all susceptibles in a relatively short time period (Figure 3A), or two outbreaks occurred (usually relatively close together in time) that also ultimately infected all hosts ( Figure 3B). Timing of the outbreaks varied within these behaviors ( Table 2). A few parameter sets showed intermediate or alternate behaviors. Changes in Tempc alone could change this behavior. For example, in one parameter set the outcome was a single epidemic fairly soon after introduction with Tempc = 22 °C ( Figure 3A). Decreasing Tempc to 16 shifted the outcome behavior to two outbreaks ( Figure 4A), while decreasing Tempc further to 10 recovered a single epidemic peak but with a delay in the time of the epidemic ( Figure 4B). The number of parameter sets exhibiting outbreaks increased with Tempc, from 75 (30%) to 178 (72%) ( Table 2).  Table  S6. Tempc = 22 °C. Blue lines, susceptible; red lines, infectious humans or mosquitoes; and green lines, recovered humans. Mosquito panels set to the same Y-axis scaling for comparison; infectious mosquitoes not always visible at this scale.  Table S6. Temp c = 22 • C. Blue lines, susceptible; red lines, infectious humans or mosquitoes; and green lines, recovered humans. Mosquito panels set to the same Y-axis scaling for comparison; infectious mosquitoes not always visible at this scale.  Logistic regression with main effects only had low predictive ability for the likelihood of epidemics ( Table 2), suggesting that interactions are important in predicting model outcome. The parameters most influencing the likelihood of epidemics changed based on the value of the central temperature for the mortality function (Tempc) ( Table 3). The annual mean temperature was a key parameter in the likelihood of epidemics at all Tempc values, but the sign of the coefficient changed (Tables 3 and S1; Figure 5). When Tempc was lower (=10 or 16), decreasing annual mean temperature increased the likelihood of epidemics. Conversely, when Tempc = 22, decreasing annual mean temperature decreased the likelihood of epidemics. Increasing Ae. albopictus population size increased the likelihood of epidemics at the lowest Tempc but decreased it at the highest Tempc, while population size in Ae. aegypti increased the likelihood at middle and high levels. Overall, Ae. aegypti parameters were more important than Ae. albopictus parameters at all Tempc levels. Only Ae. aegypti transmission rates were in the top five, and only for low and middle Tempc levels. Similar effects were observed for biting frequency, but with a negative relationship indicating that lower biting frequency of Ae. aegypti increased the likelihood of epidemics at lower Tempc levels. The slope of the mortality rate-temperature relationship for Ae. aegypti was in the top five parameters for the middle and high Tempc levels, with increasing slopes decreasing the likelihood of epidemics. Logistic regression with main effects only had low predictive ability for the likelihood of epidemics ( Table 2), suggesting that interactions are important in predicting model outcome. The parameters most influencing the likelihood of epidemics changed based on the value of the central temperature for the mortality function (Temp c ) ( Table 3). The annual mean temperature was a key parameter in the likelihood of epidemics at all Temp c values, but the sign of the coefficient changed (Table 3 and  Table S1; Figure 5). When Temp c was lower (=10 or 16), decreasing annual mean temperature increased the likelihood of epidemics. Conversely, when Temp c = 22, decreasing annual mean temperature decreased the likelihood of epidemics. Increasing Ae. albopictus population size increased the likelihood of epidemics at the lowest Temp c but decreased it at the highest Temp c , while population size in Ae. aegypti increased the likelihood at middle and high levels. Overall, Ae. aegypti parameters were more important than Ae. albopictus parameters at all Temp c levels. Only Ae. aegypti transmission rates were in the top five, and only for low and middle Temp c levels. Similar effects were observed for biting frequency, but with a negative relationship indicating that lower biting frequency of Ae. aegypti increased the likelihood of epidemics at lower Temp c levels. The slope of the mortality rate-temperature relationship for Ae. aegypti was in the top five parameters for the middle and high Temp c levels, with increasing slopes decreasing the likelihood of epidemics.

Size and Lag Time of Epidemics
One measure of epidemic activity is the number of human cases at the peak of the epidemic (MaxHi). This is influenced by the model behavior and duration of virus activity, as some simulations showed single, short outbreaks ( Figure 3) while others had prolonged, lower viral activity. The average MaxHi increased with Tempc, but the range was large in all three sets (Table 2, Figure S3). Adjusted R2 values again were low, indicating that interactions are likely to be important (Table 2). Parameters most influencing MaxHi varied between Tempc LHS sets (Tables 4 and S2). At the warmest Tempc, increased mean annual temperature and later day of introduction both decreased MaxHi, but this was not observed in the other sets. Aedes albopictus traits were important at the coldest Tempc, but became less important as Tempc increased. Aedes aegypti traits were important at all Tempc values. Unlike the likelihood of epidemics, mosquito mortality parameters did not strongly affect MaxHi. Peak activity occurred faster (lower lag between tcrit and tMaxHi) at the highest Tempc, but again

Size and Lag Time of Epidemics
One measure of epidemic activity is the number of human cases at the peak of the epidemic (MaxH i ). This is influenced by the model behavior and duration of virus activity, as some simulations showed single, short outbreaks ( Figure 3) while others had prolonged, lower viral activity. The average MaxHi increased with Temp c , but the range was large in all three sets (Table 2, Figure S3). Adjusted R 2 values again were low, indicating that interactions are likely to be important (Table 2). Parameters most influencing MaxH i varied between Temp c LHS sets (Table 4 and Table S2). At the warmest Temp c , increased mean annual temperature and later day of introduction both decreased MaxH i , but this was not observed in the other sets. Aedes albopictus traits were important at the coldest Temp c , but became less important as Temp c increased. Aedes aegypti traits were important at all Temp c values. Unlike the likelihood of epidemics, mosquito mortality parameters did not strongly affect MaxH i . Peak activity occurred faster (lower lag between t crit and tMaxH i ) at the highest Temp c , but again variation was high ( Table 2). The most influential parameters varied between Temp c LHC sets, but again the main effects regression models had low predictive power (Tables 2 and 5). The effect of mean annual temperature varied, decreasing lag at low Temp c but increasing lag at higher Temp c . At the higher Temp c values, the population size of both species was important and had a negative effect on lag, while at the coldest temperatures, transmission parameters had more effect (Table 5 and Table S3). Central Temperature, Temp c ( • C) 10 16 22 Annual mean temperature (T mean ) (-) Day of introduction (t crit ) (-) Transmission: humans to Ae.

Comparisons of Parameter Sets Across Temp c
Comparing individual parameter sets across Temp c values, the general pattern was increasing likelihood of epidemics as Temp c increases (Table S4). In 60% of parameter sets the outcome was not epidemic at Temp c = 10, then was epidemic only at Temp c = 16 (12%), only at Temp c = 22 (32%), or both (16.4%). We did not observe any parameter sets epidemic at the coldest and warmest temperatures but not epidemic with the middle Temp c temperature. Decreasing likelihood of epidemics with increasing Temp c occurred rarely, with 4.4% of parameter sets epidemic only at Temp c = 10 and 2.8% epidemic at Temp c = 10 and 16 but not 22.
The peak size of the outbreak also typically increased with Temp c . Of 57 parameter sets epidemic at all three values, the size increased as Temp c increased in 84% (Table S5). With parameter sets showing epidemics at only two values of Temp c , the peak was higher at the warmer Temp c in 77% of the sets.

Comparison of Strain-Species Combinations
Baseline model runs with species characteristics set identically (functionally one large population) showed identical outbreak sizes regardless of the vector competence assigned to each population and the relative population sizes ( Figure 6). When species-specific parameters are used, an aegypti-adapted strain shows the same outbreak sizes with only Ae. aegypti populations and with Ae. aegypti-dominated populations. Equal population sizes of both vectors reduces outbreak size slightly. However, with an albopictus-adapted strain (in transmission parameters only), a different pattern is observed. Highest outbreak sizes are lower than those observed with identical species characteristics or with aegypti-adapted virus strains, and increasing the proportion of Ae. albopictus in the population further decreases outbreak size. increasing Tempc occurred rarely, with 4.4% of parameter sets epidemic only at Tempc = 10 and 2.8% epidemic at Tempc = 10 and 16 but not 22. The peak size of the outbreak also typically increased with Tempc. Of 57 parameter sets epidemic at all three values, the size increased as Tempc increased in 84% (Table S5). With parameter sets showing epidemics at only two values of Tempc, the peak was higher at the warmer Tempc in 77% of the sets.

Comparison of Strain-Species Combinations
Baseline model runs with species characteristics set identically (functionally one large population) showed identical outbreak sizes regardless of the vector competence assigned to each population and the relative population sizes ( Figure 6). When species-specific parameters are used, an aegypti-adapted strain shows the same outbreak sizes with only Ae. aegypti populations and with Ae. aegypti-dominated populations. Equal population sizes of both vectors reduces outbreak size slightly. However, with an albopictus-adapted strain (in transmission parameters only), a different pattern is observed. Highest outbreak sizes are lower than those observed with identical species characteristics or with aegypti-adapted virus strains, and increasing the proportion of Ae. albopictus in the population further decreases outbreak size.

Discussion
Sensitivity analysis of model outcomes can include systematic variation of individual parameters or randomized methods such as the LHC method. Here, we show that systematic variation of one individual parameter could change the relative ranking of other parameters in an LHC analysis, in addition to changing model behavior. Although it is generally not possible to vary all parameters in a model analysis, this highlights the need to consider the impact of fixed parameters and to identify the values used in any analysis. Even if chosen by fitting to a particular data set, fixed

Discussion
Sensitivity analysis of model outcomes can include systematic variation of individual parameters or randomized methods such as the LHC method. Here, we show that systematic variation of one individual parameter could change the relative ranking of other parameters in an LHC analysis, in addition to changing model behavior. Although it is generally not possible to vary all parameters in a model analysis, this highlights the need to consider the impact of fixed parameters and to identify the values used in any analysis. Even if chosen by fitting to a particular data set, fixed parameters may vary in a broader context so model results should be interpreted with caution outside the conditions set for a specific analysis.
Here, we showed changes in the relative importance of other parameters with variation in one parameter, the center of the temperature range for minimal mosquito mortality. Notably, when this parameter was higher, the slope of the temperature-mortality relationship for Ae. aegypti was identified as an important factor in likelihood of epidemics. Interestingly, the same parameter for Ae. albopictus did not have a similar effect, and it did not have strong effects on the size of an outbreak for either species. We used a broad range for Temp c , likely extending beyond biological realism for Ae. aegypti. Aedes aegypti is a more tropical mosquito and mortality is likely to be high at lower temperatures. Extending the upper range of this parameter would be of interest, given the increase in the likelihood of invasion at the higher value and the relationship with the mean annual temperature. Future work will consider more specific mortality relationships and broader ranges.
Adult mosquito mortality has been identified as an important parameter in other arbovirus models, even with constant mortality rates or simple relationships with temperature. Complex relationships between mortality and temperature have been observed for Ae. aegypti, but less is known about Ae. albopictus. Clearly a better understanding of the factors influencing mortality, and improving inclusion of these complexities in arbovirus transmission models, is needed.
Our results support the general consensus that Ae. aegypti is more important as a vector in CHIKV dynamics in many environments, but Ae. albopictus can support transmission and outbreaks either alone or in mixed populations with low Ae. aegypti abundance. Interestingly, even with a viral strain adapted to Ae. albopictus, Ae. aegypti supported higher transmission ( Figure 6). However, this analysis did not include changes to the EIP known to occur with some mutations, and this should be considered further. In our analysis, transmission rates (b and β) were consistently identified as key factors, while virus development rates (γ) had limited impacts on outcome variables, suggesting that mutations affecting transmission rates may be more important. The virus may acquire adaptive mutations during an outbreak, and including variation in key parameters over time may be important in understanding how this would affect virus dynamics. In this model, the higher transmission with Ae. aegypti likely depended on the species-specific host preferences (higher biting rates on humans). The more opportunistic feeding patterns of Ae. albopictus have been reported from multiple locations, but there may also be situations where human blood meals are more frequent. In the field, host preference may be dynamic and realized feeding dependent on host availability, host preference [37,[58][59][60], and the degree to which both species exhibit gonotrophic discordance (e.g., [53,61]). Additional exploration of the model including variation in human biting rates over time would be of interest to better understand the potential role of Ae. albopictus in CHIKV local transmission.
Other models of CHIKV have shown similar results in some aspects, while diverging in others. Manore et al. [20] developed a model to compare CHIKV and DENV transmission with either Ae. aegypti or Ae. albopictus vectors, but did not include mixed mosquito populations, seasonal mosquito dynamics, or temperature dynamics. Some of the same parameters were identified as having the most influence on the outcome, including mosquito mortality, biting rate, and transmission rates. Overall their model showed much slower dynamics for CHIKV epidemic development after introduction than did ours. This is likely due to the temperature dependence in our model, leading to seasonally changing mosquito mortality rates and duration of the latency stage. Our central transmission rates were lower than the baseline transmission rates of [20] as well, likely slowing the dynamics. Adjusting parameters in our model to approximate those in the Manore et al. model led to similarly slow dynamics (not shown).
Our results support temperature effects on the likelihood, size, and structure of epidemics. However, successful invasions occurred over the full range of mean annual temperatures for Florida, and at all three values of Temp c (Figure 4). While the interaction between the annual mean temperature and the central temperature for minimum mosquito mortality was complex, the range of conditions with successful invasions implies that all of Florida is at risk of invasion by CHIKV. Other areas that have similar climates and extend over temperate-subtropical zones, such as Brazil, may have similar effects of mean annual temperatures and temperature-dependent parameters. Outbreaks in other temperate areas (e.g., Italy) further illustrate the complex interactions between temperature and transmission dynamics. Extensions of the sensitivity analysis to higher and lower T mean values would be of interest to explore the risk in other areas. This would, however, require assessment of the relationship between mean and amplitude of the temperature curve, as the relationship used here will not always be appropriate (Salami et al. unpublished). Several input parameters influenced the likelihood of epidemics, but over the ranges considered none showed threshold effects, suggesting that there is a wide range of parameter space that allows outbreaks.
The parameters identified as important in this analysis are also likely to influence other Aedes-transmitted viruses such as ZIKV and DENV, suggesting areas for research to predict the outcomes of introductions of these viruses. However, virus-specific parameter estimates, along with location-specific temperature curves, are needed to assess the likelihood of invasion and identify factors that most influence epidemic dynamics.
Adult recruitment patterns were based on rainfall and modeled similarly for both vector species. However, in reality egg diapause in Ae. albopictus reduces hatching probabilities of this species in the late fall and early winter periods in Florida [50]. By contrast, egg diapause does not occur in Ae. aegypti, which does not experience seasonal hatching interruptions. These recruitment differences could also contribute to the relative involvement of the two species in transmission dynamics. Different seasonal patterns of recruitment affected simulated dynamics of WNV [19] and should be considered in more detail for CHIKV.
With repeated introductions of arboviruses into Florida and elsewhere, we need to improve our understanding of what influences whether an introduction leads to establishment and outbreaks. Mosquito mortality has consistently been identified as important, is expected to vary considerably in the field, and is a target for novel control efforts to mitigate pathogen transmission (e.g., life-shortening microbes, Wolbachia, [62], Beauveria bassiana, [63]; gene drive systems, [64]). Further exploration of this is needed, including temperature dependency, species and individual variation, and other environmental influences. Stochastic factors in mortality and human-mosquito contact will also play a role in virus establishment, and individual based models allowing this structure would be valuable. Expanded consideration of virus lineages, genotypes, and mutations and interactions with mosquito species will be valuable to better understand the potential for virus invasion and outbreaks. Interactions between parameters are clearly important in model analysis, but it is unclear how this relates to field conditions. Additionally, spatial variation in human-related factors may influence model parameters (e.g., socioeconomic factors [65]). Incorporating interactions and increased variation into modeling studies will improve our ability to understand arbovirus dynamics and target empirical studies on critical factors.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1999-4915/12/8/830/s1, Figure S1: Relationship between mean and range of yearly temperatures for Florida weather stations. Figure S2: Phases of modeled mosquito populations. Figure S3: MaxH i histograms for each Temp c . Figure S4: Relationship between epidemic occurrence and selected parameters. Table S1, Results of logistic regression on epidemic behavior; Table S2, Results of regression on MaxH i ; Table S3, Results of regression on lag; Table S4, Increasing epidemics with Temp c ; Table S5, Increasing size of epidemics with Temp c ; Table S6