Zoonotic cutaneous leishmaniasis in northeastern Iran: a GIS-based spatio-temporal multi-criteria decision-making approach

Summary

Zoonotic cutaneous leishmaniasis (ZCL) constitutes a serious public health problem in many parts of the world including Iran. This study was carried out to assess the risk of the disease in an endemic province by developing spatial environmentally based models in yearly intervals. To fill the gap of underestimated true burden of ZCL and short study period, analytical hierarchy process (AHP) and fuzzy AHP decision-making methods were used to determine the ZCL risk zones in a Geographic Information System platform. Generated risk maps showed that high-risk areas were predominantly located at the northern and northeastern parts in each of the three study years. Comparison of the generated risk maps with geocoded ZCL cases at the village level demonstrated that in both methods more than 90%, 70% and 80% of the cases occurred in high and very high risk areas for the years 2010, 2011, and 2012, respectively. Moreover, comparison of the risk categories with spatially averaged normalized difference vegetation index (NDVI) images and a digital elevation model of the study region indicated persistent strong negative relationships between these environmental variables and ZCL risk degrees. These findings identified more susceptible areas of ZCL and will help the monitoring of this zoonosis to be more targeted.

INTRODUCTION

Leishmaniasis, one of the most neglected tropical diseases and a high-priority public health issue, is spread to humans by the bite of infected female sand flies. The World Health Organization (WHO) reported that the public health impact of leishmaniasis has been greatly underestimated for many years, and has classified the disease as an emerging and uncontrolled disease [1]. Leishmaniasis includes a wide variety of diseases which are classified into three main categories including: cutaneous leishmaniasis (CL), visceral leishmaniasis (VL) and muco-cutaneous leishmaniasis (MCL), among which CL is the most common form in Iran. The disease remains a global problem in which about 2 million new cases are believed to occur annually (1·5 million CL and 0·5 million VL cases) [2]. The geographical distribution of leishmaniasis is almost entirely restricted to tropical and subtropical regions predominantly in developing countries including Afghanistan, Algeria, Brazil, Iran, Peru, Saudi Arabia and Syria. CL is one of the most significant infectious diseases in Iran with a yearly average number of 20 000 cases [3]. Both epidemiological entities of CL including zoonotic cutaneous leishmaniasis (ZCL) caused by Leishmania major and anthroponotic cutaneous leishmaniasis (ACL) due to L. tropica are present in Iran; however, about 80% of cases are of the ZCL form [4]. Despite the fact that this common zoonotic disease rarely causes severe morbidity, the long lasting lesions leave unpleasant scars on the face or other exposed skin areas. Therefore, the disease not only causes serious psychological effects for the patients but also a great socioeconomic burden to society.

Despite all conducted efforts to control the disease, it is still endemic and has been reported in different provinces of Iran, especially Golestan province where ZCL prevention and control is highly emphasized by Ministry of Health of Iran. Incidence of ZCL in this endemic province is well correlated with the climate conditions [5]. Variations in the spatial disease pattern are influenced by environmental and landscape factors, which determine the distribution and abundance of vectors and reservoirs [6]. As diagnosed by molecular methods in Golestan province, the vector of the disease is the female sand fly of Phlebotomus papatasi species which is distributed in almost all parts of Iran, the parasitic agent of ZCL is L. major with rodents as reservoirs hosts [7].

