Abstract
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 powerlaw 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 parasitehost 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 [57]. In addition, models have also been used to study regulation of gene expression, protein synthesis, and metabolism of the parasite at the genomewide level [810]. The dynamics of parasitehost interactions in the infection process has also been studied using agentbased modeling approaches [1113]. For example, Dancik et al. [13] used an agentbased 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,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 ArraisSilva et al. [17] on the hypoxiainducible factor1 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 DeaAyuela 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.
Results
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 powerlaw formulation, the General Mass Action Power Law formalism (GMAPL) that allows for noninteger kinetic orders [21,22] with the following structure:
Figure 1. Leishmaniasis progression model. Solid arrows represent synthesis (input arrows) or degradation (output arrows) fluxes (each flux number notated by its corresponding γ_{i}); dashed arrows are the signals among processes variables which are quantified by the corresponding g_{i }value. Positive and negative signs denote activation and inhibition of the corresponding fluxes respectively. See text for the nomenclature.
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 parasitefree, 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 GMAPL formalism on the behavior of variables and parameters. The specific numerical values for the parameters σ_{ij}, γ_{j }and g_{jk }are determined using prior biological knowledge, information about the basal steadystates of the system [21], and/or dynamical experimental data [22,23]. In powerlaw models, kinetic orders can have noninteger values. One of the main advantages of powerlaw 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 saturatinglike manner. Values equal to one imply a flux that depends linearly on the variable or, in chemical terms, a first order reaction. By allowing noninteger, 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:
Figure 2 shows the model data fitting, displaying a good correlation between the experimental and estimated data.
Figure 2. Data fit and predicted model dynamics of the four model variables. The panels under Model Data Fitting shows the data fit for the time series data of the four model variables. Panels under Predicted Model dynamics show the comparison of predicted and measured system variable dynamics for an initial parasite load of 10^{6}. The continuous lines indicate the estimated dynamics while the dotted lines indicate the experimental data. Error bars indicate the standard deviation among mice.
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.
Sensitivity Analysis
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.
Figure 3. Absolute value of the system dynamic sensitivities S_{pk}^{Xi }(X_{j }p_{k}).
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 the set of parameter values (kinetic constants, g_{i }and rate constants, γ_{i}) that yields a reduced, minimum value of the parasite load, both during the time of infection evolution and at the final, 24week stage.
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 upperbound 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.
Singleparameter 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 4. Evolution of the parasite load over time for different values of g_{1}. A. Final parasite load during a time period of 24 weeks for different values of g_{1}. The dot indicates the position of the original system (g_{1 }= 0.53). B. Parasite load over time for g_{1 }= 3.
Figure 4B shows that the evolution of the parasite load decreases right from the beginning, almost linearly through the 24week 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 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 1720 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.
Figure 5. Evolution of the final parasite load over time for different values of g_{6}. A. Final parasite load during a time period of 24 weeks for different values of g_{6}. B. Comparison of the predicted evolution of the final parasite load over time for an optimized value of g_{6 }(1.02, line) with the experimentally observed parasite load values (dots). In the experimental setup the estimated value for g_{6 }is 0.022.
γ_{i }parameter scanning
We also carried out a systematic search among all rate constants. We observed that changes in all γ_{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 γ_{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 γ_{2}. Figure 6B shows the effect on the parasite load for two increased values of γ_{2}: 0.43 and 4.3.
Figure 6. Evolution of the parasite load over time for different values of γ_{2}. A. Final parasite load during a time period of 24 weeks for different values of γ_{2}. B. Parasite load over time for (γ_{2})_{original }(0.043, dotted line); 0.43 (discontinuous line); 4.3 (continuous line) and experimental data (dots).
Combined twoparameter 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). The other three variables (IgG1, IgG2a, lymphocytes) remain almost unchanged (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.
Table 1. Parameter change factors for the two parameter combinations involving g_{1.}
Figure 7. Optimized final parasite load obtained for each possible combination of two parameters of the system. The parasite load of the fitted system was 5.5.
By way of illustration, Table 1 show the parameter change factor for the twoparameter 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 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 IL12 (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 tradeoff 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. Insulinelike 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 antileishmanial effectiveness. Insulinelike growth factor also increases the number of parasites (γ_{1}) and reduces parasitetoxic 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: interleukin5,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, modelbased 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, 68 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% heatinactivated 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 28week 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% heatinactivated 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 enzymelinked 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 HRPOconjugated goat antimouse 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 Lglutamine, 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; BIORAD). Stimulation index of the lymphocyte proliferation (SI) of experimental mice was calculated according the following equation:
where Abs_{Tw }and 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 powerlaw 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 powerlaw 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 24week 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 functiondefined 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 GMAPL 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.
Authors' contributions
BV and NVT proposed the elaboration of a mathematical model based on the GMA powerlaw formalism of the progression of the disease and with the potential to be used for the design of effective therapies. They planned the work and coordinated the study. CPB carried out parasite, antigen, antibodies and lymphoproliferation quantization. BML implemented the GMAPL model.
BML, BW and CGA performed the optimum parameter searches. NVT, BV, BML and CPB and CGA defined the model and obtained the numerical parameters used in the paper. All authors read and approved the final manuscript.
Acknowledgements
The authors gratefully acknowledge Dr Carlos GonzálezAlcón and Guido Santos and for their helpful observations. This work was funded by research Grants from Spanish MICINN (Ref. BIO201129233C0202, PI081984 and RD06/0021/0005) and from Agencia Canaria de Investigación, Innovación y Sociedad de la Información (Ref. Project PIL2070901 and PIL2071001).
References

Neghina R, Neghina AM: Leishmaniasis, a global concern for travel medicine.
Scand J Infect Dis 2010, 42:563570. PubMed Abstract  Publisher Full Text

KillickKendrick R: The biology and control of phlebotomine sand flies.
Clin Dermatol 1999, 17:279289. PubMed Abstract  Publisher Full Text

Control of the leishmaniasis: report of a meeting of the WHO Expert Commitee on the Control of Leishmaniases. Geneva. 2010.

Handman E: Leishmaniasis: current status of vaccine development.
Clin Microbiol Rev 2001, 14:229243. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Chaves LF, Hernandez MJ: Mathematical modelling of American cutaneous leishmaniasis: incidental hosts and threshold conditions for infection persistence.
Acta Trop 2004, 92:245252. PubMed Abstract  Publisher Full Text

Bacaër N, Guernaoui S: The epidemic threshold of vectorborne diseases with seasonality: the case of cutaneous leishmaniasis in Chichaoua, Morocco.
J Math Biol 2006, 53:421436. PubMed Abstract  Publisher Full Text

Carneiro DD, Bavia ME, Rocha WJ, Tavares AC, Cardim LL, Alemayehu B: Application of spatiotemporal scan statistics for the detection of areas with increased risk for American visceral leishmaniasis in the state of Bahia, Brazil.
Geospat Health 2007, 2:113126. PubMed Abstract  Publisher Full Text

Leifso K, CohenFreue G, Dogra N, Murray A, McMaster WR: Genomic and proteomic expression analysis of Leishmania promastigote and amastigote life stages: the Leishmania genome is constitutively expressed.
Mol Biochem Parasitol 2007, 152:3546. PubMed Abstract  Publisher Full Text

Chavali AK, Whittemore JD, Eddy JA, Williams KT, Papin JA: Systems analysis of metabolism in the pathogenic trypanosomatid Leishmania major.
Mol Syst Biol 2008, 4:177. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Rochette A, Raymond F, Ubeda JM, Smith M, Messier N, Boisvert S, Rigault P, Corbeil J, Ouellette M, Papadopoulou B: Genomewide gene expression profiling analysis of Leishmania major and Leishmania infantum developmental stages reveals substantial differences between the two species.
BMC Genomics 2008, 9:255. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

de Almeida MC, Moreira HN: A mathematical model of immune response in cutaneous leishmaniasis.

Schweitzer N, Swinton J, Anderson R: Dynamic interaction between Leishmania infection in mice and Th1type CD4^{+ }Tcells: complexity in outcome without a requirement for Th2type responses.
Parasite Immunol 1993, 15:8599. PubMed Abstract  Publisher Full Text

Dancik GM, Jones DE, Dorman KS: Parameter estimation and sensitivity analysis in an agentbased model of Leishmania major infection.
J Theor Biol 2010, 262:398412. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

DantasTorres F, BrandãoFilho SP: Visceral leishmaniasis in Brazil: revisiting paradigms of epidemiology and control.
Rev Inst Med Trop Sao Paulo 2006, 48:151156. PubMed Abstract  Publisher Full Text

PalatnikdeSousa CB, SilvaAntunes I, Morgado Ade A, Menz I, Palatnik M, Lavor C: Decrease of the incidence of human and canine visceral leishmaniasis after dog vaccination with Leishmune^{® }in Brazilian endemic areas.
Vaccine 2009, 27:35053512. PubMed Abstract  Publisher Full Text

Courret N, Lang T, Milon G, Antoine JC: Intradermal inoculations of low doses of Leishmania major and Leishmania amazonensis metacyclic promastigotes induce different immunoparasitic processes and status of protection in BALB/c mice.
Int J Parasitol 2003, 33:13731383. PubMed Abstract  Publisher Full Text

ArraisSilva WW, Paffaro VA Jr, Yamada AT, Giorgio S: Expression of hypoxiainducible factor1alpha in the cutaneous lesions of BALB/c mice infected with Leishmania amazonensis.
Exp Mol Pathol 2005, 78:4954. PubMed Abstract  Publisher Full Text

Lira R, Doherty M, Modi G, Sacks D: Evolution of lesion formation, parasitic load, immune response, and reservoir potential in C57BL/6 mice following highand lowdose challenge with Leishmania major.
Infect Immun 2000, 68:51765182. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Requena JM, Soto M, Doria MD, Alonso C: Immune and clinical parameters associated with Leishmania infantum infection in the golden hamster model.
Vet Immunol Immunopathol 2000, 76:269281. PubMed Abstract  Publisher Full Text

DeaAyuela MA, RamaIñiguez S, Alunda JM, BolásFernandez F: Setting new immunobiological parameters in the hamster model of visceral leishmaniasis for in vivo testing of antileishmanial compounds.
Vet Res Commun 2007, 31:703717. PubMed Abstract  Publisher Full Text

Voit EO: Computational Analysis of Biochemical Sustems. A Practical Guide for Biochemists and Molecular Biologist. Cambridge University Press; 2000.

Vera J, Bachmann J, Pfeifer AC, Becker V, Hormiga JA, Darias NV, Timmer J, Klingmüller U, Wolkenhauer O: A systems biology approach to analyse amplification in the JAK2STAT5 signalling pathway.
BMC Syst Biol 2008, 2:38. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Vera J, Curto R, Cascante M, Torres NV: Detection of potential enzyme targets by metabolic modelling and optimization. Application to a simple enzymopathy.
Bioinformatics 2007, 23:22812289. PubMed Abstract  Publisher Full Text

Guebel DV, Cánovas M, Torres NV: Model Identification in Presence of Incomplete Information by Generalized Principal Component Analysis: Application to the Common and Differential Responses of Escherichia coli to Multiple Pulse Perturbations in Continuous, HighBiomass Density Culture.
Biotechnol Bioeng 2009, 104:785795. PubMed Abstract  Publisher Full Text

Danø S, Madsen MF, Schmidt H, Cedersund G: Reduction of a biochemical model with preservation of its basic dynamic properties.
FEBS J 2006, 273:48624877. PubMed Abstract  Publisher Full Text

Abbas AK, Lichtman AH, Pillai S: Cellular and molecular immunology. Elsevier Press; 2010.

Monks A, Scudiero D, Skehan P, Shoemaker R, Paull K, Vistica D, Hose C, Langley J, Cronise P, VaigroWolff A: Feasibility of a highflux anticancer drug screen using a diverse panel of cultured human tumor cell lines.
J Natl Cancer Inst 1991, 83:757766. PubMed Abstract  Publisher Full Text

Chang KP, McGwire BS: Molecular determinants and regulation of Leishmania virulence.
Kinetoplastid Biol Dis 2002, 1:17. PubMed Abstract  BioMed Central Full Text  PubMed Central Full Text

Chang KP, Reed SG, McGwire BS, Soong L: Leishmania model for microbial virulence: the relevance of parasite multiplication and pathoantigenicity.
Acta Trop 2003, 85:375390. PubMed Abstract  Publisher Full Text

Mougneau E, Bihl F, Glaichenhaus N: Cell biology and immunology of Leishmania.
Immunol Rev 2011, 240:286296. PubMed Abstract  Publisher Full Text

Kaur S, Kaur T, Garg N, Mukherjee S, Raina P, Athokpam V: Effect of dose and route of inoculation on the generation of CD4+ Th1/Th2 type of immune response in murine visceral leishmaniasis.
Parasitol Res 2008, 103:14131419. PubMed Abstract  Publisher Full Text

Scott P, Natovitz P, Coffman RL, Pearce E, Sher A: Immunoregulation of cutaneous leishmaniasis. T cell lines that transfer protective immunity or exacerbation belong to different T helper subsets and respond to distinct parasite antigens.
J Exp Med 1988, 168:16751684. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Duncan DD, Swain SL: Role of antigenpresenting cells in the polarized development of helper T cell subsets: evidence for differential cytokine production by Th0 cells in response to antigen presentation by B cells and macrophages.
Eur J Immunol 1994, 24:25062514. PubMed Abstract  Publisher Full Text

van Griensven J, Balasegaram M, Meheus F, Alvar J, Lynen L, Boelaert M: Combination therapy for visceral leishmaniasis.
Lancet Infect Dis 2010, 10:184194. PubMed Abstract  Publisher Full Text

SegoviaJuarez JL, Ganguli S, Kirschner D: Identifying control mechanisms of granuloma formation during M. tuberculosis infection using an agentbased model.
J Theor Biol 2004, 231:357376. PubMed Abstract  Publisher Full Text

Buffet PA, Sulahian A, Garin YJ, Nassar N, Derouin F: Culture microtitration: a sensitive method for quantifying Leishmania infantum in tissues of infected mice.
Antimicrob Agents Chemother 1995, 39:21672168. PubMed Abstract  Publisher Full Text  PubMed Central Full Text

Larreta R, Guzman F, Patarroyo ME, Alonso C, Requena JM: Antigenic properties of the Leishmania infantum GRP94 and mapping of linear Bcell epitopes.
Immunol Lett 2002, 80:199205. PubMed Abstract  Publisher Full Text

Hormiga JA, Vera J, Frías I, Torres Darias NV: Growth and ligninolytic system production dynamics of the Phanerochaete chrysosporium fungus. A modelling and optimization approach.
J Biotechnol 2008, 137:5058. PubMed Abstract  Publisher Full Text

Albeck JG, MacBeath G, White FM, Sorger PK, Lauffenburger DA, Gaudet S: Collecting and organizing systematic sets of protein data.
Nat Rev Mol Cell Biol 2006, 7:803812. PubMed Abstract  Publisher Full Text