San Juan is a municipality of Puerto Rico and is located in the northeastern part of the Caribbean island. The population was 395, 326 in 2010 (US Census Bureau, 2010 Census). San Juan has a humid tropical climate with minimal variation in seasonal temperature. Precipitation occurs all year, but it is notably drier during boreal winter and early spring. San Juan is an ideal study location because dengue is endemic in the municipality, and weekly dengue case data are available through CDC ArboNet (USGS, 2014) from 2010–2013. During that time, annual case numbers ranged from 500 (2011) to 919 (2012). These data are crucial for model evaluation.

The mosquito population was simulated using the general structure of the Dynamic Mosquito Simulation Model, (DyMSiM) [8] but the model was parameterized for Ae. aegypti mosquitoes with additional components added, including the epidemiological SEIR model (S2 Table). This enabled the simulation of virus transmission between the human and mosquito populations in DyMSiM. The model is deterministic and was implemented using Euler’s method of integration in Stella 10.06 software. A conceptual diagram of model processes is provided in Fig 1. Parameter value constants and equations are provided in S1 Table, governing equations for human and mosquito populations are provided in S2 Table, and further details are discussed below.

Because it is difficult to obtain precise measurements for many of the DyMSiM parameters, especially those related to immature container habitat, we used a Monte Carlo approach to assess outcomes for a range of possible parameter values. We performed 9600 simulations with DyMSiM using discrete parameter values representative of the distribution of values identified in the literature and the preliminary model simulations (Bold in Column 3 of S1 Table) to replicate DENV transmission in San Juan, Puerto Rico from 2010–2013. During each simulation, one value was changed until all parameter value combinations were used. All humans and mosquitoes were assumed to be susceptible at the beginning of the runs, however, runs were started in 2009 to provide spin-up time to build infection and immunity (human only) within the human and mosquito population. Additionally, a background infection rate was used to initiate infections and is based off of the average number of dengue cases that occur during the low transmission season according to the 2010–2013 case data for San Juan. The assumption that the population is completely susceptible is not accurate, however, given the low incidence (3192 total cases over a population of 395,326) and short time interval over which the model is run, the influence on outcomes should be minor. Weekly CDC reported case data were used to evaluate the model.

Container habitat area estimates were determined by first performing preliminary model simulations under parameter values that would produce a maximum and minimum number of dengue cases in San Juan. By comparing these runs with the reported dengue data, an estimated range of habitat area was determined for the study. Because there is great uncertainty in the container habitat area and it strongly influences mosquito population size, the number and range of values chosen greatly exceeds that of the other parameters. Fewer values and smaller ranges were chosen for parameters where more research has been performed resulting in greater confidence in the selected values.

There were 3,192 cases of DF reported during the study period with 901 cases occurring in 2010, 500 cases in 2011, 919 cases in 2012, and 872 cases in 2013. Reported cases began increasing during week 20 in 2010, 2012, and 2013 but occurred later (week 25) in 2011 which was relatively cool compared to the other years and experienced the fewest number of cases (Fig 2). During most years (2010, 2011, and 2013) the peak in reported DF cases occurred between weeks 32 and 37, however, during 2012 the increase in cases was much more gradual and the peak did not occur until week 49 and did not decline to minimum season case numbers until week 15 of 2013. Temperatures were also much warmer during the second half of 2012 compared to other years, while precipitation was notably lower, especially during summer and fall. Although 2013 experienced similar reported case numbers as 2010 and 2012, most cases occurred during the first 15 weeks of the year (Fig 2) with fewer cases occurring during its summer peak. The decline to reported minimum season case numbers also occurred earlier than the other years. Aside from 2012, increases in reported DF followed rains occurring in late spring. The earliest increase and peak in reported DF cases occurred in 2010 which also experienced the warmest spring temperatures and earliest precipitation.

