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

doi:10.1186/1752-0509-6-1

BMC Syst Biol. 2012; 6: 1.

Published online 2012 January 5. doi: 10.1186/1752-0509-6-1

PMCID: PMC3293051

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

Bettina M Länger,¹ Cristina Pou-Barreto,² Carlos González-Alcón,³ Basilio Valladares,² Bettina Wimmer,¹ and Néstor V Torres¹

¹Grupo de Tecnología Bioquímica. Departamento de Bioquímica y Biología Molecular. Universidad de La Laguna. 38206. San Cristóbal de La Laguna. Tenerife. Spain

²Instituto Universitario de Enfermedades Tropicales y Salud Pública de Canarias. Universidad de La Laguna. 38206. San Cristóbal de La Laguna. Tenerife. Spain

³Grupo de Tecnología Bioquímica. Departamento de Estadística, I.O. y Ciencias de la Computación. Universidad de La Laguna. 38206. San Cristóbal de La Laguna. Tenerife. Spain

Corresponding author.

Bettina M Länger: bettina.laenger/at/gmx.at; Cristina Pou-Barreto: cripou/at/ull.es; Carlos González-Alcón: cgalcon/at/ull.es; Basilio Valladares: bvallada/at/ull.es; Bettina Wimmer: bettina.wimmer/at/student.tuwien.ac.at; Néstor V Torres: ntorres/at/ull.edu.es

Received August 27, 2011; Accepted January 5, 2012.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

This article has been cited by other articles in PMC.

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 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-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-10]. The dynamics of parasite-host interactions in the infection process has also been studied using agent-based modeling approaches [11-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,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, host-pathogen 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 Figure11 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:

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 (more ...)

(1)

In the above expression, X_i, σ_ij, γ_jand g_jkrepresent the normalized variable set, the stoichiometric matrix, the rate constants, and the kinetic orders, respectively. The variables lymphocytes proliferation (X₂), IgG1 (X₃) and IgG2a (X₄) 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 σ_ij, γ_jand g_jkare determined using prior biological knowledge, information about the basal steady-states of the system [21], and/or dynamical experimental data [22,23]. In power-law 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 Figure11 shows the model scheme of the chosen variables and the influences among them denoted by arrows and parameters. g₁, ... , g₁₄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₁) stimulates its immune response. The parasites multiply in macrophages by binary division. The parasite load growth (V₁) has a nonlinear dependence on the parasite load through the kinetic order g₁. Increased parasite load leads to a decrease in the proliferation rate of lymphocytes (V₄); this interaction is represented by the kinetic order g₇[26]. Proliferation (or multiplication) of T lymphocytes (X₂) occurs when naïve T cells are activated by antigens of the pathogen (g₅) 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₃is also stimulated by X₂, through g₆. Cell mediated effectors enhance X₁decay (V₂); this effect is represented by the positive kinetic order g₃[28,29].

The host immune system produces IgG1 (X₃) and IgG2a (X₄) 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₂on the rate synthesis of IgG1 (V₅) through g₉, and on the rate synthesis of IgG2a (V₇) through g₁₂. These two immunoglobulins are antagonistic, so each of them has a negative influence on the generation rate of the other, namely X₄on V₅and X₃on V₇. These effects are represented by the kinetic orders g₁₀and g₁₃, respectively [32,33]. The IgG2a influences macrophage activity by stimulating the X₁rate decay, V₂. This interaction is represented in our model by the positive kinetic order g₄. It is assumed that the transformation rates V₂, V₄, V₅and V₈are proportional to X₁, X₂, X₃and X₄. These dependences are represented in the model by the positive kinetic orders g₂, g₈, g₁₁and g₁₄, 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 Figure22 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 (more ...)

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 Figure22 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⁶. 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 Figure33 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₁₀; X₂)) and 1.83 (S(γ₁; X₁)). The lower value corresponds to the influence of g₁₀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_jp_k).

In general, the higher sensitivities correspond to the variable indicated by the arrow, except for the parameters directly influencing lymphocyte proliferation (γ₃, γ₄, g₅, g₆, g₇and g₈). For the parasite load, sensitivities with absolute values greater than 1 are S_γ1^X1and 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^X1is 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 (γ₁, γ₂, g₁, g₂, g₃, and g₄) 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_iand rate constants, γ_i) that yields a reduced, minimum value of the parasite load, both during the time of infection evolution and at the final, 24-week stage.

