A Comparative Assessment of Epidemiologically Different Cutaneous Leishmaniasis Outbreaks in Madrid, Spain and Tolima, Colombia: An Estimation of the Reproduction Number via a Mathematical Model

Leishmaniasis is a neglected tropical disease caused by the Leishmania parasite and transmitted by the Phlebotominae subfamily of sandflies, which infects humans and other mammals. Clinical manifestations of the disease include cutaneous leishmaniasis (CL), mucocutaneous leishmaniasis (MCL) and visceral leishmaniasis (VL) with a majority (more than three-quarters) of worldwide cases being CL. There are a number of risk factors for CL, such as the presence of multiple reservoirs, the movement of individuals, inequality, and social determinants of health. However, studies related to the role of these factors in the dynamics of CL have been limited. In this work, we (i) develop and analyze a vector-borne epidemic model to study the dynamics of CL in two ecologically distinct CL-affected regions—Madrid, Spain and Tolima, Colombia; (ii) derived three different methods for the estimation of model parameters by reducing the dimension of the systems; (iii) estimated reproduction numbers for the 2010 outbreak in Madrid and the 2016 outbreak in Tolima; and (iv) compared the transmission potential of the two economically-different regions and provided different epidemiological metrics that can be derived (and used for evaluating an outbreak), once R0 is known and additional data are available. On average, Spain has reported only a few hundred CL cases annually, but in the course of the outbreak during 2009–2012, a much higher number of cases than expected were reported and that too in the single city of Madrid. Cases in humans were accompanied by sharp increase in infections among domestic dogs, the natural reservoir of CL. On the other hand, CL has reemerged in Colombia primarily during the last decade, because of the frequent movement of military personnel to domestic regions from forested areas, where they have increased exposure to vectors. In 2016, Tolima saw an unexpectedly high number of cases leading to two successive outbreaks. On comparing, we estimated reproduction number of the Madrid outbreak to be 3.1 (with range of 2.8–3.9), which was much higher than reproduction number estimates of the Tolima first outbreak 1.2 (with range of 1.1–1.3), and the estimate for the second outbreak in Tolima of 1.019 (with range of 1.018–1.021). This suggests that the epidemic outbreak in Madrid was much more severe than the Tolima outbreak, even though Madrid was economically better-off compared to Tolima. It indicates a potential relationship between urban development and increasing health disparities.

