Bernstein collocation method for neutral type functional differential equation

1.   Introduction

The immune system response to pathogens or abnormal exogenous stimuli is a very complex and well-orchestrated process whose aim is to keep the host safe and healthy. It is composed of two major lines of defense: innate and adaptive immune responses, to counter any and all threats. Innate immunity, which is the first line of defense, is made up of external defenses, such as the skin and mucous membranes for inner surfaces, while phagocytic cells such as neutrophils, macrophages, and natural killer cells are components of internal defenses. The innate response to an infection is immediate and passive. Adaptive immunity only kicks in when the innate response is unsuccessful in getting rid of the pathogen. Though it takes a longer time to get activated, it targets the pathogen more accurately. The adaptive response is antigen- specific, well diversified and has memory of past responses ^[1,2,3].

Despite its sophistication, some viruses, bacteria, and parasites are able to breach or escape the defenses and infect the host, amongst which are the protozoa that causes leishmaniasis ^[4]. Leishmania infects about 2 million people annually worldwide and is transmitted through the bite of infected sandfly. The transmission cycle begins when female sandflies in need of protein to develop eggs ingest the Leishmania parasites in the course of taking blood meals from infected hosts. The parasites incubate inside female sandflies, alternating through multiple developmental stages and in the process infecting them. Infected female sandflies transmit the parasite when taking fresh blood meals from naive hosts, thereby completing the transmission cycle ^[5,6,7].

The mechanism of Leishmania infection has been carefully considered and the parasite’s survival depends on the availability of macrophages ^[8,9,10]. It is not well understood how the parasite enters inside macrophages. Some studies suggest that Leishmania are able first to get shelter within neutrophils which usually arrive at the infection site before macrophages ^[11,12,13]. Once inside the macrophages, the parasite differentiates into a smaller immobile form known as an amastigote. In order to survive and replicate, the amastigote both manipulates and disables the macrophage’s ability to destroy it or demonstrate signs of cellular infection, consequently avoiding being destroyed by other components of the immune system. Through maneuvering the amastigote is also able to prolong the lifespan of the macrophage giving itself enough time to replicate while simultaneously feeding on molecules within the macrophage till the host cell is eventually destroyed releasing many more amastigotes.

To successfully evade destruction by the immune system and to proliferate, it is crucial for the amastigote to induce an overexpression of genes that codes the surface of the host macrophage such as receptors, antibodies and cytokines ^[14,15,16]. An Fc receptor is a protein found on the surface of certain immune system cells including macrophages, neutrophils and natural killer cells. They play a critical role in the protective functions of the immune system by enabling phagocytic or cytotoxic cells to bind to antibodies that are attached to the surface of microbes or microbe infected cells, ear-marking them for destruction. Fc gamma receptors (Fc$\gamma$R), divided into the subclasses Fc$\gamma$RI, -Ⅱ, -Ⅲ and -Ⅳ are the most important for stimulating phagocytosis of marked cells. Fc$\gamma$Rs are activated by immunoglobulin G (IgG), the main type of antibody for controlling infection and made up of four subclasses IgG1, 2, 3, and 4, with IgG1 being the most abundant. For example, the activation of Fc$\gamma$RIII by IgG causes the release of interferon gamma (IFN$\gamma$), a cytokine that is an important activator of macrophages and is critical for the immune system response against an infection. IFN$\gamma$ is a small protein secreted by an immune cell that signals to other immune cells it is harboring an antigen. Unfortunately, the synthesis of IFN$\gamma$ along with others can be inhibited by another cytokine, IL-10, that is primarily produced by monocytes. IL-10 has a strong ability to suppress the antigen-presentation capacity of antigen presenting cells like macrophages and consequently leads to an infection establishing itself. Additionally, IL-10 can down-regulate nitric oxide production in infected cells, thus preventing parasite death ^[17,18,19].

