Climatological, virological and sociological drivers of current and projected dengue fever outbreak dynamics in Sri Lanka

The largest ever Sri Lankan dengue outbreak of 2017 provides an opportunity for investigating the relative contributions of climatological, epidemiological and sociological drivers on the epidemic patterns of this clinically important vector-borne disease. To do so, we develop a climatologically driven disease transmission framework for dengue virus using spatially resolved temperature and precipitation data as well as the time-series susceptible-infected-recovered (SIR) model. From this framework, we first demonstrate that the distinct climatological patterns encountered across the island play an important role in establishing the typical yearly temporal dynamics of dengue, but alone are unable to account for the epidemic case numbers observed in Sri Lanka during 2017. Using a simplified two-strain SIR model, we demonstrate that the re-introduction of a dengue virus serotype that had been largely absent from the island in previous years may have played an important role in driving the epidemic, and provide a discussion of the possible roles for extreme weather events and human mobility patterns on the outbreak dynamics. Lastly, we provide estimates for the future burden of dengue across Sri Lanka using the Coupled Model Intercomparison Phase 5 climate projections. Critically, we demonstrate that climatological and serological factors can act synergistically to yield greater projected case numbers than would be expected from the presence of a single driver alone. Altogether, this work provides a holistic framework for teasing apart and analysing the various complex drivers of vector-borne disease outbreak dynamics.

