Research article

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

  • Received: 19 September 2021 Revised: 30 November 2021 Accepted: 03 December 2021 Published: 14 February 2022
  • MSC : 92F05, 34K20, 34K18

  • This study aims to propose and analyze a mathematical model of the competitive interaction of the pathogen-immune system. Some effects of the existence of natural delays and the addition of therapeutic proteins are considered in the model. A delay arises from the indirect response of the host body when a pathogen invades. The other comes from the maturation of immune cells to produce immune memory cells since the immune system and antigenic substances responsible for provoking the production of immune memory cells. Analytical investigations suggest several sufficient conditions for the existence of a positive steady-state solution. There is a critical pair of delays at which oscillatory behavior appears around the positive steady-state solution. Numerical simulations were carried out to describe the results of the analysis and show that the proposed model can describe the speed of pathogen eradication due to the addition of therapeutic proteins as antigenic substances.

    Citation: Kasbawati, Yuliana Jao, Nur Erawaty. Dynamic study of the pathogen-immune system interaction with natural delaying effects and protein therapy[J]. AIMS Mathematics, 2022, 7(5): 7471-7488. doi: 10.3934/math.2022419

    Related Papers:

    [1] Kasbawati, Mariani, Nur Erawaty, Naimah Aris . A mathematical study of effects of delays arising from the interaction of anti-drug antibody and therapeutic protein in the immune response system. AIMS Mathematics, 2020, 5(6): 7191-7213. doi: 10.3934/math.2020460
    [2] Jinting Lin, Changjin Xu, Yiya Xu, Yingyan Zhao, Yicheng Pang, Zixin Liu, Jianwei Shen . Bifurcation and controller design in a 3D delayed predator-prey model. AIMS Mathematics, 2024, 9(12): 33891-33929. doi: 10.3934/math.20241617
    [3] Wei Ou, Changjin Xu, Qingyi Cui, Yicheng Pang, Zixin Liu, Jianwei Shen, Muhammad Zafarullah Baber, Muhammad Farman, Shabir Ahmad . Hopf bifurcation exploration and control technique in a predator-prey system incorporating delay. AIMS Mathematics, 2024, 9(1): 1622-1651. doi: 10.3934/math.2024080
    [4] Qinghui Liu, Xin Zhang . Chaos detection in predator-prey dynamics with delayed interactions and Ivlev-type functional response. AIMS Mathematics, 2024, 9(9): 24555-24575. doi: 10.3934/math.20241196
    [5] Erxi Zhu . Time-delayed feedback control for chaotic systems with coexisting attractors. AIMS Mathematics, 2024, 9(1): 1088-1102. doi: 10.3934/math.2024053
    [6] Yingyan Zhao, Changjin Xu, Yiya Xu, Jinting Lin, Yicheng Pang, Zixin Liu, Jianwei Shen . Mathematical exploration on control of bifurcation for a 3D predator-prey model with delay. AIMS Mathematics, 2024, 9(11): 29883-29915. doi: 10.3934/math.20241445
    [7] San-Xing Wu, Xin-You Meng . Dynamics of a delayed predator-prey system with fear effect, herd behavior and disease in the susceptible prey. AIMS Mathematics, 2021, 6(4): 3654-3685. doi: 10.3934/math.2021218
    [8] Jorge Luis Ramos-Castellano, Miguel Angel Dela-Rosa, Iván Loreto-Hernández . Bifurcation analysis for the coexistence in a Gause-type four-species food web model with general functional responses. AIMS Mathematics, 2024, 9(11): 30263-30297. doi: 10.3934/math.20241461
    [9] Ahmed M. Elaiw, Amani S. Alsulami, Aatef D. Hobiny . Global properties of delayed models for SARS-CoV-2 infection mediated by ACE2 receptor with humoral immunity. AIMS Mathematics, 2024, 9(1): 1046-1087. doi: 10.3934/math.2024052
    [10] Xin-You Meng, Fan-Li Meng . Bifurcation analysis of a special delayed predator-prey model with herd behavior and prey harvesting. AIMS Mathematics, 2021, 6(6): 5695-5719. doi: 10.3934/math.2021336
  • This study aims to propose and analyze a mathematical model of the competitive interaction of the pathogen-immune system. Some effects of the existence of natural delays and the addition of therapeutic proteins are considered in the model. A delay arises from the indirect response of the host body when a pathogen invades. The other comes from the maturation of immune cells to produce immune memory cells since the immune system and antigenic substances responsible for provoking the production of immune memory cells. Analytical investigations suggest several sufficient conditions for the existence of a positive steady-state solution. There is a critical pair of delays at which oscillatory behavior appears around the positive steady-state solution. Numerical simulations were carried out to describe the results of the analysis and show that the proposed model can describe the speed of pathogen eradication due to the addition of therapeutic proteins as antigenic substances.



    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 [24]. 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 [59].

    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 [1117]. 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 [1115,17]. Some researchers are also concerned to study and observe the effects of delays occurring in the immune system [1829]. 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.

    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, Is which states the number of specific immune cells, Ip 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:

    Figure 1.  Compartment diagram of pathogen responses to the therapeutic protein addition to the host immune system.

    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 ω. 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 σ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 μ.

    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 τ1, where c is the rate of cellular and biochemical reactions in the body and λ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 λMand 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 τ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 α is the rate of regeneration of immune cells and β is the rate of reduction of immune cells because it loses against pathogens. The number of immune cells will increase if α > β, namely the rate of regeneration is greater than the rate of competition with pathogens and will decrease if β > α.

    4. The production of memory cells usually takes place 14 days after the first pathogen invasion, therefore a delay of τ3 occurs. Memory cells will be produced together with primary immune production at the rate of γ and decay at the rate of δ. 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 λ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:

    dP(t)dt=P(t)[abP(t)]ωP(t)(σμ)P(t)Is(t),
    dIp(t)dt=c[λ0P(tτ1)+λMM(t)P(t)Ip+kD(t)Ip(t)],
    dIs(t)dt=cIs(tτ2)(βα)P(t)Is(t), (1)
    dM(t)dt=γIp(tτ3)δM(t),
    dD(t)dt=λDkD(t)Ip(t),

    with initial conditions M(0)=M0>0,D(0)=D0>0,and the historical functions: P(t)=β1,t[τ1,0],IS(t)=β2,t[τ2,0],Ip(t)=β3,t[τ3,0], for β1,β2,β3R+.

    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, R5+={(P,Ip,Is,M,D)R5|P0,Ip0,Is0,M0,D0}. Since the vector field does not point to the exterior of R5+ then Eq (1) is defined in the region R5+. Steady state of (1) is the stationary solution E=(P,Ip,Is,M,D) that fulfills,

    P[abP]ωP(σμ)PIs=0, (2)
    c[λ0P+λMMPIp+kDIp]=0, (3)
    cIs(βα)PIs=0, (4)
    γIpδM=0,
    λDkDIp=0. (5)

    With several mathematical processes we have the following results.

    Proposition 1. System (1) has a nontrivial positive steady state solution E=(P,Ip,Is,M,D) with P=cβα, Ip=δ(cλ0+λD(βα))δ(βα)λMcγ, Is=(aω)(βα)bc(βα)(σμ), M=γ(cλ0+λD(βα))δ(βα)λMcγ, and D=λD(δ(βα)λMcγ)kδ(cλ0+λD(βα)) if it fulfills conditions β>α, δ(βα)>λMcγ, (aω)(βα)>bc and σ>μ, or (aω)(βα)<bc and σ<μ.

    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,

    dPdt=f1(P,Ip,Is,M,D,Pτ1,Isτ2,Ipτ3),
    dIpdt=f2(P,Ip,Is,M,D,Pτ1,Isτ2,Ipτ3),
    dIsdt=f3(P,Ip,Is,M,D,Pτ1,Isτ2,Ipτ3), (6)
    dMdt=f4(P,Ip,Is,M,D,Pτ1,Isτ2,Ipτ3),
    dDdt=f5(P,Ip,Is,M,D,Pτ1,Isτ2,Ipτ3).

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

    [˙P(t)˙Ip(t)˙Is(t)˙M(t)˙D(t)]=J0[P(t)Ip(t)Is(t)M(t)D(t)]+J1[P(tτ1)Ip(tτ1)Is(tτ1)M(tτ1)D(tτ1)]+J2[P(tτ2)Ip(tτ2)Is(tτ2)M(tτ2)D(tτ2)]+J3[P(tτ3)Ip(tτ3)Is(tτ3)M(tτ3)D(tτ3)], (7)

    where

    J0=[Δ10Δ400Δ2Δ50Δ8Δ9Δ30Δ700000δ00Δ600Δ10],J1=[00000cλ00000000000000000000],J2=[000000000000c000000000000],

    and J3=[0000000000000000γ00000000],

    with

    Δ1=a2bPω(σμ)Is,
    Δ2=cλMM,
    Δ3=(αβ)Is,
    Δ4=(σμ)P,
    Δ5=c+ckD,
    Δ6=kD,
    Δ7=(αβ)P,
    Δ8=cλMP,
    Δ9=ckIp,
    Δ10=kIp.

    The characteristic equation associated with (7) is given by

    λ5+r1λ4+r2λ3+r3λ2+r4λ+(r5λ4+r6λ3+r7λ2+r8λ+r9)eλτ2+(r10λ3+r11λ2+r12λ+r13)eλτ3+(r14λ2+r15λ+r16)eλ(τ2+τ3)+r17=0, (8)

    with

    r1=δΔ1Δ5Δ7,
    r2=δ(Δ1+Δ5+Δ7+Δ10)+Δ1(Δ5+Δ7+Δ10)+Δ5(Δ7+Δ10)Δ3Δ4Δ6Δ9+Δ7Δ10,
    r3=δ(Δ1Δ5+Δ1Δ7+Δ1Δ10Δ3Δ4+Δ5Δ7+Δ5Δ10Δ6Δ9+Δ7Δ10)Δ1Δ5(Δ7+Δ10)+Δ7Δ10(Δ5Δ1)+Δ6Δ9(Δ1+Δ7)+Δ3Δ4(Δ5+Δ10),
    r4=δ(Δ1Δ5Δ7+Δ1Δ5Δ10Δ1Δ6Δ9+Δ1Δ7Δ10Δ3Δ4Δ5Δ3Δ4Δ10+Δ5Δ7Δ10Δ6Δ7Δ9)+(Δ5Δ10Δ6Δ9)(Δ1Δ7Δ3Δ4),
    r5=c,
    r6=c(Δ1+Δ5+Δ10δ),
    r7=cδ(Δ1+Δ5+Δ10)c(Δ1Δ5+Δ1Δ10+Δ5Δ10Δ6Δ9),
    r8=cδ(Δ1Δ5+Δ1Δ10+Δ5Δ10Δ6Δ9)+c(Δ1Δ5Δ10Δ1Δ6Δ9),
    r9=cδ(Δ1Δ5Δ10Δ1Δ6Δ9),
    r10=γΔ8,
    r11=γ(Δ1Δ8Δ7Δ8+Δ8Δ10),
    r12=γ(Δ1Δ7Δ8Δ1Δ8Δ10+Δ3Δ4Δ8Δ7Δ8Δ10),
    r13=γΔ8Δ10(Δ1Δ7Δ3Δ4),
    r14=cγΔ8,
    r15=cγΔ8(Δ1+Δ10),
    r16=cγΔ1Δ8Δ10,
    r17=δΔ1Δ7(Δ5Δ10Δ6Δ9)δΔ3Δ4(Δ5Δ10Δ6Δ9).

    Equation (8) is a five-degree exponential polynomial in λ. 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., τ2 and τ3 with free τ1 meaning that stability of E is only influenced by τ2 and τ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 λ. Therefore, analysis of the roots of Eq (8) will be presented in some cases as follow.

    Case 1: ττ2=ττ3=0 with τ1>0

    Assume that τ2=τ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,

    λ5+a1λ4+a2λ3+a3λ2+a4λ+a5=0, (9)

    with

    a1=r1+r5,
    a2=r2+r6+r10,
    a3=r3+r7+r11+r13,
    a4=r4+r8+r12+r15,
    a5=r9+r14+r16+r17.

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

    D1=|a1|>0,
    D2=|a1a31a2|>0ora1a2>a3,
    D3=|a1a3a51a2a40a1a3|>0ora1a2a3+a1a5>a32+a12a4,
    D4=|a1a3a501a2a400a1a3a501a2a4|>0or
    a1a2a3a4+2a1a4a5+a2a3a5>a12a42+a32a4+a52+a1a22a5,
    D5=|a1a3a5001a2a4000a1a3a5001a2a4000a1a3a5|>0or
    a5(a1a2a3a4)+a25(2a1a4a5+a2a3a5)>a5(a12a42+a32a4+a52+a1a22a5).

    Case 2: ττ1,ττ2>0,ττ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 τ3=0). It means that the body can quickly respond the invasion of pathogen such that infection occurred (it also can make τ10). Then Eq (8) becomes,

    λ5+r1λ4+z1λ3+z2λ2+z3λ+z4+(r5λ4+r6λ3+z5λ2+z6λ+z7)eλτ2=0, (10)

    with

    z1=r2+r10,
    z2=r3+r11,
    z3=r4+r12,
    z4=r13+r17,
    z5=r7+r14,
    z6=r8+r15,
    z7=r9+r16.

    Since Eq (10) will become Eq (9) when τ2=0, then based on the Proposition 2, the equilibrium point E is locally asymptotically stable. When τ2>0, the roots of Eq (10) are not easy to obtain explicitly. The eigenvalues of Eq (10) will depend on τ2. Suppose that the eigenvalue of Eq (10) is in a complex form, λ(τ2)=v(τ2)±iω(τ2) with τ2>0. To find out whether E is stable when τ2>0or has a limit cycle, the analysis will be processed by assuming the eigenvalue of Eq (10) is in an imaginary form, λ=±iω,ωR. If iω is a purely imaginary root of (10), then it fulfills

    (iω)5+r1(iω)4+z1(iω)3+z2(iω)2+z3(iω)+z4+(r5(iω)4+r6(iω)3+z5(iω)2+z6(iω)+z7)eiωτ2=0.

    Simplifying the equation, we get

    [r1ω4z2ω2+z4+r5ω4cos(ωτ2)r6ω3sin(ωτ2)z5ω2cos(ωτ2)+z6ωsin(ωτ2)+z7cos(ωτ2)]+i[ω5z1ω3+z3ωr5ω4sin(ωτ2)r6ω3cos(ωτ2)+z5ω2sin(ωτ2)+z6ωcos(ωτ2)z7sin(ωτ2)]=0.

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

    r1ω4z2ω2+z4+r5ω4cos(ωτ2)r6ω3sin(ωτ2)z5ω2cos(ωτ2)+z6ωsin(ωτ2)+z7cos(ωτ2)=0, (11)

    and

    ω5z1ω3+z3ωr5ω4sin(ωτ2)r6ω3cos(ωτ2)+z5ω2sin(ωτ2)+z6ωcos(ωτ2)z7sin(ωτ2)=0. (12)

    From Eq (11) we have

    sin(ωτ2)=r1ω4z2ω2+z4+cos(ωτ2)(r5ω4z5ω2+z7)(r6ω3z6ω), (13)
    cos(ωτ2)=r1ω4z2ω2+z4+sin(ωτ2)(z6ωr6ω3)(z5ω2r5ω4z7). (14)

    While from Eq (12) we have

    sin(ωτ2)=ω5z1ω3+z3ω+cos(ωτ2)(z6ωr6ω3)(r5ω4z5ω2+z7), (15)
    cos(ωτ2)=ω5z1ω3+z3ω+sin(ωτ2)(z5ω2r5ω4z7)(r6ω3z6ω). (16)

    From Eqs (13) and (15) we get

    cos(ωτ2)=s1s3s2s4s21+s22, (17)

    with s1=r6ω3z6ω, s2=r5ω4z5ω2+z7, s3=ω5z1ω3+z3ω, and s4=r1ω4z2ω2+z4. From Eqs (14) and (16) we get

    sin(ωτ2)=s1s4+s2s3s21+s22. (18)

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

    ω10+(r12r522z1)ω8+(2r5z52r1z2+2z3+z12r62)ω6+(2r1z4+2r6z62r5z72z1z3+z22z52)ω4+(2z5z72z2z4+z32z62)ω2+(z42z72)=0. (19)

    Consider p1=r12r522z1,

    p2=2r5z52r1z2+2z3+z12r62,
    p3=2r1z4+2r6z62r5z72z1z3+z22z52,
    p4=2z5z72z2z4+z32z62,
    p5=z42z72,

    Then Eq (19) can be written as

    ω10+p1ω8+p2ω6+p3ω4+p4ω2+p5=0.

    If we consider X=ω2 then we have

    F(X)=X5+p1X4+p2X3+p3X2+p4X+p5=0. (20)

    Assume

    (H2) p5<0, p4>0, p3<0, p2>0, and p1<0.

    If condition H2 is fulfilled then all roots of F(X) are positive, say Xl,l=1,2,5 such that ωl=Xl. Therefore condition H2 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 dRe(λ)dτ2.

    Now, consider Eqs (17) and (18). For ω=ωl, we have the solutions of both equations,

    ωlτ2=arccos[(r6ωl3z6ωl)(ωl5z1ωl3+z3ωl)(r5ωl4z5ωl2+z7)(r1ωl4z2ωl2+z4)(r6ωl3z6ωl)2+(r5ωl4z5ωl2+z7)2]

    if sin(ωlτ2)>0, and

    ωlτ2=2πarccos[(r6ωl3z6ωl)(ωl5z1ωl3+z3ωl)(r5ωl4z5ωl2+z7)(r1ωl4z2ωl2+z4)(r6ωl3z6ωl)2+(r5ωl4z5ωl2+z7)2]

    if sin(ωlτ2)0. If we define two sequences {τ1,j2,l} and {τ2,j2,l} for l=1,,5 and jN then we have

    τ1,j2,l=1ωl[arccos((r6ωl3z6ωl)(ωl5z1ωl3+z3ωl)(r5ωl4z5ωl2+z7)(r1ωl4z2ωl2+z4)(r6ωl3z6ωl)2+(r5ωl4z5ωl2+z7)2)+2jπ],
    τ2,j2,l=1ωl[2πarccos((r6ωl3z6ωl)(ωl5z1ωl3+z3ωl)(r5ωl4z5ωl2+z7)(r1ωl4z2ωl2+z4)(r6ωl3z6ωl)2+(r5ωl4z5ωl2+z7)2)+2jπ].

    Lemma 1. Suppose that τ2,l is an element of either the sequence {τ1,j2,l} or {τ2,j2,l} with ωl. Then the characteristic equation (10) has a pair of conjugate pure imaginary roots λ=±iωl that fulfils

    sign{dRe(λ)dτ2|τ2=τ2,l}=sign{F'(ωl2)}.

    Proof:

    Consider the roots of (10), λ(τ2)=v(τ2)+iω(τ2) where v(τ2,l)=0 and ω(τ2,l)=ωl. By differentiating (10) with respect to τ2, we get

    [5λ4+4r1λ3+3z1λ2+2z2λ+z3+(4r5λ3+3r6λ2+2z5λ+z6τ2(r5λ4+r6λ3+z5λ2+z6λ+z7))eλτ2]dλdτ2=λ(r5λ4+r6λ3+z5λ2+z6λ+z7)eλτ2

    or

    (dλdτ2)1=τ2λ+5λ4+4r1λ3+3z1λ2+2z2λ+z3eλτ2λ(r5λ4+r6λ3+z5λ2+z6λ+z7)+(4r5λ3+3r6λ2+2z5λ+z6)λ(r5λ4+r6λ3+z5λ2+z6λ+z7).

    From Eq (10) we have eλτ2=λ5+r1λ4+z1λ3+z2λ2+z3λ+z4r5λ4+r6λ3+z5λ2+z6λ+z7 such that we have

    (dλdτ2)1=τ2λ5λ4+4r1λ3+3z1λ2+2z2λ+z3λ(λ5+r1λ4+z1λ3+z2λ2+z3λ+z4)+(4r5λ3+3r6λ2+2z5λ+z6)λ(r5λ4+r6λ3+z5λ2+z6λ+z7).

    For τ2=τ2,l, we obtain

    (dλdτ2)1τ2=τ2,l=τ2,l(iωl)5(iωl)4+4r1(iωl)3+3z1(iωl)2+2z2(iωl)+z3(iωl)((iωl)5+r1(iωl)4+z1(iωl)3+z2(iωl)2+z3(iωl)+z4)
    +(4r5(iωl)3+3r6(iωl)2+2z5(iωl)+z6)(iωl)(r5(iωl)4+r6(iωl)3+z5(iωl)2+z6(iωl)+z7),

    or

    (dλdτ2)1τ2=τ2,l=1ωl2[iτ2,lωl5ωl5i4r1ωl43z1ωl3+i2z2ω2l+z3ωl(ωl5+ir1ωl4+z1ωl3iz2ωl2z3ωl+iz4)+i4r5ωl43r6ωl3+i2z5ωl2+z6ωlir5ωl4+r6ωl3iz5ωl2z6ωl+iz7]. (21)

    Let

    A=5ωl53z1ωl3+z3ωl+i(2z2ω2l4r1ωl4)(ωl5+z1ωl3z3ωl+i(r1ωl4+z2ωl2+z4) and B=z6ωl3r6ωl3i(2z5ωl24r5ωl4)r6ωl3z6ωl+i(r5ωl4z5ωl2+z7).

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

    Re(A)=5ωl53z1ωl3+z3ωl+i(2z2ω2l4r1ωl4)(ωl5+z1ωl3z3ωl+i(r1ωl4+z2ωl2+z4).(ωl5+z1ωl3z3ωli(r1ωl4+z2ωl2+z4)(ωl5+z1ωl3z3ωli(r1ωl4+z2ωl2+z4)
    =5ω10l+(4r128z1)ω8l+(6r1z2+3z12+6z3)ω6l+(4r1z44z1z3+2z22)ω4l+(2z2z4+z32)ω2lω10l+(r122z1)ω8l+(z122r1z2+2z3)ω6l+(2r1z42z1z3+z22)ω4l+(z322z2z4)ω2l+z42,
    Re(B)=z6ωl3r6ωl3+i(2z5ωl24r5ωl4)r6ωl3z6ωl+i(r5ωl4z5ωl2+z7).r6ωl3z6ωli(r5ωl4z5ωl2+z7)r6ωl3z6ωli(r5ωl4z5ωl2+z7)
    =4r52ω8l+(6r5z53r62)ω6l+(4r5z7+4r6z62z52)ω4l+(2z5z7z62)ω2lr52ω8l+(2r5z5+r62)ω6l+(2r5z72r6z6+z52)ω4l+(2z5z7+z62)ω2l+z72.

    By rewritten Eq (21) as (dλdτ2)1τ2=τ2,l=1ωl2[A+B+iτ2,lωl], we have

    Re[(dλdτ2)1τ2=τ2,l]=1ωl2[Re(A)+Re(B)]
    =1ωl2[4r52ω8l+(6r5z53r62)ω6l+(4r5z7+4r6z62z52)ω4l+(2z5z7z62)ω2lr52ω8l+(2r5z5+r62)ω6l+(2r5z72r6z6+z52)ω4l+(2z5z7+z62)ω2l+z72]
    +5ω10l+(4r128z1)ω8l+(6r1z2+3z12+6z3)ω6l+(4r1z44z1z3+2z22)ω4l+(2z2z4+z32)ω2lω10l+(r122z1)ω8l+(z122r1z2+2z3)ω6l+(2r1z42z1z3+z22)ω4l+(z322z2z4)ω2l+z42.

    For ω=ωl, from Eq (19) we have

    ω10l+(r122z1)ω8l+(z122r1z2+2z3)ω6l+(2r1z42z1z3+z22)ω4l+(z322z2z4)ω2l+z42=r52ω8l+(2r5z5+r62)ω6l+(2r5z72r6z6+z52)ω4l+(2z5z7+z62)ω2l+z72,

    such that

    Re[(dλdτ2)1τ2=τ2,l]=ϕr52ω8l+(2r5z5+r62)ω6l+(2r5z72r6z6+z52)ω4l+(2z5z7+z62)ω2l+z72, (22)

    with

    ϕ=[5ω8l+(4r124r528z1)ω6l+(6r5z56r1z2+6z3+3z123r62)ω4l+(4r1z44z1z34r5z7+4r6z62z52+2z22)ω2l+(2z5z72z2z4+z32z62)].

    Consider Eq (20). By differentiating Eq (20) with respect to X=ωl2, we get

    F'(ωl2)=5ωl8+4p1ωl6+3p2ωl4+2p3ωl2+p4. (23)

    By substituting (23) into (22) we obtain

    Re[(dλdτ2)1τ2=τ2,l]=F'(ωl2)(r5ωl4z5ωl2+z7)2+ω2l(r6ωl2z6)2.Sincesign{Re[(dλdτ2)1τ2=τ2,l]}=sign{dRe(λ)dτ2|τ2=τ2,l},

    then we get

    sign{dRe(λ)dτ2|τ2=τ2,l}=sign{F'(ωl2)}.

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

    Theorem 1. Suppose that the conditions (H1) and (H2) are fulfilled. Let τ2=minl=1,,5;jN{τ1,j2,l,τ2,j2,l}. Then when τ2<τ2, the equilibrium point Eis locally asymptotically stable. Furthermore, if F'(ωl2)>0 then a Hopf bifurcation occurs at E when τ2=τ2.

    Case 3: τi>0,i=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 τ1,τ2,τ3>0.

    Lemma 2. If all roots of Eq (10) have negative real parts for τ2>0, then there exists a τ3(τ2)>0 such that all roots of Eq (10) have negative real parts when τ3<τ3(τ2).

    Proof:

    For τ2>0, we assume that Eq (10) has no roots with nonnegative real part. Then Eq. (10) with τ3=0 and τ2>0 has no root with nonnegative real part. Suppose τ3>0 then Eq (8) is analytic in λ and τ3. Furthermore, as τ3 increases, the stability of the equilibrium point will be lost. Since Eq (8) with τ3=0 has no root with nonnegative real part then there exists τ3(τ2)>0 such that all roots of Eq (8) with τ3<τ3(τ2) have negative real part.

    Theorem 2. Suppose that the conditions (H1) and (H2) are fulfilled. Let τ2=minl=1,,5;jN{τ1,j2,l,τ2,j2,l}. Then there exists a τ3(τ2)>0, for any τ2[0,τ2), such that equilibrium point Eis locally asymptotically stable when τ3[0,τ3(τ2)).

    Proof:

    Using Theorem 1, we have that the roots of Eq (10) have negative real parts for τ2[0,τ2). Using Lemma 2, we conclude that the equilibrium point Eis locally asymptotically stable when τ3[0,τ3(τ2)).

    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., τ1,τ2,τ3>0. The delays are chosen at the critical delays, the conditions for the existence of oscillatory behavior, i.e., τ2=2.64 days and τ3=7.776 days while τ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(0)=0,D(0)=0,P(t)=4x103,t[τ1,0],IS(t)=6x103,t[τ2,0],Ip(t)=3.5x101, t[τ3,0]. 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.

    Table 1.  Definition and value of model parameters.
    Parameter Definition Value and Unit
    a Pathogen replication rate 0.02 mol/days
    b Carrying capacity of the host body environment 0.0035 mol/days
    ω Pathogen natural death rate 0.01 mol/days
    σ Pathogen death rate due to interactions with immune cells 2.5 × 101 mol/days
    μ Pathogen evolution rate 0.0015 mol/days
    c The rate of cellular and biochemical reactions in the body 5.073 mol/days
    λ0 Immunity level from the burden of infection 1.5 mol/days
    λM Immune cell regeneration rate 1.467 × 10−3 mol/days
    τ1 Time delay for immune cell regeneration rate (after pathogen invasion) 1.2 days
    τ2 Critical delay for maturation of immune cells 2.64 days
    τ3 Critical delay for formation of memory cells 7.776 days
    k Therapeutic protein drug binding rate by immune cells 36 mol/days
    α Immune cell regeneration rate 1 × 103 mol/days
    β Immune cell reduction rate 3.8 × 103 mol/days
    γ Immune cell production rate 2.3 × 10−3 mol/days
    δ Immune cell decay rate 0.2 mol/days
    λD Protein drug input dose rate 3.5 × 10 −2days

     | Show Table
    DownLoad: CSV

    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.

    Figure 2.  The transient behavior of the: (a) pathogen; (b) primary immune cell; (c) specific immune cell; (d) memory cell in the body without intervention of protein drug. Pathogens replicate exponentially and fuse with their carrying capacity while specific immune cells as the pathogen barrier decay by following their oscillatory dynamics towards zero at the fifteen' day.

    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.

    Figure 3.  The transient behavior of the: (a) pathogen; (b) primary immune cell; (c) specific immune cell; (d) memory cell in the host body after the addition of the protein drug. The existence of protein drugs triggers the production of primary cells which directly affect the decreasing of pathogens which exponentially converge to zero in no more than fifteen days.

    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.

    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.

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



    [1] D. D. Chaplin, Overview of the immune response, J. Allergy. Clin. Immun., 125 (2010), S3–S23. https://doi.org/10.1016/j.jaci.2009.12.980 doi: 10.1016/j.jaci.2009.12.980
    [2] C. A. Jr. Janeway, P. Travers, M. Walport, et al., Immunobiology: The Immune System in Health and Disease. 5th edition. New York: Garland Science, 2001.
    [3] C. R. Maldini, G. I. Ellis, J. L. Riley, CAR T cells for infection, autoimmunity and allotransplantation, Nat. Rev. Immunol., 18 (2018), 605–616. https://doi.org/10.1038/s41577-018-0042-2 doi: 10.1038/s41577-018-0042-2
    [4] L. B. Nicholson, The immune system, Essays Biochem., 60 (2016), 275–301. https://doi.org/10.1042/EBC20160017 doi: 10.1042/EBC20160017
    [5] J. M. Carton, W. R. Strohl, Protein therapeutics (introduction to biopharmaceuticals): Introduction to Biological and Small Molecule Drug Research and Development, Elsevier, 2013,127–159, ISBN 9780123971760, https://doi.org/10.1016/B978-0-12-397176-0.00004-2.10.1016/B978-0-12-397176-0.00004-2
    [6] M. Lever, T. D. C. Hirata, P. Russo, H. I. Nakaya, Systems immunology, Theoretical and Applied Aspects of Systems Biology., Springer International Publishing (2018), 159–173. https://doi.org/10.1007/978-3-319-74974-7.
    [7] N. Chirmule, V. Jawa, B. Mibohm, Immunogenicity to therapeutic proteins: impact on PK/PD and efficacy, AAPS J., 14 (2012), 296–302. https://doi.org/10.1208/s12248-012-9340-y doi: 10.1208/s12248-012-9340-y
    [8] K. Bloem, B. Hernández-Breijo, A. Martínez-Feito, T. Rispens, Immunogenicity of therapeutic antibodies: Monitoring antidrug antibodies in a clinical context, Ther Drug Monit., 4 (2017), 327–332. https://doi.org/10.1097/FTD.0000000000000404
    [9] J. J. P. Ruixo, P. Ma, A. T. Chow, The utility of modeling and simulation approaches to evaluate immunogenicity effect on the therapeutic protein pharmacokinetics, AAPS J., 15 (2013), 172–182. https://doi.org/10.1208/s12248-012-9424-8 doi: 10.1208/s12248-012-9424-8
    [10] R. Eftimie, J. J. Gillard, D. A. Cantrell, Mathematical models for immunology: Current state of the art and future research directions, B. Math. Biol., 78 (2016), 2091–2134. https://doi.org/10.1007/s11538-016-0214-9 doi: 10.1007/s11538-016-0214-9
    [11] S. Banerjee, S. Khajanchi, S. Chaudhuri, A mathematical model to elucidate brain tumor abrogation by immunotherapy with T11 target structure, PLoS ONE, 10 (2015), e0123611, https://doi.org/10.1371/journal.pone.0123611. doi: 10.1371/journal.pone.0123611
    [12] J. Reyes-Silveyra, A. R. Mikler, Modeling immune response and its effect on infectious disease outbreak dynamic, Theor. Biol. Med. Model., 13 (2016). https://doi.org/10.1186/s12976-016-0033-6 doi: 10.1186/s12976-016-0033-6
    [13] S. Kathman, T. M. Thway, L. Zhou, S. Lee, Utility of a bayesian mathematical model to predict the impact of immunogenicity on pharmacokinetics of therapeutic proteins, AAPS J., 18 (2016), 424–431. https://doi.org/10.1208/s12248-015-9853-2 doi: 10.1208/s12248-015-9853-2
    [14] W. L. Duan, H. Fang, C. Zeng, The stability analysis of tumor-immune responses to chemotherapy system with gaussian white noises. Chaos, Soliton. Fract., 127 (2019), 96–102. https://doi.org/10.1016/j.chaos.2019.06.030. doi: 10.1016/j.chaos.2019.06.030
    [15] G. A. Bocharov, V. Volpert, B. Ludewig, A. Meyerhans, Mathematical modeling of the immune system in homeostasis, infection and disease, Front. Immunol., 10 (2020), 2944. https://doi.org/10.3389/fimmu.2019.02944 doi: 10.3389/fimmu.2019.02944
    [16] G. A. Bocharov, D. S. Grebennikov, R. S. Savinkov, Mathematical immunology: From phenomenological to multiphysics modelling, Russ. J. Numer. Anal. Mathematical M., 35 (2020), 203–213. https://doi.org/10.1515/rnam-2020-0017 doi: 10.1515/rnam-2020-0017
    [17] S. A. Alharbi, A. Z. Rambely, Dynamic behaviour and stabilisation to boost the immune system by complex interaction between tumour cells and vitamins intervention, Adv. Differ. Equ-Ny., 1 (2020), 412. https://doi.org/10.1186/s13662-020-02869-6 doi: 10.1186/s13662-020-02869-6
    [18] A. Fenton, J. Lello, M. B. Bonsall, Pathogen responses to host immunity: The impact of time delays and memory on the evolution of virulence, Proc. R. Soc. B., 273 (2006). https://doi.org/10.1098/rspb.2006.3552 doi: 10.1098/rspb.2006.3552
    [19] F. A. Rihan, D. H. A. Rahman, Delay differential model for tumour-immune dynamics with HIV infection of CD4+ T-cells, Int. J. Comput. Math., 90 (2013), 594–614, http://dx.doi.org/10.1080/00207160.2012.726354. doi: 10.1080/00207160.2012.726354
    [20] F. A. Rihan, D. H. Abdelrahman, F. Al-Maskari, F. Ibrahim, M. A. Abdeen, Delay differential model for tumour-immune response with chemoimmunotherapy and optimal control, Comput. Math. Method. M., 2014 (2014), Article ID 982978. http://dx.doi.org/10.1155/2014/982978.10.1155/2014/982978.
    [21] S. Khajanchi, S. Banerjee, Stability and bifurcation analysis of delay induced tumor immune interaction model, Appl. Math. Comput., 248 (2014), 652–671. https://doi.org/10.1016/j.amc.2014.10.009. doi: 10.1016/j.amc.2014.10.009
    [22] S. Kayan, H. Merdan, R. Yafia, S. Goktepe, Bifurcation analysis of a modified tumor-immune system interaction model involving time delay, Math. Model. Nat. Pheno., 12 (2017), 120–145. https://doi.org/10.1051/mmnp/201712508 doi: 10.1051/mmnp/201712508
    [23] F. A. Rihan, S. Lakshmanan, H. Maurer, Optimal control of tumour-immune model with time-delay and immuno-chemotherapy. Appl. Math. Comput., 353 (2019), 147–165. https://doi.org/10.1016/j.amc.2019.02.002. doi: 10.1016/j.amc.2019.02.002
    [24] J. P. Mendonça, I. Gleria, M. L. Lyra, Delay-induced bifurcations and chaos in a two-dimensional model for the immune response, Physica A, 517 (2019), 484–490. https://doi.org/10.1016/j.physa.2018.11.039 doi: 10.1016/j.physa.2018.11.039
    [25] F. Fatehi, Y. N. Kyrychko, K. B. Blyuss, Time-delayed model of autoimmune dynamics, Math. Biosci. Eng., 16 (2019), 5613–5639. https://doi.org/10.3934/mbe.2019279 doi: 10.3934/mbe.2019279
    [26] Kasbawati, Mariani, N. Erawaty, N. Aris, A mathematical study of effects of delays arising from the interaction of anti-drug antibody and therapeutic protein in the immune response system, AIMS Math., 5 (2020), 7191–7213. https://doi.org/10.3934/math.2020460 doi: 10.3934/math.2020460
    [27] P. Das, P. Das, S. Das, Effects of delayed immune-activation in the dynamics of tumor-immune interactions, Math. Model. Nat. Pheno., 15 (2020), 45. https://doi.org/10.1051/mmnp/2020001 doi: 10.1051/mmnp/2020001
    [28] Q. Tang, G. Zhang, Stability and Hopf bifurcations in a competitive tumour-immune system with intrinsic recruitment delay and chemotherapy, Math. Biosci. Eng., 18 (2021), 1941–1965. https://doi.org/10.3934/mbe.2021101 doi: 10.3934/mbe.2021101
    [29] W. L. Duan, L. Lin, Noise and delay enhanced stability in tumor-immune responses to chemotherapy system, Chaos, Soliton. Fract., 148 (2021), 111019. https://doi.org/10.1016/j.chaos.2021.111019. doi: 10.1016/j.chaos.2021.111019
  • This article has been cited by:

    1. Muhammad Rifqy Adha Nurdiansyah, Syamsuddin Toaha, Stability analysis and numerical simulation of rabies spread model with delay effects, 2024, 9, 2473-6988, 3399, 10.3934/math.2024167
    2. Xiong Zhang, Zhongyi Xiang, Piecewise immunosuppressive infection model with viral logistic growth and effector cell-guided therapy, 2024, 9, 2473-6988, 11596, 10.3934/math.2024569
  • Reader Comments
  • © 2022 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(1666) PDF downloads(82) Cited by(2)

Figures and Tables

Figures(3)  /  Tables(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog