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.

Introduction

Leishmaniases are diseases caused by Leishmania spp. parasites transmitted through the bites of infected female sand flies [1]. There are three main types of leishmaniasis: visceral (VL), the most serious form of the disease; cutaneous (CL), the most common; and mucocutaneous (MCL), the most disabling form of the disease [1]. One hundred and two countries or territories remain endemic for leishmaniasis, with more than one billion people living at risk of acquiring the disease [1]. There are about 30,000 new cases of VL and >1 million new cases of CL that occur annually [1]. Although more than 100 countries are endemic to leishmaniases, high transmission is concentrated within a few. In 2018, over 90% of global VL cases were reported from seven countries: Brazil, Ethiopia, India, Kenya, Somalia, South Sudan, and Sudan [1,2]. Over 90% of MCL cases occurred in Bolivia, Brazil, Ethiopia, and Peru, and 10 countries reported more than 5000 CL cases [1,2].

Three South Asian countries, India, Nepal, and Bangladesh, accounted for about half the global burden of VL [2]. Encouraged by the recent progress achieved in leishmaniasis control, these countries made a commitment to eliminate VL in the region by 2015 with the deadline extended thereafter [3,4]. The elimination target was set at less than one case per 10,000 people per year, an incidence rate considered to be no longer of public health concern [3,4]. However, Sri Lanka, India’s neighbor, has reported a substantial surge in clinical leishmaniasis cases in the past 20 years [5]. The reported cases of leishmaniasis in Sri Lanka increased from 22 cases in 2001 to 426 cases in 2010 and 3,271 cases in 2018, with one district had an incidence rate of >10 cases/10,000 people/year and eight additional districts (out of the total 25 districts) had an incidence rate of >1 cases/10,000 people/year in 2018 [5]. The vast majority of leishmaniasis clinical cases in Sri Lanka are CL caused by Leishmania donovani, the same parasite that causes VL elsewhere, including in India [2,6]. The probable vector of the Sri Lankan leishmaniasis parasite is the Phlebotomus glaucus sand fly, which demonstrates zoophilic behavior and differs from the Phlebotomus species of sand flies found in southern India, near the India–Sri Lanka strait [7,8]. The continuous upsurge of the disease transmission in Sri Lanka may hamper the regional effort of eliminating VL in the South Asia and calls for local interventions.

The question is where to start the interventions? Intervention strategies should be built based upon the understanding of the epidemiology of leishmaniasis. Leishmaniasis studies in Sri Lanka have demonstrated two transmission hotspots, one in the south coast and another in the north central region of the country, with a possible biannual seasonal variation [5,9]. Clinical leishmaniasis in Sri Lanka also showed a shift in sex and age distribution, from predominantly males aged 21–40 in early 2000s to a more balanced sex ratio and a 50/50 split between ages 21–40 and >40 years since 2015 [10]. Furthermore, outdoor-associated activities, including occupational exposure and living near a potential vector breeding site, were identified as key risk factors of leishmaniasis infection in Sri Lanka [11,12]. Previous studies have also explored the spatiotemporal patterns of leishmaniasis infections [13-16], but most of them, if not all, lack predictive power, with the output limited to distribution patterns. Other studies used agent-based modeling to predict the distribution of leishmaniasis [17,18]. These models were built based on the deeper understanding of the biology of sand fly, parasite, and host [19] and therefore, useful for the understanding of the interactions among vectors, parasites, and hosts, and how environmental (include climatic) and socioeconomic (include demographic) factors influence the distribution of leishmaniasis infections [17-19]. Since some of the variables used in the models such as climatic data is only available in retrospect, these models are good for static risk assessments and suffer similar limitations of lack of predictive power for future distributions. Since long-term leishmaniasis surveillance data are available at local scale in many countries, this may provide a good opportunity to better utilize these clinical data for future transmission predictions [5,20-22]. Time series models are useful tools for analyzing such long-term observational data [23-25], however that too is devoid of a spatial component. Nonetheless, the information revealed through past studies provide a good understanding of leishmaniasis epidemiology and transmission and pave the way for greater public awareness and disease control [5,19]. However, further questions remain regarding the analyses and predictions of spatiotemporal dynamics in leishmaniasis transmission and the driving forces behind it. For example, is the leishmaniasis parasite circulated locally or spread between neighboring areas? This information is important for determining control strategies. Furthermore, since the development of both the sand fly and the parasite inside its gut are affected by climatic conditions [26,27]; how do climatic factors affect leishmaniasis transmission? How can we better utilize this information to predict the future distribution of leishmaniasis in Sri Lanka? To answer these questions, spatial-temporal predictive models are required, which can predict the influence of potential parasite movement between neighboring areas that can facilitate spread. However, the persistence of local transmission may also lead to a buildup of case numbers. These components remain as gaps from previous work in this field [13-19,23-25].