Murine studies for different Leishmania species have been attributed to disease susceptibility and progression due to the parasite's ability to interfere with the normal production of cytokines ^{[20,21,22,23]}. These studies that aim at delineating the pathways of chronic infection utilize mice in which genes responsible for producing specific immune cell proteins have been knocked out (KO). Buxbaum et al. ^[24] observed that mice deficient in $\beta 2$-microglobulin ($\beta 2$m) (a protein that is found on the surface of almost all nucleated cells), or the cytokine IL-10-, or Fc$\gamma$R were able to resolve infection with Leishmania mexicana. Further investigation of the Fc$\gamma$RIII subclass by Thomas et al. ^[25], demonstrated that Fc$\gamma$RIII KO mice had a stronger immune response for resisting and resolving cutaneous Leishmania mexicana lesions. The implication being that Leishmania is able to trick Fc$\gamma$RIII working in tandem with IL-10 in its favor, among many other pathways of establishing infection. Because Fc$\gamma$RIII has highest affinity for IgG1 and IgG1 also acts almost exclusively through Fc$\gamma$RIII (Ravetch et al. ^[26]), Chu et al. ^[27] also investigated IgG1 KO mice for Leishmania mexicana and came to the conclusion that IgG1 is indeed pathogenic by inducing immunosuppressive IL-10 through Fc$\gamma$RIII.

Mathematical models studying the response of the immune system to different types of pathogens have been considered in the literature ^{[28,29,30,31,32]}. In ^[28], two ordinary differential equations are used to construct a model of the relationship between the immune system and a given target population. The model generates several types of immune responses when different nonlinear interactions between the immune system and its targets are considered. Although strongly idealized, the model is able to predict a threshold level for eliminating the target population by the immune system. Because time lags in population dynamics are common and in some cases natural, Buric et al. ^[29] updated the model in ^[28] by introducing time delays only in the immune response equation. The time delays enabled the model to qualitatively capture the irregular oscillations that are associated with the state of the immune system. Most recently, Mendonca et al. ^[32] further updated the model of Mayer et al. ^[28] by also introducing time delays in both the target population equation and the immune response equation. By so doing, they observed that time delays in the immune response do lead to stability switches of the fixed points. Also, delays in the target replication and cooperative immune response were shown to induce bifurcations and chaos. A few models for the immune system response to Leishmania parasite have also been considered ^{[33,34,35,36]}. In Siewe et al. ^[35], a model is formulated and used to simulate different treatment regimens for visceral leishmaniasis. The model tracks the density of macrophages, the density of dendritic cells, the density of certain T cells, the concentration of specific cytokines and the density of Leishmania parasites. Results from their model show a linear correlation between the ratio Leishmania/macrophages during the early stages of infection and that drugs which promote the proliferation of T cells effectively hinder disease progression. In De Almeida et al. ^[36], a model for cutaneous leishmaniasis is formulated, investigating the Th1/Th2 paradigm by tracking the densities of T helper cells, Th1 (that generate responses against intracellular parasites) and Th2 (that respond against extracellular parasites), and the density of Leishmania parasites. The model was used to predict disease resolution or progression when either Th1 or Th2 is absent, and also when both T helper cells co-exist. However, the existing models are limited to some aspects of the immune system response to Leishmania parasites. In this report, we follow ^[28] and model the collective response of the entire immune system to the Leishmania parasite. The model derived is used to study the effect that different genes contribute to the eradication or establishment of leishmaniasis within the host by focusing on the collective response of the immune system.

2.   Mathematical model

Following Mayer et al. ^[28], we consider the relationship between the parasite and the immune system as a feedback loop. We assume that the temporal change of the parasite population is described by the model equation

i.e., we assume a logistic growth model with the parasites normalized to the carrying capacity (so the maximum possible value of $P$ is 1) and reproduction constant $r$. The interaction between the parasite and the immune system is assumed to be concentration linked, with interaction constant $k$. The temporal change in the immune system response is modeled by

where $w$, $v$, $n$ are positive integers, $\rho$, $m$, $c$, and $d$ are positive constants. The first term describes the parasite competition response to the immune system; the constants in the expression vary the shape of the sigmoid response. The second term describes the autocatalytic or cooperative reinforcement of the immune activation processes; again, the constants in the expression give different types of sigmoid shapes. The last term represents the finite lifetime of the immune cells, with a positive death rate constant $d$.

We will assume that $w = v = n = 2$ and $m = c = 1$ for the current model.

Remark 1. We note here that in Mayer et al. ^[28], the equation representing the target population (parasite in our case) assumed an exponential growth term while in this work we have assumed a logistic growth term. The utilization of a logistic growth together with the choice of the parameter values $w = v = n = 2$ indeed depict the type of steady states (presented in subsequent sections) that govern leishmaniasis infection.

3.   Stability analysis

