==== Front PLoS Negl Trop Dis PLoS Negl Trop Dis plos plosntds PLoS Neglected Tropical Diseases 1935-2727 1935-2735 Public Library of Science San Francisco, CA USA 33226987 PNTD-D-19-00900 10.1371/journal.pntd.0008868 Research Article Biology and Life Sciences Population Biology Population Metrics Death Rates Medicine and Health Sciences Medical Conditions Infectious Diseases Disease Vectors Insect Vectors Mosquitoes Biology and Life Sciences Species Interactions Disease Vectors Insect Vectors Mosquitoes Biology and Life Sciences Zoology Entomology Insects Mosquitoes Biology and Life Sciences Organisms Eukaryota Animals Invertebrates Arthropoda Insects Mosquitoes Biology and Life Sciences Zoology Animals Invertebrates Arthropoda Insects Mosquitoes Medicine and Health Sciences Medical Conditions Tropical Diseases Neglected Tropical Diseases Dengue Fever Medicine and Health Sciences Medical Conditions Infectious Diseases Viral Diseases Dengue Fever Medicine and Health Sciences Medical Conditions Infectious Diseases Infectious Disease Control Physical Sciences Mathematics Differential Equations Biology and Life Sciences Immunology Immunity Medicine and Health Sciences Immunology Immunity Medicine and Health Sciences Epidemiology Engineering and Technology Signal Processing Autocorrelation Research and Analysis Methods Mathematical and Statistical Techniques Statistical Methods Autocorrelation Physical Sciences Mathematics Statistics Statistical Methods Autocorrelation Linking mosquito surveillance to dengue fever through Bayesian mechanistic modeling Mosquito surveillance and dengue riskhttps://orcid.org/0000-0003-4280-0375Leach Clinton B. ConceptualizationFormal analysisMethodologySoftwareVisualizationWriting – original draftWriting – review & editing12* https://orcid.org/0000-0003-0414-0193Hoeting Jennifer A. ConceptualizationMethodologySupervisionWriting – review & editing2 Pepin Kim M. ConceptualizationWriting – review & editing3 Eiras Alvaro E. Data curationResourcesWriting – review & editing4 Hooten Mevin B. ConceptualizationMethodologyWriting – review & editing256 Webb Colleen T. ConceptualizationMethodologySupervisionWriting – review & editing1 1 Graduate Degree Program in Ecology, Colorado State University, Fort Collins, Colorado, United States of America 2 Department of Statistics, Colorado State University, Fort Collins, Colorado, United States of America 3 National Wildlife Research Center, United States Department of Agriculture, Wildlife Services, Fort Collins, Colorado, United States of America 4 Departamento de Parasitologia, Universidade Federal de Minas Gerais, Belo Horizonte, Minas Gerais, Brazil 5 U.S. Geological Survey, Colorado Cooperative Fish and Wildlife Research Unit, Fort Collins, Colorado, United States of America 6 Department of Fish, Wildlife, and Conservation Biology, Colorado State University, Fort Collins, Colorado, United States of America Althouse Benjamin Muir Editor Institute for Disease Modeling, UNITED STATES The authors have declared that no competing interests exist. * E-mail: clint.leach@gmail.com 11 2020 23 11 2020 14 11 e000886828 5 2019 8 10 2020 This is an open access article, free of all copyright, and may be freely reproduced, distributed, transmitted, modified, built upon, or otherwise used by anyone for any lawful purpose. The work is made available under the Creative Commons CC0 public domain dedication.Our ability to effectively prevent the transmission of the dengue virus through targeted control of its vector, Aedes aegypti, depends critically on our understanding of the link between mosquito abundance and human disease risk. Mosquito and clinical surveillance data are widely collected, but linking them requires a modeling framework that accounts for the complex non-linear mechanisms involved in transmission. Most critical are the bottleneck in transmission imposed by mosquito lifespan relative to the virus’ extrinsic incubation period, and the dynamics of human immunity. We developed a differential equation model of dengue transmission and embedded it in a Bayesian hierarchical framework that allowed us to estimate latent time series of mosquito demographic rates from mosquito trap counts and dengue case reports from the city of Vitória, Brazil. We used the fitted model to explore how the timing of a pulse of adult mosquito control influences its effect on the human disease burden in the following year. We found that control was generally more effective when implemented in periods of relatively low mosquito mortality (when mosquito abundance was also generally low). In particular, control implemented in early September (week 34 of the year) produced the largest reduction in predicted human case reports over the following year. This highlights the potential long-term utility of broad, off-peak-season mosquito control in addition to existing, locally targeted within-season efforts. Further, uncertainty in the effectiveness of control interventions was driven largely by posterior variation in the average mosquito mortality rate (closely tied to total mosquito abundance) with lower mosquito mortality generating systems more vulnerable to control. Broadly, these correlations suggest that mosquito control is most effective in situations in which transmission is already limited by mosquito abundance. Author summary The contribution of the mosquito vector Aedes aegypti to the spread of dengue fever depends not only on their abundance, but also on the likelihood of an exposed mosquito living long enough to incubate the dengue virus and subsequently transmit it to a susceptible human host. We developed a mechanistic model that accounts for the role of this process in the dynamics of dengue fever and fit the model to a time series of human case reports and mosquito trap counts from the city of Vitória, Brazil. We then used this fitted model to simulate the effect of mosquito control implemented at different times of the year and found that mosquito control leads to the largest reduction in human dengue cases over the following year when implemented in early September, during the dengue off-season. Further, the effectiveness of mosquito control was strongly negatively correlated with the overall average abundance of mosquitoes. Together with the timing of effective control, these results suggest that mosquito control is most effective when mosquitoes are already limiting to transmission. http://dx.doi.org/10.13039/100000001National Science FoundationDGE-1321845https://orcid.org/0000-0003-4280-0375Leach Clinton B. http://dx.doi.org/10.13039/501100003593Conselho Nacional de Desenvolvimento Científico e Tecnológico312335/2018-0Eiras Alvaro E. CBL was supported by the NSF Graduate Research Fellowship Program under grant DGE-1321845. AEE was supported by the Brazilian National Council for Scientific and Technological Development (CNPq grant 312335/2018-0). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. PLOS Publication Stagevor-update-to-uncorrected-proofPublication Update2020-12-07Data AvailabilityAll data files are available from Figshare at https://doi.org/10.6084/m9.figshare.7905254.Data Availability All data files are available from Figshare at https://doi.org/10.6084/m9.figshare.7905254. ==== Body Introduction Dengue fever is a massive global public health burden, with millions of cases per year [1]. Because the dengue virus (DENV) is transmitted by the mosquito Aedes aegypti, dengue fever is prevented primarily through mosquito control programs [2]. Though there have been documented successes, there is limited evidence for the long-term sustainability and effectiveness of these control programs [3]. As a result, there is a growing recognition that effective control needs to be guided by high quality vector surveillance, together with quantitative tools that synthesize vector surveillance with clinical surveillance, account for local epidemiology, and facilitate local decision making [3, 4]. Moreover, mosquito control needs to be guided by an understanding of the link between mosquito abundance and disease risk so that the mosquitoes most responsible for transmission can be targeted [4, 5]. Many of the attempts to establish this link have found a weak relationship between mosquito abundance indices and incidence of disease in humans [6–8]. However, these attempts often do not account for the complex, non-linear interactions that mediate the relationship between mosquito abundance and human disease. In particular, host immunity is a key intrinsic driver of infectious disease dynamics, and conditions favorable for transmission can only lead to an outbreak of disease when there is a sufficiently large population of susceptible hosts [9, 10]. As such, the ability of mosquitoes to contribute to DENV transmission depends critically on the level of immunity in the human population [5]. Further, the cycle of transmission between humans and mosquitoes is influenced not just by mosquito abundance, but also by mosquito survival relative to the virus incubation period in mosquitoes [11]. In fact, whether or not an exposed mosquito will survive long enough to become infectious represents a critical bottleneck in the transmission process and leads to nonlinear dependence of transmission on mosquito survival [11]. The importance of intrinsic nonlinearities, potentially alongside seasonality and stochastic forcing [9, 12, 13], in governing human disease risk highlights the need to integrate mechanistic modeling into the quantitative tools used to understand the effects of control interventions. Such mechanistic models can often perform better than complex autoregressive statistical models in describing and forecasting population dynamics [14]. Moreover, in the absence of case-control studies, mechanistic models can provide scenario-based tools that can be used to predict the effect of management actions [15]. Differential equation models provide a natural way to describe mechanistic processes, but that description must also account for sources of uncertainty [16, 17]. In particular, the values of parameters (e.g., the average length of time for which a host is infectious) are often uncertain, which can lead to large uncertainty about the effects of management actions [18]. The structure of the processes themselves can be uncertain [12], and needs to be informed by available, often noisy, data. Bayesian hierarchical modeling provides a coherent framework to account for and integrate this uncertainty across the three levels of the model (data, process, and parameters [19, 20]). In what follows, we integrate these elements—a detailed mechanistic model of dengue transmission with a full Bayesian accounting of uncertainty—to better understand the interplay of forces governing dengue dynamics and their interaction with potential vector control interventions. We apply this framework to clinical and entomological surveillance data from the city of Vitória, Brazil. These data allow us to estimate a latent time series of mosquito mortality rates that modulate the transmission process and link mosquito abundance to human disease. We then use the fitted model to explore how perturbations to the mosquito population propagate and interact with the nonlinearities of dengue transmission to better inform mosquito control efforts. Methods Ethics statement We did not obtain Institutional Review Board (IRB) approval for this work as the data we worked with were received by us as aggregated data at the weekly and neighborhood level. Hence, this research does not meet the definition of human subjects research requiring IRB approval. The data were analyzed in the aggregated form, which protects the anonymity of individuals. Study system and data Vitória is a coastal city and the capital of the state of Espírito Santo, Brazil, with a population of 327,801 as of 2010 [21]. Since 2008, the company Ecovec has monitored mosquito abundance for the city using approximately 1327 sticky traps (MosquiTRAP, [22]) arranged in a roughly 250m grid across the city [7, 23]. Each trap is checked weekly and the mosquitoes inside counted and identified, with the results sent to a central database that city managers then use to map mosquito infestations and target control. These data comprise 243 weeks (week 1 of 2008 through week 34 of 2012) of total city-wide counts of trapped gravid female Aedes aegypti. It is important to note that this time series reflects both natural fluctuations in mosquito density and fluctuations driven by the city’s existing mosquito control program. In addition, dengue fever is a mandatory notifiable disease, and thus the city’s Ministry of Health Secretary maintains a database of weekly notified probable dengue cases (i.e., medical care sought for dengue-like symptoms) for the same time period. Process model Dengue epidemiology is complicated considerably by the presence of four simultaneously circulating serotypes. Infection with one serotype confers life-long immunity to that serotype, along with temporary immunity to other serotypes [24]. As this cross-immunity wanes, antibodies from the previous infection can result in antibody-dependent enhancement (ADE), wherein human hosts are more susceptible to infection with the other serotypes and more likely to develop severe symptoms (i.e., dengue hemorrhagic fever or dengue shock syndrome) [24]. The strength and duration of these different inter-serotype interactions are not well understood, although different models suggest that temporary cross-immunity alone (without ADE) is sufficient to reproduce observed multi-annual dynamics in Thailand [24, 25]. Explicitly capturing the cross-immune interactions among all four serotypes, or even only two of the four [26], leads to a large and complex mechanistic model. Moreover, because dengue case reports do not identify serotype, there is not enough information in our data to inform the dynamics of individual serotypes. As such, we captured temporary cross-immunity, and the potential for multiple sequential infections, as simply and tractably as possible in a susceptible–exposed–infectious–recovered–susceptible (SEIRS) compartment model, similar to [27–30]. Although it captures the critical influence of temporary immunity, this framework does not account for the potential relationship between an individual’s infection history and the likelihood that a new infection will be symptomatic (and thus reported). In particular, secondary infections appear more likely to be symptomatic than primary infections [31, 32], while third and fourth infections appear much less likely to be symptomatic [33] (but see [34] who found similar rates of symptomatic cases across infection number). Despite these potential differences in reporting rate, modeling work has suggested that the dynamics of primary and secondary infections are closely coupled (and thus not dynamically distinct) under many conditions [35]. Moreover, given the short time scale of our data (5 years) relative to the period of cross-immunity (roughly 2 years [25]), we expect that third and fourth infections will be relatively rare. We thus expect the SEIRS framework, and the assumption of equal symptomatic rates across infection number, to be sufficient for capturing the dengue dynamics of Vitória and the relationship between mosquito abundance and human disease. In the SEIRS framework, the total human population of Vitória (N) is divided into susceptible (S), exposed (E), infectious (I), and immune (R) classes. Susceptible humans (S) become exposed (E) through contact with infectious mosquitoes (VI). Following a latent period (1ρ), exposed humans become infectious (I) at which point they can infect susceptible mosquitoes (VS). Infectious humans recover at rate γ and subsequently remain immune (R) for a period (1δ) after which they re-enter the susceptible class. Similarly, susceptible mosquitoes (VS) become exposed (VE) by biting infectious humans and pass through a temperature-dependent incubation period (1ρv(t)) before becoming infectious (VI). Because the assumption of an exponentially distributed incubation period (implicit in the specification of a differential equation model) is a poor fit to laboratory observations [36], we instead implemented a gamma-distributed incubation period by chaining together multiple exposed classes (VEj, taking advantage of the fact that a gamma-distributed random variable can be generated through the sum of exponential random variables with the same rate parameter) [37]. Total mosquito population size (VN) is controlled by a forced, seasonally varying growth rate (r(t)), while the transmission bottleneck is captured with a forced, seasonally varying mortality rate (d(t)). Captured mosquitoes (VC) accumulate at rate ϕqτ(t), where ϕq is the per-trap capture rate, and τ(t) is the number of traps deployed in week t. We specified the differential equations governing the human population as: dSdt=bN-bS-λVINS+δR(1) dEdt=λVINS-(ρ+b)E(2) dIdt=ρE-(γ+b)I(3) dRdt=γI-(δ+b)R(4) while the equations governing the mosquito (vector) population are: dVNdt=r(t)VN-ϕqτ(t)VN(5) dVE1dt=λINVS-(4ρv(t)+d(t)+ϕqτ(t))VE1(6) dVE2dt=4ρv(t)VE1-(4ρv(t)+d(t)+ϕqτ(t))VE2(7) dVE3dt=4ρv(t)VE2-(4ρv(t)+d(t)+ϕqτ(t))VE3(8) dVE4dt=4ρv(t)VE3-(4ρv(t)+d(t)+ϕqτ(t))VE4(9) dVIdt=4ρv(t)VE4-(d(t)+ϕqτ(t))VI(10) dVCdt=ϕqτ(t)VN(11) VS=VN-VE-VI.(12) We modeled the centered and log-transformed mosquito mortality rate (ν) and the per-capita mosquito growth rate (r) as forced harmonic oscillators with natural periods of one year: d2νdt2=-ω2ν+ϵνt(13) d2rdt2=-ω2r+ϵrt,(14) where the angular frequency of the oscillator, ω = 2π/52, the mosquito death rate d(t) = d0 exp(ν(t)), and ϵνt∼Normal(0,σν2)(15) ϵrt∼Normal(0,σr2),(16) for each week t = 1, …, 243. These stochastically-forced harmonic oscillators provide a flexible framework for generating smooth seasonal oscillations in the latent mosquito processes [38]. Data model To connect the differential equation model to the observed case reports, we added an extra state, C, that collects the cumulative number of transitions from the exposed to infectious class (assuming that case reporting coincides with the onset of symptoms). We then modeled the number of new cases reported in week t (yt) as: yt∼NegBin(ϕy(C(t)-C(t-1)),ηy),(17) where ϕy is the reporting probability, C(t) − C(t − 1) is the number of new infectious humans in week t, and ηy controls the overdispersion relative to the Poisson distribution. We similarly modeled the number of mosquitoes trapped in week t (qt) as: qt∼NegBin(VC(t)-VC(t-1),ηq),(18) where VC(t) − VC(t − 1) is the number of new mosquitoes captured in week t, and ηq controls overdispersion relative to the Poisson distribution. Parameterization and priors Several of the parameters in this model are assumed to be fixed and known (Table 1). The human population size and average life span (which we use to parameterize the birth/death rate) for Vitória were taken from the 2010 census. To maintain identifiability, the transmission rate (λ) was also fixed at literature values. Lastly, the extrinsic incubation period in mosquitoes was modeled as a function of weekly mean temperature and forced with weather station data obtained from WeatherUnderground [39]. 10.1371/journal.pntd.0008868.t001Table 1 Model parameters and their values. The values in parentheses after the posterior means give the 80% credible interval. See S1 Text for a full description of all prior distributions. Parameter Description Prior mean Posterior mean Citation N Human population size in Vitória, Brazil 327801 [21] 1/d Human life-span 76 years [40] λ Transmission rate 4.87 week−1 [41] 1/ρv(t) Extrinsic incubation period 17exp(7.9-0.21T(t)) weeks [36] VE0 Initial exposed mosquitoes 0 VI0 Initial infectious mosquitoes 0 d0 Baseline mosquito mortality rate 1.47 week−1 0.88 (0.7, 1.1) [42] 1/ρ Latent period in host 0.87 weeks 1.72 (1.2, 2.3) [36] γ Rate of loss of infectiousness 3.5 week−1 3.6 (3.2, 4.1) [43] 1/δ Period of cross-immunity 97 weeks 114 (72, 160) [25] σr Standard deviation of mosquito growth rate forcing 0* 0.013 (0.01, 0.02) σν Standard deviation of mosquito mortality rate forcing 0* 0.0005 (0.0004, 0.0008) S0 Proportion initially susceptible 0.4 0.42 (0.28, 0.57) [44] E0 Number initially exposed 100 148 (104, 196) I0 Number initially infectious 60 79 (45, 116) r0 Initial mosquito population growth rate 0 0.008 (-0.01, 0.02) ν0 Initial unconstrained mosquito mortality rate 0 -0.16 (-0.4, 0.08) VN0 Initial mosquito population size 2N 1.5N (1.1N, 1.9N) ϕy Reporting probability 0.083 0.14 (0.1, 0.18) [45] log(ϕq) Log per-trap mosquito capture rate −13 -13.2 (-13.6, -13) ηy Overdispersion of case reports 0* 0.12 (0.1, 0.13) ηq Overdispersion of mosquito trap counts 0* 0.14 (0.12, 0.16) * indicates prior mode, rather than mean. The remaining parameters include the epidemiological parameters controlling the average latent, infectious, and immune periods (ρ, γ, δ) and average mosquito lifespan (d0), the initial conditions of the model (S0, E0, I0, R0, VN0, ν0, r0), the variances of the latent mosquito processes (σr2, σν2), and the remaining measurement parameters (ϕy, ϕq, ηy, ηq). Where possible, we specified informative prior distributions for these parameters based on existing laboratory and field studies (see Table 1 for means and S1 Text for detailed explanations). Implementation Combining the data, process, and parameter models [19], we summarize the full hierarchical model as: yt|·∼NegBin(ϕy(C(t)-C(t-1)),ηy),(19) qt|·∼NegBin(VC(t)-VC(t-1),ηq)(20) (C,VC)=M(ϵr,ϵν,θ)(21) ϵrt∼Normal(0,σr2)(22) ϵνt∼Normal(0,σν2)(23) θ∼[θ](24) where θ is a vector of all the model parameters and initial conditions, and M(ϵr,ϵν,θ) represents the (numeric) solution to the differential equation model (Eqs 1–14) as a function of θ and the weekly stochastic forcing terms (ϵr, ϵν). Sampling from the posterior distribution of the parameters in a mechanistic model is difficult due to multimodality, variable parameter sensitivities (e.g., small changes in one parameter may lead to large changes in output, while similar changes in another parameter may have little effect), and potentially strong posterior correlations induced by the nonlinearity of the differential equation model [14, 46, 47]. However, the variable ϵt introduces flexibility to the mechanistic model that remedies lack-of-fit when the process parameters are far from optimal [48], thereby reducing multimodality and helping to smooth the posterior surface. Gradient-based methods like Hamiltonian Monte Carlo (HMC) can then more easily and efficiently traverse the posterior. Samples from the posterior distribution were generated using HMC implemented in the rstan package [49, 50] for R [51]. We ran 3 chains with different starting values for 4,000 iterations each, discarding the first 2,000 as burn-in. Convergence diagnostics and mixing were evaluated using the shinystan package [52]. In our implementation, the solution to the differential equation model was approximated with an Euler scheme with a time step of 1 day. Timesteps as small as 1/8 of a day were explored and did not qualitatively change the modeled dynamics. Code is available from https://github.com/clint-leach/mosquito-recon. Mosquito control simulations Given a subset of the samples from the posterior distribution as obtained above (2000, taken to reduce computation time), we simulated the effects of a single pulse of mosquito control applied in each week of the first three years of the time series. Because the city already implements responsive, targeted control with the aim of reducing local mosquito density during an existing outbreak, we focused our simulations on exploring the longer-term feedbacks induced by mosquito control and the ability of an intervention to reduce the disease burden over the following year. For each week and each posterior sample, we simulated the dynamics resulting from a 5% reduction in the mosquito population implemented at the beginning of that week (affecting susceptible, exposed, and infectious mosquitoes equally). To capture the likely rapid rebound in mosquito abundance following a single pulse of control [28, 53], we also simulated a 5% increase in mosquito birth rate in the following week to return mosquito abundance to its previous trajectory (without this, the 5% reduction in abundance persists indefinitely). We then compared the number of cases produced over the year following the control intervention in the control scenario to the same number in the uncontrolled scenario. A 5% reduction in mosquito abundance was chosen to keep our simulations conservative relative to field estimates of the mortality induced by spraying [54], and to avoid pushing the model into the unrealistic range of dengue eradication. Results The model captured the observed dynamics of both case reports and mosquito trap counts (Fig 1). The estimated posterior median case reports explained 91% of the variation in the observed time series, while the posterior median mosquito trap counts explained 46% of the variation in the observed time series. In addition, posterior predictive checks showed that the model reproduced the total number of cases reported and mosquitoes captured as well as the autocorrelation structure of both time series (with the exception of slightly underestimating the autocorrelation for short lags, S6, S7 and S8 Figs). The posterior distributions of the rate of infectious decay (γ) and the period of cross-immunity (1/δ) did not differ substantially from their priors, suggesting that the Vitória data contained little additional information about these parameters (S3 Fig). The estimated latent period in a human host (1/ρ, the expected time it takes for an exposed human to become infectious to biting mosquitoes) was influenced more strongly by the data, with a posterior mean of 1.73 weeks compared to a prior mean of 0.87 weeks. Further, the posterior mean case reporting rate (ϕ) was 0.14, larger than the prior mean of 0.08. 10.1371/journal.pntd.0008868.g001Fig 1 Vitória data and model estimates. A: weekly observed case reports (points), with corresponding posterior median (black line) and 80% posterior credible interval (gray band). B: weekly mosquito trap counts (points), with posterior median (black line) and 80% posterior credible interval (gray band). C: extrinsic incubation period (EIP; weeks), computed from weekly mean temperature data. D: estimated weekly mosquito mortality rate, with the posterior median (black line) and the 80% posterior credible interval (gray band). The estimated weekly mosquito mortality rate varied seasonally, with generally high mortality early in the year and low mortality in August to October (Fig 1). This seasonal trend broadly tracked seasonal variation in temperature (Fig 2, correlation coefficient of 0.64) and mosquito trap counts (Fig 1), though the shape of the annual trajectory differed from year to year. The posterior distribution of the baseline mortality rate (d0) had a mean of 0.88/week, roughly 60% of the prior mean. The marginal posterior means of the ϵνt forcing the mosquito mortality process exhibited a higher-frequency periodic oscillation (S1 Fig), although the marginal posterior distribution of each ϵνt overlapped zero. In addition, the standard deviation of the mortality forcing terms was small relative to the weak prior (E(σν|y) = 0.01, S2 Fig). 10.1371/journal.pntd.0008868.g002Fig 2 Mosquito mortality and temperature. Posterior median mosquito mortality rate as a function of weekly mean temperature (degrees Celsius). The effect of a given mosquito control intervention (i.e., the temporary removal of 5% of the adult population in a given week) on the number of cases in the following year (relative to no control) varied both seasonally and interannually (Fig 3). This variation in the effect of control was tightly correlated with the estimated mosquito mortality, with a median posterior correlation between the two time series of 0.96. As such, the seasonal variation in the effectiveness of control followed the same trend as mosquito mortality rate (although the annual minima in the case ratio time series generally fell 1 to 2 weeks before the minima in the mosquito mortality time series), with the largest reductions in case load resulting from interventions during the dengue off-season (early September for 2008 and 2009, and mid-July for 2010). Summing over the interannual variation to compute the overall effect of control implemented in a given week of the year, we found that mosquito control was most effective when implemented around week 34 (late August/early September), reducing the estimated case load by roughly 14% (Fig 3). This broadly corresponds to the end of the dry season in Vitória, when both dengue case reports and mosquito trap counts are low. 10.1371/journal.pntd.0008868.g003Fig 3 The effect of mosquito control as a function of the week in which it was applied. Summary of the posterior predicted effect of mosquito control implemented in a given week of the year on the number of cases in the following year (relative to the number of cases expected without control). Black lines indicate the posterior median, while gray ribbons indicate the 80% credible interval. The first three panels show the results for control implemented in the years 2008-2010, and the last panel shows the overall effect of control impelmented in a given week of the year, summing over all three years. For example, mosquito control applied in week 37 of 2008 would have prevented about 13% of the human cases over the following year (i.e., the caseload would have been 87% of the expectation without control). The variation associated with the posterior predicted effect of control implemented in week 34 of the year (i.e., the width of ribbon in Fig 3) was correlated with the baseline mosquito mortality rate (d0, posterior correlation coefficient of 0.65) and the case reporting probability (ϕ, posterior correlation coefficient of -0.19). Simulated mosquito control created the largest reduction in case reports in posterior samples with low baseline mosquito mortality rate and/or high reporting probability, while control was relatively less effective in simulations from samples with high mosquito mortality or low reporting probability (Fig 4). Thus mosquito control was more effective at reducing disease burden in simulations with long average mosquito lifespans (i.e., low mosquito mortality rates) or low overall prevalence (i.e., fewer undetected cases). 10.1371/journal.pntd.0008868.g004Fig 4 Posterior correlation between system parameters and the effectiveness of control. The y-axis represents the effectiveness of optimally timed control, i.e., the effect of control implmented in the 34th week of the year on the relative number of cases in the following year, summed over 2008, 2009, and 2010. Each point represets a single sample from the posterior distribution, giving the number of cases in the controlled simulation (relative to the number of cases expected without control) as a function of A: the mean mosquito mortality rate, d0, and B: the case reporting probability, ϕ, from that posterior sample. Discussion Processes driving effect of mosquito control The dynamics of dengue fever, like those of many infectious diseases [9, 12] and ecological systems [55], are driven by the combined efforts of intrinsic non-linearities, seasonality, and stochasticity. Seasonality, in particular, is an important factor in capturing the annual cycle of dengue outbreaks [24, 25, 56]. However, the observed seasonality in transmission likely emerges from the combined effects of multiple seasonally-varying components that may be driven by different environmental factors that oscillate in different phases (e.g., mosquito abundance seems to lag slightly behind temperature-driven variation in extrinsic incubation period, Fig 1). Integrating these seasonally-varying components into a synthetic measure of transmission potential (e.g., a temperature-dependent effective reproduction number, [57]), or more specifically, a measure of the transmission potential of mosquitoes, is difficult. We positioned the latent mosquito mortality as the link between mosquito abundance, the extrinsic incubation period, and human cases. In this way, the estimated mosquito mortality rate serves as an index of transmission potential, opening or closing the mosquito life history bottleneck [11] as necessary to fit to the case reports data. The resulting seasonality in the estimated trajectory suggests that the seasonality in mosquito abundance and the extrinsic incubation period was not sufficient to capture the observed case reports. In particular, when mosquitoes were relatively scarce and transmission limited, we estimated a relatively low mosquito mortality, suggesting that long-lived mosquitoes were required to maintain observed levels of transmission through the off-season. On the other hand, when mosquitoes were abundant, we estimated relatively high mosquito mortality rates, suggesting that transmission needed to be damped. The importance of mosquito longevity in driving disease dynamics highlights the potential effectiveness of control efforts that target adult mosquitoes and disrupt transmission by preventing mosquitoes from living long enough to progress through the extrinsic incubation period to the infectious state and bite a susceptible human. In fact, this forms the basis for much of the theory of adult mosquito control [3, 11, 28]. The high correlation between our estimated mosquito mortality and the effect of control confirms this theory, suggesting that control is most effective when it targets long-lived mosquitoes. Specifically, our simulated mosquito control interventions were most effective at reducing the disease burden when applied around week 34 (i.e., early September), in the dengue off-season. The effectiveness of this control was likely driven by the fact that transmission during the off-season was already limited by low mosquito abundance and a relatively high extrinsic incubation period. Given that limitation, transmission was maintained by relatively few long-lived mosquitoes, making the system vulnerable to perturbation. On the other hand, we found that a single pulse of control was relatively less effective when implemented during an outbreak, when mosquitoes were abundant (and unlikely to be limiting transmission) but short-lived. Given the relatively high mosquito mortality rates during this time, exposed and infectious mosquitoes were already fairly ephemeral, such that the relatively small disruption induced by control likely made little difference. Moreover, due to the large number of infectious human hosts available to transmit to the remaining (and rapidly rebounding) mosquito population, the population of exposed mosquitoes likely recovered quickly [27, 28]. Barsante et al. [58] and Oki et al. [59] similarly found that control was most effective when applied well before peak prevalence, either during the dry season [58] or early in the rainy season [59] (September is near the end of the dry season in Vitória). While the immediate effects of mosquito control implemented during the decline phase of an outbreak may be masked by the natually fading transmission intensity [60], our results nonetheless indicated that disrupting inter-seasonal transmission can be an effective longer-term strategy [61]. Further, although large pulses of imported cases could potentially swamp the effects of early control, in additional simulations we found that our results were robust to the import of 10 infectious humans (roughly the same order as the number of locally reported cases) just before the annual outbreak. The Bayesian framework allowed us to account for uncertainty across the data, process, and parameter levels of our model [19]. We carried this uncertainty through to our simulations of control interventions [18] and found that there was substantial uncertainty in the proportion of cases prevented by a control intervention (i.e., the width of the ribbons in Fig 3). Much of this uncertainty could be attributed to posterior uncertainty in the case reporting rate (ϕ) and the average mosquito mortality rate (d0). Specifically, we found that control implemented at the overall optimum (week 34) had the largest impact (i.e., the lowest case ratio) in simulations with a low average mosquito mortality and/or a high case reporting rate (Fig 4). Mosquito mortality and the case reporting rate were correlated with the overall level of mosquito abundance and the overall size of the susceptible population, respectively, suggesting that control was most effective in simulations with fewer mosquitoes and more susceptible humans. Hladish et al. [61] similarly found that simulated indoor residual spraying was more effective when the modeled mosquito abundance was already low. This suggests that control efforts that reduce the ability of mosquitoes to transmit DENV are most effective for situations in which mosquito abundance is already the limiting component to maintaining transmission (relative to other factors like human immunity). The posterior correlation between the effectiveness of control, the case reporting rate, and the average mosquito mortality rate emphasizes that the impact of mosquito control is jointly regulated by both mosquito population dynamics and human immune processes. Similar observations were made by ten Bosch et al. [62] who found that models with longer periods of cross-immunity (such that susceptibles replenished more slowly) generated systems in which transmission was more difficult to disrupt with control actions. As a result of these relationships, efficient deployment of mosquito control, and accurate prediction of its effects, is likely to depend in part on our ability to monitor and predict the dynamics of human immunity. To meet these needs, existing mosquito monitoring efforts need to be paired with more detailed clinical surveillance [3] and tighter estimates of the period of cross-immunity [62] and the number of unreported cases [45]. In the absence of, or as a supplement to, such data, mechanistic models like the one developed here, or so-called TSIR (Time-series Susceptiple-Infected-Recovered) frameworks that reconstruct the dynamics of the susceptible class [25, 63], need to be further developed to better inform and understand mosquito control efforts. Interpretation of the estimated mosquito mortality rate The positive correlation between the estimated mosquito mortality rate and the simulated effect of control suggests that mosquito control was most effective when it targeted long-lived mosquitoes during the inter-epidemic periods when mosquito abundance was low. However, the fact that the mosquito mortality forcing terms (ϵν, Eqs 13 and 16) were the only source of variability in the transmission process implies that the estimated mosquito mortality time series could have absorbed other sources of stochasticity or model misspecification. Hooker and Ellner [64] provide a framework for diagnosing such model misspecification in differential equation models using forcing functions similar to our implementation of the ϵν. In that framework, Hooker and Ellner [64] estimate nonparametric forcing functions that modify a fitted differential equation model to provide a good fit to the data. These forcing functions serve as residuals on the time derivatives, and can be more readily interpreted as indicators of lack-of-fit than residuals on the state variables [16, 64]. We do not employ the same explicit goodness-of-fit testing framework as [64], but we can inspect our estimated ϵν forcing terms in the same spirit. The periodic structure in the time series of the posterior means of the ϵν (S1 Fig) suggests that these terms were accounting for more than just noise, and there may have been some unmodeled process influencing fluctuations in mosquito mortality and/or transmission. Following Hooker and Ellner, we can explore whether this process is likely to result from misspecification of the rates of change of the existing state variables (indicated by a dependence of ϵνt on other state variables), or from missing state variables altogether (indicated by an additional dependence of ϵνt on its own lagged values). The lack of any apparent relationship between the forcing terms and any of the estimated state variables, combined with the dependence of ϵνt on previous values (as apparent through the periodic structure), suggest that unmodeled state variables may be the more likely driver of model misspecification. These unmodeled components could include additional mosquito population dynamic processes (e.g., aquatic stage dynamics, environmental drivers, or control interventions), or epidemiological processes (e.g., multiple circulating serotypes of the dengue virus, subsets of the population with different mixing or risk levels). Despite these potential sources of model misspecification, our estimated mosquito mortalities nonetheless fell within the reasonable range from the literature [42, 65]. Moreover, the fact that our estimated mortality rate increased with temperature also broadly agrees with the empirical literature on mosquito survival [42, 66, 67]. This suggests that regardless of unexplained structure in the forcing terms, the pattern of case reports was still very well described by realistic seasonal fluctuations in the mosquito mortality rate. As demonstrated by Reiner et al. [68] for malaria transmission, estimates of transmission potential can be sensitive to fluctuations in mosquito abundance and age structure. Moreover, given the broad importance of seasonality in understanding dengue epidemiology [62], and the role of mosquito mortality and age in driving the effect of control interventions, future work should focus on developing a more complete, predictive understanding of the seasonal drivers of mosquito mortality (including, potentially, control itself). Additional considerations and extensions In addition to epidemiological complexity, dengue dynamics are further complicated by the reporting process. We estimated a relatively long average intrinsic incubation period (1/ρ, posterior mean of 1.7 weeks) relative to our prior mean (0.87 weeks), suggesting possible reporting delays [69]. Moreover, the case data to which we fit the model represent reports of “dengue-like illness,” without laboratory confirmation, and as such could include cases of other diseases with similar symptoms (e.g., chikungunya or Zika). However, neither chikungunya nor Zika had emerged as substantial public health threats in Brazil by the end of our time series in late 2012 [70, 71]. In addition, given the high underreporting rate expected for dengue fever [45], and the uncertainty incorporated into the measurement model, we expect misreported cases to have a small effect on our analyses. Mosquito control interventions can prevent cases by acting on any of the components of vectorial capacity. We focused on the direct effect of killing adult mosquitoes on transmission, but adult control can also act by reducing egg laying and the number of mosquitoes in the next generation [72]. Given the relative simplicity of our mosquito model, we were unable to explore the feedbacks that adult control may induce in mosquito population dynamics, and instead assumed that mosquito populations quickly rebound from any perturbations [53]. We expect that capturing these feedbacks would likely reinforce our conclusions about the utility of off-season control, as the disruption to mosquito population dynamics would more strongly limit the ability of the mosquito population to maintain transmission through the off-season. Although our results suggest that a pulse of adult control in the off-season may be an effective tool for preventing human cases, achieving particular control thresholds or policy goals will likely require deploying a combination of control interventions [72]. In fact, it is important to note that the data to which we fit our model implicitly reflect the control efforts already enacted by the city (the effects of which may have influenced our estimates of mosquito demographic rates). Thus, our mosquito control simulations should be interepreted as exploring the effect of an additional pulse of city-wide control in addition to the existing control activities. These control efforts are guided by the MI-Dengue system, which uses the trap-level mosquito surveillance data to target areas of high mosquito infestation for control (including source reduction, larvacide, and adulticide) [22]. Although these targeted, reactive interventions are necessary to help reduce local disease risk at times and locations of high mosquito abundance [73], our results suggest that an additional pulse of proactive control in the off-season when mosquitoes are less abundant would minimize human cases. The frequent use of spatially-targeted mosquito control highlights the potential for spatial heterogeneities in disease risk within a city. In particular, structured human movement within the city is likely to induce heterogeneous human-mosquito mixing [74–76]. In addition, spatial variability in socioeconomic factors within the city may also modulate the extent to which mosquitoes in different parts of the city contribute to disease spread [77–80]. Although our model does not account for these heterogeneities, it was nonetheless able to capture the city-wide dynamics well, suggesting that there may be sufficient mixing to appear homogeneous at the city scale. Further, although our mechanistic model may be able to suggest when a city-wide intervention is likely to be effective, to make the best use of limited resources, spatial prioritization may still be necessary. Conclusions Efforts to connect mosquito abundance to human disease are often hampered by the confounding influences of human immunity and mosquito survival. The challenges presented by these confounding factors highlight the value of mechanistic information in studying the effect of mosquito control on disease spread. In the Bayesian context we deployed, this mechanistic knowledge, as formalized in the specification of a differential equation model, can be viewed as part of the prior knowledge on the relationship between mosquito abundance and human disease [12, 17]. As such, we should neither ignore this mechanistic information, nor encode it so rigidly that it overwhelms the signal in our data. We developed a simple yet realistic mechanistic model of dengue fever spread that represents the fundamental elements of our prior understanding of dengue epidemiology, while also allowing for uncertainty and flexibility in the fluctuations of mosquito demographic rates. This mechanistic framework allowed us to capture the critical contribution of long-lived, off-season mosquitoes to the maintenance of transmission and to identify critical intervention points that would not be apparent otherwise. The fully hierarchical Bayesian framework in which we embedded the mechanistic model allowed for a thorough accounting of uncertainty that was carried through to the evaluation of different control strategies. This combination of model features helps to meet the need for more effective, biologically grounded, and data-driven dengue control policies and offers a building block on which these tools can be further developed in the future. Supporting information S1 Fig Posterior estimates of harmonic oscillator forcing functions. A: mosquito mortality forcing. B: mosquito growth forcing. Median posterior estimate (black line) and 80% credible interval (gray band). (EPS) Click here for additional data file. S2 Fig Posterior distribution of variances on forcing terms. A: mosquito mortality forcing. B: mosquito growth forcing. (EPS) Click here for additional data file. S3 Fig Posterior distribution of epidemiological parameters. The black line indicates the prior density. A: scaled latent period in human host (1/ρ). B: scaled rate of infectious decay in human host (γ). C: scaled period of cross-immunity (1/δ). D: scaled baseline mosquito mortality rate (d0). (EPS) Click here for additional data file. S4 Fig Posterior distribution of initial conditions. The black line indicates the prior density. A: initial proportion of the human population that is susceptible. B: initial number of exposed individuals. C: initial number of infectious individuals. D: log initial number of mosquitoes-per-person. E: initial value of centered mosquito mortality rate oscillator. F: initial value of mosquito population growth rate oscillator. (EPS) Click here for additional data file. S5 Fig Posterior distribution of measurement parameters. The black line shows the prior density. A: log of the mosquito trap capture rate (which is strongly correlated with initial mosquito abundance). B: human case reporting probability. C: mosquito trap count measurement over-dispersion (prior density too low to be visible). D: case report measurement over-dispersion (prior density too low to be visible). (EPS) Click here for additional data file. S6 Fig Posterior check of autocorrelation structure in data. A: autocorrelation in case report time series. B: autocorrelation in trapped mosquito time series. Points indicate the observed autocorrelation at the given lag, and vertical lines give the 80% posterior credible interval of the autocorrelation in the estimated time series. (EPS) Click here for additional data file. S7 Fig Posterior distribution of the total number of cases reported (A) and the total number of mosquitoes captured (B). The thick line indicates the observed totals. (EPS) Click here for additional data file. S8 Fig Posterior distributions of the minimum and maximum weekly reported cases (A and C) and trapped mosquitoes (B and D). The thick line indicates the observed minima and maxima. (EPS) Click here for additional data file. S9 Fig Posterior estimates of the mosquito mortality rate and the probability of a mosquito surving the external incubation period. The black line indicates the posterior median, while the gray ribbon shows the 80% credible interval. The probability of surviving the extrinsic incubation period (bottom panel) was computed by taking the ratio of the number of mosquitoes that entered the infectious class in a week over the total number of mosquitoes that exited the exposed class that week. That is, the probability of surviving the extrinsic incubatio period represents the fraction of mosquitoes leaving the exposed class that enter the infectious class. (EPS) Click here for additional data file. S1 Text Full specification of the model and prior distributions. (PDF) Click here for additional data file. CBL would like to thank members of the Webb Lab for feedback and support through innumerable previous iterations of this project. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government. ==== Refs References 1 Bhatt S , Gething PW , Brady OJ , Messina JP , Farlow AW , Moyes CL , et al The global distribution and burden of dengue . Nature . 2013 ;496 :504 –507 . 10.1038/nature12060 23563266 2 Achee NL , Gould F , Perkins TA , Reiner RC , Morrison AC , Ritchie Sa , et al A critical assessment of vector control for dengue prevention . PLoS Negl Trop Dis . 2015 ;9 (5 ):e0003655 10.1371/journal.pntd.0003655 25951103 3 Morrison AC , Zielinski-Gutierrez E , Scott TW , Rosenberg R . Defining challenges and proposing solutions for control of the virus vector Aedes aegypti . PLoS Med . 2008 ;5 (3 ):0362 –0366 . 10.1371/journal.pmed.0050068 18351798 4 Scott TW , Morrison AC . Longitudinal field studies will guide a paradigm shift in Dengue prevention In: Vector Biology, Ecology, and Control ; 2010 p. 139 –161 . http://link.springer.com/10.1007/978-90-481-2458-9. 5 Scott TW , Morrison AC . Vector dynamics and transmission of dengue virus: implications for dengue surveillance and prevention strategies . Dengue Virus . 2010 ;(338 ):161 19802582 6 Bowman LR , Runge-Ranzinger S , McCall PJ . Assessing the relationship between vector indices and dengue transmission: A systematic review of the evidence . PLoS Negl Trop Dis . 2014 ;8 (5 ):e2848 10.1371/journal.pntd.0002848 24810901 7 Pepin KM , Leach CB , Marques-Toledo C , Laass KH , Paixao KS , Luis AD , et al Utility of mosquito surveillance data for spatial prioritization of vector control against dengue viruses in three Brazilian cities . Parasites and Vectors . 2015 ;8 (1 ):1 –15 . 10.1186/s13071-015-0659-y 25561160 8 Cromwell EA , Stoddard ST , Barker CM , Van Rie A , Messer WB , Meshnick SR , et al The relationship between entomological indicators of Aedes aegypti abundance and dengue virus infection . PLoS Negl Trop Dis . 2017 ;11 (3 ):e0005429 10.1371/journal.pntd.0005429 28333938 9 Koelle K , Pascual M . Disentangling extrinsic from intrinsic factors in disease dynamics: A nonlinear time series approach with an application to cholera . Am Nat . 2004 ;163 (6 ):901 –913 . 10.1086/420798 15266387 10 Koelle K , Rodó X , Pascual M , Yunus M , Mostafa G . Refractory periods and climate forcing in cholera dynamics . Nature . 2005 ;436 :696 –700 . 10.1038/nature03820 16079845 11 Smith DL , Battle KE , Hay SI , Barker CM , Scott TW , Mckenzie FE . Ross, Macdonald, and a theory for the dynamics and control of mosquito-transmitted pathogens . PLoS Pathog . 2012 ;8 (4 ). 10.1371/journal.ppat.1002588 12 Ellner SP , Bailey Ba , Bobashev GV ,Gallant AR , Grenfell BT , Nychka DW . Noise and nonlinearity in measles epidemics: combining mechanistic and statistical approaches to population modeling . Am Nat . 1998 ;151 (5 ):425 –440 . 10.1086/286130 18811317 13 Grenfell BT , Bjørnstad ON , Finkenstädt BF . Dynamics of measles epidemics: Scaling noise, determinism, and predictability with the TSIR model . Ecol Monogr . 2002 ;72 (2 ):185 –202 . 10.1890/0012-9615(2002)072[0185:DOMESN]2.0.CO;2 14 Reilly C , Zeringue A . Improved predictions of lynx trappings using a biological model In: Gelman A , Meng XL , editors. Applied Bayesian Modeling and Causal Inference from Incomplete-Data Perspectives . John Wiley & Sons Ltd ; 2005 p. 297 –308 . 15 Buckland ST , Newman KB , Fernández C , Thomas L , Harwood J . Embedding population dynamics models in inference . Stat Sci . 2007 ;22 (1 ):44 –58 . 10.1214/088342306000000673 16 Hotelling H . Differential equations subject to error, and population estimates . J Am Stat Assoc . 1927 ;22 (159 ):283 –314 . 10.1080/01621459.1927.10502963 17 Wikle CK , Hooten MB . A general science-based framework for dynamical spatio-temporal models . Test . 2010 ;19 :417 –451 . 10.1007/s11749-010-0214-2 18 Elderd BD , Dukic VM , Dwyer G . Uncertainty in predictions of disease spread and public health responses to bioterrorism and emerging diseases . Proc Natl Acad Sci U S A . 2006 ;103 (42 ):15693 –15697 . 10.1073/pnas.0600816103 17030819 19 Berliner LM . Hierarchical Bayesian time series models In: Hanson KM , Silver RN , editors. Maximum entropy and Bayesian methods. Fundamental Theories of Physics (An international book series on the fundamental theories of physics: Their clarification, development and application) . vol. 79 Springer , Dordrecht ; 1996 . 20 Cressie N , Calder Ca , Clark JS , Ver Hoef JM , Wikle CK . Accounting for uncertainty in ecological analysis: the strengths and limitations of hierarchical statistical modeling . Ecol Appl . 2009 ;19 (3 ):553 –70 . 10.1890/07-0744.1 19425416 21 Demographic Census 2000 and 2010 and population estimates 2001 to 2009 and 2011 to 2014. IBGE;. http://legado.vitoria.es.gov.br/regionais/Censo_2010/Tab17_crespopulacional.asp. 22 Eiras Á , Resende M . Preliminary evaluation of the “Dengue-MI” technology for Aedes aegypti monitoring and control . Cadernos de Saúde Pública . 2009 ; p. 45 –58 . 10.1590/S0102-311X2009001300005 19287866 23 Lana RM , Morais MM , de Lima TFM , Carneiro TGdS , Stolerman LM , dos Santos JPC , et al Assessment of a trap based Aedes aegypti surveillance program using mathematical modeling . PLoS ONE . 2018 ;13 (1 ):1 –16 . 10.1371/journal.pone.0190673 29304070 24 Wearing HJ , Rohani P . Ecological and immunological determinants of dengue epidemics . Proc Natl Acad Sci U S A . 2006 ;103 (31 ):11802 –7 . 10.1073/pnas.0602960103 16868086 25 Reich NG , Shrestha S , King Aa , Rohani P , Lessler J , Kalayanarooj S , et al Interactions between serotypes of dengue highlight epidemiological impact of cross-immunity . J R Soc Interface . 2013 ;10 (86 ):20130414 10.1098/rsif.2013.0414 23825116 26 Aguiar M , Kooi BW , Rocha F , Ghaffari P , Stollenwerk N . How much complexity is needed to describe the fluctuations observed in dengue hemorrhagic fever incidence data? Ecol Complex . 2013 ;16 :31 –40 . 10.1016/j.ecocom.2012.09.001 27 Newton EAC , Reiter A . A model of the transmission of dengue fever with an evaluation of the impact of ultra-low volume (ULV) insecticide applications on dengue epidemics . Am J Trop Med Hyg . 1992 ;47 (6 ):709 –720 . 10.4269/ajtmh.1992.47.709 1361721 28 Burattini MN , Chen M , Chow A , Coutinho FAB , Goh KT , Lopez LF , et al Modelling the control strategies against dengue in Singapore . Epidemiol Infect . 2008 ;136 (03 ):309 –319 . 10.1017/S0950268807008667 17540051 29 Pinho STR , Ferreira CP , Esteva L , Barreto FR , Morato e Silva VC , Teixeira MGL . Modelling the dynamics of dengue real epidemics . Philos Trans R Soc A Math Phys Eng Sci . 2010 ;368 (1933 ):5679 –5693 . 10.1098/rsta.2010.0278 21078642 30 Hooten MB , Anderson J , Waller La . Assessing North American influenza dynamics with a statistical SIRS model . Spatial and Spatio-temporal Epidemiology . 2010 ;1 :177 –85 . 10.1016/j.sste.2010.03.003 21218144 31 Imai N , Dorigatti I , Cauchemez S , Ferguson NM . Estimating dengue transmission intensity from case-notification data from multiple countries . PLoS Neglected Tropical Diseases . 2016 ;10 (7 ):e0004833 10.1371/journal.pntd.0004833 27399793 32 Clapham HE , Cummings DAT , Johansson MA . Immune status alters the probability of apparent illness due to dengue virus infection: Evidence from a pooled analysis across multiple cohort and cluster studies . PLoS Neglected Tropical Diseases . 2017 ;11 (9 ):e0005926 10.1371/journal.pntd.0005926 28953902 33 Olkowski S , Forshey BM , Morrison AC , Rocha C , Vilcarromero S , Halsey ES , et al Reduced risk of disease during postsecondary dengue virus infections . The Journal of Infectious Diseases . 2013 ;208 :1026 –33 . 10.1093/infdis/jit273 23776195 34 Montoya M , Gresh L , Mercado JC , Williams KL , Vargas MJ , Gutierrez G , et al Symptomatic versus inapparent outcome in repeat dengue virus infections is influenced by the time interval between infections and study year . PLoS Neglected Tropical Diseases . 2013 ;7 (8 ):e2357 10.1371/journal.pntd.0002357 23951377 35 Schwartz IB , Shaw LB , Cummings DAT , Billings L , Mccrary M , Burke DS . Chaotic desynchronization of multistrain diseases . Physical Review E—Statistical, Nonlinear, and Soft Matter Physics . 2005 ;72 (6 ):066201 10.1103/PhysRevE.72.066201 16486034 36 Chan M , Johansson MA . The incubation periods of dengue viruses . PLoS One . 2012 ;7 (11 ):1 –7 . 10.1371/journal.pone.0050972 23226436 37 Lloyd AL . Realistic distributions of infectious periods in epidemic models: Changing patterns of persistence and dynamics . Theor Popul Biol . 2001 ;60 (1 ):59 –71 . 10.1006/tpbi.2001.1525 11589638 38 Ramsay J , Hooker G . Dynamic Data Analysis: Modeling Data with Differential Equations . Springer ; 2017 . 39 Historial weather data for station IESVITOR2, 2008-2012. WeatherUnderground;. https://www.wunderground.com/dashboard/pws/IESVITOR2. 40 Longevity Subindex. IDHM;. http://legado.vitoria.es.gov.br/regionais/indicadores/esperancadevidanascer.asp. 41 Scott TW , Amerasinghe PH , Morrison AC , Lorenz LH , Clark GG , Strickman D , et al Longitudinal studies of Aedes aegypti (Diptera: Culicidae) in Thailand and Puerto Rico: blood feeding frequency . J Med Entomol . 2000 ;37 (1 ):89 –101 . 10.1603/0022-2585-37.1.77 15218911 42 Brady OJ , Johansson MA , Guerra CA , Bhatt S , Golding N , Pigott DM , et al Modelling adult Aedes aegypti and Aedes albopictus survival at different temperatures in laboratory and field settings . Parasit Vectors . 2013 ;6 :351 10.1186/1756-3305-6-351 24330720 43 Nguyet MN , Duong THK , Trung VT , Nguyen THQ , Tran CNB , Long VT , et al Host and viral features of human dengue cases shape the population of infected and infectious Aedes aegypti mosquitoes . Proc Natl Acad Sci U S A . 2013 ;110 (22 ):9072 –7 . 10.1073/pnas.1303395110 23674683 44 Cardoso I , Cabidelle A . Dengue: clinical forms and risk groups in a high incidence city in the southeastern region of Brazil . Revista da Sociedade Brasileira de Medicina Tropical . 2011 ;44 (4 ):430 –435 . 10.1590/S0037-86822011005000044 21789354 45 Silva MMO , Rodrigues MS , Paploski IAD , Kikuti M , Kasper AM , Cruz JS , et al Accuracy of dengue reporting by national surveillance system, Brazil . Emerg Infect Dis . 2016 ;22 (2 ):336 –339 . 10.3201/eid2202.150495 26812472 46 Girolami M . Bayesian inference for differential equations . Theor Comput Sci . 2008 ;408 (1 ):4 –16 . 10.1016/j.tcs.2008.07.005 47 Calderhead B , Girolami M . Statistical analysis of nonlinear dynamical systems using differential geometric sampling methods . Interface Focus . 2011 ;1 (6 ):821 –835 . 10.1098/rsfs.2011.0051 23226584 48 Leander J , Lundh T , Jirstrand M . Stochastic differential equations as a tool to regularize the parameter estimation problem for continuous time dynamical systems given discrete time measurements . Math Biosci . 2014 ;251 (1 ):54 –62 . 10.1016/j.mbs.2014.03.001 24631177 49 Carpenter B , Gelman A , Hoffman M , Lee D , Goodrich B , Betancourt M , et al Stan: A Probabilistic Programming Language . J Stat Softw . 2016 ;VV (Ii ). 50 Stan Development Team. Rstan: the R interface to Stan, R package version 2.16.2; 2017. 51 R Core Team . R: A Language and Environment for Statistical Computing ; 2016 https://www.R-project.org/. 52 Gabry J. shinystan: Interactive visual and numerical diagnostics and posterior analysis for Bayesian models; 2018. https://CRAN.R-project.org/package=shinystan. 53 Focks DA , Kloter KO , Carmichael GT . The impact of sequential ultra-low volume ground aerosol applications of malathion on the population dynamics of Aedes aegypti (L.) . American Journal of Tropical Medicine and Hygiene . 1987 ;36 (3 ):639 –647 . 10.4269/ajtmh.1987.36.639 3578659 54 Esu E , Lenhart A , Smith L , Horstick O . Effectiveness of peridomestic space spraying with insecticide on dengue transmission; Systematic review . Trop Med Int Heal . 2010 ;15 (5 ):619 –631 . 20214764 55 Bjørnstad ON , Grenfell BT . Noisy clockwork: time series analysis of population fluctuations in animals . Science . 2001 ;293 (5530 ):638 –43 . 10.1126/science.1062226 11474099 56 Aguiar M , Ballesteros S , Kooi BW , Stollenwerk N . The role of seasonality and import in a minimalistic multi-strain dengue model capturing differences between primary and secondary infections: Complex dynamics and its implications for data analysis . J Theor Biol . 2011 ;289 :181 –196 . 10.1016/j.jtbi.2011.08.043 21907213 57 Codeço CT , Villela DAM , Coelho FC . Estimating the effective reproduction number of dengue considering temperature-dependent generation intervals . Epidemics . 2018 ;25 (May ):101 –111 . 29945778 58 Barsante LS , Cordeiro FS , Cardoso RTN , Acebal JL , Paixão KS ,Eiras ÁE . Controle da dengue na estação seca da cidade de Curvelo-MG através de um modelo preditivo . Proceedings of the Brazilian Society of Applied and Computational Mathematics . 2015 ;3 (1 ):1 –7 . 59 Oki M , Sunahara T , Hashizume M , Yamamoto T . Optimal timing of insecticide fogging to minimize dengue cases: Modeling dengue transmission among various Seasonalities and transmission intensities . PLoS Neglected Tropical Diseases . 2011 ;5 (10 ). 10.1371/journal.pntd.0001367 22039560 60 Stoddard ST , Wearing HJ , RC R Jr, Morrison AC , Astete H , Vilcarromero S , et al Long-term and seasonal dynamics of dengue in Iquitos, Peru . PLoS Neglected Tropical Diseases . 2014 ;8 (7 ):19 –21 . 10.1371/journal.pntd.0003003 25033412 61 Hladish TJ , Pearson CAB , Patricia Rojas D , Gomez-Dantes H , Halloran ME , Vazquez-Prokopec GM , et al Forecasting the effectiveness of indoor residual spraying for reducing dengue burden . PLoS Neglected Tropical Diseases . 2018 ;12 (6 ):1 –16 . 10.1371/journal.pntd.0006570 29939983 62 ten Bosch QA , Singh BK , Hassan MRA , Chadee DD , Michael E . The role of serotype interactions and seasonality in dengue model selection and control: Insights from a pattern matching approach . PLoS Neglected Tropical Diseases . 2016 ;10 (5 ):e0004680 10.1371/journal.pntd.0004680 27159023 63 Finkenstadt BF , Grenfell BT . Time series modelling of childhood diseases: a dynamical systems approach . Applied Statistics . 2000 ;49 :187 –205 . 64 Hooker G , Ellner SP . Goodness of fit in nonlinear dynamics: Misspecified rates or misspecified states? Ann Appl Stat . 2015 ;9 (2 ):754 –776 . 65 Maciel-de-Freitas R . Calculating the survival rate and estimated population density of gravid Aedes aegypti (Diptera, Culicidae) in Rio de Janeiro, Brazil . Cadernos de Saúde Pública . 2008 ;24 (12 ):2747 –2754 . 10.1590/S0102-311X2008001200003 19082265 66 Yang HM , Macoris MLG , Galvani KC , Andrighetti MTM , Wanderley DMV . Assessing the effects of temperature on the population of Aedes aegypti, the vector of dengue . Epidemiol Infect . 2009 ;137 (08 ):1188 10.1017/S0950268809002052 19192322 67 Morin CW , Comrie AC , Ernst K . Climate and dengue transmission: Evidence and implications . Environmental Health Perspectives . 2013 ;121 (11-12 ):1264 –1272 . 10.1289/ehp.1306556 24058050 68 Reiner RC , Guerra C , Donnelly MJ , Bousema T , Drakeley C , Smith DL . Estimating malaria transmission from humans to mosquitoes in a noisy landscape . Journal of the Royal Society Interface . 2015 ;12 (20150478 ). 69 Reich NG , Lauer SA , Sakrejda K , Iamsirithaworn S , Hinjoy S , Suangtho P , et al Challenges in Real-Time Prediction of Infectious Disease: A Case Study of Dengue in Thailand . PLoS Neglected Tropical Diseases . 2016 ;10 (6 ):1 –17 . 10.1371/journal.pntd.0004761 27304062 70 Pan American Health Organization, World Health Organization. Epidemiological alert: Chikungunya and dengue fever in the Americas; 2014. August. 71 Pan American Health Organization, World Health Organization. Epidemiological update: Zika virus infection; 2015. October. 72 Brady OJ , Godfray HCJ , Tatem AJ , Gething PW , Cohen JM , McKenzie FE , et al Vectorial capacity and vector control: Reconsidering sensitivity to parameters for malaria elimination . Transactions of the Royal Society of Tropical Medicine and Hygiene . 2016 ;110 (2 ):107 –117 . 10.1093/trstmh/trv113 26822603 73 Pepin KM , Marques-toledo C , Scherer L , Morais MM , Ellis B , Eiras AE . Cost-effectiveness of novel system of mosquito surveillance . Emerg Infect Dis . 2013 ;19 (4 ). 10.3201/eid1904.120117 23628282 74 Adams B , Kapan DD . Man bites mosquito: Understanding the contribution of human movement to vector-borne disease dynamics . PLoS One . 2009 ;4 (8 ):e6763 10.1371/journal.pone.0006763 19707544 75 Cosner C , Beier JC , Cantrell RS , Impoinvil D , Kapitanski L , Potts MD , et al The effects of human movement on the persistence of vector-borne diseases . J Theor Biol . 2009 ;258 (4 ):550 –60 . 10.1016/j.jtbi.2009.02.016 19265711 76 Stoddard ST , Morrison AC , Vazquez-Prokopec GM , Paz Soldan V , Kochel TJ , Kitron U , et al The role of human movement in the transmission of vector-borne pathogens . PLoS Negl Trop Dis . 2009 ;3 (7 ):e481 10.1371/journal.pntd.0000481 19621090 77 Mondini A , Chiaravalloti-Neto F . Spatial correlation of incidence of dengue with socioeconomic, demographic and environmental variables in a Brazilian city . Sci Total Environ . 2008 ;393 (2-3 ):241 –8 . 10.1016/j.scitotenv.2008.01.010 18262225 78 Honório NA , Nogueira RMR , Codeço CT , Carvalho MS , Cruz OG , Magalhães MdAFM , et al Spatial evaluation and modeling of Dengue seroprevalence and vector density in Rio de Janeiro, Brazil . PLoS Negl Trop Dis . 2009 ;3 (11 ):e545 10.1371/journal.pntd.0000545 19901983 79 Hu W , Clements A , Williams G , Tong S , Mengersen K . Spatial patterns and socioecological drivers of dengue fever transmission in Queensland, Australia . Environ Health Perspect . 2012 ;120 (2 ):260 –6 . 10.1289/ehp.1003270 22015625 80 de Mattos Almeida MC , Caiaffa WT , Assunção RM , Proietti FA . Spatial vulnerability to dengue in a Brazilian urban area during a 7-year surveillance . J Urban Health . 2007 ;84 (3 ):334 –45 . 10.1007/s11524-006-9154-2 17243024