Modeling of leishmaniasis infection dynamics: novel application to the design of effective therapies

Background The WHO considers leishmaniasis as one of the six most important tropical diseases worldwide. It is caused by parasites of the genus Leishmania that are passed on to humans and animals by the phlebotomine sandfly. Despite all of the research, there is still a lack of understanding on the metabolism of the parasite and the progression of the disease. In this study, a mathematical model of disease progression was developed based on experimental data of clinical symptoms, immunological responses, and parasite load for Leishmania amazonensis in BALB/c mice. Results Four biologically significant variables were chosen to develop a differential equation model based on the GMA power-law formalism. Parameters were determined to minimize error in the model dynamics and time series experimental data. Subsequently, the model robustness was tested and the model predictions were verified by comparing them with experimental observations made in different experimental conditions. The model obtained helps to quantify relationships between the selected variables, leads to a better understanding of disease progression, and aids in the identification of crucial points for introducing therapeutic methods. Conclusions Our model can be used to identify the biological factors that must be changed to minimize parasite load in the host body, and contributes to the design of effective therapies.


Background
The WHO considers leishmaniasis as one of the six most important tropical diseases worldwide [1]. It is caused by parasites of the genus Leishmania that are passed to humans and animals by sandflies of the subfamily Phlebotominae [2]. Leishmaniasis, which is endemic in 88 countries, has an annual incidence of two million cases and is estimated to cause over 50,000 deaths per year [3]. The disease has three main forms: cutaneous leishmaniasis, mucocutaneous leishmaniasis and visceral leishmaniasis. Visceral leishmaniasis, the most severe form of the disease, is also known as "kala azar", "black fever" or "dumdum fever". It especially affects hosts with weak immune systems, such as children or adults suffering from malnutrition or HIV. The human immune response that limits leishmaniasis is mediated by Th1 cells that activate macrophages to kill the parasite (cellular immunity). When cellular immunity is deficient, an expansion of Th2 cells occurs which allows the parasite to survive within the monocytes and fosters disease development [4]. After an incubation period that varies from ten days to two years [3], typical symptoms are fever, diarrhea, body weight loss, lymphadenopathy, hepatomegaly and splenomegaly. Despite all of the current research, there is still a lack of understanding about the metabolism of the parasite and disease progression.
Mathematical modeling of the processes involved in parasite-host interactions has become a necessary tool in the study of diseases, leishmaniasis being no exception. A significant part of the modeling work in this field is epidemiological [5][6][7]. In addition, models have also been used to study regulation of gene expression, protein synthesis, and metabolism of the parasite at the genome-wide level [8][9][10]. The dynamics of parasite-host interactions in the infection process has also been studied using agent-based modeling approaches [11][12][13]. For example, Dancik et al. [13] used an agent-based stochastic model of the immune response of mice to L. major infection to identify parameters that are important in changing the dynamics of the infection process, and to quantify the influence of those parameters. The authors showed that increasing parasite growth rate decreases pathogen load in some circumstances.
There are many studies regarding the biology, epidemiology and immunology of leishmaniasis [5,6,14,7,15], yet there are fewer studies related to the evolution of the infection in animal models. A significant reference for the latter is the work of Courret et al. [16] where lesion development, cellular response, expression of cytokines, as well as parasite load in the spleen of BALB/c mice infected with L. amazonensis is described. In this vein there are also the works of Arrais-Silva et al. [17] on the hypoxia-inducible factor-1 from L. amazonensis infection; of Lira et al. [18] on BALB/c mice symptoms, parasite load and immune response in C57BL/6 mice infected with L. major; and the work of Requena et al. [19] and of Dea-Ayuela et al. [20] that explores the clinical symptoms, parasite loads and antibody levels in susceptible, oligosymptomatic and resistant hamsters. The study of Requena et al. [19] compared these parameters together with lymphocyte population and proliferation, in two groups infected with different amounts of parasites and a control group.
In the present work we adopted a systems biology approach for understanding disease evolution, hostpathogen interactions, and immune response function. We performed this task by using experimental time series measurements in BALB/c mice infected with L. amazonensis to parameterize a mathematical model that accounts for immune response and parasite load. Based on this model we were able to quantify the biological interrelations between variables, perform predictive simulations, carry out sensitivity analysis to evaluate the significance of the system parameters, and solve an optimization problem for minimizing the parasite load. This analysis contributes valuable information to the drug discovery pipeline for developing effective therapeutic methods against leishmaniasis.