most vulnerable population, due to the continuous deployment of troops to forested areas of high endemicity and high density of the insect vector [15].
Health disparities: A central aspect of disparities is to identify and study differences in health status between groups, which negatively impact less advantaged groups. These differences could be because of socioeconomic status, gender and ethnicity disparities, and accessibility to health care and interventions. For example, despite the United States' economic dominance and status as one of the most developed countries, an estimated 12 million Americans living in poverty suffer from at least one neglected tropical disease (NTD) [16]. While there are immense challenges to systematically investigating the potential impact of health disparities on an outbreak, overall potential of transmission and size of an infection for distinct populations can be estimated via metrics such as reproduction number, inoculation rate, epidemic size, vectorial capacity, and so forth. In this study, we modeled leishmaniasis outbreaks in the regions of Madrid, Spain and in Tolima, Colombia, as two 'distinct' populations to study differences in transmission potential. Spain's healthcare system is regularly rated among the world's best (with~90% of patients accessing public healthcare and~20% accessing some part of private healthcare) and spends about 10% of its gross domestic product (GDP) on healthcare. State healthcare guarantees universal coverage, although one may have to travel far to find, or wait a significant time to access, a public healthcare facility. On the other hand, Colombia ranks 22nd on the WHO's list of the best healthcare systems, with private healthcare establishments accounting for around 57% of establishments (thus, relatively rapid access to healthcare). Total expenditures on health constitute around 7% of Colombia's GDP; however, urban and rural areas have significant differences in access to health care (see Table in the Supplementary Material, Section S.1, for details on potential health disparity between Spain and Colombia). Here, we do not aim to identify (or study) specific factors for health disparity leading to a leishmaniasis epidemic. However, a general discussion of the potential impact of health disparity on the disease outbreak is provided. Significant social and environmental data are needed to truly capture the differences and study the role of health disparity in the transmission dynamics of CL.
Mathematical modeling study of leishmaniasis and reproduction number: Initial mathematical models of CL transmission dynamics were developed by Dye et al. [17,18], in which they analyzed simple discrete-time epidemic models to study the mechanism behind observed inter-epidemic periods and the intensity of infection in dog reservoirs. Other studies included models with heterogeneous biting among age-structured dog populations and used serological data for the dog population in Gozo, Malta to estimate the basic reproduction number [19], a quantity that measures the intensity of an outbreak. Formally, the basic reproduction number (R 0 ) of an infection can be interpreted as the average number of new cases generated by a typical infectious individual over the course of its infectious period, in an otherwise uninfected population. It is a key parameter, the value of which characterizes the transmission potential of an epidemic and hence, is often used to inform the potential effectiveness of intervention strategies. There are various ways to estimate R 0 [20,21]. Here, we derive three novel methods to estimate R 0 for a CL outbreak via a mathematical model with region-dependent features. The methods are tested using data from two ecologically distinct regions-Madrid (a city in a 'developed' country, Spain) and Tolima (a city in a 'developing' country, Colombia)-as a case study. In the literature, models have especially been used to estimate R 0 [18] (see Table 1 for review on R 0 estimates); however, such studies have used data primarily from different unrelated studies to collect point estimates of model parameters instead of applying a rigorous parameter estimation procedure. In the present study, we developed mathematical procedures for the estimation of model parameters via fitting the model to temporal incidence data, using three different techniques. Research focus of the study: This study attempts to understand three major CL outbreaks, the first one in Madrid, Spain from 2009 to 2012 and the second and third outbreaks in Tolima, Colombia-both occurred in 2016. The outbreak in Madrid was started mainly by dogs, which are reservoir hosts of the disease, and the outbreaks in Tolima initiated because of the movement of soldiers, particularly those coming from the jungle to urban areas after staying in forestlands for long periods of time. In this work, we study the dynamics of CL in Spain and Colombia using a simple vector-borne disease model, while incorporating local characteristics and data on the disease. We used the country-specific model to estimate the transmission potential of each of three outbreaks via three novel parameter estimation procedures. These two regions were selected because of their distinct characteristics related to the disease and to highlight the comparison of the transmission potential between the ecologically and economically different regions. The estimated model parameters were then used to estimate the local reproduction number for each outbreak and each region. The difference in the basic reproduction number between regions could shed light on potential differences in health inequality, population immunity and transmissibility of leishmaniasis. This information is also important for designing effective control policies.

Materials and Methods
We used outbreak specific-models for the two regions-Madrid and Tolima-and estimated their parameters using three different estimation methods (methodologies shown in Figure 1; details of the methods are provided in later sections).

Data Sources for Madrid, Spain
The data corresponding to the 2009-2012 CL outbreak in communities of Madrid, Spain were collected ( Figure 2) for the study by Mar Noguerol Álvarez et al. [30]. These reported data consisted of a monthly collection of new cases over the four-year period. The epidemic started in July 2009 (7th month of the year), ended in March 2012 and resulted in 156 total cases. The patients were primarily reported from four communities, namely, Fuenlabrada, Leganés, Getafe, and Humanes de Madrid. The mean incidence was found to be 14 cases per 100,000 inhabitants during this period. The incidence rate during the 2009 to 2012 epidemic was much higher than the corresponding rate during the 2000 to 2009 period from the same regions (between 1 and 6 new cases per year or incidence rate of less than 1 per 100,000 inhabitants). More than 60% of reported cases were male with a mean age of 46 years. L. infantum was identified as the causative pathogen species. It is important to note that around one quarter of cases had contact with dogs in one or more places in the domestic or peridomestic  The mean incidence was found to be 14 cases per 100,000 inhabitants during this period. The incidence rate during the 2009 to 2012 epidemic was much higher than the corresponding rate during the 2000 to 2009 period from the same regions (between 1 and 6 new cases per year or incidence rate of less than 1 per 100,000 inhabitants). More than 60% of reported cases were male with a mean age of 46 years. L. infantum was identified as the causative pathogen species. It is important to note that around one quarter of cases had contact with dogs in one or more places in the domestic or peridomestic The mean incidence was found to be 14 cases per 100,000 inhabitants during this period. The incidence rate during the 2009 to 2012 epidemic was much higher than the corresponding rate during the 2000 to 2009 period from the same regions (between 1 and 6 new cases per year or incidence rate of less than 1 per 100,000 inhabitants). More than 60% of reported cases were male with a mean age of 46 years. L. infantum was identified as the causative pathogen species. It is important to note that around one quarter of cases had contact with dogs in one or more places in the domestic or peridomestic environment. The present study uses these 2009-2012 epidemic data from Madrid and estimates the basic reproduction number for the outbreak using a mathematical model.