susceptible mosquito becomes infectious and able to transmit DENV to future hosts following a delay time known as the extrinsic incubation period (EIP) [5]. The EIP is dependent on the environmental temperature, along with the life cycle duration and survival of the mosquito [5,6]. In addition, precipitation and the local hydrological conditions can affect the quantity of pooled surface water, which serves as the breeding sites in which mosquitoes deposit their eggs. In particular, for dengue, Aedes mosquitoes are primarily container breeders (natural and artificial) and thrive in both clean and organically rich water [7]. Consequently, climate variables interact in complex ways with the biological parameters of the vector to determine not only the viable geographical range for the transmission of the disease [1,6,8] but also the timing and seasonality of infection in regions where the disease is endemic.
In 2017, Sri Lanka experienced its largest recorded dengue outbreak, with 186 101 patient hospitalizations reported [9]. DENV has been endemic in Sri Lanka since the mid-1960s, and the four primary virus serotypes have been co-circulating there for more than three decades [7]. Furthermore, topographical influences as well as the seasonality of the monsoon cycles result in highly variable rainfall characteristics across the compact island nation [10], providing a unique opportunity to analyse localized interactions between climate and disease transmission. Indeed, several studies have determined a role for climatological drivers in the temporal dynamics of dengue incidence in Sri Lanka [11,12] (and more broadly [1,[13][14][15]). In Sri Lanka, a variety of other factors have also been implicated in the epidemic and transmission dynamics of dengue, including changes in viral serotype [9], high population densities in affected areas [16], and poor hygiene conditions in urban environments [17]. This complex dynamical network of potential drivers is depicted schematically in figure 1. Recent work related to a dengue epidemic in Guangzhou, China has demonstrated that dengue epidemic periods in geographical regions with inter-annual transmission variability can serve as important case studies for assessing the relative importance of different drivers on DENV dynamics [19].
In this work, we systematically assess potential roles for climatological drivers including seasonal patterns in average temperature and precipitation as well as extreme weather events, epidemiological drivers including serotype change, and sociological drivers including human movement patterns in response to a severe drought, on the outbreak dynamics of dengue in Sri Lanka. We consider four distinct time periods in this analysis: the years 2010-2015 when endemic transmission patterns persisted, the transition year of 2016, the epidemic year of 2017, and finally future transmission patterns. We find that although seasonal climatological conditions play an important role in establishing monthly dengue transmission patterns in Sri Lanka, they alone are unable to account for the severity of the 2017 epidemic. Using a simplified two-strain susceptible-infected-recovered (SIR) model, we show that the re-introduction of a DENV serotype may have played an important role in driving the epidemic, and that in combination, climatological and serological drivers can act synergistically to yield greater projected case numbers than the presence of a single driver alone.

Results and Discussion
In figure 2, an overview of the dengue incidence time series as well as relevant demographic and climatological information for Sri Lanka is provided. Dengue case data at a national level between the years 2010 and 2018 is plotted in figure 2a, showing a striking peak in case numbers in 2017. Population levels by district, averaged over the years 2010-2018, are presented in figure 2b.  Figure 1. Schematic depicting the various potential drivers of the transmission dynamics of dengue as well as the relationships between these drivers. The depicted phylogeny of DENV is from [18].
In addition to substantial demographic variation by district, the tropical monsoon climate of Sri Lanka is characterized by elevated minimally varying average temperatures as well as highly variable rainfall characteristics [10]. Two precipitationbased climate zones can be distinguished within the island: the wet zone lying to the west of the central highlands, and the relatively drier zone lying to the east and north [10]. In general, the wet zone is exposed to the southwest monsoon from the months of May through to September as well as the first intermonsoon (lasting from March through to April), while the dry zone experiences substantial rainfall during the second intermonsoon (October through to November) and the northeast monsoon of December through to February [10].
In figure 2c,d, these climatic features for the pre-epidemic years of 2010-2016 are summarized and contrasted with the seasonality of dengue cases. In figure 2c, the mean monthly number of dengue cases, precipitation, and average temperature between the years 2010 and 2016 averaged across all districts are plotted as the solid lines, with the shaded areas indicating the standard error across districts. One immediately apparent feature is the relative uniformity of temperature both temporally over the year and across districts, as opposed to the clear seasonality in both precipitation and dengue cases, with larger standard errors indicative of greater variation between districts. Furthermore, two peaks in precipitation and dengue can be distinguished at this macroscopic level, with the former leading the latter by a period of two months.
In order to explore regional variation in seasonality of precipitation and dengue, a monthly relative entropy metric S [20,21] is calculated by district for these same data spanning the years 2010-2016, where In the centre panel of figure 2d, the entropy metric S corresponding to the monthly mean precipitation between the years 2010 and 2016 is presented for each district. A notable geographical trend is observed, with low levels of precipitation entropy, or minimal seasonality, in the wet zone or southwest region of the country, and increasing precipitation entropy, or increasing seasonality, towards the northeast and the climatic drier zone. To illustrate the precipitation patterns corresponding to these values of S, time series of mean monthly precipitation from 2010 to 2016 with standard deviations   Dashed grey lines in all charts are aligned with local maxima in precipitation, while the solid grey lines are aligned with the peaks in dengue cases. Together, these illustrate a two month lag between peak rainfall and peak dengue incidence.
(d ) Monthly entropy metric S by district for precipitation (centre panel) and dengue cases (right panel) from 2010 to 2016. Darker colours indicate districts with greater seasonal variability. Time series of monthly averaged precipitation (solid line) and standard deviation across years (shaded regions) for this same time period corresponding to the districts with the minimum (Nuwara Eliya, S ¼ 0:05), and maximum (Jaffna, S ¼ 0:67) precipitation entropy are also shown in the left panel. The approximate location of the national capital city of Colombo is indicated by a star in (b) and in the centre and right panels of (d).
across years indicated corresponding to the districts with the minimum (Nuwara Eliya, S ¼ 0:05), and maximum (Jaffna, S ¼ 0:67) precipitation entropy are shown in the left-hand panel. In the right-hand panel of figure 2d, this entropy calculation is repeated for mean monthly dengue case incidence between the years 2010 and 2016, and plotted for each district. Although less clear than that for precipitation, a similar trend of increasing seasonality in dengue case numbers is observed towards the northeast drier zone of Sri Lanka. Motivated by the apparent correlation between precipitation and dengue cases observed in figure 2d as well as previous evidence of a role for climatological and other factors in establishing the temporal dynamics of dengue in the region [11,12], we analyse these epidemiological data within the framework of various disease drivers during four distinct periods: the years 2010-2015 when endemic transmission patterns persisted, the transition year of 2016, the epidemic year of 2017, as well as future years extending until 2099. We develop and use a climatologically informed time-series susceptible-infected-recovered (TSIR) model [22], in which temperature and precipitation levels or approximated serotype susceptibility play a role in determining either the dynamical transition of individuals within a given population from the susceptible to the infected category or the relative sizes of these subpopulations. The details of the development of this model are provided in the Material and methods section.

2010-2015: endemic dynamics
We begin by studying the ability of the TSIR model [22] to capture the more regular endemic dynamics of dengue incidence In order to evaluate the predictive ability of the TSIR model during these endemic years, we train the model on district-level case data from 2010 to 2014, and assess its ability to predict the 2015 case numbers. In figure 3b, the number of cases from 2010 to 2015 adjusted by the reporting rate ρ obtained from the TSIR model fitted to the case data from 2010 to 2014, the TSIR fitted from 2010 to 2014, and the predicted case numbers for 2015 are plotted as the solid blue, solid black and dashed red curves, respectively, for the districts in which the ratio C 2015 /C 2010-2014 deviates the most (Ampara, C 2015 /C 2010-2014 = 0.36), and least (Vavuniya, C 2015 /C 2010-2014 = 0.99) from unity. Here, we consider only districts included in the climatological regression model, thus excluding (for example) Mullaitivu, which has a ratio of C 2015 /C 2010-2014 = 2.06 (as described in the Material and methods section, Mullaitivu is excluded from this analysis as an outlier, consistent with its exclusion from the climatological regression model). As can be seen, the TSIR model produces a satisfactory prediction of the 2015 case numbers in both districts, with cyclic seasonal trends well-captured despite disagreement in quantitative values. This is also summarized in the scatter plot of the number of dengue cases in 2015 adjusted by the reporting rate ρ obtained from the TSIR model fitted to the data from the years 2010 to 2014 versus the TSIR predictions for the number of infected individuals in 2015 for the districts of Vavuniya ( purple colour scheme, circles) and Ampara (yellow colour scheme, triangles) presented in figure 3c. These data are further organized by season, with the lightest hues corresponding to winter months (December, January, February (DJF)), then spring months (March, April, May (MAM)), summer months (June, July August (JJA)) and finally the darkest hues correspond to autumn months (September, October, November (SON)). In Vavuniya, the observed and predicted case numbers agree well (lie near the dotted black line indicating parity) apart from in the winter months when the observed case numbers increase sharply. By contrast, the case numbers are over-predicted in Ampara for all seasons.

2016: transitional dynamics
In certain districts, the onset of the epidemic period occurred as early as partway through 2016, thus warranting a separate consideration of this transitional year. To illustrate this, in figure 3d, the ratio C 2016 /C 2010-2015 of dengue cases by district in 2016 to the average of those reported in the years 2010-2015 is shown. As can be seen, in contrast with figure 3a, the value of this ratio is now greater than unity in nearly all districts. As in §2.1., we next assess the ability of the TSIR model to predict the 2016 case numbers in the district for which C 2016 /C 2010-2015 deviates the most (Galle, C 2016 / C 2010-2015 = 2.76), and least (Puttalam, C 2016 /C 2010-2015 = 1.04) from unity. Here, we consider only districts included in the climatological regression model, thus excluding (for example) Mannar, which has a ratio of C 2015 /C 2010-2014 = 1.00 (as described in the Material and methods section, Mannar is excluded from this analysis as an outlier, consistent with its exclusion from the climatological regression model).
In figure 3e, the number of cases from 2010 to 2016 adjusted by the reporting rate ρ obtained from the TSIR model fitted to the case data from 2010 to 2015, the TSIR fitted from 2010 to 2015, and the predicted case numbers for 2016 are plotted as the solid blue, solid black, and dashed red curves, respectively, for the districts of Puttalam and Galle. As can be seen, for the district of Puttalam where the epidemic had not begun by the end of 2016, the TSIR model produces a satisfactory prediction of the 2016 case numbers, with the climatological seasonal cycle well-captured despite disagreement in quantitative values. Conversely, in Galle, where the epidemic begins midway through the year 2016, there is substantial disagreement between the prediction of the TSIR model and the reported case numbers.

2017: epidemic dynamics
Next, we consider the dynamics of the epidemic year of 2017. As can be seen in figure 4a, dengue incidence in 2017 was higher in all districts compared to the previous 7 years, ranging from a 10.7-fold increase in the district of Trincomalee to a 2.3-fold increase in the district of Mannar. In this section, we seek to explore the possible drivers that may have contributed to these epidemic dynamics.  for initial number of susceptibles (10, 11) (9, 10) (8, 9) prediction for 2016 was satisfactory) as well as the district of Colombo which contains the national capital city of the same name. These results are shown in figure 4b. We additionally perform the model prediction for 2017 using values for the transmission rate β derived from the climatological regression model (see the Material and methods section) that explicitly accounts for the climate of 2017, and plot these results as the solid red curve in figure 4b. As can be seen, although the timing of dengue case seasonality is relatively well captured by the TSIR predictions in both districts, even when the climate from 2017 is explicitly accounted for, the reported case numbers are substantially under-predicted. Consequently, it would appear that although climate plays an important role in establishing the timing of annual dengue transmission patterns in Sri Lanka, the severity of the 2017 epidemic cannot be explained by an anomalous climate that year. Figure 3b is reproduced for all districts in the electronic supplementary material, figure S8.