Mathematical Model
Experimental measurements obtained in BALB/c mice were used to fit the parameters of the mathematical model shown in Figure 1 that accounts for the progression of Leishmaniasis (see Methods for details). In the present work we used a general power-law formulation, the General Mass Action Power Law formalism (GMA-PL) that allows for non-integer kinetic orders [21,22] with the following structure: In the above expression, X i , σ ij , γ j and g jk represent the normalized variable set, the stoichiometric matrix, the rate constants, and the kinetic orders, respectively. The variables lymphocytes proliferation (X 2 ), IgG1 (X 3 ) and IgG2a (X 4 ) were normalized with respect to the respective value in the control group of mice. Because the control group is parasite-free, the same approach could not be used to normalize parasite load. In this case the variable was normalized with respect to its own mean value. This standardization reduces the range of variation of the parameters and computation time, and also exploits various properties of the GMA-PL formalism on the behavior of variables and parameters. The specific numerical values for the parameters s ij , γ j and g jk are determined using prior biological knowledge, information about the basal steady-states of the system [21], and/or dynamical experimental data [22,23]. In powerlaw models, kinetic orders can have non-integer values. One of the main advantages of power-law models is that they allow for the condensation of several steps into simplified representations [21,24,25]. The parameters of the model are kinetic orders and rate constants. Negative values for the kinetic order represent inhibition, that is, an increase in its variable leads to a diminution of the rate involved, while a zero indicates that the variable does not affect the described process. When positive values are considered for a kinetic order, several alternatives are possible: values between zero and one imply a flux that depends on the variable in a saturating-like manner. Values equal to one imply a flux that depends linearly on the variable or, in chemical terms, a first order reaction. By allowing non-integer, positive or negative, kinetic orders, we are able to consider a larger class of kinetic models from which we can select a suitable candidate without changing the (original) model structure. Figure 1 shows the model scheme of the chosen variables and the influences among them denoted by arrows and parameters. g 1 , ... , g 14 stand for kinetic orders representing influences on the creation or degradation fluxes (V i ) of the four variables.
The total parasite load in the host (X 1 ) stimulates its immune response. The parasites multiply in macrophages by binary division. The parasite load growth (V 1 ) has a nonlinear dependence on the parasite load through the kinetic order g 1 . Increased parasite load leads to a decrease in the proliferation rate of lymphocytes (V 4 ); this interaction is represented by the kinetic order g 7 [26]. Proliferation (or multiplication) of T lymphocytes (X 2 ) occurs when naïve T cells are activated by antigens of the pathogen (g 5 ) and then differentiated into effector cells (Th1 or Th2) and memory cells. The activation of lymphocytes is an essential event in the production of specific immune responses (both humoral and cellular) against pathogens. Proliferation was measured following the protocol by Monks et al. [27] using the Stimulation Index (see Methods). The lymphocyte proliferation V 3 is also stimulated by X 2 , through g 6 . Cell mediated effectors enhance X 1 decay (V 2 ); this effect is represented by the positive kinetic order g 3 [28,29]. The host immune system produces IgG1 (X 3 ) and IgG2a (X 4 ) antibodies which could be linked to the Th2 and Th1 mechanisms respectively [30,31]. This is represented in our model through a positive influence of X 2 on the rate synthesis of IgG1 (V 5 ) through g 9 , and on the rate synthesis of IgG2a (V 7 ) through g 12 . These two immunoglobulins are antagonistic, so each of them has a negative influence on the generation rate of the other, namely X 4 on V 5 and X 3 on V 7 . These effects are represented by the kinetic orders g 10 and g 13 , respectively [32,33]. The IgG2a influences macrophage activity by stimulating the X 1 rate decay, V 2 . This interaction is represented in our model by the positive kinetic order g 4 . It is assumed that the transformation rates V 2 , V 4 , V 5 and V 8 are proportional to X 1 , X 2 , X 3 and X 4 . These dependences are represented in the model by the positive kinetic orders g 2 , g 8 , g 11 and g 14 , respectively.
Given that every variable has an influx and an outflow, the stoichiometric coefficients are 1 and -1 for the synthesis and transformation processes respectively. Model parameters were determined by fitting the model to experimental data from mice using a genetic algorithm as described in the methods section.
Accordingly, the power law model derived from the above scheme is given by: (2) Figure 2 shows the model data fitting, displaying a good correlation between the experimental and estimated data.