Data Sources for Tolima, Colombia
The CL epidemic was observed in Tolima, Colombia in 2016 ( Figure 3). The epidemic data representing new cases per week were obtained from the National Institute of Health's Weekly Bulletin [10]. The outbreak started during the first week of January of 2016 (epidemiological week 1) and ended during the first week of January of 2017 (epidemiological week 53). There were two consecutive outbreaks resulting in 3223 total reported cases. The second outbreak might have been a result of heavy movement of troops back to the city from forest areas as a result of peace deal signed between government and Revolutionary Armed Forces of Colombia (FARC) rebels. The first outbreak occurred between week 1 and week 35 and the second outbreak occurred between week 36 and week 53. The highest number of reported cases (~149) occurred on 15 May (~week 20) and 15 November (~week 46). In this study, we use these two outbreaks in Tolima and estimate the corresponding basic reproduction numbers using a mathematical model. The CL epidemic was observed in Tolima, Colombia in 2016 ( Figure 3). The epidemic data representing new cases per week were obtained from the National Institute of Health's Weekly Bulletin [10]. The outbreak started during the first week of January of 2016 (epidemiological week 1) and ended during the first week of January of 2017 (epidemiological week 53). There were two consecutive outbreaks resulting in 3223 total reported cases. The second outbreak might have been a result of heavy movement of troops back to the city from forest areas as a result of peace deal signed between government and Revolutionary Armed Forces of Colombia (FARC) rebels. The first outbreak occurred between week 1 and week 35 and the second outbreak occurred between week 36 and week 53. The highest number of reported cases (~149) occurred on 15 May (~week 20) and 15 November (~week 46). In this study, we use these two outbreaks in Tolima and estimate the corresponding basic reproduction numbers using a mathematical model.

Mathematical Model for CL Epidemic in Madrid, Spain
Since dogs are the major reservoir of the disease in Spain, we considered a transmission dynamics model consisting of human (represented by subscript 'h'), dogs (represented by subscript 'A'), and sandfly vector (represented by subscript 'v') populations [16]. The susceptible, infectious and recovered subcategories for the three populations are represented by X, Y, and Z, respectively. The flowchart of the epidemic model for the Madrid outbreak is illustrated in Figure 4. The model equations are given in supplementary material (Section S.2) and state variables and parameters are explained in Table 2.

Mathematical Model for CL Epidemic in Madrid, Spain
Since dogs are the major reservoir of the disease in Spain, we considered a transmission dynamics model consisting of human (represented by subscript 'h'), dogs (represented by subscript 'A'), and sandfly vector (represented by subscript 'v') populations [16]. The susceptible, infectious and recovered subcategories for the three populations are represented by X, Y, and Z, respectively. The flowchart of the epidemic model for the Madrid outbreak is illustrated in Figure 4. The model equations are given in Supplementary Material (Section S.2) and state variables and parameters are explained in Table 2.     Here, N h , N v and N A represent the total population of humans, vectors and dogs, respectively. We assume these populations to be constant by taking equal natural birth and death rates (where, µ h and µ A represented both birth and death per capita rates for the three populations).
respectively. Similarly, the probability that a vector bites a dog is N A respectively. The corresponding model's equations are included in the Supplementary Material (Section S.2.1).

