Periodic synchronization of dengue epidemics in Thailand: the roles played by temperature and immunity

The spatial distribution of dengue and its vectors (spp. Aedes) may be the widest it has ever been, and projections suggest that climate change may allow the expansion to continue. However, the largest impacts of climate change on dengue might be in regions where the pathogen is already endemic. In these areas, the waxing and waning of immunity has a large impact on temporal dynamics of cases of dengue haemorrhagic fever. Here, we use 51 years of data across 72 provinces and characterise spatio-temporal patterns of dengue in Thailand, where dengue has caused almost 1.5 million cases over the last thirty years, and examine the roles played by temperature and dynamics of immunity in giving rise to those patterns. We find that timescales of multiannual oscillations in dengue vary in space and time and uncover an interesting spatial phenomenon: Thailand has experienced multiple, periodic synchronization events. We show that patterns in synchrony of dengue are consistent with those observed in temperature. Applying a temperature-driven dengue model, we explore how dynamics of immunity interact with temperature to produce the observed multiannual dynamics and patterns in synchrony. While multiannual oscillations are readily produced by immunity in absence of multiannual timescales in temperature, synchrony in temperature can synchronise dengue dynamics in different locations. However, at higher mean temperatures and lower seasonal variation, immune dynamics become more predominant, and dengue dynamics become more insensitive to multiannual fluctuations in temperature. These findings can help underpin predictions of disease patterns as global temperatures rise. Author summary

The spatio-temporal dynamics of animal populations, including fluctuations and 2 correlations in amplitudes and phases, has been an important area of research in ecology 3 and the physical sciences for decades. Empirical observations of patterns in population 4 dynamics have propelled theory forward leading to a better understanding of 5 predator-prey interactions [1,2], allee effects [3] and the interactions between 6 deterministic and stochastic elements of systems [4,5]. While a rich literature exists on 7 synchronization of systems and their causes [6], there is little empirical evidence for 8 periodic oscillations or fluctuations in spatial synchrony, particularly in epidemiology, 9 where identification of long-term patterns in infectious disease dynamics could be 10 critical for the success of health interventions. For instance, knowing when pathogen 11 population levels are particularly low concurrently across a region presents improved 12 opportunities for pathogen elimination [7]. Likewise, anticipation of global epidemics 13 may assist in the structured allocation of resources, such as vaccines, across space and 14 time. While regular spatial synchrony in other wildlife species has been described 15 (e.g., [8][9][10][11][12]), how the degree of synchrony may vary over time has received less 16 attention [13][14][15]. Here, we describe regular periodic synchronization in the dynamics of 17 dengue haemorrhagic fever (DHF) in Thailand over a 51-year interval. in which dengue is efficiently transmitted [42][43][44]. However, the impact of changing 50 temperature regimes on dengue in locations where the pathogen is already endemic is 51 less well studied. It is in these areas which, even in projections to 2050 or 2080, 52 comprise a majority of the world's population at risk of dengue, that the impacts of 53 climate change might be greatest [42][43][44]. 54 Here, we examine the spatio-temporal dynamics of dengue in Thailand between 1968 55 and 2018. We characterise multiannual cycles and uncover an interesting spatial 56 phenomenon in which Thailand has experienced periodic synchronizations of dengue 57 incidence. In contrast to previous studies, which tended to focus on either extrinsic or 58 intrinsic drivers to explain observed patterns, we hypothesise that immunity constitutes 59 a strong dynamical filter necessary to understand impacts of temperature on dengue. 60 For this reason, we adapt a mechanistic, temperature-dependent, four-serotype dengue 61 model to disentangle how temperature and dynamics of immunity interact to generate 62 periodic synchronizations. We focus on temporary cross-protection between serotypes 63 because of the clear empirical support, although there are other mechanisms, such as 64 antibody-dependent enhancement, that could conceivably also play a role.  Supplementary information (SI)). The presence of multiannual timescales has been 86 previously observed [27,34], but here, with the benefit of longer time series, we show 87 that these change over time and space.