The basic elements of the qualitative behavior of the model given by the system of equations (2.1)–(2.2) are the equilibrium solutions. In this section, we establish the existence of equilibrium solutions and analyze them for linear stability. We start by delineating a positive solution space for the system (2.1)–(2.2). The following result shows that a positive solution of the system that starts in the positively defined closed region will continue to propagate in the closed region for all time $t\geq 0$.

Lemma 2. The biologically relevant part of the phase space, the rectangular region $0\leq P\leq 1$, $0\leq I\leq (\rho+s)/d$ is positively invariant.

Proof. The system of equations (2.1)–(2.2) shows that $P = 0$ is invariant, and $dP/dt\leq 0$ when $P\geq 1$; also, $dI/dt\geq 0$ for $I = 0$, $P\geq 0$ and $dI/dt\leq 0$ when $I\geq (\rho+s)/d$. These show that the statement of the lemma holds.

In the next result, we identify the disease-free equilibria of the system.

Lemma 3. If $s < 2d$, the only disease-free equilibrium of the system (2.1)–(2.2) is $P = 0$, $I = 0$; if $s > 2d$, there are two additional parasite-free, immune response positive equilibria given by $P = 0$, $I_1 = (s+\sqrt{s^2-4d^2})/2d$ and $P = 0$, $I_2 = (s-\sqrt{s^2-4d^2})/2d$.

Proof. From $f(P, I) = 0$ we obtain that $P = 0$ or $P = 1-kI/r$. The disease-free equilibria then will have to satisfy the equation $sI^2/(1+I^2)-dI = 0$, which gives $I = 0$ (for all values of $s$ and $d$); the other solutions of this equation then satisfy $dI^2-sI+d = 0$, and are identical to the solutions given in the statement of the lemma when $s > 2d$.

The linearized stability results for these disease-free equilibria are summarized in the following result.

Theorem 4. The disease-free, immune-response free equilibrium solution is unstable. When $s > 2d$, one of the disease-free, non-trivial immune response equilibria, $I_1$, is locally asymptotically stable if $r/k < (s+\sqrt{s^2-4d^2})/2d$, while the other one, $I_2$, is unstable.

Proof. The Jacobian of the system is given by

At the parasite-free, $I = 0$ equilibrium, $J = \text{diag}(r, -d)$, thus the system has a saddle at $(0, 0)$ and this equilibrium is unstable as claimed. The stable manifold is given by $P = 0$, extending up until the above $(0, I_2) = (0, (s- \sqrt{s^2-4d^2})/2d)$ equilibrium point (if exists, i.e., when $s > 2d$), as equations (2.1)–(2.2) show.

For $s < 2d$, there are no other equilibrium on $P = 0$; when $s\geq 2d$, there exists the pair $(0, I_1)$, $(0, I_2)$ of equilibrium points.

For $i = 1, 2$, the Jacobian at the equilibria $(0, I_i)$ is

After simplification and using that $sI_i = d(1+I_i^2)$, we obtain that $J = \text{diag}(r-kI_i, d(1-I_i^2)/(1+I_i^2))$. Now the second eigenvalue of this diagonal matrix is clearly negative for $I_1$ ($I_1 > 1$, because $s > 2d$), and it is positive for $I_2$ because a simple computation shows that in this case $I_2^2 < 1$. This means that the stability of $I_1$ will depend on whether $r-kI_1$ is positive or negative; it will be locally asymptotically stable if $r/k < (s+\sqrt{s^2-4d^2})/2d$ as claimed.

Remark 5. Lemma 3 and Theorem 4 show that there are three disease-free equilibria, $(0, 0)$, $(0, I_1)$ and $(0, I_2)$ of the system (2.1)–(2.2). The first one which is the origin is always unstable. In the context under consideration, it is understandable since there is no time when the entire immune system is absent or non-responsive. Indeed, in the absence of parasites or pathogens, aspects of the immune system are always active to a certain degree. One of the remaining two disease-free equilibria is stable while the other is unstable. The stable equilibrium $(0, I_1)$ can be attributed to the baseline immune system components that are always available to kickup (the external barrier of the skin being an example). When the immune system is activated by a pathogen, the generated response is only on need basis. That is, it is dependent on the quantity, severity and longevity of the infection. As a result, we can associate the unstable equilibrium $(0, I_2)$ with the state achieved immediately after the infection is resolved. This state then changes with time and is therefore unstable. The latter two disease-free equilibria can also represent aspects of the innate and adaptive immune system response.

