Predicting re-emergence times of dengue epidemics at low reproductive numbers: DENV1 in Rio de Janeiro, 1986-1990

5 Predicting arbovirus re-emergence remains challenging in regions with limited off-6 season transmission and intermittent epidemics. Current mathematical models treat the 7 depletion and replenishment of susceptible (non-immune) hosts as the principal drivers 8 of re-emergence, based on established understanding of highly transmissible childhood 9 diseases with frequent epidemics. We extend an analytical approach to determine the 10 number of ‘skip’ years preceding re-emergence for diseases with continuous seasonal 11 transmission, population growth and under-reporting. Re-emergence times are shown 12 to be highly sensitive to small changes in low R0 (secondary cases produced from a 13 primary infection in a fully susceptible population). We then fit a stochastic SIR 14 (Susceptible-Infected-Recovered) model to observed case data for the emergence of 15 dengue serotype DENV1 in Rio de Janeiro. This aggregated city-level model 16 substantially over-estimates observed re-emergence times either in terms of skips or 17 outbreak probability under forward simulation. The inability of susceptible depletion and 18 replenishment to explain re-emergence under ‘well-mixed’ conditions at a city-wide 19 scale demonstrates a key limitation of SIR aggregated models including those applied 20 to other arboviruses. The predictive uncertainty and high skip sensitivity to 21 epidemiological parameters suggest a need to investigate the relevant spatial scales of 22 susceptible depletion and the scaling of microscale transmission dynamics to formulate 23 simpler models that apply at coarse resolutions. 24 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint NOTE: This preprint reports new research that has not been certified by peer review and should not be used to guide clinical practice. 2

25 Epidemics of arboviruses such as dengue (1), Zika (2, 3), and chikungunya (4) 26

in substantial global morbidity. Over the past decade, invasions of several 27 arboviruses have triggered large outbreaks in the Western Hemisphere. In Brazil, these 28 invasions include dengue serotype DENV4 in 2012 (5) as well as Zika (2, 6) and 29 chikungunya (7) between 2014-2016. Predicting and understanding the re-emergence 30 of arboviruses after these invasions has important consequences for epidemic 31 preparedness, particularly in regions where climate factors limit mosquito transmission 32 in the off-season. These regions typically exhibit highly intermittent seasonal 33 epidemics, lasting one to three years with long periods of no, or low, reported cases in 34 between, and low mean reproductive numbers (the number of secondary cases arising 35 from each primary case in a completely susceptible population, R0) (5, 8-10). Several 36 proposed explanations include the depletion of susceptible individuals following initial 37 epidemics (11) and the time required for their replenishment via population growth (12), 38 inter-annual variation in climate (13-17), and antigenic interactions between strains of 39 different serotypes (18-21). These temporal patterns contrast with the recurrent 40 seasonal outbreaks observed in childhood diseases with high reproductive numbers, 41 whose extensive study has provided the basis for our theoretical understanding of SIR 42 (Susceptible-Infected-Recovered) dynamics in infections that confer lifelong or lasting 43 immune protection (22-29). 44 Statistical models of dengue transmission that take into account climate 45 dependencies can be used to make short-term re-emergence forecasts on the order of 46 4 months (30) or 16 weeks (15). Many epidemiological models that predict the re-47 emergence of arboviruses such as Zika (11, 31) on longer time-scales of a year (11) or 48 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 3 several decades (31) rely however on compartmental formulations such as SIR-type 49 approaches (11) or Ross-McDonald equations that explicitly incorporate vector 50 transmission (31). Both formulations assume transmission between any two individuals 51 in the population (‘well-mixed’ conditions), typically at aggregated spatial scales. These 52 process-based formulations, for example those recently applied to Zika, represent the 53 acquisition of immunity in the population and its loss via demographic growth and 54 turnover. These models do take into account seasonality of transmission and spatial 55 heterogeneity in the intensity of transmission due to climate at coarse resolutions (at 56 large city, state, or country-level scales). Nevertheless, the replenishment of a well-57 mixed susceptible population is frequently assumed to be the principal driver 58 determining when the disease will re-emerge given a particular seasonal pattern for R0 59 at a particular location(31). Stochasticity can also play an important role in long-term 60 models of re-emergence (31). Variation in reporting rates of arboviruses between 61 locations (32) can add further complexity. 62 Although childhood diseases with high reproductive numbers display different 63 dynamics from emergent arboviruses (22-26), their compartmental models share a 64 basic SIR structure given the acquisition of long-term immunity after infection. The 65 resulting depletion and replenishment of the susceptible population is known to clearly 66 drive inter-annual variability and re-emergence in the former (25, 27, 28). In particular, 67 recent theory (29) has derived analytical expressions for the number of “skip” years for 68 a measles-like disease in the pre-vaccine era, where “skips” are defined as seasons 69 when transmission occurs but does not cause susceptible depletion. In other words, 70 although the number of infections increases in such seasons, it is not large enough to 71 offset the growth in the susceptible population due to demography. The resulting 72 expressions specifically provide a threshold condition for the number of skips expected 73 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 4 following an initial invasion as a function of R0. Their derivation did not include under-74 reporting and assumed a closed-population SIR model with ‘school-term’ seasonality, 75 alternating two different rates for low and high transmission. 76 We examine in this work whether replenishment of susceptible individuals under 77 the typical ‘well-mixed’ assumption explains dengue (DENV1) re-emergence at the 78 whole-city aggregated level. We specifically address the uncertainty inherent in such 79 predictions at the low reproductive numbers characteristic of arboviruses, not previously 80 considered in applications of the analytical approach. To this end, we first extend the 81 threshold derivation to take into account population growth, continuous (sinusoidal) 82 seasonality, and under-reporting of cases. We then fit a stochastic SIR model to 83 observed monthly dengue case counts from the DENV1 invasion in Rio de Janeiro, 84 Brazil from 1986-1988 (8, 10, 33) and numerically predict expected times to re-85 emergence. We describe high uncertainty in re-emergence times for these seasonal, 86 low transmission regions, and show the insufficiency of susceptible replenishment in a 87 simple SIR model to explain the short periods observed in DENV1 re-emergence. We 88 discuss possible explanations and the need for model formulations that would scale to 89 coarse spatial resolutions. 90 91