The aim of this study was to further explore spatiotemporal leishmaniasis transmission patterns; to examine the contributions of climatic factors, local carryover of infections, and inter-region transmission to an epidemic situation; and to develop space-time models for the prediction of future clinical case distribution in Sri Lanka. This study can serve as a prototype of an early warning model for the prediction of future epidemics.

Methods

Ethics statement

Ethics approval for the study was obtained from the ethics review committee, Faculty of Medicine, University of Colombo (Reference number EC-17-062). This study only used aggregated data collected from health facilities, without any individual patient information, therefore, no consent/assent were required.

Study area

This study covered the entire Sri Lanka. It’s a country located in South Asia on the southeastern tip of the Indian sub-continent between latitudes 5° and 10° N, and longitudes 79° and 82° E. It is an island with an area of 65,610 km2 and a population of 21,803,000 (2019 estimate; 2012 census population 20,277,597) (http://www.statistics.gov.lk). Administratively, it has 25 districts and 357 divisions (vary over time due to the creation of new divisions). It consists mostly of flat to rolling coastal plains, with mountains rising in the south-central part with the highest elevation of 2,524 meters above sea level. The climate in Sri Lanka is tropical and warm with a mean temperature ranging from 17°C in the highlands to 33°C in the coastal areas. The rainfall patterns are influenced by the seasonal monsoon winds, with an annual rainfall of 2,500 mm in the central highlands down to around 1,200 mm in some of the dry zones of east, southeast, and northern parts of Sri Lanka. The annual rainfall is from 800–1,200 mm in the semiarid northwest and southeast coasts.

Clinical and climatic data collection

This study focused only on anthroponotic cutaneous leishmaniasis. Clinical data collection has been described in a previous study [5]. Briefly, clinically-confirmed leishmaniasis case data was obtained from the records maintained at the diagnostic and research laboratory at the Faculty of Medicine, University of Colombo, the Epidemiology Unit of the Ministry of Health, Sri Lanka and from small health administrative units (the medical officer of health branches) in each district. Case counts of CL at each district covered the entire country from 2001 to 2019; division-level case counts covered the period from 2015 to 2019. The population in each district was obtained from census data of 2001 and 2012 with projections for each year as given by the government of Sri Lanka, the population in each division was obtained from 2012 Sri Lanka census (http://www.statistics.gov.lk). For the district level data analysis, annual total population was used to calculate the disease incidence rate in a given district and given year from 2001 to 2019. For the division level data analysis, 2012 census projected population at each division was used to calculate incidence rate at division level for years from 2015 to 2019. The climatic data was collected from 20 meteorological stations where long-term climatic data are available in the country. Climatic variables included monthly mean of maximum, minimum, and mean temperature and monthly accumulative precipitation from 2001 to 2019. Climatic data was only used for district level disease data analysis.

Spatial distribution and temporal trend classification

The distribution of district-level annual incidence rates was used to determine the temporal trend in leishmaniasis transmission, which was analyzed using hierarchical clustering to determine the similarity in incidence rates among different years (spatial similarity) [28]. The purpose of this step was to analyze how disease transmission had evolved over the past 19 years. The district-level incidence dynamics and hierarchical clustering was used to determine the similarity in incidence dynamics among different districts (temporal similarity) [28]. The last step was to measure space-time correlation, i.e., to determine transmission spatial synchrony [29]. Spatial synchrony refers to coincident changes in the incidence or other time-varying characteristics of geographically adjacent areas [29]. Ecologically, disease spatial synchrony may arise from three primary mechanisms, i.e., dispersal among areas, congruent dependence of incidence dynamics on exogenous factors such as climatic parameters, and host-vector-parasite interactions in which vectors are spatially synchronous [29]. This analysis will help to explain the mechanisms behind disease transmission, e.g., spatial dispersal, locally correlated climatic variables, and/or how vectors might interact with local parasites to produce different patterns of transmission. The cross-correlation coefficient was used to measure region-wide transmission synchrony, and a modified Mantel correlogram to measure local transmission synchrony [29,30]. The correlations were only measured between neighboring districts and correlation with farther distance was not measured due to limited number of districts.

Risk factor analysis

This step was included to identify the driving forces behind transmission dynamics in different districts. Three types of risk factors were tested in this study: climatic variables, inter-district dispersal, and localized transmission. The following mixed regression-autoregression model incorporated with a spatial component was used to assess the risk factors [31-33]:

Where I(D) stands for incidence rate at district i in year t; ND is all neighboring districts of district i; Σ is the sum of all incidence in neighboring districts; MaxT, MeanT, MinT and Precip represent the annual average maximum, mean, and minimum temperature and the annual cumulative precipitation at district i, respectively; and α, β, γ, and β, i = 1–4, are constant coefficients.

The climatic effect was tested to measure if the dynamics of infections were affected by local climatic factors. The neighboring-district effect was tested to determine if the infections in a given district were associated with the dispersal of parasites from neighboring districts. Infections from the previous year were used to test the local infection carryover effect, i.e., the localized transmission trend carried over from the previous year. Here the carryover effect means the effect of the current infections in a given area on the future infections in the same area. In other words, if a district has high number of confirmed infections this year, one would expect a similarly high number of infections in that district by next year because parasites may well be maintained within local vectors or other reservoir host populations. Risk factor analysis was conducted only in districts where mean annual incidence rate > 5 cases/1,000 population and where climatic observations were available.

The risk factors were tested in an additive fashion; i.e., we first tested the climatic effect, with subsequent additions of the neighboring effect, and thereafter, the carryover effect. The stepwise multiple regression analysis and least square method for variable selection were used to determine if one factor was more significantly correlated with leishmaniasis incidence dynamics than the others. The correlation coefficient was used to measure the overall goodness-of-fit of the model.

Space-time prediction

If localized spatiotemporal correlation existed as determined from the previous two steps, and if the climatic effect was not strong, it was possible to predict the future case distribution based on the carryover and neighboring-district dispersal effects. While this approach may underutilize the climatic information, the model can predict future case numbers without climatic data (since future climatic data is unavailable) to make predictions. Since local transmission and spatial dispersal may occur at a finer scale, we used division-level data to build the following predictive model [33,34]:

Where, I(D) is the number of cases in division i at time t; ND represents neighboring divisions of division i; and α, β, and γ are constant coefficients. The least-square multiple regression analysis was used for model parameter estimation. However, directly observed climatic data at division level was not available for use in Sri Lanka.

The division-level case numbers from 2015–2018 were used to build the model and 2019 data was used for validation. We used the model to predict the case distribution for 2020 and total number of cases from 2020 up to 2025 together with their 95% confidence intervals.

The maps and spatial data were generated using ArcGIS 10.0 (Redlands, CA, USA). Other data analyses were conducted using R 3.6.3. with four R packages, viz. car, lme4, MASS, and synchrony.

Results

Leishmaniasis epidemiology: Spatiotemporal dynamics

There were 19,361 clinically-confirmed leishmaniasis cases from 2001 to 2019. Case incidences were low from 2001 to 2008, followed by a slow but steady increase in case incidence from 2009 to 2017, until a sudden outbreak occurred in 2018 (Fig 1A) with case numbers increasing from 1,508 cases in 2017 to 3,271 in 2018, and then 4,061 cases in 2019. District-level average incidence rate increased about 4-fold from 2017 to 2018. Four districts had average annual incidence rates ≥ 10 cases/1,000 people, and another nine districts had annual incidence rates ≥ 1 case/1,000 people (Fig 1B and 1C).

Fig 1

a) Nationwide mean annual incidence rate (cases/1,000 people/year) from 2001 to 2019 (left), b) the distribution of average incidence rate in each district over the study period (middle), and c) the corresponding names and incidence rate in each district (right). Source: Authors’ own map using data from https://data.humdata.org/dataset/sri-lanka-administrative-levels-0-4-boundaries.