Statistical measures of accuracy between the reported data and ensemble averages using the top 1% (n = 96) best-fit simulations, for both the entire time period (S2 Dataset) and individual years (S3–S6 Datasets), are reported in Table 1. Although inter-annual variability in DF cases was well-simulated in the ensemble average from the 96 model simulations over the entire time period (r² = 0.44, p = 0.002), the annual epidemic curves for individual years often lacked precision (Fig 2). When evaluating ensemble averages from the simulations over individual years (same simulations but the statistics are only performed on reported case data for individual years), however, intra-annual variability was replicated with much higher accuracy for 2010 (r² = 0.90, p = 0.03), 2011 (r² = 0.83, p = 0.05), and 2012 (r² = 0.94, p = 0.01) (Fig 3). Simulations for 2013 were markedly less accurate (r² = 0.25, p = 0.25) compared to the other years (Fig 3). By contrast, the yearly model accuracy assuming no a-priori knowledge to guide parameter selection was lower and more variable: the r² values for 2010, 2011, 2012, and 2013 were 0.77, 0.64, 0.05, and 0.12 respectively.

Fig 4 displays the number of times each parameter value was used in the top 1% of simulations. Notably, in 2010 and 2011 the number of simulations retained that used a high proportion of open containers (i.e., more precipitation dependent sources) vs. a low proportion of open containers (i.e., more human-managed sources) is much larger than in 2012 and 2013 (76 and 59 vs. 0 and 0). Additionally, patterns in container habitat area, background infection rate, and carrying capacity vary considerably between years. The retained simulations for most years used the higher adult daily survival rate, and the shorter length of infectious period (except 2012). With the exception of 2011, the range of values tested for container height and host infection probability are used with equal frequency between years indicating that they did not have a strong impact on temporal variation in dengue transmission dynamics.

Fig 5 displays the average rank of simulations using each parameter value among all 9600 simulations. Simulations that used a high proportion of open containers ranked slightly better for 2011, while the opposite is true for 2012 and all years. This is consistent with the trend in the top 1% of simulations, however, there is no strong difference for 2010 and 2013. Simulations using higher values for container habitat area were markedly worse, though this effect was lesser for 2012. Additionally, the lower carrying capacity value tended to produce better simulations. Average ranks among different values did not vary markedly for the other parameters.

Precipitation was highest in 2010 (227.5 cm) and 2011 (224.0 cm), with the onset of the rainy season starting a few weeks earlier in 2010 (Fig 6). Precipitation was by far the lowest in 2012 (140.3 cm), however, temperatures were warmest in 2012, especially during late fall (~1°C warmer in November). Precipitation was 216.3 cm in 2013 –only slightly drier than 2010 and 2011—but the timing of the precipitation was different, exhibiting a marked drop off during late summer/early fall.

While we used the best available information from the literature to specify values of the DyMSiM parameters, some values are based on the results of only a few studies or were reasonable estimates when quantitative results were unavailable. Papers that consolidated and considered data from multiple studies proved especially useful because they accounted for variations arising from conducting research on different colonies of mosquitoes and strains of virus or using different methods to answer the same question [17,39]. Often averages were used to enhance simplicity and interpretability of the model; yet, many parameters may have values across an unknown distribution surrounding their mean due to genetic and environmental variability. For example, genetic variations have been found to alter temperature regulation of West Nile virus transmission [63]. If a similar variation in temperature sensitivity exists between or within dengue serotypes it could be an important component to the severity of DF epidemics. We attempted to mitigate our limited knowledge of some parameter values by using a Monte Carlo approach when they were particularly difficult to ascertain. The model showed that some parameters, such as container habitat area and composition of containers, were particularly important for simulating transmission. The carrying capacity density was also important as it influenced the mosquito productivity within containers, and was directly related to their area. For example, simulations were degraded when using larger container areas, but this effect was lessened when using the lower maximum larval density value. The value of the background infection rate was also important when selecting the best fit models. This may indicate the importance of varying virus introduction or virulence between years. Vertical transmission of the virus was not considered in this model though it has been observed in the laboratory and field setting and could help maintain the virus during unfavorable conditions [64,65]. However, Adams et al. [66] found it did not have considerable importance in their model and the frequency of its occurrence in nature is still uncertain [67]. As more evidence is generated establishing links among meteorological factors, vectors, virus, and humans, parameter values can be better defined.