Model validation
We validated the model by using it to make predictions about the way the system would behave under initial parasite loads that were different from those used to calculate parameter models (106 as compared to 103, see Methods for details). We then performed the corresponding experiments in vivo (see Methods), measuring  the four variables described by the model and comparing their observed behavior to model predictions. This initial number of parasites (which mimics a severe leishmaniasis condition) was chosen to check the model capacity to correctly describe the behavior in extreme and differing conditions of initial parasite load. Since the model's main purpose is for the design of therapeutic strategies, a model able to describe a wide range of parasite load dynamics is of foremost interest. Figure 2 shows the results obtained. There it can be seen that the deviation of model prediction from the experimental data is reasonably good during the initial 20 weeks after infection (which means that new treatments are applicable in this time period) with a parasite load of 10 6 .
Since the model describes the evolution over the first 20 weeks after infection, the observed discrepancies in the two experimental conditions considered (model fitting and validation) can be deemed as reasonable in light of the associated experimental error. In this regard, we want to stress the fact that in the experimental data used for model verification, other elements of the immune system may be playing a significant role not addressed by the model, but which could be relevant in conditions of massive infection. Figure 3 shows the dynamic sensitivity analysis for the kinetic orders (g i ) and rate constants (γ j ). The System Parameter Dynamic Sensitivity is noted as S pk Xi ; pk is the parameter under scrutiny and Xi is the considered variable (see equation 6 in Methods). As is shown, all sensitivities have absolute values between 4·10 -5 (S(g 10 ; X 2 )) and 1.83 (S(γ 1 ; X 1 )). The lower value corresponds to the influence of g 10 on the lymphocyte proliferation, and the higher value measures the influence of the rate constant associated to the parasite multiplication rate on the parasite load. This range of values, together with the observation that the median of the absolute values of sensitivities is 0.067, indicates a robust model.

Sensitivity Analysis
In general, the higher sensitivities correspond to the variable indicated by the arrow, except for the parameters directly influencing lymphocyte proliferation (γ 3 , γ 4 , g 5 , g 6 , g 7 and g 8 ). For the parasite load, sensitivities with absolute values greater than 1 are S γ1 X1 and S g1 X . This implies that the generation rate of parasites, and the effect of parasites on their own generation, strongly influence parasite load. Sensitivity S γ2 X1 is much lower than S γ1 X1 . All the other parameters yield sensitivities with absolute values of less than 0.1 for parasite load. The sensitivity of IgG1 and IgG2a to most of the parameters is higher than the sensitivity of lymphocyte proliferation or parasite load to the same parameters. This could be a consequence of the fact that most of the values for parameters directly influencing the parasite load (γ 1 , γ 2 , g 1 , g 2 , g 3 , and g 4 ) are small (< 1), as opposed to those directly influencing lymphocyte proliferation.

Systematic search of parameter profiles for the minimization of parasite load
In order to apply the model for therapeutic purposes, we carried out a systematic search of parameter values that minimize the parasite load. The aim was to discover Figure 3 Absolute value of the system dynamic sensitivities S pk Xi (X j p k ). The search program was organized in two phases. In the first phase, we changed only one of the parameters at a time (g i or γ i ), with the others maintaining their original values. In this case the value of the candidate parameter is initially set to 10% of the model estimated parameter, the following to 20%, 30% and so on, until the parameter reaches the upper-bound range that was assumed feasible and physiologically relevant. Then, for each changing factor, the model solutions were calculated. In order to evaluate the effectiveness of the change in parameter value, the mean, maximal and final parasite loads were calculated. The mean parasite load reflects the average severity of the disease, the maximal value accounts for the maximal number of parasites along the infection dynamics (which has to be lower than the maximum number of parasites the organism can bear), and the final parasite load represents the final outcome of the disease.

Single-parameter search for identification of optimum parameter values g i parameter scanning
After a systematic search among all kinetic constants, we found that g 1 and g 6 were the most suitable parameters to be changed for reducing parasite load. g 1 , describes the influence of parasites on their own proliferation and is the most significant in this regard, since changes in its value causes the greatest reduction. This is achieved by increasing g 1 value from 0.53 to 3. Figure  4A shows the progression of parasite load during a time period of 24 weeks for values of g 1 ranging from 0.01 to 3. Our model suggests that increases in g 1 from its initial value (0.53) have a therapeutic effect, because they lead to a decrease in parasite load and therefore to healing. It is important to note that for values of g 1 < 0.5, final parasite load is proportional to g 1 . However, for values of g 1 > 0.5, final parasite load becomes inversely proportional to the value of this parameter. This latter region includes the actual g 1 value. This implies that therapeutic strategies should aim to increase g 1 . If decreases are sought, such decreases must be well below 0.5 in order to have a similar effect. Figure 4B shows that the evolution of the parasite load decreases right from the beginning, almost linearly through the 24-week period.
At first glance this result could appear paradoxical and certainly counterintuitive. But, in fact, the sensitivity of g 1 with respect to the parasite load (X 1 ), S g1 X1 , is negative (see Figure 2), therefore an increase in g 1 should lead to a decrease in parasite load. This prediction holds true for the mean parasite load over a time period of 24 weeks. However, it has been observed that both the maximum and the final value of parasite load increase in the beginning until they reach a threshold value and decrease from that point on. Again, such a decrease in parasite load indicates that increasing g 1 produces a therapeutic effect. The same result has recently been observed by Dancik et al. [13]. Their model showed that increasing parasite growth rate (V 1 in our model, which is influenced by g 1 ) impairs the pathogen load in certain stages of the disease. Since the modeling strategy used by Dancik et al. [13] is different from our approach, there is not a straightforward translation of the parameters. However, there are some correlations between their findings and our model predictions. In particular, their analysis [13] indicates that increasing growth rate A B can, in some circumstances, suppress pathogen loads, which is what our model also predicts. Due to the particular formalism we have used, we are able to point to the increase of g 1 as the mechanism that increases parasite growth rate. They reported the same evolution pattern in the parasite load that we found: a higher parasite growth rate yields a higher increase in pathogen load in the beginning but also a higher decrease afterwards in such a way that, as a whole, parasites are eliminated earlier. They observed that infection was cleared after eight weeks versus the 17-20 weeks in our model, but this can be attributed to the different Leishmania strain and the lower initial parasite load (1000 versus 50). The authors explain this behavior by stating that a pathogen that proliferates rapidly is more likely to be detected by the immune system. Therefore, pathogen load decreases as growth rate increases, with slowly replicating pathogens persisting longer than fast growing ones. Another g i parameter with similarly minimizing effects on the parasite load is g 6 . g 6 stands for the influence of lymphocytes on their own production. It has been observed that the parasite load (final, maximal and mean) can be reduced by increasing g 6 from its original value of 0.02 to 1.02 ( Figure 5). This figure shows that parasite load first increases until reaching a maximum after 18 to 20 weeks and then abruptly decreases. The effect of the augmentation of g 6 could be related to system immune response enhancement. Lymphocytes need some time to identify the pathogens, thus there is a time lag between the start of immune response, the identification of parasites, and their elimination.
Accordingly, the parasite load augments until reaching a point where it suddenly decreases. g i parameter scanning We also carried out a systematic search among all rate constants. We observed that changes in all g i parameters reduce the final parasite load (see Figure 6). The rate constant that yields the minimum for final, maximal and mean parasite load is g 2 (after slightly increasing its value), which is the rate constant for parasite degradation (see Figure 1). Figure 6A shows final parasite load for different values of g 2 . Figure 6B shows the effect on the parasite load for two increased values of g 2 : 0.43 and 4.3.