There are many epidemiological studies conducted throughout Iran that have concentrated on the zoonotic features of ZCL and the role of reservoir hosts for its transmission, but little attention has been paid to environmental and geographical characteristics affecting the disease spread. Since 1990, Geographic Information System (GIS) has become a powerful tool in epidemiological investigations and disease surveillance providing a potential tool for monitoring ZCL epidemics in tropical countries. Although this tool is well-known in health systems in developed countries, its application is very limited in developing countries such as Iran [8]. Several studies regarding the use of GIS and spatial analysis of CL have been reported from different parts of the world. In this region, Mollalo et al. [9] showed that the geography of the area has played an important role in ZCL distribution. Moreover, the spatial distribution of ZCL due to environmental factors was confined to the northern and north-eastern low-lying regions of the province. In another study conducted in Fars province located in southern Iran, a spatial environmentally based model using a geographic weighted regression model was developed. The results showed that minimum temperature, maximum relative humidity, population density, days of rainfall and wind velocity were the most significant risk factors explaining 0·388 of the associated factors of CL [10]. In central Tunisia, Salah et al. [11] used the scan statistics technique to identify spatial, temporal and spatio-temporal clusters of ZCL to visualize the abnormally high incidence rates. Their results demonstrated that the most likely spatial clusters were located close to a dam. Seid et al. [12] used GIS and statistical analysis to develop a risk map for CL based on environmental factors in Ethiopia. Their results indicated that slope, elevation and annual rainfall were good predictors of CL presence based on the probabilistic and weighted overlay approaches.

Although GIS-based multi-criteria decision-making (MCDM) methods have been occasionally used for zoonosis diseases, according to Clements et al. [13], it provides a great understanding about the uncertainty related to the geographical distribution of diseases. Modelling of a disease based on its cause-and-effect parameters can support public health welfare decision makers in monitoring, controlling and preventing diseases. This study thus concentrated on the identifying high- and low-risk areas of ZCL using the combination of GIS and MCDM analyses. Making proper decisions based on analyses of such information, would facilitate reaching desired results in a shorter time with less costs.

MATERIAL AND METHODS

Study area

This analytical cross-sectional investigation was conducted over a period of 3 years from January 2010 to December 2012 throughout Golestan province, northeastern Iran. As shown in Figure 1, this province is located between latitudes 36° 30′ to 38° 8′ N of equator and longitudes 53° 57′ to 56° 22′ E of the Greenwich meridian. The province with the centre of Gorgan is one of the 31 provinces of Iran, bordering Turkmenistan, occupying an area of 20 893 km² with a total population of about 1 750 000 people. The data concerning administrative boundaries of province and counties, locations of villages as well as population statistics were obtained from the Ministry of Interior of Iran for 2011 under the Shapefile format (the Ministry of Interior of Iran, unpublished data).

Fig. 1. Geographical location of Golestan province and its counties, North-east Iran.

This region is characterized by widely diverse regional and/or seasonal variations, ranging from extremely harsh conditions to very hot, dry summers. Northern regions are located in arid and semi-arid climates with hot and dry summers and cool and rainy winters. The southern regions represent a mountainous climate, and central and southwest parts of the province lie under the influence of the Mediterranean climate (Weather Centre, Hashemabad of Gorgan, 2007, unpublished data).

Data collection and preparation

Various data in different formats (attribution, vector and raster formats) and different scales received from multiple data sources [Centre for Disease Control and Prevention (CDC), meteorological stations, Ministry of Interior]. The data were prepared in a GIS environment so as to have the same coordinate system, spatial extent, and spatial and temporal resolution. Disease data together with various environmental variables which directly/indirectly effect on the ecology of ZCL (mainly vectors and reservoirs) were studied and selected. For the disease cases, suitable spatial and attribution data were collected and manipulated, having village as the spatial unit.

Epidemiological data

This study is based on the passive data of 2893 indigenous ZCL cases from 2010 to 2012. Cases were confirmed by a positive skin test and/or parasitological examination and were officially provided by Golestan CDC. Data containing monthly reports of disease onset, and place of residence at the village level were checked meticulously to prevent any possible duplicate records and were linked with their corresponding geographical location at the village level. Detailed information with regard to demographic characteristics of ZCL cases and identification methods have been presented in the paper by Mollalo et al. [9].

Climate data

In addition to ZCL data, climate characteristics including the minimum, maximum and mean temperature (°C); minimum, maximum and mean relative humidity (%); mean evaporation (mm), accumulative precipitation (mm), and number of freezing days (number of days with minimum temperatures ⩽0° C) and rainy days were retrieved from Golestan province meteorological centres and also neighbouring provinces for a more accurate interpolation process. Taking into account spatial autocorrelations, raster maps were produced using the inverse distance weighting (IDW) method for each of the aforementioned parameters in a GIS framework with a yearly interval at 100 m spatial resolution. Based on the reported data of synoptic stations of the study area, during the study period, the annual average of minimum, maximum, and mean temperature were 12·6 °C, 22·9 °C and 17·8 °C, respectively; the yearly average of minimum, maximum and mean humidity were 53%, 88% and 71%, respectively; yearly total rainfall ranged between 69·8 and 103·2 mm; and the yearly number of rainy and freezing days were about 76–108 days and 10–26 days, respectively.