The next result shows that the model contains an endemic equilibrium under appropriate conditions.

Lemma 6. If $s < d(1+r^2/k^2)/(r/k)$, then there exists at least one endemic equilibrium of the system (1)-(2).

Proof. Assume that $f(P, I) = 0$; we obtain that any non-parasite free equilibrium has to satisfy $P = 1-kI/r$. If $P = 1-kI/r$, then the equation $g(P, I) = 0$ gives, after simplification, that

This is a quintic equation with a negative leading coefficient, which satisfies $h(0) = \rho > 0$, thus it is guaranteed that at least one nonnegative root exists. Because the corresponding $P$ value is given by $P = 1-kI/r$, in order to get a non-parasite-free equilibrium, we need that at least one root of $h(I)$ has to satisfy $I < r/k$. A sufficient assumption for this to happen is $h(r/k) < 0$, which is the same as $-dr^3/k^3+sr^2/k^2-dr/k < 0$ or, after simplification, $s < d(1+r^2/k^2)/(r/k)$ as claimed.

The next result provides sufficient conditions for the local asymptotic stability of any endemic equilibrium solution whose existence is guaranteed by the preceding lemma.

Lemma 7. If $s < d(1+r^2/k^2)/(r/k)$ and $s < 8d/(3\sqrt{3})$, then any endemic equilibrium $(P^*, I^*)$ corresponding to a positive root $I^* < r/k$ of $h(I)$ is locally asymptotically stable.

Proof. By Lemma 4, at least one positive $I^* < r/k$ exists. The Jacobian in this case becomes

The function $2x/(1+x^2)^2$ has a maximum of $3\sqrt{3}/8$ at $x = 1/\sqrt{3}$, thus the element in the second row and second column of the Jacobian is guaranteed to be negative in case $s < 8d/(3\sqrt{3})$. This then implies that the trace of $J$ is negative, while the determinant is positive, so the equilibrium is locally asymptotically stable as claimed.

Lemma 8. Lemma 6 and Lemma 7 establish the existence of a biologically relevant positive endemic equilibrium $(P^*, I^*)$ that results in the resolution or persistence of an infection. It is important to note here that in the case of leishmaniasis, infection resolution does not necessarily imply complete eradication of the parasites as there is always a residual amount left. Rather, resolution implies that the immune system has generated a sufficient response that reduces the infection to the point where it is not harmful to the host.

4.   Numerical results and discussion

The model presented in this work does not have a closed form solution using elementary functions. In order to consider the time evolution of the model, we numerically investigate it for a chosen set of parameter values. We begin by carrying out a phase plane analysis and investigate the system for bifurcations, and afterward compare the model with experimental data for leishmaniasis.

4.1.   Bifurcation analysis

In the analysis below, we will assume that our bifurcation parameter is the autocatalytic response size parameter $s$; we will fix the other parameters and investigate how the phase plane is changing if we vary the value of $s$.

Supposition 9. Assume that $r = 1/2$, $k = 1$, $\rho = 1$ and $d = 2$. By Lemma 3, the parasite-free non-trivial equilibria exist when $s > 2d = 4$. The sufficient condition of Lemma 6 for $I^* < r/k = 1/2$ is given by $s < (5/2)d = 5$. Local asymptotic stability of the endemic equilibrium for $(P^*, I^*)$ is guaranteed by Lemma 7 for $s < 16/(3\sqrt{3})$.

For these fixed values, the roots $I^*$ satisfying $I^* < 1/2$, the corresponding $P^*$ values, and the eigenvalues of the Jacobian at $(P^*, I^*)$ are as follows:

for $s = 3.0$: $I^*\approx 0.19098$, $P^*\approx 0.61803$, $\lambda_{1, 2}\approx -0.6212\pm 0.6959i$;

for $s = 3.6$: $I^*\approx 0.20161$, $P^*\approx 0.59677$, $\lambda_{1, 2}\approx -0.4789\pm 0.7031i$; and

for $s = 4.2$: $I^*\approx 0.21548$, $P^*\approx 0.56904$, $\lambda_{1, 2}\approx -0.3157\pm 0.6779i$.