Combined two-parameter searches for identification of optimum parameter values
We carried out a systematic scanning of all the combinations of two parameters that yielded the minimum final parasite load. The rationale is that a combination of drugs makes it possible to reduce the parasite load in greater quantity, more quickly, and with lighter dosage than using only one drug. The search was limited to smaller parameter changes in the range of 60% -180% of a parameter's original estimated value.
We found that a total reduction of the observed final (as well as the maximal and mean) parasite load can be attained by simultaneously increasing g 1 from its original value by a factor of 1.6 (approximately), and by changing any other of the remaining 21 parameters by different factors (Table 1) (results not shown). Figure 7 shows all possible combinations of two parameters and the optimum predicted final parasite load. Aside from the best combinations (see Figure 7, black column and file) g 1 and g 2 are also good choices for the minimization of parasite load. The remaining combinations also produce a reduction of parasite load, but to a lesser degree, and are considered to be of minor interest. By way of illustration, Table 1 show the parameter change factor for the two-parameter combinations involving g 1 (black column and file in Figure 7). Interestingly, the solutions that lead to the lowest parasite loads have values for most parameters that are approximately 80% higher than their basal values. Exceptions include γ 1 , γ 4 , γ 6 , g 2 , g 3 , g 5 , and g 14 . Their values are reduced by factors ranging from 5 to 85% (see Table 1).

Discussion
The standard leishmaniasis treatments are chemotherapy based, though some new treatments are based on the use of immunotherapy. In our model, the chemotherapeutical agents are those that target parasite destruction (g 2 ) or inhibit proliferation (g 1 ), whereas immunotherapeutic treatment implies changing parameters g 3 , ... g 8 and g 3 , g 4 , g 6 , g 8 , g 9 , ...g 12 . In most cases the exact interaction mechanism of the drug is not yet known, though it is possible to associate them to the corresponding parameters that are being influenced. It is important to mention that if a given therapeutic agent has an influence that is not represented by any of our model's parameters but corresponds with the in-or outflux of a model variable, the effect of this agent can be translated in our model by a change in the respective rate constant γ i .
Regarding drug therapy, we have found three parameters which cause parasite load reduction: g 1 , which describes the influence of parasites on their own proliferation; g 6 , which represents the influence of lymphocytes on their own proliferation; and γ 2 , the rate constant for parasite degradation.
Examination of the standard drugs used for leishmaniasis treatment shows that most are aimed at parasite destruction. In our model that translates a an increase in γ 2 , the rate of parasite destruction [34], an A B  The Parameter combination columns stand for the combination of parameter while the FC(g 1 )/FCp k columns indicate the ratio between the factor change (FC) of the original value or g1 over the FC of the corresponding parameter (p k ) observation that is coherent with our model's predictions. This is the case for several standard treatments such as amphotericin B (partially inhibits the completion of the parasite's membrane), antimonials (decrease biosynthesis of energy in the amastigote), and itraconazole and pyrazolopyrimidines (inhibit the parasite growth). Other substances currently under evaluation, such as betle leaves extract (reduces viability of promastigotes), interferon (actives macrophages that reduce parasite load), and IL-12 (stimulates Th1) also increase γ 2 . These observations constitute a pragmatic, a posteriori verification of our model's predictions.
Most of the therapeutic drugs used also seem to inhibit, albeit through different mechanisms, parasite proliferation: aminoglycosides alter parasite messenger RNA, pentamidine inhibits polyamine and DNA synthesis in the parasite, imidazole and itraconazole inhibit demethylation of membrane, and pyrazolopyrimidines block protein synthesis and destroy parasite RNA. All these effects can be interpreted, in terms of our model, as a decrease in g 1 . The discrepancy in our model's predictions can be explained by several facts. First, in all cases where a decrease in g 1 could be assumed, there is also the concomitant effect of increasing γ 2 , as noted above. Thus, a trade-off of these two actions should be previously evaluated in order to have an accurate account of the whole drug effect. Second, it should be taken into account that the effect of a g 1 modulation could be different depending on the stage of the disease. It has also been shown that if parasites replicate quickly, the immune system is able to recognize them more easily [13]. Parasites use mechanisms like inhibition of antigen presentation to escape immune response, however, a high growth rate induces massive macrophage recruitment [13]. At this point it should be stressed that our model considers the infection from the very initial stage. A third explanation could be that the parasites produce a certain molecule that stimulates an immune response of the body.
(While investigating a model of tuberculosis infection [35], it was also found that the partial rank correlation between growth rate and extracellular bacteria is negative in a certain time interval.) The factor that increases the influence of parasites on their own proliferation (g1) is crucial according to our model's results; and currently, no pharmaceuticals that increase g1 have been tested against Leishmania. Insuline-like growth factor 1, interferon, and possibly TNF-α cytikine could be considered as potential targets for stimulating parasite replication inside macrophages, and it would be of great interest to test their anti-leishmanial effectiveness. Insuline-like growth factor also increases the number of parasites (γ 1 ) and reduces parasite-toxic production of nitric oxide (γ 2 ). Furthermore, no existing drug is known to have an effect on g 6 , which, in our analysis, is also seen as a possible effective pharmaceutical target. This clearly points to the new, potential application of existing and current therapeutic strategies.
The approach used for detecting key processes that must be regulated in order to reduce parasite load also allowed us to identify combinations of two drugs that would eventually be more effective than a single drug treatment. As is showed in Figure 7, combinations of drugs able to increase g 1 or g 2 and simultaneously change any other parameter, or, alternatively, combinations of drugs that decrease g 1 together with the change in another parameter, would cause significant reduction in the final parasite load. These findings greatly amplify the number of therapeutic options available, although they still remain to be tested. By way of illustration, we could suggest the combination of any of the available drugs that increase g 2 (amphotericin B, aminoglycosides, antimonials, pentamidine imidazole, itraconazole, and pyrazolopyrimidines) together with any of the following: interleukin-5,6,13 and MHC class II molecules (both increasing g 3 ), rLmSTI1 (increase in g 4 ), and chemokines (that increase g 3 and g 2 simultaneously). In the same mouse model we will test the effects on the variables of different drug combinations to verify the model's predictions and to eventually refine and extend the model by including new variables and mechanisms.
A limitation of the present approach is that our model is a simplification and does not include a detailed description of all the factors involved in the interaction mechanism of the drug in the body. However, given that these mechanisms are often not known, the modeling approach constitutes an approximation to the understanding of a complex dynamic system based on available information and informed hypothesis.