The search program was organized in two phases. In the first phase, we changed only one of the parameters at a time (g_ior γ_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_iparameter scanning

After a systematic search among all kinetic constants, we found that g₁and g₆were the most suitable parameters to be changed for reducing parasite load. g₁, 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₁value from 0.53 to 3. Figure Figure4A4A shows the progression of parasite load during a time period of 24 weeks for values of g₁ranging from 0.01 to 3. Our model suggests that increases in g₁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₁< 0.5, final parasite load is proportional to g₁. However, for values of g₁> 0.5, final parasite load becomes inversely proportional to the value of this parameter. This latter region includes the actual g₁value. This implies that therapeutic strategies should aim to increase g₁. 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₁. A. Final parasite load during a time period of 24 weeks for different values of g₁. The dot indicates the position of the original system (g₁= 0.53). B. Parasite load over time for (more ...)

Figure Figure4B4B 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₁with respect to the parasite load (X₁), S_g1^X1, is negative (see Figure Figure2),2), therefore an increase in g₁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₁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₁in our model, which is influenced by g₁) 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₁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_iparameter with similarly minimizing effects on the parasite load is g₆. g₆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₆from its original value of 0.02 to 1.02 (Figure (Figure5).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₆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₆. A. Final parasite load during a time period of 24 weeks for different values of g₆. B. Comparison of the predicted evolution of the final parasite load over time for an optimized (more ...)

γ_iparameter scanning

We also carried out a systematic search among all rate constants. We observed that changes in all γ_iparameters reduce the final parasite load (see Figure Figure6).6). The rate constant that yields the minimum for final, maximal and mean parasite load is γ₂(after slightly increasing its value), which is the rate constant for parasite degradation (see Figure Figure1).1). Figure Figure6A6A shows final parasite load for different values of γ₂. Figure Figure6B6B shows the effect on the parasite load for two increased values of γ₂: 0.43 and 4.3.

Figure 6

Evolution of the parasite load over time for different values of γ₂. A. Final parasite load during a time period of 24 weeks for different values of γ₂. B. Parasite load over time for (γ₂)_original(0.043, dotted line); 0.43 (discontinuous (more ...)

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₁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 (Table1).1). The other three variables (IgG1, IgG2a, lymphocytes) remain almost unchanged (results not shown). Figure Figure77 shows all possible combinations of two parameters and the optimum predicted final parasite load. Aside from the best combinations (see Figure Figure7,7, black column and file) g₁and g₂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 Table11 show the parameter change factor for the two-parameter combinations involving g₁(black column and file in Figure Figure7).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 γ₁, γ₄, γ₆, g₂, g₃, g₅, and g₁₄. Their values are reduced by factors ranging from 5 to 85% (see Table Table11).

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₂) or inhibit proliferation (g₁), whereas immunotherapeutic treatment implies changing parameters g₃, ... g₈and g₃, g₄, g₆, g₈, g₉, ...g₁₂. 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₁, which describes the influence of parasites on their own proliferation; g₆, which represents the influence of lymphocytes on their own proliferation; and γ₂, 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 γ₂, 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 IL-12 (stimulates Th1) also increase γ₂. 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₁. The discrepancy in our model's predictions can be explained by several facts. First, in all cases where a decrease in g₁could be assumed, there is also the concomitant effect of increasing γ₂, 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₁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 (γ₁) and reduces parasite-toxic production of nitric oxide (γ₂).

Furthermore, no existing drug is known to have an effect on g₆, 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 Figure7,7, combinations of drugs able to increase g₁or g₂and simultaneously change any other parameter, or, alternatively, combinations of drugs that decrease g₁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₂(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₃), rLmSTI1 (increase in g₄), and chemokines (that increase g₃and g₂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 parasitic-related 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³(10⁶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^-6were prepared twice for each homogenate in wells of 96-well 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.

(3)

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₄Cl, 10 mM KHCO₃, 1 mM EDTA, pH 7.4) and remaining cells were finally resuspended to a density of 2.5 × 10⁶cells/ml in DMEM containing 10% FCS, 2 mM L-glutamine, 0.05 mM 2-mercaptoethanol, 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⁵cells) in 96-well plates and they were allowed to proliferate for 3 days at 37°C in an atmosphere containing 5% CO₂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:

(4)

where

Abs_Twand equation M6

Abs_Cwstand 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:

(5)