These are all locally asymptotically stable with illustrations given in Figure 1, Figure 2 and Figure 3 for values of $s = 3.0$, $s = 3.6$ and $s = 4.2$, respectively. The only noticeable change in the system as the value of $s$ increases from $3.0$ to $3.6$ (see Figure 1 and Figure 2) is an increase in the magnitude of the endemic equilibrium. However, as the value of $s$ transitions between $3.6$ and $4.2$, a significant change occurs. For the last value of $s$, we have another locally asymptotically stable solution and a nonzero disease-free equilibrium solution $(0, I_1)\approx (0, 1.37016)$ that is also locally asymptotically (see Figure 3). The basins of the two equilibria in Figure 3 are separated from each other by the unstable manifold of $(0, I_2)$.

Remark 10. The phase plane results in Figure 1 and Figure 2 show that irrespective of the strength (or size) of the initial infection and the state of the immune response, the system will always converge to the endemic equilibrium solution. This paradigm is no longer maintained as the value of $s$ increases from $3.6$ to $4.2$. As indicated in Figure 3 a bifurcation occurs leading to two basins of attraction delineated by a separatrix. The separatrix divides the phase plane into two sections, a bottom part containing the endemic equilibrium and a top section containing the nonzero disease free equilibrium. An initial infection originating within the lower basin converges to the endemic equilibrium while an infection that begins in the top section converges to the nonzero disease free equilibrium. All infections originating in the basin of the nonzero disease free equilibrium are eradicated. We note here that an infection that originates in the basin of the endemic equilibrium can only be eradicated if it is coupled with treatment that pushes it across the separatrix into the basin of the nonzero disease-free equilibrium.

Another bifurcation that we did observe is that raising the value of $s$ from $4.2$ to $5.0$, we obtain more possible endemic equilibria as the second positive root $I^{**}$ of the quintic becomes less than $r/k$. Also, the new pair of nonzero parasite-free equilibria is created by a saddle-node bifurcation, and the endemic equilibria lose stability. This shows that in case of larger $s$ values (the autocatalytic reaction is strong), one of the nonzero parasite-free equilibria becomes locally asymptotically stable, and the long term behavior indicates the decline of the parasite $P$.

The parameter values were chosen in Supposition 9 so that the ratio $r/k < 1$. In the next supposition, we choose the parameter values such that the ratio $r/k > 1$. Notice that all parameters values have been maintained except the value of $k$.

Supposition 11. Assume that $r = 1/2$, $k = 1/3$, $\rho = 1$ and $d = 2$. By Lemma 3, the parasite-free non-trivial equilibria exist when $s > 2d = 4$. The sufficient condition of Lemma 6 for $I^* < r/k = 3/2$ is given by $s < (13/6)d = 13/3$. Local asymptotic stability of the endemic equilibrium for $(P^*, I^*)$ is guaranteed by Lemma 7 for $s < 16/(3\sqrt{3})$.

In a similar way, for these fixed values, the roots $I^*$ satisfying $I^* < 1/2$, the corresponding $P^*$ values, and the eigenvalues of the Jacobian at $(P^*, I^*)$ are as follows:

for $s = 3.0$: $I^*\approx 0.34675$, $P^*\approx 0.76884$, $\lambda_{1, 2}\approx -0.3633\pm 0.5995i$;

for $s = 3.6$: $I^*\approx 0.73034$, $P^*\approx 0.51311$, $\lambda_{1, 2}\approx -0.0101\pm 0.2323i$; and

for $s = 4.2$: $I^*\approx 1.38038$, $P^*\approx 0.07975$, $\lambda_{1, 2}\approx -0.6243, -0.0420004$.

Notice also that, these are all locally asymptotically stable. The numerical simulations are given in Figure 4, Figure 5 and Figure 6 for values of $s = 3.0$, $s = 3.6$ and $s = 4.2$, respectively. Observe that the changes associate with this case on raising the value of $s$ show up early compared to the case when $r/k < 1$. The first bifurcation turns the locally asymptotically stable equilibrium point from a focus into a node as $s$ transitions from $s = 3.0$ to $s = 3.6$ (see Figure 4 and Figure 5). As the value of $s$ increases from $s = 3.6$ to $s = 4.2$, another bifurcation takes place yielding an endemic equilibrium solution that is locally asymptotically stable and a nonzero disease-free equilibrium solution $(0, I_1)\approx (0, 1.37016)$ that is unstable within the same basin of attraction (see Figure 6).