Conclusions
In the present work we have illustrated a novel approach for the design of effective therapeutic strategies for leishmaniasis treatment. The approach is based on the integration of experimentally available information on infection development in an animal model using a mathematical model that describes the system dynamics observed. Many of the predictions concur with the standard practice, while others remain to be explored. We are confident that this rational, model-based approach is of great interest given that it overcomes the limitations of a trial and error strategy, and provides an extra layer of rationality in the search for new therapeutic formulas. This approach is also readily applicable to other parasiticrelated illnesses.

Methods
Mice BALB/c mice, 6-8 weeks old, were obtained from the animal breeding facilities of the Universidad de La Laguna. The experimental protocols used were approved by the Animal Care and Use Committee of the University of La Laguna (Approval ID number 132).

Parasites and experimental infection of mice
Amastigotes of MPRO/BR/77/LTB0016 strain of L. amazonensis were prepared from infected BALB/c mice, retaining full virulence. Promastigotes derived from tissue amastigotes were then grown in Schneider's insect medium, pH 7.2, supplemented with 10% heat-inactivated fetal bovine serum, 100 U/ml penicillin and 100 μg/ml streptomycin at 26°C.
In order to follow up the evolution of the infection, 10 3 (10 6 in case of the verification experiments) stationary phase L. amazonensis promastigotes contained in 30 μl of PBS, were injected subcutaneously into the tarsi of right hind leg of 20 BALB/c. In 8 of the 20 mice the response of IgG1 and IgG2a antibodies was evaluated every 2 weeks during the 28-week period of study. The cellular response and parasite load were also measured, for which we had to sacrifice four mice at 8, 16 and 24 weeks of infection.

Parasite quantitation
Estimates of the parasite number present in infected organs were done as described in Buffet et al. [36] which allowed for quantifying the parasite load from tissue homogenates.
Serial threefold dilutions ranging from 1 to 1/3 × 10 -6 were prepared twice for each homogenate in wells of 96well plates, containing 200 μl of Schneider's insect medium, pH 7.2, supplemented with 10% heat-inactivated fetal bovine serum, 200 U/ml penicillin and 200 μg/ml streptomycin. After 7 and 10 days at 26°C, each well was examined and defined as positive or negative based on the presence or absence of viable promastigotes. A limiting dilution analysis was applied to the data to estimate the number of viable Leishmania expressed as Limiting Dilution Assay Unit. The total number of parasites per gram (parasite load) was calculated by equation 4, where D+ is the last dilution positive and Po is the weight of the piece of tissue.

