Spatiotemporal distribution of cutaneous leishmaniasis in Sri Lanka and future case burden estimates

Background Leishmaniasis is a neglected tropical vector-borne disease, which is on the rise in Sri Lanka. Spatiotemporal and risk factor analyses are useful for understanding transmission dynamics, spatial clustering and predicting future disease distribution and trends to facilitate effective infection control. Methods The nationwide clinically confirmed cutaneous leishmaniasis and climatic data were collected from 2001 to 2019. Hierarchical clustering and spatiotemporal cross-correlation analysis were used to measure the region-wide and local (between neighboring districts) synchrony of transmission. A mixed spatiotemporal regression-autoregression model was built to study the effects of climatic, neighboring-district dispersal, and infection carryover variables on leishmaniasis dynamics and spatial distribution. Same model without climatic variables was used to predict the future distribution and trends of leishmaniasis cases in Sri Lanka. Results A total of 19,361 clinically confirmed leishmaniasis cases have been reported in Sri Lanka from 2001–2019. There were three phases identified: low-transmission phase (2001–2010), parasite population buildup phase (2011–2017), and outbreak phase (2018–2019). Spatially, the districts were divided into three groups based on similarity in temporal dynamics. The global mean correlation among district incidence dynamics was 0.30 (95% CI 0.25–0.35), and the localized mean correlation between neighboring districts was 0.58 (95% CI 0.42–0.73). Risk analysis for the seven districts with the highest incidence rates indicated that precipitation, neighboring-district effect, and infection carryover effect exhibited significant correlation with district-level incidence dynamics. Model-predicted incidence dynamics and case distribution matched well with observed results, except for the outbreak in 2018. The model-predicted 2020 case number is about 5,400 cases, with intensified transmission and expansion of high-transmission area. The predicted case number will be 9115 in 2022 and 19212 in 2025. Conclusions The drastic upsurge in leishmaniasis cases in Sri Lanka in the last few year was unprecedented and it was strongly linked to precipitation, high burden of localized infections and inter-district dispersal. Targeted interventions are urgently needed to arrest an uncontrollable disease spread.

Describe where the data may be found in full sentences. If you are copying our sample text, replace any instances of XXX with the appropriate details.
If the data are held or will be held in a public repository, include URLs, accession numbers or DOIs. If this information will only be available after acceptance, indicate this by ticking the box below. For example: All XXX files are available from the XXX database (accession number(s) XXX, XXX. Leishmaniases are neglected tropical diseases caused by Leishmania parasites and transmitted 84 through the bites of vector sand flies [1]. There are three main types of leishmaniasis: visceral 85 (VL), the most serious form of the disease; cutaneous (CL), the most common; and 86 mucocutaneous (MCL), the most disabling form of the disease. Ninety two countries or 87 territories remain endemic for leishmaniasis, with more than one billion people living at risk of 88 acquiring the disease [1]. There are about 30,000 new cases of VL and >1 million new cases of 89 CL that occur annually [1]. Although nearly 100 countries are endemic to leishmaniases, major 90 transmissions are concentrated within a few. In 2018, over 90% of global VL cases were 91 reported from seven countries: Brazil, Ethiopia, India, Kenya, Somalia, South Sudan and Sudan. 92 Over 90% of MCL cases occurred in Bolivia, Brazil, Ethiopia and Peru, and 10 countries 93 reported more than 5000 CL cases [1,2]. 94 Three South Asian countries, India, Nepal, and Bangladesh, accounted for about half the global 95 burden of VL [2]. Encouraged by the recent progress achieved in leishmaniasis control, these 96 countries made a commitment to eliminate VL in the region by 2015 with the deadline extended 97 thereafter [3, 4]. The elimination target was set at less than one case per 10,000 people per 98 year, an incidence rate considered to be no longer of public health concern [3, 4]. However, Sri 99 Lanka, India's neighbor, has reported a substantial surge in clinical leishmaniasis cases in the 100 past 20 years [5]. The vast majority of leishmaniasis clinical cases in Sri Lanka are CL caused 101 by Leishmania donovani, the same parasite that causes VL elsewhere, including in India [2,6].

Clinical and climatic data collection 129
Clinical data collection has been described in a previous study [5]. Briefly, we collected clinically

Spatial distribution and temporal trend classification 141
To determine the temporal trend in leishmaniasis transmission, we used the distribution of 142 district-level annual incidence rates and hierarchical clustering to analyze the similarity in 143 incidence rates among different years (spatial similarity) [15]. The purpose of this step was to 144 analyze how transmission had evolved over the past 19 years. We used district-level incidence 145 dynamics and hierarchical clustering to determine the similarity in incidence dynamics among 146 different districts (temporal similarity) [15]. The last step was to measure space-time correlation, 147 i.e., to determine transmission synchrony [16]. This analysis will help to explain the mechanisms 148 behind transmission, e.g., spatial dispersal, locally correlated climatic variables, and/or how 149 vectors might interact with local parasites to produce different patterns of transmission. We used 150 the cross-correlation coefficient to measure region-wide transmission synchrony, and a modified 151 164 The climatic effect was tested to measure if the dynamics of infections were affected by local 165 climatic factors. The neighboring-district effect was tested to determine if the infections in a 166 given district were associated with the dispersal of parasites from neighboring districts. 167 Infections from the previous year were used to test the local infection carryover effect, i.e., the 168 localized transmission trend carried over from the previous year. Risk factor analysis was 169 conducted only in districts where mean annual incidence rate > 5 cases/1,000 population and 170 where climatic observations were available. 171 We tested each risk factor in an additive fashion; i.e., we first tested the climatic effect, then 172 added the neighboring effect, then added the carryover effect. We used stepwise analysis for 173 variable selection to determine if one factor was more significantly correlated with leishmaniasis 174 incidence dynamics than the other factors. 175

Space-time prediction 176
If localized spatiotemporal correlation existed as determined from the previous two steps, and if 177 the climatic effect was not strong, we were able to predict future case distribution based on the 178 carryover and neighboring-district dispersal effects. While this approach may underutilize the 179 climatic information, the model can predict future case numbers without climatic data; we simply 180 cannot use future climatic data (not yet available) to make predictions. Since local transmission 181 and spatial dispersal may occur at a finer scale, we used division-level data to build the  The data we used to build the model were division-level case numbers from 2015-2018. We 187 used 2019 data for validation. We used the model to predict the case distribution for 2020. 188 The maps and spatial data were generated using ArcGIS 10.0 (Redlands, CA, USA). Other data 189 analyses were conducted using R 3.6.3. 190

Leishmaniasis epidemiology: Spatiotemporal dynamics 193
Case incidences were low from 2001 to 2008, followed by a slow but steady increase in case 194 incidence from 2009 to 2017, until a sudden outbreak occurred in 2018 (Fig 1a). District-level 195 average incidence rate increased about 4-fold from 2017 to 2018. Four districts had average 196 annual incidence rates ≥ 10 cases/1,000 people, and another nine districts had annual 197 incidence rates ≥ 1 case/1,000 people (Fig 1b). 198 Districts could be classified into three groups based on similarity in incidence dynamics (Fig 3). 214 The majority of the districts with low incidence rates were in one group (gray on the map, Fig 3). 215 Four districts with intermediate incidence rates had similar temporal dynamics (yellow on the 216 map). The remaining seven districts (green on the map) were not necessarily classified into one 217 group; however, they had temporal dynamics that differed from the other two groups, and they 218 were also the seven districts with the highest average annual incidence rates (Fig 3). Incidence 219 rates in these seven districts were subjected to further risk factor analysis. The heat map in Figure 4 illustrates the distribution of incidence rates over time (Fig 4). Results 224 from the space-time joining classification indicated that the study areas could be roughly 225 classified into three groups: A) continuously low incidence until 2019 (top section); B) relatively 226 high incidence before 2008, followed by a decline and then a resurgence starting in 2012 227 (middle section); and C) continuously increasing incidence (bottom section) (Fig 4).