Classification of district-level incidence distribution found that the study period could be classified into three based on similar distribution patterns: 2001–2010, 2011–2017, and 2018–2019 (Fig 2). All average incidence rates were < 1 case/1,000 people from 2001 to 2010, between 1 and 5 from 2011 to 2017, and >10 cases/1,000 people after 2017 (Fig 2). The similarity in incidence distributions did not follow chronological order, as indicated in branches of the classification tree (Fig 2). Overall, 2001–2010 could be seen as a low-transmission period, 2011–2017 as a parasite population build-up period, and 2018–2019 as an epidemic period (Figs 1A and 2).

Fig 2

Classification of study period by disease distribution at each district.

Districts could be classified into three groups based on similarity in incidence dynamics (Fig 3). The majority of the districts with low incidence rates were in one group (gray on the map, Fig 3). Four districts with intermediate incidence rates had similar temporal dynamics (yellow on the map). The remaining seven districts (green on the map) were not necessarily classified into one group; however, they had temporal dynamics that differed from the other two groups, and they were also the seven districts with the highest average annual incidence rates (Fig 3). Incidence rates in these seven districts were subjected to further risk factor analysis.

Fig 3

Classification of study areas by disease dynamics at each district.

Source: Authors’ own map using data from https://data.humdata.org/dataset/sri-lanka-administrative-levels-0-4-boundaries.