Antigen
The soluble antigen of the parasite (LaAgS) used for enzyme-linked immunosorbent assay (ELISA) determination and splenocytes proliferation assay was obtained from stationary phase cultures of MHOM/BR/77/ LTB0016 strain promastigotes of L. amazonensis and according to Larreta et al. [37].

Antibodies
The specific antibody response levels of IgG1 and IgG2a against LaAgS were determined by indirect ELISA in serum of BALB/c mice. ELISA assays were carried out using standard conditions. Microtiter plates were coated with 0.8 μg per well of the antigen. The sera from the mice were assayed at 1:80 dilutions. As secondary antibody the HRPO-conjugated goat anti-mouse IgG1 and IgG2a were used to 1:8000 and 1:1000, respectively. We used two groups of 8 mice, one experimental and one control.

Lymphoproliferation
In vitro lymphoproliferation assays were carried out to measure the capacity of the LaAgS to stimulate the lymphocyte multiplication or proliferation, as indicative of the capacity of the parasite to produce a specific immune response. The lymphocytes were aseptically removed from spleens of experimental and control BALB/c mice and disrupted in PBS with 1% FCS. The cells were centrifuged, the erythrocytes were lysed in a lysis solution (150 mM NH 4 Cl, 10 mM KHCO 3 , 1 mM EDTA, pH 7.4) and remaining cells were finally resuspended to a density of 2.5 × 10 6 cells/ml in DMEM containing 10% FCS, 2 mM L-glutamine, 0.05 mM 2mercaptoethanol, 12 mM HEPES, pH 7.1, 100 IU/ml penicillin and 100 μg/ml streptomycin. Lymphocytes were divided into 100 μl aliquots (2.5 × 10 5 cells) in 96well plates and they were allowed to proliferate for 3 days at 37°C in an atmosphere containing 5% CO 2 and 95% humidity in the presence or absence of LaAgS (final concentration 40 μg/ml). Proliferation was measured by SRB assay, following the protocol by Monks et al. [27]. Absorbance due to the incorporation of the SRB coloring to anionic proteins of viable lymphocytes, as an index of proliferation stimulation, was measured using a micro plate reader at 570 nm (Model 680; BIO-RAD). Stimulation index of the lymphocyte proliferation (SI) of experimental mice was calculated according the following equation: where X Abs Tw and X Abs Cw stand for the average absorbance value for the antigen treated wells and the average of the average absorbance value for the control wells (only DMEM), respectively. The Cutoff point was calculated as the average difference between the treated wells absorbance minus the control wells absorbance plus 3 SD (standard deviation) of the absorbance of lymphocytes of control mice. The Cutoff was established for each assay. We used two groups of 12 mice, one experimental and one control. All assays were performed in triplicate with four mice representing each group.