Mathematical Model for CL Epidemic in Tolima, Colombia
Since the movement of military personal from the forest to civilian regions resulted in an unprecedented higher number of cases in Colombia, we considered the transmission dynamical model consisting of civilian (represented by subscript 'c'), military (represented by subscript 'm') and sandfly vector (represented by subscript 'v') populations interacting with each other [16]. The susceptible, infectious and recovered subcategories are represented by S, i, and R, respectively. The model framework that is used for the two Tolima, Colombia outbreaks is illustrated in the flow chart in Figure 5. The model state variables and parameters are explained in Table 3.

Mathematical Model for CL Epidemic in Tolima, Colombia
Since the movement of military personal from the forest to civilian regions resulted in an unprecedented higher number of cases in Colombia, we considered the transmission dynamical model consisting of civilian (represented by subscript 'c'), military (represented by subscript 'm') and sandfly vector (represented by subscript 'v') populations interacting with each other [16]. The susceptible, infectious and recovered subcategories are represented by S, i, and R, respectively. The model framework that is used for the two Tolima, Colombia outbreaks is illustrated in the flow chart in Figure 5. The model state variables and parameters are explained in Table 3.  The civilian or military susceptible population moves to the class of infected population by the bite of an infected vector biting at the rate b. Similarly, the susceptible vector population moves to the infectious class by the bite of a female vector to an infected civilian or military. For each population, the  The civilian or military susceptible population moves to the class of infected population by the bite of an infected vector biting at the rate b. Similarly, the susceptible vector population moves to the infectious class by the bite of a female vector to an infected civilian or military. For each population, the total recruitment rate is Λ c , Λ m and Λ v and the per capita death rates are µ c , µ m and µ v . The recovery rate for the civilian population is γ c . It is assumed that proportion q of the traveling military individuals move to the town are infected and proportion (1-q) are assumed to be susceptible. The corresponding model's differential equations are included in the Supplementary Material (Section S.2.2).
We assumed two different new incoming rate parameters each for civilian (Λ c ) and military (Λ m ) populations. The incoming rate for military incorporate two constant rates: the net recruitment rate in military from civilian population and the net movement rate from other surrounding areas to the modeled region of Tolima. The recruitment into military occurs at national (country) level. In urban areas of Tolima, Colombia, military populations are stationed in battalion camps, which are typically on the outskirts of the urban areas. These military camps have all relevant facilities but only limited interactions with civilians occur. This is to ensure safety of military individuals (as they fight with local guerillas, who are often friendly with civilians) and to maintain secrecy of military operations.

Analysis
The well-known Ross-Macdonald model, developed in 1911, formally initiated the field of modeling of complex transmission cycles of vector borne diseases [31] and provided a theoretical The key quantities such as the basic reproduction number, vectorial capacity and inoculation rates derived from the analysis of vector-borne models become the central to the quantification of transmission. Once the model parameters are estimated, these quantities can be easily computed. Here, we develop the enhanced expression and estimate of some of these quantities using region-specific characteristics and data.