Remark 12. The phase plane analysis in Figure 4, Figure 5 and Figure 6 clearly show that all initial infections will always converge to the endemic equilibrium state that is locally asymptotically stable. The analysis also show that the autocatalytic reaction is weak for small $s$ values and strong for large $s$ values. A more careful observation of the time domain solutions in Figure 4, Figure 5 and Figure 6 show that the locally asymptotically stable endemic equilibrium correspond to parasite persistence for small $s$ values and parasite resolution for large $s$ values. This case renders the model applicable to Leishmania infection because a residue of the parasite is always retained when an infection is classified as resolved.

4.2.   Application to leishmaniasis

We will now apply the model to cutaneous leishmaniasis infection. In particular, we focus on Leishmania mexicana (L. mexicana) and compare the model prediction with experimental data. Before proceeding, we comment here that the model presented is not limited to Leishmania mexicana. With appropriate choices of the parameter values, it can be used to study other species of Leishmania and also other infectious diseases. Also, as mentioned earlier, the motivation to apply the model to Leishmania infection is due to the nature of the endemic equilibria involved. One of the endemic equilibria leads to disease persistence while the other leads to infection resolution. An infection is considered resolved if a growing lesion reverses course and in the process gets eliminated. Alternatively, an infection is considered resolved when the parasite load falls below the initial amount. In practice, there is always a residue of parasites left after an infection is deemed as resolved.

The data used below were obtained from different murine studies ^[24,25,27] where C57BL/6(B6) mice lacking specific genes (also known as knocked out (KO)) and normal B6 mice were infected with L. mexicana. In each experiment, the mice were inoculated with $5$ million parasites in the hind footpad and lesions were monitored for several weeks. Experimental data was reported using the size of the lesion with the first measurement taken $2$ weeks after the inoculation.

We present two sets of results below related to the parasite load. The first set of results given in Figures 7-10, represent cases in which the lesion size is used as a proxy for parasite load. The second set of results is given in Figures 11-12 where the actual parasite load is considered. We do not present corresponding results for the immune response because no experimental data is available that we can use for comparison. In all simulations reported below, parameter values were chosen to explore the behavior of the model, guided by values reported in the published literature ^[35]. In Siewe et al. ^[35], the parasite growth rate in pro-inflammatory macrophages is given as $3.09$ per day and the parasite growth rate in anti-inflammatory macrophages is given as $3.82$ per day. Taking the average of these values, we set the parasite growth rate at $r = 1/2$ per week. Also in ^[35], the parasite death rate in pro-inflammatory macrophages attributed to Nitric Oxide and IFN-$\gamma$ (which are aspects of the immune response) is given as $1.85$ per day, and the parasite death rate in anti-inflammatory macrophages is given as $2.22$ per day. Guided by the average of these values, we set the parasite clearance rate by the entire immune response at $k = 1/3$ per week. Further, in Siewe et al. ^[35], the immune system strength factor is given by the range $(0-150)$. Because we are considering the entire response of the immune system which is determine by the parameters $\rho$, $s$ and $d$, we choose them within the range $(0-150)$. The choices we made present numerical solutions that are in close agreement with experimental data.

In the first set of results, unless otherwise stated, we fixed the parameter values at $r = 1/2$, $\rho = 1$, $k = 1/3$, and $d = 2$. We also fixed the autocatalytic response parameter at $s = 3.0$ in simulating normal B6 mice in all cases, while prescribing a different $s$ value for each KO case. In all simulations reported below, we assumed that the initial condition of the immune response is $I(0) = 0$. Initial conditions for the parasite loads for both wild-type and KO mice are also assigned for each case reported. We ran several simulations with the initial conditions for the parasites loads perturbed and the results were similar to the ones presented in this work.

Using the experimentally observed lesion sizes as proxy for parasite loads, we simulate the murine studies in ^[24,25,27]. Before we continue, it is necessary to pause and remark here that the nature of the relationship between lesion size and parasite load is indeed not well established. This is because apart from parasites, lesions are composed of immune cells like T cells, infected macrophages, and some other material/debris. Indeed, some infections do lead to parasite clearance but not to lesion resolution. However, the cases reported in this paper do have chronic disease in wild-type B6 mice with plateauing lesions and parasite loads vs knock-out mice that heal lesions and clear parasites. Therefore, it is important to emphasize that the simulations are meant to predict the overall general behavior of the parasite load, not to recreate the exact lesion size. The main point is that the change in the value of the autocatalytic response parameter changes the behavior of the solution of the system the exact same way the lesion size is changing in the different cases. Thus, the model includes the possibility of different behavior in the different cases.

As our starting point, we compare the model simulations to the experimental results of Buxbaum et al. ^[24]. It was observed in ^[24] that B6 mice deficient in IL-$10$, Fc$\gamma$R and $\beta 2$m were able to resolve L. mexicana infection while normal B6 mice containing these genes did not. Below, we only present model simulations for IL-$10$ and $\beta 2$m. To apply the model in the case of $\beta 2$m KO and normal B6 mice, we set the autocatalytic response parameter at $s = 6.0$ for the $\beta 2$m KO mice. The initial conditions for the wild-type mice and the KO mice were set at $P(0) = 0.1$ and $P(0) = 0.02$, respectively. A comparison of our numerical simulations in time domain and experimental data for $\beta 2$m KO and normal B6 mice is given in Figure 7. We see a strong agreement between the model results and experimental data. The parasite load in the wild-type mice approaches the stable endemic equilibrium with $P^* = 0.769$ indicating infection persistence while the parasite load in KO mice approaches the stable non-endemic equilibrium with $P^* = 0$, signifying infection resolution.

Next, we apply the model to experimental data for IL-$10$ KO and normal B6 mice. We maintain all parameter values as earlier defined and set the autocatalytic response parameter at $s = 4.5$ for IL-$10$ KO mice. The initial conditions for the wild-type mice and the KO mice were set at $P(0) = 0.2$ and $P(0) = 0.02$, respectively. The numerical results in time domain are given in Figure 8. We also see a very good correlation between model simulations and experimental data in this case. Similar to the preceding case, here, the parasite load in the wild-type mice also approaches the stable endemic equilibrium with $P^* = 0.769$ indicating infection persistence while the parasite load in KO mice approaches the stable non-endemic equilibrium with $P^* = 0$, leading to infection resolution.

Although not considered in this work, experimental observations in ^[24] showed that B6 mice deficient in Fc$\gamma$R protein were also able to resolve infection with L. mexicana. This led to further studies by Thomas et al. ^[25] where Fc$\gamma$RIII, a subclass of Fc$\gamma$R was investigated. It was also established in ^[25] that mice deficient in Fc$\gamma$RIII did resolve infection with L. mexicana while normal B6 mice did not. To apply the model to this case, retaining all other parameter values, we set the autocatalytic response parameter at $s = 4.1$ for Fc$\gamma$RIII KO mice. The initial conditions for the wild-type mice and the KO mice were set at $P(0) = 0.2$ and $P(0) = 0.1$, respectively. The results of our computation in time domain and experimental data is presented in Figure 9. We see here as in previous cases that the time evolution of the infection closely matches with experimental observations. The parasite load in the wild-type mice still approaches the stable endemic equilibrium with $P^* = 0.769$ because the parameter values remain the same as in previous cases while the parasite load in KO mice approaches the stable endemic equilibrium with $P^* = 0.145$. Note that because $P^* < P(0)$ we conclude that the infection is resolved as the parasite load approaches the stable endemic equilibruim $P^* = 0.145$.

Finally, we apply the model to experimental data for IgG1 reported in Chu et al. ^[27]. Considering the fact that IgG1 binds mostly to Fc$\gamma$RIII and because the presence of Fc$\gamma$RIII had been implicated in infection progression, Chu et al. ^[27] investigate the role of IgG1 protein. They came to the conclusion that IgG1 KO mice were able to resolve infection with L. mexicana while the normal B6 mice did not. In order to apply the model to this case, we again maintain all parameter values as earlier defined with the exception of the autocatalytic response parameter set at $s = 4.2$ for IgG1 KO mice and the value of $\rho = 2$. The initial conditions for the wild-type mice and the KO mice were set at$P(0) = 0.2$ and $P(0) = 0.15$, respectively. The results in time domain for this case are given in Figure 10. Clearly, the model simulations do mimic experimental observations depicting infection resolution in IgG1 KO mice and infection persistence in normal B6 mice. By observing the differences in the experimental results for wild-type mice from the three studies (all carried out in a similar way), we altered the value of $\rho$ from $1$ to $2$ in this case. With this change, the parasite load in the wild-type mice now approaches the stable endemic equilibrium with $P^* = 0.513$ while the parasite load in KO mice approaches the stable endemic equilibrium with $P^* = 0.075$. We also say here that, beacuse $P^* < P(0)$ the infection is resolved as the parasite load approaches the stable endemic equilibruim $P^* = 0.075$.

