Dynamic study of the pathogen-immune system interaction with natural delaying effects and protein therapy

1.   Introduction

The immune system is a complex system consisting of many interrelated processes that form a coordinating system that aims to protect the integrity and identity of the host. This system has a role in preventing the invasion of harmful substances such as pathogens in the environment that can self-destruct ^[1,2]. The major role of the immune system is to effectively protect the host against pathogens by identifying and destroying pathogen-infected cells. Identification and destruction processes of pathogens involve two stages, namely innate immunity and adaptive immunity which influence each other because the innate immune response plays an important role in infection control during periods of delayed activation of the adaptive immune response. In addition, cells of the innate immune system initiate and determine the subsequent direction of the adaptive immune response and participate in the elimination of pathogens that are targets of the adaptive immune response ^[1,3].

All cells in the immune system are produced and reproduced in the bone marrow as immature cells. The immature cells are then educated to become mature cells of the immune system which can be categorized as lymphocytes, neutrophils, and macrophages ^[4]. T cells are a subset of lymphocytes which are the main cellular component of the adaptive immune response because they can directly attack infected cells. T cells consist of three types of cells namely cytotoxic, helper, and regulatory that react strongly to foreign antigens to be effective for immunity. As part of the immune system, T cells focus on identifying certain foreign particles and circulating them until they encounter specific antigens ^[3]. Since T cells are self-tolerant, they have the ability to recognize harmful antigens and self-produced antigens by identifying molecules on the cell surface. In addition to T cells, lymphocytes also have a subset called B cells, specialized cells with the main function of producing antibodies, i.e. the main proteins that are produced when cells of the immune system react with foreign protein antigens. B cells also mature in memory cells so that the system can respond quickly when a similar attack is encountered. Together with proteins, cells collaborate to recognize and react to foreign substances in providing defense against invading pathogens. Therefore the immune system becomes a wonderful collaboration between cells and proteins which also results in a very complex system with inherent nonlinear properties ^[2–4]. The development of science and technology, especially in recombinant technology, has made great progress in the development of proteins such as laboratory-engineered therapeutic proteins which showed highly effective in replacing abnormal or deficient proteins in the host's body as well as increasing the body's supply of beneficial proteins to help reduce the effects of disease or chemotherapy ^[5–9].

Over the past decade, theoretical studies have been continuously developed to capture and understand the complex multiscale dynamics of the immune system and its immunogenicity that may not be captured directly in experiments. Eftimie et al. in ^[10] presented a study review regarding the development of mathematical model tools as a theoretical approach in studying immunology systems qualitatively and quantitatively. Moreover, some mathematical approaches were also developed and analyzed to study the complexities of the immune system involving many interacting processes ^[11–17]. Some of them are focused on the innate immunity stage and some on the adaptive immune stage or both which include humoral aspects such as antibodies, cell-mediated aspects, and immunity effects and their responses to some deadly diseases such as tumors and other infectious diseases ^[11–15,17]. Some researchers are also concerned to study and observe the effects of delays occurring in the immune system ^[18–29]. Mathematical models were developed as a series of structured population models to describe interactions between infectious agents and immune cells with a series of delay differential equations. However, for simplicity, some assumptions and limitations about the nature of various processes were made such as regarding the time delay between pathogen injection and immune activation. For example, Fenton et al. ^[18] considered a structured population model showing the presence of specific immune cells that were not directly activated and reached their maximal efficacy. Gradual enhancement of immune efficacy was considered to be a major aspect in changing the interaction dynamics of the pathogen immune system and the ability of the immune response to eliminate pathogens ^[18]. However, the distributed immune system model (immune pathway) produces a structured population model chain with high-dimensional variables that are not easy to analyze analytically. In this paper, we intend to extend the study of Fenton et al. by proposing another approach and developing a mathematical model which considers the distributed immune system (the immune pathway) as a lumped process. This extension considers only major parts of the immune pathway called primary and specific immune cells, and replaces long distributed responses with specific natural delays i.e. a delay between early pathogen infection and immune response, and a delay of immune cell maturation. Other delay is also considered from the first invasion of pathogens until the initial formation of memory cells. We consider a competitive interaction between the pathogen and the host immune system under the influence of adding therapeutic proteins to trigger the immune system. The model is generated as a system of five nonlinear delay differential equations involving three discrete delays. We believe that the proposed model will complement the results of the previous theoretical studies on the complex dynamics of the immune response system with inherent nonlinear nature.