Reproduction Numbers
The type-reproduction number (R T ) for a specific host type is interpreted as the average number of secondary cases of that type produced by the primary cases of the same host type during the entire course of infection. It takes into account not only the secondary cases directly transmitted from the specific host but also the cases indirectly transmitted by way of other types, who were infected from the primary cases of the specific host with no intermediate cases of the target host. It is a useful measure when a particular single host type is targeted in the disease control effort in a community with various types of host [32]. R T can be seen as an extension of R 0 in a sense that the threshold condition of the total population growth can be formulated by the reproduction process of the target type only.
The next generation matrix NGM using Figure 6 can be computed as and hence, the basic reproduction number is where K ij represents average number of new infections among the susceptible of type i, generated by an infected of type j. Note, for vector borne diseases, infected host of type 1 cannot directly infect susceptible host of type 3 and vice versa, and type 2 represents the vector. The type reproduction number for Type One is Trop. Med. Infect. Dis. 2018, 2, 9 of 23 total recruitment rate is , and and the per capita death rates are , and . The recovery rate for the civilian population is . It is assumed that proportion q of the traveling military individuals move to the town are infected and proportion (1-q) are assumed to be susceptible. The corresponding model's differential equations are included in the supplementary material (Section S.2.2).
We assumed two different new incoming rate parameters each for civilian ( ) and military ( ) populations. The incoming rate for military incorporate two constant rates: the net recruitment rate in military from civilian population and the net movement rate from other surrounding areas to the modeled region of Tolima. The recruitment into military occurs at national (country) level. In urban areas of Tolima, Colombia, military populations are stationed in battalion camps, which are typically on the outskirts of the urban areas. These military camps have all relevant facilities but only limited interactions with civilians occur. This is to ensure safety of military individuals (as they fight with local guerillas, who are often friendly with civilians) and to maintain secrecy of military operations.

Analysis
The well-known Ross-Macdonald model, developed in 1911, formally initiated the field of modeling of complex transmission cycles of vector borne diseases [31] and provided a theoretical support for understanding the dynamics of those infections. The key quantities such as the basic reproduction number, vectorial capacity and inoculation rates derived from the analysis of vector-borne models become the central to the quantification of transmission. Once the model parameters are estimated, these quantities can be easily computed. Here, we develop the enhanced expression and estimate of some of these quantities using region-specific characteristics and data.

Reproduction Numbers
The type-reproduction number (RT) for a specific host type is interpreted as the average number of secondary cases of that type produced by the primary cases of the same host type during the entire course of infection. It takes into account not only the secondary cases directly transmitted from the specific host but also the cases indirectly transmitted by way of other types, who were infected from the primary cases of the specific host with no intermediate cases of the target host. It is a useful measure when a particular single host type is targeted in the disease control effort in a community with various types of host [32]. RT can be seen as an extension of R0 in a sense that the threshold condition of the total population growth can be formulated by the reproduction process of the target type only.  If K 32 K 23 > 1 then the series fails to converge and Type Three hosts, N H2 , is a reservoir of infection. The type reproduction numbers of Type Two and Type Three, respectively, are R T2 = K 21 K 12 + K 32 K 23 and R T1 = K 23 K 32 The next-generation matrix (NGM), introduced by Diekmann et al. [33], provides a procedure to derive the basic reproduction number, R 0 . This matrix (often denoted by K = [k ij ]) gives the average number of new infections among the susceptible individuals of type i, generated by an infected individual of type j and R 0 is identified as its dominant eigenvalue (that is, R 0 = ρ(K)). In some special models, K = F × V −1 , where F is the new generation matrix and V represents the transition matrix [34]. However, this is not true for vector-borne models with multiple hosts, as is the case in this study. Note, R 0 < 1 i f f R Ti for all host type i and if R Ti > 1 then host type i is a reservoir of infection.
For the K = [k ij ], one identifies the set of targeted entries S, that is, the set of entries in K that are subject to change in control. The target matrix K S is identified as [K S ] ij = k ij if (i,j) ∈ S, and zero otherwise. The target reproduction number, R S is defined as R S = ρ(K S ·(I − K + K S ) −1 ) provided that ρ(K − K S ) < 1, where I is the identity matrix [35]. The last condition can be referred to as the condition for controllability, since if the spectral radius is greater than 1 then the disease cannot be eliminated by targeting only S (in such case, R S is not defined [36]). The controlled NGM, K c , is formulated by replacing the entry k ij in K by k ij /R S whenever (i,j) ∈ S.
Typically, in the case of the simple one-host vector-borne epidemic model, the computed basic reproduction number is given by R 2 0 = R h 0 R v 0 , where R h 0 represents average number of human cases generated by one vector and R v 0 is the average number of vector cases generated by one host. Therefore, the basic reproduction number R 0 gives the average number of secondary infectious hosts (or vectors) produced by one primary infectious hosts (or vector) introduced in completely susceptible populations of hosts and vectors. The effective reproduction number R e f f (t) can be defined as product of partial effective reproduction numbers as where S h (t), N h (t), S v (t), and N v (t) represents number of susceptible hosts, total size of host population, density of susceptible vectors, total density of vector population, respectively. At the time of the beginning of epidemics, R e f f = R 0 because all hosts and vectors in their respective populations are susceptible. Moreover, for large time, the epidemics reaches a steady state, which occurs due to R e f f = 1. Also, R eff formula could be used to numerically see the difference in Madrid and Tolima outbreaks over time.

Mathematical Computations for Spain Model
The basic reproductive number for the model (S1)-(S7) has the form We transform the complex model (S1)-(S7) and we obtain the following effective simple epidemic model where γ e f f ,h = γ h + µ h and The basic reproduction number for this model is

Mathematical Computations for Colombia Model
The basic reproduction number for the model (S25)-(S31) has the form We transform the complex model (S25)-(S31) into an effective simple SIR model and obtain and for the basic reproduction number we obtain Note, Equations (3)-(5) and (9)-(11) are similar.

Parameter Estimation Procedure
Three different methods were used to estimate model parameters and hence R 0 . Supplementary Material (Section S.3) contains the detailed explanation of the parameter estimation procedures and the corresponding codes. (5) can be written only in terms of single variable Z h (t). Expanding right hand side of this new equation using its Taylor series and approximating only up to quadratic terms. Finally, integrating it to obtain the solution of Z h (t).

Method 0 (Cumulative Incidence Technique) Dividing Equation (3) by (5) and then integrating we can compute expression of X h (t) in terms of
That is, the SIR model (3)-(5) can be solved approximately (see detailed computations in Section 2.3.2 of [37]) as Similarly, the cumulative incidence formula corresponding to the Equations (9)-(11)can be derived.

Method 1 (Incidence Technique)
Taking the temporal derivative of (13) we obtain the theoretical incidence curve given by

Method 2 (Incidence Technique)
Another method to obtain the theoretical incidence is to use (13), discretize time and take the difference between the two successive times, namely This equation can be simplified as

Epidemiological Evaluation Metrics Using Components of R 0
Once R 0 is estimated, we can derive other epidemiological metrics for evaluation of an outbreak. However, additional data may be also needed to compute some of these metrics. We do not use these other metrics for comparison of the outbreaks in the two regions considered in this study.
Vectorial Capacity (VC) describes the potential for a vector population to transmit a parasite and can be interpreted as the total number of potentially infectious bites that would eventually arise from all the vectors biting a single perfectly infectious (i.e., all vector bites result in infection) human on a single day. For the Madrid model and under the assumption of presence of only human host, where m is the number of vectors per hosts (i.e., N v /N h ), n is the extrinsic incubation period and p = e −µ v is the probability of a mosquito surviving a day. In fact, can be written as a function of R 0 [38].
The Entomological Inoculation Rate (EIR) is in general defined as the number of infectious bites per person per unit time and can be computed by multiplying the human biting rate with the fraction of infectious vector [39]. In other words, Often, defining the risk associated with a host or with a region is needed. There are two types of risk indexes commonly used for the vector borne diseases-the Transmission Risk Index and the Vulnerability Risk Index. The transmission risk index for host i, TR i , is defined as the probability P i that a host i gets infected, multiplied by the secondary cases it generate. Thus, it is given by TR i = P i R hi 0 , where P i = N i ∑ N k and summation is over only host populations. This index indicates the risk of the host i becoming the main transmitter of the infection at the beginning of an epidemic. On the other hand, we define the Vulnerability Risk index for host j, VR j , as the secondary infections in host j when another host of type i ( =j) becomes infected first at the beginning of an epidemic. Thus, VR j = R hj 0 ∑ i =j P i R vi 0 .

Results
Note, the next generation matrix NGM for the model in Figure 6 is The components of the K matrix for our two models are given by:

Estimates of R 0 for Madrid-Spain
We use equation (13) to fit the observed curve of accumulated cases to estimate the model parameters (that is, using Method 0; Figure 7) hence, obtaining and we use equation (17) to fit the observed curve of incidence obtaining the following estimated values of the parameters (see Figure 8).
We use equation (19) to fit the observed curve of incidence obtaining the following estimated values of the parameters (see Figure 9)       The details of the estimation procedure and confidence interval are given in the Supplementary Material, Section S.3.1. The results suggest that around 626 individuals (with range of 606-711) were atrisk of CL in Madrid before 2009-2012 outbreak (see Table 4). The estimated reproduction number of the Madrid CL outbreak was found to be 3.1 with range from 2.8 to 3.9 (see Table 4).  Table 4 summarizes the estimates obtained in this case. The details of the estimation procedure and confidence interval are given in the Supplementary Material, Section S.3.1. The results suggest that around 626 individuals (with range of 606-711) were at-risk of CL in Madrid before 2009-2012 outbreak (see Table 4). The estimated reproduction number of the Madrid CL outbreak was found to be 3.1 with range from 2.8 to 3.9 (see Table 4).

Tolima First Outbreak
In this case, we use similar methods to estimate the parameters of the respective model (see Supplementary Material Section S.3.2) obtaining the following results for the two outbreaks. The parameter estimation is carried out using three methods (the fitting are shown in Figures 10-12).     The details of the estimation procedure and confidence interval are given in the Supplementary Material, Section S.3.2. The results suggest that around 6656 individuals (with range of 5288-7301) were at-risk of CL in Tolima-Colombia before 2016 outbreak (see Table 5). The estimated reproduction number of the Tolima CL outbreak was found to be 1.2 (see Table 5). In this case also the parameter estimation is carried out using the same three methods (the fitting are shown in Figures 13-15). The details of the estimation procedure and confidence interval are given in the Supplementary Material, Section S.3.2. The results suggest that around 6656 individuals (with range of 5288-7301) were at-risk of CL in Tolima-Colombia before 2016 outbreak (see Table 5). The estimated reproduction number of the Tolima CL outbreak was found to be 1.2 (see Table 5).

Tolima Second Outbreak
In this case also the parameter estimation is carried out using the same three methods (the fitting are shown in Figures 13-15).      The details of the estimation procedure and confidence interval are given in the Supplementary Material, Section S.3.2. The results suggest that the estimated reproduction number of the Tolima CL second outbreak was found to be 1.02 (see Table 6). The details of the estimation procedure and confidence interval are given in the Supplementary Material, Section S.3.2. The results suggest that the estimated reproduction number of the Tolima CL second outbreak was found to be 1.02 (see Table 6). 1 Estimation using the observed curve of accumulated cases (EAC or M0), 2 Estimation using the observed incidence curve (EIC or M1), 3 Estimation using the observed incidence curve (EIC or M2).

Comparison of Estimates of R 0 between Two Regions
The estimates of the basic reproduction number for Madrid outbreak was around three times more than the estimates of the reproduction number for Tolima outbreaks (see Table 7). The first outbreak of Tolima has larger reproduction number as compared to the corresponding number for the second outbreak. Table 7. Comparison of estimates of the basic reproduction number for the outbreaks.