The heat map in Fig 4 illustrates the distribution of incidence rates over time (Fig 4). Results from the space-time joining classification indicated that the study areas could be roughly classified into three groups: A) continuously low incidence until 2019 (top section); B) relatively high incidence before 2008, followed by a decline and then a resurgence starting in 2012 (middle section); and C) continuously increasing incidence (bottom section) (Fig 4).

Fig 4

Joining classification of study areas and study period by incidence rate at each district.

The global mean correlation among district incidence dynamics was 0.30 (95% CI 0.25–0.35). However, the localized Mantel correlation between neighboring districts was on average 0.58 (95% CI 0.42–0.73), indicating a rather high correlation between neighboring districts.

The seven districts for which we built risk analysis models were all classified into one group based on transmission dynamics (Fig 3), and they all had average annual incidence rates > 5 cases/1,000 people. The climatic variables–only model showed that three out of seven models yielded a correlation coefficient R2 > 0.5, indicating good correlation with climatic factors (R2 values 0.711–0.829), while two districts showed no significant correlation between disease dynamics and temperature or precipitation variables (Table 1). Without considering the neighboring-district dispersal and carryover effects, precipitation alone was not significantly correlated (P < 0.05) with leishmaniasis transmission dynamics (Table 1).

Table 1

Risk factors of clinical leishmaniasis in different districts.

District	Hambantota	Matara	Kurunegala	Polonnaruwa	Matale	Anuradhapura	Vavuniya
Mean annual incidence rate (cases/ 1,000 population)	41.817	19.264	6.972	21.779	7.505	23.532	9.550
Climate model
Adjusted R2	0.829	N/A	N/A	0.784	0.425	0.711	0.185
ANOVA F- value	30.127			22.812	7.648	23.196	5.079
P-value	<0.0001			<0.0001	0.0047	<0.0001	0.0377
Coefficients
Minimum temperature	-68.499***					-78.453***	6.547*
Mean temperature	64.465 **			83.635***	64.183**
Maximum temperature	-35.656 *			-37.536 **	-39.912 *	69.320***
Precipitation
Climate + Neighbor-district model
Adjusted R2	0.838	0.808	0.873	0.930	0.832	0.726	0.185
ANOVA F- value	47.675	38.967	62.606	119.738	45.612	24.861	5.079
P value	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001	0.0377
Coefficients
Minimum temperature	-35.213 **	22.311***					6.547*
Mean temperature
Maximum temperature
Precipitation			-0.005 **	0.004 *	-0.009***	0.014 **
Neighboring-district mean	3.340 ***	1.262***	1.815***	3.061***	1.551***	2.444***
Climate + Neighbor-district+Carryover model
Adjusted R2	0.915	0.846	0.937	0.930	0.832	0.816	0.316
ANOVA F- value	65.449	33.858	89.557	119.738	45.612	40.796	9.310
P-value	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001	<0.0001	0.0072
Coefficients
Minimum temperature
Mean temperature
Maximum temperature
Precipitation	0.015*	-0.009*	-0.003*	0.004*	-0.009*

Note: Stepwise regression analysis was used, insignificant variables were removed at level of 0.05

*, ** and*** represent significant level of <0.05, <0.01 and <0.001, respectively.

When the neighboring-district effect was included, six of the seven models had R2 > 0.7 (R2 values 0.726–0.930) and the neighboring-district effect appeared in six models (P < 0.05), precipitation appeared in four models (P < 0.05), and one model still showed very weak correlation (Vavuniya District, R2 = 0.18) between transmission dynamics and risk factors (Table 1).

When climatic, neighboring-district dispersal, and carryover effects were all included, precipitation appeared in five models, neighboring-district effect in six models, and carryover effect in five models (R2 values 0.816–0.930). Temperature did not appear in any of the models tested (Table 1). Once again, models for disease dynamics in Vavuniya District showed weak correlation (R2 = 0.32) with the risk factors investigated; only the carryover effect was selected by the stepwise process, indicating probable locally isolated transmission (Table 1).

