Spatial-time Visualization of Cutaneous Leishmaniasis in an Urban Area: Informing Health Policy and Practice

Cutaneous leishmaniasis (CL) is an important public health concern worldwide. Iran is among the most CL-affected countries, being listed as one of the rst six endemic countries in the world. In order to develop targeted interventions, we performed a spatial-time visualization of CL cases in an urban area to identify high-risk and low-risk areas during 2016-2019. Methods: This cross-sectional study was conducted in the city of Mashhad. Patient data were gathered from Mashhad health centers. All cases (n=2425) were diagnosed in two stages; the initial diagnosis was based on clinical ndings. Subsequently, clinical manifestation was conrmed by parasitological tests. The data were aggregated at the neighborhood and district levels and smoothed CL incidence rates per 100,000 individuals were calculated using the spatial empirical Bayesian approach. Furthermore, we used the Anselin Local Moran’s I statistic to identify clusters and outliers of CL distribution during 2016-2019 in Mashhad.

The southwestern area of Mashhad had the highest CL incidence rates. This piece of information might be of value to design tailored interventions such as running effective resource allocation models, informed control plans and implementation of e cient surveillance systems. Furthermore, this study generates new hypotheses to test potential relationships between socio-economic and environmental risk factors and incidence of CL in areas with higher associated risks. Background Leishmaniasis, a serious worldwide public health problem caused by intracellular protozoan parasites, was prevalent in 92 countries in 2018 with an estimated 700,000 to 1 million new annual infections (1,2).
It is a vector-borne disease with a wide range of symptoms, which is transmitted to humans by phlebotomine sand ies (3). There are different variants of this disease, the most common of which is the cutaneous form. Over 85% of new identi ed cutaneous leishmaniasis (CL) cases originated in Afghanistan, Algeria, Brazil, Colombia, Iran (Islamic Republic of), Pakistan, Peru, and the Syrian Arab Republic in 2018 (2,4). Recently, the number of CL cases has signi cantly increased in the endemic countries (5). This increase is attributed to changes in both natural and human-made environments such as fast and unplanned urbanization, increase in agricultural development, large migrations, and deforestation (6,7). CL causes psychological, social, and economic problems (8), and is a major public health challenge in the affected areas.
Iran is among the top six countries in the world with the highest annual rate of CL incidence (2,4), with 50-250 cases per 100,000 individuals (9). CL, the second most prevalent arthropod-transmitted disease (10), is a complex public health problem and its distribution follows a spatial pattern in Iran (2,4,11). The Khorasan-Razavi province, located in northeastern Iran, has the highest prevalence of CL in the country (12,13). Mashhad, the capital of Khorasan-Razavi province, is the second most populous city in Iran and is bordered by two neighboring CL-prevalant countries, Turkmenistan and Afghanistan. Moreover, Mashhad, one of the main CL endemic cities (14,15), has different CL incidence rates across different geographical areas (16).
Identi cation of geographical distribution and spatio-temporal patterns of disease can help to implement effective preventive strategies. Geographic information system (GIS) has been previously applied to visualize the geographical distributions of disease and the spatial modelling of risk factors, and in particular environmental factors (17,18). GIS is a powerful tool to gather and analyze spatial and nonspatial data regarding the epidemiology of disease and understanding its underlying ecology. Spatial and spatial-time analyses can provide insight into disease occurrence and guide potential tailored disease control interventions. As such, identifying high-risk areas can be useful to implement surveillance programs and come up with effective disease control strategies (19,20). Because environmental factors can potentially impact the geographical distribution of leishmaniasis, using spatial analysis can foster valuable insight for health policymakers (21).
According to previous studies, GIS-based techniques led to the detection of countries with the highest rate of CL cases between 2001 and 2011 in Latin America (22). In a study conducted in Brazil, the areas with high CL incidence rate were identi ed using spatial analysis (23). A number of studies used spatial analysis to assess spatial distribution and pattern of CL in some areas of Iran at provincial level (24)(25)(26)(27)(28).
Furthermore, other studies have attempted to determine relations between epidemiological and environmental factors and CL incidence in endemic areas (24,27,29). A few studies have been carried out in Khorasan provinces (North, Razavi and South) to assess the spatial distribution of CL at provincial and national levels. In these studies, southern areas of North Khorasan province and northern areas of Razavi Khorasan province were identi ed as hot areas (28,30). Another study revealed the high prevalence of CL in southwestern areas of Mashhad using molecular methods (14). An epidemiological study assessed CL prevalence in Mashhad over a period of 20 years (1992-2014) and identi ed climaterelated factors with remarkable in uence on the prevalence of CL (31). The previous studies con rmed the high prevalence of CL in Mashhad, but to the best of our knowledge, no study has investigated the spatial-time pattern of CL incidence at district and neighborhood levels in Mashhad. This underscores the need to explore the geographical distribution and high-risk disease areas. Exploring high-risk areas is an important step in developing strategies to effectively allocate healthcare resources, so as to overcome the disease and its associated complications (32). The objective of the present study is to apply a spatial analysis of CL cases in the city of Mashhad between 2016 and 2019, and to identify the spatial-time pattern of CL cases and identify the critically affected areas.