Remote sensing and topographic data

To reflect vegetation cover, mean annual normalized difference vegetation index (NDVI) values, a dimensionless index, which are calculated from the red (R) and near infrared (NIR) bands as follows were originated from MODIS sensor data (http://modis.gsfc.nasa.gov/):

$${\rm NDVI} = {\rm} \left( {{\rm NIR} -{\rm R}} \right)/\left( {{\rm NIR} + {\rm R}} \right).$$

According to the MODIS satellite data, during the study period, the average range of NDVI was somewhere between −0·15 and 0·80 with a mean value of 0·14. In addition, locations of dams in Golestan province were extracted from Google Earth satellite images. Figure 2 shows generated maps of the climate factors by the IDW method together with the average NDVI satellite image of MODIS sensor, during the study period. Moreover, assuming topography may influence on the ZCL distribution, topographical data providing elevation data at 90 m spatial resolution were derived from a STRM digital elevation model.

Fig. 2. Generated climate maps using inverse distance weighting method as well as MODIS normalized difference vegetation index (NDVI) image in Golestan province, Iran, 2010–2012.

Spatial analyses

A three-stage approach was adopted to investigate ZCL risk degrees in Golestan province. First, by considering the knowledge of local health authorities, Pearson's correlation analysis for each year was conducted to choose the independent variables. Then, linear multivariate regression analyses were developed to determine which independent environmental factors plays more important role in ZCL incidence distribution. Second, both analytical hierarchy process (AHP) and fuzzy AHP methodologies were employed to produce the ZCL risk map for each year, separately. Finally, the obtained results were compared with geo-referenced ZCL data. Furthermore, frequency of ZCL cases in each risk zone was calculated and compared with corresponding environmental factors.

AHP

Each independent affecting environmental factor has a different level of importance and weight in modelling of ZCL risk. To obtain the mapping weight or importance of each ZCL risk factor, the AHP method proposed by Saaty [14] was employed. This technique is a powerful, flexible, and systematic tool for creating a hierarchical model of the spatial decision problem, assisting policy makers to set priorities and make better decisions [15]. AHP decomposes an unstructured and complex problem to a hierarchy so that it is easier to study and weight each criteria and sub-criteria individually.

To designate the importance of each factor, it is necessary to specify the weight of each. One of the advantages of MCDM is that it is based on the pairwise comparison that facilitates the decision-making because it reveals the compatibility/incompatibility of the decisions. In the pairwise comparison weighting method, criteria are compared two by two and their importance is determined towards each other. A rank matrix is created where its inputs are the same determined weights and its outputs are the relative weights related to criteria. For example, in the pairwise matrix the element a _ij  = 3 indicates that the ith element is slightly more important than jth element in ZCL risk, thus the knowledge of users can be applied in these methods. Detailed information in terms of a pairwise matrix has been presented in Tables 1 and 2. In this research, in order to weight the factors using pairwise comparison method, the Expert Choice version 11 software (expertchoice.com) was used. Table 3 presents the assigned relative weights of each input parameter for each year, separately.

Table 1. Saaty's [14] pairwise comparison table with 9 degrees

Table 2. Triangular fuzzy conversion scale

Table 3. Environmental variables affecting zoonotic cutaneous leishmaniasis incidence rate utilized in the models along with their initial weights and positive/negative signs in 2010, 2011 and 2012, Golestan province, Iran

Moreover, one of the most important advantages of this process is the ability to estimate the inconsistency rate of pairwise comparison matrix for more accurate judgments. For this purpose, the consistency ratio (CR) was used. This index reflects the consistency of one’s judgment so that a CR of ≤ 0·1 is considered acceptable. Conversely, any higher value at any level is considered an inconsistent decision indicating that the judgement requires revision.

Fuzzy AHP

One limitation of AHP is that this method does not consider the uncertainty and imprecision associated with the mapping and nature of the variables [16]. To reduce the accuracy of this major drawback, an alternative to the AHP method, designated FAHP, was employed. This method uses linguistic expression for modelling uncertainty in the decision-making process [17]. Each effective factor in the prediction model was fuzzified using triangular fuzzy numbers according to Chang's extent analysis method [18]. Similar to the stages of the AHP method the risk maps in all the 3 years were provided by running MATLAB codes with the exception of fuzzy operators and rules.

After obtaining the weights of all contributing factors by MCDM analyses, pixel values of corresponding factors, multiplied by their respective weights, and the final value of each pixel (i.e. produced risk maps) calculated by weighted arithmetic mean formula in a GIS environment for further spatial analyses. Based on the distribution of data (pixel values of output raster map) the natural-break classification method which identifies real classes was used to classify the ZCL risk map. This method finds points that minimize within-class variance and maximizes variances between classes. The obtained risk maps were then classified into five classes of risk categories ranging from ‘very low’ to ‘very high’ risk areas for prioritizing control interventions. The range of values of the ZCL risk map, produced by both AHP and FAHP are presented in Table 7, for each year separately. The accuracy of the produced risk maps was assessed by comparing them with the spatial distribution of ZCL in the study area. Moreover, the spatially averaged environmental factors for each risk category were calculated and evaluated.

RESULTS

Table 3 summarizes the independent explanatory variables together with their weights resulting from both Pearson's correlation and multivariate regression analyses of the relationship between ZCL incidence and different environmental variables. Then, the weights were used as initial weights in AHP and FAHP methods. In total, for all of the 3 years of the study, topography and vegetation cover had a major inverse relationship in ZCL distribution, while temperature variables had a direct relationship.

In this study, CRs of 0·01 and 0 were obtained which are below the threshold of 0·10 suggesting a reasonable level of consistency in the pairwise comparisons or acceptable judgement in identifying ZCL risk zones.

After conducting both modelling methods with the information of the pairwise comparison matrices presented in Tables 4 and 5, the final weights of each criteria were calculated (Table 6) and the layers were overlaid based on their obtained weights within the GIS environment. As the generated risk maps, which were classified based on the natural-break method for both methods show (Fig. 3), high-risk areas were mainly located at the northern and northeastern parts of the province while low-risk areas were located at the southern parts. Comparison of the generated risk maps with geocoded ZCL cases at the village level demonstrated that in both methods more than 90%, 70% and 80% of the cases occurred in high and very high risk areas for the years 2010, 2011, and 2012, respectively (Fig. 4). The results indicates the capability of both models to predict susceptible ZCL areas with an accuracy exceeding 70%.

Fig. 3. AHP (top row) and FAHP (bottom row) derived zoonotic cutaneous leishmaniasis risk zones in Golestan province, Iran, 2010–2012.

Fig. 4. Frequency of zoonotic cutaneous leishmaniasis (ZCL) occurrence in different level of risks using (a) AHP and (b) FAHP methods in Golestan province, Iran, 2010–2012.

Table 4. Criteria pairwise comparison matrix, Golestan province, Iran, in 2010, 2011 and 2012

Alt, Altitude from the sea level (m); Frz, number of freezing days; Hum, humidity (%); NDVI, normalized difference vegetation index; Max temp, maximum temperature (°C); Evap, evaporation (mm); Max wind, maximum direction of wind velocity; Temp, temperature (°C); Prc, precipitation (mm)

Pairwise comparison matrixes are symmetric, thus real numbers like 0·33 indicates that Evap is somewhat more important than NDVI (the weight of Evap is three times more than NDVI).

Table 5. Fuzzified criteria pairwise comparison matrix, Golestan province, Iran, in 2010, 2011 and 2012

Alt, Altitude from the sea level (m); Frz, number of freezing days; Hum, humidity (%); NDVI, normalized difference vegetation index; Max temp, maximum temperature (°C); Evap, evaporation (mm); Max wind, maximum direction of wind velocity; Temp, temperature (°C); Prc, precipitation (mm).

CR_m and CR_g are consistency ratios of fuzzy pairwise comparison matrix which were defuzzified to crisp numbers.

Table 6. Final weights calculated from AHP and FAHP approaches for each factor in 2010, 2011 and 2012, Golestan province, Iran

AHP, Analytical hierarchy process; FAHP, fuzzy analytical hierarchy process.

Table 7. The range values of zoonotic cutaneous leishmaniasis risk map produced by AHP and FAHP, 2010–2012

AHP, Analytical hierarchy process; FAHP, fuzzy analytical hierarchy process; VL, very low degree of risk; L, low degree of risk; M, moderate degree of risk; H, high degree of risk; VH, very high degree of risk.

Moreover, using several spatial analyses, including zonal statistics, each category of risk map was compared in relation to their corresponding spatially averaged environmental variables. Based on the results of AHP and FAHP presented in Table 8, comparison of the risk categories with spatially averaged NDVI images for each of the 3 years shows that the majority of cases occurred at low and very low vegetation areas. As depicted in Figure 5, it can be clearly seen that there is a persistent strong negative relationship between the vegetation cover and level of ZCL risk in Golestan province, implying the relationship remains somewhat similar. Similarly, in the 3 years of study, the effect of altitude on ZCL occurrence is apparent and persistent as well. High-risk areas were situated at altitudes between 500 m above mean sea level suggesting ZCL occurrence tended to be more prevalent in the plains or at relatively low altitudes, while low-risk areas were located at high altitudes (Fig. 6).

Fig. 5. Level of zoonotic cutaneous leishmaniasis risks concerning mean of normalized difference vegetation index (NDVI) using (a) AHP and (b) FAHP methods, Golestan province, Iran, 2010–2012.

Fig. 6. Level of zoonotic cutaneous leishmaniasis risks concerning mean of altitudes using (a) analytical hierarchy process (AHP) and (b) fuzzy AHP methods, Golestan province, Iran, 2010–2012.

Table 8. ZCL level of risks expressed in terms of frequency of cases and environmental variables extracted from AHP and FAHP risk models, Golestan province, Iran (2010-2012)

ZCL, Zoonotic cutaneous leishmaniasis; NDVI, normalized difference vegetation index; AHP, analytical hierarchy process; FAHP, fuzzy analytical hierarchy process.

DISCUSSION

ZCL risk maps play a key role in public health and epidemiology of zoonoses in that they highlight areas which are more susceptible and more suitable for breeding and maintenance of sand flies and reservoirs with a high incidence rate. Thus, visualization of high- and low-risk areas can provide valuable information for public health decision makers in geographical management of ZCL occurrence and give direction regarding where their control efforts, such as prioritizing proper allocation of the budget, personnel and equipment, should be targeted.

This study confirmed the capabilities of decision-making methods in monitoring and prediction of ZCL occurrence. Other studies have been conducted using these decision-making techniques. For example, one spatial study in Thailand developed AHP and FAHP models for hand, foot and mouth disease. Their results showed that FAHP performed more accurately than AHP in defining high-risk areas [19]. However, in the present study, both methods successfully predict high and very high risk areas in all 3 years of study, the models were also in agreement with each other in defining high risk areas with no significant differences.

It should be noted that MCDA provides a great opportunity for using the knowledge and experience of the users or experts in weighting factors by using the abilities of the pairwise matrix, this way it is possible to test various scenarios that meet the requirement of the CR. In reality, it is not true to assume the whole study area as spatially homogeneous because the criteria vary across space [20]. Regression analysis was employed based on non-spatial to provide global, simple and quick initial weights for running AHP in Expert Choice software as well as to investigate the independence of parameters. In other words, because of insufficient knowledge of environmental factors, regression analysis was used to determine significant factors and their initial weights through a significance test of R ², while it was possible to gain the weights based on the experience of the expertise or even trial-and-error method, as well.

The public health surveillance system of Iran has been well-founded, especially in rural areas, with more than 95% coverage of the entire country. One of the basic duties of the system is to provide primary healthcare and to register health records including ZCL cases [21]. However, in reality, most ZCL cases are not observed, some are observed but not recognized; some are recognized but not reported [22]. According to Alvar and colleagues [23], the estimated degree of under-reporting of cutaneous leishmaniasis cases in Iran is somewhere between 2·8- and 4·8-fold. On the other hand, pure statistical models require accurate and reliable data preferably for a long period of time in order to reach a solid conclusion. Thus, in such poor conditions of observed data, AHP and FAHP models are superior to pure statistical models to fill the gap of reliable and historical data.

The approximate location of identified high-risk zones, based on the used factors in the present study, both support and extend the findings of previous work of cluster detection analysis in the same study area. The findings of Mollalo et al. [9, 24], using spatial scan statistics cluster detection technique in this endemic area, illustrated that the most likely spatial clusters were located in northern and northeastern parts of the study area, with arid and semi-arid climates and low vegetation cover, supporting the view that this areas contains potential high-risk populations and warrants closer consideration. However, the methods used in this study to define high-risk areas are more robust than cluster detection methods, which were merely based on case incidence rates due to the fact that results of clustering methods cannot be extrapolated to other areas, while risk factors can be used to identify high-risk zones in other areas where risk factor data are available [25].

Regardless of investigating the stability of the results, the most important reason for conducting analyses for each year separately was based on the descriptive statistics of a previous study in this region by Mollalo et al. [9] who observed that the frequency of ZCL cases in 2010 (1660 cases) significantly decreased in the next years (660 cases in 2011); thus we supposed that environmental changes in the area may be responsible for such a huge difference. It can be a reason indicating why the set of environmental variables differed from one year to another. It is obvious that some phenomena like earthquakes, flood, or even climate change can change the behaviour and pattern of infected vectors and consequently the frequency of human ZCL cases.

Results of the risk models were further compared with discrete correlation analysis between environmental variables and ZCL incidence rate for the same study region [5], where dynamic monthly significant associations were observed between environmental variables (including vegetation cover, topography and climate factors) and ZCL incidence rate in Golestan province. Both studies showed strong negative influence of altitude and NDVI factors, as proxies of environmental changes, on the ZCL incidence rate providing an excellent niche for ZCL transmission. While temperature and relative humidity variables in discrete correlation analyses showed strong significant association with ZCL incidence rate, the current study signified poor relationships between these variables. This might be due to the fact that in the previous study the associations between each variable and ZCL occurrence were analysed individually, regardless of influence or dependency between parameters. However, the advantage of such modelling techniques is that not only is it based on independent constituting variables but also the interaction between the variables can be reflected in the results. It should be noted that the results of above models might vary by different geographical areas due to different seasonal patterns in different ecological zones of the province, thus the results of this study are not plausible for other endemic areas. For instance, visual comparison between location of dams and produced risk maps of the study area shows that almost all of the dams are located in low and very low risk areas, which is in contrast with the study of Salah et al. [11] who observed location of spatial hotspots close to dams in Tunisia.

It should acknowledged that the current study has shortcomings from two perspectives. The first weakness is short data length (3 years) which might not lead to robust and reliable results and the second limitation is related to absence of biological conditions and socioeconomic status in the models which are important determinants of ZCL risk. Therefore, future research should consider other neglected influencing factors and/or even culture and life-style of the population at risk to better describe the epidemiology of the disease. Despite existing limitations encountered in the current study, these findings identified more susceptible areas and will help to make the control and monitoring of ZCL more targeted in Golestan province.

CONCLUSION

By identifying the independent environmental factors associated with the ZCL incidence rate, the present work emphasized the spatial characteristics of this common zoonosis in an endemic focus of Iran. The developed risk maps provide an effective visual tool for public health policy makers with regard to where the control programmes must be targeted, assisting more optimal allocation of budget and health facilities under future environmental changes. Integration of popular decision-making methods (i.e. AHP and FAHP modelling techniques) and powerful GIS can help the results to be more precise, knowledge-based, and cost-effective in defining high-risk areas of ZCL.