Overall, the risk model–predicted disease dynamics matched well with the observed dynamics in six of the seven districts, except Vavuniya (Fig 5 and Table 1). Leishmaniasis transmission in Vavuniya was likely to be an outlier.

Fig 5

Observed and model-predicted temporal changes in incidence rate in the seven districts with the highest case numbers.

Climate, neighbor, and carryover represent climatic, neighboring-district dispersal, and carryover effects, respectively. The results of goodness-of-fit of each model were presented in Table 1, predictions was only made with significantly fitted models.

Spatial distribution and future trend prediction

Using cases from the previous year in the same division and in neighboring divisions, we built a prediction model with R2 = 0.823 and adjusted R2 = 0.823 (ANOVA F2,1321 = 3075.71, P < 0.0001, Table 2). Model-predicted cases matched well with observed cases for all years except 2018, which was an epidemic year (Figs 6 and 7 and S1). The key mismatch between observed and predicted cases in 2018 was the intensity in the north-central and southwestern coastal areas, where the model underestimated the observations in some high-transmission divisions (Figs 6 and S1). The model predicted both the expansion of the epidemic area and the increased transmission intensity, as well as the expansion of the southern hotspot from coastal to inland areas (Fig 6). The predicted case number in 2020 is 5,379 (95% CI [4739, 6017]), which is about a 30% increase from 2019 observations (4,064 cases) (Fig 7). If the current trend holds, the predicted case number will be 9,115 (95% CI [8,157, 10,070]) by year 2022 and close to 20,000 by year 2025 (19,212 cases, 95% CI [17,389, 21,31]).

Table 2

Division-level predictive model summary.

Variable	Coefficient (95% CI)	t-value	p value
Intercept	0.8102 [0.1921, 1.4284]	2.57	0.0102
Cases from previous year	1.0679 [1.0273, 1.1086]	51.52	< 0.0001
Cases in neighboring divisions from previous year	0.0331 [0.0215, 0.0446]	5.61	< 0.0001

Fig 6

Observed and model-predicted leishmaniasis cases by division and year.

Source: Authors’ own map using data from https://data.humdata.org/dataset/sri-lanka-administrative-levels-0-4-boundaries.

Fig 7

Observed vs. model-predicted total cases from 2015 to 2025.

Discussion

Leishmaniasis clinical cases are still on the rise as reported by the Sri Lankan Ministry of Health, in contrast to the substantial decline in case numbers recorded in other South Asian countries including India, Sri Lanka’s nearest neighbor [2,5]. Although leishmaniasis clinical cases were first reported from Sri Lanka in 1992 [35], the case numbers have remained low and sporadic prior to 2001 [36]. Reported case numbers increased 15-fold from 2001 to 2010, followed by a 10-fold increase from 2010 to 2019, with a clear outbreak in 2017–2018. In addition to localized transmission, inter-area dispersal seems to have played an important role in spreading the parasite, and precipitation is probably the key climatic factor fueling the increasing trends in infections. The spatiotemporal model predicted that cases will continue to rise unless viable interventions are implemented. There are no formal leishmaniasis interventions implemented either locally or nationwide in Sri Lanka, except for diagnosis and treatment of patients who report to hospitals. With due consideration given to the highly localized transmission with spread among neighboring areas, introduction of focal vector control measures may help to effectively reduce the local transmission and parasite spread, since factors such as the climate, population movement and vector dispersal are beyond control., The potential risk of uncontrollable epidemic calls for urgent introduction of effective intervention strategies.