88
Synchrony in dengue 89 We estimate spatial synchrony in dengue cases using a range of methods; we here focus 90 on one approach, wavelet mean fields (WMFs; [45,46]), but provide details and results 91 for the others in Section "Perspectives on synchrony" in SI. WMFs measure synchrony 92 as a function of both timescale and time; they indicate timescales and time points at 93 which both phases and magnitude of oscillations are consistent (or more synchronous) 94 across provinces. The WMF for dengue haemorrhagic fever cases across Thailand 95 describes a system that appears to fluctuate in and out of synchrony (Fig 2b), results reconstructions of time series using multiannual components only, per province arranged from north (top) to south (bottom). To improve clarity, values for each province were normalised to standard scores (µ = 0, σ = 1) in all panels. Blues (respectively reds) are lower (respectively higher) numbers, and whites correspond to a standard score of zero (the mean for each province). Edge effects in the wavelet transforms may influence results before and after the vertical dashed lines in (b) (see Materials and Methods).
that are consistent with the travelling waves observed across Thailand [27] and 97 Southeast Asia [34]. The synchrony in dengue is statistically significant at all times, 98 meaning that even when synchrony is lower (the whiter areas in Fig 2b), the degree of 99 synchrony is still greater than would be expected in an asynchronous system.

100
Furthermore, moments of greater synchrony take place at different timescales. For 101 example, while the synchrony event of the late 1980s occurs with a timescale of ≈ 2.2 102 years, the one taking place around the year 2000 has a timescale of ≈ 4.1 years. These 103 results are further corroborated using four additional approaches to estimate synchrony 104 (Fig 3); all methods describe a system that appears to fluctuate in and out of synchrony 105 (Section "Perspectives on synchrony" in SI). In all cases, periods of greater synchrony 106 appear to coincide with larger outbreaks (e.g., comparing Fig S6a and c in SI), as also 107 previously observed [34,47].

108
Synchrony in temperature 109 We calculated the temperature WMF across provinces for the same time period as the 110 dengue passive surveillance data . As expected, synchrony in temperature 111 was significant at all points in time. The WMF for temperature describes a system 112 similar to that found in dengue, where temperature fluctuates in and out of synchrony 113 at multiannual timescales (Fig 2a; for results using alternative approaches see Figs S10 114     and S11 in SI). Furthermore, the moments of greater synchrony appear to align, both in 115 timing and timescale, with those of dengue (Fig 2a,b). 116 If the patterns in synchrony in dengue were to be directly linked to patterns in 117 synchrony in temperature, we would expect the magnitudes of the two respective 118 WMFs (in the multiannual time scales, and for all points in time) to correlate positively 119 and significantly. We found a Pearson correlation of 0.21 (P = 0.060; see Material and 120 Methods for significance testing) and a Spearman correlation of 0.20 (P = 0.069). Thus, 121 while correlations were positive, they were not statistically significant. However, see the 122 result for WMFs of dengue data and a mechanistic model based on temperature in 123 'Driving the model with actual temperature time series' below.

124
To address whether the patterns of synchrony in temperature were consistent with 125 those observed in the dengue cases in greater detail, we extended the methods of 126 Sheppard et al. [45,46] to calculate cross-wavelet mean fields (CWMFs). Cross-wavelets 127 quantify the similarity in two time series' (temperature and dengue cases) wavelet power 128 Our analyses describe a system that oscillated in and out of synchrony (Figs 2b, 3). To 137 describe the periodicity of the degree of synchrony itself, we collapsed WMFs to a time 138 series by taking the mean WMF across multiannual timescales for each point in time

139
(red line in Figure 3b). We then normalised the resulting time series to standard scores, 140 and applied WTs.

141
Because these time series of synchrony contain, at most, four distinct cycles (Fig 3), 142 we chose to characterise their periodicity by averaging the wavelet power over time, by 143 cross-protective interactions between serotypes of varying mean lengths of time.