Virological drivers
As illustrated in figure 4c, according to hospital surveillance reports, DENV1 was the predominant serotype in circulation in Sri Lanka from 2009 until mid-2016, during which time DENV2 and DENV3 were not detected [9]. DENV2 and DENV3 were the predominant serotypes in circulation prior to 2009, and indeed all four virus serotypes have been co-circulating in Sri Lanka for more than three decades [7,9]. In mid-2016, the re-emergence of DENV2 coincided with the onset of the epidemic [9]. Sequencing studies of the 2016 Sri Lankan DENV2 strain obtained from patients who were treated for acute dengue infection at the National Institute of Infectious Disease confirmed it to be from the same clade as those circulating in 2015-2016 in Singapore, and indicated that it shared a common ancestor with the 2014 Malaysian strains [9]. In fact, these studies found that the 2016 Sri Lankan DENV2 strain was genetically more distant from the DENV2 strains that circulated in Sri Lanka from 1981 to 2004 than these Singaporean and Malaysian strains [9], suggesting case importation. The phylogeny of DENV from [18] is included in the electronic supplementary material, figure S1a. In order to estimate the possible role for serotype reemergence on the epidemic dynamics, we model the predicted effect of an effective increase in population-level susceptibility at the start of the epidemic. This increase can be interpreted as representing individuals in the population who had previously been infected with DENV1 only (and were hence still susceptible to DENV2 infection), or an increase in the severity of secondary DENV infections relative to primary ones [2] following replacement of the dominant strain. Indeed, over short time periods, the TSIR yields similar predicted increases in case numbers when either the initial number of susceptible individuals or the transmission rate β is increased. It is worth noting, however, that in one clinical report [9], all patients with acute dengue infection recruited after January 2017 by the National Institute of Infectious Disease of Sri Lanka were infected with DENV2, and that high levels of patients presenting with both DENV1 (pre-epidemic) and DENV2 (post-epidemic) were experiencing secondary infections. Whether this extends to the district-level data used in this paper is unknown. We emphasize then that with these clinical observations in mind, the analysis in this section serves as a reasonable starting point for considering the role of virological drivers given the absence of serotype information. Further caveats are discussed in §2.5.
In figure 4d, the multiplicative factor on the initial value of susceptible individuals at the start of 2017 that minimizes the difference in the peak number of reported and predicted cases in 2017 is plotted for each district, and ranges between a value of 1.11 and 2.61. In figure 4e, the number of cases from 2010 to 2017 adjusted by the reporting rate ρ obtained from the TSIR model fitted to the case data from 2010 to 2016, the TSIR fitted from 2010 to 2016, and the predicted case numbers for 2017 are plotted as the solid blue, solid black and dashed red curves, respectively, for the district of Colombo. Additionally, we plot the TSIR prediction for 2017 following re-initialization of the initial number of susceptible individuals by the multiplicative factor obtained in figure 4d as the solid red curve.
In order to contextualize these determined increases in the number of susceptible individuals, we consider the simplified scenario of a two-strain system, in which we assume that a single DENV serotype has achieved steady-state transmission dynamics prior to the arrival by exogenous importation of a new viral strain. As described above, this approximation is consistent with sequencing results of the genetic origins of the 2016 Sri Lankan DENV2 strain. However, this scenario ignores the complexities of crossimmunity associated with the long-term co-circulation of all four DENV strains in Sri Lanka. Further discussion of the caveats of this two-strain system are discussed in §2.5. Nevertheless, in this simplified scenario, the number of individuals susceptible to the newly imported strain and/or the previously circulating one would exceed the steady-state number of susceptible individuals to the previously established strain only. A schematic of this simplified two-strain system based on that developed in [23] is shown in the electronic supplementary material, figure S1b and the mathematical details are described in the Material and methods section.
At the time of importation of serotype 2 infections, we approximate the actual number of susceptible individuals S 0 as , where S is the number of individuals susceptible to infection by both strains, S* = N/β(γ + μ) is the steady-state value of this quantity when a single strain is in circulation (see the Material and methods section), I 1 and R 1 are, respectively, the number of individuals presently infected with or recovered from serotype 1 infections only, and N is the total number of individuals. Furthermore, μ, γ, β and ρ 2 denote the rates of birth, death, transmission and importation of exogenous infections of serotype 2, respectively. The analytic expression for the factorial increase in the susceptible population in equation (  The 2017 drought had a severe impact on the agricultural sector, causing a 40% decrease in the cultivated area and production of rice as depicted in figure 5c. Economic theory and empirical research in a number of developing economies including Senegal have shown that increases in the ratio of urban-to-rural wages, possibly driven by decreases in agricultural output, provide the foundation of rural-urban migration [26,27], and news reports suggest that this phenomenon occurred in Sri Lanka during this period [28]. Changes in human social ecology, including population size and density and urbanization, have been implicated in altering the probabilities of infectious disease emergence and transmission by creating highways for 'microbial traffic' and providing opportunities for rapid dissemination of emerging infections [29]. These various factors suggest that indepth studies into the effect of human mobility patterns on the dengue epidemic of 2017 may be warranted, such as those conducted using mobile phone data in Pakistan [30] and for rubella in Kenya [31]. It is nevertheless interesting to consider the magnitude of the population that would be required to relocate in order to achieve the fractional increases in the susceptible population calculated in §2.3.2. In the district of Colombo, for instance, a 27% increase in the susceptible population (i.e. a fractional increase of S 0 = 1.27S*) corresponds to an influx of 256 000 susceptible individuals; approximately 11% of its population in 2017.

Other drivers
Lastly, we note that the length of the approximately 7.5 year period from 2009 until mid-2016 during which DENV1 was the primary serotype in circulation prior to the re-emergence of DENV2 is broadly consistent with the periodicity of serotype prevalence of 8-9 years determined both theoretically and from epidemiological data in [32], although shorter cycles of 3-5 years have also been observed in both data and simulations [2]. Depending on the modelling framework, a number of factors have been found to contribute to the maintenance of these multi-year cycles, including virological cross-protection [2] as well as the segregation of hosts via spatially arranged communities with the possibility of importation of infections via temporal human movements [32]. Therefore, although these complex multi-annual patterns are undoubtedly driven by a number of factors, it is possible that these periods of serotype prevalence define an upper bound on the allowed build-up of the susceptible cohort, beyond which the re-introduction of a non-predominant viral strain may result in the onset of an epidemic. royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 17: 20200075 the panel regression model with all districts included. The coefficients on the climate variables used for the transmission rate projections for these three districts, however, were set to those obtained by the regression model when these three districts were excluded. Furthermore, in obtaining the future values, it is important to note that only climatological variables were modified in the panel regression model and all location-and time-specific dummy variables were held constant. In order to better interpret these predictions for the empirical transmission rate, we re-run the TSIR model for 2017, as was done in §2.3.1., using the projected values for the transmission rate obtained from the panel regression model with climate data for 2080-2099. Furthermore, because serotype invasion was found in §2.3.2. to be a potentially important driver of outbreak dynamics, we also consider the compound scenario in which the same climatological factors are considered and in addition the baseline number of susceptible individuals is increased to S 0 (see equation (2.2)). Consequently, these simulations are relevant for considering the effect of altered climate variables and/or the introduction of a previously absent serotype when all other demographic variables including birth rates and population levels are held constant.

Projections for future transmission dynamics
In figure 6b, the total number of cases for the year 2017 under these various scenarios are reported for each district. The blue circles and solid blue line correspond to the actual case numbers adjusted for the reporting rate ρ obtained from the TSIR model fitted to the case data from the years 2010 to 2016; the yellow circles and dashed-dotted yellow line correspond to the predicted case numbers when the transmission rate obtained from the panel regression model with 2080-2099 climate data is used; the red circles and dashed red lines correspond to the scenario of simulated serotype invasion shown in figure 4e, in which only the initial number of susceptible individuals is adjusted by the district-specific multiplicative factor reported in figure 4d; and the green circles and solid green line correspond to the scenario of both 2080-2099 climate and serotype invasion. Once again, substantial variation is observed across districts, but critical increases in case numbers under the compound climatological and serological driver scenario are observed in a number of districts, including the most populous district of Colombo in which the national capital city of the same name is located. In quantitative terms, in Colombo the total predicted case numbers under the compound driver scenario are 2.84-fold and 2.44-fold higher than those predicted by the scenarios of 2080-2099 climate and serotype invasion alone, respectively. Dengue incidence across the year is shown for Colombo in figure 6c under these same conditions, as well as the predictions of the TSIR model for 2017 following training on the 2010-2016 data only (the dashed black curve). Figure 6c is reproduced for all districts in the electronic supplementary material, figure S9.

Caveats
This work provides an interdisciplinary analysis of the potential drivers of dengue dynamics in Sri Lanka, and consequently it is important to explicitly mention a number of caveats that should be considered when interpreting the findings presented herein. Firstly, as described in the Material and methods section, the panel regression model assumes a Brière functional form for the average temperature in determining the empirical transmission rate. This was done to capture the non-monotonic thermal responses generally observed for mosquito and viral traits [33]. One immediate issue with this approach, however, is the shortage of data points on the decreasing portion of the Brière curve (see the electronic supplementary material, figure S5c). Although we lack observational data to populate this curve owing to the limited range of average temperatures contained within the regression panel, we hypothesize that above a critical value, increasing values of average temperature would correlate negatively with log-transmission rate. Indeed, several of the critical temperature values estimated in [33] lie above the range of average temperature values measured in our climatological dataset. Somewhat reassuringly, however, because the predicted future values of average temperature (see the electronic supplementary material, figure S2) generally lie below the critical thermal response temperature obtained for the panel regression model (see the Material and methods section), the assumption of either a Brière or linear functional form for average temperature yielded similar predictions for log-transmission rates. One interesting exception is the district of Nuwara Eliya, whose yearly average temperature values in the present day as well as projected future values generally lie below the value of T min determined for the Brière function (see the electronic supplementary material, figure S2). The assumption of a Brière functional form for average temperature in the regression model consequently results in minimal projected changes in log transmission (figure 6a), while a much greater increase is projected when a linear functional form is used (electronic supplementary material, figure S6). Significant variation across districts was obtained for the lag times between dengue case numbers and both average temperature and precipitation (see the electronic supplementary material, figure S4), and these values were also sensitive to the number of years of data considered by the panel regression model. In order to capture relevant long-term trends, we chose to develop the panel regression model using the maximum number of years of epidemiological data available prior to the epidemic year (2010-2016), and assumed a single lag value between the case data and climate variables for the entire country by selecting that which occurred most frequently across districts. Nevertheless, if more accurate dengue case predictions based on future climate projections for a specific district are desired, it may be more appropriate to use district-specific as opposed to national lag values for the climate variables, and to adjust the years considered in the panel regression model accordingly.
Furthermore, given the absence of serotype data, it is important to carefully consider the limitations of the serologically-motivated analysis in § §2.3.2. and 2.4. As previously mentioned, modelling an effective increase in populationlevel susceptibility or an increase in the transmission rate β yield similar increases in predicted case numbers using the TSIR model over short-time periods. These scenarios may represent additional cases owing to secondary infections and potential increases in the severity of secondary DENV infections relative to primary ones [2]. While both of these scenarios are relevant to multi-strain pathogenic systems, they highly simplify the reality of DENV infections in Sri Lanka, which include endemic co-circulation of four viral strains and temporary periods of cross-immunity between heterologous strains [2]. There is evidence, however, that for multi-year periods, one DENV strain dominates and is eventually replaced by another serotype (see case data from Thailand [2] and Puerto Rico [32], for instance). In that sense, at specific periods in time corresponding to strain replacement, it may be justified to think about DENV infection as a two-strain system, although in order to model long-term serotype interaction dynamics this approach is clearly not justified. Therefore, given the absence of serotype-specific case data but clinical evidence of serotype invasion, our approach provides a starting point for analysis despite its numerous limitations. Furthermore, in our analytical consideration of the potential magnitude of case importation, it was assumed that R 0 was identical for primary and secondary dengue infections. Although this may be correct to first order, it is certainly possible that transmission rates may vary between serotypes. Additionally, in the derivation of the multiplicative increase factor for the number of susceptible individuals (see the Material and methods section and equation (2.2)), it was assumed that the number of imported cases was much smaller than the number of autochthonous ones. Although this assumption may hold true in regions where dengue is endemic, it may not be valid in non-endemic areas.
Finally, as described in §2.4., projections of future transmission rates and dengue case numbers were obtained by considering modifications to the climatological and serological drivers of the TSIR model only. In other words, inherent in these projections are assumptions of unchanging location-and time-specific dummy variables in the panel regression model, as well as unchanging demographic patterns including birth rates and population levels in the TSIR model. Consequently, our simulations are designed to examine the effect of a changing climate and serotype invasion on dengue transmission all other things being equal, and care should be taken when basing claims about future dynamical trends on these results.

Conclusion
Dengue is a geographically widespread disease of major clinical significance, whose transmission dynamics are well known to be modulated by a number of drivers including the local climate, human ecology and behaviour, and biological interactions between the four distinct primary serotypes. Particularly in light of the absence of vaccines or antivirals for this disease, there remain important opportunities to analyse the transmission dynamics of dengue in a holistic fashion by integrating these distinct and complex drivers. In this work, we do so in the context of the largest ever dengue epidemic in Sri Lanka during 2017, specifically considering climate, serotype invasion, and human mobility as epidemic drivers.
Sri Lanka provides a particularly powerful location for studying the effect of climate on dengue transmission as a result of highly variable rainfall patterns within the island. By developing a climatologically informed TSIR model, we found, in agreement with previous work, that seasonal climatological conditions play an important role in establishing yearly dengue transmission patterns in Sri Lanka, however they were unable to independently account for the severity of the 2017 epidemic. Motivated by hospital reports of a shift in predominant serotype from DENV1 to DENV2 midway through 2016, we modified a simplistic two-strain model for dengue transmission and demonstrated that increasing the number of susceptible individuals reconstructed by the TSIR model at the start of the epidemic year in order to account for secondary infections, was sufficient to capture the epidemic case numbers of 2017. Additionally, evidence of significantly reduced agricultural output during 2017 as a result of a series of severe climatic events suggests that a careful analysis of internal migration during this year as a driver of increased dengue transmission may be warranted. Critically, our simulations also suggest that the combination of climatological changes and serotype invasion can create compound effects on dengue transmission, resulting in even greater predicted case numbers than when only a single driver is assumed to be present. This possibility of synergistic interactions between disease drivers further reinforces the need to consider the epidemiology of dengue holistically.
To conclude, a number of important studies have demonstrated a role for the various drivers presented in figure 1 in the outbreak dynamics of dengue [1,2,[13][14][15]23,[30][31][32]. Here, we establish a holistic framework for teasing apart and analysing these various drivers, thus bringing together much of this prior work in the context of the 2017 Sri Lankan dengue epidemic. Such a multifaceted analysis may prove to be increasingly important in the face of changes in human ecology driven by urbanization and development, as well as climate change.

Dengue case data
Dengue case data were obtained from the publicly available disease surveillance dataset provided online by the Epidemiology Unit of the Sri Lankan Ministry of Health [34]. Cases are defined based on the presentation of a required number of symptoms. Cases are reported based on three types of surveillance: passive surveillance (submission of a notification card to a Medical Officer of Health), sentinel site surveillance (cases entered into webbased system by hospital officer), and special surveillance (used to identify the dynamics of dengue fever and dengue haemorrhagic fever). Monthly case data ranging from January 2010 to December 2018 were recorded for each of the 25 districts.

Demographic data
District-level demographic data for Sri Lanka were obtained from the CEIC database [35] as well as the Department of Census and Statistics of Sri Lanka [36].

Climatological data
Daily surface air temperature data were obtained from the ERA-5 global atmospheric reanalysis dataset [37]. The native 0.25 degree gridded dataset was interpolated bilinearly to a 0.05 degree grid. Temperature data for a given district were obtained by averaging the values at all grid points located within its geographical limits.
Daily precipitation data were obtained from the Climate Hazards Group InfraRed Precipitation with Station (CHIRPS) 0.05 degree gridded dataset [38,39]. As with temperature, precipitation data for a given district was obtained by averaging the values at all grid points located within its geographical limits.
Monthly SPEI data with a three-month time scale were obtained from the 1 degree gridded SPEI Global Drought Monitor dataset [40]. The time series of SPEI was computed by averaging SPEI values of all pixels within the island boundaries.
Climate projection data were obtained from the CMIP5 multi-model means dataset MMM-v2 which averages over 38 models under Representative Concentration Pathway 8.5 [41]. The native 2 degree gridded dataset was interpolated bilinearly to a 0.05 degree grid. As for the reanalysis data, the value for a given district was obtained by averaging the values at all grid points located within its geographical limits. Projections for temperature T fut (i, m) during month m in district i in the future period of interest were obtained using the formula

Agricultural data
Data for annual district level cultivated area and rice production were obtained from the Department of Census and Statistics of Sri Lanka [36]. The data included two cultivation seasons: Maha from September to March, and Yala from May to August. Annual values were obtained by summing the cultivated area and production data for the two cultivation seasons every year.

Development of the climatologically informed transmission model
SIR models form a family of simple mass-action epidemic models in which individuals within an assumed well-mixed population transition between the categories of susceptible (S), infected (I) or recovered (R) depending on their disease status to a given pathogen [42]. In this framework, transitions between categories are governed by the demographic processes of birth (which feeds the population of susceptible individuals) and death, as well as the disease transmission rate β driving susceptible individuals into the infected category and the recovery rate γ at which individuals move from the infected to recovered class. For both directly-and vector-transmitted diseases, the transmission rate β is frequently seasonal, reflecting yearly variation in disease drivers such as climatic conditions or demographic considerations such as the timing of school terms [22]. The stochastic, discrete-time TSIR model [22] provides a computationally inexpensive, highly tractable method for calibrating these seasonally varying, continuous time SIR models against time-series case data [42,43], and has previously been used to model the temporal dynamics of DENV [2,8]. From time series of disease case counts, births and total population numbers, the TSIR model constructs estimates for the susceptible population and transmission rate at discrete time points corresponding to the generation time of the disease in question. In this work, numerical implementation of the TSIR model was performed using the R package tsiR [42]. Following [43], we seek to formulate a regression model to estimate the effect of the climatological variables of average temperature and precipitation on the transmission rate β constructed by the TSIR model.

Determination of the empirical transmission rate
The empirical transmission rate β t for each generation time point by district was used as the dependent variable in a panel regression. From the output of the TSIR model, the empirical transmission rate is calculated as where I denotes the number of infected individuals, S denotes the number of susceptible individuals, N denotes the district population and α is a constant that captures heterogeneities in population mixing and the discretization of a continuous time process. S, α and the reporting rate ρ, which allows for a corrected estimate of the true case numbers, are reconstructed from the district-specific output of the TSIR model. The time period separating time points t and t + 1 is the generation time of dengue, which is taken to be two weeks [2]. For certain time steps, the estimated values of β t were undefined depending on the corresponding values of infected individuals I at the relevant time points t and t + 1. In a first instance to mitigate this effect, the transmission rate was assigned a value of β t = 0 when the number of infected individuals at the next time point was zero (I t+1 = 0). Furthermore, following the determination of the empirical transmission rate time series, undefined values of β t were replaced with either the mean value of surrounding finite values of β or simply β t−1 , depending on the position of t in the panel time series.
When constructing the panel regression model, districts for which the R 2 value associated with the linear regression of the transmission rate time series derived from the TSIR model against the empirical transmission rate time series was R 2 ≤ 0.5 were omitted (three districts). These districts (Kilinochchi, Mannar and Mullaitivu; indicated as the shaded areas in the electronic supplementary material, figure S4) share the common feature of low case numbers, and resulted in empirical transmission rates with high variance and irregular temporal dynamics.

Determination of the lag period for climate variables
Unlike the directly transmitted infection varicella considered in [43], the vector-borne nature of dengue implies that its transmission dynamics depend on prior as opposed to contemporary climatic conditions. In order to estimate the correct lag for the climate variables considered here, for each of the districts retained in the regression model, the empirical transmission rate time series were cross-correlated against each the average temperature and precipitation time series (interpolated to also be structured at time points coinciding with the generation time of dengue) for the pre-epidemic years of 2010-2016. Lag times of up to one year were considered in each cross-correlation, and only such that the climate variables led the empirical transmission rate. Using a peak finding algorithm, the local maxima and minima of the cross correlation function were subsequently located, and statistically significant values retained. Finally, from this subset, the shortest lag time corresponding to a local cross-correlation maximum or minimum was selected as the appropriate lag for that district and climate variable. These lag values are shown for all districts retained in the regression model in the electronic supplementary material, figure S4, for both average temperature (figure S4a) and precipitation (figure S4b). Note that the same three districts omitted from the panel regression model were also excluded from these calculations as a result of the irregular temporal dynamics of their empirical transmission rate time series.
Once the appropriate district-level lag value was determined for each climate variable, the national lag value was chosen as the most frequently occurring district-level lag. For the subsequent regression model, the climate variable-specific national lag value was used for all districts. In practice, for the present-day reanalysis climate data, this lag was implemented by shifting the time series in order to include data from the previous year. For the future climate projections over the two 20 year time periods, the lag was implemented by rotating the annual cycle by the appropriate delay period. In the case of average temperature, a 10 week lag was found. This value is in reasonable agreement with the 10 week lag between temperature and incidence assumed in [33], which was selected to account for a six week lag between temperature and oviposition based off findings for Aedes aegypti in Ecuador [13] as well as an assumed four week period to account for transmission, disease development, medical-care seeking, and case reporting. In the case of precipitation, a four week lag was found, which again is in reasonable agreement with the lag of three weeks between precipitation and ovitrap data in [13]. Although beyond the scope of the present work, we note that several disease transmission models explicitly include climate-driven vector life cycle dynamics, and for examples of such work the interested reader is referred to [44,45].

Development of the panel regression model
Finally, the regression panel was constructed with log(β t + 1) as the dependent variable, and the appropriately lagged climate parameters as the dependent variables, yielding the regression model where t denotes the time and i denotes the district. Location-bymonth dummies (γ i,m ) are included to remove location-specific seasonal variation in transmission which may be confounded by other seasonally varying factors. Location-by-year dummies (δ i,y ) are included to remove location-specific trends in transmission that may be spuriously correlated with climate. Standard errors were clustered at the district level. In order to determine appropriate functional forms for the climate variables in equation (4.4), non-parametric binned models of the lagged average temperature and precipitation were constructed. The results from these models with 95% confidence intervals for precipitation and average temperature are shown in the electronic supplementary material, figures S5a and S5b, respectively.
Although precipitation is critical in the mosquito life cycle for the provision of breeding sites as well as the aquatic larval and pupal stages, its quantitative effect on vector viability is less well established. By contrast, numerous studies have experimentally determined the effect of average temperature on various biological traits of mosquitoes and, when appropriate for the disease in question, parasites or viruses [6,33,[46][47][48][49]. In general, mosquito and virus traits that drive disease transmission exhibit non-monotonic thermal responses, increasing initially as a function of temperature before decreasing, following peak trait performance at a critical temperature value [33].
Motivated by this, we fitted the temperature response to a Brière function: f Bri ere (Tavg) ¼ cTavg(Tavg À T min )(T max À Tavg) 1=2 , (4:5) following [33]. Caveats to this fitted selection are discussed in §2.5. Next, we fitted equation (4.5) to the data in the electronic supplementary material, figure S5b, assuming that above the fitted maximum and minimum temperatures there is no effect of average temperature on log transmission. We impose a value of T max = 38.63°C, the average value of the fitted maximum temperature values for all A. aegpyti traits reported in the electronic supplementary material, table B of [33], and obtain parameter estimates of c = 4.99 × 10 −4 and T min = 25.27°C from this fit. The resulting curve is shown in the electronic supplementary material, figure S5c along with the filled circles corresponding to the temperature response points from the binned model. Finally, we obtain two possible regression models for log-transmission rate on the lagged climate variables of interest by considering either linear responses in both average temperature and precipitation (i.e. f (Tavg lag,t,i ) = Tavg lag,t,i and f (Precip lag,t,i ) = Precip lag,t,i ), or a Brière response for average temperature and a linear response for precipitation (i.e. f(Tavg lag,t,i ) ¼ f Bri ere (Tavg lag ) (equation (4.5)) and f (Precip lag,t,i ) = Precip lag,t,i ). The regression tables for both models are shown in the electronic supplementary material, table S1.
Because the predicted future values of average temperature generally lie below the critical thermal response temperature (see the electronic supplementary material, figures S2 and S5c), both models yield similar predictions for the effect of future average temperature on log transmission (see the electronic supplementary material, figure S6). The predictions for the district of Nuwara Eliya provide an interesting exception to this pattern, as discussed in §2.5. Ultimately, because of the more realistic prediction of the thermal response of the biological traits provided by the Brière temperature function, the first regression model of a linear response in precipitation and a Brière function for average temperature is used in the analysis of the main text.

Estimation of the fractional increase in susceptible
individuals at the onset of serotype invasion for the two-strain susceptible-infected-recovered model We consider a simplified two-strain system based on that proposed in [23], shown schematically in the electronic supplementary material, figure S1b, in order to study the possibility of a population-level susceptibility increase owing to serotype invasion. Let S denote the category of individuals susceptible to infection by both strains, I 1 and I 2 denote individuals with primary infections of serotype 1 or 2, respectively, and R 1 and R 2 denote individuals recovered from primary infections with serotypes 1 or 2, respectively. I 12 and I 21 denote individuals previously recovered from infections with serotypes 1 or 2, now experiencing a secondary infection with serotypes 2 or 1, respectively, and R denotes the class of individuals recovered from secondary infections. The total number of individuals is N. The birth and death rate are denoted by μ, and γ is the recovery rate from infection. Finally, ρ 1 and ρ 2 denote the importation rate of exogenous infections of serotype 1 and 2, respectively. Additionally, we make a few simplifying assumptions compared to the model described in [23]. First, we assume that the transmission rate β is identical for primary and secondary infections, implying that R 0 , or the number of secondary infections arising from a single primary infection in an entirely susceptible population, is identical in both cases [50]. This assumption still allows for the possibility of antibody dependent enhancement, or facilitated within-host replication of a heterologous viral strain upon secondary infection as a result of cross-reactive antibodies generated by a previous exposure [51]. Furthermore, the assumption of a uniform R 0 has been made previously in two-strain models of DENV infection [51], and it has been demonstrated that complex transmission dynamics can be achieved without asymmetries in R 0 when effects of temporary cross-immunity between heterologous strains are considered [52]. As an additional modification to the model in [23], we ignore short periods of cross-immunity as we are primarily concerned with the initial conditions at serotype invasion as opposed to the longer-term dynamics.
Prior to the re-emergence of serotype 2, the ordinary differential equations governing this mass action model are At steady state, these equations yield S* = N/β(γ + μ), At the moment of importation of ρ 2 N serotype 2 infections, the governing differential equations become Consequently, gathering terms from equations (4.9)-(4.13) associated with the generation of newly infected individuals, we see that the total rate of generation of new infections at this time point is given by which can be re-expressed as the product of the populationnormalized transmission rate and the total number of susceptible and infected (I 1 + ρ 2 N) individuals at the previous time step. Doing so we obtain b(I 1 þ r 2 N) N S þ R 1 1 À Assuming the number of imported cases is small relative to the number of autochthonous infections (i.e. ρ 2 N ≪ I 1 ), the expression in equation (4.15) simplifies to b(I 1 þ r 2 N) N S þ R 1 r 2 N I 1 : (4: 16) royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 17: 20200075