Study area
This cross-sectional study was conducted in the city of Mashhad, located in northeastern Iran ( Fig. 1) with an approximate population of 3,000,000. Mashhad is situated at a latitude of 35'N − 37'N and longitude of 59'E − 60'E, with a total area of 328 km 2 .

Cutaneous leishmaniasis de nition and diagnosis methods
Cutaneous leishmaniasis is the most common form of leishmaniasis, which causes skin ulcers and serious disabilities persisting for months and sometimes for years. The skin lesions develop in the exposed parts of the body (2,33). CL is caused by parasites, transmitted by the bite of phlebotomine sand ies, and affects humans (3). Diagnosis of CL is through a combination of clinical symptoms and laboratory validation (2). In this study, all cases were diagnosed in two stages. This includes an initial diagnosis based on clinical examination followed by a validation using parasitological tests.

Data sources
Data were obtained from different sources at different spatial scales. Individual CL data were obtained from ve Mashhad health Centers, including patient data collected between 2016 and 2019. In the past years, CL data were generally gathered at the municipal level. However, during the study period, data collection was performed at district and neighborhood scales. The base maps for conducting spatial analyses were obtained from Mashhad municipality.

Data preparation
All patients' addresses were geocoded manually using Google's MyMaps software (https://www. google.com/maps/about/mymaps). In order to generate the point map ( Fig. 1) and protect patients' privacy, we used a geomasking technique. Each patient's address was randomly assigned to a location within a radius of half a kilometer around their actual home address.

Descriptive mapping
To generate descriptive maps, we used the natural break classi cation with ve classes. This classi cation method produces the best arrangement of values into different classes by minimizing the class's average deviation from the class mean. This method aims to reduce the variance within classes and increase the variance between classes (34).

Empirical Bayesian Smoothing
The relevant population at risk typically varies across the areas, so the precision of raw incidence rate will vary as well. This variance instability requires the smoothing of rates. The main objective in CL mapping was to estimate the incidence rate for each geographical unit using the observed number of patients.
Empirical Bayes is one of the spatial statistical smoothing methods which is used to reduce the random uctuations due to population size. In this model, the estimation of underlying rates is attained by making an assumption about prior disease distribution using observed data (35)(36)(37). The empirical Bayes technique computes the risk as a weighted sum of the raw rate for each unit and a prior mean (38,39). The incidence rate of CL during 2016-2019 was adapted for smoothing by the Bayesian models.

Spatial cluster analysis
Mashhad has 172 neighborhoods and 13 districts. The population of all neighborhoods and districts was obtained from the Statistics Center of Iran. The Local Moran's I statistic (40) was performed to identify potential spatial-time clusters in terms of CL incidence rate during the 2016-2019 period at the neighborhood level. We used two incidence rates for spatial clustering: raw incidence rate and empirical Bayesian incidence rate. We computed the leishmaniasis raw incidence rate for each geographical unit by dividing leishmaniasis cases in each unit by its total population at risk in that unit. The use of local Moran's I clustering method has been suggested in many studies (24,41,42). This model calculates a zscore and p-value to determine whether the apparent similarity (a spatial clustering of either high or low values) or dissimilarity (a spatial outlier) is more pronounced than the one expected in a random distribution. The null hypothesis states that CLs are randomly distributed across the study area. In this method, the areas with similar high or low CL rates are indicated by High-High (HH) and Low-Low (LL) areas, while the regions with dissimilar rates are described as High-Low (HL) and Low-High (LH) regions (43). In other words, the HH and LL areas indicate the presence of CL clusters, while HL and LH areas indicate outliers.

Moran Scatter plot
The Moran scatter plot consists of a plot with the original variable and spatially lagged variable. In this plot the spatial autocorrelation is classi ed into four categories that is visualized as four quadrants.
Positive spatial autocorrelation (similar values in nearby or neighboring areas), known as HH and LL areas, are shown in the upper-right quadrant and the lower-left quadrants, respectively. In contrast, the lower-right and upper-left quadrants are related to negative spatial autocorrelation (unlike neighboring areas), known as HL and LH, respectively. The slope of scatter plot's line corresponds to Moran's I as the most common indicator of spatial autocorrelation. In this plot, the original variable has been standardized (zero mean and unit standard deviation). Also, the spatial lag variable is considered for those standardized values (44). In this study, the original variable was CL incidence and the spatially lagged variable was lagged incidence.  The temporal uctuation of CL-positive cases is shown in Fig. 2 The highest incidence rates of CL cases at district levels were concentrated in the northwestern and southwestern areas of the city from 2016 to 2019. The reduction of CL incidence during the three years is clearly seen in the aforementioned areas (Fig. 3).
The descriptive maps of CL crude incidence at the neighborhood level revealed that the southwestern area of the city had the highest incidence rate. The local Moran's I revealed a statistically signi cant spatial autocorrelation regarding CL incidence at the neighborhood level. The cluster maps revealed HH clusters of CL in the southwestern area from March 2016 to March 2019. Furthermore, there was an LL cluster in the northeastern area of the city from March 2016 to March 2019 (Fig. 4). All clusters and outliers were statistically signi cant (p < 0.05).
As outlined in the Methods section, rates have an intrinsic variance instability, which may lead to the identi cation of spurious outliers. In order to correct for this, we can use smoothing approaches (also called shrinkage estimators), which improve the precision of crude rate by borrowing strength from other observations. The descriptive maps of smoothed incidence CL cases revealed that the southwestern area of the city was the area most affected by CL incidence. Also, the cluster maps con rm these results and show the HH clusters of CL cases in the southwestern area of city from March 2016 to March 2019. Furthermore, there were LL clusters in the northeastern area of the city from March 2016 to March 2019, which had the smallest areas with high priority (Fig. 5).
The rst scatter plot shows the CL incidence for 2016-2017 with a Moran's value of 0.236 that should be interpreted as a spatial pattern with a cluster tendency (HH and LL). The second and third scatter plots revealed the CL incidence by Moran's value of 0.296 and 0.231 during 2017-2018 and 2018-2019, respectively. On the lower end of the spectrum, there is a much smaller spread in the CL incidence over the three-year period (Fig. 6).

Discussion
The main aim of this study was to explore the spatial-time patterns of CL incidence using two spatial analysis techniques: local Moran's clustering analysis on crude incidence rate and spatial empirical Bayesian smoothed rate. To the best of our knowledge, this is the rst study in the city of Mashhad to assess spatial patterns of CL incidence using these two approaches at the nest geographical level (neighborhoods). In the present study, our objective was to understand the spatial and temporal trends of CL. Hence, we evaluated the incidence rates of CL, and designed risk maps for the 2016-2019 period using GIS techniques. The results identi ed the southwestern areas of the city as high-risk areas.
Furthermore, northeastern areas were identi ed as low-risk CL areas. Also, temporal variation of CL revealed that CL incidence rates have steadily declined during the study period. The declining trend of CL is also reported in a previous study, which assessed CL incidence during the last four years in Mashhad (45). An important factor that may have contributed to the steady decline of CL infection rates over last years is the new municipal waste management system. This is further explained by the fact that this new system has facilitated the collection of construction waste in urban areas, which is an attraction site for sand ies breeding. Also, a recent change in climate in Mashhad, which has resulted in heavy rains and oods, may have been another contributing factor in the decline of CL infections. This is because the heavy precipitation can wash away sand ies' eggs and larva.
Khorasn-Razavi province is an area of long-standing hyper-endemic CL with a reported prevalence rate of 64/100,000 in 2013 (46). Mashhad, the province capital, is one of the main CL endemic areas with a high prevalence rate (45). Despite previous attempts to prevent and control leishmaniasis, new cases of CL are reported every year, highlighting the need for further interventions (47). The rst step in the design of CL control and prevention programs is the identi cation of high-risk areas. Taking advantage of GIS techniques and spatial analysis plays an important role in the identi cation and visualization of highly affected regions. Previous studies assessed the spatial distribution of CL in some parts of Iran (24,25,27,29) and also other endemic countries (22,23), in an attempt to identify areas of high priority. In some studies, the clustering method has been used for spatial analysis of CL incidence in different countries and provinces (23,24,27,48), while other studies have used different approaches such as geographic weighted regression (GWR), ordinary least square (OLS) regression, bivariate and multivariate logistic regression to evaluate the relative environmental factors (29,49). The present study revealed a high incidence rate of CL in southwestern Mashhad using both crude incidence rate and spatial empirical Bayes rate analysis.
Previous epidemiological studies in Mashhad revealed the regions with high incidence. The Ab-o-Bargh region had the highest prevalence of CL in 2001. This nding is in line with the results of the present study. Some socio-epidemiological risk factors were identi ed for high incidence in this region, which include density of population, stray dogs, and some cultural and economic factors (13). Some studies also revealed a high frequency of CL in the cities of Neishabour (50) and Torghabeh-Shandiz (16). Both of these cities are located to the west of Mashhad, and these results reinforce the ndings of the present study. A previous study indicated that the high frequency of CL was due to deforestation or human migration (51). Our ndings strengthen this hypothesis since there has been fast and uncontrolled migration to southwestern areas of Mashhad, where we found a high rate of CL incidence. Migration is potentially a contributing factor since it can provide the means for disease transmission. Another important factor, which was also identi ed in the previous studies and contributes to the increase in CL incidence, is urbanization (52). Human activities such as construction of dams, roads, and new residential complexes can cause high CL incidence rates. Also, living in the marginalized areas of the city, which suffer from poor health services, can be regarded as another risk factor for CL (31,53).
Ab-o-Bargh and Torghabeh-Shandiz regions are mountainous areas. As such, they are home to rodent nests and carnivores, which are the main reservoirs of leishmaniasis. Moreover, the two aforementioned neighborhoods are humid, which provide a contributing environment for sand y breeding (54). These may be regarded as considerable environmental factors that can be taken into account for future implementation of control strategies.

Limitations
There are some limitations for this study. Iran has a passive system for CL surveillance, which may contribute to underreporting of CL cases. Also, there may exist cases with small lesions who might not refer to health centers. This can affect the number of diagnosed and reported cases in some areas. Furthermore, the contribution of social and environmental factors was not investigated in this study, which can be the topic of future studies. Despite these limitations, the data used in this study were su ciently robust for spatial-time analysis of CL cases, and for identifying the priority areas for surveillance in Mashhad.

Conclusions
The present study provides a vivid depiction of the spatial-time spread of CL, with most cases concentrated in the southwestern region of Mashhad during the three-year study period. In our study, high-risk and low-risk areas were identi ed using Local Moran's I on crude incidence rate and empirical Bayesian smoothed rate. These ndings indicate that GIS is a feasible approach for identifying spatialtime patterns of particular infectious diseases. The presented approach can be applied to other geographical areas, or other forms of disease. Our results are designed to help policymakers implement more targeted interventions for disease prevention and control. Furthermore, these ndings can be used to hold community training programs in high-risk areas so as to improve the awareness and skills of individuals who are at the higher risk of disease for using personal protective equipment and home spraying of sand ies. Furthermore, identifying high risk areas can help in effective health source allocations such as distribution of free window grids.

Consent for publication
Not applicable.

Availability of data and materials
All data presented in the study are available from the corresponding authors.

Competing interests
The authors declare that they have no competing interests.

Funding
This study received funding from Mashhad University of Medical Sciences (number= 971558). The funder provided the cost of geocoding and cleaning the data. E.M and M.F provided the data. A.R.R geocoded the data. R.S prepared and cleaned the data. B.K, N.F and A.M conducted spatial analyses. B.K and N.F drafted the manuscript. M.M, E.M and D.H revised the text. B.K designed the study and was the research leader. The authors read and approved the nal manuscript.