166
The patterns quantified so far describe a system that fluctuates in and out of 167 synchrony. When seeking potential drivers for the observed patterns, these would need 168 to account for both asynchrony, and thus the ability to produce a range of spatially 169 heterogeneous dynamics, and synchrony, during which dynamics are more homogeneous. 170 We ran three sets of simulations to see whether we could better understand how 171 temperature might drive patterns in synchrony, using both real temperature time series 172 and simplified experiments (  locations would be sufficient to produce synchrony in an otherwise asynchronous system 192 (simulation 3). Table 1. Simulations used in this study. We use the mechanistic dengue model, with different temperature inputs and numbers of locations N , to address specific but related questions (column Aim). Each location was run as an independent simulation (no host movement between locations), using the same host birth and death rates and population sizes and starting conditions, such that the only difference across locations was temperature (see section "Further details on simulation studies" in SI). In all cases, simulations were run assuming mean cross-protections between dengue serotypes of six months, one year, and two years. Outputs had a monthly resolution.
Sim. N Temperature input Aim 1 72 Real temperature time series Reproduce empirical patterns using temperature as the input.   (Table S1 in SI). The phase angles between temperature and model output 207 across provinces were, as expected, consistent and statistically significant over most 208 time points (Fig 2e). Both the WMF and CWMF for the model output contain a 209 notable synchrony event in the early-1970s, reproducing that found in temperature; this 210 feature is missing from the WMF and CWMF estimated using the data (Fig 2b,d).

211
WMFs and CWMFs for model outputs using different assumptions on cross-protection 212 are shown in Fig S13 in SI, and show that with longer mean cross-protections (e.g., two 213 years), the degree of synchrony in the model output, and consistency in the phase angles 214 between modelled dengue and temperature, is lower than with shorter cross-protections. 215 A temperature gradient can generate asynchronous dengue dynamics 216 For the idealised experiments of simulation 2, we found that with at least a mean of one 217 year temporary cross-protection between serotypes, the resulting dengue dynamics 218 include multiannual cycles (Fig 4a-c) components in dengue go down (i.e., multiannual cycles are more frequent), and the 225 power at the multiannual timescales increases. At the highest mean temperature, with a 226 mean of two years of cross-protection, seasonal cycles all but disappear (Fig 4c) due to 227 two main trends. On the one hand, as cross-protection becomes longer, the power of the 228 multiannual cycles (and therefore their amplitude) increases (Fig 4a-c). On the other, 229 the seasonal variation in transmission is most limited at the highest mean temperatures 230 (Fig S37 in SI), thus producing more limited seasonal variation in dengue cases. Taken 231 together, these two trends mean that dynamics become increasingly dominated by the 232 multiannual cycles produced through cross-protection, and less so by seasonal variation 233 in temperature. While the periodicity of the multiannual cycles for each temperature 234 regime changes little between means of one and two years of cross-protection, their 235 prominence does: with longer cross-protection, these cycles clearly increase in mean 236 wavelet power. The range of dynamics produced across locations under all assumptions 237 of cross-protection are sufficiently diverse: no synchrony is detected (e.g., Fig 4d).  Because different temperature regimes produce distinct dengue dynamics (Fig 4a-c), 241 one simple explanation for synchrony could be that temperatures across the country 242 converge to a more similar temperature regime during those periods, resulting in more 243 synchronous dengue dynamics. However, this has not been the case in Thailand 244 (Fig S32 in SI). The next hypothesis is that a common multiannual fluctuation in 245 temperature can synchronise dengue dynamics (simulation 3). The result for simulation 246 3 was that when the amplitude of the multiannual fluctuation in temperature was at 247 least 20% of the seasonal cycle, the five locations were temporarily synchronised 248 (Fig 4e), under all assumptions on cross-protection. Synchrony in dengue occurred at 249 approximately the same timescale as that of the multiannual fluctuation in temperature, 250 so, for example, when the multiannual fluctuation in temperature was a single four-year 251 cycle, synchrony in dengue was detected at a four-year timescale. However, when the 252 multiannual fluctuation had an amplitude 10% that of the seasonal cycle and mean 253 cross-protection was two years, synchrony was not detectable (Fig S28 in SI).  Specifically, as the duration of cross-protection and/or mean temperatures increase (and 269 as seasonality in temperature is reduced), the power of the multiannual intrinsic 270 dynamics increases, such that the effect of the multiannual fluctuation in temperature is 271 overwhelmed. In our simulations, this is particularly the case for mean temperatures of 272 30 • C and a mean cross-protection of two years (Fig 5).

273
It is also worth noting that the multiannual fluctuation in temperature is not between simulations with a multiannual fluctuation in temperature, and those without 278 (Fig S29 in SI). These show both that the effect of the multiannual fluctuation in 279 temperature can be lasting, and also that the effect of the multiannual fluctuation is 280 reduced with longer cross-protections and higher mean temperatures.

281
This interaction between temperature and cross-protection explains why with longer 282 cross-protection, the degree of synchrony appears to go down: as cross-protection and by multiannual fluctuations in the environment using simulations. We compare the mean wavelet power of dengue at the timescale of the environmental multiannual fluctuation, with the mean wavelet power of dengue at the timescale at which wavelet power is maximised without a multiannual fluctuation in temperature (the multiannual timescale attributable to cross-protection alone). Positive values mean that the multiannual fluctuation in temperature plays a more important role driving dengue multiannual dynamics, and negative values mean that dengue dynamics are less sensitive to the multiannual fluctuation in temperature. For example, for a mean temperature of 30 • C, a multiannual fluctuation with a timescale of four years, and a mean cross-protection of two years, we compared the mean wavelet power at a timescale of four years (that of the multiannual fluctuation), with the mean wavelet power at a timescale of 2.7 years (the multiannual timescale at which wavelet power was maximised with no multiannual fluctuation in temperature, i.e., the peak on the top row of Fig 4c). In this example, the mean wavelet power was distinctly greater at the 2.7 year timescale than at the 4 year timescale, meaning dengue dynamics were insensitive to the multiannual fluctuation in temperature. Differences in wavelet power shown here have been normalised against the maximum value for each timescale of the multiannual fluctuation, to improve clarity. Overall, longer cross-protection (the right column of points) and mean temperatures (purple colors) are more likely to lead to a system dominated by intrinsic dynamics, and less sensitive to multiannual fluctuations in temperature. Mean wavelet power is estimated for the duration of the multiannual fluctuation only (so for a four-year multiannual fluctuation in temperature, mean wavelet power is estimated during those four years only).
Here, the multiannual fluctuation has an amplitude 0.2 times that of the seasonal cycle, but the patterns are qualitatively similar for the other amplitudes of the multiannual fluctuations.

287
Using multiple approaches, we have described an interesting spatial phenomenon in 288 dengue: the country has repeatedly and periodically moved in and out of synchrony in 289 dengue cases. In general, periods of greater synchrony coincided with larger outbreaks. 290 When analysing temperature time series across Thailand, we found a similar pattern in 291 spatial synchrony. Although the correlation between the overall patterns in synchrony 292 in temperature and dengue cases was not statistically significant, the relationship 293 between temperature and dengue cases was particularly consistent across provinces 294 during times of greater synchrony. Accounting for the dynamics of immunity is essential 295 for understanding dengue dynamics; any effects of temperature on dynamics are 296 modulated by immunity. To gain a more mechanistic understanding of how temperature 297 interacts with the dynamics of immunity to produce the observed patterns in synchrony, 298 we adapted a mechanistic dengue model with temporary cross-protection. When 299 running the model using the observed temperatures across Thai provinces, the patterns 300 in synchrony in the resultant model output dengue time series correlated positively and 301 significantly with observed patterns in dengue cases from the passive surveillance data. 302 The fact that this correlation was statistically significant, but the one between the 303 patterns in synchrony in dengue cases and temperature was not, further highlights the 304 importance of nonlinearities in the way temperature affects dengue and dynamics of 305 immunity. While (a sufficiently long) cross-protection alone can produce multiannual 306 cycles, we found that their periodicity was modulated by temperature; different 307 temperature regimes produced characteristically different dengue dynamics. However, 308 although the timescales of the multiannual periodicities of dengue varied little, their 309 power became more prominent with longer cross-protections. This result sheds some 310 light on the roles played by both temperature and dynamics of immunity, and suggest 311 that both are essential to understand patterns in dengue dynamics. We also found that 312 a range of temperature regimes produced asynchronous dengue dynamics, but the

327
Most locations' temperature appear to be increasing in mean and decreasing in 328 seasonal variation (Fig S33 and S36 in SI). If these trends were to be maintained, our 329 results suggest that dengue dynamics could slowly shift too, and approach a regime 330 characterised by higher average incidence (Figs S15-S17 in SI), and dynamics 331 increasingly dominated by multiannual, rather than seasonal, cycles. Multiannual cycles 332 may also be, to greater extents, dictated by intrinsic factors over the multiannual 333 timescales present in temperature. These trends may also suggest that in the future, 334 synchrony may be driven by increasing similarity in temperatures across Thailand, and 335 less so by synchronous multiannual fluctuations in temperature. However, these 336 hypotheses are contingent on a more precise understanding of the interactions between 337 serotypes and contributions from other intrinsic factors (such as changes in demography 338 and movement of hosts between locations).

339
There are several caveats to our analyses. We chose to focus on temperature and its 340 interaction with the dynamics of immunity. However, we could not ascribe every 341 pattern in synchrony in dengue to synchrony in temperature. For example, Thailand 342 experienced a strongly synchronous two-year fluctuation in temperature around 1972, 343 and while this fluctuation was reproduced in our model output, it was conspicuously 344 absent from the observed dengue dynamics. There are multiple reasons that could 345 explain this discrepancy. During the earlier part of the records, dengue cases might have 346 been detected less effectively, thus potentially obscuring multiannual patterns. On the 347 other hand, other environmental variables (notably precipitation) and their interaction 348 with temperature are likely important for the dynamics of dengue. For instance, the 349 extent to which a synchronous event in temperature can synchronise multiannual cycles 350 in dengue might depend on the amount of precipitation or timing of the rainy seasons 351 during those years. Additionally, other non-environmental mechanisms might also 352 influence synchrony in dengue. Complex multiannual cycles, synchrony, and variations 353 in the degree of synchrony over time can in theory be produced via, for example, the 354 movement of hosts between locations (weakly-coupled oscillators; [15]). There are also a 355 range of concurrently ongoing changes in Thailand that may also influence multiannual 356 dynamics and synchrony. For instance, the birth rates across the country have been 357 declining [48], and the rates of decline may be spatially heterogeneous. Similarly,   Wavelet transforms 391 We applied continuous wavelet transforms (WTs) to explore how the oscillatory 392 behaviour of dengue cases has changed over time. WTs have now been applied 393 extensively in the study of infectious disease dynamics [31,34,[50][51][52] and more broadly 394 in ecology [45,46]. WTs decompose time series into frequency components, but, as 395 opposed to Fourier transforms, WTs are also localised in time. Thus, WTs can be used 396 to characterise non-stationary time series and how the relative importance of different 397 frequencies change over time. The basis for the WT is a "mother" wavelet, a wave 398 localised in both frequency and time (i.e., of limited duration), which is then scaled or 399 stretched (for the frequency component) and shifted along the time axis (for the 400 temporal component) to derive a set of "daughter" wavelets. The WT can then be 401 understood as a correlation between (or more specifically the convolution of) the time 402 series and the set of daughter wavelets [53]. 403 We use the WT as implemented in R package 'WaveletComp' v1.1 -details are 404 provided in the package documentation [54] and code, and a summary is given here.

446
To test whether synchrony is statistically significant in the multiannual range 447 against a null hypothesis of no synchrony, we focus on the consistency of phases across 448 locations. We do so by estimating the wavelet phasor mean field (WPMF), defined as which retains information on the complex phases of the transforms [46]. When different 450 locations have a similar phase, the value WPMF will be large, while when the phase in 451 each location is independent from all other locations, the value of the WPMF will be obtained, as above, by preserving the same Fourier spectrum as the original time series 479 of dengue cases, but assuming asynchrony. The Pearson and Spearman correlations 480 between each of these surrogate WMFs and the temperature WMF or the model output 481 WMF produce the null distributions against which we compare the observed 482 correlations.

483
To determine in greater detail whether the patterns of synchrony in temperature 484 were consistent with those observed in the dengue cases, we extended the method above 485 for WMFs [46] to calculate cross-wavelet mean fields (CWMFs). Cross-wavelets 486 quantify the similarity in two time series' (temperature and dengue cases) wavelet power 487 at each scale and time point, so the CWMFs quantify the consistency of this similarity 488 across locations. CWMFs are an average of the power-normalised cross-wavelets at each 489 location between temperature and dengue cases. A cross-wavelet X between two time 490 series is defined as where the cross-wavelet is corrected following Veleda et al. [57]. The argument of X using temperatures for Thailand we were to produce a CWMF using temperatures from 500 a different location, with similar temporal and spatial autocorrelation, we could 501 potentially observe similar patterns in synchrony, despite the two variables being clearly 502 unrelated. For this reason, we test significance in the multiannual timescales using the 503 same three-step approach described above, but here, surrogate dengue datasets are 504 produced that not only preserve the autocorrelation structure of the time series, but 505 also the cross-correlation structure across locations (i.e., synchrony-preserving 506 surrogates). We do this by generating time series with the same Fourier spectrum as the 507 original time series (as above), but then adding the same random uniformly distributed 508 phase at each frequency across all time series [46]. As above, we test significance by 509 focussing on the consistency of phase angles between dengue and temperature. Using 510 these surrogate dengue datasets, cross-wavelets are estimated relative to the real 511 temperature dataset, and these in turn are used to produce the surrogate cross-wavelet 512 phasor mean fields, which we use for comparison. Material has been reviewed by the Walter Reed Army Institute of Research. There is no 521 objection to its presentation and/or publication. The opinions or assertions contained 522 herein are the private views of the author, and are not to be construed as official, or as 523 reflecting true views of the Department of the Army or the Department of Defense.