2.   Mathematical model formulation

A mathematical model of the immune system is formulated to study the interactions between pathogens and immune cells with the addition of antigens contained in a therapeutic protein that triggers antibody production. Our model follows the line of Fenton et al. ^[18] with some extensions in their model. As we have stated before that the distributed immune system model in Fenton et al. is considered as a lumped pathway such that we consider only the main part of the immune system. Here, we focus on the five compartments namely P which states as the number of pathogens in the host body, I_s which states the number of specific immune cells, I_p which states the number of primary immune cells, M which determines the number of memory cells, and D which states as the number of therapeutic proteins injected into the host body to enhance the production of primary cells (see Figure 1 for illustration of the interaction between compartment). Moreover, our mathematical model was formulated based on the following general assumptions:

1. Pathogens have the ability to reproduce themselves in the cell's body by using facilities and nutrients from the host body.

2. Host body is unable to respond directly to the pathogen injection such that a delay occurs on the first injection until the primary immune cells are activated.

3. On the other hand, the primary immune cells cannot directly act to attack the pathogens. It requires a certain time for maturation and growth to become specific immune cells that are ready to operate the immunity function. On the other hand, the specific immune cells are produced after pathogens have passed through a non-specific immune system.

4. Pathogens can quickly evolve after an interaction with the immune cells so that they can pass through immune cells. This is a part of pathogen self-defense through producing low protein so that it can not be detected by the immune cells.

5. Therapeutic proteins are not produced by the host body. It is injected into the body as antigens which are needed by the body and stimulate the production of primary immune cells. Doses of therapeutic proteins are assumed constant over the period of taking.

6. Immune cells have the ability to remember antigens entered the body by producing memory cells such that when similar antigens reenter the body, the primary immune response will be faster. Therefore, there is a lag by the inclusion of an explicit delay term before initiation of the immune memory.

7. Interaction between pathogens and cells is a competitive interaction where the pathogens are stronger than the immune cells. The number of immune cells will be reduced but the brain will instruct the body to produce more optimal immune cells to destroy the pathogens.

More specifically, assumptions for deriving the mathematical model are explained as follows:

1. Pathogens can replicate to multiply by using the facilities needed by the host such as nutrition so that the pathogen will grow exponentially with the replication rate $a$. Pathogen cells will die due to age or other factors with natural death rate $\omega$. Because pathogens need nutrition in their development there will be competition between fellow pathogens so that there is a logistical growth model by including environmental carrying capacity as $b$. The existence of immune cells will cause a number of pathogens because the immune cells will destroy objects that are considered foreign to the body with $\sigma$state the reduction rate of pathogens. The pathogen evolution process causes pathogens to avoid detection of immune cells so the number of pathogens will increase at the rate $\mu$.

2. The presence of pathogens in the body will stimulate the production of primary immune cells. However, there is a delay between the injection of pathogens and the initiation of primary immune cell production symbolized by ${\tau }_{1}$, where $c$ is the rate of cellular and biochemical reactions in the body and ${\lambda }_{0}$ is the level of immunity depending on the infection. The addition of therapeutic protein drugs will stimulate the formation of primary immune cells because they are considered as antigens by the body, so the number of primary immune cells produced will increase by $k$ which states the binding rate of therapeutic protein drugs by immune cells. When the same pathogen injection occurs in the future, the components of memory cells will be activated to stimulate the formation of immune cells again at the rate of ${\lambda }_{M}$and will destroy pathogens as before.

3. Specific immune cells need time to mature their cells so that their functions become more effective in destroying pathogens and symbolized as ${\tau }_{2}$. The body produces specific immune cells to fight pathogens and will be reduced by competition but the brain will send signals so that regeneration will occur and new stronger immune cells will be produced where $\alpha$ is the rate of regeneration of immune cells and $\beta$ is the rate of reduction of immune cells because it loses against pathogens. The number of immune cells will increase if $\alpha$ > $\beta$, namely the rate of regeneration is greater than the rate of competition with pathogens and will decrease if $\beta$ > $\alpha$.

4. The production of memory cells usually takes place 14 days after the first pathogen invasion, therefore a delay of ${\tau }_{3}$ occurs. Memory cells will be produced together with primary immune production at the rate of $\gamma$ and decay at the rate of $\delta$. The doses of therapeutic protein drugs in the body depend on two things, increasing because the rate of input of therapeutic protein drug doses which are symbolized as ${\lambda }_{D}$, and decreasing because it is bound and fused with primary immune cells which is symbolized as $k$.