Discussion
Proper surveillance is crucial for controlling leishmaniasis in endemic countries; however, there is a need to develop methods that can measure disease transmission rates effectively using existing limited data [29] and can be used to evaluate control programs [40]. Leishmaniasis-affected regions are primarily resource-constrained and hence face various challenges to gathering regular comprehensive data. In such scenarios, model-driven decisions might be helpful and can provide understanding of region-specific transmission dynamics [37]. In this study, we provide methodologies to estimate the basic reproduction number, R 0 , for CL with regional dependent factors. The estimation methods were tested using case studies from the two economically contrasting regions, Madrid, Spain and Tolima, Colombia. The Madrid model considers dog reservoir hosts (since most cases had contact with infected dogs) whereas the Tolima model takes into account the movement of military personnel on the transmission dynamics of CL. The three estimation procedures were developed for the two models to estimate their parameters using reported incidence data. Unlike the traditionally used estimating process, in which point estimation of model parameters is taken directly from independent studies reported in the literature, the methods in this research provide a simple but consistent way to estimate model parameters. The estimation of model parameters is followed by the estimation of the basic reproduction number, R 0 , and the computation of various epidemiologically important quantities, the type reproduction number (RT) vectorial capacity (VC) and entomological inoculation rates (EIR).
Prior to the estimation of model parameters, outbreak-related CL incidence data from the two ecologically and epidemiologically different regions (Tolima, Colombia and Madrid, Spain) are first analyzed. Various differences are found in the outbreaks: (i) The Madrid incidence data were in months whereas the Tolima data were collected every week (this suggests a difference in the reporting systems of the two countries and potentially different infectious and latent periods between the regions), (ii) The outbreak in Madrid peaked in winters (Dec and Jan) as compared to the outbreak in Tolima, where the highest incidence was observed in Spring (April) and fall (in October), (iii) Each of the Tolima outbreaks was short lived (1/2 year) whereas the outbreak in Madrid lasted for 3 years (the Madrid outbreak was from 2009 to 2012 whereas the two Tolima outbreaks both occurred during 2016), (iv) there were two successive outbreaks in Tolima whereas there was a single outbreak in Madrid (the first outbreak in Tolima was much more lethal than its second outbreak), and (v) Dog reservoirs were important in the Madrid transmission cycle but in the Tolima outbreak, the frequently moving military population played a critical role in its spread.
The key parameter describing the spread of an infection is the basic reproduction numbers, R 0 , which is defined as the number of secondary infections generated by an infected index case in otherwise susceptible population. This study uses mathematical models to estimate R 0 for the 2009-2012 CL outbreak in Madrid and the two CL outbreaks in Tolima during 2016. The mean estimates of R 0 are found to be 3.1 for Madrid, 1.2 for the first outbreak of Tolima and 1.01 for the second outbreak of Tolima. The R 0 estimate for Madrid seems to be significantly higher than corresponding estimate for Tolima. This could be a result of differences in the population density (60 persons/km 2 in Tolima vs. 5400 persons/km 2 in Madrid), climatic factors, human mobility, and/or health disparities in sub-communities [41]. In the Madrid outbreak, dogs were the main reservoir host, P. perniciosus was the principal vector of leishmania and L. infantum was the primary parasite species [5]. Epidemic outbreaks of CL in Tolima were caused by L. braziliensis, L. guyanensis and L. panamensis, with intra-and peri-domiciliary transmission. The military showed the highest incidence of CL due to the continuous deployment of troops to areas of high endemicity.
Risk factors for CL such as urbanization, malnutrition, health seeking behaviors and disparity have been reported in the literature [42]; however, their impact on the dynamics of CL is less known [43]. This study attempts to provide a simple framework by which the impact of risk factors can be captured using limited reported data. We made some simplifying assumptions to reduce the dimension of the models and to obtain an explicit analytical formula for estimating R 0 . The data used here for fitting the model were obtained through passive case detection and therefore may be prone to high underreporting. Nevertheless, the results in this study suggest a One Health perspective, for example, if animals are key reservoirs of CL, interventions should not be human focused only, instead control programs should be heterogeneous, focusing on both human and animal hosts [44]. It suggests identification of disease-affected communities also in Blue Marble Health countries, wealthy nations with high GDP but also high endemicity of neglected diseases in hidden pockets [41]. Thus, the current estimation study needs to also include proper cost analysis in order to study transmission dynamics comprehensively [45]. As more data accumulate in the future, a more thorough analysis will allow for more accurate estimates of R 0 , together with less uncertainty around them and greater understanding of impact of socio-economic conditions on its estimates.