Treatment of leishmaniasis with chemotherapy and vaccine: a mathematical model

Abstract

Leishmaniasis, an infectious disease, manifests itself mostly in two forms, cutaneous leishmaniasis (CL) and, a more severe and potentially deadly form, visceral leishmaniasis (VL). The current control strategy for leishmaniasis relies on chemotherapy drugs such as sodium antimony gluconate (SAG) and meglumine antimoniate (MA). However, all these chemotherapy compounds have poor efficacy, and they are associated with toxicity and other adverse effects, as well as drug resistance. While research in vaccine development for leishmaniasis is continuously progressing, no vaccine is currently available. However, some experimental vaccines such as LEISH-F1+MPL-SE (V) have demonstrated some efficacy when used as drugs for CL patients. In this paper we use a mathematical model to address the following question: To what extent vaccine shots can enhance the efficacy of standard chemotherapy treatment of leishmaniasis? Starting with standard MA treatment of leishmaniasis and combining it with three injections of V , we find, by Day 84, that efficacy increased from 29% to 65–91% depending on the amount of the vaccine. With two or just one injection of V , efficacy is still very high, but there is a definite resurgence of the disease by end-time.

Keywords:

• Leishmaniasis
• vaccine LEISH-F1+
• MPL-SE
• sodium stibugluconate
• sodium antomony gluconate
• meglumine antimoniate
• immune response

1. Introduction

Leishmaniases are a group of diseases caused by protozoan obligate intracellular parasites from more than 20 Leishmania species. The parasites are transmitted to humans by a bite of an infected female sand fly, which injects into the body the flagellated form of the parasite, the promastigote. Promastigotes are endocytozed into phagocytic cells (macrophages) and quickly transform into amastigotes (non-flagellated Leishmania). The amastigotes then mature and multiply, and eventually cause the infected cells to burst.

The main two forms of the disease are cutaneous leishmaniasis (CL) and visceral leishmaniasis (VL). CL is the more common form of the disease, while VL is the more severe form of the disease. If not treated, VL is expected to cause death in more than 90% of the cases. Incidence of 0.7–1 million new cases of leishmaniasis occur annually, while 310 million people are vulnerable to infection.

The immune response to leishmaniasis is taken by CD4${}^{+}$ and CD8${}^{+}$ T cells. The T cells produce inflammatory cytokines that kill the intracellular parasites [53]. But in most cases, the immune system is unable to resolve the disease on its own. Currently, the control strategy for leishmaniasis relies on pharmacological treatment. But only a handful of chemotherapeutic agents are in use, and each of the options has several limitations, such as adverse effects (e.g. cardiotoxicity), poor efficacy, and drug resistance [73, 75, 86]. Over the last decades, only a few drugs have been approved, and currently, combinations of drugs are recommended by WHO.

Several vaccines have been developed to activate the immune system against the parasite in both CL and VL patients. Although work on vaccine development is continuously progressing [64], no vaccine is currently available to prevent any form of the disease.

In this paper, we focus on Leishmania vaccine LEISH-F1+MPL-SE, a recombinant protein through genetically engineered-cells. LEISH-F1 acts as an antigen, and MPL-SE is an adjuvant [90]. In clinical trials [25, 65] it was shown that as immunotherapy LEISH-F1+MPL-SE confers benefits in CL and VL patients, and the benefits increase when used in combination with meglumine antimoniate, one of the antimony-carbohydrate complexes.

A review article [23] suggests that chemotherapy along with immunotherapy can elicit a protective immune response and clear infection more effectively. More specifically, we consider the following question: To what extent can one of the experimental vaccines enhance the efficacy of chemotherapy treatment of leishmaniasis? In this paper, we address this question in a case where LEISH-F1+MPL-SE (as an immunotherapy drug) is combined with one of two pentavalent antimonial drugs: meglumine antimoniate (MA) and stibugluconate (also called sodium antimony gluconate (SAG)). Both drugs are first-line chemotherapy drugs for the treatment of leishmaniasis [28, 37, 86]. Using a mathematical model, we evaluate the efficacy of various combinations of the vaccine with one of the two drugs. The model is based on an earlier mathematical model by Siewe et al. [79] on the immune response to infection by Leishmania, which was extended in order to include the effects of the drugs.

2. Mathematical model

The model variables are listed in Table 1 in units of g/cm${}^3$. The immune system includes pro-inflammatory macrophages M1, anti-inflammatory macrophages M2, dendritic cells D, and T cells, which are taken to represent CD4${}^{+}$ Th1 cells and CD8${}^{+}$ T cells combined. Cytokines produced by the inflammatory immune cells are IL-2, IL-12, IFN-γ, TNF-α and nitric oxide (NO); IL-10 and TGF-β are produced by M2 macrophages. T cells are activated by dendritic cells presenting pieces of the parasites by MHC class I. T cells are also activated by IL-12, a process blocked by IL-10, and IL-2 increases the proliferation of T cells. IFN-γ is produced by T cells, TNF-α and NO are produced by macrophages, and macrophages kill the intracellular parasites. TNF-α mediates the polarization of M2 into M1 macrophages.

Table 1. List of variables of the model.

Variables	Descriptions	Variables	Descriptions
D	Density of dendritic cells	$M_1$	Density of pro-inflammatory
$M_2$	Density of anti-inflammatory macrophages	T	Density of combined CD4${}^{+}$ and CD8${}^{+}$ T cells
$P_1$	Density of Leishmania in M1 macrophages	$P_2$	Density of Leishmania in M2 macrophages
$I_2$	Concentration of IL-2	$I_{10}$	Concentration of IL-10; IL-10 is combined with IL-13
$I_{12}$	Concentration of IL-12	$I_{\gamma }$	Concentration of IFN-γ
NO	Concentration of Nitric Oxide	$T_{\alpha }$	Concentration of TNF-α
$T_{\beta }$	Concentration of TGF-β
A	Concentration of meglumine antimoniate (MA)	S	Concentration of sodium antimony gluconate (SAG)
V	Concentration of vaccine LEISH-F1+MPL-SE

Macrophage polarization plays a role in leishmaniasis [69]. TNF-α and IFN-γ induce polarization from M2 to M1, while IL-10 and TGF-β induce polarization from M1 to M2.

Meglumine antimoniate (A) enters macrophages and enhances their production of TNF-α and NO by M1 macrophages and IL-10 production by M2 macrophages [34]. SAG (S) induces macrophages to produce NO [37]; triggers MHC class I in dendritic cells, hence the production of IL-12 by the increased activities of M1 macrophages [7, 31, 37], and stimulates T cells without antigen, i.e. by the production of IL-2 by T cells [7, 37].

Leishmania vaccine LEISH-F1+MPL-SE (V) increases the production of IFN-γ by the CD4${}^{+}$ T cells (both in patients and in healthy individuals) [14, 16, 90], and thereby increases the activity of M1 macrophages in their production of TNF-α [7]; it also increases the IL-2 production by CD4${}^{+}$ T cells [16].

Figure 1 is a network of interactions among the model variables; the parasites within M1 and M2 macrophages are denoted by $P_1$ and $P_2$, respectively. The mathematical model is based on Figure 1, and it is represented by a system of ODEs. The equations are given below with full explanations and references.

Figure 1. Effects of vaccine LEISH-F1+MPL-SE (V), and drugs sodium antimony gluconate (SAG, S) and meglumine antimoniate (MA, A) on the immune response to Leishmania infection.

2.1. Equation for dendritic cells D

Dendritic cells are activated by ingesting parasites $P_1+P_2$. We denote the average number of parasites in a bursting dendritic cell, or a macrophage, by N, and assume, as in [79], that the mass of 150 parasites is equal to the mass of one macrophage. By the Michaelis-Menten law, the activation of D is proportional to $f(P_1+P_2)\equiv \frac{P_1+P_2}{(P_1+P_2)+K_P}$where we take the Michaelis-Menten parameter $K_P$ such that $f(P_1+P_2)=\frac{1}{2}f(\mathrm{\infty })$where $P_1+P_2$ is the average bursting pressure of the parasites; hence $K_P=\frac{N(M_1+M_2)}{150}.$The density of dendritic cells satisfies the equation: (1) $\frac{\mathrm{d}D}{\mathrm{d}t}=\underset{\mathrm{activation}}{\underbrace{{\lambda }_D{D}_0\frac{P_1+P_2}{(P_1+P_2)+K_P}}}-\underset{\mathrm{loss}\mathrm{due}\mathrm{to}\mathrm{death}}{\underbrace{{\mu }_{D}D}},$(1) where ${\lambda }_D$ is proportional to the density of inactive dendritic cells $D_0$. We assume, for simplicity, that the number of inactive dendritic cells remains constant during infection.

2.2. Equations for macrophages $M_1$ and $M_2$

In the sequel, we take the activation rate of cells by a cytokine X to be proportional to $X/(X+K_X)$, and, the resistance to activation of cells, to be proportional to $1/(1+X/K_X)$. We take $K_X$ to be the steady-state, or average, of X and refer to it as the half-saturation of X.

The dynamics of the density of $M_1$ is given by: (2) $\begin{array}{rl}\frac{\mathrm{d}{M}_1}{\mathrm{d}t}& =\underset{\mathrm{recruitment}}{\underbrace{{\lambda }_{M_1}}}+\underset{I_{\gamma },T_{\alpha }\text{-induced transition from }{M}_2\mathrm{to}{M}_1}{\underbrace{M_2({\lambda }_{M_2{M}_1{I}_{\gamma }}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{M_2{M}_1{T}_{\alpha }}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & -\underset{\text{loss due to transition to }{M}_2}{\underbrace{M_1({\lambda }_{M_1{M}_2{I}_{10}}\frac{I_{10}}{I_{10}+K_{I_{10}}}+{\lambda }_{M_1{M}_2{T}_{\beta }}\frac{T_{\beta }}{T_{\beta }+K_{T_{\beta }}})\frac{1}{1+I_{\gamma }/K_{I_{\gamma }}}}}\\ & -\underset{\text{loss due to bursting}}{\underbrace{{\tilde{\mu }}_{M_1}{M}_1\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}}}-\underset{\text{loss due to death}}{\underbrace{{\mu }_{M_1}{M}_1}}.\end{array}$(2) The first term on the right-hand side is the recruitment of new macrophages from the host immune system. The second term on the right-hand side is the gain in $M_1$ due to $I_{\gamma }$-induced and $T_{\alpha }$-induced transition from $M_2$ to $M_1$ [6, 17, 50, 94], respectively. The third term is the loss of $M_1$ due to the transition to $M_2$, a process stimulated by $I_{10}$ and $T_{\beta }$ [32], and inhibited by $I_{\gamma }$ [46, 59, 84]. The fourth term is the loss of $M_1$ due to bursting [9, 21, 38], modelled by the Hill dynamics, where $K=\frac{N{M}_1}{150}$ is the parasite pressure at the bursting of $M_1$.

The dynamics of the density of $M_2$ satisfies the equation: (3) $\begin{array}{rl}\frac{\mathrm{d}{M}_2}{\mathrm{d}t}& =\underset{\mathrm{recruitment}}{\underbrace{{\lambda }_{M_2}}}+\underset{\text{gain due to transition from}{M}_1}{\underbrace{M_1({\lambda }_{M_1{M}_2{I}_{10}}\frac{I_{10}}{I_{10}+K_{I_{10}}}+{\lambda }_{M_1{M}_2{T}_{\beta }}\frac{T_{\beta }}{T_{\beta }+K_{T_{\beta }}})\frac{1}{1+I_{\gamma }/K_{I_{\gamma }}}}}\\ & -\underset{I_{\gamma },T_{\alpha }\text{-induced transition from}{M}_2\mathrm{to}{M}_1}{\underbrace{M_2({\lambda }_{M_2{M}_1{I}_{\gamma }}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{M_2{M}_1{T}_{\alpha }}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & -\underset{\text{loss due to bursting}}{\underbrace{{\tilde{\mu }}_{M_2}{M}_2\frac{P_2^2}{P_2^2+{\left(\frac{N{M}_2}{150}\right)}^2}}}-\underset{\text{loss due to death}}{\underbrace{{\mu }_{M_2}{M}_2}},\end{array}$(3) where the fourth term on the right-hand side accounts for bursting. The remaining terms are similar or complementary to the terms in the right-hand side of Equation (2(2) $\begin{array}{rl}\frac{\mathrm{d}{M}_1}{\mathrm{d}t}& =\underset{\mathrm{recruitment}}{\underbrace{{\lambda }_{M_1}}}+\underset{I_{\gamma },T_{\alpha }\text{-induced transition from }{M}_2\mathrm{to}{M}_1}{\underbrace{M_2({\lambda }_{M_2{M}_1{I}_{\gamma }}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{M_2{M}_1{T}_{\alpha }}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & -\underset{\text{loss due to transition to }{M}_2}{\underbrace{M_1({\lambda }_{M_1{M}_2{I}_{10}}\frac{I_{10}}{I_{10}+K_{I_{10}}}+{\lambda }_{M_1{M}_2{T}_{\beta }}\frac{T_{\beta }}{T_{\beta }+K_{T_{\beta }}})\frac{1}{1+I_{\gamma }/K_{I_{\gamma }}}}}\\ & -\underset{\text{loss due to bursting}}{\underbrace{{\tilde{\mu }}_{M_1}{M}_1\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}}}-\underset{\text{loss due to death}}{\underbrace{{\mu }_{M_1}{M}_1}}.\end{array}$(2) ).

2.3. Equations for the parasites $P_1$ and $P_2$

The parasites $P_1$ reside in $M_1$ and the parasites $P_2$ reside in $M_2$. The Leishmania density $P_1$ satisfies the equation: (4) $\begin{array}{rl}\frac{\mathrm{d}{P}_1}{\mathrm{d}t}& =\underset{\text{growth in}{M}_1}{\underbrace{{\lambda }_{P_1}{P}_1{(1-\frac{P_1}{\frac{N_{\ast }{M}_1}{150}})}^{+}}}+\underset{I_{\gamma },T_{\alpha }\text{-induced transition from}{M}_2\mathrm{to}{M}_1}{\underbrace{P_2({\lambda }_{M_2{M}_1{I}_{\gamma }}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{M_2{M}_1{T}_{\alpha }}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & +\underset{\text{gain after burst of}{M}_2}{\underbrace{{\tilde{\mu }}_{M_2}\frac{N{M}_2}{150}\frac{P_2^2}{P_2^2+{\left(\frac{N{M}_2}{150}\right)}^2}\theta }}+\underset{\text{gain after burst of}{M}_1}{\underbrace{{\tilde{\mu }}_{M_1}\frac{N{M}_1}{150}\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}\theta }}\\ & -\underset{\text{loss due to burst of}{M}_1}{\underbrace{{\tilde{\mu }}_{M_1}\frac{N{M}_1}{150}\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}}}\\ & -\underset{\text{loss due to transition from}{M}_1\text{to}{M}_2}{\underbrace{P_1({\lambda }_{M_1{M}_2{I}_{10}}\frac{I_{10}}{I_{10}+K_{I_{10}}}+{\lambda }_{M_1{M}_2{T}_{\beta }}\frac{T_{\beta }}{T_{\beta }+K_{T_{\beta }}})\frac{1}{1+I_{\gamma }/K_{I_{\gamma }}}}}\\ & -\underset{\text{death augmented by}NO,I_{\gamma }\text{and}{T}_{\alpha }}{\underbrace{{\mu }_{P_1}{P}_1(1+{\lambda }_{NO{P}_1}\frac{NO}{NO+K_{NO}}+{\lambda }_{I_{\gamma }{P}_1}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{T_{\alpha }{P}_1}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & -\underset{\text{death by}{M}_1\text{apoptosis}}{\underbrace{{\mu }_{M_1}{N}_1{P}_1}},\end{array}$(4) where $0<\theta <1,0<{\tilde{\mu }}_{M_1}<{\tilde{\mu }}_{M_2}$.

The first term on the right-hand side of Equation (4(4) $\begin{array}{rl}\frac{\mathrm{d}{P}_1}{\mathrm{d}t}& =\underset{\text{growth in}{M}_1}{\underbrace{{\lambda }_{P_1}{P}_1{(1-\frac{P_1}{\frac{N_{\ast }{M}_1}{150}})}^{+}}}+\underset{I_{\gamma },T_{\alpha }\text{-induced transition from}{M}_2\mathrm{to}{M}_1}{\underbrace{P_2({\lambda }_{M_2{M}_1{I}_{\gamma }}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{M_2{M}_1{T}_{\alpha }}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & +\underset{\text{gain after burst of}{M}_2}{\underbrace{{\tilde{\mu }}_{M_2}\frac{N{M}_2}{150}\frac{P_2^2}{P_2^2+{\left(\frac{N{M}_2}{150}\right)}^2}\theta }}+\underset{\text{gain after burst of}{M}_1}{\underbrace{{\tilde{\mu }}_{M_1}\frac{N{M}_1}{150}\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}\theta }}\\ & -\underset{\text{loss due to burst of}{M}_1}{\underbrace{{\tilde{\mu }}_{M_1}\frac{N{M}_1}{150}\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}}}\\ & -\underset{\text{loss due to transition from}{M}_1\text{to}{M}_2}{\underbrace{P_1({\lambda }_{M_1{M}_2{I}_{10}}\frac{I_{10}}{I_{10}+K_{I_{10}}}+{\lambda }_{M_1{M}_2{T}_{\beta }}\frac{T_{\beta }}{T_{\beta }+K_{T_{\beta }}})\frac{1}{1+I_{\gamma }/K_{I_{\gamma }}}}}\\ & -\underset{\text{death augmented by}NO,I_{\gamma }\text{and}{T}_{\alpha }}{\underbrace{{\mu }_{P_1}{P}_1(1+{\lambda }_{NO{P}_1}\frac{NO}{NO+K_{NO}}+{\lambda }_{I_{\gamma }{P}_1}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{T_{\alpha }{P}_1}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & -\underset{\text{death by}{M}_1\text{apoptosis}}{\underbrace{{\mu }_{M_1}{N}_1{P}_1}},\end{array}$(4) ) is a logistic growth of Leishmania within $M_1$, where $N_{\ast }$ is the carrying capacity for $P_1$ in $M_1$; $N_{\ast }\ge N$ [9, 21, 38]. Note that if $P_1$ exceeds $N{M}_1/150$, there is no growth.

When $M_2$ polarizes to $M_1$ (the second term on the right-hand side of Equation (2(2) $\begin{array}{rl}\frac{\mathrm{d}{M}_1}{\mathrm{d}t}& =\underset{\mathrm{recruitment}}{\underbrace{{\lambda }_{M_1}}}+\underset{I_{\gamma },T_{\alpha }\text{-induced transition from }{M}_2\mathrm{to}{M}_1}{\underbrace{M_2({\lambda }_{M_2{M}_1{I}_{\gamma }}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{M_2{M}_1{T}_{\alpha }}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & -\underset{\text{loss due to transition to }{M}_2}{\underbrace{M_1({\lambda }_{M_1{M}_2{I}_{10}}\frac{I_{10}}{I_{10}+K_{I_{10}}}+{\lambda }_{M_1{M}_2{T}_{\beta }}\frac{T_{\beta }}{T_{\beta }+K_{T_{\beta }}})\frac{1}{1+I_{\gamma }/K_{I_{\gamma }}}}}\\ & -\underset{\text{loss due to bursting}}{\underbrace{{\tilde{\mu }}_{M_1}{M}_1\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}}}-\underset{\text{loss due to death}}{\underbrace{{\mu }_{M_1}{M}_1}}.\end{array}$(2) )), the parasites $P_2$ in $M_2$ become parasites in $M_1$, and this is accounted for by the second term on the right-hand side of Equation (4(4) $\begin{array}{rl}\frac{\mathrm{d}{P}_1}{\mathrm{d}t}& =\underset{\text{growth in}{M}_1}{\underbrace{{\lambda }_{P_1}{P}_1{(1-\frac{P_1}{\frac{N_{\ast }{M}_1}{150}})}^{+}}}+\underset{I_{\gamma },T_{\alpha }\text{-induced transition from}{M}_2\mathrm{to}{M}_1}{\underbrace{P_2({\lambda }_{M_2{M}_1{I}_{\gamma }}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{M_2{M}_1{T}_{\alpha }}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & +\underset{\text{gain after burst of}{M}_2}{\underbrace{{\tilde{\mu }}_{M_2}\frac{N{M}_2}{150}\frac{P_2^2}{P_2^2+{\left(\frac{N{M}_2}{150}\right)}^2}\theta }}+\underset{\text{gain after burst of}{M}_1}{\underbrace{{\tilde{\mu }}_{M_1}\frac{N{M}_1}{150}\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}\theta }}\\ & -\underset{\text{loss due to burst of}{M}_1}{\underbrace{{\tilde{\mu }}_{M_1}\frac{N{M}_1}{150}\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}}}\\ & -\underset{\text{loss due to transition from}{M}_1\text{to}{M}_2}{\underbrace{P_1({\lambda }_{M_1{M}_2{I}_{10}}\frac{I_{10}}{I_{10}+K_{I_{10}}}+{\lambda }_{M_1{M}_2{T}_{\beta }}\frac{T_{\beta }}{T_{\beta }+K_{T_{\beta }}})\frac{1}{1+I_{\gamma }/K_{I_{\gamma }}}}}\\ & -\underset{\text{death augmented by}NO,I_{\gamma }\text{and}{T}_{\alpha }}{\underbrace{{\mu }_{P_1}{P}_1(1+{\lambda }_{NO{P}_1}\frac{NO}{NO+K_{NO}}+{\lambda }_{I_{\gamma }{P}_1}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{T_{\alpha }{P}_1}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & -\underset{\text{death by}{M}_1\text{apoptosis}}{\underbrace{{\mu }_{M_1}{N}_1{P}_1}},\end{array}$(4) ). When $M_2$ bursts, $N{M}_2$ parasites are released, with a total mass of $N{M}_2/150$. We assume that θ-fraction of them is ingested by $M_1$ macrophages, and this is represented by the third term on the right-hand side of Equation (4(4) $\begin{array}{rl}\frac{\mathrm{d}{P}_1}{\mathrm{d}t}& =\underset{\text{growth in}{M}_1}{\underbrace{{\lambda }_{P_1}{P}_1{(1-\frac{P_1}{\frac{N_{\ast }{M}_1}{150}})}^{+}}}+\underset{I_{\gamma },T_{\alpha }\text{-induced transition from}{M}_2\mathrm{to}{M}_1}{\underbrace{P_2({\lambda }_{M_2{M}_1{I}_{\gamma }}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{M_2{M}_1{T}_{\alpha }}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & +\underset{\text{gain after burst of}{M}_2}{\underbrace{{\tilde{\mu }}_{M_2}\frac{N{M}_2}{150}\frac{P_2^2}{P_2^2+{\left(\frac{N{M}_2}{150}\right)}^2}\theta }}+\underset{\text{gain after burst of}{M}_1}{\underbrace{{\tilde{\mu }}_{M_1}\frac{N{M}_1}{150}\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}\theta }}\\ & -\underset{\text{loss due to burst of}{M}_1}{\underbrace{{\tilde{\mu }}_{M_1}\frac{N{M}_1}{150}\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}}}\\ & -\underset{\text{loss due to transition from}{M}_1\text{to}{M}_2}{\underbrace{P_1({\lambda }_{M_1{M}_2{I}_{10}}\frac{I_{10}}{I_{10}+K_{I_{10}}}+{\lambda }_{M_1{M}_2{T}_{\beta }}\frac{T_{\beta }}{T_{\beta }+K_{T_{\beta }}})\frac{1}{1+I_{\gamma }/K_{I_{\gamma }}}}}\\ & -\underset{\text{death augmented by}NO,I_{\gamma }\text{and}{T}_{\alpha }}{\underbrace{{\mu }_{P_1}{P}_1(1+{\lambda }_{NO{P}_1}\frac{NO}{NO+K_{NO}}+{\lambda }_{I_{\gamma }{P}_1}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{T_{\alpha }{P}_1}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & -\underset{\text{death by}{M}_1\text{apoptosis}}{\underbrace{{\mu }_{M_1}{N}_1{P}_1}},\end{array}$(4) ); the remaining $(1-\theta )$-fraction are ingested by $M_2$ macrophages (the third term on the right-hand side of Equation (5(5) $\begin{array}{rl}\frac{\mathrm{d}{P}_2}{\mathrm{d}t}& =\underset{\text{growth in}{M}_2}{\underbrace{{\lambda }_{P_2}{P}_2{(1-\frac{P_2}{\frac{N_{\ast }{M}_2}{150}})}^{+}}}\\ & +\underset{\text{gain from transition from}{M}_1\text{to}{M}_2}{\underbrace{P_1({\lambda }_{M_1{M}_2{I}_{10}}\frac{I_{10}}{I_{10}+K_{I_{10}}}+{\lambda }_{M_1{M}_2{T}_{\beta }}\frac{T_{\beta }}{T_{\beta }+K_{T_{\beta }}})\frac{1}{1+I_{\gamma }/K_{I_{\gamma }}}}}\\ & +\underset{\text{gain after burst of}{M}_2}{\underbrace{{\tilde{\mu }}_{M_2}\frac{N{M}_2}{150}\frac{P_2^2}{P_2^2+{\left(\frac{N{M}_2}{150}\right)}^2}(1-\theta )}}+\underset{\text{gain after burst of}{M}_1}{\underbrace{{\tilde{\mu }}_{M_1}\frac{N{M}_1}{150}\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}(1-\theta )}}\\ & -\underset{\text{loss due to burst of}{M}_2}{\underbrace{{\tilde{\mu }}_{M_2}\frac{N{M}_2}{150}\frac{P_2^2}{P_2^2+{\left(\frac{N{M}_2}{150}\right)}^2}}}-\underset{I_{\gamma },T_{\alpha }\text{-induced transition from}{M}_2\text{to}{M}_1}{\underbrace{P_2({\lambda }_{M_2{M}_1{I}_{\gamma }}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{M_2{M}_1{T}_{\alpha }}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & -\underset{\text{death augmented by NO},I_{\gamma }\text{and}{T}_{\alpha }}{\underbrace{{\mu }_{P_2}{P}_2(1+{\lambda }_{NO{P}_2}\frac{NO}{NO+K_{NO}}+{\lambda }_{I_{\gamma }{P}_2}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{T_{\alpha }{P}_2}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & -\underset{\text{death by}{M}_2\text{apoptosis}}{\underbrace{{\mu }_{M_2}{N}_2{P}_2}}.\end{array}$(5) )). The fourth term on the right-hand side of Equation (4(4) $\begin{array}{rl}\frac{\mathrm{d}{P}_1}{\mathrm{d}t}& =\underset{\text{growth in}{M}_1}{\underbrace{{\lambda }_{P_1}{P}_1{(1-\frac{P_1}{\frac{N_{\ast }{M}_1}{150}})}^{+}}}+\underset{I_{\gamma },T_{\alpha }\text{-induced transition from}{M}_2\mathrm{to}{M}_1}{\underbrace{P_2({\lambda }_{M_2{M}_1{I}_{\gamma }}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{M_2{M}_1{T}_{\alpha }}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & +\underset{\text{gain after burst of}{M}_2}{\underbrace{{\tilde{\mu }}_{M_2}\frac{N{M}_2}{150}\frac{P_2^2}{P_2^2+{\left(\frac{N{M}_2}{150}\right)}^2}\theta }}+\underset{\text{gain after burst of}{M}_1}{\underbrace{{\tilde{\mu }}_{M_1}\frac{N{M}_1}{150}\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}\theta }}\\ & -\underset{\text{loss due to burst of}{M}_1}{\underbrace{{\tilde{\mu }}_{M_1}\frac{N{M}_1}{150}\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}}}\\ & -\underset{\text{loss due to transition from}{M}_1\text{to}{M}_2}{\underbrace{P_1({\lambda }_{M_1{M}_2{I}_{10}}\frac{I_{10}}{I_{10}+K_{I_{10}}}+{\lambda }_{M_1{M}_2{T}_{\beta }}\frac{T_{\beta }}{T_{\beta }+K_{T_{\beta }}})\frac{1}{1+I_{\gamma }/K_{I_{\gamma }}}}}\\ & -\underset{\text{death augmented by}NO,I_{\gamma }\text{and}{T}_{\alpha }}{\underbrace{{\mu }_{P_1}{P}_1(1+{\lambda }_{NO{P}_1}\frac{NO}{NO+K_{NO}}+{\lambda }_{I_{\gamma }{P}_1}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{T_{\alpha }{P}_1}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & -\underset{\text{death by}{M}_1\text{apoptosis}}{\underbrace{{\mu }_{M_1}{N}_1{P}_1}},\end{array}$(4) ) accounts for the gain of $P_1$ from bursting $M_1$ macrophages.

Recall that the density of parasites $P_1$ reduces when $M_1$ bursts or when $M_1$ polarizes to $M_2$. The fifth term of Equation (4(4) $\begin{array}{rl}\frac{\mathrm{d}{P}_1}{\mathrm{d}t}& =\underset{\text{growth in}{M}_1}{\underbrace{{\lambda }_{P_1}{P}_1{(1-\frac{P_1}{\frac{N_{\ast }{M}_1}{150}})}^{+}}}+\underset{I_{\gamma },T_{\alpha }\text{-induced transition from}{M}_2\mathrm{to}{M}_1}{\underbrace{P_2({\lambda }_{M_2{M}_1{I}_{\gamma }}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{M_2{M}_1{T}_{\alpha }}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & +\underset{\text{gain after burst of}{M}_2}{\underbrace{{\tilde{\mu }}_{M_2}\frac{N{M}_2}{150}\frac{P_2^2}{P_2^2+{\left(\frac{N{M}_2}{150}\right)}^2}\theta }}+\underset{\text{gain after burst of}{M}_1}{\underbrace{{\tilde{\mu }}_{M_1}\frac{N{M}_1}{150}\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}\theta }}\\ & -\underset{\text{loss due to burst of}{M}_1}{\underbrace{{\tilde{\mu }}_{M_1}\frac{N{M}_1}{150}\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}}}\\ & -\underset{\text{loss due to transition from}{M}_1\text{to}{M}_2}{\underbrace{P_1({\lambda }_{M_1{M}_2{I}_{10}}\frac{I_{10}}{I_{10}+K_{I_{10}}}+{\lambda }_{M_1{M}_2{T}_{\beta }}\frac{T_{\beta }}{T_{\beta }+K_{T_{\beta }}})\frac{1}{1+I_{\gamma }/K_{I_{\gamma }}}}}\\ & -\underset{\text{death augmented by}NO,I_{\gamma }\text{and}{T}_{\alpha }}{\underbrace{{\mu }_{P_1}{P}_1(1+{\lambda }_{NO{P}_1}\frac{NO}{NO+K_{NO}}+{\lambda }_{I_{\gamma }{P}_1}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{T_{\alpha }{P}_1}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & -\underset{\text{death by}{M}_1\text{apoptosis}}{\underbrace{{\mu }_{M_1}{N}_1{P}_1}},\end{array}$(4) ) is the loss of parasites $P_1$ due to the bursting of $M_1$; this term is the parasite-equivalence of the third term in the right-hand side of Equation (2(2) $\begin{array}{rl}\frac{\mathrm{d}{M}_1}{\mathrm{d}t}& =\underset{\mathrm{recruitment}}{\underbrace{{\lambda }_{M_1}}}+\underset{I_{\gamma },T_{\alpha }\text{-induced transition from }{M}_2\mathrm{to}{M}_1}{\underbrace{M_2({\lambda }_{M_2{M}_1{I}_{\gamma }}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{M_2{M}_1{T}_{\alpha }}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & -\underset{\text{loss due to transition to }{M}_2}{\underbrace{M_1({\lambda }_{M_1{M}_2{I}_{10}}\frac{I_{10}}{I_{10}+K_{I_{10}}}+{\lambda }_{M_1{M}_2{T}_{\beta }}\frac{T_{\beta }}{T_{\beta }+K_{T_{\beta }}})\frac{1}{1+I_{\gamma }/K_{I_{\gamma }}}}}\\ & -\underset{\text{loss due to bursting}}{\underbrace{{\tilde{\mu }}_{M_1}{M}_1\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}}}-\underset{\text{loss due to death}}{\underbrace{{\mu }_{M_1}{M}_1}}.\end{array}$(2) ) where the average mass of parasites in macrophage $M_1$ at bursting time is $N{M}_1/150$. The sixth term is the loss in $P_1$ due to the transition from $M_1$ to $M_2$. Similarly, this sixth term is the parasite-equivalence of the second term in Equation (2(2) $\begin{array}{rl}\frac{\mathrm{d}{M}_1}{\mathrm{d}t}& =\underset{\mathrm{recruitment}}{\underbrace{{\lambda }_{M_1}}}+\underset{I_{\gamma },T_{\alpha }\text{-induced transition from }{M}_2\mathrm{to}{M}_1}{\underbrace{M_2({\lambda }_{M_2{M}_1{I}_{\gamma }}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{M_2{M}_1{T}_{\alpha }}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & -\underset{\text{loss due to transition to }{M}_2}{\underbrace{M_1({\lambda }_{M_1{M}_2{I}_{10}}\frac{I_{10}}{I_{10}+K_{I_{10}}}+{\lambda }_{M_1{M}_2{T}_{\beta }}\frac{T_{\beta }}{T_{\beta }+K_{T_{\beta }}})\frac{1}{1+I_{\gamma }/K_{I_{\gamma }}}}}\\ & -\underset{\text{loss due to bursting}}{\underbrace{{\tilde{\mu }}_{M_1}{M}_1\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}}}-\underset{\text{loss due to death}}{\underbrace{{\mu }_{M_1}{M}_1}}.\end{array}$(2) ). The seventh term on the right-hand-side of Equation (4(4) $\begin{array}{rl}\frac{\mathrm{d}{P}_1}{\mathrm{d}t}& =\underset{\text{growth in}{M}_1}{\underbrace{{\lambda }_{P_1}{P}_1{(1-\frac{P_1}{\frac{N_{\ast }{M}_1}{150}})}^{+}}}+\underset{I_{\gamma },T_{\alpha }\text{-induced transition from}{M}_2\mathrm{to}{M}_1}{\underbrace{P_2({\lambda }_{M_2{M}_1{I}_{\gamma }}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{M_2{M}_1{T}_{\alpha }}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & +\underset{\text{gain after burst of}{M}_2}{\underbrace{{\tilde{\mu }}_{M_2}\frac{N{M}_2}{150}\frac{P_2^2}{P_2^2+{\left(\frac{N{M}_2}{150}\right)}^2}\theta }}+\underset{\text{gain after burst of}{M}_1}{\underbrace{{\tilde{\mu }}_{M_1}\frac{N{M}_1}{150}\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}\theta }}\\ & -\underset{\text{loss due to burst of}{M}_1}{\underbrace{{\tilde{\mu }}_{M_1}\frac{N{M}_1}{150}\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}}}\\ & -\underset{\text{loss due to transition from}{M}_1\text{to}{M}_2}{\underbrace{P_1({\lambda }_{M_1{M}_2{I}_{10}}\frac{I_{10}}{I_{10}+K_{I_{10}}}+{\lambda }_{M_1{M}_2{T}_{\beta }}\frac{T_{\beta }}{T_{\beta }+K_{T_{\beta }}})\frac{1}{1+I_{\gamma }/K_{I_{\gamma }}}}}\\ & -\underset{\text{death augmented by}NO,I_{\gamma }\text{and}{T}_{\alpha }}{\underbrace{{\mu }_{P_1}{P}_1(1+{\lambda }_{NO{P}_1}\frac{NO}{NO+K_{NO}}+{\lambda }_{I_{\gamma }{P}_1}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{T_{\alpha }{P}_1}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & -\underset{\text{death by}{M}_1\text{apoptosis}}{\underbrace{{\mu }_{M_1}{N}_1{P}_1}},\end{array}$(4) ) accounts for the killing of parasites within $M_1$, which is enhanced by NO [42, 82], $I_{\gamma }$ [1, 2, 6], and $T_{\alpha }$ [52, 56, 58, 96]. We assume that when a macrophage of $M_1$ dies, all the Leishmania inside it also die. Accordingly, we get the corresponding death rate of $P_1$ to be ${\mu }_{M_1}{N}_1{P}_1$, where $N_1$ is the average number of Leishmania in one $M_1$ macrophage at the time when the macrophage undergoes apoptosis.

The equation for Leishmania density $P_2$ is the following: (5) $\begin{array}{rl}\frac{\mathrm{d}{P}_2}{\mathrm{d}t}& =\underset{\text{growth in}{M}_2}{\underbrace{{\lambda }_{P_2}{P}_2{(1-\frac{P_2}{\frac{N_{\ast }{M}_2}{150}})}^{+}}}\\ & +\underset{\text{gain from transition from}{M}_1\text{to}{M}_2}{\underbrace{P_1({\lambda }_{M_1{M}_2{I}_{10}}\frac{I_{10}}{I_{10}+K_{I_{10}}}+{\lambda }_{M_1{M}_2{T}_{\beta }}\frac{T_{\beta }}{T_{\beta }+K_{T_{\beta }}})\frac{1}{1+I_{\gamma }/K_{I_{\gamma }}}}}\\ & +\underset{\text{gain after burst of}{M}_2}{\underbrace{{\tilde{\mu }}_{M_2}\frac{N{M}_2}{150}\frac{P_2^2}{P_2^2+{\left(\frac{N{M}_2}{150}\right)}^2}(1-\theta )}}+\underset{\text{gain after burst of}{M}_1}{\underbrace{{\tilde{\mu }}_{M_1}\frac{N{M}_1}{150}\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}(1-\theta )}}\\ & -\underset{\text{loss due to burst of}{M}_2}{\underbrace{{\tilde{\mu }}_{M_2}\frac{N{M}_2}{150}\frac{P_2^2}{P_2^2+{\left(\frac{N{M}_2}{150}\right)}^2}}}-\underset{I_{\gamma },T_{\alpha }\text{-induced transition from}{M}_2\text{to}{M}_1}{\underbrace{P_2({\lambda }_{M_2{M}_1{I}_{\gamma }}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{M_2{M}_1{T}_{\alpha }}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & -\underset{\text{death augmented by NO},I_{\gamma }\text{and}{T}_{\alpha }}{\underbrace{{\mu }_{P_2}{P}_2(1+{\lambda }_{NO{P}_2}\frac{NO}{NO+K_{NO}}+{\lambda }_{I_{\gamma }{P}_2}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{T_{\alpha }{P}_2}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & -\underset{\text{death by}{M}_2\text{apoptosis}}{\underbrace{{\mu }_{M_2}{N}_2{P}_2}}.\end{array}$(5) The various terms are similar, or complementary, to the terms in Equation (4(4) $\begin{array}{rl}\frac{\mathrm{d}{P}_1}{\mathrm{d}t}& =\underset{\text{growth in}{M}_1}{\underbrace{{\lambda }_{P_1}{P}_1{(1-\frac{P_1}{\frac{N_{\ast }{M}_1}{150}})}^{+}}}+\underset{I_{\gamma },T_{\alpha }\text{-induced transition from}{M}_2\mathrm{to}{M}_1}{\underbrace{P_2({\lambda }_{M_2{M}_1{I}_{\gamma }}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{M_2{M}_1{T}_{\alpha }}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & +\underset{\text{gain after burst of}{M}_2}{\underbrace{{\tilde{\mu }}_{M_2}\frac{N{M}_2}{150}\frac{P_2^2}{P_2^2+{\left(\frac{N{M}_2}{150}\right)}^2}\theta }}+\underset{\text{gain after burst of}{M}_1}{\underbrace{{\tilde{\mu }}_{M_1}\frac{N{M}_1}{150}\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}\theta }}\\ & -\underset{\text{loss due to burst of}{M}_1}{\underbrace{{\tilde{\mu }}_{M_1}\frac{N{M}_1}{150}\frac{P_1^2}{P_1^2+{\left(\frac{N{M}_1}{150}\right)}^2}}}\\ & -\underset{\text{loss due to transition from}{M}_1\text{to}{M}_2}{\underbrace{P_1({\lambda }_{M_1{M}_2{I}_{10}}\frac{I_{10}}{I_{10}+K_{I_{10}}}+{\lambda }_{M_1{M}_2{T}_{\beta }}\frac{T_{\beta }}{T_{\beta }+K_{T_{\beta }}})\frac{1}{1+I_{\gamma }/K_{I_{\gamma }}}}}\\ & -\underset{\text{death augmented by}NO,I_{\gamma }\text{and}{T}_{\alpha }}{\underbrace{{\mu }_{P_1}{P}_1(1+{\lambda }_{NO{P}_1}\frac{NO}{NO+K_{NO}}+{\lambda }_{I_{\gamma }{P}_1}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{T_{\alpha }{P}_1}\frac{T_{\alpha }}{T_{\alpha }+K_{T_{\alpha }}})}}\\ & -\underset{\text{death by}{M}_1\text{apoptosis}}{\underbrace{{\mu }_{M_1}{N}_1{P}_1}},\end{array}$(4) ).

2.4. Equation for T cells T

The density of T cells satisfies the equation: (6) $\begin{array}{rl}\frac{\text{d}T}{\text{d}t}& =\underset{I_2\text{-induced proliferation of T cells}}{\underbrace{{\lambda }_{T}T\frac{I_2}{I_2+K_{I_2}}}}+\underset{\text{naive T cells activation}}{\underbrace{{\lambda }_{DT}{T}_0\frac{D}{D+K_D}\frac{I_{12}}{I_{12}+K_{I_{12}}}}}\\ & +\underset{I_{12}\text{-induced proliferation of T cells}}{\underbrace{{\lambda }_{MT}{T}_0\frac{M_1}{M_1+K_M}\frac{I_{12}}{I_{12}+K_{I_{12}}}\frac{1}{1+I_{10}/K_{I_{10}}}}}+\underset{\text{to become memory T cells}}{\underbrace{{\lambda }_{VT}TV}}-\underset{\text{loss due to death}}{\underbrace{{\mu }_{T}T}}.\end{array}$(6) The first term on the right-hand side represents the $I_2$-induced proliferation of T cells [22, 88]. The second term is the activation of T cells by contact with dendritic cells in the $I_{12}$ environment. The third term represents activation of naive T cells ($T_0$ cells) by contact with the inflammatory M1 macrophages in the $I_{12}$ environment. The fourth term is the production of T cells that will become memory T cells.

2.5. Equations for cytokines

Equation for interleukin-2 ($I_2$): $I_2$ production by T cells [1] is enhanced by SAG [7, 37] and LEISH-F1+MPL-SE [16]. Hence, (7) $\frac{\text{d}{I}_2}{\text{d}t}=\underset{\text{secretion by}T}{\underbrace{{\lambda }_{T{I}_2}T(1+{\lambda }_{S{I}_2}S+{\lambda }_{V{I}_2}V)}}-\underset{\text{decay}}{\underbrace{{\mu }_{I_2}{I}_2}}.$(7)

Equation for interleukin-10 ($I_{10}$): $I_{10}$ is produced primarily by $M_2$ [1, 61], and this production is enhanced by MA [34], so that (8) $\frac{\text{d}{I}_{10}}{\text{d}t}=\underset{\text{secretion by}{M}_2}{\underbrace{{\lambda }_{M_2{I}_{10}}{M}_2(1+{\lambda }_{A{I}_{10}}A)}}-\underset{\text{decay}}{\underbrace{{\mu }_{I_{10}}{I}_{10}}}.$(8)

Equation for interleukin-12 ($I_{12}$): $I_{12}$ is produced primarily by $M_1$, and this process is inhibited by $I_{10}$ (and $I_{13}$, which we combined with $I_{10}$) [1, 10, 70]. SAG enhances the production of IL-12 [7, 31, 37]. Hence, (9) $\frac{\text{d}{I}_{12}}{\text{d}t}=\underset{\text{secretion by}{M}_1}{\underbrace{{\lambda }_{M_1{I}_{12}}{M}_1\frac{1}{1+I_{10}/K_{I_{10}}}(1+{\lambda }_{S{I}_{12}}S)}}-\underset{\text{decay}}{\underbrace{{\mu }_{I_{12}}{I}_{12}}}.$(9)

Equation for interferon-γ ($I_{\gamma }$):

$I_{\gamma }$ is secreted by the T cells [2, 10, 41], a process that is enhanced in presence of vaccine LEISH-F1+MPL-SE [25], so that (10) $\frac{\text{d}{I}_{\gamma }}{\text{d}t}=\underset{\text{secretion by}T}{\underbrace{{\lambda }_{T{I}_{\gamma }}(1+c_{I_{\gamma }\text{adj}})T(1+{\lambda }_{V{I}_{\gamma }}V)}}-\underset{\text{decay}}{\underbrace{{\mu }_{I_{\gamma }}{I}_{\gamma }}},$(10) where $c_{I_{\gamma }\text{adj}}$ is the augmented secretion rate of IFN-γ due to the MPL-SE adjuvant. We set $c_{I_{\gamma }\text{adj}}=0$ if the MPL-SE adjuvant is not given, or if the full vaccine LEISH-F1+MPL-SE adjuvant is given.

Equation for nitric oxide (NO): Nitric oxide is produced by $M_1$ and $M_2$ macrophages, and this production is enhanced by $I_{\gamma }$ [13, 76, 82], SAG [37], and MA [34]. Hence, (11) $\frac{\text{d}NO}{\text{d}t}=\underset{\text{secretion by}{M}_1\text{and}{M}_2}{\underbrace{{\lambda }_{MNO}(M_1+M_2)(1+{\lambda }_{I_{\gamma }NO}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{SNO}S+{\lambda }_{ANO}A)}}-\underset{\text{decay}}{\underbrace{{\mu }_{NO}NO}}.$(11)

Equations for TNF-α ($T_{\alpha }$) and TGF-β ($T_{\beta }$): TNF-α is produced by T cells [62] and by M1 macrophages [55, 71], and V [16, 25] and A [34] increase the production of TNF-α by macrophages. Hence, (12) $\frac{\text{d}{T}_{\alpha }}{\text{d}t}=\underset{\text{secretion by}{M}_1\text{augmented by}V\text{and}A}{\underbrace{{\lambda }_{M_1{T}_{\alpha }}{M}_1(1+{\lambda }_{V{T}_{\alpha }}V+{\lambda }_{A{T}_{\alpha }}A)}}+\underset{\text{secretion by}T}{\underbrace{{\lambda }_{T{T}_{\alpha }}T}}-\underset{\text{decay}}{\underbrace{{\mu }_{T_{\alpha }}{T}_{\alpha }}}.$(12) TGF-β is secreted by M2 macrophages [19, 54, 66], so that (13) $\frac{\text{d}{T}_{\beta }}{\text{d}t}=\underset{\text{secretion by}{M}_2}{\underbrace{{\lambda }_{M_2{T}_{\beta }}{M}_2}}-\underset{\text{decay}}{\underbrace{{\mu }_{T_{\beta }}{T}_{\beta }}}.$(13)

 

2.6. Equations for MA (A), SAG (S) and LEISH-F1+MPL-SE vaccine (V)

We use a PK/PD (pharmacokinetic/pharmacodynamic) model for the drugs. We denote by ${\gamma }_A$ an amount of MA administered at times $t_0=0,t_1,t_2,\dots ,$ and by $f_A(t)$ the kinetic profile of the drug as it decays in time. We assume that the decay is exponential, so that $f_A(t)=\sum _{j=0}^k{\gamma }_A{e}^{-{\beta }_A(t-t_j)}\text{for}{t}_{k-1}<t<t_k,k=1,2,\dots ,$where ${\beta }_A$ is some positive parameter. The PD term accounts for depletion of the drug through its effect on M1 macrophages, namely, inducing the secretion of NO. We assume that this term has the form ${\mu }_{MA}(M_1+M_2)A$ where ${\mu }_{MA}$ is a positive parameter. Then the dynamic of A takes the following form: (14) $\frac{\mathrm{d}A}{\mathrm{d}t}=f_A(t)-\underset{\text{depletion through acting on M1}}{\underbrace{{\mu }_{MA}(M_1+M_2)A}}-\underset{\text{degradation}}{\underbrace{{\mu }_{A}A}},$(14) where ${\mu }_{A}A$ is the intrinsic degradation of the drug.

The equations for SAG and vaccine LEISH-F1+MPL-SE have similar forms: (15) $\begin{array}{rl}\frac{\text{d}S}{\text{d}t}& =f_S(t)-\underset{\text{depletion through acting on M1 and M2}}{\underbrace{{\mu }_{SM}(M_1+M_2)S}}-\underset{\text{depletion through acting on}T}{\underbrace{{\mu }_{ST}T}}-\underset{\text{degradation}}{\underbrace{{\mu }_{S}S}},\end{array}$(15) (16) $\begin{array}{rl}\frac{\text{d}V}{\text{d}t}& =f_V(t)-\underset{\text{depletion through acting on}T}{\underbrace{{\mu }_{VT}TV}}-\underset{\text{degradation}}{\underbrace{{\mu }_{V}V}},\end{array}$(16) where $f_S(t)$ and $f_V(t)$ have the same form as $f_A(t)$ with some parameters ${\beta }_S,{\beta }_V,{\gamma }_S,{\gamma }_V$.

3. Model simulations

All computations are done using the Python ODE solver odeint(), which uses a fourth-order Runge–Kutta scheme. In all the simulations we take an initial load of Leishmania parasites, while all other initial values are taken to be close, but not necessarily identical, to their steady state estimated in Section 5 (in units of g/cm${}^3$):

Table

$D(0)=4\times {10}^{-4}$	$M_1(0)=1.1\times {10}^{-3}$	$M_2=4.65\times {10}^{-4}$	$T(0)=3.5\times {10}^{-4}$
$P_1(0)=4.65\times {10}^{-5}$	$P_2(0)=3.5\times {10}^{-5}$	$I_2(0)=7.5\times {10}^{-11}$	$I_{10}(0)=2.4\times {10}^{-12}$
$I_{12}(0)=6.4\times {10}^{-12}$	$I_{\gamma }(0)=2.3\times {10}^{-11}$	$NO(0)=4.73\times {10}^{-8}$	$T_{\alpha }(0)=1.7\times {10}^{-11}$
$T_{\beta }(0)=2.4\times {10}^{-12}$	$A(0)=0$	$S(0)=0$	$V(0)=0$

3.1. Simulations with no drugs

Figures 2 and 3 show simulations of all the model variables in a no-drug case, when a leishmaniasis patient does not heal (Figure 2) and when a patient heals (Figure 3); correspondingly, the killing rates of parasites by $I_{\gamma },T_{\alpha }$ and NO are larger in Figure 3 by a factor of 3.

Figure 2. No self-healing; ${\lambda }_{NO{P}_1}={\lambda }_{I_{\gamma }{P}_1}={\lambda }_{I_{\gamma }{P}_2}={\lambda }_{T_{\alpha }{P}_1}={\lambda }_{T_{\alpha }{P}_2}=0.062$, ${\lambda }_{NO{P}_2}=0.56$. (a) All variables with no self-healing and (b) all macrophages and parasites with no self-healing.

Figure 3. Self-healing; (factor of 3) ${\lambda }_{NO{P}_1}={\lambda }_{NO{P}_2}={\lambda }_{I_{\gamma }{P}_1}={\lambda }_{I_{\gamma }{P}_2}={\lambda }_{T_{\alpha }{P}_1}={\lambda }_{T_{\alpha }{P}_2}=0.186$, ${\lambda }_{NO{P}_2}=1.67$. (a) All variables with self-healing and (b) all macrophages and parasites with self-healing.

In Figures 2 and 3 we see that, as t increases, cells and cytokines converge to approximately the half-saturation values listed in Table 2, which were estimated in Section 5. This shows consistency in the parameter estimates. We also note that the steady state of $I_{\gamma }$ is larger than that of $T_{\alpha }$ ($I_{\gamma }\sim 1.7T_{\alpha }$), which is in agreement with [16] (Figure 3(a)).

Table 2. Descriptions, values, and sources of the model parameters.

Parameter	Description	Value/range	Reference
$K_{I_2}$	Half-saturation of IL-2	$14\times {10}^{-11}$ g/cm${}^3$	[8]est.
$K_{I_{10}}$	Half-saturation of IL-10	$1.53\times {10}^{-12}$ g/cm${}^3$	[26]est.
$K_{I_{12}}$	Half-saturation of IL-12	$9\times {10}^{-12}$ g/cm${}^3$	[3]est.
$K_{I_{\gamma }}$	Half-saturation of IFN-γ	$4.37\times {10}^{-11}$ g/cm${}^3$	[26]est.
$K_{T_{\alpha }}$	Half-saturation of TNF-α	$21.83\times {10}^{-12}$ g/cm${}^3$	[48]est.
$K_{T_{\beta }}$	Half-saturation of TGF-β	$1.53\times {10}^{-12}$ g/cm${}^3$	fitted
$K_{NO}$	Half-saturation of NO	$5.36\times {10}^{-8}$ g/cm${}^3$	[4, 30]est.
$K_D$	Half-saturation of DCs	$7\times {10}^{-4}$ g/cm${}^3$	[79]
$K_{M_1}$	Half-saturation of M1 macrophages	$1.2\times {10}^{-3}$ g/cm${}^3$	[79, 80]est.
$K_{M_2}$	Half-saturation of M2 macrophages	$3\times {10}^{-4}$ g/cm${}^3$	[79, 80]est.
$K_P$	Half-saturation of parasites	$5\times {10}^{-4}$ g/cm${}^3$	[79, 80]est.
$K_T$	Half-saturation of T cells	$6.5\times {10}^{-4}$ g/cm${}^3$	[49]est.
${\lambda }_{M_2{M}_1{I}_{\gamma }}$	IFN-γ-induced transition rate from $M_2$ to $M_1$	0.69 d${}^{-1}$	[39, 67]
${\lambda }_{M_2{M}_1{T}_{\alpha }}$	TNF-α-induced transition rate from $M_2$ to $M_1$	0.69 d${}^{-1}$	est.
${\lambda }_{M_1{M}_2{I}_{10}}$	$I_{10}$-induced transition rate from $M_1$ to $M_2$	0.27 d${}^{-1}$	fitted
${\lambda }_{M_1{M}_2{T}_{\beta }}$	$T_{\beta }$-induced transition rate from $M_1$ to $M_2$	0.27 d${}^{-1}$	fitted
${\lambda }_{T{I}_2}$	Production rate of $I_2$ by T cells	$7.1\times {10}^{-5}$ d${}^{-1}$	est.
${\lambda }_{M_1{I}_{12}}$	Production rate of IL-12 by $M_1$	$2.1\times {10}^{-8}$ d${}^{-1}$	est.
${\lambda }_{M_2{I}_{10}}$	Production rate of IL-10 by $M_2$ macrophages	$7.4\times {10}^{-8}$ d${}^{-1}$	est.
${\lambda }_{M_1{T}_{\alpha }}$	Production rate of TNF-α by M1 macrophages	$2.35\times {10}^{-6}$ d${}^{-1}$	est.
${\lambda }_{T{T}_{\alpha }}$	Production rate of TNF-α by T cells	$2.35\times {10}^{-6}$ d${}^{-1}$	est.
${\lambda }_{T{I}_{\gamma }}$	Production rate of IFN-γ by T cells	$2.22\times {10}^{-6}$ d${}^{-1}$	est.
${\lambda }_{M_2{T}_{\beta }}$	production rate of TGF-β by M2 macrophages	$2.04\times {10}^{-6}$ d${}^{-1}$	est.
${\lambda }_{MNO}$	Production rate of NO by macrophages	$2.4\times {10}^{-3}$ d${}^{-1}$	est.

In both figures, the initial values of the inflammatory cytokines are below their steady states (or half-saturation) and this results in the initial increase in $P_1+P_2$ (and their corresponding hosts $M_1+M_2$). But, following this initial increase, $P_1+P_2$ starts to continuously decrease; in the healing case of Figure 3(b), $P_1+P_2\to 0$ while $M_1+M_2$ is continuously increasing, whereas in the non-healing case of Figure 2(b), $P_1+P_2$ remains positive, with $P_1+P_2>0.43\times {10}^{-4}$ g/cm${}^3$ and $M_1+M_2$ decreases below the half-saturation level.

Figure 4 shows additional profiles of $M_1+M_2$ and $P_1+P_2$ for some intermediate values of killing rates. The competition between the killing of parasites and their growth within macrophages results in oscillations in $P_1+P_2$, which eventually subside. The oscillations in $M_1+M_2$ are associated with their bursting by $P_1+P_2$.

Figure 4. All macrophages and parasites with no healing: (factor of 1.7) ${\lambda }_{NO{P}_1}={\lambda }_{I_{\gamma }{P}_1}={\lambda }_{I_{\gamma }{P}_2}={\lambda }_{T_{\alpha }{P}_1}={\lambda }_{T_{\alpha }{P}_2}=0.122$, ${\lambda }_{NO{P}_2}=1.106$.

3.2. Simulations with drugs in clinical trials

3.2.1. Meglumine antimoniate and vaccine

Nascimento et al. [65] conducted clinical trials with a combination of leishmania vaccine LEISH-F1+MPL-SE, or its adjuvant MPL-SE, and chemotherapy compound MA. The effect of MPL-SE is to increase the production of IFN-γ [74]. The clinical trials included three groups of CL patients:

1. All groups were treated with MA (A) in 4 cycles of 21 days; in the first 10 days they were injected with 10 mg/kg daily, and the remaining 11 days, they were at rest.

2. Group 1 receives shots of the vaccine LEISH-F1+MPL-SE (V) at Days 0, 28, and 56, with LEISH-F1 at 510or 20 µ g, and MPL-SE at 25 µ g.

3. Group 2 receives adjuvant MPL-SE shots at Days 0, 28, and 56, at 25 µ g.

4. Group 3 receives only MA.

Figure 5 represents the schedule of the treatments.

Figure 5. Schedule for drug administration. As in [65], shots of the vaccine LEISH-F1+MPL-SE adjuvant are given at Days 0, 28 and 56 for the LEISH-F1 and the MPL-SE adjuvant; meglumine antimoniate is administered daily, in cycles of 21 days, for the first 10 days followed by 11 days of rest.

The results reported in [65] are the following: At Day 84 of treatment, 80%, 50% and 38% of patients in Group 1, Group 2 and Group 3, respectively, had recovered (i.e. $P_1+P_2\to 0$). In our simulations, we represent the recovery rate by the relative difference between the parasite load without drug and the parasite load with drug at Day 84, that is, (17) $\text{Recovery}={\frac{(P_1+P_2)(\text{no drug})-(P_1+P_2)(\text{drug})}{(P_1+P_2)(\text{no drug})}|}_{\text{at Day 84}}(\times 100).$(17) Figure 6 shows our simulations of the experiments in [65]. In the simulations, we represented the effect of the MPL-SE adjuvant alone by a constant augmented secretion rate of IFN-γ, $c_{I_{\gamma }\text{adj}}$, in Equation (10(10) $\frac{\text{d}{I}_{\gamma }}{\text{d}t}=\underset{\text{secretion by}T}{\underbrace{{\lambda }_{T{I}_{\gamma }}(1+c_{I_{\gamma }\text{adj}})T(1+{\lambda }_{V{I}_{\gamma }}V)}}-\underset{\text{decay}}{\underbrace{{\mu }_{I_{\gamma }}{I}_{\gamma }}},$(10) ): we took $c_{I_{\gamma }\text{adj}}=0.35$. All other parameters in the drug equations are estimated, or fitted, in Section 5.

Figure 6. Treatment of leishmaniasis with MA (A), vaccine LEISH-F1+MPL-SE (V) and MPL-SE adjuvant as in [65]. The numbers in parentheses represent the recovery rates.

Figure 6 shows the profiles of $M_1+M_2$ and $P_1+P_2$ during 84 days. From Equation (17(17) $\text{Recovery}={\frac{(P_1+P_2)(\text{no drug})-(P_1+P_2)(\text{drug})}{(P_1+P_2)(\text{no drug})}|}_{\text{at Day 84}}(\times 100).$(17) ) we computed the recovery rates of Group 1 at 79.80%, Group 2 at 49.60%, and Group 3 at 37.83%, in agreement with the results of the clinical trials in [65].

In Figure 7, we show the profiles of the total parasite load under the same treatment protocol as in Figure 6, but with alternative doses of vaccine LEISH-F1, namely, 5 and 20 µ g. The simulations show that doubling the dose of LEISH-F1 significantly increases the recovery rate of the combination MA+LEISH-F1+MPL-SE.

Figure 7. Treatment of leishmaniasis with MA (A), vaccine LEISH-F1+MPL-SE (V) and MPL-SE adjuvant as in [65]. Profiles of total parasites load with alternative doses of V . The numbers in parentheses represent the recovery rates.

3.2.2. Sodium antimony gluconate and vaccine

We next proceed to simulate treatment with the combination of LEISH-F1+MPL-SE and sodium antimony gluconate (SAG), a first-choice compound for the treatment of leishmaniasis.

Figure 8 shows the simulations of various combinations of SAG, vaccine LEISH-F1+MPL-SE, and MPL-SE adjuvant with schedules similar to those in Figure 5. We note that our simulations, in Figure 8, with SAG only and in Figure 6 with MA only, show that the latter is more effective, which is in agreement with experimental results by Yesilova et al. [95]. We also observe a similar pattern as with MA in Figures 6 and 7, that doubling the dose of vaccine LEISH-F1+MPL-SE yields increased recovery rates.

Figure 8. Treatment of leishmaniasis with SAG, vaccine LEISH-F1+MPL-SE (V) and MPL-SE adjuvant. The numbers in parentheses represent the recovery rates.

3.3. Simulations with drugs under standard protocol

Meglumine antimoniate (MA) and sodium antimony gluconate (SAG) are first-line drugs in the treatment of leishmaniasis [28]. A standard protocol of treatment is injection of 20 mg/kg daily for 20 days (in some cases up to 28 days) [15]. Each of the drugs has its own set of adverse effects; see [27, 68] for MA, and [89, 92] for SAG.

In the following simulations, we combine the vaccine LEISH-F1+MPL-SE (V) with MA and evaluate the efficacy of the combinations at Day 84; drugs are administered for 20 days and the vaccine shots are given three times, twice or just once.

Figures 9 and 10 show the profiles of the total parasite loads and recovery rates of treatment with meglumine antimoniate in a single regimen and in combination with vaccine LEISH-F1+MPL-SE, where we give 3 vaccine shots (at Days 0, 28 and 56) in Figure 9, and 2 vaccine shots (at Days 0 and 28) in Figure 10(a), and 1 vaccine shot (at Day 0) in Figure 10(b). We see in all cases that combination therapy yields a much better recovery rate than chemotherapy alone, and doubling the dose of vaccine results in a significant increase in the recovery rate. Figure 10 also shows, in comparison with Figure 9, that reducing the number of vaccine shots decreases efficacy and accentuates disease resurgence irrespective of the vaccine dose. Similar simulations have been carried out with the chemotherapy SAG (not shown here).

Figure 9. Standard protocol treatment of leishmaniasis with MA and vaccine LEISH-F1+MPL-SE. The vaccine is given at Days 0 and 28 and 56, and the chemotherapy is administered at Days 0–20, daily. The numbers in parentheses are the recovery rates at Day 84.

Figure 10. Standard protocol treatment of leishmaniasis with MA and vaccine LEISH-F1+MPL-SE. The vaccine is given at (a) Days 0 and 28 and (b) Day 0; the chemotherapy is administered at Days 0–20, daily. The numbers in parentheses are the recovery rates at Day 84. (a) 2 vaccine dose and (b) 1 vaccine dose.

Figure 11 shows the effect on $(P_1+P_2)$ at day 84 of combining the drug MA with the LEISH-F1+MPL-SE vaccine (Figure 11(a)), and SAG with the LEISH-F1+MPL-SE vaccine (Figure 11(b)), for various doses of ${\gamma }_A,{\gamma }_S$ and ${\gamma }_V$. In both figures, increasing the doses ${\gamma }_A$ of MA and ${\gamma }_V$ of the vaccine, and of ${\gamma }_S$ of SAG and ${\gamma }_V$, simultaneously, yields a better reduction in parasite load. However, increasing ${\gamma }_A$ or ${\gamma }_S$ while keeping ${\gamma }_V$ fixed yields a much better reduction in parasite load than increasing ${\gamma }_V$ only.

Figure 11. Long-run total parasite load ($P_1+P_2$ at day 84), $P_1+P_2$, combination of (a) MA with LEISH-F1+MPL-SE and (b) SAG with LEISH-F1+MPL-SE. The parasite loads are in units of g/cm${}^3$.

We conclude this section by noting that the model system (1(1) $\frac{\mathrm{d}D}{\mathrm{d}t}=\underset{\mathrm{activation}}{\underbrace{{\lambda }_D{D}_0\frac{P_1+P_2}{(P_1+P_2)+K_P}}}-\underset{\mathrm{loss}\mathrm{due}\mathrm{to}\mathrm{death}}{\underbrace{{\mu }_{D}D}},$(1) )–(16(16) $\begin{array}{rl}\frac{\text{d}V}{\text{d}t}& =f_V(t)-\underset{\text{depletion through acting on}T}{\underbrace{{\mu }_{VT}TV}}-\underset{\text{degradation}}{\underbrace{{\mu }_{V}V}},\end{array}$(16) ) has a unique solution for all $0\le t<\mathrm{\infty }$. Indeed this follows for the following general theorem:

Theorem 3.1

Consider a system of differential equations (18) $\frac{\text{d}{x}_i}{\text{d}t}=f_i(x_1,\dots ,x_n),1\le i\le n,$(18) and assume that the $f_i$'s are continuously differentiable functions satisfying the following conditions in the nonnegative quadrant $x_1\ge 0,x_2\ge 0,\dots ,x_n\ge 0:$

1. $f_i(x_1,\dots ,x_n)<A+B\sum _{j=1}^n{x}_j$,

2. $f_i(x_1,\dots ,x_n)=g_i(x_1,\dots ,x_n)+h_i(x_1,\dots ,x_n)x_i$, where

3. $g_i(x_1,\dots ,x_n)\ge 0$ and $h_i(x_1,\dots ,x_n)\ge -A_1-B_1\sum _{j=1}^n{x}_j$, with some positive constants $A,B,A_1$ and $B_1$.

   If $x_i(0)>0$ for $1\le i\le n$, then

4. $0<x_i(t)<\frac{A}{B}+\hat{X}{e}^{nBt}$, for $1\le i\le n$ and all t>0, where $\hat{X}=\sum _{j=1}^n{x}_i(0)$.

Proof.

The inequalities in (iv) certainly hold if t is small. Hence, if the assertion (iv) is not true, then there is a smallest time τ such that (iv) holds for $t<\tau$ but not for $t=\tau$. To derive a contradiction, we get from Equation (16(16) $\begin{array}{rl}\frac{\text{d}V}{\text{d}t}& =f_V(t)-\underset{\text{depletion through acting on}T}{\underbrace{{\mu }_{VT}TV}}-\underset{\text{degradation}}{\underbrace{{\mu }_{V}V}},\end{array}$(16) ) and (i), the inequality $\frac{\text{d}z}{\text{d}t}\le nA+nBz,\text{for}t<\tau ,$where $z=\sum _{j=1}^n{x}_j$. Hence (19) $\begin{array}{rl}z(\tau )& \le e^{nB\tau }z(0)+\frac{A}{B}(1-e^{-nB\tau })\end{array}$(19) (20) $\begin{array}{rl}& <\frac{A}{B}+\hat{X}{e}^{nB\tau },\end{array}$(20) so that the second inequality in (iv) follows with $t=\tau$.

We next use the conditions (ii) and (iii) to get $\frac{\text{d}{x}_i}{\text{d}t}\ge -(A_1+B_1\sum _{j=1}^n{x}_j)x_i$and, by the second inequality in (iv), $A_1+B_1\sum _{j=1}^n{x}_j(t)\le A_1+n{B}_1(\frac{A}{B}+\hat{X}{e}^{nB\tau })\equiv C,\text{for}t\le \tau .$Hence $\frac{\text{d}{x}_i}{\text{d}t}\ge -C{x}_i,\text{or}{x}_i(\tau )>e^{-C\tau }{x}_i(0)>0.$This completes the proof of (iv) by contradiction.

As easily seen, the system (1(1) $\frac{\mathrm{d}D}{\mathrm{d}t}=\underset{\mathrm{activation}}{\underbrace{{\lambda }_D{D}_0\frac{P_1+P_2}{(P_1+P_2)+K_P}}}-\underset{\mathrm{loss}\mathrm{due}\mathrm{to}\mathrm{death}}{\underbrace{{\mu }_{D}D}},$(1) )–(16(16) $\begin{array}{rl}\frac{\text{d}V}{\text{d}t}& =f_V(t)-\underset{\text{depletion through acting on}T}{\underbrace{{\mu }_{VT}TV}}-\underset{\text{degradation}}{\underbrace{{\mu }_{V}V}},\end{array}$(16) ) satisfies the conditions (i)–(iii) of Theorem 3.1. Hence, the system has a unique solution for all t>0, and each variable satisfies the inequality in (iv).

4. Conclusion and discussion

Leishmaniasis is a disease caused by a parasite from one of 20 Leishmania species. The most common forms of the disease are cutaneous leishmaniasis (CL) and visceral leishmaniasis (VL); the latter is likely to be deadly if not treated. The current control strategy for leishmaniasis relies on chemotherapy drugs such as sodium antimony gluconate (SAG) or meglumine antimoniate (MA). However, all these chemotherapy compounds have poor efficacy, and they are associated with toxicity and other adverse effects, as well as drug resistance.

A number of experimental vaccines have been developed, although no vaccine is currently available. However, some studies show that vaccines such as LEISH-F1+MPL-SE (V) have some efficacy as drugs for CL patients. This suggests that treatment of leishmaniasis can be improved by combining such a vaccine with one of the chemotherapy agents. Indeed, clinical trials, reported in [65], show significant improvement when V was combined with MA.

In this paper, we study, by a mathematical model, the efficacy of various combinations of V with MA, or with SAG. The model is an extension of [79, 80] that enables the inclusion of these drugs. We first demonstrated that the model predicts, what is well known, that even with no treatment, some patients recover, or partially recover; this we show by modifying some ‘personal’ parameters in the model equations (Figures 2–4).

We next repeated with simulations the clinical trials in [65] with MA and MA$+V$, and showed, in Figures 6 and 7, that the efficacy of treatments is in agreement with the percentage of recovery reported in [65]. In Figure 8 we showed similar results in simulating clinical trials with SAG and SAG$+V$.

We next considered the current standard treatment with MA, where the drug is administered daily for 20 days [15]. In order to improve efficacy we added to this treatment 3 injections of V in the same schedule used in clinical trials [16, 65], namely: at Days 0, 28, 56, and in the same amounts. In Figure 9 we found that, by Day 84, the efficacy of treatment with MA alone was just 29%, but in combination with vaccine V , the efficacy increased significantly, to 65%, 78% and 91% when the amount of injection doubled from 5 to 10, and from 10 to 20 µ g. Furthermore, these recovery rates already occurred by Day 42. Similar results were obtained with SAG$+V$.

We next asked the question of whether we could get such improvements in MA$+V$ if we reduce the number of vaccine injections, administering it at Days 0 and 28, or just at Day 0. We found, as in Figure 10(a), that with 2 shots the efficacy is still high (with 62%, 76% and 89%) but the density of the parasites begins to grow after the first 40–50 days, which indicates eventual resurgence of leishmaniasis. With one vaccine shot, given at Day 0, we see, in Figure 10(b), that the profiles for the total parasite load are qualitatively similar to those in Figure 10(a), but with smaller efficacy (58%, 71% and 85%) and higher tendency toward disease resurgence.

The model has as many variables (species) as were needed in order to be able to evaluate the efficacy of drugs and vaccines. This, necessarily, resulted in introducing many unknown parameters. We estimated some of these parameters making assumptions, such as ‘steady-state’ in some of the dynamical equations. Other parameters were ‘fitted’ using a few reported clinical trials. We addressed this shortcoming of the model by performing sensitivity analysis, showing that some randomized changes in parameters do not qualitatively change model's predictions.

In this paper, we demonstrated that by combining vaccine LEISH-F1+MPL-SE with a standard treatment with chemotherapy drug, such as meglumine antimoniate, we can improve very significantly the efficacy of treatment. But these conclusions will need to be confirmed in actual clinical trials, with attention also to potential side effects.

5. Parameter estimations

We denote by $X^0$ the average density/concentration of species X in patients before the beginning of therapy and take $K_X=X^0\text{for all species.}$

5.1. Estimate for $K_T$

The level of T cells (CD4${}^{+}$ T cells+CD8${}^{+}$ T cells) was measured in non-relapsing and relapsing VL patients immediately after treatment and 6–12 months post-treatment [49]. We take the 12-month post-treatment relapsing VL case to correspond to an active VL infection: $T_4=800{\text{cells/mm}}^3(\text{Figure 2(A)}),T_8=500{\text{cells/mm}}^3(\text{Figure 2(C)}).$Since the mass of one cell is approximately $5\times {10}^{-10}$ g, we take $T^0=K_T=(800\times {10}^3+500\times {10}^3)\times 5\times {10}^{-10}=6.5\times {10}^{-4}{\text{g/cm}}^3.$

5.2. Estimates for $M_1^0$ and $M_2^0$

The level of macrophages was estimated in [79, 80] as $K_M=1.5\times {10}^{-3}$ g/cm${}^3$. As in [79, 80], we take the ratio $M_1/M_2=4$, so that $M_1^0=1.2\times {10}^{-3}{\text{g/cm}}^3,M_2^0=3\times {10}^{-4}{\text{g/cm}}^3.$

5.3. Estimates for $P_1^0,P_2^0$ and $K_P$

We have $K_P=\frac{N(M_1^0+M_2^0)}{150}=5\times {10}^{-4}{\text{g/cm}}^3.$We have the ratio $P_2/P_1=N_2/N_1=10$, as in [79, 80]. Hence, $P_1^0=4.55\times {10}^{-5}{\text{g/cm}}^3,P_2^0=4.55\times {10}^{-4}{\text{g/cm}}^3.$

5.4. Estimate for $I_2^0,I_{10}^0,I_{12}^0,I_{\gamma }^0,N{O}^0,T_{\alpha }^0$ and $T_{\beta }^0$

The level of circulating TNF-α in CL patients without treatment is $21.83\pm 4.92$ pg/ml [48] (Table 1, Figure 1). We assume that it is the same for VL patients without treatment so that $T_{\alpha }^0=21.83\times {10}^{-12}{\text{g/cm}}^3.$In [26], the lymphoproliferative response level to drug-sensitive and drug resistance VL patients was measured for TNF-α, IFN-γ, IL-10 and TGF-β; Figs. 1–4 in [26] show that $I_{\gamma }\simeq 2T_{\alpha },I_{10}\simeq 0.07T_{\alpha }$and $T_{\beta }$ is significantly higher than $I_{10}$. But, as pointed out in [26], the high level of $T_{\beta }$ was mostly associated with treatment failure. Assuming that the ratios of $I_{\gamma }$ and $I_{10}$ to $T_{\alpha }$ remain the same in the infected tissue, and that $T_{\beta }\simeq I_{10}$, we get $I_{\gamma }^0=4.37\times {10}^{-11}{\text{g/cm}}^3,T_{\beta }^0=1.53\times {10}^{-12}{\text{g/cm}}^3,I_{10}^0=1.53\times {10}^{-12}{\text{g/cm}}^3.$In tuberculosis patients, IL-2 concentration before therapy was 164.5 pg/ml and after 2-month therapy, it was 92.11 pg/ml [8], while the concentration of IL-12 in tuberculosis patients was 11.86 pg/ml compared to 5.1 pg/ml for healthy individuals [3]. We take $I_2^0=140\text{pg/ml}=14\times {10}^{-11}{\text{g/cm}}^3,I_{12}^0=9\times {10}^{-12}{\text{g/cm}}^3.$The density of NO is 1.34 g/mole. The serum concentration in health is $25\text{μ}\text{Mol}/L=25\times 1.34\times {10}^{-9}{\text{g/cm}}^3$ [30]; in sickle cell disease, it increases to $59\times 1.34\times {10}^{-9}$ g/cm${}^3$ [4]. We take $N{O}^0=40\times 1.34\times {10}^{-9}=5.36\times {10}^{-8}{\text{g/cm}}^3.$For each cytokine X, we take $K_X=X^0$.

5.5. Estimate for ${\mu }_{NO}$

The half-life of nitric oxide is approximately 445 seconds (0.00515 days) [36]. Hence, ${\mu }_{NO}=\frac{\ln 2}{0.00515}=134.6{\text{d}}^{-1}.$

5.6. Estimate for ${\mu }_A$

The half-life of meglumine antimoniate is approximately 3.42 hours (0.14 days) [35] (Figure 2(a)). Hence, ${\mu }_A=\frac{\ln 2}{0.14}=4.95{\text{d}}^{-1}.$

5.7. Estimate for ${\mu }_S$

The half-life of SAG is approximately 10 hours (0.42 days) [57]. Hence, ${\mu }_S=\frac{\ln 2}{0.42}=1.65{\text{d}}^{-1}.$

5.8. Estimate for ${\gamma }_A$

Meglumine antimoniate, in standard treatment, is administered by intravenous or injection route with a usual dose of 20 milligrams (mg) per kilogram (kg) of body weight per day [12]; in a clinical trial by Nascimentao et al. [65], meglumine antimoniate was given at 10 mg/kg per day. We estimate the dose to be per 1000 cm${}^3$. We accordingly take $\begin{array}{rl}{\gamma }_A& =20\text{mg}/1000{\text{cm}}^3{\text{d}}^{-1}=2\times {10}^{-5}{\text{g/cm}}^3{\text{d}}^{-1}\text{in standard treatment, and}\\ {\gamma }_A& =1\times {10}^{-5}{\text{g/cm}}^3{\text{d}}^{-1}\text{in clinical trials}.\end{array}$

5.9. Estimate for ${\gamma }_S$

SAG (or sodium stibogluconate) is administered by injectable solutions, intravenously or intramuscular, at a dose of 20 mg/kg per day [11]. We accordingly take, in standard treatment, ${\gamma }_S=20\text{mg}/1000{\text{cm}}^3{\text{d}}^{-1}=2\times {10}^{-5}{\text{g/cm}}^3{\text{d}}^{-1},$and ${\gamma }_S=1\times {10}^{-5}{\text{g/cm}}^3{\text{d}}^{-1}\text{in our clinical trials}.$

5.10. Estimate for ${\gamma }_V$

Chakravarty et al. [16] and Nascimento et al. [65] performed clinical trials with LEISH-F1+MPL-SE vaccine in the prevention of VL and suggested the following doses: 510and 20 µ g of LEISH-F1 antigen + 25 µ g MPL-SE adjuvant. We assume that the vaccine works through the immune system more effectively than chemotherapy drugs, and take, in the case of 10 µ g of LEISH-F1 + 25 µ g of MPL-SE adjuvant, ${\gamma }_V=2\times {10}^{-5}{\text{g/cm}}^3{\text{d}}^{-1}.$

5.11. Estimate for ${\lambda }_{M_2{M}_1{T}_{\alpha }}$

We assume that the rate of transition for $M_2$ to $M_1$ by $T_{\alpha }$ is the same as by $I_{\gamma }$ and take, as in [79], ${\lambda }_{M_2{M}_1{I}_{\gamma }}=0.69{\text{d}}^{-1}.$

5.12. Estimates for ${\lambda }_{I_{\gamma }{P}_1},{\lambda }_{I_{\gamma }{P}_2},{\lambda }_{T_{\alpha }{P}_1}$ and ${\lambda }_{T_{\alpha }{P}_2}$

We assume that the strength of IFN-γ and TNF-α in augmenting the death of parasites is the same as that by NO, and take, as in [79], ${\lambda }_{I_{\gamma }{P}_1}={\lambda }_{I_{\gamma }{P}_2}={\lambda }_{T_{\alpha }{P}_1}={\lambda }_{T_{\alpha }{P}_2}=\mathrm{0.062.}$

5.13. Estimate by equations (without drugs)

• Equation (1(1) $\frac{\mathrm{d}D}{\mathrm{d}t}=\underset{\mathrm{activation}}{\underbrace{{\lambda }_D{D}_0\frac{P_1+P_2}{(P_1+P_2)+K_P}}}-\underset{\mathrm{loss}\mathrm{due}\mathrm{to}\mathrm{death}}{\underbrace{{\mu }_{D}D}},$(1) ): By the steady-state equation ${\lambda }_D{D}_0/2-{\mu }_{D}D=0,$where ${\mu }_D=0.1$ d${}^{-1}$, $D=D^0=7\times {10}^{-4}$ g/cm${}^3$ and $D_0=7\times {10}^{-4}$ g/cm${}^3$; hence, ${\lambda }_D=0.2{\text{d}}^{-1}.$

• Equation (6(6) $\begin{array}{rl}\frac{\text{d}T}{\text{d}t}& =\underset{I_2\text{-induced proliferation of T cells}}{\underbrace{{\lambda }_{T}T\frac{I_2}{I_2+K_{I_2}}}}+\underset{\text{naive T cells activation}}{\underbrace{{\lambda }_{DT}{T}_0\frac{D}{D+K_D}\frac{I_{12}}{I_{12}+K_{I_{12}}}}}\\ & +\underset{I_{12}\text{-induced proliferation of T cells}}{\underbrace{{\lambda }_{MT}{T}_0\frac{M_1}{M_1+K_M}\frac{I_{12}}{I_{12}+K_{I_{12}}}\frac{1}{1+I_{10}/K_{I_{10}}}}}+\underset{\text{to become memory T cells}}{\underbrace{{\lambda }_{VT}TV}}-\underset{\text{loss due to death}}{\underbrace{{\mu }_{T}T}}.\end{array}$(6) ): We assume that ${\lambda }_{DT}={\lambda }_{MT}=2{\lambda }_T$, and using the steady-state equation ${\lambda }_{T}T/2+{\lambda }_{DT}{T}_0/4+{\lambda }_{MT}{T}_0/8-{\mu }_{T}T=0,$where ${\mu }_T=0.197$ d${}^{-1}$, $T^0=6.5\times {10}^{-4}$ g/cm${}^3$ and $T_0=0.9$ g/cm${}^3$, we get, ${\lambda }_T=1.9\times {10}^{-4}{\text{d}}^{-1},{\lambda }_{MT}=3.8\times {10}^{-4}{\text{d}}^{-1},{\lambda }_{DT}=3.8\times {10}^{-4}{\text{d}}^{-1}.$We combine ${\lambda }_{DT}$ and $T_0$ and write ${\lambda }_{DT{T}_0}=3.42\times {10}^{-4}$ g/cm${}^3$ d${}^{-1}$.

• Equation (7(7) $\frac{\text{d}{I}_2}{\text{d}t}=\underset{\text{secretion by}T}{\underbrace{{\lambda }_{T{I}_2}T(1+{\lambda }_{S{I}_2}S+{\lambda }_{V{I}_2}V)}}-\underset{\text{decay}}{\underbrace{{\mu }_{I_2}{I}_2}}.$(7) ): We use the steady-state equation ${\lambda }_{T{I}_2}T-{\mu }_{I_2}{I}_2=0,$where ${\mu }_{I_2}=330$ d${}^{-1}$, $T^0=6.5\times {10}^{-4}$ g/cm${}^3$ and $K_{I_2}=14\times {10}^{-11}$ g/cm${}^3$; hence, ${\lambda }_{T{I}_2}=7.1\times {10}^{-5}{\text{d}}^{-1}.$

• Equation (8(8) $\frac{\text{d}{I}_{10}}{\text{d}t}=\underset{\text{secretion by}{M}_2}{\underbrace{{\lambda }_{M_2{I}_{10}}{M}_2(1+{\lambda }_{A{I}_{10}}A)}}-\underset{\text{decay}}{\underbrace{{\mu }_{I_{10}}{I}_{10}}}.$(8) ): Using the steady-state equation (without drug) ${\lambda }_{M_2{I}_{10}}{M}_2-{\mu }_{I_{10}}{I}_{10}=0,$where ${\mu }_{I_{10}}=14.5$ d${}^{-1}$, $M_2^0=3\times {10}^{-4}$ g/cm${}^3$ and $K_{I_{10}}=1.53\times {10}^{-12}$ g/cm${}^3$, we get ${\lambda }_{M_2{I}_{10}}=7.4\times {10}^{-8}{\text{d}}^{-1}.$

• Equation (9(9) $\frac{\text{d}{I}_{12}}{\text{d}t}=\underset{\text{secretion by}{M}_1}{\underbrace{{\lambda }_{M_1{I}_{12}}{M}_1\frac{1}{1+I_{10}/K_{I_{10}}}(1+{\lambda }_{S{I}_{12}}S)}}-\underset{\text{decay}}{\underbrace{{\mu }_{I_{12}}{I}_{12}}}.$(9) ): We use the steady-state equation (without drug) ${\lambda }_{M_1{I}_{12}}{M}_1/2-{\mu }_{I_{12}}{I}_{12}=0,$where ${\mu }_{I_{12}}=1.39$ d${}^{-1}$, $M_1^0=1.2\times {10}^{-3}$ g/cm${}^3$ and $K_{I_{12}}=9\times {10}^{-12}$ g/cm${}^3$; hence, ${\lambda }_{M_1{I}_{12}}=2.1\times {10}^{-8}{\text{d}}^{-1}.$

• Equation (10(10) $\frac{\text{d}{I}_{\gamma }}{\text{d}t}=\underset{\text{secretion by}T}{\underbrace{{\lambda }_{T{I}_{\gamma }}(1+c_{I_{\gamma }\text{adj}})T(1+{\lambda }_{V{I}_{\gamma }}V)}}-\underset{\text{decay}}{\underbrace{{\mu }_{I_{\gamma }}{I}_{\gamma }}},$(10) ): By the steady-state equation ${\lambda }_{T{I}_{\gamma }}T-{\mu }_{I_{\gamma }}{I}_{\gamma }=0,$where ${\mu }_{I_{\gamma }}=33$ d${}^{-1}$, $T^0=6.5\times {10}^{-4}$ g/cm${}^3$ and $K_{I_{\gamma }}=4.37\times {10}^{-11}$ g/cm${}^3$, we get ${\lambda }_{T{I}_{\gamma }}=2.22\times {10}^{-6}{\text{d}}^{-1}.$

• Equation (11(11) $\frac{\text{d}NO}{\text{d}t}=\underset{\text{secretion by}{M}_1\text{and}{M}_2}{\underbrace{{\lambda }_{MNO}(M_1+M_2)(1+{\lambda }_{I_{\gamma }NO}\frac{I_{\gamma }}{I_{\gamma }+K_{I_{\gamma }}}+{\lambda }_{SNO}S+{\lambda }_{ANO}A)}}-\underset{\text{decay}}{\underbrace{{\mu }_{NO}NO}}.$(11) ): From the steady-state equation ${\lambda }_{MNO}(M_1+M_2)(1+{\lambda }_{I_{\gamma }NO}/2)-{\mu }_{NO}NO=0,$with ${\lambda }_{I_{\gamma }NO}=2$, ${\mu }_{NO}=134.6$ d${}^{-1}$, $K_{NO}=5.36\times {10}^{-8}$ g/cm${}^3$, $M_1^0=1.2\times {10}^{-3}$ g/cm${}^3$ and $M_2^0=3\times {10}^{-4}$ g/cm${}^3$, we find that ${\lambda }_{MNO}=2.4\times {10}^{-3}{\text{d}}^{-1}.$

• Equation (12(12) $\frac{\text{d}{T}_{\alpha }}{\text{d}t}=\underset{\text{secretion by}{M}_1\text{augmented by}V\text{and}A}{\underbrace{{\lambda }_{M_1{T}_{\alpha }}{M}_1(1+{\lambda }_{V{T}_{\alpha }}V+{\lambda }_{A{T}_{\alpha }}A)}}+\underset{\text{secretion by}T}{\underbrace{{\lambda }_{T{T}_{\alpha }}T}}-\underset{\text{decay}}{\underbrace{{\mu }_{T_{\alpha }}{T}_{\alpha }}}.$(12) ): We assume that ${\lambda }_{M_1{T}_{\alpha }}={\lambda }_{T{T}_{\alpha }}$ and use the steady-state equation ${\lambda }_{M_1{T}_{\alpha }}{M}_1+{\lambda }_{T{T}_{\alpha }}T-{\mu }_{T_{\alpha }}{T}_{\alpha }=0,$where ${\mu }_{T_{\alpha }}=199$ d${}^{-1}$, $K_{T_{\alpha }}=21.83\times {10}^{-12}$ g/cm${}^3$, $M_1^0=1.2\times {10}^{-3}$ g/cm${}^3$ and $T^0=6.5\times {10}^{-4}$ g/cm${}^3$; hence, ${\lambda }_{M_1{T}_{\alpha }}=2.35\times {10}^{-6}{\text{d}}^{-1}\text{and}{\lambda }_{T{T}_{\alpha }}=2.35\times {10}^{-6}{\text{d}}^{-1}.$

• Equation (13(13) $\frac{\text{d}{T}_{\beta }}{\text{d}t}=\underset{\text{secretion by}{M}_2}{\underbrace{{\lambda }_{M_2{T}_{\beta }}{M}_2}}-\underset{\text{decay}}{\underbrace{{\mu }_{T_{\beta }}{T}_{\beta }}}.$(13) ): Using the steady-state equation ${\lambda }_{M_2{T}_{\beta }}{M}_2-{\mu }_{T_{\beta }}{T}_{\beta }=0,$where ${\mu }_{T_{\beta }}=399.25$ d${}^{-1}$, $K_{T_{\beta }}=1.53\times {10}^{-12}$ g/cm${}^3$ and $M_2^0=3\times {10}^{-4}$ g/cm${}^3$, we get ${\lambda }_{M_2{T}_{\beta }}=2.04\times {10}^{-6}{\text{d}}^{-1}.$

6. Parameter sensitivity analysis

We performed a sensitivity analysis with respect to the total parasite load at Day 84, $(P_1+P_2)(84)$, without self-healing with no drugs, for the transition/activation parameters ${\lambda }_{M_2{M}_1{I}_{\gamma }}$, ${\lambda }_{M_2{M}_1{T}_{\alpha }}$, ${\lambda }_{M_1{M}_2{I}_{10}}$, ${\lambda }_{M_1{M}_2{T}_{\beta }}$, ${\lambda }_{NO{P}_1}$, ${\lambda }_{NO{P}_2}$, ${\lambda }_{I_{\gamma }{P}_1}$, ${\lambda }_{I_{\gamma }{P}_2}$, ${\lambda }_{T_{\alpha }{P}_1}$ and ${\lambda }_{T_{\alpha }{P}_2}$, the bursting rates ${\tilde{\mu }}_{M_1}$, ${\tilde{\mu }}_{M_2}$, and the fraction of parasites phagocyted by M1 macrophages after burst, θ (Figure 12).

Figure 12. Parameter sensitivity analysis for the total parasite load after 84 days. The p-values are less than .01.

The computations were done using Latin Hypercube Sampling/Partial Rank Correlation Coefficient (LHS/PRCC) with a Matlab package by [47, 60]. The range for the parameters in the sensitivity analysis was between $\pm 50\mathrm{\%}$ of their baseline values in Tables 2 and 3.

Table 3. Descriptions, values, and sources of the model parameters.

Parameter	Description	Value/range	Reference
${\lambda }_{NO{P}_1}$	Strength of NO-effect on parasites $P_1$	0.062	Fitted
${\lambda }_{NO{P}_2}$	Strength of NO-effect on parasites $P_2$	0.56	Fitted
${\lambda }_{I_{\gamma }NO}$	Strength of $I_{\gamma }$-effect on NO production by $M_1$	2	Est.
${\lambda }_{I_{\gamma }{P}_1}$	Strength of $I_{\gamma }$-effect on parasites $P_1$	0.062	Fitted
${\lambda }_{I_{\gamma }{P}_2}$	Strength of $I_{\gamma }$-effect on parasites $P_2$	0.062	Fitted
${\lambda }_{T_{\alpha }{P}_1}$	Strength of $T_{\alpha }$-effect on parasites $P_1$	0.062	Fitted
${\lambda }_{T_{\alpha }{P}_2}$	Strength of $T_{\alpha }$-effect on parasites $P_2$	0.062	Fitted
${\lambda }_{A{I}_{10}}$	Effect of SAG on IL-10 secretion	$5\times {10}^{10}$	Fitted
${\lambda }_{S{I}_2}$	Effect of SAG on IL-2 secretion	$2\times {10}^4$	Fitted
${\lambda }_{V{I}_2}$	Effect of LEISH-F1+MPL-SE on IL-2 secretion	${10}^{10}$ cm${}^3$/g	Fitted
${\lambda }_{S{I}_{12}}$	Effect of SAG on IL-12 secretion	$4\times {10}^2$	Fitted
${\lambda }_{V{I}_{\gamma }}$	Effect of LEISH-F1+MPL-SE on IFN-γ secretion	${10}^{10}$ cm${}^3$/g	Fitted
${\lambda }_{A{T}_{\alpha }}$	Effect of meglumine antimoniate on TNF-α secretion	${10}^{10}$ cm${}^3$/g	Fitted
${\lambda }_{V{T}_{\alpha }}$	Effect of LEISH-F1+MPL-SE on TNF-α secretion	${10}^{10}$ cm${}^3$/g	Fitted
${\lambda }_{VT}$	Transition rate of T cells to memory cells	$8\times {10}^3$ cm${}^3$/g d${}^{-1}$	Fitted
${\lambda }_{SNO}$	Effect of SAG on NO secretion	$1.7\times {10}^3$	Fitted
${\lambda }_{ANO}$	Effect of meglumine antimoniate on NO secretion	${10}^7$ cm${}^3$/g	Fitted
${\lambda }_{M_1}$	Basal production rate of M1 macrophages	$4.94\times {10}^{-5}$ g/cm${}^3$ d${}^{-1}$	Fitted
${\lambda }_{M_2}$	Basal production rate of M2 macrophages	$2.625\times {10}^{-5}$ g/cm${}^3$ d${}^{-1}$	Fitted
${\lambda }_{P_1}$	Intrinsic growth rate of parasites within $M_1$	2.21 d${}^{-1}$	[18, 29, 79]
${\lambda }_{P_2}$	Intrinsic growth rate of parasites within $M_2$	3.82 d${}^{-1}$	[18, 29, 79]
${\lambda }_T$	$I_2$-induced rate of proliferation of T cells	$1.9\times {10}^{-4}$ d${}^{-1}$	Est.
${\lambda }_D$	Activation rate of DCs	0.66 d${}^{-1}$	Est.
${\lambda }_{DT{T}_0}$	Activation rate of T cells by naive $T_0$ cells	$3.42\times {10}^{-4}$ g/cm${}^3$ d${}^{-1}$	Est.
${\lambda }_{MT}$	Activation rate of T cells by M1 macrophages	$3.8\times {10}^{-4}$ d${}^{-1}$	Est.
${\beta }_A,{\beta }_S,{\beta }_V$	PK exponents of A, S and V	0.016, 0.01, 0.01	Fitted
θ	Fraction of parasites phagocyted by $M_1$ after burst	0.8	Fitted
N	Average number of parasites for macrophage burst	50 (30–50)	[18]
$N_1$	Average number of parasites in M1 macrophages	2	[18]
$N_2$	Average number of parasites in M2 macrophages	20	[18]
$N_{\ast }$	Maximum number of parasites in macrophages	50	[18, 79]
$c_{I_{\gamma }\text{adj}}$	Augmented rate of IFN-γ secretion by MPL-SE adjuvant	0.35	Fitted
${\mu }_D$	Death rate of DCs	0.1 d${}^{-1}$	[39]
${\mu }_{NO}$	Decay rate of NO	134.6 d${}^{-1}$	[36]Est.
${\mu }_{M_1}$	Death rate of M1 macrophages	0.056 d${}^{-1}$	[79, 91]
${\mu }_{M_2}$	Death rate of M2 macrophages	0.014 d${}^{-1}$	[79, 91]
${\mu }_{P_1}$	Death rate of $P_1$ parasites	1.85 d${}^{-1}$	[33, 79, 87]
${\mu }_{P_2}$	Death rate of $P_2$ parasites	2.22 d${}^{-1}$	[33, 79, 87]
${\mu }_T$	Death rate of T cells	0.197 d${}^{-1}$	[20, 40, 93]
${\mu }_{I_2}$	Decay rate of IL-2	330 d${}^{-1}$	[24]
${\mu }_{I_{10}}$	Decay rate of IL-10	14.5 d${}^{-1}$	[45, 51, 63]
${\mu }_{I_{12}}$	Decay rate of IL-12	1.39 d${}^{-1}$	[5]
${\mu }_{I_{\gamma }}$	Decay rate of IFN-γ	33 d${}^{-1}$	[44]
${\mu }_{T_{\alpha }}$	Decay rate of TNF-α	199 d${}^{-1}$	[77, 81]
${\mu }_{T_{\beta }}$	Decay rate of TGF-β	399.25 d${}^{-1}$	[78, 85]
${\mu }_A$	Decay rate of meglumine antimoniate	4.95 d${}^{-1}$	[35]Est.
${\mu }_S$	Decay rate of SAG	1.65 d${}^{-1}$	[57]Est.
${\mu }_V$	Decay rate of vaccine LEISH-F1+MPL-SE	${10}^{-3}$ d${}^{-1}$	Fitted
${\tilde{\mu }}_{M_2}$	Bursting rate of $M_2$	0.17 d${}^{-1}$	[79, 83]
${\tilde{\mu }}_{M_1}$	Bursting rate of $M_1$	0.017 d${}^{-1}$	[72, 79]
${\mu }_{MA}$	Depletion rate of meglumine antimoniate by M1 macrophages	${10}^{1}0$ cm${}^3$/g d${}^{-}1$	Fitted
${\mu }_{SM}$	Depletion rate of SAG by M1 and M2	1 cm${}^3$/g d${}^{-1}$	Fitted
${\mu }_{ST}$	Depletion rate of SAG by T cells	1 cm${}^3$/g d${}^{-1}$	Fitted
${\mu }_{VT}$	Depletion rate of vaccine LEISH-F1+MPL-SE by T cells	100 cm${}^3$/g d${}^{-}1$	Fitted
$D_0$	Density of DC in healthy	$0.7\times {10}^{-3}$ g/cm${}^3$	[43, 79]

Expectedly, increasing any one of the strength of NO-, $I_{\gamma }$- or $T_{\alpha }$-effect on parasites $P_1$ and $P_2$, namely ${\lambda }_{NO{P}_1},{\lambda }_{NO{P}_2},{\lambda }_{I_{\gamma }{P}_1},{\lambda }_{I_{\gamma }{P}_2},{\lambda }_{T_{\alpha }{P}_1}$ and ${\lambda }_{T_{\alpha }{P}_2}$ result in decreased total parasite load. It is worth noting that the NO-effect on parasites is the most negatively correlated with respect to $P_1+P_2$.

Increased transition from M1 to M2 macrophages increases $P_1+P_2$, since the burst rate of $M_1$ is smaller than the burst rate of $M_2$, which enables additional growth time of the parasites in $M_1$. This explains why ${\lambda }_{M_2{M}_1{I}_{\gamma }}$ and ${\lambda }_{M_2{M}_2{T}_{\alpha }}$ are positively correlated, while ${\lambda }_{M_1{M}_2{I}_{10}}$ and ${\lambda }_{M_1{M}_2{T}_{\beta }}$ are negatively correlated.

Data availability statement

All relevant data are within the manuscript and its Supporting information files.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Funding

The author(s) reported there is no funding associated with the work featured in this article.