where n_varis the number of variables monitored and n_tpthe number of time points when each variable was measured. In turn, X_j(t_i) is the predicted value for the j^thvariable at the i^thtime point obtained after numerical integration of the solution, while X_j^exp(t_i) is the value of the j^thobservable variable at the i^thtime 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 γ_iwere permitted to vary within the range [0,10] and the g_iwithin the range [0,3], except the parameters g₂, g₈, g₁₁and g₁₄which were set to 1 since they relate variables with their own outflow (V₂, V₄, V₆and V₈, 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 γ_ivalues, 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^-10and 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^-6as the initial value for parasite load in the group infected with 10³parasites. Therefore 10^-6was 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 Figure1.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:

(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^Xivalue 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_kbetween the original curve and the abscissa, the area A_k2between the new curve and the abscissa and the area D_kbetween 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 power-law 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. CP-B carried out parasite, antigen, antibodies and lymphoproliferation quantization. BML implemented the GMA-PL model.

BML, BW and CG-A performed the optimum parameter searches. NVT, BV, BML and CP-B and CG-A 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ález-Alcón and Guido Santos and for their helpful observations. This work was funded by research Grants from Spanish MICINN (Ref. BIO2011-29233-C02-02, 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:563–570. doi: 10.3109/00365541003789473. [PubMed] [Cross Ref]
Killick-Kendrick R. The biology and control of phlebotomine sand flies. Clin Dermatol. 1999;17:279–289. doi: 10.1016/S0738-081X(99)00046-2. [PubMed] [Cross Ref]
World Health Organization. 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:229–243. doi: 10.1128/CMR.14.2.229-243.2001. [PMC free article] [PubMed] [Cross Ref]
Chaves LF, Hernandez MJ. Mathematical modelling of American cutaneous leishmaniasis: incidental hosts and threshold conditions for infection persistence. Acta Trop. 2004;92:245–252. doi: 10.1016/j.actatropica.2004.08.004. [PubMed] [Cross Ref]
Bacaër N, Guernaoui S. The epidemic threshold of vector-borne diseases with seasonality: the case of cutaneous leishmaniasis in Chichaoua, Morocco. J Math Biol. 2006;53:421–436. doi: 10.1007/s00285-006-0015-0. [PubMed] [Cross Ref]
Carneiro DD, Bavia ME, Rocha WJ, Tavares AC, Cardim LL, Alemayehu B. Application of spatio-temporal scan statistics for the detection of areas with increased risk for American visceral leishmaniasis in the state of Bahia, Brazil. Geospat Health. 2007;2:113–126. [PubMed]
Leifso K, Cohen-Freue 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:35–46. doi: 10.1016/j.molbiopara.2006.11.009. [PubMed] [Cross Ref]
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. [PMC free article] [PubMed]
Rochette A, Raymond F, Ubeda JM, Smith M, Messier N, Boisvert S, Rigault P, Corbeil J, Ouellette M, Papadopoulou B. Genome-wide gene expression profiling analysis of Leishmania major and Leishmania infantum developmental stages reveals substantial differences between the two species. BMC Genomics. 2008;9:255. doi: 10.1186/1471-2164-9-255. [PMC free article] [PubMed] [Cross Ref]
de Almeida MC, Moreira HN. A mathematical model of immune response in cutaneous leishmaniasis. Journal of Biological Systems. 2007;3:313–354.
Schweitzer N, Swinton J, Anderson R. Dynamic interaction between Leishmania infection in mice and Th1-type CD4⁺T-cells: complexity in outcome without a requirement for Th2-type responses. Parasite Immunol. 1993;15:85–99. doi: 10.1111/j.1365-3024.1993.tb00588.x. [PubMed] [Cross Ref]
Dancik GM, Jones DE, Dorman KS. Parameter estimation and sensitivity analysis in an agent-based model of Leishmania major infection. J Theor Biol. 2010;262:398–412. doi: 10.1016/j.jtbi.2009.10.007. [PMC free article] [PubMed] [Cross Ref]
Dantas-Torres F, Brandão-Filho SP. Visceral leishmaniasis in Brazil: revisiting paradigms of epidemiology and control. Rev Inst Med Trop Sao Paulo. 2006;48:151–156. doi: 10.1590/S0036-46652006000300007. [PubMed] [Cross Ref]
Palatnik-de-Sousa CB, Silva-Antunes 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:3505–3512. doi: 10.1016/j.vaccine.2009.03.045. [PubMed] [Cross Ref]
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:1373–1383. doi: 10.1016/S0020-7519(03)00179-6. [PubMed] [Cross Ref]
Arrais-Silva WW, Paffaro VA Jr, Yamada AT, Giorgio S. Expression of hypoxia-inducible factor-1alpha in the cutaneous lesions of BALB/c mice infected with Leishmania amazonensis. Exp Mol Pathol. 2005;78:49–54. doi: 10.1016/j.yexmp.2004.09.002. [PubMed] [Cross Ref]
Lira R, Doherty M, Modi G, Sacks D. Evolution of lesion formation, parasitic load, immune response, and reservoir potential in C57BL/6 mice following high-and low-dose challenge with Leishmania major. Infect Immun. 2000;68:5176–5182. doi: 10.1128/IAI.68.9.5176-5182.2000. [PMC free article] [PubMed] [Cross Ref]
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:269–281. doi: 10.1016/S0165-2427(00)00221-X. [PubMed] [Cross Ref]
Dea-Ayuela MA, Rama-Iñiguez S, Alunda JM, Bolás-Fernandez F. Setting new immunobiological parameters in the hamster model of visceral leishmaniasis for in vivo testing of antileishmanial compounds. Vet Res Commun. 2007;31:703–717. doi: 10.1007/s11259-007-0040-5. [PubMed] [Cross Ref]
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 JAK2-STAT5 signalling pathway. BMC Syst Biol. 2008;2:38. doi: 10.1186/1752-0509-2-38. [PMC free article] [PubMed] [Cross Ref]
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:2281–2289. doi: 10.1093/bioinformatics/btm326. [PubMed] [Cross Ref]
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, High-Biomass Density Culture. Biotechnol Bioeng. 2009;104:785–795. [PubMed]
Danø S, Madsen MF, Schmidt H, Cedersund G. Reduction of a biochemical model with preservation of its basic dynamic properties. FEBS J. 2006;273:4862–4877. doi: 10.1111/j.1742-4658.2006.05485.x. [PubMed] [Cross Ref]
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, Vaigro-Wolff A. Feasibility of a high-flux anticancer drug screen using a diverse panel of cultured human tumor cell lines. J Natl Cancer Inst. 1991;83:757–766. doi: 10.1093/jnci/83.11.757. [PubMed] [Cross Ref]
Chang KP, McGwire BS. Molecular determinants and regulation of Leishmania virulence. Kinetoplastid Biol Dis. 2002;1:1–7. doi: 10.1186/1475-9292-1-1. [PMC free article] [PubMed] [Cross Ref]
Chang KP, Reed SG, McGwire BS, Soong L. Leishmania model for microbial virulence: the relevance of parasite multiplication and pathoantigenicity. Acta Trop. 2003;85:375–390. doi: 10.1016/S0001-706X(02)00238-3. [PubMed] [Cross Ref]
Mougneau E, Bihl F, Glaichenhaus N. Cell biology and immunology of Leishmania. Immunol Rev. 2011;240:286–296. doi: 10.1111/j.1600-065X.2010.00983.x. [PubMed] [Cross Ref]
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:1413–1419. doi: 10.1007/s00436-008-1150-x. [PubMed] [Cross Ref]
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:1675–1684. doi: 10.1084/jem.168.5.1675. [PMC free article] [PubMed] [Cross Ref]
Duncan DD, Swain SL. Role of antigen-presenting 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:2506–2514. doi: 10.1002/eji.1830241037. [PubMed] [Cross Ref]
van Griensven J, Balasegaram M, Meheus F, Alvar J, Lynen L, Boelaert M. Combination therapy for visceral leishmaniasis. Lancet Infect Dis. 2010;10:184–194. doi: 10.1016/S1473-3099(10)70011-6. [PubMed] [Cross Ref]
Segovia-Juarez JL, Ganguli S, Kirschner D. Identifying control mechanisms of granuloma formation during M. tuberculosis infection using an agent-based model. J Theor Biol. 2004;231:357–376. doi: 10.1016/j.jtbi.2004.06.031. [PubMed] [Cross Ref]
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:2167–2168. [PMC free article] [PubMed]
Larreta R, Guzman F, Patarroyo ME, Alonso C, Requena JM. Antigenic properties of the Leishmania infantum GRP94 and mapping of linear B-cell epitopes. Immunol Lett. 2002;80:199–205. doi: 10.1016/S0165-2478(01)00331-5. [PubMed] [Cross Ref]
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:50–58. doi: 10.1016/j.jbiotec.2008.07.1814. [PubMed] [Cross Ref]
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:803–812. doi: 10.1038/nrm2042. [PubMed] [Cross Ref]

Articles from BMC Systems Biology are provided here courtesy of
BioMed Central