Based on these general and specific assumptions, we formulate a nonlinear delay differential system as follows:

with initial conditions $M\left(0\right) = {M}_{0} > 0,D\left(0\right) = {D}_{0} > 0,$and the historical functions: $P\left(t\right) = {\beta }_{1},\forall t\in \left[-{\tau }_{1},0\right],{I}_{S}\left(t\right) = {\beta }_{2},\forall t\in \left[-{\tau }_{2},0\right],{I}_{p}\left(t\right) = {\beta }_{3},\forall t\in \left[-{\tau }_{3},0\right]$, for ${\beta }_{1},{\beta }_{2},{\beta }_{3}\in {R}^{+}$.

3.   Steady state solution and asymptotic stability analysis

Now, let us turn our consideration to the existence of steady states for the system (1). We admit steady state solutions defined in the restriction region, ${\mathbb{R}}_{+}^{5} = \left\{\left(P,{I}_{p},{I}_{s},M,D\right)\in {\mathbb{R}}_{}^{5}|P\ge 0,{I}_{p}\ge 0,{I}_{s}\ge 0,M\ge 0,D\ge 0\right\}$. Since the vector field does not point to the exterior of ${\mathbb{R}}_{+}^{5}$ then Eq (1) is defined in the region ${\mathbb{R}}_{+}^{5}$. Steady state of (1) is the stationary solution ${E}^{*} = \left({P}^{*},{{I}_{p}}^{*},{{I}_{s}}^{*},{M}^{*},{D}^{*}\right)$ that fulfills,

With several mathematical processes we have the following results.

Proposition 1. System (1) has a nontrivial positive steady state solution ${E}^{*} = \left({P}^{*},{{I}_{p}}^{*},{{I}_{s}}^{*},{M}^{*},{D}^{*}\right)$ with ${P}^{*} = \frac{c}{\beta -\alpha }$, ${{I}_{p}}^{*} = \frac{\delta \left(c{\lambda }_{0}+{\lambda }_{D}(\beta -\alpha )\right)}{\delta \left(\beta -\alpha \right)-{\lambda }_{M}c\gamma }$, ${{I}_{s}}^{*} = \frac{\left(a-\omega \right)\left(\beta -\alpha \right)-bc}{(\beta -\alpha )\left(\sigma -\mu \right)}$, ${M}^{*} = \frac{\gamma \left(c{\lambda }_{0}+{\lambda }_{D}(\beta -\alpha )\right)}{\delta \left(\beta -\alpha \right)-{\lambda }_{M}c\gamma },$ and ${D}^{*} = \frac{{\lambda }_{D}\left(\delta \left(\beta -\alpha \right)-{\lambda }_{M}c\gamma \right)}{k\delta \left(c{\lambda }_{0}+{\lambda }_{D}(\beta -\alpha )\right)}$ if it fulfills conditions $\beta > \alpha$, $\delta \left(\beta -\alpha \right) > {\lambda }_{M}c\gamma$, $\left(a-\omega \right)\left(\beta -\alpha \right) > bc$ and $\sigma > \mu$, or $\left(a-\omega \right)\left(\beta -\alpha \right) < bc$ and $\sigma < \mu$.

Next, we focus on the analysis of the asymptotic stability of ${E}^{*}$. To this aim, we linearize (1) about ${E}^{*}$ and determine the associated characteristic equation. System (1) can be written as,

By linearizing system (6) near the equilibrium point ${E}^{*}$, we get a linear form of (6) as follow:

where

and ${J}_{3} = \left[\begin{array}{ccccc}0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0\\ 0& \gamma & 0& 0& 0\\ 0& 0& 0& 0& 0\end{array}\right]$,

with

The characteristic equation associated with (7) is given by

with

Equation (8) is a five-degree exponential polynomial in $\lambda$. The local asymptotic stability analysis of the steady state ${E}^{*}$ can be performed by identifying the sign of the real parts of the roots of (8). The steady state ${E}^{*}$ is locally asymptotically stable if and only if all roots of (8) have negative real parts, and its stability can only be lost if purely imaginary roots appear. Note that Eq (8) depends only on two discrete delays, i.e., ${\tau }_{2}$ and ${\tau }_{3}$ with free ${\tau }_{1}$ meaning that stability of ${E}^{*}$ is only influenced by ${\tau }_{2}$ and ${\tau }_{3}$. Nevertheless, it is not easy to perform the stability analysis directly via investigating the coefficients of the polynomial (8) due to its implicit relation with $\lambda$. Therefore, analysis of the roots of Eq (8) will be presented in some cases as follow.