In the second set of results, we use the model to simulate the actual parasite load for the murine studies in Thomas et al. ^[25]. We note here that although the studies ran for about $24$ weeks, parasite loads were only obtained for weeks $4$, $8$, $12$ and $23$. The process involved grinding a severed infected footpad of a sacrificed mice and analyzing collected samples. The parasite load was then calculated based on the sample amount. As with the first set of results, we also choose parameter values so as to get the best possible outcome. Although each animal was inoculated with $5$ million parasites at the onset, we assume that the majority of these were destroyed by the immune system before they were able to successfully obtain shelter within macrophages or other phagocytes and replicate. We assume that about $40,000$ and $50,000$ parasites were able to get into macrophages in the Fc$\gamma$RIII KO mice and wild-type B6 mice, respectively. We also set the carrying capacities at $10$ million and $25$ million for the Fc$\gamma$RIII KO mice and wild-type B6 mice, respectively. The parameter values associated with the parasite equation were maintained at $r = 1/2$, $k = 1/3$, $d = 2$, while those for the immune response equation, $\rho$ and $s$, were altered in each case. The values of $\rho$ and $s$ needed by the model to adequately simulate the actual parasite load were significantly higher for the knock-out mice than those for the wild-type mice.

The results in Figure 11 illustrate the time evolution of the actual parasite load for Fc$\gamma$RIII KO mice vs measured data. With the exception of the fourth week, the simulation was within the standard errors associated with the experimental data. Similarly, the results in Figure 12 show the time progression of the parasite load in wild-type B6 mice vs measured data. The simulation was also mostly within bounds of observed data.

5.   Conclusion

In this paper we developed a model of the immune response to threats from parasitic organisms. The model constitutes a coupled system of differential equations in which one of the equations represents the parasite density and the other equation is a simplified representation of the immune response. The parasitic equation follows logistic growth while the immune response equation is determined by parasitic activation and cooperative reinforcement.

The model possesses a zero equilibrium solution which is also a disease-free equilibrium that is unstable. In addition to the zero solution, the model also contains two other disease-free equilibria, one stable and the other unstable. Further, the model possesses a positive endemic equilibrium that is stable. A numerical bifurcation analysis for chosen parameter values show that the endemic equilibrium transitions between disease resolution and transition as the autocatalytic response of the immune system changes. The behavior of the endemic equilibrium makes the model particularly suitable for application to leishmaniasis.

Leishmaniasis, of which are many species remains endemic in many regions of the world. While it can be successfully treated, there is yet to be an effective vaccine for controlling the disease. Clinical observations aimed at understanding the parasite's pathway of establishing infection have been undertaken in recent decades. The majority of these observations are murine based and have demonstrated the elusiveness of the parasite through manipulation of the immune response to its advantage. As an example, studies of L. mexicana have implicated aspects of the immune response such as the cytokine IL-$10$, antileishmanial antibody IgG1, Fc$\gamma$ receptor, and the subclass receptor Fc$\gamma$RIII in disease persistence.

The model was used to study different cases of L. mexicana infections in which wild type (normal) B6 mice where compared to IL-$10$ KO or $\beta 2$m KO or IgG1 KO or Fc$\gamma$RIII KO mice. In experimental studies ^[24,25,27], the infection was observed to persist in the normal B6 mice while it was resolved in mice with KO genes. The numerical simulations of the model formulated in this work are in strong agreement with experimental observations. Furthermore, the bifurcations associated with the endemic equilibrium as the autocatalytic response changes correspond to disease resolution in the KO mice and disease persistence in the normal B6 mice for each infection case.

In conclusion, we remark that the accurate predictions of the model are mostly obtained by varying the value of the parameter that determines the strength of the autocatalytic response. This makes the model more appealing for investigating L. mexicana. We postulate that apart from the KO genes reported in this work, other genes with implication in disease persistence can be predicted using formulated model.

Acknowledgments

This work was supported by a summer Faculty Education and Development (FEAD) grant from the College of Science, Rochester Institute of Technology. The experimental studies reported in this work were supported by a Veterans Affairs Merit Review Grant, NIH K08-AI01805, NIH R01AI81717.

Conflict of interest

No conflict of interest.