Leishmaniasis is an infectious disease caused by parasites from Leishmania species. It is considered a neglected tropical disease, affecting 96 countries worldwide (Alvar et al., 2012). More than 700.000 new cases are reported by WHO every year (WHO, 2017). The disease presents in two main different clinical forms, visceral and cutaneous leishmaniasis, and their outcome ranges from self-healing cutaneous lesions to disseminated lesions or to visceral dissemination of the parasite, possibly leading to death if not properly treated (Carvalho et al., 1985; Alvar et al., 2012). The clinical form and severity of leishmaniasis depend on the parasite species/strain involved and on the host immune response mounted (Sacks and Noben-Trauth, 2002; Hurdayal and Brombacher, 2014).

The immune response involved in leishmaniasis has mainly been studied in mice. The animal models most frequently used are BALB/c and C57BL/6 mice infected with L. major, which represent typical models for susceptibility and resistance respectively. The susceptibility in the BALB/c strain is related to development of (Th2)-polarized immune response, genetically characterized by a high production of interleukin (IL)-4 by CD4⁺ T cells (Sacks and Noben-Trauth, 2002; Mougneau et al., 2011). Conversely, the C57BL/6 strain represents a model for resistance to L. major, due to (Th1) polarized immune response, characterized by the production of high amounts of interferon-gamma (IFN-γ) by CD4⁺ T cells (Sacks and Noben-Trauth, 2002; Mougneau et al., 2011). The differentiation of naïve helper T (Th) cells into Th1 or Th2 cells depends on antigen presentation, timing, and cytokines produced by dendritic cells (DCs) after contact with the parasite (Vieira et al., 1994; Hurdayal and Brombacher, 2014). When DCs express nuclear factors such as STAT6 and secrete IL-4 and IL-10, they instruct naïve Th cells to differentiate into the Th2 phenotype (Dent et al., 1999; Sacks and Noben-Trauth, 2002). DCs expressing the nuclear factors STAT4 (Buxbaum et al., 2002), STAT1, and IL-12, induce the differentiation of naïve Th cells into the Th1 subtype, in turn producing IFN-γ (Stamm et al., 1999; Sacks and Noben-Trauth, 2002; Jayakumar et al., 2008).

The immune responses to other species of Leishmania such as L. braziliensis and L. amazonensis are less well studied than that to L. major, and susceptibility mechanisms are not as clearly linked to a Th2 immune response as they are in the L. major model. However, it seems to be clear that control of parasite growth is always dependent on IFN-γ and other inflammatory cytokines. In the L. braziliensis BALB/c infection model, secretion of IFN-γ and TNF-α, is observed but no Th2 immune response, is induced even in the absence of IL-12 (Souza-Neto et al., 2004; Vargas-Inchaustegui et al., 2008). Also, susceptibility to L. amazonensis is related to very low levels of IFN-γ production and low cell proliferation rate in response to its expression of serine phosphate on its membrane and to a high expression of ecto-nucleotidases (Maioli et al., 2004; Franca-Costa et al., 2012).

Parasites including the Leishmania species employ mechanisms to escape the immune response, interfering with signaling pathways of antigen-presenting cells (APCs) and with the differentiation of Th cells (Mougneau et al., 2011). An important escape mechanism employed by Leishmania may be related to the conversion of trinucleotides to adenosine (Ado). It has been reported that ATP leads mostly to pro-inflammatory signals while Ado acts by limiting the inflammation (Bours et al., 2006; Cekic and Linden, 2016). Our group and others have shown that increased ecto-nucleotidase activity on the surface of these parasites correlates with different virulence levels of the cutaneous form of the disease (Maioli et al., 2004; Marques-da-Silva et al., 2008). L. amazonensis has the highest activity of ecto-nucleotidases, which leads to a higher concentration of Ado, decreasing the capacity of DC to present antigens and induce differentiation of Th cells leading to less Th cell proliferation and cytokine production (de Souza et al., 2010; Leite et al., 2012; Figure 1).