92 We start with the analytical approach for a seasonally forced SIR system with in-93 termittent outbreaks and population turnover, to consider general features of re-emer-94 gence at low R0. In such a system, the onset of the off-season can bring an end to an 95 initial outbreak, and the replenishment of susceptible individuals due to births and popu-96 lation turnover can be a major determinant of recurrence times. Let S represent the 97 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 5 number of susceptible individuals in a population and s0, the fraction of the population 98 still susceptible at the end of an initial epidemic, t0, when a prediction for the time to 99 the next outbreak will be made. If there are enough susceptible individuals left in the 100 population (i.e. if s0 is large), another outbreak will occur in the following year once the 101 on-season resumes. However, if the initial outbreak was very large, s0 may be too small, 102 and the outbreak may “skip” one or more years. A skip year is defined as a year in 103 which the susceptible population does not decrease, whether or not infections increase. 104 The smaller the fraction of the susceptible population at the time of prediction (s0), the 105 longer it will take for the susceptible population to replenish, and the larger the number 106 of skips that will occur. Previous theory(29) allows prediction of the number of skips 107 that will occur given s0. Specifically, it demonstrated that s0 must fall below some thresh-108 old sc (k) for k skips to occur. An analytical expression was provided for sc(k) in terms of 109 the reproductive number and population turnover rate for a closed-population SIR model 110 with school-term seasonality (29). The derivation of the threshold presented in (29) re-111 quires the assumption that the transmission rate or reproductive number of the disease 112 is high and that the fraction of the population susceptible at the time of prediction (s0) is 113 small. 114 We extend this approach to take into account population growth and sinusoidal 115 seasonality (which de scribes the transmission rate of dengue more accurately than a 116 discrete high -low representation). Our derivation does not require as suming that the 117 transmission rate or reproductive number are high or that the fraction of the population 118 susceptible at the time of prediction is small. We follow the criteria developed in (29) (see 119 details in (34)), which essentially consider the sign of the logarithm of the ratio between 120 the respective number of infections at two times, t0 and tn>t0. A positive value indicates 121 that an outbreak will still occur at tn; conversely a negative value indicates no outbreak at 122 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 6 that time. By setting the logarithm of this ratio to zero, the threshold sc is obtained (See 123 Section 1 of the Supporting Information for details). 124 The resulting expression for sc(n), the critical fraction of susceptible individuals 125 required at the time of prediction for n or more skip seasons to occur, is 126 s!(𝑛) = 1 + "($%&')(') ! "# ))$* +, (+,*,.,%) (1) 127 where 𝑓(𝜔, 𝛿, 𝑟, 𝑛) = (1 + 𝑒). $ %($%&'))𝜔 δ (𝜔$ + 𝑟$)⁄ − (1 − 𝑒). $ %($%&')) 𝑟1 , R0 is the 128 annual mean of the reproductive number, d, the amplitude of seasonal transmission (as 129 infectious contacts per person per day), 𝜔, the transmission frequency (in days-1) and r, 130 the population growth rate (also in days-1). The full expression for the seasonal 131 transmission rate is given by 𝛽(𝑡) = 𝛽/(1 + 𝛿𝑠𝑖𝑛(𝜔𝑡 + 𝜙)), where f corresponds to the 132 phase (in radians) and b0, to the mean seasonal transmission rate (infectious contacts 133 per person per day). The quantity b0 is related to the annual mean reproductive number 134 R0 via the expression 𝑅/ = 𝛽//𝛾, where g is the recovery rate (in days-1). 135 Figure 1 illustrates the implications of this formula. As before, t0 corresponds to 136 the time of prediction, in practice usually after a large initial epidemic or invasion. Like-137 wise, s0 represents the fraction of the population susceptible at the time of prediction. In-138 tuitively, the smaller the fraction of the population susceptible at the time of prediction 139 (s0), the longer it will take for the susceptible population to replenish, and the larger the 140 number of skips that will occur. In practice, as we will illustrate below, values of s0 can 141 be computed from surveillance data provided one has an estimate of the reporting rate. 142 For n skips to occur, the fraction of the population susceptible at the time of pre-143 diction (s0) must fall below the susceptibility threshold sc(n). Figure 1A shows that the 144 larger the number of skips n one is considering, the smaller the threshold sc(n) that s0 145 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 7 must fall below for at least n skips to occur. Let nc denote the critical skip number cor-146 responding to the number of skips expected at the time of prediction (t0). We use the 147 fraction of the population susceptible at the time of prediction (s0) and identify the maxi-148 mum value of n for which s0 is smaller than sc(n). In the example shown in Fig. 1A, this 149 fraction s0 = 0.7 is smaller than sc(n = 6) and bigger than sc(n = 7), which means nc = 6. 150 We therefore expect six years of skips followed by re-emergence in the seventh year. 151 Formally, for a given value of s0 at the end of the transmission season, we define the 152 critical skip number nc as the value of n for which sc(nc) > s0 > sc(nc + 1). 153 154 155 Fig 1. A) Graphical illustration of how the expected number of skips ( nc) is 156 calculated. The black dots represent the threshold fraction of the population susceptible 157 at the time of prediction required for n skips to occur (sc(n)). The plot shows (sc(n)) as a 158 function of n (the number of skips) obtained from Equation 1 with seasonality amplitude 159 d=0.2 ( contacts per person per day) and reproductive number R0=1.4. In this example, 160 the red line represents the fraction of the population susceptible at the time of prediction 161 (s0). If s0 is smaller than sc(n), at least n skips will occur. To find the expected number of 162 skips (nc), we identify the largest number of skips n such that s0 is smaller than the 163 susceptibility threshold required for those skips sc(n). In this example, t he red line 164 intersects the sc(n) curve between sc(n=6) and sc(n=7). Therefore, a critical skip number 165 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 8 of nc=6 is obtained. B) and C) The critical skip value n c as a function of R0 for (B) 166 different values of the amplitude of seasonal transmission d with s0=0.7 and (C) different 167 values of the fraction of the population susceptible at the time of prediction (s0) with 168 d=0.70. In all three panels, the frequency of transmission w, the population turnover rate 169 µ, and population growth rate r are fixed at respective values w = (2p/365) day-1 170 corresponding to a n annual periodicity, µ= 1/ (74.46*365)) day-1 corresponding to an 171 average lifespan of ~75 years, and r=1.55µ day-1 consistent with the growth of the city of 172 Rio de Janeiro. These values were chosen for the purpose of illustration, based on the 173 inverse of the average life expecta ncy in Brazil in 2012 according to the 2010 census 174 (35), and the interpolation of population estimates for the resident population of the 175 municipality of Rio de Janeiro from the 199 1 (36) and 2000 (37) censuses assuming 176 exponential growth. 177 178 With this general approach at hand, w e explored the effects of the reproductive 179 number R0, amplitude of seasonal transmission d and fraction of the population 180 susceptible at time of prediction s0, on the critical number of skips nc (Figure 1 Panels B 181 and C). Consideration of both the variation of the reproductive number R0 and fraction of 182 the population susceptible at time of prediction s0 is relevant here. Different combinations 183 of transmission ra te (β0) and duration of the infection (1/ γ) can yield the same R0 but 184 different fractions of the population susceptible at the time of prediction (Supplemental 185 Figure S16). Importantly, Fig. 1 p anels B and C show that the time to re -emergence is 186 very sensitive to R0. A singularity is observed as R0 approaches 1 where the expected 187 number of skips goes to infinity. The approach to that singularity can be very steep, 188 meaning that small cha nges in R0 can result in large increases in the expected re -189 emergence time. The obtained values of nc are not as sensitive to the amplitude of 190 seasonal transmission (Fig. 1 Panel B) but are sensitive to the fraction of the population 191 susceptible at the time of prediction (Fig. 1 Panel C). The shift of the curve in Fig 1 Panel 192 C for small values of the fraction of the population susceptible at time of prediction s0 193 means that, for a given R0, more time is required to replenish the susceptible population 194 and therefore to observe a re-emergence. 195 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 9 We next apply this approach to the surveillance data from the 1986 invasion of 196 DENV1 in Rio de Janeiro (Figure 2). The initial DENV1 invasion in Rio de Janeiro is an 197 ideal initial test case for this technique given the lack of widespread prior immunity from 198 to prior dengue epidemics, vaccination campaigns, cross-immunity from other disease 199 outbreaks. Specifically, the 1986 invasion occurred prior to the development of dengue 200 vaccines. The outbreak was the first dengue invasion in the area since the initial eradi-201 cation of the Aedes aegypti mosquito in Brazil in the 1950s (38-41) following a sus-202 tained intervention program that began in the 1930s and 1940s in Rio de Janeiro and 203 other cities (39). Cross-immunity from yellow fever vaccination appears to be very lim-204 ited (42). Given the young age distribution of the population in 1986 (43), most individu-205 als were not alive during the period when mass yellow fever vaccination or prior dengue 206 epidemics occurred. 207 We let our tim e of prediction t0 be equal to September 1, 1987, corresponding to 208 the end of the initial DENV1 invasion (see panel A of Figure 2). In panel B of Figure 2, we 209 evaluate the number of expected skips expected in Rio de Janeiro, nc, on the basis of a 210 range of R0 values from 1.18 to 2.02 from the literature (44, 45). The critical susceptibility 211 threshold for n skips to occur ( sc(n)) is calculated using Equation 1 with an annual 212 seasonality, a population growth rate interpolated fr om the census (see Materials and 213