Spatiotemporal analysis is a powerful tool, not only for detecting transmission hotspots but also for understanding the driving forces behind transmission, e.g., how spatial dispersal and locally correlated climatic variables contribute to transmission [11,28,29,37-39]. The spatiotemporal and risk factor analyses of epidemiological studies conducted in the past have yielded limited outcome [11,37-42] and they lacked predictive power [43,44] and failed to optimally utilize spatial information [20-22]. The multi-step approach adopted in this study included not only spatiotemporal clustering but also spatiotemporal risk modeling. Study of disease dynamics in Sri Lanka based on spatiotemporal classification analysis described in this paper demonstrated that the study period from 2001 to 2019 may be divided into three phases: the low-transmission phase from 2001 to 2010, the parasite population build-up phase from 2011 to 2017, and the sudden epidemic starting in 2018. From 2001 to 2010, parasites circulated locally with a few hotspots [5]. Prior to 2017, no more than three districts had relatively high incidence rates (>50 cases/1,000 people-years) each year. However, this number increased to five districts in 2018 and seven in 2019. Spatial correlation analysis indicated that transmission dynamics between neighboring districts were positively correlated, suggesting inter-area spreading. Risk factor analysis in the heavily hit districts suggested that precipitation was likely the key climatic factor driving the increase in cases. Persisting CL lesions due to poor responsiveness to routine therapy [45] may have contributed to the sudden leishmaniasis surge in Sri Lanka in 2018. However, further investigations will be required for its confirmation. The results from this study provide a better understanding of the epidemiology of leishmaniasis and the factors related to the emergence of the disease. The temporal prediction is useful for epidemic preparedness, and the spatial prediction can help to guide targeted interventions for control.

Leishmaniasis transmission in Vavuniya District was a special case. Vavuniya and Kilinochchi districts were in the front lines of the 30-year armed conflict in Sri Lanka, which ended in 2009 [45,46]. With the existence of indigenous circulating leishmaniasis parasites, the transmission of leishmaniasis in Vavuniya District seemed focal, and might be due to the displacement of locals during the conflict and the postwar return of individuals who had fled, coupled with the influx of armed forces personnel who assisted with resettlement and reconstruction of infrastructure [46,47]. It may be described as a test cage, where parasites circulated freely. Once resettlement was completed, local infections may have been mass tested, clinical cases treated and armed forces personnel redeployed elsewhere [46,47], which may have influenced the control of parasite transmission. Such factors may explain the failure to reliably predict the transmission dynamics in this district in relation to the climatic variables and neighboring-district dispersal.

A limitation of this study is the spatiotemporal predictive model, which lacked the power to predict the 2018 leishmaniasis outbreak. Since long-term dynamics of climatic variables cannot be predicted in a futuristic manner, any climate-based model will lack predictive power. Climatic factors were not included in the predictive model in this study; however, if climatic factors were indeed the driving force behind the 2018 outbreak, model predictions excluding climatic variables would surely miss the epidemic. Studies in other countries have linked climatic variability with outbreaks of leishmaniasis [23-25]. Using long-term mean climatic data to predict future disease transmission is an option. However, it lacks variability and thus likely to miss the abnormal changes; in addition, it creates another dimension of uncertainty for the predictions. Nonetheless, we speculate that climatic variation and/or persistent CL lesions due to poor drug response may have been the causes of the 2018 outbreak in Sri Lanka; however, further evidence is required. Another limitation is the short temporal period used for building the predictive model (the model without climatic variables). We collected disease data from some 350 divisions for each year, which provided a good sample size in terms of the goodness-of-fit of the model; still, four years is a short period. Adding data observed from other years might improve model performance, although collecting accurate long-term nationwide division-level clinical case data has proven difficult. There are other minor limitations in this study, for example, the neighboring area effect and annual climatic data. Since monthly confirmed case numbers in each district are available since 2014, using monthly climatic data may fine tune the modeling results at least to a certain degree. At divisional level, infection may spread beyond the neighboring divisions because some divisions are small in size. Although the second-order neighborhood or distance weighted neighborhood effect could have been analyzed, this may have diverged the focus of the study.

In conclusion, leishmaniasis transmission in Sri Lanka is still spreading and has recently intensified. Spatiotemporal analysis incorporating a neighboring-area dispersal effect provides a new framework for analyzing the spatiotemporal dynamics of the disease, detecting the risk factors that affect the dynamics, and predicting potential future distribution of infections. Such information on future transmission distribution predictions are useful for disease control programs and policy makers for planning and prioritizing of targeted areas.

Correlation between observed and predicted cases at each division from 2016 to 2019.

(TIF)

Click here for additional data file.

2 Feb 2021

Dear Prof. Karunaweera,

Thank you very much for submitting your manuscript "Spatiotemporal transmission of cutaneous leishmaniasis in Sri Lanka and future case burden estimates" for consideration at PLOS Neglected Tropical Diseases. As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Alberto Novaes Ramos Jr

Associate Editor

PLOS Neglected Tropical Diseases

Epco Hasker

Deputy Editor

PLOS Neglected Tropical Diseases

***********************

Reviewer's Responses to Questions

Key Review Criteria Required for Acceptance?

As you describe the new analyses required for acceptance, please consider the following:

Methods

-Are the objectives of the study clearly articulated with a clear testable hypothesis stated?