Some of the inconsistencies between modeled and reported case data may also be attributed to the nature of surveillance data. Because the surveillance system is passive, the reported case data are only an approximation of total DF incidence for the municipality. Incomplete or inaccurate reporting can arise due to sub-clinical infections, misdiagnosis, failure to report cases, differences between where virus transmission is reported and where it was acquired, and variability in the time between transmission of the virus, development of symptoms, and clinical visits. Relative increases in dengue reporting may be seen if there is increased media attention, circulation of a more severe strain, or increases in laboratory testing ability, etc. The serotype of the virus is often not reported but could be of particular importance if a high proportion of the population is immune to the circulating serotype. It is likely that not all cases are being reported to the CDC and, therefore, the reported case data is only a subset of the actual number of cases occurring in San Juan though it is difficult to estimate the amount of underreporting. The influence of under-reporting should be minimal if the patterns of reported cases represent the yearly variability even if the overall magnitude of reported cases is inaccurate. The largest bias in the calculation of the correlation between modeled output and reported case numbers would occur if there were significant changes to the reporting system or if there are interannual differences in reporting. There were no significant changes to the surveillance system protocol during this time period.

It is difficult to determine inter-annual variability in reporting, but there is the possibility for bias. For example, if individuals at the beginning of the season were more likely to seek care and be reported through the surveillance system, or if physicians were less likely to diagnose and report cases during typically lower periods of dengue, the shape of the curve might shift. Given the uncertainty of these types of scenarios, trying to determine the relative impact they may have on the correlation between model output and reported case numbers would be highly speculative. We acknowledge, however, that under-reporting of dengue is an important issue and therefore, the model output should not be interpreted as an estimate of actual dengue case numbers but as a representation of inter-annual variability of DENV transmission. If more reliable estimates of reporting bias and information on the temporal variability in asymptomatic cases or the reporting of symptomatic cases become available, it would be useful to compare numbers of modeled and reported dengue case data. Despite these limitations, the reported case data are currently the best metric available for tracking DENV transmission and have been used in other studies to represent levels of DENV transmission in Puerto Rico [4,61,68].

The results of this study provide important information for dengue control. First, epidemics can occur during both wet and dry years. Dengue cases were comparatively high in 2010, one of the wettest years on record for San Juan, and in 2012, a drier-than-normal year in which the rains came later than average. The warmer winter and spring temperatures are likely responsible for the greater number of DF cases during 2010 while warmer fall and winter conditions in 2012 likely helped propagate dengue transmission later in the year. It is common for dry years to be hotter due to a combination of factors that include: 1) fewer clouds facilitating higher daytime temperatures and 2) less evapotranspiration freeing a greater amount of solar radiation for heating. Second, the timing of the dengue cases varies substantially from year-to-year and appears to be linked with the timing and magnitude of rainfall. In the wetter 2010 year, DF case numbers rose sharply after the onset of the rains (Week 10) and peaked near Week 30, whereas in the drier 2012 year the cases began later (Week 20) after a highly sporadic beginning to the rainy season, and gradually climbed to a much later peak near Week 48. In both instances, the onset of DF cases coincides with or begins shortly after the onset of the rainy season, and rises steadily due to the time required for the mosquito population to development in newly formed habitats and completion of the EIP. Lastly, despite a tropical climate, small fluctuations in temperature can have considerable effects on DENV transmission. Often tropical environments are thought to be warm enough for year-round transmission with precipitation being the limiting factor. Our results indicate that nonlinear relationships result in greater temperature impacts on DENV transmission than would be expected. Monitoring climate as a proxy of dengue risk may provide public health workers with a simple tool to prepare and execute transmission intervention methods and arrange testing and treatment protocols before the DF season peaks. Early warning systems can be developed based on the identification of rainfall and temperature patterns that promote rapid dengue transmission [69–72].