Risk factor analysis 236
The seven districts for which we built risk analysis models were all classified into one group 237 based on transmission dynamics (Fig 3), and they all had average annual incidence rates > 5 238 cases/1,000 people. The climatic variables-only model showed that three out of seven models 239 yielded a correlation coefficient R 2 > 0.5, indicating good correlation with climatic factors, while 240 two districts showed no significant correlation between disease dynamics and temperature or 241 precipitation variables (Table 1). Without considering the neighboring-district dispersal and 242 carryover effects, precipitation alone was not significantly correlated (P < 0.05) with 243 leishmaniasis transmission dynamics ( When the neighboring-district effect was included, six of the seven models had R 2 > 0.7 and the 253 neighboring-district effect appeared in six models (P < 0.05), precipitation appeared in four 254 models (P < 0.05), and one model still showed very weak correlation (Vavuniya District, R 2 = 255 0.18) between transmission dynamics and risk factors (Table 1). 256 When climatic, neighboring-district dispersal, and carryover effects were all included, 257 precipitation appeared in five models, neighboring-district effect in six models, and carryover 258 effect in five models. Temperature did not appear in any of the models tested (Table 1). Once 259 again, models for disease dynamics in Vavuniya District showed weak correlation (R 2 = 0.32) 260 with the risk factors investigated; only the carryover effect was selected by the stepwise 261 process, indicating likely locally isolated transmission (Table 1). 262 Overall, the risk model-predicted disease dynamics matched well with the observed dynamics 263 in six of the seven districts, all but Vavuniya (Fig 5). Leishmaniasis transmission in Vavuniya 264 was likely to be an outlier.  Table 2). Model-predicted cases matched well with observed cases for all years except 274 2018, which was an epidemic year (Fig 6 and 7, Supplemental Fig S1). The key mismatch 275 between observed and predicted cases in 2018 was the intensity in the north-central and 276 southwestern coastal areas, where the model underestimated the observations in some high-277 transmission divisions (Fig 6, Supplemental Fig S1). The model predicted both the expansion of 278 the epidemic area and the increased transmission intensity, as well as the expansion of the 279 southern hotspot from coastal to inland areas (Fig 6). The predicted case number in 2020 is 280 5,379 (95% CI [4739,6017]), which is about a 30% increase from 2019 observations (4,064 281 cases) (Fig 7). If the current trend holds, the predicted case number will be 9,115 (95% CI