Parameter estimation
In the power-law model used, the parameters of the model, kinetic orders (g i ) and rate constants (γ i ) were estimated from experimental data using a genetic algorithm adapted for power-law models [22,38]. The initial population of solutions is generated through a random exploration of the search space, which is defined using feasible intervals of values for the parameters. The best individuals of the population are selected in the considered iteration based on the value of the following objective function: where n var is the number of variables monitored and n tp the number of time points when each variable was measured. In turn, X j (t i ) is the predicted value for the j th variable at the i th time point obtained after numerical integration of the solution, while X j exp (t i ) is the value of the j th observable variable at the i th time point measured in the experiment. This objective function has been tested in other case studies [22,38] showing that it is independent of time and therefore will not be affected by the measurements schedule. We consider it to be well suited for this type of fitting problem. The objective function takes into account the 24-week period for which there is available experimental data of all variables; in the fitting process all variables are weighted equally. The fitting process stops if either the objective function converges below a threshold value or the maximum number of iterations is exceeded. The first criterion is fulfilled when it reaches objective function values smaller than a previously defined one. The second occurs when it reaches a certain number of generations. However, it is possible to set up a combination of both so that parameter estimation stops when a solution with an objective function value below a predetermined parameter is reached, or a certain number of generations have been explored. Every few generations a controllable number of individuals of the total population are selected, and the objective function-defined surface is submitted to local scrutiny taking as starting points selected individuals. For each starting point, parameters are changed in a small and random amount and the new objective function is evaluated. If the new value is lower than its initial value, the new set of parameters is adopted, otherwise the original parameters are kept. In our present case we used a combination stopping criterion, namely 1000 generations, as the maximal number of iterations, and an objective function value smaller than 0.2.
Parameter estimation was executed from the previously obtained time series of experimental data after normalization [39]. All γ i were permitted to vary within the range [0,10] and the g i within the range [0,3], except the parameters g 2 , g 8 , g 11 and g 14 which were set to 1 since they relate variables with their own outflow (V 2 , V 4 , V 6 and V 8 , respectively). This model hypothesis is a biologically plausible one that permits reduction in terms of the number of parameters. These ranges come from physiological as well as kinetic considerations. Processes in which kinetic orders are greater than 3 are rarely observed since this corresponds to processes or reactions with extreme sensitivity to changes in one of the reactants. A similar reasoning applies to higher γ i values, although in this case the admissible range is wider. However, in the optimization search we used greater ranges. The rationale is that in the optimization search we can assume possible wider ranges as physiologically acceptable because the optimized system would correspond with altered kinetic behavior with the use of drugs or other therapeutic agents.
Parameter estimation using 1000 iterations which required a computation time of 96 hours. The objective function value finally reached was 0.0517, with a maximal absolute error between interpolated data and model of 0.5828. In accordance with the objective function definition, this measures the average distance between the experimental data and the model. Since the standardized experimental variation values range between 0 and 2.5, an objective function value of 0.0517 means an average relative error lower than 5%, which is less than the experimental error and therefore enough to ensure that the model represents the experimental data.
All variables in the model are normalized. Parasite load is normalized with respect to its proper mean and the remaining variables are normalized with respect to the control group of mice. This standardization reduces the range of variation of the parameters and computation time and also exploits various properties of the GMA-PL formalism concerning the behavior of variables and parameters. Initial values of parasite load between 10 -10 and 0.01 were implemented, comparing them in terms of the objective function. The smaller the minimum of the objective function, the better the model (with the respective initial value) fits the experimental data. Of the initial values for which simulations meet the criterion f obj < 0.2, we found that the best fit and lowest objective function was attained with 10 -6 as the initial value for parasite load in the group infected with 10 3 parasites. Therefore 10 -6 was used as the initial value in that group.

Model selection strategies
The topology of the model is assumed to be as it appears in Figure 1. There are different models according to the number of parameters used; the most general model has 22 parameters. In other simplified models the value of certain parameters were assumed to be zero or one [21]. Among the different strategies of model selection for finding the one that yielded the "best" fit, we finally chose to specify a fixed value of the objective function and, accordingly, selected models that fell below this value. Out of these, we chose the simplest model, that is, the model with the least number of parameters.

Dynamic Sensitivities
Sensitivity analysis is a tool useful for model robustness evaluation and system dynamics characterization. Since our model studies the system dynamics, this tool enables us to identify the parameters with major influence on the transient dynamics. We have used the System Parameter Dynamic Sensitivity, S pk Xi , defined in equation 6: In the above expression, Xi is the considered variable, while pk is the parameter under scrutiny. Ak,2 is the area below the solution curve after a change of pk to pk,2. Ak is the area under the solution curve using the original value pk. Dk is the area between the two solution curves, using pk and pk,2, respectively. In our analysis we have considered changes in the two kinds of parameters, kinetic orders and rate constants. The S pk Xi value corresponds to the variation of the area under the variable time course after perturbation in parameter space. For each model variable the absolute values of the area A k between the original curve and the abscissa, the area A k2 between the new curve and the abscissa and the area D k between the original and the new curve are calculated using the trapezoidal method. Positive sensitivity means that the area under the "new" curve is greater than the area under the original curve, i.e., that the parameter pk has a positive influence on variable Xi. Negative sensitivity means the opposite. Zero sensitivity means that small changes in the parameter have no influence on the variable. All sensitivities were computed for standardized variables.