Case 1: ${\bf{ \pmb{\mathsf{ τ}} }}_{2} = {\bf{ \pmb{\mathsf{ τ}} }}_{3} = 0$ with ${\mathit{\tau }}_{1} > 0$

Assume that ${\tau }_{2} = {\tau }_{3} = 0$ (there is no delays on the maturation of primary cells and activation of memory cells or the delays can be ignored), then the characteristic Eq (8) can be written as a five-degree polynomial equation,

with

Proposition 2. Based on the Routh-Hurwitz stability criteria, Eq (9) has roots with negative real part if and only if it meets the following conditions (H1):

Case 2: ${\bf{ \pmb{\mathsf{ τ}} }}_{1},{\bf{ \pmb{\mathsf{ τ}} }}_{2} > 0,{\bf{ \pmb{\mathsf{ τ}} }}_{3} = 0$

Assumed that the pathogen has invaded the host body meaning that the memory cells already exist and can quickly stimulate the formation of primary immunity (so that ${\tau }_{3} = 0$). It means that the body can quickly respond the invasion of pathogen such that infection occurred (it also can make ${\tau }_{1}\to 0)$. Then Eq (8) becomes,

with

Since Eq (10) will become Eq (9) when ${\tau }_{2} = 0$, then based on the Proposition 2, the equilibrium point ${E}^{*}$ is locally asymptotically stable. When ${\tau }_{2} > 0$, the roots of Eq (10) are not easy to obtain explicitly. The eigenvalues of Eq (10) will depend on ${\tau }_{2}$. Suppose that the eigenvalue of Eq (10) is in a complex form, $\lambda \left({\tau }_{2}\right) = v\left({\tau }_{2}\right)\pm i\omega \left({\tau }_{2}\right)$ with ${\tau }_{2} > 0$. To find out whether ${E}^{*}$ is stable when ${\tau }_{2} > 0$or has a limit cycle, the analysis will be processed by assuming the eigenvalue of Eq (10) is in an imaginary form, $\lambda = \pm i\omega ,\omega \in \mathbb{R}$. If $i\omega$ is a purely imaginary root of (10), then it fulfills

Simplifying the equation, we get

By taking its real part and imaginary part become zero, we have

and

From Eq (11) we have

While from Eq (12) we have

From Eqs (13) and (15) we get

with ${s}_{1} = {r}_{6}{\omega }^{3}-{z}_{6}\omega$, ${s}_{2} = {r}_{5}{\omega }^{4}-{z}_{5}{\omega }^{2}+{z}_{7}$, ${s}_{3} = {\omega }^{5}-{z}_{1}{\omega }^{3}+{z}_{3}\omega$, and ${s}_{4} = {r}_{1}{\omega }^{4}-{z}_{2}{\omega }^{2}+{z}_{4}$. From Eqs (14) and (16) we get

If $\omega$ is the solution of Eqs (17) and (18) then $\text{-}\omega$ is also the solution of both equations. Therefore, in the following, we only look for positive solutions $\omega$ of (17) and (18). Adding the squares of Eqs (17) and (18) we obtain,

Consider ${p}_{1} = {{r}_{1}}^{2}-{{r}_{5}}^{2}-2{z}_{1}$,

Then Eq (19) can be written as

If we consider $X = {\omega }^{2}$ then we have

Assume

$\left({\boldsymbol{H}_{2}}\right)$ ${p}_{5} < 0$, ${p}_{4} > 0$, ${p}_{3} < 0$, ${p}_{2} > 0$, and ${p}_{1} < 0$.

If condition ${\boldsymbol{H}}_{2}$ is fulfilled then all roots of $F\left(X\right)$ are positive, say ${X}_{l},l = \mathrm{1,2}\dots ,5$ such that ${\omega }_{l} = \sqrt{{X}_{l}}$. Therefore condition ${\boldsymbol{H}}_{2}$ guarantees the existence of imaginary eigenvalues of (10) where a Hopf bifurcation will probably arise. This presumption can be adjusted by investigating the sign of $\frac{dRe\left(\lambda \right)}{d{\tau }_{2}}$.