The complex dynamics of the Leishmania-host interaction can be addressed with a certain degree of success by using mathematical and computational modeling. A system of ordinary differential equations (ODEs) was proposed (Nelson and Velasco-Hernandez, 2002) and later expanded (Biswas et al., 2016) to describe the dynamics of macrophages and parasites in the early phase of infection prior to the development of the adaptive immune response. Another approach—agent-based modeling (i.e., a class of discrete computational models; Castiglione, 2009)—has been implemented with L. major infection data to describe the dynamics of parasites and macrophages in the later phase of infection (Dancik et al., 2010). The authors found that the decrease in the number of macrophages following peak infection could be explained by their uptake of necrotic tissues. Paradoxically, they also found that a decrease in the parasite reproduction rate might eventually lead to more parasites. In line with these results, a different ODE model showed a negative correlation between parasite load in the initial stage of the infection and the overall number of parasites at the end of the observed time window (Länger et al., 2012). The hypothesis arising from both models is that a smaller reproduction rate elicits a weaker immune response, resulting in higher survival rates of the parasite.

Despite some progress, comprehensive models for Leishmania-host interactions and leishmaniasis progression have not been yet fully implemented: most models so far are based only on specific aspects of the disease such as the interaction between macrophages and parasites. To our knowledge, this is the first model of leishmaniasis that covers the general dynamics of the infection and takes into account the importance of purinergic (i.e., adenosine- and ATP-based) signaling. Here we propose a model of cutaneous leishmaniasis caused by different Leishmania species, L. major, L. braziliensis, and L. amazonensis, which aims to establish a minimum set of rules that can describe the development of infection for each species and to test the importance of Ado release as a virulence factor. We created a model composed of the key immune competent cell types (CD4⁺ Th cells, macrophages, DCs and epithelial cells) and molecules (IL-4, IL-12, IFN-γ, Ado, and nitric oxide, NO) reported in the literature, and the parasite. By modeling such key aspects together, we were able to demonstrate the effect of Ado on the number of parasites and the lesion formation process, thus showing the importance of Ado in the inhibition of inflammatory Th cells in Leishmania infection, as well as in the lesion formation processes.

The model we present here was created with the aim to reproduce the behavior of Leishmania infection from experimental parasite burden curve related to the immune response developed by the host in the mouse model.

In this model, the time-step length was not predefined by the parameters in Equations (2)–(10) but defined a posteriori by superimposing simulated dynamics with real data. In this sense, comparisons with experimental data (Belkaid et al., 2000) serve as a qualitative validation of the model and as a means of extracting values such as time-step length and the starting point of simulation (Figure 3). During the fine-tuning, it was found that the key parameter controlling step-size was Th-cell activation probability (k₁). Any value of k₁ in the range of 0.001–100% probability reproduced published data but with markedly differences in the step-size. With a probability of activation of 100% the time step is close to ~80 h (~half a week), while when such probability is 0.001% the time step is 2 h (80 steps/week). The inverse reasoning can be applied to check that it can make perfect sense. In a long period, such as ~80 h of antigen exposure, the probability of finding a Th-cell specific for leishmania antigen will be close to 100%, while for a smaller period (e.g., 2 h) this probability will be very low at least for low dose inoculations. The other four parameters had an important bias in the outcome of the disease but they were found not to influence the time-step length. It is interesting to notice that k₂ and k₄ are related to levels of inhibitors (IL-4 and Ado) which do not depend on time while k₁, k₃, and k₅ are related to Th-cell activation, and Th-cell survival probabilities and leishmania replication rate., All these depend on time, but only k₁ affects time-step length.

The parameter Th-cell half-life (k₃) is sensitive to Th-cell activation (k₁). With k₁ high (>10%) the k₃ value is nearly irrelevant but for k₁ low (<0.01) k₃ must have a value greater than zero for a cure to be possible. Biologically, this finding makes sense since Th-cells have half-lives of days (Murphy et al., 2008) and k₁ low and k₃ close to zero would lead to an unrealistically short half-life of hours or minutes. It can be also observed that k₃ correlates negatively with the number of parasites throughout the simulation (Table 4). Most experiments used a value of 18% for k₃ (Table 3). This value will give a Th-cell half-live, in the absence of antigen, of about 11.5 steps (1 day), which is roughly close to the half-life of these cells in vivo. Throughout the simulation, the survival of Th cells will also be influenced by parasite presence and Ado production. Entities representing molecules (cytokines, Ado and NO) were modeled with a lifespan of 1 time step (2 h), in agreement with the fact that these substances have half-lives of minutes to hours (Bocchi, 1991; Loffler et al., 2007).