-Is the study design appropriate to address the stated objectives?

-Is the population clearly described and appropriate for the hypothesis being tested?

-Is the sample size sufficient to ensure adequate power to address the hypothesis being tested?

-Were correct statistical analysis used to support conclusions?

-Are there concerns about ethical or regulatory requirements being met?

Reviewer #1: Yes

Reviewer #2: my comments are attached.

Reviewer #3: Partially.

Reviewer #4: Major Revision

Reviewer #5: Please mention from where did you obtain meteorological data.

How many clinically confirmed cases you included for this study?

What are the tests you applied to analyze risk factors

--------------------

Results

-Does the analysis presented match the analysis plan?

-Are the results clearly and completely presented?

-Are the figures (Tables, Images) of sufficient quality for clarity?

Reviewer #1: Yes

Reviewer #2: my comments are attached.

Reviewer #3: Yes

Reviewer #4: Minor Revision

Reviewer #5: fig 1 (b) incidence rate - which year? ?

Add p values for each clinical factors in every districts

--------------------

Conclusions

-Are the conclusions supported by the data presented?

-Are the limitations of analysis clearly described?

-Do the authors discuss how these data can be helpful to advance our understanding of the topic under study?

-Is public health relevance addressed?

Reviewer #1: Yes

Reviewer #2: my comments are attached.

Reviewer #3: Yes

Reviewer #4: Major Revision

Reviewer #5: conclusions supported by the data presented

limitations of analysis clearly described

--------------------

Editorial and Data Presentation Modifications?

Use this section for editorial suggestions as well as relatively minor modifications of existing data that would enhance clarity. If the only modifications needed are minor and/or editorial, you may wish to recommend “Minor Revision” or “Accept”.

Reviewer #1: Minor Revision

Reviewer #2: my comments are attached.

Reviewer #3: (No Response)

Reviewer #4: (No Response)

Reviewer #5: (No Response)

--------------------

Summary and General Comments

Use this section to provide overall comments, discuss strengths/weaknesses of the study, novelty, significance, general execution and scholarship. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. If requesting major revision, please articulate the new experiments that are needed.

Reviewer #1: In this study the authors addressed spatiotemporal transmission of cutaneous leishmaniasis and predicted future case burden in Sri Lanka. The objective of this study is clear, the data are well-analyzed, and the results are well-presented. This analytic model will be a powerful tool for the prediction of future risk of leishmaniasis, and contribute to the control of this disease I enjoyed reading this manuscript. Following is a minor comment.

1. The data of figure 5 is not explained well in the text. I recommend mentioning about it more detail in the text, or move this figure into a supplemental figure.

Reviewer #2: my comments are attached.

Reviewer #3: The paper ‘Spatiotemporal transmission of cutaneous leishmaniasis in Sri Lanka and future case burden estimates’ covers a very relevant and urgent topic in the Sri Lank context.

Below you can find recommendations to further improve this investigation:

L99 – Please provide a detailed information about this upsurge in leishmaniasis cases. This is critical to justify the need for the present study.

L128-138 – The sections about the study area and data sources is very incomplete. The following information needs to be included:

a. Description of the territorial units of analysis, more precisely population size (mean and SD or median and IQR)

b. Description of the information systems that collect and provide data on leishmaniosis cases (flow of information from diagnosis until register in the database; coverage; existence of changes in data quality through time, etc.)

c. Description of the climatic data, more precisely the data source, temporal resolution and number of meteorological stations in the climate network from Sri Lanka.

If data sources are prone to any type of information and selection bias, that needs to be discussed in the limitations section of the manuscript.

L148 – Please clarify the concept of transmission synchrony. The readers may not know this term.

L190 – Please indicate the R packages used for the analysis.

L180 – The authors do not use climatic variable to make the space-time predictions. However, the authors found that climatic factors played a role in the infections. Could it be possible to make predictions based on meteorological scenarios (high, low, medium precipitation year; high, low, medium temperature year)? I believe it would be an interesting add-on to the present or future research.

Figure 1 – To make the figure more informative I recommend numbering the districts and add the names as a legend of the figure.

Reviewer #4: (No Response)

Reviewer #5: Don't start sentences as "we.................. Please write in a scientific way.

Please give your suggestions more about how a high correlation between neighboring districts and incidence rate of leishamaniasis in endemic regions.

--------------------

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: Yes: Ana Isabel Ribeiro

Reviewer #4: No