section), and d =0.7 (44). The fraction of the population susceptible at the time 214 of the prediction (s0) is estimated as the difference between the total population N0 (total 215 population N at (t0=Sep. 1987)) and the total number of people infected between the start 216 of the invasion and the time of prediction ( September 1, 1987). The total number of 217 infected people during the outbreak is computed by summing the ratio between the 218 observed monthly cases and the reporting rate for DENV1 in the city. Literature estimates 219 from serology during the DENV1 invasion in Rio de Janeiro indicate a reporting rate of 220 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 10 around 3% (33) which we use and fix for this analysis. For comparison purposes, we also 221 include the number of skips expected under a higher reporting rate of 10%. These curves 222 show that the expected re-emergence could be very sensitive to small variation in R0 and 223 ρ, two quantities that are difficult to estimate with precision in the absence of serology. In 224 particular, assuming a reporting rate of 3%, a reproductive number of 1.2 with 20% 225 uncertainty can yield large changes in the expected re-emergence time. We highlight the 226 potential sensitivity of the expected number of skips to the repo rting rate as well to 227 illustrate the importance of uncertainty in th is parameter in cities or epidemics where its 228 value is unknown. 229 230 Fig 2. (A) Observed dengue case data. Monthly reported dengue cases in the city of 231 Rio de Janeiro, Brazil from April 1986-1995. The grey shaded region denotes 232 observations that were included in the fitting of the stochastic model from May 1, 1986 to 233 July 1, 1988 inclusive. Serotype DENV1 re-emerged in 1990. DENV2 was first detected 234 in the state of Rio de Janeiro in 1990 but did not become dominant until 1991 (8, 9). Both 235 co-circulated afterwards. We focus on the invasion of DENV1 from 1986 -1987 and its 236 initial re-emergence in DENV1 in 1990 using a single serotype transmission model. This 237 allows us to evaluate this transmission model in a region where only one serotype was 238 circulating, where cross-immunity could not easily be invoked to explain the absence or 239 reduction of dengue in a given year. (B) Deterministic critical number of DENV1 skips 240 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 11 nc for Rio de Janeiro from September 1988 . Expected number of skips nc with 241 amplitude of seasonal transmission d=0.7 and the fraction of the population susceptible 242 after the first DENV1 invasion as of September 1, 1987 (s0) calculated from the data (A). 243 We use a reporting rate r of 3% when calculating s0 , consistent with serological 244 estimates from the literature (33). For comparison purposes, we also include the expected 245 number of skips nc assuming a reporting rate of 10%. 246 247 248 Replenishment of Susceptible Individuals is Insufficient to Explain Re-Emergence 249 To obtain more precise bounds for the reporting rate and R0 and to determine if 250 the depletion and replenishment of susceptible individuals could explain the rapid re-251 emergence of dengue in Rio de Janeiro, we fit a stochastic aggregate SIR model to 252 case data from the DENV1 invasion from 1986-1988. The stochastic SIR model 253 assumes that the underlying deterministic transmission rate varies seasonally as a 254 sinusoidal function with annual mean b0, seasonal transmission amplitude d, frequency 255 w (equal to 2p/365), and phase f. The model takes into account demographic 256 stochasticity, environmental stochasticity in the transmission rate, and measurement 257 error due to under-reporting and variation in reporting of cases (See Materials and 258

and the Supporting Information). Panels A,B, and C of Figure 3 show the 259 likelihood profile of the annual mean transmission rate, b0, the amplitude of seasonal 260 transmission d, and the reporting rate r, respectively. In particular, our estimate of the 261 reporting rate matches that from serology in the literature (Panel C). 262 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 12 263 Fig 3. A-C) Selected parameter profiles for the stochastic model. Profiles of the mean 264 annual transmission rate b0 (A), seasonal transmission amplitude d (B), and reporting rate 265 r (C). The red curve is a polynomial fit to the subset of the profile points shown on the 266 figure. The single dashed grey horizontal line represents the likelihood value 2 log 267 likelihood units below the maximum likelihood estimate. This line provides an estimate of 268 confidence intervals for the given parameter. The grey vertical line denotes the parameter 269 value of the maximum likelihood estimate . The maximum likelihood estimate for the 270 reporting rate in panel C is very close to the literature value obtained from serology 271 (approximately 3 percent). (33). 272 273 Overall, the model is able to capture key dynamics of the DENV1 invasion 274 including the two peaks of incidence in 1986 and 1987 and the subsequent reduction of 275 transmission in 1988. This is shown by comparing the trajectories for an ensemble of 276 simulations with the fitted model to the observed values of cases (Fig. 4). Estimated 277 values for the transmission rate indicate a low value for R0 (Figure 4 Panel C). Both of 278 these conclusions generally hold even if one takes into account uncertainty in 279 parameter estimates by examining all parameter combinations with log likelihoods 280 within 2 log likelihood units of the maximum likelihood estimate (the grey region in 281 Figure 4 Panel C as well as Supplemental Figure S1), although some parameter 282 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 13 combinations (not the maximum likelihood estimate) have substantial process noise 283 (Supplemental Figure S1). 284 285 286 287 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 14 288 Fig 4. A-B) Comparison of simulated values with the fitted model and observed 289 data on a log (A) and regular (B) scale. Observed monthly cases from April 1986 to 290 June 1988 are shown in blue. Median values from 100 simulations with the maximum 291 likelihood parameter combination are shown in red. The shaded red region denotes the 292 2.5% and 97.5%th quantile boundaries from those simulations. C) Estimates for R0(t). 293 The black line denotes the trajectory of R0(t) for the maximum likelihood estimate. The 294 shaded grey region represents the 2.5% and 97.5%th quantile boundaries for 295 trajectories from all parameter combinations within 2 log likelihood units of the maximum 296 likelihood estimate. Each parameter combination has only one seasonal trajectory for 297 R0(t) since R0(t) is a deterministic quantity. R0(t) for all parameter estimates ranges from 298 1.79-2.09 in the on season to 0.31-0.52 in the off-season. 299 300 We now apply the obtained parameter estimates from the fitted model to 301 address the expected re-emergence time on the basis of, first, the analytical expression 302 for the skip calculation (Equation 1), and then the stochastic simulations of the fitted 303 model. The parameter estimates used here are those for the reporting rate r, the 304 reproductive number R0, and the amplitude of seasonal transmission d from all 305 combinations within 2 log likelihood units of the MLE. The expected number of skips 306 following the DENV1 invasion in 1986-1988 is considerably higher than the observed 2 307 years. Depending on the parameter combination used, we obtain anywhere from 27 to 308 100 skips (Panel A of Figure 5). Even the fastest estimated return from the skip analysis 309 (27 years) is much slower than the observed re-emergence time. 310 Forward simulation of the stochastic model likewise does not predict the rapid re-311 emergence of DENV1 (Panel B of Figure 5). Under a pulse of 20 infected individuals 312 arriving per day, there was a low probability of re-emergence for parameter 313 combinations with low process noise (Panel B of Figure 5). Only parameter 314 combinations with high amounts of process noise (which have limited predictive value) 315 had a non-zero emergence probability. We consider alternate pulse rates in 316 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 15 Supplemental Figure S14. Re-emergence probabilities under forward simulation of the 317 stochastic model thus corroborated the deterministic skip findings. The depletion of 318 susceptible individuals from 1986-1988 and their replenishment via population growth 319 from 1989-1990 under an aggregate SIR model was unable to explain the rapid re-320 emergence of DENV1 in 1990. 321 322 323 Fig 5. A) Expected number of skips (nc) calculated using parameters obtained 324 from the fitted stochastic model. The open circles show the expected number of 325 skips nc from Equation 1 using parameters and the fraction of the population susceptible 326 after the initial DENV1 invasion (s0) estimated from the fitted stochastic model. Each 327 circle corresponds to one parameter combination, and we included here all parameter 328 combinations for the fitted model with a seasonal transmission amplitude (d) of 0.7 ( 329 contacts per person per day) and a likelihood value within two log-likelihood units of the 330 maximum likelihood estimate (MLE). See Figure S15 for expected skips from 331 parameter combinations with different values of d, and Figure S10 for parameter 332 combinations from the profile of the recovery rate, g. For comparison purposes, the 333 black line shows the expected number of skips for the deterministic skip calculation from 334 panel B of Figure 2 with the reporting rate r fixed at the literature value of 3%. B) 335 Probability of epidemic in 1990 under forward stochastic simulation of fitted 336 model. The fitted stochastic model was simulated forward in time from 1986-1990 with 337 population growth. A pulse of 20 infected individuals were assumed to arrive each day 338 in January 1990. Each parameter combination within 2 log likelihood units of the 339 maximum likelihood estimate was simulated 100 times. The re-emergence probability 340 was calculated by determining the number of simulations in which the susceptible 341 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 16 population decreased in 1990. The plot shows re-emergence probability as a function of 342 the process noise intensity sP. Each point represents a single parameter combination. 343 The maximum likelihood estimate parameter combination is circled in red. 344 345 Sensitivity Analysis: 346 To examine the robustness of our findings to adding an incubation period or 347 altering the form of seasonality, we conducted a sensitivity analysis by considering both 348 SIR and SEIR models with spline seasonality. The results are presented and discussed 349 in the Supporting Information and show that our conclusions remain unchanged. (See 350 the Supporting Information including Supplemental Figures S2-S7 and Supplemental 351 Tables ST2 and ST3). 352 Comparison with Vector Model and literature R0 353 The fitted stochastic SIR model uses a cosine function as a simplification to rep-354 resent the seasonal forcing that would be created by climate variation (temperature 355 (46)) via the changes in infected mosquitos. To evaluate whether this simplification is 356 realistic, we take two approaches. The first one compares the mean seasonal R0 result-357 ing from our model to values of this reproductive number directly estimated from time 358 series data in the literature for DENV1 and DENV4 in Rio de Janeiro from 2010-2016. 359 There is a close match between these very different ways to estimate R0, and in particu-360 lar the shape of the seasonality produced by our model is realistic (Supplemental Figure 361 S18). 362 The second approach considers a simple temperature-driven vector model. To 363 this end, we initially show that the seasonal variation in temperature in Rio de Janeiro 364 can be approximated via a cosine function (Panel A of Supplemental Figure S19 and 365 use this approximation to drive a transmission rate that includes the vector explicitly. 366 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 17 To obtain an expression for the seasonal transmission rate we consider an ex-367 plicit mosquito model with compartments for infectious and susceptible mosquitoes in 368 which a number of parameters depend on temperature (T) (see Section 4 of the Sup-369 porting Information). By assuming fast dynamics of the mosquito (so that levels of in-370 fection in the mosquito population quickly equilibrate to the dynamics of infection in the 371 human population), we derive the following expression for the effective transmission 372 rate in the mosquito-human model in terms of the biting rate a(T), probability of human 373 infection given an infectious bite b(T), probability of mosquito infection given biting of an 374 infectious human pMI(T), adult mosquito mortality rate µM, carrying capacity K of the 375 mosquito population, human population size N, and mosquito demographic function 376 g(T): 377 𝛽0,, = 𝑎(𝑇)$𝑏(𝑇)(𝑝𝑀𝐼(𝑇)) 𝜇1 𝐾 𝑁 D1 − 𝜇1 𝑔(𝑇)F (2) 378 379 The function g(T) is the product of the eggs laid per female mosquito per gono-380 trophic cycle, the mosquito egg-to-adult survival probability, and the mosquito egg-adult 381 development rate divided by the adult mosquito mortality rate µM. The temperature-de-382 pendence of these components was borrowed from the literature (47, 48) (see Support-383 ing Information Section 4 for details). 384 Under the fast dynamics assumption, this effective transmission rate beff is an 385 implicit representation of the force of infection inflicted on humans by the vectors of the 386 coupled human-vector model. When re-scaled between 0 and 1, beff corresponds 387 closely with bMLE, the transmission rate from the fitted stochastic SIR cosine model 388 (Panel B of Supplemental Figure S19). This close correspondence indicates that the 389 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 18 SIR cosine model is able to capture the shape of the seasonality of DENV1 in Rio de 390 Janeiro. 391 392

393 We developed two lines of evidence regarding the uncertainty and predictability 394 of the time to re-emergence for diseases with low reproductive numbers, on the basis of 395 a seasonally forced SIR model under the ‘well-mixed’ assumption at aggregated, city-396 wide, scales. We showed with an analytical approach that the time to re-emergence 397 (expressed as the number of “skip” years) was highly sensitive to small changes in R0 398 and the fraction of the population still susceptible s0 at the time of prediction (e.g. at 399 the end of the initial outbreak). This sensitivity applies to dengue in Rio de Janeiro 400 where re-emergence times can vary on the order of decades based on literature 401 parameters. This uncertainty contrasts with previous applications of this analytical 402 approach to SIR dynamics in childhood diseases such as measles with much higher R0 403 values where accurate predictions of much shorter skip times have been made (29). We 404 also showed with a stochastic SIR model with seasonal transmission fit to DENV1 405 observed case data for Rio de Janeiro from 1986-1988 that susceptible depletion and 406 replenishment are insufficient to explain dengue re-emergence. The fitted model failed 407 to predict by far the re-emergence of DENV1 in 1990 in terms of either the number of 408 skips expected or the outbreak probability under forward simulation. 409 Transmission parameters like R0 are generally defined with respect to a particular 410 model. Given that we aggregated cases at the city level and used a short time series, 411 care should be taken in interpreting parameter values. Nevertheless, fitted transmission 412 parameters correspond well with literature values and exhibit well-defined confidence 413 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 19 intervals. Estimates of the reporting rate in particular closely match the 3% value (8) 414 obtained via a serological study conducted during the 1986 invasion (8, 33). Reporting 415 rates during the onset of an epidemic may be much lower in regions that have not 416 recently experienced transmission (33, 49) than in those with re-occurring outbreaks 417 and an established surveillance network. This may explain why serological studies of 418 the 1986 invasion (8, 33) and our results, estimate a lower reporting rate for dengue 419 than studies conducted in later years in Brazil (50). Even though different combinations 420 of the transmission rate and duration of infection can yield the same reproductive 421 number, the parameter estimates that compose R0 across all models considered in the 422 sensitivity analysis (which take into account those different combinations) are relatively 423 well-defined. These values are also consistent with the effective reproductive number 424 estimated for local dengue epidemics from 2012-2016 (44) and 1996-2014 (45) taking 425 into account differences in serotype circulation and population size during those 426 periods. 427 More complex model structures are possible and often used for arboviruses that 428 include an explicit representation of the vector. We expect our results to hold as this 429 vector component should largely affect the phase and shape of seasonality in the 430 transmission between human hosts, which we have modeled phenomenologically as a 431 cosine wave. With two typical successive epidemic years from an emergent virus, 432 parameter inference from such short observation period is unlikely to justify a more 433 complex model. Nevertheless, to examine transmission seasonality further, we 434 compared the seasonal R0 resulting from the fitted model to the seasonal R0 directly 435 estimated from time series of cases in the literature (44). We also considered the 436 transmission rate experienced by humans in a simple vector-human model forced by 437 the typical seasonality of temperature in Rio de Janeiro. The shape and timing of the 438 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 20 vector-human model’s transmission rate was comparable to that of the cosine 439 transmission rate we employed. More complex models that do not assume fast 440 dynamics of infection in the vector relative to epidemic spread would likely exhibit a 441 difference relative to our transmission rate, especially a delayed phase, whose 442 consequences should be examined in future work. We posit that this difference would 443 not influence our results on the predictability and uncertainty of re-emergence, since the 444 values of other parameters (such as the length of infection in humans) can compensate 445 for it. 446 Factors that could explain the observed rapid re-emergence include inter-annual 447 climate anomalies, antigenic evolution, or micro-scale spatial heterogeneity in 448 transmission intensity and associated susceptible depletion. Larvae washout following 449 flooding coupled with temperature-driven seasonality in transmission could have 450 temporarily halted the invasion in 1988 and delayed the epidemic in 1989. Widespread 451 flooding was reported in February 1988 (51). Large amounts of rainfall washed away 452 mosquito larvae in lab and field studies (52). High rainfall negatively affected dengue 453 transmission in Singapore (53, 54) and India (55). The impact could be compounded in 454 Rio de Janeiro if the high rainfall occurs during the transmission season. If the larvae 455 population has not fully recovered before the start of the off-season, the impact of the 456 rainfall anomaly could extend to the subsequent season. 457 The large amount of process noise observed in the aggregate model would be 458 consistent with this effect, given that the process noise parameter sP represents random 459 variation in the transmission rate due to environmental factors. However, the model’s 460 inherent structure limits its ability to take into account flooding events via sP, since the 461 magnitude of the process noise does not change between years. Incorporating an inter-462 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 21 annual climate driver could provide more accurate re-emergence predictions. The 463 response to rainfall would be nonlinear: positive at low to moderate levels and negative 464 at higher ones. 465 Intra-serotype antigenic evolution from 1986-1990 could also facilitate faster re-466 emergence. Many models focus on inter-serotype variation and assume long-lasting 467 homosubtypic immunity (18, 19, 21). However, the antigenic variation within and across 468 dengue serotypes is comparable (56), and antigenic differences between strains of the 469 same serotype influence overall dengue evolution(57) . Sequences associated with 470 case data were unavailable, making direct analysis challenging. We cannot rule out the 471 possibility that genetic differences between the circulating strains enabled re-infection. A 472 future SIRS-type model (Susceptible-Infected-Recovered-Susceptible) could 473 incorporate this intra-serotype antigenic evolution. 474 Micro-scale spatial heterogeneity in transmission intensity and the effects of 475 human movement between neighborhoods could also explain the rapid re-emergence. 476 Small-scale differences in socioeconomic status and population density between 477 neighborhoods in a large city can result in different relationships between mosquito and 478 human population sizes, resulting in widespread heterogeneity in R0 across 479 neighborhoods (58). Previous studies of mosquito trap data in the city have 480 demonstrated that neighborhoods with differing socioeconomic characteristics have 481 different vector population patterns (46). In fact, schoolchildren from neighborhoods 482 with divergent socioeconomic characteristics had varying levels of seroconversion 483 during the 1986 invasion (33). Human movement between neighborhoods may also 484 influence transmission within (59) and between (60) those neighborhoods, potentially 485 resulting in non-uniform depletion of susceptible populations between highly connected 486 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 22 and isolated areas of a city. Whether arising through the effects of spatial heterogeneity 487 in transmission or intra-city movement, non-uniform levels of herd immunity could 488 enable faster re-emergence. 489 Our findings reveal the uncertainty of re-emergence predictions with the simplest 490 SIR models, those that would be most useful at times of emergent public health threats. 491 Consideration of the above factors in transmission models whose goal is to inform 492 public health over large regions, and to do so soon after, if not during, an emergent 493 outbreak, is clearly a challenge. For example, coarse resolutions are typically used 494 because of the scales at which the observed cases are reported, the scales at which the 495 climate covariates are available, and the difficulties inherent in incorporating microscale 496 variation including connectivity. Our results should motivate further research into the 497 central question of how we can scale microscale heterogeneity to formulate aggregated 498 models that include it implicitly. It should also motivate the related further understanding 499 of how such microscale heterogeneity influences susceptible depletion and 500 replenishment in particular case studies. From such efforts, we should be able to 501 evaluate whether the increasing availability of high-resolution data makes it feasible to 502 parameterize transmission models at higher resolutions, or to inform new model 503 formulations at coarser resolutions. 504 The inability of susceptible depletion and replenishment in a simple seasonal SIR 505 formulation at a large, city-wide scale, to explain DENV1 re-emergence has potential 506 implications for other arboviruses. Recent long-term Zika forecasts (31) assume that 507 susceptible depletion and replenishment brought an end to the 2015-2017 epidemics 508 and will determine when re-emergence occurs. DENV1 and Zika share the same vector 509 and invaded a completely susceptible population (not accounting for pre-existing cross-510 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 23 immunity from dengue). If factors absent from the basic model were key drivers of 511 DENV1 inter-annual variability, it would not be unreasonable to infer that similar types of 512 factors could have played a major role in the Zika dynamics observed from 2015-2017. 513 Zika re-emergence could similarly occur much earlier than expected. 514 With changing temperature patterns due to climate change, cities in Asia, 515 Europe, and the western hemisphere that currently do not have recurrent local 516 transmission may transition in the near future to the kinds of dynamics studied here. Our 517

suggest that estimates should be interpreted in the context of this sensitivity to 518 small changes in the reporting rate and reproductive number. Factors like variation in 519 reporting rates, micro-scale transmission heterogeneity and inter-annual climate drivers 520 that are often ignored in long-term forecasts may thus become critical in determining re-521 emergence times. Overall, the large uncertainty in re-mergence times may be 522 unavoidable for these regions. Improved models are needed together with richer data 523 than currently used, to address the question of the relevant spatial scales of susceptible 524 depletion. 525 526

527 The derivation of the expression for the number of skip years (Equation 1) is 528 included in Section 1 of the Supporting Information. We fitted a stochastic version of the 529 SIR model to observed monthly case counts in Rio de Janeiro from 1986-1988 to 530 estimate parameters needed to apply this expression, and also to separately predict in 531 parallel the time to re-emergence via numerical simulation Expected re-emergence 532 times were then compared for the two approaches. 533 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 24 Data Description 534 We used monthly dengue case estimates in the city of Rio de Janeiro, Brazil 535 from 1986-1990. Cases were reported to the local public health surveillance system (9, 536 10). The case counts did not contain serotype information, but prior studies indicated 537 that the dengue serotype DENV1 invaded the city of Rio de Janeiro in 1986 (10) and 538 was the dominant serotype in circulation in the state of Rio de Janeiro from 1986-1990 539 (8) prior to the arrival of DENV2 in 1990. DENV2 did not become dominant until 1991 540 (9). 541 Basic Model Formulation 542 Because dengue infection confers full immunity to the same serotype, we 543 considered an SIR (Susceptible-Infected-Recovered) model. The deterministic model 544 for the number of individuals in the Susceptible (S), Infected (I), or Recovered (R) class 545 is given by the following system of ordinary differential equations: 546 𝑑𝑆 𝑑𝑡 = 𝑟𝑁 − 𝜆(𝑡)𝑆 − 𝜇2𝑆 (3) 547 𝑑𝐼 𝑑𝑡 = 𝜆(𝑡)𝑆 − 𝛾𝐼 − 𝜇2𝐼 (4) 548 𝑑𝑅 𝑑𝑡 = 𝛾𝐼 − 𝜇2𝑅 (5) 549 550 𝜆(𝑡) = 𝛽(𝑡)( 𝐼 𝑁) (6) 551 552 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 25 553 𝛽(𝑡) = 𝛽/(1 + 𝛿𝑠𝑖𝑛(𝜔𝑡 + 𝜙)) (7) 554 Deaths occur at rate (µH) given by the inverse of the life expectancy of Brazil in 555 2012 (74.49 years(35)). All individuals are born susceptible. The term r represents 556 population growth. The human population growth rate was estimated from census 557 resident population estimates in 1991 (36) and 2000 (37) assuming exponential growth. 558 This rate was used to interpolate the estimated population in 1986 (See Supporting 559 Information Section 2.1.1 for details). 560 The per capita rate at which susceptible individuals become infected was given 561 by the force of infection l(t) (Equation 6). Individuals recovered at per-capita rate g 562 whose inverse is the duration of infection. Estimates of the duration of infection in den-563 gue vary. One analysis estimated that symptoms of dengue infection last 2-7 days fol-564 lowing an incubation period of 4-10 days (61, 62). For our analysis, we fixed the recov-565 ery rate 𝛾 to be 1/17, assuming an exponentially distributed duration of infection with 566 mean of 17 days encapsulating the maximum extent of the combined incubation and 567 symptomatic period in humans. We take into account the possibility that duration of in-568 fection could vary by profiling over the duration of infection in the sensitivity analysis. 569 The short duration of the available time series meant that fitting a formal vector model 570 could prove difficult and could require additional assumptions in terms of which parame-571 ters could be fitted or fixed from existing formulations in the literature. We therefore 572 used an SIR framework in which the infected stage served as a proxy for the exposed 573 and infected human compartments in a vector model of dengue transmission, and we 574 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 26 assumed infection levels in vector rapidly equilibrate with those in humans (as de-575 scribed later when we consider the vector explicitly (See Section 4 of the Supporting In-576 formation). A duration of infection was thus chosen that corresponds to the upper bound 577 of the estimated pre-infectious period (4-10 days) and infectious period (2-7 days) in hu-578 mans (61, 62). We profiled over the duration of infection in the sensitivity analysis to 579 verify that this parameterization is reasonable. 580 This transmission rate b(t) was represented as a cosine function with mean b0, 581 (units of contacts per person per day) and seasonal oscillations of amplitude d (same 582 units as b0) and frequency w, which was assumed to be annual (w = 2p/365) days-1. The 583 annual mean R0 was thus given by: 584 𝑅/ = 𝛽/ 𝛾 + 𝜇 (8) 585 The observed dengue data in Rio de Janeiro consisted of monthly case counts. 586 Serological studies of the DENV1 invasion in Rio de Janeiro also indicated substantial 587 under-reporting (8, 33). Let C represent the true number of monthly cases that would be 588 obtained by summing the number of individuals entering the infected class (I) over the 589 course of a month. For the purposes of the skip analysis, we assume that a fixed 590 fraction r of the true cases C are observed, where r is the reporting rate. 591 The stochastic model is an approximation of the deterministic one used for the 592 skip analysis. For simplicity and given the short time interval, we assumed that there 593 was no population growth over the two and half years of the DENV1 invasion (r =µH) 594 and that births and deaths occurred at rate µH = (1/(74.9*365)), which is equal to the 595 inverse of the average life expectancy in Brazil from the 2010 census (35). However, 596 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 27 population growth is taken into account when simulating forward in time from the fitted 597 stochastic model. We also assumed that there were no recovered individuals at the start 598 of the epidemic, so all other individuals in the population not initially infected were 599 susceptible. We considered time in units of days and used a time step Dt of 1 day. 600 The stochastic model is a discrete-time model with fixed time step Dt and a discrete 601 state space (i.e. the number of people in each compartment S, I, R, and C, at any point 602 in time must be integers). The number of individuals who moved from one compartment 603 to another over the course of each day was calculated via Euler simulation from the 604 deterministic equations (See Supporting Information). Demographic stochasticity was 605 then incorporated into the Euler approximations to obtain integer state variable values 606 after each time step. We specifically assumed that the number of individuals making 607 each state transition was drawn from a binomial distribution with exponentially decaying 608 probability (See Supporting Information). Environmental noise (variation in the 609 transmission rate b(t) due to random environmental variation) was captured via 610 multiplicative gamma white noise in the transmission rate as described by (63, 64). On 611 time step size D t, we multiplied the transmission rate by DG / D t, where DG / D t was 612 drawn from a Gamma distribution with mean 1 and variance sP2 / D t. 613 The measurement model assumed that the observed number of monthly dengue 614 cases (Y(t)) at time t were drawn from a negative binomial distribution with mean equal 615 to the true number of monthly cases C multiplied by a reporting rate r, with dispersion 616 parameter sM. More details of the measurement model can be found in Section 2.4 of 617 the Supporting Information. 618 Fitting the stochastic model 619 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 28 We fitted the transmission parameters (b0 and d), reporting rate (r), process 620 noise parameter (sP), measurement noise parameter (sM), and the number of infected 621 individuals at the start of the outbreak in May 1986 (I0). While the first cases of DENV1 622 were reported in April 1986, we started the model fitting in May 1986 to avoid 623 complications from changes in the reporting rate as the surveillance system was 624 established during the start of the DENV1 invasion. We used in an interpolated initial 625 population size of 5,281,842 for Rio de Janeiro in May 1986.The model was fit using the 626 mif2 method in the R-package pomp. The model fitting method is described further in 627 the Supporting Information and in (65). 628 Calculating expected skips using parameter estimates from stochastic model 629 Following the completion of the Monte Carlo Profiles, a maximum likelihood 630 estimate (MLE) parameter combination was obtained from the Monte Carlo Profiles of 631 the fitted model by selecting the parameter combination with the highest likelihood 632 across all profiles. The table of MLE parameter values is shown in Supplemental Table 633 ST1. All sets of parameter combinations within 2 log-likelihood units of the maximum 634 likelihood estimate (from all profiles) were used for the expected skip calculation. The 635 reporting rate (r), b0, and d value of each parameter combination within 2 log likelihood 636 units of the maximum likelihood estimate were applied to a finer gridded version of the 637 deterministic skip calculation described earlier. A distribution for the number of skips 638 expected in Rio de Janeiro following the DENV1 invasion from 1986-1988 was 639 obtained. 640 Stochastic Simulation 641 We then simulated re-emergence probabilities under the stochastic model. Each 642 parameter combination within 2 log likelihood units of the MLE estimate from the 643 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 29 stochastic fit was simulated again without any immigration from 1986 until 1990 but with 644 population growth. During January 1990, “sparks” of infectious individuals were 645 assumed to have arrived in the city at some fixed rate. There were low but non-zero 646 levels of DENV1 incidence from 1988-1989. We chose to wait until January 1990 before 647 introducing new DENV1 infections to be conservative, as this is when an uptick in 648 DENV1 incidence was first observed. Had we introduced sparks earlier in 1988-1989, 649 we would likely have observed even earlier re-emergence times. We explored rates 650 from 5 to 100 infected individuals per day. This process was repeated 100 times, and 651 the probability of an epidemic occurring in 1990 was calculated. An epidemic 652 occurrence in this situation was defined as a net decrease in the susceptible population 653 over the course of the year (after taking into account population growth), to best match 654 the definition of an epidemic used in the skip analysis. 655 Sensitivity Analysis 656 We assessed how parameter estimates of R0 and r may depend on the model 657 formulation by fitting several more complex SIR-type models to the same data using the 658 fitting procedure described in the Methods section: an SIR Spline Model and SEIR 659 Spline Model. As an additional sensitivity analysis, we profiled over the recovery rate for 660 the SIR Cosine Model (Supplemental Figure S9). For details, see the Supporting 661 Information. 662 Comparison with Vector Model and literature R0 663 For a full description of the explicit coupled human-mosquito model with 664 compartments for infectious and susceptible mosquitoes and comparison of 665 transmission rates between this model and the simpler seasonally forced SIR , see the 666 Supporting Information. 667 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 30 668 Author Contributions 669 RS, VRA, MP designed the experiments. RS and VRA conducted the experi-670 ments. RS and EI developed the stochastic model fitting pipeline. CC provided 671 the data and expertise on dengue in Rio de Janeiro. RS, VRA, and MP drafted 672 the manuscript. All authors contributed to the writing of the manuscript. 673 674

675 The authors would like to thank Aaron King for his advice , and two referees 676 for their insightful comments. This work was completed with resources and 677 support provided by the University of Chicago’s Research Computing Cen-678 ter. 679 680 Data Accessibility 681 Data and code used to fit the stochastic model is accessible via a public 682 GitHub repository (https://github.com/pascualgroup/JRSI_DENV1_skips_Rio 683 ). 684 685 Ethics 686 This article does not present research with ethical considerations. The authors 687 declare that they have no competing interests. 688 689 Funding Statement 690 RS was supported by a National Science Foundation Research Train-691 eeship (# 1735359: NRT-INFEWS: Computational data science to advance 692 research at the energy environment nexus). MP and EI were supported by 693 a collaborative grant from the National Science Foundation’s Division of 694 Mathematical Sciences and the National Institutes of Health (# 1761612: 695 Collaborative Research: Urban Vector-Borne Disease Transmission De-696 mands Advances in Spatiotemporal Statistical Inference). VRA was jointly 697 supported by this grant and by the Mansueto Institute for Urban Innovation 698 through a Mansueto Institute Postdoctoral Fellowship. 699 700 The funders had no role in study design, data collection and analysis, deci-701 sion to publish, or preparation of the manuscript. 702 703 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 31 The content is solely the responsibility of the authors and does not neces-704 sarily represent the official views of the National Institutes of Health or the 705 National Science Foundation. 706 707 708

709 710 1. Bhatt S, Gething PW, Brady OJ, Messina JP, Farlow AW, Moyes CL, et al. The global 711 distribution and burden of dengue. Nature. 2013;496(7446):504. 712 2. Baud D, Gubler DJ, Schaub B, Lanteri MC, Musso D. An update on Zika virus infection. 713 The Lancet. 2017;390(10107):2099-109. 714 3. Petersen LR, Jamieson DJ, Powers AM, Honein MA. Zika Virus. New England Journal 715 of Medicine. 2016;374(16):1552-63. 716 4. Nasci RS. Movement of Chikungunya Virus into the Western Hemisphere. Emerging 717 Infectious Diseases. 2014;20(8):1394-5. 718 5. Nogueira RM, Eppinghaus AL. Dengue virus type 4 arrives in the state of Rio de Janeiro: 719 a challenge for epidemiological surveillance and control. Memórias do Instituto Oswaldo Cruz. 720 2011;106:255-6. 721 6. Zanluca C, de Melo VCA, Mosimann ALP, dos Santos GIV, dos Santos CND, Luz K. 722 First report of autochthonous transmission of Zika virus in Brazil. Memórias do Instituto 723 Oswaldo Cruz. 2015;110(4):569-72. 724 7. Nunes MRT, Faria NR, de Vasconcelos JM, Golding N, Kraemer MUG, de Oliveira LF, 725 et al. Emergence and potential for spread of Chikungunya virus in Brazil. BMC Medicine. 726 2015;13(1):102. 727 8. Nogueira RMR, Miagostovich MP, Schatzmayr HG, Santos FBd, Araújo ESMd, Filippis 728 AMBd, et al. Dengue in the State of Rio de Janeiro, Brazil, 1986-1998. Memórias do Instituto 729 Oswaldo Cruz. 1999;94:297-304. 730 9. Nogueira RMR, Miagostovich MP, Lampe E, Souza RW, Zagne SMO, Schatzmayr HG. 731 Dengue Epidemic in the State of Rio de Janeiro, Brazil, 1990-1: Co-Circulation of Dengue 1 and 732 Dengue 2 Serotypes. Epidemiol Infect. 1993;111(1):163-70. 733 10. Miagostovich MP, Nogueira RM, Cavalcanti S, Marzochi KB, Schatzmayr HG. Dengue 734 epidemic in the state of Rio de Janeiro, Brazil: virological and epidemiological aspects. Revista 735 do Instituto de Medicina Tropical de São Paulo. 1993;35(2):149-54. 736 11. O'Reilly KM, Lowe R, Edmunds WJ, Mayaud P, Kucharski A, Eggo RM, et al. 737 Projecting the end of the Zika virus epidemic in Latin America: a modelling analysis. BMC 738 medicine. 2018;16(1):180-. 739 12. Struchiner CJ, Rocklöv J, Wilder-Smith A, Massad E. Increasing Dengue Incidence in 740 Singapore over the Past 40 Years: Population Growth, Climate and Mobility. PLOS ONE. 741 2015;10(8):e0136286. 742 13. Stewart-Ibarra AM, Lowe R. Climate and non-climate drivers of dengue epidemics in 743 southern coastal ecuador. Am J Trop Med Hyg. 2013;88(5):971-81. 744 14. Vincenti-Gonzalez MF, Tami A, Lizarazo EF, Grillet ME. ENSO-driven climate 745 variability promotes periodic major outbreaks of dengue in Venezuela. Scientific reports. 746 2018;8(1):5727-. 747 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 32 15. Hii YL, Rocklöv J, Ng N, Tang CS, Pang FY, Sauerborn R. Climate variability and 748 increase in intensity and magnitude of dengue incidence in Singapore. Global Health Action. 749 2009;2(1):2036. 750 16. Cazelles B, Chavez M, McMichael AJ, Hales S. Nonstationary Influence of El Niño on 751 the Synchronous Dengue Epidemics in Thailand. PLOS Medicine. 2005;2(4):e106. 752 17. Harris M, Caldwell JM, Mordecai EA. Climate drives spatial variation in Zika epidemics 753 in Latin America. Proceedings of the Royal Society B: Biological Sciences. 754 2019;286(1909):20191578. 755 18. Recker M, Blyuss Konstantin B, Simmons Cameron P, Hien Tran T, Wills B, Farrar J, et 756 al. Immunological serotype interactions and their effect on the epidemiological pattern of 757 dengue. Proceedings of the Royal Society B: Biological Sciences. 2009;276(1667):2541-8. 758 19. Reich NG, Shrestha S, King AA, Rohani P, Lessler J, Kalayanarooj S, et al. Interactions 759 between serotypes of dengue highlight epidemiological impact of cross-immunity. Journal of 760 The Royal Society Interface. 2013;10(86). 761 20. Cummings DAT, Schwartz IB, Billings L, Shaw LB, Burke DS. Dynamic effects of 762 antibody-dependent enhancement on the fitness of viruses. Proceedings of the National Academy 763 of Sciences of the United States of America. 2005;102(42):15259. 764 21. Ferguson N, Anderson R, Gupta S. The effect of antibody-dependent enhancement on the 765 transmission dynamics and persistence of multiple-strain pathogens. Proceedings of the National 766 Academy of Sciences of the United States of America. 1999;96(2):790-4. 767 22. Bjørnstad ON, Finkenstädt BF, Grenfell BT. Dynamics of measles epidemics: estimating 768 scaling of transmission rates using a time series SIR model. Ecological Monographs. 769 2002;72(2):169-84. 770 23. Earn DJD, Rohani P, Bolker BM, Grenfell BT. A Simple Model for Complex Dynamical 771 Transitions in Epidemics. Science. 2000;287(5453):667-70. 772 24. Finkenstädt B, Grenfell B. Empirical determinants of measles metapopulation dynamics 773 in England and Wales. Proc Biol Sci. 1998;265(1392):211-20. 774 25. Graham M, Winter AK, Ferrari M, Grenfell B, Moss WJ, Azman AS, et al. Measles and 775 the canonical path to elimination. Science. 2019;364(6440):584-7. 776 26. Grenfell BT, Bjørnstad ON, Finkenstädt BF. DYNAMICS OF MEASLES EPIDEMICS: 777 SCALING NOISE, DETERMINISM, AND PREDICTABILITY WITH THE TSIR MODEL. 778 Ecological Monographs. 2002;72(2):185-202. 779 27. Domenech de Cellès M, Magpantay FMG, King AA, Rohani P. The impact of past 780 vaccination coverage and immunity on pertussis resurgence. Science Translational Medicine. 781 2018;10(434):eaaj1748. 782 28. Lavine JS, King AA, Bjørnstad ON. Natural immune boosting in pertussis dynamics and 783 the potential for long-term vaccine failure. Proceedings of the National Academy of Sciences. 784 2011;108(17):7259-64. 785 29. Stone L, Olinky R, Huppert A. Seasonal dynamics of recurrent epidemics. Nature. 786 2007;446(7135):533. 787 30. Lowe R, Bailey TC, Stephenson DB, Jupp TE, Graham RJ, Barcellos C, et al. The 788 development of an early warning system for climate-sensitive disease risk with a focus on 789 dengue epidemics in Southeast Brazil. Statistics in medicine. 2012;32(5):864-83. 790 31. Ferguson NM, Cucunubá ZM, Dorigatti I, Nedjati-Gilani GL, Donnelly CA, Basáñez M-791 G, et al. Countering the Zika epidemic in Latin America. Science. 2016;353(6297):353-4. 792 32. Imai N, Dorigatti I, Cauchemez S, Ferguson NM. Estimating Dengue Transmission 793 Intensity from Case-Notification Data from Multiple Countries. PLOS Neglected Tropical 794 Diseases. 2016;10(7):e0004833. 795 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 33 33. Figueiredo LTM, Cavalcante SMB, Simoes MC. Dengue serologic survey of 796 schoolchildren in Rio de Janeiro, Brazil, in 1986 and 1987. 1990. 797 34. Olinky R, Huppert A, Stone L. Seasonal dynamics and thresholds governing recurrent 798 epidemics. Journal of Mathematical Biology. 2008;56(6):827-39. 799 35. (IBGE) BIoGaS. Brazil Demographic Census 2010. . Rio de Janeiro, Brazil: Brazilian 800 Institute of Geography and Statistics (IBGE). 2012. 801 36. (IBGE) BIoGaS. "Censo Demographico- 1991-Rio de Janeiro". 1991;"Censo 802 demográfico : 1991 : resultados do universo relativos as características da população e dos 803 domicílios"(Table 1.4: "População residente, por grupos de idade, segundo ti lolesorregiàes, as 804 Microrregiões, os Municípios,os Distritos e o sexo"):32-41. 805 37. (IBGE) BIoGaS. Censo Demográfico. 2000(População residente, sexo e situação do 806 domicílio; Total column). 807 38. Dick OB, San Martín JL, Montoya RH, del Diego J, Zambrano B, Dayan GH. The 808 history of dengue outbreaks in the Americas. The American journal of tropical medicine and 809 hygiene. 2012;87(4):584-93. 810 39. Severo OP. Eradication of the Aedes aegypti mosquito from the Americas. 1955. 811 40. Löwy I. Leaking Containers: Success and Failure in Controlling the Mosquito Aedes 812 aegypti in Brazil. Am J Public Health. 2017;107(4):517-24. 813 41. PAHO. The feasibility of erradicating Aedes aegypti in the Americas. Pan Am J Public 814 Health. 1997;1. 815 42. Souza NCSe, Félix AC, de Paula AV, Levi JE, Pannuti CS, Romano CM. Evaluation of 816 serological cross-reactivity between yellow fever and other flaviviruses. International Journal of 817 Infectious Diseases. 2019;81:4-5. 818 43. de Carvalho JAM. Demographic dynamics in Brazil: recent trends and perspectives. 819 Brazilian journal of population studies. 1997;1(1):5-23. 820 44. Codeço CT, Villela DA, Coelho FC. Estimating the effective reproduction number of 821 dengue considering temperature-dependent generation intervals. Epidemics. 2018;25:101-11. 822 45. Coelho FC, Carvalho LMd. Estimating the Attack Ratio of Dengue Epidemics under 823 Time-varying Force of Infection using Aggregated Notification Data. Scientific Reports. 824 2015;5:18455. 825 46. Honório NA, Codeço CT, Alves FC, Magalhães MAFM, Lourenço-de-Oliveira R. 826 Temporal Distribution of Aedes aegypti in Different Districts of Rio De Janeiro, Brazil, 827 Measured by Two Types of Traps. Journal of Medical Entomology. 2009;46(5):1001-14. 828 47. Mordecai EA, Cohen JM, Evans MV, Gudapati P, Johnson LR, Lippi CA, et al. 829 Detecting the impact of temperature on transmission of Zika, dengue, and chikungunya using 830 mechanistic models. PLOS Neglected Tropical Diseases. 2017;11(4):e0005568. 831 48. Otero M, Schweigmann N, Solari HG. A Stochastic Spatial Dynamical Model for Aedes 832 Aegypti. Bulletin of Mathematical Biology. 2008;70(5):1297. 833 49. Musa SS, Zhao S, Chan H-S, Jin Z, He D. A mathematical model to study the 2014–2015 834 large-scale dengue epidemics in Kaohsiung and Tainan cities in Taiwan, China. Mathematical 835 biosciences and engineering. 2019. 836 50. Imai N, Dorigatti I, Cauchemez S, Ferguson NM. Estimating dengue transmission 837 intensity from sero-prevalence surveys in multiple countries. PLoS neglected tropical diseases. 838 2015;9(4):e0003719-e. 839 51. Ap. 127 Die in Floods and Mud Slides in Brazil. New York Times. 1988 02/08/ 840 1988 Feb 08. 841 52. Fouque F, Carinci R, Gaborit P, Issaly J, Bicout DJ, Sabatier P. Aedes aegypti survival 842 and dengue transmission patterns in French Guiana. Journal of Vector Ecology. 2006;31(2):390-843 400. 844 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint 34 53. Benedum CM, Seidahmed OME, Eltahir EAB, Markuzon N. Statistical modeling of the 845 effect of rainfall flushing on dengue transmission in Singapore. PLOS Neglected Tropical 846 Diseases. 2018;12(12):e0006935. 847 54. Seidahmed OME, Eltahir EAB. A Sequence of Flushing and Drying of Breeding Habitats 848 of Aedes aegypti (L.) Prior to the Low Dengue Season in Singapore. PLOS Neglected Tropical 849 Diseases. 2016;10(7):e0004842. 850 55. Kakarla SG, Caminade C, Mutheneni SR, Morse AP, Upadhyayula SM, Kadiri MR, et al. 851 Lag effect of climatic variables on dengue burden in India. Epidemiol Infect. 2019;147:e170. 852 56. Katzelnick LC, Fonville JM, Gromowski GD, Arriaga JB, Green A, James SL, et al. 853 Dengue viruses cluster antigenically but not as discrete serotypes. Science. 854 2015;349(6254):1338. 855 57. Bell SM, Katzelnick L, Bedford T. Dengue genetic divergence generates within-serotype 856 antigenic variation, but serotypes dominate evolutionary dynamics. eLife. 2019;8:e42496. 857 58. Romeo-Aznar V, Paul R, Telle O, Pascual M. Mosquito-borne transmission in urban 858 landscapes: the missing link between vector abundance and human density. Proceedings of the 859 Royal Society B: Biological Sciences. 2018;285(1884):20180826. 860 59. Stoddard ST, Forshey BM, Morrison AC, Paz-Soldan VA, Vazquez-Prokopec GM, 861 Astete H, et al. House-to-house human movement drives dengue virus transmission. Proceedings 862 of the National Academy of Sciences. 2013;110(3):994-9. 863 60. Guzzetta G, Marques-Toledo CA, Rosà R, Teixeira M, Merler S. Quantifying the spatial 864 spread of dengue in a non-endemic Brazilian metropolis via transmission chain reconstruction. 865 Nature Communications. 2018;9(1):2837. 866 61. Organization WH. Dengue and severe dengue. Fact sheet Updated April 2017. 2017. 867 62. Organization WH. Dengue and severe dengue. World Health Organization. Regional 868 Office for the Eastern Mediterranean; 2014. 869 63. Bretó C, Ionides EL. Compound markov counting processes and their applications to 870 modeling infinitesimally over-dispersed systems. Stochastic Processes and their Applications. 871 2011;121(11):2571-91. 872 64. He D, Ionides EL, King AA. Plug-and-play inference for disease dynamics: measles in 873 large and small populations as a case study. Journal of the Royal Society Interface. 874 2009;7(43):271-83. 875 65. Ionides EL, Nguyen D, Atchadé Y, Stoev S, King AA. Inference for dynamic and latent 876 variable models via iterated, perturbed Bayes maps. Proceedings of the National Academy of 877 Sciences. 2015;112(3):719-24. 878 879 . CC-BY-NC-ND 4.0 International licenseIt is made available under a is the author/funder, who has granted medRxiv a license to display the preprint in perpetuity. (which was not certified by peer review) The copyright holder for this preprint this version posted June 22, 2020. ; https://doi.org/10.1101/2020.05.02.20074104doi: medRxiv preprint