A broad range of values of k₁ could reproduce qualitatively the profile published by Belkaid (Belkaid et al., 2000): indeed Th-cells activation probability in the range of 10-100% leads to a growth in parasite number of the order of 10–10². The same study reported that the number of parasites, following an injection of 10² early developmental stage (amastigotes) parasites, has grown up to 10⁵ late developmental stage (promastigotes) parasites within 4.5 weeks (Belkaid et al., 2000). It is to be noted that, according to von Stebut (2007), 90% of parasites injected are killed by the complement within 3 min: this means a 10³–10⁴-fold increase in the number of parasites. These figures are better reproduced by an activation probability k₁ ~0.001%. The parameter controlling parasite growth rate (k₅) also plays a role in tuning the results to better reproduce literature data. In this work, the pair of values k₁ = 0.001% and k₅ = 35% was very successful in reproducing published data.

The comparisons made in Figure 3 also showed that the simulation starts at 2.5 weeks post-infection. The key reason why it could not simulate the first 2.5 weeks of the disease is that this model is a quadratic approximation of the real phenomenon. As an approximation, it cannot reproduce the whole dynamics but only a certain range, similar to what one would expect with a Taylor linearization over a critical point. Besides that, 3 lattice-points containing parasites can be translated to 10² parasites in a mouse which is not very far from the expected value after 2.5 weeks post-infection given that only about 10 of 100 injected parasites survive after 3 min (von Stebut, 2007). A similar approach of only simulating the late phase of Leishmania infection was used by Dancik et al. (2010); the simulations in his work start at 3.5 weeks post-infection. Nevertheless, the present model simulates part of the early silent phase of infection that according to Belkaid et al. (2000) lasts for the first 3 or 4 weeks.

It was possible to expand a L. major model to simulate other cutaneous leishmaniasis (L. braziliensis and L. amazonensis). As it can be seen in Table 3 there is no IL-4 induced production of IL-12 for these other leishmaniasis (for a review on this see Hochrein et al., 2000). This adequately models the finding of Maioli et al. (2004) that no IL-4 producing cells could be detected during L. braziliensis or L. amazonensis C57BL/6 infection. An increase in Ado concentration to inhibit the activation of Th-cells was used to simulate L. amazonensis infection, in agreement with literature reports of incurable lesions in C57BL/6 mice infected with L. amazonensis, and with higher ecto-nucleotidase activity in the surface of this parasite (Berredo-Pinho et al., 2001; Ji et al., 2003; Maioli et al., 2004; Marques-da-Silva et al., 2008; de Souza et al., 2010; Gomes et al., 2015). This higher activity may lead to higher Ado and lower ATP concentration at the site of infection.

Notice that for simulation of L. braziliensis infection in C57BL/6 mice it was not necessary to tune the parameters of the model, but only to turn ϕ_IL12 on; this change was nevertheless extracted directly from literature (Souza-Neto et al., 2004; Vargas-Inchaustegui et al., 2008). In the case of L. amazonensis only Ado inhibition over Th-cells (k4) was changed.

Comparison of the lesion size modeled in these simulations with experimental data from Marques-da-Silva et al. (2008), Maioli et al. (2004), Ji et al. (2003) (Figures 4A–G) showed that the model is capable of simulating different mouse models of cutaneous leishmaniasis. These comparisons not only validate the new simulations (L. braziliensis and L. amazonensis) but they also serve as a further validation for the L. major model. It is interesting to notice that the simulations were fitted to these other published data (Figures 4A–G) using the step-size found with Belkaid et al. (2000) comparison (80 steps/week). The appearance of the lesion takes a few weeks in mice. Thus, the starting point had to be changed from 2.5 to zero weeks post-infection. A justification for that is that in Belkaid's article (Belkaid et al., 2000), the infection was induced with a low dose of inoculum (10² parasites injected in the ear), while in the other studies in our comparison the infection was induced with high dose inoculum (10⁵–10⁶ parasites injected in the footpad). It is therefore reasonable to assume a different starting point for the appearance of the lesion. It is interesting that, according to the simulations, the inoculum size does not change the dynamics of late infection.