Now, consider Eqs (17) and (18). For $\omega = {\omega }_{l}$, we have the solutions of both equations,

if $\mathrm{sin}\left({\omega }_{l}{\tau }_{2}\right) > 0$, and

if $\mathrm{sin}\left({\omega }_{l}{\tau }_{2}\right)\le 0$. If we define two sequences $\left\{{\tau }_{2,l}^{1,j}\right\}$ and $\left\{{\tau }_{2,l}^{2,j}\right\}$ for $l = 1,\dots ,5$ and $j\in \mathbb{N}$ then we have

Lemma 1. Suppose that ${\tau }_{2,l}^{*}$ is an element of either the sequence $\left\{{\tau }_{2,l}^{1,j}\right\}$ or $\left\{{\tau }_{2,l}^{2,j}\right\}$ with ${\omega }_{l}$. Then the characteristic equation (10) has a pair of conjugate pure imaginary roots $\lambda = \pm i{\omega }_{l}$ that fulfils

Proof:

Consider the roots of (10), $\lambda \left({\tau }_{2}\right) = v\left({\tau }_{2}\right)+i\omega \left({\tau }_{2}\right)$ where $v\left({\tau }_{2,l}^{*}\right) = 0$ and $\omega \left({\tau }_{2,l}^{*}\right) = {\omega }_{l}$. By differentiating (10) with respect to ${\tau }_{2}$, we get

or

From Eq (10) we have ${e}^{-\lambda {\tau }_{2}} = -\frac{{\lambda }^{5}+{r}_{1}{\lambda }^{4}+{z}_{1}{\lambda }^{3}+{z}_{2}{\lambda }^{2}+{z}_{3}\lambda +{z}_{4}}{{r}_{5}{\lambda }^{4}+{r}_{6}{\lambda }^{3}+{z}_{5}{\lambda }^{2}+{z}_{6}\lambda +{z}_{7}}$ such that we have

For ${\tau }_{2} = {\tau }_{2,l}^{*}$, we obtain

or

Let