Reviewer #5: Yes: Lahiru Sandaruwan Galgamuwa

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org.

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, PLOS recommends that you deposit laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see https://journals.plos.org/plosntds/s/submission-guidelines#loc-methods

Submitted filename: PNTD-D-20-02228_reviewer-1.pdf

Click here for additional data file.

Submitted filename: Rewiever Comments.pdf

Click here for additional data file.

12 Mar 2021

Submitted filename: Responses to reviewers comments.docx

Click here for additional data file.

30 Mar 2021

Dear Prof. Karunaweera,

We are pleased to inform you that your manuscript 'Spatiotemporal distribution of cutaneous leishmaniasis in Sri Lanka and future case burden estimates' has been provisionally accepted for publication in PLOS Neglected Tropical Diseases.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Neglected Tropical Diseases.

Best regards,

Alberto Novaes Ramos Jr

Associate Editor

PLOS Neglected Tropical Diseases

Epco Hasker

Deputy Editor

PLOS Neglected Tropical Diseases

***********************************************************

Reviewer's Responses to Questions

Key Review Criteria Required for Acceptance?

As you describe the new analyses required for acceptance, please consider the following:

Methods

-Are the objectives of the study clearly articulated with a clear testable hypothesis stated?

-Is the study design appropriate to address the stated objectives?

-Is the population clearly described and appropriate for the hypothesis being tested?

-Is the sample size sufficient to ensure adequate power to address the hypothesis being tested?

-Were correct statistical analysis used to support conclusions?

-Are there concerns about ethical or regulatory requirements being met?

Reviewer #1: Yes

Reviewer #3: The authors addressed my previous concerns in the revision

Reviewer #4: Yes

Reviewer #5: yes

**********

Results

-Does the analysis presented match the analysis plan?

-Are the results clearly and completely presented?

-Are the figures (Tables, Images) of sufficient quality for clarity?

Reviewer #1: Yes

Reviewer #3: The authors addressed my previous concerns in the revision

Reviewer #4: Yes

Reviewer #5: yes

**********

Conclusions

-Are the conclusions supported by the data presented?

-Are the limitations of analysis clearly described?

-Do the authors discuss how these data can be helpful to advance our understanding of the topic under study?

-Is public health relevance addressed?

Reviewer #1: Yes

Reviewer #3: (No Response)

Reviewer #4: Yes

Reviewer #5: yes

**********

Editorial and Data Presentation Modifications?

Use this section for editorial suggestions as well as relatively minor modifications of existing data that would enhance clarity. If the only modifications needed are minor and/or editorial, you may wish to recommend “Minor Revision” or “Accept”.

Reviewer #1: (No Response)

Reviewer #3: no

Reviewer #4: Accept

Reviewer #5: Accept

**********

Summary and General Comments

Use this section to provide overall comments, discuss strengths/weaknesses of the study, novelty, significance, general execution and scholarship. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. If requesting major revision, please articulate the new experiments that are needed.

Reviewer #1: The paper was revised properly.

Reviewer #3: The authors addressed my previous concerns in the revision

Reviewer #4: I have read with interest your manuscript as well as reviewers’ concerns and your response. The revisions that you made to the manuscript are very effective in addressing the remaining concerns. The revised manuscript is overall well-written, well-structured and clear. So, in my opinion, the manuscript is well-suited for publication in PLOS Neglected Tropical Diseases.

Reviewer #5: Accept

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #3: Yes: Ana Isabel Ribeiro

Reviewer #4: No

Reviewer #5: Yes: Lahiru Sandaruwan Galgamuwa

14 Apr 2021

Dear Prof. Karunaweera,

We are delighted to inform you that your manuscript, "Spatiotemporal distribution of cutaneous leishmaniasis in Sri Lanka and future case burden estimates," has been formally accepted for publication in PLOS Neglected Tropical Diseases.

We have now passed your article onto the PLOS Production Department who will complete the rest of the publication process. All authors will receive a confirmation email upon publication.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any scientific or type-setting errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript. Note: Proofs for Front Matter articles (Editorial, Viewpoint, Symposium, Review, etc...) are generated on a different schedule and may not be made available as quickly.

Soon after your final files are uploaded, the early version of your manuscript will be published online unless you opted out of this process. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting open-access publishing; we are looking forward to publishing your work in PLOS Neglected Tropical Diseases.

Best regards,

Shaden Kamhawi

co-Editor-in-Chief

PLOS Neglected Tropical Diseases

Paul Brindley

co-Editor-in-Chief

PLOS Neglected Tropical Diseases