It is important to highlight that the comparisons in Figure 4 are only meant to be qualitative, since different lesion metrics (lesion area vs. thickness) are used. The aim of these comparisons is to see if the model can capture key aspects of lesion dynamics such as when they appear and disappear and the time when they reach a peak. Data in Figures 4A–G were normalized by the size of L. amazonensis infection lesion in C57BL/mice at 6 weeks post-infection. This approach made it possible to compare the relative sizes of different cutaneous leishmaniasis lesions. We observed a good agreement between the model and literature with L. amazonensis lesions, i.e., the largest lesion, followed by L. major, and by L. braziliensis, which induces the smallest wounds in C57BL/6 mice.

While k₁ is the key value for setting the time-step length, Ado inhibition over effector cells (k₄) is the most important parameter for simulating the outcome of the disease. This can be observed in the sensitivity analysis (Table 4) and in Figures 5, 6. This shows that this parameter has the most influence over the number of parasites throughout the infection and it also controls whether the infection will be cured. These results agree with Marques-da-Silva et al. (2008) and correlate the different outcomes between the three leishmaniasis in this study with ecto-nucleotidase activity on the surface of parasites and a correlation of ecto-nucleotidase activity with the virulence of different strains of L. amazonensis (de Souza et al., 2010).

Lesion growth is influenced by Ado concentration as is shown in Figures 6, 7. These results show that this purine has a non-linear relationship with growth. Moderate inhibition of Th-cells activation may lead to the largest lesions, in agreement with Marques-da-Silva et al. (2008), Maioli et al. (2004) and Ji et al. (2003). But high Ado concentration would have a stronger inhibition of Th-cell activation and it may lead to smaller and more diffuse lesions (Supplementary Material). These data may correlate with the different outcomes of L. amazonensis infection in humans (Silveira et al., 2009; Hombach and Clos, 2014). It may be the case that for some genetic backgrounds, the immune system responds strongly to Ado with more secretion of IL-10 and TGF-β leading to diffuse cutaneous leishmaniasis outcome. Nevertheless, such effects remain speculative and more studies are needed to fully decipher the influence of Ado in human leishmaniosis.

This model also includes the role of IL-4. This cytokine will induce the production of IL-12 (Hochrein et al., 2000) and at the same time inhibits the action of IL-12 in Th-cells (Szabo et al., 1997). It is reported in the literature that if IL-4 is presented in the earlier phase of leishmaniasis it leads to the production of IL-12 and prime Th1 cells, while if it is presented during Th-cells priming it will lead to Th2 production and will suppress the expression of the IL-12 receptor (Biedermann et al., 2001). Furthermore, Belkaid et al. (2000) showed that IL-4 is produced during the whole infection (weeks 1–22) and that it peaks during the Th-cell priming phase (weeks 4–8) without impairing the Th1 response. The apparent contradiction between these two results can most likely be explained by the dose of IL-4. Biedermann and coworkers (Biedermann et al., 2001) injected 1 μg of IL-4 into the mice, which in a conservative calculation will lead to 10⁴ pg/ml of that cytokine, while Belkaid (Belkaid et al., 2000) measures physiological doses of 10–100 pg/ml in the lymph node. The current work simulates physiological conditions that are more in line with what Belkaid (Belkaid et al., 2000) showed. Nevertheless, sensitivity analysis showed that IL-4, despite stimulating the production of IL-12, promotes susceptibility to infection, in accordance with the literature (Biedermann et al., 2001; Cangussú et al., 2009).

It was observed in the animations of L. major infection (Supplementary Material) that there is a spatial negative correlation between L. major parasites and NO production. This finding agrees with results from (Cangussú et al., 2009).

We have shown how a robust model of cutaneous leishmaniasis has been built and validated. It consists of a minimum set of rules that can describe and differentiate L. amazonensis, L. braziliensis, and L. major infection. The simulator agrees with key experimental results published in the literature such as the antagonism between IL-4 and IFN-γ and brings new lights to the influence of Ado signaling in these infections. We could observe a non-linear relationship between this purine and lesion formation and to confirm that it is the key factor in differentiating cutaneous leishmaniasis. Higher Ado concentration can inhibit Th-cell activation leading to infection outcomes that have been seen in experimental models of leishmaniasis. This model can therefore represent a valuable tool for answering questions regarding cellular and molecular players, inflammatory processes related to these infections and finally exploration of combination therapies involving drugs such as sodium stibogluconate, topical paromomycin preparations, etc.

HR, TM, LF, PT, and FC conceived of the study, HR and FC ran the simulations, all authors contributed to the analysis, wrote the manuscript and finally read and approved the final manuscript.