Dominant Ae. aegypti container habitat appeared to vary with precipitation patterns. During the wetter years, 2010 and 2011, simulations using a smaller container area and more rain dependent containers were ranked better; however, during 2012 simulations using a greater amount of human managed water sources ranked better. This suggests that Ae. aegypti inhabit rain filled containers when available, in addition to human managed water sources, but drier conditions may limit immature mosquito habitats to permanent (septic tanks, large water tanks, etc.) or human managed (animal drinking pans, plant trivets) water sources. This pattern has been observed in field studies [58]. From a public health perspective, the magnitude of precipitation in a given year may inform how to refine the message conveyed to the public about how to limit mosquito habitat. During wet years, mosquito control strategies should focus on eliminating old tires, buckets, and other items that can be filled via precipitation. Drier years should shift the focus to treating permanent or human managed water sources and covering stored water to prevent mosquito infestation. Further investigation of how weather and rainfall patterns influence human behaviors such as water storage are needed to determine if this could be incorporated into the model.

Understanding weather influences on DENV transmission is important in the context of climate change especially as DENV range expands in the Americas. Outbreaks in Hawaii, along the Texas/Mexico border, and in southern Florida indicate a reemergence of dengue in the U.S. [73–75]. Additionally, sero-prevalence studies indicate that dengue cases are likely underreported in the U.S. [73,75] due to subclinical infections, insufficient resources for testing, lack of sensitization of the medical community and patients seeking medical attention across the border in Mexico [75]. Although infrastructure differences, such as the prevalence of air conditioning, were shown to reduce transmission on the US side of the Mexico border in Laredo, autochthonous cases still occur [18]. Viral introduction in response to frequent travel across the border by residents may be a large source of infections [76] and weather may still have an important influence on transmission risk [20]. Whether DENV is endemic to an area or currently only poses a threat, understanding the influence of weather and climate on DENV ecology can facilitate strategies to prevent or mitigate transmission.

Using known relationships between climate variables, Ae. aegypti dynamics, DENV replication and transmission, and a basic SEIR model within a modeling framework known as DyMSiM we were able to simulate inter-annual variability in DF cases in San Juan County, PR for the years 2010–2013. When optimizing the DyMSiM parameters for a single year, we were able to simulate intra-annual variability well for three of the four years. Meteorological factors were important influences on dengue transmission for 2010, 2011, and 2012 but not as important in 2013. Parameter values changed in response to meteorological conditions, illustrating the complexity of DENV ecology.

The results of this study quantify meteorological impacts on DENV transmission and demonstrate the variance in parameter values that result from changing environmental conditions. Rain patterns were especially important for determining the timing of epidemics and the primary habitat for immature Ae. aegypti vectors. Abnormally high temperatures have the potential to extend the transmission season (even during dry years) as exemplified in 2012. In general, temperature had the greatest influence on annual case numbers. The sensitivity of DENV ecology to meteorological variables and their interactions underscores the utility of process-based modeling when studying the impacts of climate variability and climate change on vector-borne disease. Information from this study can be used by public health officials to improve DENV transmission intervention strategies in San Juan, PR by way of advanced preparedness through model predictions and creation of more targeted vector control campaigns. As our understanding of the ecology of the virus and its thresholds improves, so will our ability to implement effective mitigation strategies against DENV transmission and reduce the disease burden on human populations.

Future research endeavors will seek to enhance model performance and utility. Validation of DyMSiM in other areas with longer term datasets would be useful in determining its ability to simulate DENV transmission in other environments and its ability to simulate DF incidence over long time periods. The deficit of long, reliable dengue surveillance records has forced training of dengue models over periods of less than 10 years and validation on one or only a few years [70,77,78]. Inclusion of serotypes would also enhance the model given that immunity to specific virus serotypes can affect transmission dynamics by limiting the susceptible population. The current SEIR model is simple but can still be used to examine the dynamics of dengue over short time periods. Incorporating surveillance of multiple virus serotypes, including short term cross-immunity, will be a priority in future research given the importance of immunological status on virus transmission dynamics. Additionally, the inclusion of social factors that influence population risk, such as migration, age structure, personal and household level prevention behaviors and other factors influencing vector-human contact, may provide further insights into the nature of DF epidemics. Finally, future work will concentrate on the possible impacts of climate change on dengue transmission. For instance, under climate change, winter temperatures may become suitable for completion of the EIP allowing DENV transmission to continue. Understanding the effects of future weather and climate conditions on DENV transmission is vital for limiting increases in incidence of disease and range expansion.