$A = -\frac{5{{\omega }_{l}}^{5}-3{z}_{1}{{\omega }_{l}}^{3}+{z}_{3}{\omega }_{l}+i(2{z}_{2}{\omega }_{l}^{2}-4{r}_{1}{{\omega }_{l}}^{4})}{\left({{-\omega }_{l}}^{5}+{{z}_{1}{{\omega }_{l}}^{3}-{z}_{3}{\omega }_{l}+i(r}_{1}{{\omega }_{l}}^{4}+-{z}_{2}{{\omega }_{l}}^{2}+{z}_{4}\right)}$ and $B = \frac{{z}_{6}{\omega }_{l}-3{r}_{6}{{\omega }_{l}}^{3}i(2{z}_{5}{{\omega }_{l}}^{2}-4{r}_{5}{{\omega }_{l}}^{4})}{{r}_{6}{{\omega }_{l}}^{3}-{z}_{6}{\omega }_{l}+{i(r}_{5}{{\omega }_{l}}^{4}-{z}_{5}{{\omega }_{l}}^{2}+{z}_{7})}$.

By multiplying A and B with their respectively conjugate of the denominator, we obtain

By rewritten Eq (21) as ${\left(\frac{d\lambda }{d{\tau }_{2}}\right)}_{{\tau }_{2} = {\tau }_{2,l}^{*}}^{-1} = \frac{1}{{{\omega }_{l}}^{2}}\left[A+B+i{\tau }_{2,l}^{*}{\omega }_{l}\right],$ we have

For $\omega = {\omega }_{l},$ from Eq (19) we have

such that

with

Consider Eq (20). By differentiating Eq (20) with respect to $X = {{\omega }_{l}}^{2}$, we get

By substituting (23) into (22) we obtain

then we get

Summarizing our analysis results for this case, we have the following theorem.

Theorem 1. Suppose that the conditions $\left({H}_{1}\right)$ and $\left({H}_{2}\right)$ are fulfilled. Let ${{\tau }_{2}}^{*} = \underset{l = 1,\cdots ,5;j\in \mathbb{N}}{\mathit{min}}\left\{{\tau }_{2,l}^{1,j},{\tau }_{2,l}^{2,j}\right\}$. Then when ${\tau }_{2} < {{\tau }_{2}}^{*},$ the equilibrium point ${E}^{*}$is locally asymptotically stable. Furthermore, if ${F}^{\text{'}}\left({{\omega }_{l}}^{2}\right) > 0$ then a Hopf bifurcation occurs at ${E}^{*}$ when ${\tau }_{2} = {{\tau }_{2}}^{*}$.

Case 3: ${\mathit{\tau }}_{\mathit{i}} > 0,i = \mathrm{1,2},3$

For this case, we assume that a new type of pathogen has entered the host body. It is assumed that there has been no previous contact between the pathogen and the host immune system. As a result, the body cannot directly respond to pathogens such that we have ${\tau }_{1},{\tau }_{2},{\tau }_{3} > 0$.

Lemma 2. If all roots of Eq (10) have negative real parts for ${\tau }_{2} > 0$, then there exists a ${{\tau }_{3}}^{*}\left({\tau }_{2}\right) > 0$ such that all roots of Eq (10) have negative real parts when ${\tau }_{3} < {{\tau }_{3}}^{*}\left({\tau }_{2}\right)$.

Proof:

For ${\tau }_{2} > 0$, we assume that Eq (10) has no roots with nonnegative real part. Then Eq. (10) with ${\tau }_{3} = 0$ and ${\tau }_{2} > 0$ has no root with nonnegative real part. Suppose ${\tau }_{3} > 0$ then Eq (8) is analytic in $\lambda$ and ${\tau }_{3}$. Furthermore, as ${\tau }_{3}$ increases, the stability of the equilibrium point will be lost. Since Eq (8) with ${\tau }_{3} = 0$ has no root with nonnegative real part then there exists ${\tau }_{3}^{\mathrm{*}}\left({\tau }_{2}\right) > 0$ such that all roots of Eq (8) with ${\tau }_{3} < {\tau }_{3}^{\mathrm{*}}\left({\tau }_{2}\right)$ have negative real part.

Theorem 2. Suppose that the conditions $\left({H}_{1}\right)$ and $\left({H}_{2}\right)$ are fulfilled. Let ${{\tau }_{2}}^{*} = \underset{l = 1,\cdots ,5;j\in \mathbb{N}}{\mathit{min}}\left\{{\tau }_{2,l}^{1,j},{\tau }_{2,l}^{2,j}\right\}$. Then there exists a ${\tau }_{3}^{\mathrm{*}}\left({\tau }_{2}\right) > 0$, for any ${\tau }_{2}\in [0,{\tau }_{2}^{\mathrm{*}})$, such that equilibrium point ${E}^{*}$is locally asymptotically stable when ${\tau }_{3}\in [0,{\tau }_{3}^{\mathrm{*}}\left({\tau }_{2}\right))$.

Proof:

Using Theorem 1, we have that the roots of Eq (10) have negative real parts for ${\tau }_{2}\in [0,{\tau }_{2}^{\mathrm{*}})$. Using Lemma 2, we conclude that the equilibrium point ${\mathrm{E}}^{\mathrm{*}}$is locally asymptotically stable when ${\mathrm{\tau }}_{3}\in [0,{\mathrm{\tau }}_{3}^{\mathrm{*}}\left({\mathrm{\tau }}_{2}\right))$.

4.   Numerical simulation

We present a numerical simulation for the system by assuming that the pathogen invades the host for the first time so that no memory cells are preformed. As a result, there is a certain time for the immune system to respond by producing primary immune cells along with the formation of memory cells after pathogen injection, i.e., ${\tau }_{1},{\tau }_{2},{\tau }_{3} > 0$. The delays are chosen at the critical delays, the conditions for the existence of oscillatory behavior, i.e., ${\tau }_{2} = 2.64$ days and ${\tau }_{3} = 7.776$ days while ${\tau }_{1}$ is assumed 1.2 days. The parameter values for the numerical simulations are presented in Table 1. The simulations are presented in two cases, with and without the addition of therapeutic protein. This is intended to study how the therapeutic proteins affect the dynamics of pathogens and immune cells systems. The initial values and the historical functions are chosen as follows: $M\left(0\right) = 0,D\left(0\right) = 0,P\left(t\right) = 4x{10}^{-3}, \forall t\in \left[-{\tau }_{1},0\right],{I}_{S}\left(t\right) = 6x{10}^{-3}, \forall t\in \left[-{\tau }_{2},0\right],{I}_{p}\left(t\right) = 3.5x{10}^{-1},$ $\forall t\in \left[-{\tau }_{3},0\right]$. These initial conditions show that the immune system could not destroy the pathogen such that assistance in the addition of the therapeutic proteins was needed.

We first present a simulation of the system without the addition of the protein drug to study how the state of the system was initially, and how the protein drug affects the dynamics of the system. The simulation results are given in Figure 2. We can observe that the number of pathogens in the body increase as time increase in days without intervention of the drug therapy. It means that the body's immune system cannot reduce the number of pathogens by itself due to the highest of pathogenesis in the system. This is also due to the existence of delays in producing primary immune cells and maturing specific immune cells. As shown in Figure 2-a, the pathogen will continue to multiply, increase exponentially, and converge to its carrying capacity even though occasionally it decreases due to the interactions with the specific immune cells. Around 17 days, the pathogen reaches a saturation point and does not increase again. In Figure 2-b, we can observe that the number of primary immune cells increases as a result of pathogen infection. However, primary cells cannot stop the invasion of pathogens due to their inefficiency. Since there exists a delay in the maturation of primary cells into specific cells, the existing specific cells cannot avoid the increasing of pathogens due to the higher number of pathogenic replications compared to the number of the existing specific cells. In Figure 2-c, we show that the oscillatory behavior that appears in the system is clearly observed in the dynamic of the specific immune cells which means that natural delays greatly affect the number of specific immune cells that play an important role in fighting pathogens. The specific immune cell solution shows a fluctuation at the beginning of time observation. However, for a long time, the number of specific immune cells goes to zero caused by the high pathogenesis of the pathogen. Even though the host body produces many immune cells, the cells will die after their interaction with the pathogens. On the other hand, initially, the number of the memory cells is very small due to the existence of a time delay in the generation of the new memory cells. After the delay, the number of memory cells increases as time increase and converge to its stationary point.

After observing the condition of the system without protein drug intervention, we next present the simulation results for the system with therapeutic protein intervention. Figure 3 shows the dynamics of the pathogen which is decreases as time increase. It decreases from the initial observation time until it reaches zero at 18 days. This behavior is quite different from the previous simulation in which the number of pathogens increases as time increase. This means that the addition of therapeutic proteins has succeeded in reducing the invasion of pathogens and even eliminating pathogens in the host's body. Due to the decreasing of the pathogen, the number of primary immune cells is also decreasing at the final time. The host body will reduce the production of primary immune cells when the number of pathogens decrease. In addition, the presence of protein drugs also affects the dynamics of specific immune cells where the addition of therapeutic proteins has reduced fluctuations in specific immune solutions with relatively small oscillatory effects. For the memory cells, the addition of protein drug increases the number of memory cells as a result of the initial increase in primary cells that triggered by the presence of protein drug.

5.   Conclusions

A modified mathematical model of immune system was proposed by considering the immune long pathway as a simple lumped pathway. It was assumed that the interaction of pathogens and immune cells is a competitive interaction, and under the influences of therapeutic protein, the production of primary immune cells was triggered. Some natural discrete delays were introduced into the model to accommodate the slow response of primary cells to the pathogen attack, and to consider the slow maturation rate of the primary cells. We found that there exists a pair of critical delays for which the system appearance oscillatory behavior. However, the presence of natural delays has the greatest effect on the dynamics of specific immune cells, i.e. immune cells that are responsible for attacking pathogens directly. It was numerically observed that the fluctuation of specific immune cells was quite high when the protein drug was not injected into the host body. Moreover, the number of pathogens increases when the number of specific immune cells decreases, and converges to their stationary point when the number of specific immune cells converges to zero. After the addition of therapeutic proteins, the invasion of pathogens can be reduced or even eliminated from the system, and protein drugs were also successful in reducing the fluctuations of specific immune cells. Based on the presented numerical simulations, it concluded that the proposed model can capture the dynamic of the delayed system with or without drug intervention. The existence of therapeutic protein was successfully minimizing the number of pathogens and accelerated the elimination of pathogens in the host body. It also observed that the existence of therapeutic protein minimized the oscillation in the system especially in the specific immune cells dynamic. It meant that the protein drugs successfully speed up the specific immune performance in eliminating pathogens invasion in the host body.

Acknowledgments

This work was supported by the Institute for Research and Community Service, Hasanuddin University [grant number 1516/UN4.22/PT.01.03/2020]. We thank the anonymous reviewers for their valuable suggestions helped improve this manuscript.

Conflict of interest

The authors declare that there is no conflict of interest regarding the publication of this paper.