# Mathematical Modelling of the Inhibitory Role of Regulatory T Cells in Tumor Immune Response

Authors: Zhongtao Yang; Cuihong Yang; Yueping Dong; Yasuhiro Takeuchi
Journal: Complexity
 (2020)

---

## Abstract

The immune system against tumors acts through a complex dynamical process showing a dual role. On the one hand, the immune system can activate some immune cells to kill tumor cells (TCs), such as cytotoxic T lymphocytes (CTLs) and natural killer cells (NKs), but on the other hand, more evidence shows that some immune cells can help tumor escape, such as regulatory T cells (Tregs). In this paper, we propose a tumor immune interaction model based on Tregs-mediated tumor immune escape mechanism. When helper T cells’ (HTCs) stimulation rate by the presence of identified tumor antigens is below critical value, the coexistence (tumor and immune) equilibrium is always stable in its existence region. When HTCs stimulation rate is higher than the critical value, the inhibition rate of effector cells (ECs) by Tregs can destabilize the coexistence equilibrium and cause Hopf bifurcations and produce a limit cycle. This model shows that Tregs might play a crucial role in triggering the tumor immune escape. Furthermore, we introduce the adoptive cellular immunotherapy (ACI) and monoclonal antibody immunotherapy (MAI) as the treatment to boost the immune system to fight against tumors. The numerical results show that ACI can control TCs more, while MAI can delay the inhibitory effect of Tregs on ECs. The result also shows that the combination of both immunotherapies can control TCs and reduce the inhibitory effect of Tregs better than a single immunotherapy can control.


---

## Body

## 1. Introduction

Tumors can be benignant (not cancerous), premalignant (precancerous), and malignant (cancerous). Every year millions of people suffer with cancer and die from this disease throughout the world [1]. It is important to understand tumor’s mechanisms of establishment and destruction, cell-mediate immunity with cytotoxic T lymphocytes (CTLs), and natural killer cells (NKs), generally called effector cells (ECs) that are cytotoxic to tumor cells (TCs), and play a basic role in immune response against tumors [2, 3]. Moreover, efficient antitumor immunity requires the action of helper T cells (HTCs), which can directly activate naive CD8+ T cells to differentiate into CTLs [4–6]. Recently, it has been reported that regulatory T cells (Tregs) can inhibit CTLs and promote the escape of TCs [7]. Tregs suppress immune cells, and when the war between T cells and infection is over, the Tregs signal to stop [8]. Cancer immunotherapy fights against cancer by strengthening the body’s immune system, but the involvement of Tregs inhibits the immune response and turns off the anticancer effect. Tregs inhibition is important in the dynamics of the tumor immune system, which is one motivation of this work.Adoptive T cell immunotherapy (ACI) as a common immunotherapy involves injecting adoptive T cells directly into tumor patients [9–11]. Its advantages are good destruction of tumor and persistence, while its disadvantages are serious toxic and side effects. The monoclonal antibody immunotherapy (MAI) is the immune checkpoint inhibitor [12, 13], which has the advantage of removing the suppression state of the immune system and restoring the immune function of the body to TCs and the disadvantage of having serious immune-related adverse reactions.Tregs have become an important target in tumor immunotherapy because of their contribution to tumor immune escape. Cytotoxic T lymphocyte antigen 4 (CTLA-4) is a marker that is expressed on the surface of activated T cells and transmits inhibitory signals in the immune response [14–16]. Blocking CTLA-4 can reduce the inhibitory activity of Tregs, and the anti-CTLA-4 humanized monoclonal antibody Ipilimumab and Tremelimumab are used to treat advanced melanoma and malignant mesothelioma, respectively [17]. Similar to CTLA-4, programmed death receptor 1 (PD-1) can also promote the activation and development of Tregs [18–21]. Blocking PD-1 can prevent the development of Tregs and prevent the conversion of HTCs into Tregs [22]. Currently, OPDIVO (Nivolumab), an anti-PD-1 monoclonal antibody, has been approved by the US FDA for the treatment of melanoma, renal cell carcinoma, and non-small cell lung cancer [23–25]. Establishment of a mathematical model to study the immunotherapy on the reduction of Tregs inhibition has both theoretical and practical significance.In order to describe the mechanisms of host’s own immune response to against TCs, various types of mathematical models have been proposed [26–43]. The modelling of the tumor immune system described by ordinary differential equations (ODEs) has a long history, which can be traced back to the classic research of Stepanova in 1980 [26]. In 1994, Kuznetsov et al. established the famous two-dimensional ODEs model, postulating that tumor growth follows the Logistic growth pattern. They evaluated the parameters of the model by fitting experimental data from mice [27]. In 2003, Stolongo-Costa et al. assumed that TCs follows the exponential growth pattern and constructed a two-dimensional ODEs model. They analyzed the basic properties of the model and provided conditions for stability of the tumor-free equilibrium, explaining its epidemiological significance [28]. In 2004, Galach simplified Kuznetsov’s system to account for the effect of immune delay on the tumor immune system [29]. In 2014, Dong et al. constructed a three-dimensional ODEs model focusing on the effects of HTCs on the tumor immune system [4].In 1998, Kirschner and Panetta generalized Kuznetsov–Taylor model and illustrated the dynamics between TCs, ECs, and IL-2. They firstly introduced ACI into their model which can explain both short-term tumor oscillations in tumor sizes as well as long-term tumor relapse [11]. In 2003, in order to study the role of cytokine therapy in the activation of the immune system, Stolongo-Costa et al. introduced cycle therapy term Fcos2ωt and established a cycle immunotherapy model. They obtained some thresholds of the frequency and intensity of immunotherapy [28]. In 2006, de Pillis et al. constructed the six-dimensional ODEs model to investigate the effects of combined chemotherapy and immunotherapy on tumor control. They briefly analyzed the nature of the model and discussed the optimal treatment using optimal control theory [30]. In 2008, Bunimovich-Mendrazitsky et al. established a pulsed differential equation model with Bacillus Calmette-Guerin tumor immunotherapy. They obtained the critical threshold and pulse frequency of BCG injection dose that could successfully treat superficial bladder cancer [31]. In 2012, Wilson and Levy established a mathematical model containing Tregs. They studied the absence of treatment, vaccine treatment, anti-TGF treatment, and combination vaccine and anti-TGF treatment, as well as sensitivity analysis of some important parameters [8]. In 2018, Radunskaya et al. established a mathematical model with blood, spleen, and tumor compartments to study PD-L1 inhibitors in the role of tumor immunotherapy. The model was used to fit parameters with the experimental data. The results showed that increasing the resistance of PD-L1 doses can greatly improve the clearance rate of tumor [32].This paper investigates the role of Tregs in the tumor immune system. Therefore, we incorporate the fourth population of Tregs into the previous system in [8]. For the mathematical simplicity, a bilinear term also has been used to describe the interactions between immune response and tumor. To our knowledge, HTCs can recognize TCs and promote the growth of ECs. And ECs can provide direct protective immunity by attacking TCs. When there are more HTCs and ECs, in order to maintain immune homeostasis, the body will produce corresponding Tregs to suppress ECs, and Tregs originating from both HTCs and ECs. Then, we establish a four-dimensional ODEs model described as below:(1)dTtdt=aTt1−bTt−nEtTt,dEtdt=pEtHt−qRtEt−d1Et,dHtdt=s2+k2TtHt−d2Ht,dRtdt=r1Et+r2Ht−d3Rt,where Tt, Et, Ht, and Rt represent the populations of TCs, ECs, HTCs, and Tregs, respectively. The first equation describes the rate change of TCs population. The tumor follows logistic growth dynamics with growth rate a, and the maximum capacity is 1/b. n represents the loss rate of TCs by ECs interaction. The second equation describes the rate change of the ECs population. d1 is the mortality rate of ECs. p is the activation rate of ECs by HTCs, and q is the inhibition rate of Tregs on ECs. The third equation describes the rate change of the HTCs population, s2 is birth rate of HTCs produced in the bone marrow, and HTCs have a natural lifespan of an average 1/d2 days. k2 is HTCs stimulation rate by the presence of identified tumor antigens. The fourth equation gives the rate change of the Tregs population, r1 and r2 are the activation rates of Tregs by ECs and HTCs, respectively. d3 represents per capita decay rate of Tregs. A diagram of the various interactions between these cell populations is shown in Figure 1.Figure 1
A diagram of interactions among the different cell populations in model (1).We nondimensionalize model (1) by taking the following scaling:(2)t=τnT0,Tt=T0xτ,Et=E0yτ,Ht=H0zτ,Rt=R0uτ,α=anT0,β=bT0,ρ=pn,ω2=k2n,θ=qR0nT0,δ1=d1nT0,δ2=d2nT0,δ3=d3nT0,σ2=s2nT0H0,γ1=r1nR0,γ2=r2nR0,and we choose the scaling T0=E0=H0=R0=106. By replacing τ by t, we obtain the following scaled model:(3)dxtdt=αxt1−βxt−xtyt,dytdt=ρytzt−θytut−δ1yt,dztdt=σ2+ω2xtzt−δ2zt,dutdt=γ1yt+γ2zt−δ3ut,with initial conditions(4)x0=x0≥0,y0=y0≥0,z0=z0≥0,u0=u0≥0.Here,x,y,z, and u denote the dimensionless densities of TCs, ECs, HTCs, and Tregs populations, respectively.


## 2. Model Analysis

### 2.1. Well Posedness of Model (3)

The following proposition establishes the well posedness of model (3) with initial conditions (4).Proposition 1.
The solutionsxt,yt,zt,ut of model (3) with initial conditions (4) are existent, unique, and nonnegative on the interval 0,+∞.Proof.
Since the right-hand side of model (3) is completely continuous and locally Lipschitz on the interval 0,+∞, there exists a constant δ>0 such that the solutions xt,yt,zt,ut of model (3) with initial conditions (4) are existent and unique on the interval 0,δ, where 0<δ≤+∞ [44].
From model (3) with initial conditions (4), we obtain(5)xt=x0e∫0tα1−βs−ysds≥0,yt=y0e∫0tρzs−θus−δ1ds≥0,zt=z0e∫0tω2xτ−δ2dτ+σ2∫0te∫stω2xτ−δ2dτds≥0,ut=u0e−δ3t+∫0tγ1ys+γ2zse−δ3t−sds≥0.
It is obvious to see that the solutionsxt,yt,zt,ut of model (3) with initial conditions (4) are nonnegative for all t≥0. The proof is complete.


### 2.2. Existence of Equilibria

In this section, we will study the existence of various equilibria for system (3). We set dx/dt=0,dy/dt=0,dz/dt=0,anddu/dt=0 in system (3), and we have(6)αx1−βx−xy=0,ρyz−θyu−δ1y=0,σ2+ω2xz−δ2z=0,γ1y+γ2z−δ3u=0.Puttingx=0andy=0 yields the tumor-free and ECs-free equilibrium, namely,(7)P0=x0,y0,z0,u0=0,0,σ2δ2,γ2σ2δ2δ3,which always exists.Puttingx=0andy≠0 yields the tumor-free equilibrium, namely,(8)P1=x1,y1,z1,u1=0,δ3ρσ2−δ1δ2−γ2θσ2θγ1δ2,σ2δ2,ρσ2−δ1δ2δ2θ,which exists when ρ>δ1δ2/σ2 and θ<δ3ρσ2−δ1δ2/γ2σ2≜θ1.Puttingx=1/β yields the tumor-dominant equilibrium, namely,(9)P2=x2,y2,z2,u2=1β,0,βσ2βδ2−ω2,βγ2σ2δ3βδ2−ω2,which exists when ω2<βδ2≜ω2∗.Puttingx≠0andx≠1/β and eliminating y,z,andu in (6), we have(10)Fx=Ax2+Bx+C=0,where A=−αβθγ1ω2<0,B=θαγ1ω2∗+ω2+δ1δ3ω2>0,andC=δ3ρσ2−δ1δ2−θαγ1δ2+γ2σ2.Then, we have the coexistence equilibrium(11)P∗=x∗,y∗,z∗,u∗=x∗,α1−βx∗,σ2δ2−ω2x∗,αγ11−βx∗δ2−ω2x∗+γ2σ2δ3δ2−ω2x∗,which exists when x∗<μ≜min1/β,δ2/ω2, where x∗ is the positive root of (10).Below, we consider the existence condition ofP∗:Case (1):C≥0, i.e., θ≤δ3ρσ2−δ1δ2/αγ1δ2+γ2σ2≜θ0.SinceF1/β=θαβγ1δ2+δ1δ3ω2+Cβ/β>0 and Fδ2/ω2=θαγ1δ2+δ1δ2δ3+C>0, we have 1/β<x∗andδ2/ω2<x∗, which contradicts with x∗<μ. Hence, P∗ does not exist in this case.Case (2):C<0, i.e., θ>θ0.LetΔ=B2−4AC=4αβγ1θω2δ3ρσ2−δ1δ2−θαγ1δ2+γ2σ2+θαγ1ω2∗+ω2+δ1δ3ω22=θαγ1ω2−ω2∗+δ1δ3ω22+4αβθγ1ω2σ2ρδ3−θγ2.IfΔ<0, then (10) does not have a positive root. Hence, P∗ does not exist in this case.IfΔ=0, then (10) has one positive root x∗, where x∗=θαγ1ω2∗+ω2+δ1δ3ω2/2αβθγ1ω2>ω2∗+ω2/2βω2=1/2δ2/ω2+1/β≥μ, which contradicts with x∗<μ. Hence, P∗ does not exist in this case.IfΔ>0, then (10) has two positive roots, where x+∗=θαγ1ω2∗+ω2+δ1δ3ω2+Δ/2αβθγ1ω2 and x−∗=θαγ1ω2∗+ω2+δ1δ3ω2−Δ/2αβθγ1ω2. Since x+∗>θαγ1ω2∗+ω2+δ1δ3ω2/2αβθγ1ω2>ω2∗+ω2/2βω2=1/2δ2/ω2+1/β≥μ, which contradicts with x∗<μ. If x−∗<μ, i.e., F1/β=δ1δ3ω2+βδ2ρσ2−δ1δ2−θβγ2σ2/β>0⇒θ<δ1δ3ω2+βδ3ρσ2−δ1δ2/βγ2σ2≜θ2 and Fδ2/ω2=σ2ρδ3−θγ2>0⇒θ<ρδ3/γ2≜θ3, there exists a unique coexistence equilibrium P∗=x∗,y∗,z∗,u∗=x∗,α1−βx∗,σ2/δ2−ω2x∗,αγ11−βx∗δ2−ω2x∗+γ2σ2/δ3δ2−ω2x∗, where x∗=x−∗.The existence conditions for each equilibrium are given in Table1.Table 1
The existence and stability conditions of each equilibrium.
Existence conditionsStability conditionsP0AlwaysUnstableP1θ<θ1θ<θ0P2ω2<ω2∗θ>θ2P∗θ>θ0i:ω2≤ω2∗,θ<θ2<θ3ii:ω2>ω2∗,θ<θ3<θ2(20)


### 2.3. Stability of Equilibria

In order to investigate the local stability of the above equilibriaP0,P1,P2,andP∗ of system (3), we linearize the system and obtain Jacobian matrix at each equilibrium P¯x¯,y¯,z¯,u¯:(12)JP¯=α−2αβx¯−y¯−x¯000ρz¯−θu¯−δ1ρy¯−θy¯ω2z¯0ω2x¯−δ200γ1γ2−δ3.The corresponding characteristic equation is(13)detJP¯−λI=α−2αβx¯−y¯−λ−x¯000ρz¯−θu¯−δ1−λρy¯−θy¯ω2z¯0ω2x¯−δ2−λ00γ1γ2−δ3−λ=0.Theorem 1.
System (3) always has one tumor-free and ECs-free equilibrium P0, which is unstable.Proof.
AtP0, characteristic (13) becomes(14)λ−αλ+δ2λ+δ3λ+θγ2σ2+δ1δ2δ3−ρσ2δ3δ2δ3=0.
It can be seen that one eigenvalueα is positive. Hence, P0 is unstable.Theorem 2.
System (3) has one tumor-free equilibrium P1, when ρ>δ1δ2/σ2 and θ<θ1, which is locally asymptotically stable (LAS) if the inequality θ<θ0 holds.Proof.
AtP1, characteristic (13) becomes(15)λ−θαγ1δ2+γ2σ2−δ3ρσ2−δ1δ2θγ1δ2λ+δ2λ2+δ3λ+θγ1y1=0.
Then, one root of characteristic equation is−δ2<0. It is easily noted that as δ3>0,θγ1y1>0, so λ2+δ3λ+θγ1y1=0 has solutions with negative real parts. Therefore, P1 is LAS when θαγ1δ2+γ2σ2−δ3ρσ2−δ1δ2/θγ1δ2<0, that is, θ<θ0.Theorem 3.
System (3) has one tumor-dominant equilibrium P2, when ω2<ω2∗, which is LAS if the inequality θ>βδ3ρσ2−δ1δ2+δ1δ3ω2/βγ2σ2≡θ2 holds.Proof.
AtP2, characteristic (13) becomes(16)λ+αλ+δ3λ−ω2−ω2∗βλ−ρz2−θu2−δ1=0.
The three roots of characteristic equation areλ1=−α<0,λ2=−δ3<0, and λ3=ω2−ω2∗/β<0, if ω2<ω2∗. So, P2 is LAS when λ4=ρz2−θu2−δ1=ρβσ2δ3−δ1δ3ω2∗−ω2−θβγ2σ2/δ3ω2∗−ω2<0, which is equivalent to θ>θ2.
AtP∗, characteristic (13) becomes(17)λ4+A1λ3+A2λ2+A3λ+A4=0,where(18)A1=αβx∗+δ3+δ2−ω2x∗>0,A2=θγ1y∗+αβδ2−ω2x∗x∗+αβx∗+δ2−ω2x∗δ3>0,A3=θγ1αβx∗+δ2−ω2x∗y∗+αβδ3δ2−ω2x∗x∗+ρω2x∗y∗z∗>0,A4=θαβγ1δ2−ω2x∗x∗y∗+ρδ3−θγ2ω2x∗y∗z∗.
Note thatδ2>ω2x∗ and θ<minθ2,θ3 are necessary for the existence of P∗, then we have A4>0. By the Routh–Hurwitz criterion, the roots of (17) have only negative real parts if and only if(19)A1>0,A2>0,A3>0,A4>0,A1A2−A3>0,A1A2A3−A1A4−A32>0.
Hence, we obtain the sufficient conditions for stability ofP∗:(20)H1:δ2>ω2x∗,θ<minθ2,θ3,H2:A1A2−A3=θγ1αβx∗+δ2−ω2x∗y∗+αβδ3δ2−ω2x∗x∗−ρx∗y∗z∗>0,H3:A1A2A3−A1A4−A32=B1θ+B2>0,where(21)B1=αβx∗+δ2−ω2x∗2δ3+αβx∗+δ2−ω2x∗δ32−ρω2x∗y∗z∗αβx∗+δ2−ω2x∗γ1y∗−αβx∗+δ2−ω2x∗+δ32αβγ1δ2−ω2x∗x∗y∗−ω2γ2x∗y∗z∗,B2=αβx∗+δ2−ω2x∗2δ3+αβx∗+δ2−ω2x∗δ32−ρω2x∗y∗z∗αβδ3δ2−ω2x∗x∗+ρω2x∗y∗z∗−αβx∗+δ2−ω2x∗+δ32ρω2δ3x∗y∗z∗.Theorem 4.
System (3) has a unique coexistence equilibrium P∗, when ω2≤ω2∗,θ0<θ<θ2 or ω2>ω2∗,θ0<θ<θ3 hold, and P∗ is LAS, when conditions (20) are satisfied.
We summarize the above results in Table1, where ω2∗=βδ2,θ0=δ3ρσ2−δ1δ2/αγ1δ2+γ2σ2, θ1=δ3ρσ2−δ1δ2/γ2σ2, θ2=βδ3ρσ2−δ1δ2+δ1δ3ω2/βγ2σ2, and θ3=ρδ3/γ2.


### 2.4. Numerical Simulations

In this section, we choose some suitable parameters in (3) to simulate numerically the theoretical conclusions obtained in the previous sections by using the Matlab software package MATCONT [45]. We select the following parameters’ set [4]:(22)α=1.636,β=0.002,δ1=0.3743,σ2=0.38,δ2=0.055,ρ=0.48,γ1=0.15,γ2=0.2,δ3=0.25.Note thatω2∗=βδ2=0.00011. By calculations, we have(23)θ0=δ3ρσ2−δ1δ2αγ1δ2+γ2σ2=0.452,θ1=δ3ρσ2−δ1δ2γ2σ2=0.5323,θ3=ρδ3γ2=0.6,θ2=βδ3ρσ2−δ1δ3+δ1δ3ω2βγ2σ2=615.625ω2+0.5323.And we find the stability region ofP∗=x∗,y∗,z∗,u∗ (see Figure 2). P∗ is LAS in regions I and II. P∗ is unstable in region III.Figure 2
Stability region ofP∗ and Hopf bifurcation curve in θ−ω2 parameter plane.Let us denote the point inθ−ω2 plane as Qi=θ,ω2.Case (a): we choose a pointQ1=0.48,0.0001 in the region I; then, system (3) has one interior equilibrium:(24)P1∗=77.6449,1.38195,8.0448,7.265.The eigenvalues of Jacobian matrix of (12) are(25)−0.4623,−0.049151,−0.01991−0.378iand−0.01991+0.378i,soP1∗ is stable, as shown in Figure 3(a).Case (b): we choose a pointQ2=0.455,0.0004 in the region II; then, system (3) has one interior equilibrium:(26)P2∗=2.84485,1.62669,7.05506,6.62006.The eigenvalues of Jacobian matrix of (12) are(27)−0.0931+0.301i,−0.0931−0.301i,−0.0633+0.0188iand−0.0633−0.0188i,soP2∗ is stable, as shown in Figure 3(b).Case (c): we choose a pointQ3=0.48,0.000111 in the region II, then system (3) has one interior equilibrium:(28)P3∗=72.569,1.39855,8.094,7.314.The eigenvalues of Jacobian matrix of (12) are(29)−0.458,−0.0491,−0.0132+0.383iand−0.0132−0.383i,soP3∗ is stable, as shown in Figure 3(c).Case (d): we choose a pointQ4=0.48,0.0004 in the region III; then, system (3) has one interior equilibrium:(30)P4∗=26.365,1.54973,8.54816,7.76837.Figure 3
The other parameter values are given in (22). (a) P1∗77.644,1.381,8.044,7.265 is LAS, θ=0.48andω2=0.0001. (b) P2∗2.844,1.626,7.055,6.62 is LAS, θ=0.455andω2=0.0004. (c) P3∗72.569,1.398,8.094,7.314 is LAS, θ=0.48andω2=0.000111. (d) P4∗26.365,1.549,8.548,7.768 is unstable, θ=0.48andω2=0.0004. (e) The 3D phase portrait depicts tumor cell, ECs, and Tregs, θ=0.48andω2=0.0004.
(a)(b)(c)(d)(e)The eigenvalues of Jacobian matrix of (12) are(31)−0.4328,−0.0497,0.0509+0.416iand0.0509−0.416i,so P4∗ is unstable, as shown in Figure 3(d).Below, we perform numerically bifurcation analysis ofx against θ for different values of ω2.Case (1): we chooseω2=0.0001<ω2∗=0.00011 (see Figure 4(a)).We obtain the following result.Figure 4
The bifurcation diagrams ofx with respect to θ for different ω2=0.0001 and 0.0004, respectively. The stable state is represented by the blue curve, and the unstable one corresponds to the red curve. The equilibria are P00,0,z0,u0,P10,y1,z1,u1,P2500,0,z2,u2, and P∗x∗,y∗,z∗,u∗ of system (3) (in order to distinguish P0 from P1, we move P1 up a little bit). (a) ω2=0.0001. (b) ω2=0.0004.
(a)(b)Proposition 2.
When0<θ<θ0, P0 exists and is unstable, P1 exists and is LAS, P2 exists and is unstable, and P∗ does not exist. When θ0<θ<θ1, P0 exists and is unstable, P1 exists and is unstable, P2 exists and is unstable, and P∗ exists and is LAS. When θ1<θ<θ2, P0 exists and is unstable, P1 does not exist, P2 exists and is unstable, and P∗ exists and is LAS. When θ>θ2, P0 exists and is unstable, P1 does not exist, P2 exists and is LAS, and P∗ does not exist.Case (2): we chooseω2=0.0004>ω2∗=0.00011 (see Figure 4(b)).
We obtain the following result.Proposition 3.
When0<θ<θ0, P0 exists and is unstable, P1 exists and is LAS, and P∗ does not exist. When θ0<θ<θ4, where θ4=0.4657 is a Hopf bifurcation, P0 exists and is unstable, P1 exists and is LAS, and P∗ exists and is LAS. When θ4<θ<θ1, P0 exists and is unstable, P1 exists and is unstable, and P∗ exists and is unstable. When θ1<θ<θ3, P0 exists and is unstable, P1 does not exist, and P∗ exists and is unstable. When θ>θ3, P0 exists and is unstable, P1 does not exist, and P∗ does not exist.Consider the case where HTCs stimulation rate (ω2) is low by the presence of identified tumor antigens (ω2<ω2∗). When the inhibition rate of Tregs to ECs θ is lower than θ0, the solution of system (3) approaches P1 implying that ECs can still effectively remove TCs. When θ0<θ<θ2, the solution of system (3) approaches P∗ showing the coexistence of TCs and immune cells, which means that the patient can survive with tumors. When θ>θ2, TCs escape the control of the immune system and develop into malignant tumors.Next consider the case where HTCs stimulation rate (ω2) is high by the presence of identified tumor antigens (ω2>ω2∗). When θ<θ0, the solution of system (3) approaches P1 implying that TCs can be effectively removed by ECs. When θ0<θ<θ4, the solution of system (3) approaches P∗ implying that TCs can still be controlled by the immune system. When θ>θ4, the system has a Hopf bifurcation point and induces a limit cycle (see Figure 3(e)). In the biological sense, it can be understood that the number of TCs presents a periodic change.


## 2.1. Well Posedness of Model (3)

The following proposition establishes the well posedness of model (3) with initial conditions (4).Proposition 1.
The solutionsxt,yt,zt,ut of model (3) with initial conditions (4) are existent, unique, and nonnegative on the interval 0,+∞.Proof.
Since the right-hand side of model (3) is completely continuous and locally Lipschitz on the interval 0,+∞, there exists a constant δ>0 such that the solutions xt,yt,zt,ut of model (3) with initial conditions (4) are existent and unique on the interval 0,δ, where 0<δ≤+∞ [44].
From model (3) with initial conditions (4), we obtain(5)xt=x0e∫0tα1−βs−ysds≥0,yt=y0e∫0tρzs−θus−δ1ds≥0,zt=z0e∫0tω2xτ−δ2dτ+σ2∫0te∫stω2xτ−δ2dτds≥0,ut=u0e−δ3t+∫0tγ1ys+γ2zse−δ3t−sds≥0.
It is obvious to see that the solutionsxt,yt,zt,ut of model (3) with initial conditions (4) are nonnegative for all t≥0. The proof is complete.


## 2.2. Existence of Equilibria

In this section, we will study the existence of various equilibria for system (3). We set dx/dt=0,dy/dt=0,dz/dt=0,anddu/dt=0 in system (3), and we have(6)αx1−βx−xy=0,ρyz−θyu−δ1y=0,σ2+ω2xz−δ2z=0,γ1y+γ2z−δ3u=0.Puttingx=0andy=0 yields the tumor-free and ECs-free equilibrium, namely,(7)P0=x0,y0,z0,u0=0,0,σ2δ2,γ2σ2δ2δ3,which always exists.Puttingx=0andy≠0 yields the tumor-free equilibrium, namely,(8)P1=x1,y1,z1,u1=0,δ3ρσ2−δ1δ2−γ2θσ2θγ1δ2,σ2δ2,ρσ2−δ1δ2δ2θ,which exists when ρ>δ1δ2/σ2 and θ<δ3ρσ2−δ1δ2/γ2σ2≜θ1.Puttingx=1/β yields the tumor-dominant equilibrium, namely,(9)P2=x2,y2,z2,u2=1β,0,βσ2βδ2−ω2,βγ2σ2δ3βδ2−ω2,which exists when ω2<βδ2≜ω2∗.Puttingx≠0andx≠1/β and eliminating y,z,andu in (6), we have(10)Fx=Ax2+Bx+C=0,where A=−αβθγ1ω2<0,B=θαγ1ω2∗+ω2+δ1δ3ω2>0,andC=δ3ρσ2−δ1δ2−θαγ1δ2+γ2σ2.Then, we have the coexistence equilibrium(11)P∗=x∗,y∗,z∗,u∗=x∗,α1−βx∗,σ2δ2−ω2x∗,αγ11−βx∗δ2−ω2x∗+γ2σ2δ3δ2−ω2x∗,which exists when x∗<μ≜min1/β,δ2/ω2, where x∗ is the positive root of (10).Below, we consider the existence condition ofP∗:Case (1):C≥0, i.e., θ≤δ3ρσ2−δ1δ2/αγ1δ2+γ2σ2≜θ0.SinceF1/β=θαβγ1δ2+δ1δ3ω2+Cβ/β>0 and Fδ2/ω2=θαγ1δ2+δ1δ2δ3+C>0, we have 1/β<x∗andδ2/ω2<x∗, which contradicts with x∗<μ. Hence, P∗ does not exist in this case.Case (2):C<0, i.e., θ>θ0.LetΔ=B2−4AC=4αβγ1θω2δ3ρσ2−δ1δ2−θαγ1δ2+γ2σ2+θαγ1ω2∗+ω2+δ1δ3ω22=θαγ1ω2−ω2∗+δ1δ3ω22+4αβθγ1ω2σ2ρδ3−θγ2.IfΔ<0, then (10) does not have a positive root. Hence, P∗ does not exist in this case.IfΔ=0, then (10) has one positive root x∗, where x∗=θαγ1ω2∗+ω2+δ1δ3ω2/2αβθγ1ω2>ω2∗+ω2/2βω2=1/2δ2/ω2+1/β≥μ, which contradicts with x∗<μ. Hence, P∗ does not exist in this case.IfΔ>0, then (10) has two positive roots, where x+∗=θαγ1ω2∗+ω2+δ1δ3ω2+Δ/2αβθγ1ω2 and x−∗=θαγ1ω2∗+ω2+δ1δ3ω2−Δ/2αβθγ1ω2. Since x+∗>θαγ1ω2∗+ω2+δ1δ3ω2/2αβθγ1ω2>ω2∗+ω2/2βω2=1/2δ2/ω2+1/β≥μ, which contradicts with x∗<μ. If x−∗<μ, i.e., F1/β=δ1δ3ω2+βδ2ρσ2−δ1δ2−θβγ2σ2/β>0⇒θ<δ1δ3ω2+βδ3ρσ2−δ1δ2/βγ2σ2≜θ2 and Fδ2/ω2=σ2ρδ3−θγ2>0⇒θ<ρδ3/γ2≜θ3, there exists a unique coexistence equilibrium P∗=x∗,y∗,z∗,u∗=x∗,α1−βx∗,σ2/δ2−ω2x∗,αγ11−βx∗δ2−ω2x∗+γ2σ2/δ3δ2−ω2x∗, where x∗=x−∗.The existence conditions for each equilibrium are given in Table1.Table 1
The existence and stability conditions of each equilibrium.
Existence conditionsStability conditionsP0AlwaysUnstableP1θ<θ1θ<θ0P2ω2<ω2∗θ>θ2P∗θ>θ0i:ω2≤ω2∗,θ<θ2<θ3ii:ω2>ω2∗,θ<θ3<θ2(20)


## 2.3. Stability of Equilibria

In order to investigate the local stability of the above equilibriaP0,P1,P2,andP∗ of system (3), we linearize the system and obtain Jacobian matrix at each equilibrium P¯x¯,y¯,z¯,u¯:(12)JP¯=α−2αβx¯−y¯−x¯000ρz¯−θu¯−δ1ρy¯−θy¯ω2z¯0ω2x¯−δ200γ1γ2−δ3.The corresponding characteristic equation is(13)detJP¯−λI=α−2αβx¯−y¯−λ−x¯000ρz¯−θu¯−δ1−λρy¯−θy¯ω2z¯0ω2x¯−δ2−λ00γ1γ2−δ3−λ=0.Theorem 1.
System (3) always has one tumor-free and ECs-free equilibrium P0, which is unstable.Proof.
AtP0, characteristic (13) becomes(14)λ−αλ+δ2λ+δ3λ+θγ2σ2+δ1δ2δ3−ρσ2δ3δ2δ3=0.
It can be seen that one eigenvalueα is positive. Hence, P0 is unstable.Theorem 2.
System (3) has one tumor-free equilibrium P1, when ρ>δ1δ2/σ2 and θ<θ1, which is locally asymptotically stable (LAS) if the inequality θ<θ0 holds.Proof.
AtP1, characteristic (13) becomes(15)λ−θαγ1δ2+γ2σ2−δ3ρσ2−δ1δ2θγ1δ2λ+δ2λ2+δ3λ+θγ1y1=0.
Then, one root of characteristic equation is−δ2<0. It is easily noted that as δ3>0,θγ1y1>0, so λ2+δ3λ+θγ1y1=0 has solutions with negative real parts. Therefore, P1 is LAS when θαγ1δ2+γ2σ2−δ3ρσ2−δ1δ2/θγ1δ2<0, that is, θ<θ0.Theorem 3.
System (3) has one tumor-dominant equilibrium P2, when ω2<ω2∗, which is LAS if the inequality θ>βδ3ρσ2−δ1δ2+δ1δ3ω2/βγ2σ2≡θ2 holds.Proof.
AtP2, characteristic (13) becomes(16)λ+αλ+δ3λ−ω2−ω2∗βλ−ρz2−θu2−δ1=0.
The three roots of characteristic equation areλ1=−α<0,λ2=−δ3<0, and λ3=ω2−ω2∗/β<0, if ω2<ω2∗. So, P2 is LAS when λ4=ρz2−θu2−δ1=ρβσ2δ3−δ1δ3ω2∗−ω2−θβγ2σ2/δ3ω2∗−ω2<0, which is equivalent to θ>θ2.
AtP∗, characteristic (13) becomes(17)λ4+A1λ3+A2λ2+A3λ+A4=0,where(18)A1=αβx∗+δ3+δ2−ω2x∗>0,A2=θγ1y∗+αβδ2−ω2x∗x∗+αβx∗+δ2−ω2x∗δ3>0,A3=θγ1αβx∗+δ2−ω2x∗y∗+αβδ3δ2−ω2x∗x∗+ρω2x∗y∗z∗>0,A4=θαβγ1δ2−ω2x∗x∗y∗+ρδ3−θγ2ω2x∗y∗z∗.
Note thatδ2>ω2x∗ and θ<minθ2,θ3 are necessary for the existence of P∗, then we have A4>0. By the Routh–Hurwitz criterion, the roots of (17) have only negative real parts if and only if(19)A1>0,A2>0,A3>0,A4>0,A1A2−A3>0,A1A2A3−A1A4−A32>0.
Hence, we obtain the sufficient conditions for stability ofP∗:(20)H1:δ2>ω2x∗,θ<minθ2,θ3,H2:A1A2−A3=θγ1αβx∗+δ2−ω2x∗y∗+αβδ3δ2−ω2x∗x∗−ρx∗y∗z∗>0,H3:A1A2A3−A1A4−A32=B1θ+B2>0,where(21)B1=αβx∗+δ2−ω2x∗2δ3+αβx∗+δ2−ω2x∗δ32−ρω2x∗y∗z∗αβx∗+δ2−ω2x∗γ1y∗−αβx∗+δ2−ω2x∗+δ32αβγ1δ2−ω2x∗x∗y∗−ω2γ2x∗y∗z∗,B2=αβx∗+δ2−ω2x∗2δ3+αβx∗+δ2−ω2x∗δ32−ρω2x∗y∗z∗αβδ3δ2−ω2x∗x∗+ρω2x∗y∗z∗−αβx∗+δ2−ω2x∗+δ32ρω2δ3x∗y∗z∗.Theorem 4.
System (3) has a unique coexistence equilibrium P∗, when ω2≤ω2∗,θ0<θ<θ2 or ω2>ω2∗,θ0<θ<θ3 hold, and P∗ is LAS, when conditions (20) are satisfied.
We summarize the above results in Table1, where ω2∗=βδ2,θ0=δ3ρσ2−δ1δ2/αγ1δ2+γ2σ2, θ1=δ3ρσ2−δ1δ2/γ2σ2, θ2=βδ3ρσ2−δ1δ2+δ1δ3ω2/βγ2σ2, and θ3=ρδ3/γ2.


## 2.4. Numerical Simulations

In this section, we choose some suitable parameters in (3) to simulate numerically the theoretical conclusions obtained in the previous sections by using the Matlab software package MATCONT [45]. We select the following parameters’ set [4]:(22)α=1.636,β=0.002,δ1=0.3743,σ2=0.38,δ2=0.055,ρ=0.48,γ1=0.15,γ2=0.2,δ3=0.25.Note thatω2∗=βδ2=0.00011. By calculations, we have(23)θ0=δ3ρσ2−δ1δ2αγ1δ2+γ2σ2=0.452,θ1=δ3ρσ2−δ1δ2γ2σ2=0.5323,θ3=ρδ3γ2=0.6,θ2=βδ3ρσ2−δ1δ3+δ1δ3ω2βγ2σ2=615.625ω2+0.5323.And we find the stability region ofP∗=x∗,y∗,z∗,u∗ (see Figure 2). P∗ is LAS in regions I and II. P∗ is unstable in region III.Figure 2
Stability region ofP∗ and Hopf bifurcation curve in θ−ω2 parameter plane.Let us denote the point inθ−ω2 plane as Qi=θ,ω2.Case (a): we choose a pointQ1=0.48,0.0001 in the region I; then, system (3) has one interior equilibrium:(24)P1∗=77.6449,1.38195,8.0448,7.265.The eigenvalues of Jacobian matrix of (12) are(25)−0.4623,−0.049151,−0.01991−0.378iand−0.01991+0.378i,soP1∗ is stable, as shown in Figure 3(a).Case (b): we choose a pointQ2=0.455,0.0004 in the region II; then, system (3) has one interior equilibrium:(26)P2∗=2.84485,1.62669,7.05506,6.62006.The eigenvalues of Jacobian matrix of (12) are(27)−0.0931+0.301i,−0.0931−0.301i,−0.0633+0.0188iand−0.0633−0.0188i,soP2∗ is stable, as shown in Figure 3(b).Case (c): we choose a pointQ3=0.48,0.000111 in the region II, then system (3) has one interior equilibrium:(28)P3∗=72.569,1.39855,8.094,7.314.The eigenvalues of Jacobian matrix of (12) are(29)−0.458,−0.0491,−0.0132+0.383iand−0.0132−0.383i,soP3∗ is stable, as shown in Figure 3(c).Case (d): we choose a pointQ4=0.48,0.0004 in the region III; then, system (3) has one interior equilibrium:(30)P4∗=26.365,1.54973,8.54816,7.76837.Figure 3
The other parameter values are given in (22). (a) P1∗77.644,1.381,8.044,7.265 is LAS, θ=0.48andω2=0.0001. (b) P2∗2.844,1.626,7.055,6.62 is LAS, θ=0.455andω2=0.0004. (c) P3∗72.569,1.398,8.094,7.314 is LAS, θ=0.48andω2=0.000111. (d) P4∗26.365,1.549,8.548,7.768 is unstable, θ=0.48andω2=0.0004. (e) The 3D phase portrait depicts tumor cell, ECs, and Tregs, θ=0.48andω2=0.0004.
(a)(b)(c)(d)(e)The eigenvalues of Jacobian matrix of (12) are(31)−0.4328,−0.0497,0.0509+0.416iand0.0509−0.416i,so P4∗ is unstable, as shown in Figure 3(d).Below, we perform numerically bifurcation analysis ofx against θ for different values of ω2.Case (1): we chooseω2=0.0001<ω2∗=0.00011 (see Figure 4(a)).We obtain the following result.Figure 4
The bifurcation diagrams ofx with respect to θ for different ω2=0.0001 and 0.0004, respectively. The stable state is represented by the blue curve, and the unstable one corresponds to the red curve. The equilibria are P00,0,z0,u0,P10,y1,z1,u1,P2500,0,z2,u2, and P∗x∗,y∗,z∗,u∗ of system (3) (in order to distinguish P0 from P1, we move P1 up a little bit). (a) ω2=0.0001. (b) ω2=0.0004.
(a)(b)Proposition 2.
When0<θ<θ0, P0 exists and is unstable, P1 exists and is LAS, P2 exists and is unstable, and P∗ does not exist. When θ0<θ<θ1, P0 exists and is unstable, P1 exists and is unstable, P2 exists and is unstable, and P∗ exists and is LAS. When θ1<θ<θ2, P0 exists and is unstable, P1 does not exist, P2 exists and is unstable, and P∗ exists and is LAS. When θ>θ2, P0 exists and is unstable, P1 does not exist, P2 exists and is LAS, and P∗ does not exist.Case (2): we chooseω2=0.0004>ω2∗=0.00011 (see Figure 4(b)).
We obtain the following result.Proposition 3.
When0<θ<θ0, P0 exists and is unstable, P1 exists and is LAS, and P∗ does not exist. When θ0<θ<θ4, where θ4=0.4657 is a Hopf bifurcation, P0 exists and is unstable, P1 exists and is LAS, and P∗ exists and is LAS. When θ4<θ<θ1, P0 exists and is unstable, P1 exists and is unstable, and P∗ exists and is unstable. When θ1<θ<θ3, P0 exists and is unstable, P1 does not exist, and P∗ exists and is unstable. When θ>θ3, P0 exists and is unstable, P1 does not exist, and P∗ does not exist.Consider the case where HTCs stimulation rate (ω2) is low by the presence of identified tumor antigens (ω2<ω2∗). When the inhibition rate of Tregs to ECs θ is lower than θ0, the solution of system (3) approaches P1 implying that ECs can still effectively remove TCs. When θ0<θ<θ2, the solution of system (3) approaches P∗ showing the coexistence of TCs and immune cells, which means that the patient can survive with tumors. When θ>θ2, TCs escape the control of the immune system and develop into malignant tumors.Next consider the case where HTCs stimulation rate (ω2) is high by the presence of identified tumor antigens (ω2>ω2∗). When θ<θ0, the solution of system (3) approaches P1 implying that TCs can be effectively removed by ECs. When θ0<θ<θ4, the solution of system (3) approaches P∗ implying that TCs can still be controlled by the immune system. When θ>θ4, the system has a Hopf bifurcation point and induces a limit cycle (see Figure 3(e)). In the biological sense, it can be understood that the number of TCs presents a periodic change.


## 3. Treatment Model

In order to investigate well the effect of Tregs in tumor immune response under the treatment, we follow the way in [11] to introduce the constant treatment term s1 into the second equation of (1). Since CTLA-4 and PD-1 can inhibit the development of Tregs and prevent the transformation of HTCs into Tregs [12], the effect of CTLA4 or PD-1 on Tregs can be considered. Then, we shall establish a five-dimensional ODEs model described as below:(32)dTtdt=aTt1−bTt−nEtTt,dEtdt=s1+pEtHt−qRtEt−d1Et,dHtdt=s2+k2TtHt−d2Ht,dRtdt=r1Et+r2Ht−d3R−mRtWt,dWtdt=s3−d4Wt,where Wt represents the concentration of monoclonal antibody in human body at time t, s1 represents the treatment term of introducing LAK and TIL into the region of tumor localization, m represents the inhibition rate of monoclonal antibody on Tregs, s3 represents the amount of monoclonal antibody entering human body at time t, and d4 represents the attenuation coefficient of monoclonal antibody.We scale those new parameters in model (32) as follows:(33)Wt=W0vτ,σ1=s1nE0T0,ξ=mW0nT0,σ3=s3nW0T0,δ4=d4nT0,and we choose the scaling W0=106. By replacing τ by t, we obtain the following scaled model with treatment:(34)dxtdt=αxt1−βxt−xtyt,dytdt=σ1+ρytzt−θytut−δ1yt,dztdt=σ2+ω2xtzt−δ2zt,dutdt=γ1yt+γ2zt−δ3ut−ξutvt,dvtdt=σ3−δ4vt,with initial conditions(35)x0=x0≥0,y0=y0≥0,z0=z0≥0,u0=u0≥0,v0=v0≥0.Here,v denotes the dimensionless concentration of monoclonal antibody.


### 3.1. Model Analysis

By a similar proof of Proposition1, we can obtain the well posedness of model (34) with initial conditions (35) as follows.Proposition 4.
The solutionsxt,yt,zt,ut of model (34) with initial conditions (35) are existent, unique, and nonnegative on the interval 0,+∞.Next, we discuss the effect of three types of immunotherapy on tumors such asCase I: ACI model (σ1>0,σ3=0)Case II: MAI model (σ1=0,σ3>0)Case III: combined immunotherapy model (σ1>0,σ3>0)We setdx/dt=0,dy/dt=0,dz/dt=0,du/dt=0,anddv/dt=0 in (34), and we have(36)αx1−βx−xy=0,ρyz−θyu−δ1y+σ1=0,σ2+ω2xz−δ2z=0,γ1y+γ2z−δ3u−ξuv=0,σ3−δ4v=0.


#### 3.1.1. ACI Model

Whenσ1>0andσ3=0, we put v=0 and (36) becomes(37)αx1−βx−xy=0,ρyz−θyu−δ1y+σ1=0,σ2+ω2xz−δ2z=0,γ1y+γ2z−δ3u=0.Puttingx=0 yields the tumor-free equilibrium, namely,(38)E0=x0,y0,z0,u0=0,y0,σ2δ2,γ1δ2y0+γ2σ2δ2δ3,where y0=−B+B2−4AC/2A, A=θγ1δ2>0,B=δ1δ2δ3+θγ2σ2−ρσ2δ3,andC=−σ1δ2δ3<0. E0 always exists.Puttingx≠0 and eliminating y,z,andu in (37), we have(39)Fx=C3x3+C2x2+C1x+C0=0,C3=θα2β2γ1ω2>0,C2=−θα2βγ1ω2−αβθαγ1βδ2+ω2+δ1δ3ω2<0,C1=θα2αβγ1δ2+αγ1ω2+βγ2σ2+αδ3δ1ω2−βρσ2−δ1δ2−σ1ω2δ3,C0=σ1δ2δ3+αδ3ρσ2−δ1δ2−θααγ1δ2+γ2σ2.Case (a): ifF0Fμ=C0Fμ<0 (where μ=min1/β,δ2/ω2, due to F+∞>0 and F−∞<0, then (39) has one positive root x∗ in the interval 0,μ. Therefore, system (34) has a unique coexistence equilibrium:(40)E∗=x∗,y∗,z∗,u∗=x∗,α1−βx∗,σ2δ2−ω2x∗,αγ11−βx∗δ2−ω2x∗+γ2σ2δ3δ2−ω2x∗.We show the existence conditions ofE∗ in the following:Case (1):C0>0,Fμ<0.SinceC0>0, we have θ<σ1δ2δ3+αδ3ρσ2−δ1δ2/ααγ1δ2+γ2σ2≜θ5. Since Fμ<0, we have the following two cases:Case where1/β<δ2/ω2, i.e., ω2<βδ2. Since μ=1/β, we have Fμ=δ3σ1βδ2−ω2/β>0. That is a contradiction.Case where1/β>δ2/ω2, i.e., ω2>βδ2. Since μ=δ2/ω2, we have Fμ=ασ2βδ2−ω2γ2θ−δ3ρ/ω2<0 implying θ>ρδ3/γ2=θ3. If θ5≥θ3, i.e., σ1δ2δ3+αδ3ρσ2−δ1δ2/ααγ1δ2+γ2σ2≥ρδ3/γ2, we have σ1≥ααργ1+δ1γ2/γ2, there exists E∗.Case (2):C0<0,Fμ>0.SinceC0<0, we have θ>θ5. Since Fμ>0, we have the following two cases:Case where1/β<δ2/ω2, i.e., ω2<βδ2. Since μ=1/β, we have Fμ=δ3σ1βδ2−ω2/β>0, there exists E∗.Case where1/β>δ2/ω2, i.e., ω2>βδ2. Since μ=δ2/ω2, we have Fμ=ασ2βδ2−ω2γ2θ−δ3ρ/ω2>0 implying θ<θ3. If θ3≥θ5, i.e., ρδ3/γ2≥σ1δ2δ3+αδ3ρσ2−δ1δ2/ααγ1δ2+γ2σ2, we have σ1≤ααργ1+δ1γ2/γ2, there exists E∗.Case (b): IfF0Fμ=C0Fμ>0, then (39) can have two positive roots x1∗ and x2∗ in the interval 0,μ, where x1∗ and x2∗ are the roots of (39). Therefore, system (34) can have two coexistence equilibria:(41)E1∗=x1∗,y1∗,z1∗,u1∗=x1∗,α1−βx1∗,σ2δ2−ω2x1∗,αγ11−βx1∗δ2−ω2x1∗+γ2σ2δ3δ2−ω2x1∗,(42)E2∗=x2∗,y2∗,z2∗,u2∗=x2∗,α1−βx2∗,σ2δ2−ω2x2∗,αγ11−βx2∗δ2−ω2x2∗+γ2σ2δ3δ2−ω2x2∗.In the following, we show the existence condition ofE1∗ and E2∗.Differentiating (39), we have(43)F′x=3C3x2+2C2x+C1=0,where C3>0 and C2<0. If C1>0 and Δ≜C22−3C1C3>0, then (43) has two positive roots x1+=−C2−Δ/3C3 and x2+=−C2+Δ/3C3. And if C1<0, then (43) has one positive root x2+=−C2+Δ/3C3.Case (1):C0>0,Fμ>0.SinceC0>0, we have θ<θ5. Since Fμ>0, we have the following two cases:Case where1/β<δ2/ω2, i.e., ω2<βδ2. Since μ=1/β, we have Fμ=δ3σ1βδ2−ω2/β>0. If Fx2+<0, then there exist E1∗ and E2∗.Case where1/β>δ2/ω2, i.e., ω2>βδ2. Since μ=δ2/ω2, we have Fμ=ασ2βδ2−ω2γ2θ−δ3ρ/ω2>0 implying θ<θ3. Then, we have θ<minθ3,θ5. If Fx2+<0, then there exist E1∗ and E2∗.Case (2):C0<0,Fμ<0.SinceC0<0, we have θ>θ5. Since Fμ>0, we have the following two cases:Case where1/β<δ2/ω2, i.e., ω2<βδ2. Since μ=1/β, we have Fμ=δ3σ1βδ2−ω2/β>0. That is a contradiction.Case where1/β>δ2/ω2, i.e., ω2>βδ2. Since μ=δ2/ω2, we have Fμ=ασ2βδ2−ω2γ2θ−δ3ρ/ω2<0 implying θ>θ3. Then, we have θ>maxθ3,θ5. If Fx1+>0, then there exist E1∗ and E2∗.Next, we show the stability of the above two equilibriaE0 and E∗ of system (34) with v=0. We linearize the system and obtain the characteristic equation whose expression is the same as (13).Theorem 5.
System (34) with v=0 always has one tumor-free equilibrium E0, which is LAS if the inequality θ<θ5 holds.Proof.
AtE0=x0,y0,z0,u0=0,y0,σ2/δ2,γ1δ2y0+γ2σ2/δ2δ3, the characteristic equation becomes(44)λ−α+y0λ+δ2λ2+σ1y0+δ3λ+θγ1y0+σ1δ3y0=0,and one root of characteristic equation is −δ2<0. It is easily noted that as σ1/y0+δ3>0,θγ1y1+σ1δ3/y0>0, so the solutions of λ2+σ1/y0+δ3λ+θγ1y1+σ1δ3/y0=0 have always negative real parts. Therefore, E0 is LAS if and only if α−y0=2αA+B−B2−4AC/2A<0, that is, if and only if θ<θ5.
AtE∗=x∗,y∗,z∗,u∗, characteristic (13) becomes(45)λ4+B1λ3+B2λ2+B3λ+B4=0,where(46)B1=αβx∗+δ3+σ1y∗+σ2z∗>0,B2=αβx∗+σ2z∗σ1y∗+δ3+αβx∗σ2z∗+δ3σ1y∗+θγ1y∗>0,B3=θγ1y∗+δ3σ1y∗αβx∗+σ2z∗+αβx∗σ2z∗δ3+σ1y∗+ρω2x∗y∗z∗>0,B4=αβx∗σ2z∗θγ1y∗+δ3σ1y∗+ρδ3−θγ2ω2x∗y∗z∗.
By the Routh–Hurwitz criterion, the roots of (45) have negative real parts if and only if B1>0,B2>0,B3>0,B4>0,B1B2−B3>0,andB1B2B3−B1B4−B32>0.
Hence, we obtain the sufficient condition for the stability ofP∗:(47)H4:B4=αβx∗σ2z∗θγ1y∗+δ3σ1y∗+ρδ3−θγ2ω2x∗y∗z∗>0,H5:B1B2−B3=K1θ+K2>0,H6:B1B2B3−B1B4−B32=B3K1θ+K2−B12B4=K3θ2+K4θ+K5>0,where(48)K1=σ1+δ3y∗γ1,K2=δ3+σ1y∗+σ2z∗α2β2x∗2+αβx∗+σ2z∗σ1y∗+δ3+αβx∗σ2z∗+δ3σ1y∗σ1y∗+δ3−ρω2x∗y∗z∗,K3=K1αβx∗+σ2z∗γ1y∗,K4=K2γ1y∗+K1σ1δ3y∗αβx∗+σ2z∗+K1αβx∗σ2z∗δ3+σ1y∗−αβx∗+δ3+σ1y∗+σ2z∗2γ1y∗+ρ+γ2ω2x∗y∗z∗,K5=K2αβx∗+σ2z∗σ1δ3y∗−αβx∗σ2z∗K2δ3+σ1y∗−αβx∗+δ3+σ1y∗+σ2z∗2σ1δ3y∗+K2−δ3ρω2x∗y∗z∗.Theorem 6.
System (34) with v=0 has a coexistence equilibrium E∗ when ω2>βδ2,σ1≥ααργ1+δ1γ2/γ2,θ3<θ<θ5 or ω2>βδ2, σ1≤ααργ1+δ1γ2/γ2,θ5<θ<θ3 or ω2<βδ2,θ>θ5 hold. Further E∗ is LAS when condition (47) is satisfied.


#### 3.1.2. MAI Model

Whenσ1=0,σ3>0, (36) becomes(49)αx1−βx−xy=0,ρyz−θyu−δ1y=0,σ2+ω2xz−δ2z=0,γ1y+γ2z−δ3u−ξuv=0,σ3−δ4v=0.Equation (49) can be simplified as follows:(50)αx1−βx−xy=0,ρyz−θyu−δ1y=0,σ2+ω2xz−δ2z=0,γ1y+γ2z−δ31u=0,where δ31=δ3δ4+ξσ3/δ4.The analysis of (50) is similar to that of (12), and the equilibria of (34) can be obtained as follows:(51)F0=0,0,σ2δ2,γ2σ2δ2δ31,σ3δ4=0,δ31ρσ2−δ1δ2−θγ2σ2θγ1δ2,σ2δ2,ρσ2−δ1δ2θδ2,σ3δ4,F2=1β,0,βσ2βδ2−ω2,βγ2σ2δ31βδ2−ω2,σ3δ4,F∗=x∗,α1−βx∗,σ2δ2−ω2x∗,αγ11−βx∗δ2−ω2x∗+γ2σ2δ31δ2−ω2x∗,σ3δ4.Here,x∗∈0,μ and satisfies the equation a2x2+a1x+a0=0, where coefficients are defined as a2=−θαβγ1ω2<0,a1=θαγ1βδ2+ω2+δ1δ31ω2>0,anda0=δ31ρσ2−δ1δ2−θαγ1δ2+γ2σ2.The Jacobian matrix of system (34) at any equilibrium F¯x¯,y¯,z¯,u¯,v¯ is as follows:(52)JF¯=α−2αβx¯−y¯−x¯0000ρz¯−θu¯−δ1ρy¯−θy¯0ω2z¯0ω2x¯−δ2000γ1γ2−δ3−ξv¯ξu¯0000−δ4=α−2αβx¯−y¯−x¯0000ρz¯−θu¯−δ1ρy¯−θy¯0ω2z¯0ω2x¯−δ2000γ1γ2−δ31ξu¯0000−δ4.The corresponding characteristic equation is(53)JF¯−λI=0.Substitutingδ31 for δ3 in (12), we have(54)JF¯−λI=λ+δ4JP¯−λI=0.Thus, the stability analysis ofF¯ is similar to that of P¯.Substitutingδ31 for δ3 in (20) to get new conditions for the stability of F∗:(55)h1:A4=θαβγ1δ2−ω2x∗x∗y∗+ρδ31−θγ2ω2x∗y∗z∗>0,h2:A1A2−A3=θγ1αβx∗+δ2−ω2x∗y∗+αβδ31δ2−ω2x∗x∗−ρx∗y∗z∗>0,h3:A1A2A3−A1A4−A32=B1θ+B2>0,where(56)B1=αβx∗+δ2−ω2x∗2δ31+αβx∗+δ2−ω2x∗δ312−ρω2x∗y∗z∗αβx∗+δ2−ω2x∗γ1y∗−αβx∗+δ2−ω2x∗+δ312αβγ1δ2−ω2x∗x∗y∗−ω2γ2x∗y∗z∗,B2=αβx∗+δ2−ω2x∗2δ31+αβx∗+δ2−ω2x∗δ312−ρω2x∗y∗z∗αβδ31δ2−ω2x∗x∗+ρω2x∗y∗z∗−αβx∗+δ2−ω2x∗+δ312ρω2δ31x∗y∗z∗.We set(57)θ70=δ31ρσ2−δ1δ2αγ1δ2+γ2σ2,θ71=δ31ρσ2−δ1δ2γ2σ2,θ72=βδ31ρσ2−δ1δ2+δ1δ31ω2βγ2σ2,θ73=ρδ31γ2,and we obtain the following results.Theorem 7.
System (34) always has one tumor-free and ECs-free equilibrium F0, which is unstable.Theorem 8.
System (34) has one tumor-free equilibrium F1 when ρ>δ1δ2/σ2 and θ<θ71, which is LAS if the inequality θ<θ70 holds.Theorem 9.
System (34) has one tumor-dominant equilibrium F2 when ω2<ω2∗, which is LAS if the inequality θ>θ72 holds.Theorem 10.
System (34) has a unique coexistence equilibrium F∗ when ω2<ω2∗,θ0<θ<θ2 or ω2>ω2∗,θ0<θ<θ3 hold. Further F∗ is LAS, when condition (55) is satisfied.


#### 3.1.3. Combined Immunotherapy Model

Whenσ1>0andσ3>0, (36) can be simplified as follows:(58)αx1−βx−xy=0,ρyz−θyu−δ1y+σ1=0,σ2+ω2xz−δ2z=0,γ1y+γ2z−δ31u=0,where δ31=δ3δ4+ξσ3/δ4.The analysis of (58) is similar to that of (37), and the equilibria of (34) can be obtained as follows:(59)N0=x0,y0,z0,u0,v0=0,y0,σ2δ2,γ1δ2y0+γ2σ2δ2δ31,σ3δ4,where y0=−B+B2−4AC/2A with A=θγ1δ2>0,B=δ31δ1δ2−ρσ2+θγ2σ2,andC=−σ1δ2δ31<0.(60)N∗=x∗,y∗,z∗,u∗,v∗=x∗,α1−βx∗,σ2δ2−ω2x∗,αγ11−βx∗δ2−ω2x∗+γ2σ2δ31δ2−ω2x∗,σ3δ4,where x∗∈0,μ and satisfies the equation c3x3+c2x2+c1x+c0=0 with coefficients(61)c3=θα2β2γ1ω2>0,c2=−θα2βγ1ω2−αβθαγ1βδ2+ω2+δ1δ31ω2<0,c1=θα2αβγ1δ2+αγ1ω2+βγ2σ2+αδ31δ1ω2−βρσ2−δ1δ2−σ1ω2δ31,c0=σ1δ2δ31+αδ31ρσ2−δ1δ2−θααγ1δ2+γ2σ2.For system (34) at N0=x0,y0,z0,u0,v0, Jacobian matrix is expressed as follows:(62)JN0=α−y000000−σ1y0ρy0−θy00ω2z00−δ2000γ1γ2−δ31ξu00000−δ4.The corresponding characteristic equation is(63)JN0−λI=λ+δ4JE0−λI=0.Thus, the stability analysis ofN0 is similar to that of E0. Set θ8=αδ31ρσ2−δ1δ2+σ1δ2δ31/ααγ1δ2+σ2γ2, and we can get the following result.Theorem 11.
System (34) always has one tumor-free equilibrium N0, which is LAS if the inequality θ<θ8 holds.
For system (34) at N∗=x∗,y∗,z∗,u∗,v∗, Jacobian matrix is given as follows:(64)JN∗=−αβx∗−x∗0000−σ1y∗ρy∗−θy∗0ω2z∗0ω2x∗−δ2000γ1γ2−δ31ξu∗0000−δ4.
The corresponding characteristic equation is(65)JN∗−λI=λ+δ4JE∗−λI=0.
Thus, the stability analysis ofN∗ is similar to that of E∗. Set θ9=ρδ31/γ2 and substitute δ31 for δ3 in (47) to get new conditions:(66)h4:B4=αβx∗σ2z∗θγ1y∗+δ31σ1y∗+ρδ31−θγ2ω2x∗y∗z∗>0,h5:B1B2−B3=k1θ+k2>0,h6:B1B2B3−B1B4−B32=B3k1θ+k2−B12B4=k3θ2+k4θ+k5>0,where(67)k1=σ1+δ31y∗γ1,k2=δ31+σ1y∗+σ2z∗α2β2x∗2+αβx∗+σ2z∗σ1y∗+δ31+αβx∗σ2z∗+δ31σ1y∗σ1y∗+δ31−ρω2x∗y∗z∗,k3=k1αβx∗+σ2z∗γ1y∗,k4=k2γ1y∗+k1σ1δ31y∗αβx∗+σ2z∗+k1αβx∗σ2z∗δ31+σ1y∗−αβx∗+δ31+σ1y∗+σ2z∗2γ1y∗+ρ+γ2ω2x∗y∗z∗,k5=k2αβx∗+σ2z∗σ1δ31y∗−αβx∗σ2z∗k2δ31+σ1y∗−αβx∗+δ31+σ1y∗+σ2z∗2σ1δ31y∗+k2−δ31ρω2x∗y∗z∗.
Therefore, we can obtain the following result.Theorem 12.
System (34) has a coexistence equilibrium N∗ when ω2>βδ2,σ1≥ααργ1+δ1γ2/γ2,θ9<θ<θ8 or ω2>βδ2, σ1≤ααργ1+δ1γ2/γ2,θ8<θ<θ9 or ω2<βδ2,θ>θ8 hold. Further N∗ is LAS when condition (66) is satisfied.


### 3.2. Numerical Simulations

#### 3.2.1. ACI Model (σ1>0,σ3=0)

We conduct numerical simulations of the ACI model. We chooseω2=0.0004 and the other parameter values are given in (22). We find the stability region of E∗=x∗,y∗,z∗,u∗, as shown in Figure 5(a). E∗ is stable in region I and unstable in region II.Case (a): we choose a pointG1=σ1,θ=0.4,0.5 in the region I, then system (34) has one coexistence equilibrium:(68)E∗=11.560,1.598,7.543,6.993.The eigenvalues of characteristic equation of (47) are(69)−0.269,−0.041,−0.1388+0.334iand−0.1388−0.334i,soE∗ is stable, as shown in Figure 6(b).Case (b): we choose a pointG2=σ1,θ=0.4,0.55 in the region II; then, system (34) has one coexistence equilibrium:(70)E∗=69.131,1.409,13.895,11.962.The eigenvalues of characteristic equation of (47) are(71)−0.830,−0.021,0.0326+0.6039iand0.0326−0.6039i,soE∗ is unstable, as shown in Figures 6(c) and 6(d).Figure 5
(a) Stability region ofE∗ and Hopf bifurcation curve in θ−σ1 parameter plane. E∗ is stable in region I and unstable in region II. (b) ACI curve. E0 is stable in region III and unstable in region IV. (c) The bifurcation diagram of x with respect to θ for fixed σ1=0.4. Hopf bifurcation point appears at θACI∗=0.5354. (d) The bifurcation diagram of x with respect to θ for different σ1=0,0.1,0.3,0.5,0.7,0.9 (the blue curves represent the stable steady states of system (34), while the red curves represent the unstable steady states).
(a)(b)(c)(d)Figure 6
Chooseσ1=0.4, and the other parameter values are given in (22). (a) When θ=0.48, E00,1.782,6.909,6.596 is LAS. (b) When θ=0.50, E∗11.560,1.598,7.543,6.993 is LAS. (c) When θ=0.55, E∗69.131,1.409,13.895,11.962 is unstable. (d) When θ=0.55, the 3D phase portrait depicts TCs, ECs, and Tregs.
(a)(b)(c)(d)According to the stability condition of tumor-free equilibriumE0 in Theorem 5, ACI curve can be obtained as(72)θ=αδ3ρσ2−δ1δ2+σ1δ2δ3ααγ1δ2+γ2σ2=0.09391σ1+0.452.We find the stability region ofE0=0,y0,z0,u0 shown in Figure 5(b). E0 is stable in the region III and unstable in the region IV. With the increase of σ1, region III gradually increases, while region IV gradually decreases.Case (c): we choose a pointG3=σ1,θ=0.4,0.48 in the region III; then, system (34) has a tumor-free equilibrium:(73)E0=0,1.782,6.909,6.596.The eigenvalues of Jacobian matrix ofE0 are(74)−0.055,−0.146,−0.237+0.358iand−0.237−0.358i,soE0 is stable, as shown in Figure 6(a).Case (d): we choose a pointG4=σ1,θ=0.4,0.53 in the region IV; then, system (34) has a tumor-free equilibrium:(75)E0=0,1.141,6.909,6.212.The eigenvalues of Jacobian matrix ofE0 are(76)0.494,−0.055,−0.3+0.297iand−0.3−0.297i,so E0 is unstable.We choose ACI parameterσ1=0.4 to study the relationship between the number of TCs x and the parameter θ (see Figure 5(c)). By some calculations, we have θ3=ρδ3/γ2=0.6,θ5=αδ3ρσ2−δ1δ2+σ1δ2δ3/ααγ1δ2+γ2σ2=0.489. And Hopf bifurcation point appears at θACI∗=0.5354. We can obtain the following result.Proposition 5.
When0<θ<θ5, E0 exists and is LAS, and E∗ is nonexistent. When θ5<θ<θACI∗, E0 exists and is unstable, and E∗ exists and is LAS. When θACI∗<θ<θ3, E0 exists and is unstable, and E∗ exists and is unstable. When θ>θ3, E0 exists and is unstable, and E∗ is nonexistent.
Next, we choose different ACI parametersσ1=0,0.1,0.3,0.5,0.7,0.9, to study the relationship between the number of TCs x and the parameter θ (see Figure 5(d)).
By comparing the curves in Figure5(d), we see that, as σ1 increases gradually, the stable region of tumor-free equilibrium E0 of system (34) gradually increases (the intersection point of the curve and the x-coordinate gradually moves to the right) and the stability region of equilibrium E∗ in system (34) gradually increases (the blue curve gradually moves upward). This indicates that increasing the injection volume of adoptive T cells can not only delay the inhibitory effect of Tregs on tumor immune response but also help the immune system to remove more TCs. It also helps the immune system to control more TCs, keeping them at a stable state.


#### 3.2.2. MAI Model (σ1=0,σ3>0)

We conduct numerical simulations of the MAI model. We chooseω2=0.0004,ξ=0.05,andδ4=0.25 and the other parameter values are given in (22). We find the stability region of F∗=x∗,y∗,z∗,u∗, as shown in Figure 7(a). F∗ is stable in region I and unstable in region II.Case (a): we choose a pointM1=σ3,θ=0.4,0.6 in region I; then, system (34) has one coexistence equilibrium:(77)F∗=2.413,1.628,7.032,5.002,1.6.The eigenvalues of Jacobian matrix ofF∗ are(78)−0.25,−0.146+0.336i,−0.146−0.336i,−0.0499+0.033iand−0.0499−0.033i,soF∗ is stable, as shown in Figure 8(b).Case (b): we choose a pointM2=σ3,θ=0.4,0.635 in region II; then, system (34) has one coexistence equilibrium:(79)F∗=27.345,1.546,8.624,5.93,1.6.The eigenvalues of Jacobian matrix ofF∗ are(80)−0.4334,−0.25,−0.0644,0.0172+0.429iand0.0172−0.429i,soF∗ is unstable, as shown in Figures 8(c) and 8(d).According to the stability condition of tumor-free equilibriumF1 in Theorem 7, MAI curve can be obtained as(81)θ=δ3δ4+ξσ3ρσ2−δ1δ2/αγ1δ2+σ2γ2δ4=0.3616σ3+0.452.We find the stability region ofF1, as shown in Figure 7(b).F0 is stable in region III and unstable in region IV. With the increase of σ3, region III gradually increases, while region IV gradually decreases.Case (c): we choose a pointM3=σ3,θ=0.4,0.5 in region III; then, system (34) has a tumor-free equilibrium:(82)F1=0,2.556,6.909,5.349,1.6.The eigenvalues of Jacobian matrix ofF1 are(83)−0.92,−0.25,−0.055,−0.164+0.428iand−0.164−0.428i,soF1 is stable, as shown in Figure 8(a).Case (d): we choose a pointM4=σ3,θ=0.4,0.6 in region IV; then, system (34) has a tumor-free equilibrium:(84)F1=0,1.575,6.909,4.903,1.6.Figure 7
(a) Stability region ofF∗ and Hopf bifurcation curve in θ−σ3 parameter plane. F∗ is stable in region I and unstable in region II. (b) MAI curve. F1 is stable in region III and unstable in region IV. (c) The bifurcation diagram of x with respect to θ for fixed σ3=0.4. Hopf bifurcation point appears at θMAI∗=0.627. (d) The bifurcation diagrams of x with respect to θ for different σ3=0,0.1,0.3,0.5,0.7,0.9 (the blue curves represent the stable steady states of system (34), while the red curves represent the unstable steady states).
(a)(b)(c)(d)Figure 8
Chooseσ3=0.4, and the other parameter values are given in (22). (a) When θ=0.5, F10,2.556,6.909,5.349,1.6 is LAS. (b) When θ=0.6, F∗2.413,1.628,7.032,5.002,1.6 is LAS. (c) When θ=0.635, F∗27.345,1.546,8.624,5.93,1.6 is unstable. (d) When θ=0.635, the 3D phase portrait depicts TCs, ECs, and Treg.
(a)(b)(c)(d)The eigenvalues of Jacobian matrix ofF1 are(85)0.061,−0.25,−0.055,−0.164+0.338iand−0.164−0.338i,so F1 is unstable.We choose MAI parameterσ3=0.4 to study the relationship between the number of TCs x and the parameter θ (see Figure 7(c)). By some calculations, we have θ70=δ31ρσ2−δ1δ2/αγ1δ2+γ2σ2=0.597, θ71=δ31ρσ2−δ1δ2/γ2σ2=0.703, and θ73=ρδ31/γ2=0.792. And Hopf bifurcation point appears at θMAI∗=0.627. We can obtain the following result.Proposition 6.
When0<θ<θ70, F0 exists and is unstable, F1 exists and is LAS, and F∗ exists and is LAS. When θ70<θ<θMAI∗, F0 exists and is unstable, F1 exists and is unstable, and F∗ exists and is unstable. When θMAI∗<θ<θ71, F0 exists and is unstable, F1 exists and is unstable, and F∗ exists and is unstable. When θ71<θ<θ73, F0 exists and is unstable, F1 is nonexistent, and F∗ exists and is unstable. When θ>θ3, F0 exists and is unstable, F1 is nonexistent, and F∗ is nonexistent.
Next, we choose different MAI parametersσ3=0,0.1,0.3,0.5,0.7,0.9, to study the relationship between the number of TCs x and the parameter θ (see Figure 7(d)).
By comparing the curves in Figure7(d), we can find that, with the gradual increase of σ3, the stability region of the tumour-free equilibrium F1 of system (34) gradually increases (the intersection point of curves and x-coordinate gradually moves to the right). This shows that the increase of antibody injection quantity can help to slow down Tregs inhibition of tumor immune responses. The stability region of equilibrium F∗ of system (34) gradually increases (the blue curve gradually moves upward), which means that increasing the amount of antibody injected can help the immune system to control more TCs. By comparing Figures 5 and 7, we find that the effect of MAI is better than that of AIC in delaying the inhibitory effect of Tregs on tumor immune response (at the same injection dose, the intersection point of curves and x-coordinate in Figure 7 moves to the right more widely than that in Figure 5). AIC is more effective than MAI in controlling TCs (the blue curve in Figure 7 moves up more than the blue curves in Figure 5 at the same injection dose).


#### 3.2.3. Combined Immunotherapy Model (σ1>0,σ3>0)

We conduct numerical simulations of the combined immunotherapy model. We chooseω2=0.0004,ξ=0.05,andδ4=0.25, and the other parameter values are given in (22). We find the stability region of N∗=x∗,y∗,z∗,u∗,v∗, as shown in Figure 9(a). When σ3 increases gradually, the stability region of N∗ also increases gradually.Figure 9
(a) Choose differentσ3=0.1,0.5,0.8. Stability region of N∗ and Hopf bifurcation curve in θ−σ1 parameter plane. N∗ is stable to the left of the curve and unstable to the right of the curve. (b) Combined immunotherapy surface. N0 is stable below the surface and unstable above the surface. (c) When σ1=0.4andσ3=0.4, the bifurcation diagram of x with respect to θ for system (34). Hopf bifurcation point appears at θCI∗=0.727. (d) Choose different σ1,σ3=0.4,0,0,0.4,0.4,0.4,0.7,0.7. The bifurcation diagram of x with respect to θ (the stable steady state is represented by the blue curve, and the unstable one corresponds to the red curve).
(a)(b)(c)(d)Figure 10
(a) Whenσ1=σ3=0.5,θ=0.4, N00,8.3,6.909,7.505,2 is LAS. (b) When σ1=0.4,σ3=0.5,θ=0.7, N∗11.56,1.59,7.54,4.99,2 is LAS. (c) When σ1=0.4,σ3=0.5,θ=0.78, N∗77.887,1.381,15.936,9.698,2 is unstable. (d) When σ1=0.4,σ3=0.5,θ=0.78, the 3D phase portrait depicts TCs, ECs, and Tregs.
(a)(b)(c)(d)Case (a): we choose a pointL1=θ,σ1=0.7,0.4 on the left of the curve σ3=0.5; then, system (34) has one coexistence equilibrium:(86)N∗=11.56,1.59,7.54,4.99,2.The eigenvalues of Jacobian matrix of (64) are(87)−0.06,−0.176,−0.25,−0.22−0.37iand−0.22+0.37i,soN∗ is stable, as shown in Figure 10(b).Case (b): we choose a pointL2=θ,σ1=0.78,0.4 on the right of the curve σ3=0.5; then, system (34) has one coexistence equilibrium:(88)N∗=77.887,1.381,15.936,9.698,2.The eigenvalues of Jacobian matrix of (64) are(89)−0.906238,−0.0248,−0.25,0.00641−0.66iand0.00641+0.66i,soN∗ is unstable, as shown in Figures 10(c) and 10(d).According to the stability condition of tumor-free equilibriumN0 in Theorem 11, combined immunotherapy surface can be obtained as(90)θ=αδ31ρσ2−δ1δ2+σ1δ2δ31ααγ1δ2+σ2γ2=0.055σ1+0.26471.366σ3+1.707.Therefore, the stability region of tumor-free equilibriumN0 of system (34) is obtained (see Figure 9(b)).Case (c): we choose a pointL3=σ1,σ2,θ=0.5,0.5,0.4 below surface (90); then, system (34) has a tumor-free equilibrium:(91)N0=0,8.3,6.909,7.505,2.The eigenvalues of Jacobian matrix of (64) are(92)−6.664,−0.055,−0.25,−0.204+0.69iand−0.204−0.69i,soN0 is stable, as shown in Figure 10(a).Case (d): we choose a pointL4=σ1,σ2,θ=0.5,0.5,0.8 above surface (90); then, system (34) has a tumor-free equilibrium:(93)N0=0,0.886,6.909,4.32,2.The eigenvalues of Jacobian matrix of (64) are(94)0.75,−0.055,−0.25,−0.453+0.311iand−0.453−0.311i,so N0 is unstable.We choose combined immunotherapy parameterσ1=0.4andσ3=0.4 to study the relationship between the number of TCs x and the parameter θ (see Figure 9(c)). By some calculations, we have θ8=αδ31ρσ2−δ1δ2+σ1δ2δ31/ααγ1δ2+σ2γ2=0.646,andθ9=ρδ31/γ2=0.791. And Hopf bifurcation point appears at θCI∗=0.727. We can obtain the following result.Proposition 7.
When0<θ<θ8, N0 exists and is LAS, and N∗ is nonexistent. When θ8<θ<θCI∗, N0 exists and is unstable, and N∗ exists and is LAS. When θCI∗<θ<θ9, N0 exists and is unstable, and N∗ exists and is unstable. When θ>θ9, N0 exists and is unstable, and N∗ is nonexistent.
For combined immunotherapy, by comparing the curves in Figure9(d), we know that, with the increase of σ1 and σ3, the stability region of the tumour-free equilibrium of N0 of system (34) gradually increases (the intersection point of curves and x-coordinate gradually moves to the right), and the stability region of equilibrium N∗ of system (34) gradually increases (the blue curve gradually moves upward). By comparing Figure 9 with Figures 5 and 7, it can be seen that combined immunotherapy has better effects on delaying the inhibitory effect of Tregs on tumor immune response and helps the immune system to control more TCs than ACI or MAI.


## 3.1. Model Analysis

By a similar proof of Proposition1, we can obtain the well posedness of model (34) with initial conditions (35) as follows.Proposition 4.
The solutionsxt,yt,zt,ut of model (34) with initial conditions (35) are existent, unique, and nonnegative on the interval 0,+∞.Next, we discuss the effect of three types of immunotherapy on tumors such asCase I: ACI model (σ1>0,σ3=0)Case II: MAI model (σ1=0,σ3>0)Case III: combined immunotherapy model (σ1>0,σ3>0)We setdx/dt=0,dy/dt=0,dz/dt=0,du/dt=0,anddv/dt=0 in (34), and we have(36)αx1−βx−xy=0,ρyz−θyu−δ1y+σ1=0,σ2+ω2xz−δ2z=0,γ1y+γ2z−δ3u−ξuv=0,σ3−δ4v=0.


### 3.1.1. ACI Model

Whenσ1>0andσ3=0, we put v=0 and (36) becomes(37)αx1−βx−xy=0,ρyz−θyu−δ1y+σ1=0,σ2+ω2xz−δ2z=0,γ1y+γ2z−δ3u=0.Puttingx=0 yields the tumor-free equilibrium, namely,(38)E0=x0,y0,z0,u0=0,y0,σ2δ2,γ1δ2y0+γ2σ2δ2δ3,where y0=−B+B2−4AC/2A, A=θγ1δ2>0,B=δ1δ2δ3+θγ2σ2−ρσ2δ3,andC=−σ1δ2δ3<0. E0 always exists.Puttingx≠0 and eliminating y,z,andu in (37), we have(39)Fx=C3x3+C2x2+C1x+C0=0,C3=θα2β2γ1ω2>0,C2=−θα2βγ1ω2−αβθαγ1βδ2+ω2+δ1δ3ω2<0,C1=θα2αβγ1δ2+αγ1ω2+βγ2σ2+αδ3δ1ω2−βρσ2−δ1δ2−σ1ω2δ3,C0=σ1δ2δ3+αδ3ρσ2−δ1δ2−θααγ1δ2+γ2σ2.Case (a): ifF0Fμ=C0Fμ<0 (where μ=min1/β,δ2/ω2, due to F+∞>0 and F−∞<0, then (39) has one positive root x∗ in the interval 0,μ. Therefore, system (34) has a unique coexistence equilibrium:(40)E∗=x∗,y∗,z∗,u∗=x∗,α1−βx∗,σ2δ2−ω2x∗,αγ11−βx∗δ2−ω2x∗+γ2σ2δ3δ2−ω2x∗.We show the existence conditions ofE∗ in the following:Case (1):C0>0,Fμ<0.SinceC0>0, we have θ<σ1δ2δ3+αδ3ρσ2−δ1δ2/ααγ1δ2+γ2σ2≜θ5. Since Fμ<0, we have the following two cases:Case where1/β<δ2/ω2, i.e., ω2<βδ2. Since μ=1/β, we have Fμ=δ3σ1βδ2−ω2/β>0. That is a contradiction.Case where1/β>δ2/ω2, i.e., ω2>βδ2. Since μ=δ2/ω2, we have Fμ=ασ2βδ2−ω2γ2θ−δ3ρ/ω2<0 implying θ>ρδ3/γ2=θ3. If θ5≥θ3, i.e., σ1δ2δ3+αδ3ρσ2−δ1δ2/ααγ1δ2+γ2σ2≥ρδ3/γ2, we have σ1≥ααργ1+δ1γ2/γ2, there exists E∗.Case (2):C0<0,Fμ>0.SinceC0<0, we have θ>θ5. Since Fμ>0, we have the following two cases:Case where1/β<δ2/ω2, i.e., ω2<βδ2. Since μ=1/β, we have Fμ=δ3σ1βδ2−ω2/β>0, there exists E∗.Case where1/β>δ2/ω2, i.e., ω2>βδ2. Since μ=δ2/ω2, we have Fμ=ασ2βδ2−ω2γ2θ−δ3ρ/ω2>0 implying θ<θ3. If θ3≥θ5, i.e., ρδ3/γ2≥σ1δ2δ3+αδ3ρσ2−δ1δ2/ααγ1δ2+γ2σ2, we have σ1≤ααργ1+δ1γ2/γ2, there exists E∗.Case (b): IfF0Fμ=C0Fμ>0, then (39) can have two positive roots x1∗ and x2∗ in the interval 0,μ, where x1∗ and x2∗ are the roots of (39). Therefore, system (34) can have two coexistence equilibria:(41)E1∗=x1∗,y1∗,z1∗,u1∗=x1∗,α1−βx1∗,σ2δ2−ω2x1∗,αγ11−βx1∗δ2−ω2x1∗+γ2σ2δ3δ2−ω2x1∗,(42)E2∗=x2∗,y2∗,z2∗,u2∗=x2∗,α1−βx2∗,σ2δ2−ω2x2∗,αγ11−βx2∗δ2−ω2x2∗+γ2σ2δ3δ2−ω2x2∗.In the following, we show the existence condition ofE1∗ and E2∗.Differentiating (39), we have(43)F′x=3C3x2+2C2x+C1=0,where C3>0 and C2<0. If C1>0 and Δ≜C22−3C1C3>0, then (43) has two positive roots x1+=−C2−Δ/3C3 and x2+=−C2+Δ/3C3. And if C1<0, then (43) has one positive root x2+=−C2+Δ/3C3.Case (1):C0>0,Fμ>0.SinceC0>0, we have θ<θ5. Since Fμ>0, we have the following two cases:Case where1/β<δ2/ω2, i.e., ω2<βδ2. Since μ=1/β, we have Fμ=δ3σ1βδ2−ω2/β>0. If Fx2+<0, then there exist E1∗ and E2∗.Case where1/β>δ2/ω2, i.e., ω2>βδ2. Since μ=δ2/ω2, we have Fμ=ασ2βδ2−ω2γ2θ−δ3ρ/ω2>0 implying θ<θ3. Then, we have θ<minθ3,θ5. If Fx2+<0, then there exist E1∗ and E2∗.Case (2):C0<0,Fμ<0.SinceC0<0, we have θ>θ5. Since Fμ>0, we have the following two cases:Case where1/β<δ2/ω2, i.e., ω2<βδ2. Since μ=1/β, we have Fμ=δ3σ1βδ2−ω2/β>0. That is a contradiction.Case where1/β>δ2/ω2, i.e., ω2>βδ2. Since μ=δ2/ω2, we have Fμ=ασ2βδ2−ω2γ2θ−δ3ρ/ω2<0 implying θ>θ3. Then, we have θ>maxθ3,θ5. If Fx1+>0, then there exist E1∗ and E2∗.Next, we show the stability of the above two equilibriaE0 and E∗ of system (34) with v=0. We linearize the system and obtain the characteristic equation whose expression is the same as (13).Theorem 5.
System (34) with v=0 always has one tumor-free equilibrium E0, which is LAS if the inequality θ<θ5 holds.Proof.
AtE0=x0,y0,z0,u0=0,y0,σ2/δ2,γ1δ2y0+γ2σ2/δ2δ3, the characteristic equation becomes(44)λ−α+y0λ+δ2λ2+σ1y0+δ3λ+θγ1y0+σ1δ3y0=0,and one root of characteristic equation is −δ2<0. It is easily noted that as σ1/y0+δ3>0,θγ1y1+σ1δ3/y0>0, so the solutions of λ2+σ1/y0+δ3λ+θγ1y1+σ1δ3/y0=0 have always negative real parts. Therefore, E0 is LAS if and only if α−y0=2αA+B−B2−4AC/2A<0, that is, if and only if θ<θ5.
AtE∗=x∗,y∗,z∗,u∗, characteristic (13) becomes(45)λ4+B1λ3+B2λ2+B3λ+B4=0,where(46)B1=αβx∗+δ3+σ1y∗+σ2z∗>0,B2=αβx∗+σ2z∗σ1y∗+δ3+αβx∗σ2z∗+δ3σ1y∗+θγ1y∗>0,B3=θγ1y∗+δ3σ1y∗αβx∗+σ2z∗+αβx∗σ2z∗δ3+σ1y∗+ρω2x∗y∗z∗>0,B4=αβx∗σ2z∗θγ1y∗+δ3σ1y∗+ρδ3−θγ2ω2x∗y∗z∗.
By the Routh–Hurwitz criterion, the roots of (45) have negative real parts if and only if B1>0,B2>0,B3>0,B4>0,B1B2−B3>0,andB1B2B3−B1B4−B32>0.
Hence, we obtain the sufficient condition for the stability ofP∗:(47)H4:B4=αβx∗σ2z∗θγ1y∗+δ3σ1y∗+ρδ3−θγ2ω2x∗y∗z∗>0,H5:B1B2−B3=K1θ+K2>0,H6:B1B2B3−B1B4−B32=B3K1θ+K2−B12B4=K3θ2+K4θ+K5>0,where(48)K1=σ1+δ3y∗γ1,K2=δ3+σ1y∗+σ2z∗α2β2x∗2+αβx∗+σ2z∗σ1y∗+δ3+αβx∗σ2z∗+δ3σ1y∗σ1y∗+δ3−ρω2x∗y∗z∗,K3=K1αβx∗+σ2z∗γ1y∗,K4=K2γ1y∗+K1σ1δ3y∗αβx∗+σ2z∗+K1αβx∗σ2z∗δ3+σ1y∗−αβx∗+δ3+σ1y∗+σ2z∗2γ1y∗+ρ+γ2ω2x∗y∗z∗,K5=K2αβx∗+σ2z∗σ1δ3y∗−αβx∗σ2z∗K2δ3+σ1y∗−αβx∗+δ3+σ1y∗+σ2z∗2σ1δ3y∗+K2−δ3ρω2x∗y∗z∗.Theorem 6.
System (34) with v=0 has a coexistence equilibrium E∗ when ω2>βδ2,σ1≥ααργ1+δ1γ2/γ2,θ3<θ<θ5 or ω2>βδ2, σ1≤ααργ1+δ1γ2/γ2,θ5<θ<θ3 or ω2<βδ2,θ>θ5 hold. Further E∗ is LAS when condition (47) is satisfied.


### 3.1.2. MAI Model

Whenσ1=0,σ3>0, (36) becomes(49)αx1−βx−xy=0,ρyz−θyu−δ1y=0,σ2+ω2xz−δ2z=0,γ1y+γ2z−δ3u−ξuv=0,σ3−δ4v=0.Equation (49) can be simplified as follows:(50)αx1−βx−xy=0,ρyz−θyu−δ1y=0,σ2+ω2xz−δ2z=0,γ1y+γ2z−δ31u=0,where δ31=δ3δ4+ξσ3/δ4.The analysis of (50) is similar to that of (12), and the equilibria of (34) can be obtained as follows:(51)F0=0,0,σ2δ2,γ2σ2δ2δ31,σ3δ4=0,δ31ρσ2−δ1δ2−θγ2σ2θγ1δ2,σ2δ2,ρσ2−δ1δ2θδ2,σ3δ4,F2=1β,0,βσ2βδ2−ω2,βγ2σ2δ31βδ2−ω2,σ3δ4,F∗=x∗,α1−βx∗,σ2δ2−ω2x∗,αγ11−βx∗δ2−ω2x∗+γ2σ2δ31δ2−ω2x∗,σ3δ4.Here,x∗∈0,μ and satisfies the equation a2x2+a1x+a0=0, where coefficients are defined as a2=−θαβγ1ω2<0,a1=θαγ1βδ2+ω2+δ1δ31ω2>0,anda0=δ31ρσ2−δ1δ2−θαγ1δ2+γ2σ2.The Jacobian matrix of system (34) at any equilibrium F¯x¯,y¯,z¯,u¯,v¯ is as follows:(52)JF¯=α−2αβx¯−y¯−x¯0000ρz¯−θu¯−δ1ρy¯−θy¯0ω2z¯0ω2x¯−δ2000γ1γ2−δ3−ξv¯ξu¯0000−δ4=α−2αβx¯−y¯−x¯0000ρz¯−θu¯−δ1ρy¯−θy¯0ω2z¯0ω2x¯−δ2000γ1γ2−δ31ξu¯0000−δ4.The corresponding characteristic equation is(53)JF¯−λI=0.Substitutingδ31 for δ3 in (12), we have(54)JF¯−λI=λ+δ4JP¯−λI=0.Thus, the stability analysis ofF¯ is similar to that of P¯.Substitutingδ31 for δ3 in (20) to get new conditions for the stability of F∗:(55)h1:A4=θαβγ1δ2−ω2x∗x∗y∗+ρδ31−θγ2ω2x∗y∗z∗>0,h2:A1A2−A3=θγ1αβx∗+δ2−ω2x∗y∗+αβδ31δ2−ω2x∗x∗−ρx∗y∗z∗>0,h3:A1A2A3−A1A4−A32=B1θ+B2>0,where(56)B1=αβx∗+δ2−ω2x∗2δ31+αβx∗+δ2−ω2x∗δ312−ρω2x∗y∗z∗αβx∗+δ2−ω2x∗γ1y∗−αβx∗+δ2−ω2x∗+δ312αβγ1δ2−ω2x∗x∗y∗−ω2γ2x∗y∗z∗,B2=αβx∗+δ2−ω2x∗2δ31+αβx∗+δ2−ω2x∗δ312−ρω2x∗y∗z∗αβδ31δ2−ω2x∗x∗+ρω2x∗y∗z∗−αβx∗+δ2−ω2x∗+δ312ρω2δ31x∗y∗z∗.We set(57)θ70=δ31ρσ2−δ1δ2αγ1δ2+γ2σ2,θ71=δ31ρσ2−δ1δ2γ2σ2,θ72=βδ31ρσ2−δ1δ2+δ1δ31ω2βγ2σ2,θ73=ρδ31γ2,and we obtain the following results.Theorem 7.
System (34) always has one tumor-free and ECs-free equilibrium F0, which is unstable.Theorem 8.
System (34) has one tumor-free equilibrium F1 when ρ>δ1δ2/σ2 and θ<θ71, which is LAS if the inequality θ<θ70 holds.Theorem 9.
System (34) has one tumor-dominant equilibrium F2 when ω2<ω2∗, which is LAS if the inequality θ>θ72 holds.Theorem 10.
System (34) has a unique coexistence equilibrium F∗ when ω2<ω2∗,θ0<θ<θ2 or ω2>ω2∗,θ0<θ<θ3 hold. Further F∗ is LAS, when condition (55) is satisfied.


### 3.1.3. Combined Immunotherapy Model

Whenσ1>0andσ3>0, (36) can be simplified as follows:(58)αx1−βx−xy=0,ρyz−θyu−δ1y+σ1=0,σ2+ω2xz−δ2z=0,γ1y+γ2z−δ31u=0,where δ31=δ3δ4+ξσ3/δ4.The analysis of (58) is similar to that of (37), and the equilibria of (34) can be obtained as follows:(59)N0=x0,y0,z0,u0,v0=0,y0,σ2δ2,γ1δ2y0+γ2σ2δ2δ31,σ3δ4,where y0=−B+B2−4AC/2A with A=θγ1δ2>0,B=δ31δ1δ2−ρσ2+θγ2σ2,andC=−σ1δ2δ31<0.(60)N∗=x∗,y∗,z∗,u∗,v∗=x∗,α1−βx∗,σ2δ2−ω2x∗,αγ11−βx∗δ2−ω2x∗+γ2σ2δ31δ2−ω2x∗,σ3δ4,where x∗∈0,μ and satisfies the equation c3x3+c2x2+c1x+c0=0 with coefficients(61)c3=θα2β2γ1ω2>0,c2=−θα2βγ1ω2−αβθαγ1βδ2+ω2+δ1δ31ω2<0,c1=θα2αβγ1δ2+αγ1ω2+βγ2σ2+αδ31δ1ω2−βρσ2−δ1δ2−σ1ω2δ31,c0=σ1δ2δ31+αδ31ρσ2−δ1δ2−θααγ1δ2+γ2σ2.For system (34) at N0=x0,y0,z0,u0,v0, Jacobian matrix is expressed as follows:(62)JN0=α−y000000−σ1y0ρy0−θy00ω2z00−δ2000γ1γ2−δ31ξu00000−δ4.The corresponding characteristic equation is(63)JN0−λI=λ+δ4JE0−λI=0.Thus, the stability analysis ofN0 is similar to that of E0. Set θ8=αδ31ρσ2−δ1δ2+σ1δ2δ31/ααγ1δ2+σ2γ2, and we can get the following result.Theorem 11.
System (34) always has one tumor-free equilibrium N0, which is LAS if the inequality θ<θ8 holds.
For system (34) at N∗=x∗,y∗,z∗,u∗,v∗, Jacobian matrix is given as follows:(64)JN∗=−αβx∗−x∗0000−σ1y∗ρy∗−θy∗0ω2z∗0ω2x∗−δ2000γ1γ2−δ31ξu∗0000−δ4.
The corresponding characteristic equation is(65)JN∗−λI=λ+δ4JE∗−λI=0.
Thus, the stability analysis ofN∗ is similar to that of E∗. Set θ9=ρδ31/γ2 and substitute δ31 for δ3 in (47) to get new conditions:(66)h4:B4=αβx∗σ2z∗θγ1y∗+δ31σ1y∗+ρδ31−θγ2ω2x∗y∗z∗>0,h5:B1B2−B3=k1θ+k2>0,h6:B1B2B3−B1B4−B32=B3k1θ+k2−B12B4=k3θ2+k4θ+k5>0,where(67)k1=σ1+δ31y∗γ1,k2=δ31+σ1y∗+σ2z∗α2β2x∗2+αβx∗+σ2z∗σ1y∗+δ31+αβx∗σ2z∗+δ31σ1y∗σ1y∗+δ31−ρω2x∗y∗z∗,k3=k1αβx∗+σ2z∗γ1y∗,k4=k2γ1y∗+k1σ1δ31y∗αβx∗+σ2z∗+k1αβx∗σ2z∗δ31+σ1y∗−αβx∗+δ31+σ1y∗+σ2z∗2γ1y∗+ρ+γ2ω2x∗y∗z∗,k5=k2αβx∗+σ2z∗σ1δ31y∗−αβx∗σ2z∗k2δ31+σ1y∗−αβx∗+δ31+σ1y∗+σ2z∗2σ1δ31y∗+k2−δ31ρω2x∗y∗z∗.
Therefore, we can obtain the following result.Theorem 12.
System (34) has a coexistence equilibrium N∗ when ω2>βδ2,σ1≥ααργ1+δ1γ2/γ2,θ9<θ<θ8 or ω2>βδ2, σ1≤ααργ1+δ1γ2/γ2,θ8<θ<θ9 or ω2<βδ2,θ>θ8 hold. Further N∗ is LAS when condition (66) is satisfied.


## 3.1.1. ACI Model

Whenσ1>0andσ3=0, we put v=0 and (36) becomes(37)αx1−βx−xy=0,ρyz−θyu−δ1y+σ1=0,σ2+ω2xz−δ2z=0,γ1y+γ2z−δ3u=0.Puttingx=0 yields the tumor-free equilibrium, namely,(38)E0=x0,y0,z0,u0=0,y0,σ2δ2,γ1δ2y0+γ2σ2δ2δ3,where y0=−B+B2−4AC/2A, A=θγ1δ2>0,B=δ1δ2δ3+θγ2σ2−ρσ2δ3,andC=−σ1δ2δ3<0. E0 always exists.Puttingx≠0 and eliminating y,z,andu in (37), we have(39)Fx=C3x3+C2x2+C1x+C0=0,C3=θα2β2γ1ω2>0,C2=−θα2βγ1ω2−αβθαγ1βδ2+ω2+δ1δ3ω2<0,C1=θα2αβγ1δ2+αγ1ω2+βγ2σ2+αδ3δ1ω2−βρσ2−δ1δ2−σ1ω2δ3,C0=σ1δ2δ3+αδ3ρσ2−δ1δ2−θααγ1δ2+γ2σ2.Case (a): ifF0Fμ=C0Fμ<0 (where μ=min1/β,δ2/ω2, due to F+∞>0 and F−∞<0, then (39) has one positive root x∗ in the interval 0,μ. Therefore, system (34) has a unique coexistence equilibrium:(40)E∗=x∗,y∗,z∗,u∗=x∗,α1−βx∗,σ2δ2−ω2x∗,αγ11−βx∗δ2−ω2x∗+γ2σ2δ3δ2−ω2x∗.We show the existence conditions ofE∗ in the following:Case (1):C0>0,Fμ<0.SinceC0>0, we have θ<σ1δ2δ3+αδ3ρσ2−δ1δ2/ααγ1δ2+γ2σ2≜θ5. Since Fμ<0, we have the following two cases:Case where1/β<δ2/ω2, i.e., ω2<βδ2. Since μ=1/β, we have Fμ=δ3σ1βδ2−ω2/β>0. That is a contradiction.Case where1/β>δ2/ω2, i.e., ω2>βδ2. Since μ=δ2/ω2, we have Fμ=ασ2βδ2−ω2γ2θ−δ3ρ/ω2<0 implying θ>ρδ3/γ2=θ3. If θ5≥θ3, i.e., σ1δ2δ3+αδ3ρσ2−δ1δ2/ααγ1δ2+γ2σ2≥ρδ3/γ2, we have σ1≥ααργ1+δ1γ2/γ2, there exists E∗.Case (2):C0<0,Fμ>0.SinceC0<0, we have θ>θ5. Since Fμ>0, we have the following two cases:Case where1/β<δ2/ω2, i.e., ω2<βδ2. Since μ=1/β, we have Fμ=δ3σ1βδ2−ω2/β>0, there exists E∗.Case where1/β>δ2/ω2, i.e., ω2>βδ2. Since μ=δ2/ω2, we have Fμ=ασ2βδ2−ω2γ2θ−δ3ρ/ω2>0 implying θ<θ3. If θ3≥θ5, i.e., ρδ3/γ2≥σ1δ2δ3+αδ3ρσ2−δ1δ2/ααγ1δ2+γ2σ2, we have σ1≤ααργ1+δ1γ2/γ2, there exists E∗.Case (b): IfF0Fμ=C0Fμ>0, then (39) can have two positive roots x1∗ and x2∗ in the interval 0,μ, where x1∗ and x2∗ are the roots of (39). Therefore, system (34) can have two coexistence equilibria:(41)E1∗=x1∗,y1∗,z1∗,u1∗=x1∗,α1−βx1∗,σ2δ2−ω2x1∗,αγ11−βx1∗δ2−ω2x1∗+γ2σ2δ3δ2−ω2x1∗,(42)E2∗=x2∗,y2∗,z2∗,u2∗=x2∗,α1−βx2∗,σ2δ2−ω2x2∗,αγ11−βx2∗δ2−ω2x2∗+γ2σ2δ3δ2−ω2x2∗.In the following, we show the existence condition ofE1∗ and E2∗.Differentiating (39), we have(43)F′x=3C3x2+2C2x+C1=0,where C3>0 and C2<0. If C1>0 and Δ≜C22−3C1C3>0, then (43) has two positive roots x1+=−C2−Δ/3C3 and x2+=−C2+Δ/3C3. And if C1<0, then (43) has one positive root x2+=−C2+Δ/3C3.Case (1):C0>0,Fμ>0.SinceC0>0, we have θ<θ5. Since Fμ>0, we have the following two cases:Case where1/β<δ2/ω2, i.e., ω2<βδ2. Since μ=1/β, we have Fμ=δ3σ1βδ2−ω2/β>0. If Fx2+<0, then there exist E1∗ and E2∗.Case where1/β>δ2/ω2, i.e., ω2>βδ2. Since μ=δ2/ω2, we have Fμ=ασ2βδ2−ω2γ2θ−δ3ρ/ω2>0 implying θ<θ3. Then, we have θ<minθ3,θ5. If Fx2+<0, then there exist E1∗ and E2∗.Case (2):C0<0,Fμ<0.SinceC0<0, we have θ>θ5. Since Fμ>0, we have the following two cases:Case where1/β<δ2/ω2, i.e., ω2<βδ2. Since μ=1/β, we have Fμ=δ3σ1βδ2−ω2/β>0. That is a contradiction.Case where1/β>δ2/ω2, i.e., ω2>βδ2. Since μ=δ2/ω2, we have Fμ=ασ2βδ2−ω2γ2θ−δ3ρ/ω2<0 implying θ>θ3. Then, we have θ>maxθ3,θ5. If Fx1+>0, then there exist E1∗ and E2∗.Next, we show the stability of the above two equilibriaE0 and E∗ of system (34) with v=0. We linearize the system and obtain the characteristic equation whose expression is the same as (13).Theorem 5.
System (34) with v=0 always has one tumor-free equilibrium E0, which is LAS if the inequality θ<θ5 holds.Proof.
AtE0=x0,y0,z0,u0=0,y0,σ2/δ2,γ1δ2y0+γ2σ2/δ2δ3, the characteristic equation becomes(44)λ−α+y0λ+δ2λ2+σ1y0+δ3λ+θγ1y0+σ1δ3y0=0,and one root of characteristic equation is −δ2<0. It is easily noted that as σ1/y0+δ3>0,θγ1y1+σ1δ3/y0>0, so the solutions of λ2+σ1/y0+δ3λ+θγ1y1+σ1δ3/y0=0 have always negative real parts. Therefore, E0 is LAS if and only if α−y0=2αA+B−B2−4AC/2A<0, that is, if and only if θ<θ5.
AtE∗=x∗,y∗,z∗,u∗, characteristic (13) becomes(45)λ4+B1λ3+B2λ2+B3λ+B4=0,where(46)B1=αβx∗+δ3+σ1y∗+σ2z∗>0,B2=αβx∗+σ2z∗σ1y∗+δ3+αβx∗σ2z∗+δ3σ1y∗+θγ1y∗>0,B3=θγ1y∗+δ3σ1y∗αβx∗+σ2z∗+αβx∗σ2z∗δ3+σ1y∗+ρω2x∗y∗z∗>0,B4=αβx∗σ2z∗θγ1y∗+δ3σ1y∗+ρδ3−θγ2ω2x∗y∗z∗.
By the Routh–Hurwitz criterion, the roots of (45) have negative real parts if and only if B1>0,B2>0,B3>0,B4>0,B1B2−B3>0,andB1B2B3−B1B4−B32>0.
Hence, we obtain the sufficient condition for the stability ofP∗:(47)H4:B4=αβx∗σ2z∗θγ1y∗+δ3σ1y∗+ρδ3−θγ2ω2x∗y∗z∗>0,H5:B1B2−B3=K1θ+K2>0,H6:B1B2B3−B1B4−B32=B3K1θ+K2−B12B4=K3θ2+K4θ+K5>0,where(48)K1=σ1+δ3y∗γ1,K2=δ3+σ1y∗+σ2z∗α2β2x∗2+αβx∗+σ2z∗σ1y∗+δ3+αβx∗σ2z∗+δ3σ1y∗σ1y∗+δ3−ρω2x∗y∗z∗,K3=K1αβx∗+σ2z∗γ1y∗,K4=K2γ1y∗+K1σ1δ3y∗αβx∗+σ2z∗+K1αβx∗σ2z∗δ3+σ1y∗−αβx∗+δ3+σ1y∗+σ2z∗2γ1y∗+ρ+γ2ω2x∗y∗z∗,K5=K2αβx∗+σ2z∗σ1δ3y∗−αβx∗σ2z∗K2δ3+σ1y∗−αβx∗+δ3+σ1y∗+σ2z∗2σ1δ3y∗+K2−δ3ρω2x∗y∗z∗.Theorem 6.
System (34) with v=0 has a coexistence equilibrium E∗ when ω2>βδ2,σ1≥ααργ1+δ1γ2/γ2,θ3<θ<θ5 or ω2>βδ2, σ1≤ααργ1+δ1γ2/γ2,θ5<θ<θ3 or ω2<βδ2,θ>θ5 hold. Further E∗ is LAS when condition (47) is satisfied.


## 3.1.2. MAI Model

Whenσ1=0,σ3>0, (36) becomes(49)αx1−βx−xy=0,ρyz−θyu−δ1y=0,σ2+ω2xz−δ2z=0,γ1y+γ2z−δ3u−ξuv=0,σ3−δ4v=0.Equation (49) can be simplified as follows:(50)αx1−βx−xy=0,ρyz−θyu−δ1y=0,σ2+ω2xz−δ2z=0,γ1y+γ2z−δ31u=0,where δ31=δ3δ4+ξσ3/δ4.The analysis of (50) is similar to that of (12), and the equilibria of (34) can be obtained as follows:(51)F0=0,0,σ2δ2,γ2σ2δ2δ31,σ3δ4=0,δ31ρσ2−δ1δ2−θγ2σ2θγ1δ2,σ2δ2,ρσ2−δ1δ2θδ2,σ3δ4,F2=1β,0,βσ2βδ2−ω2,βγ2σ2δ31βδ2−ω2,σ3δ4,F∗=x∗,α1−βx∗,σ2δ2−ω2x∗,αγ11−βx∗δ2−ω2x∗+γ2σ2δ31δ2−ω2x∗,σ3δ4.Here,x∗∈0,μ and satisfies the equation a2x2+a1x+a0=0, where coefficients are defined as a2=−θαβγ1ω2<0,a1=θαγ1βδ2+ω2+δ1δ31ω2>0,anda0=δ31ρσ2−δ1δ2−θαγ1δ2+γ2σ2.The Jacobian matrix of system (34) at any equilibrium F¯x¯,y¯,z¯,u¯,v¯ is as follows:(52)JF¯=α−2αβx¯−y¯−x¯0000ρz¯−θu¯−δ1ρy¯−θy¯0ω2z¯0ω2x¯−δ2000γ1γ2−δ3−ξv¯ξu¯0000−δ4=α−2αβx¯−y¯−x¯0000ρz¯−θu¯−δ1ρy¯−θy¯0ω2z¯0ω2x¯−δ2000γ1γ2−δ31ξu¯0000−δ4.The corresponding characteristic equation is(53)JF¯−λI=0.Substitutingδ31 for δ3 in (12), we have(54)JF¯−λI=λ+δ4JP¯−λI=0.Thus, the stability analysis ofF¯ is similar to that of P¯.Substitutingδ31 for δ3 in (20) to get new conditions for the stability of F∗:(55)h1:A4=θαβγ1δ2−ω2x∗x∗y∗+ρδ31−θγ2ω2x∗y∗z∗>0,h2:A1A2−A3=θγ1αβx∗+δ2−ω2x∗y∗+αβδ31δ2−ω2x∗x∗−ρx∗y∗z∗>0,h3:A1A2A3−A1A4−A32=B1θ+B2>0,where(56)B1=αβx∗+δ2−ω2x∗2δ31+αβx∗+δ2−ω2x∗δ312−ρω2x∗y∗z∗αβx∗+δ2−ω2x∗γ1y∗−αβx∗+δ2−ω2x∗+δ312αβγ1δ2−ω2x∗x∗y∗−ω2γ2x∗y∗z∗,B2=αβx∗+δ2−ω2x∗2δ31+αβx∗+δ2−ω2x∗δ312−ρω2x∗y∗z∗αβδ31δ2−ω2x∗x∗+ρω2x∗y∗z∗−αβx∗+δ2−ω2x∗+δ312ρω2δ31x∗y∗z∗.We set(57)θ70=δ31ρσ2−δ1δ2αγ1δ2+γ2σ2,θ71=δ31ρσ2−δ1δ2γ2σ2,θ72=βδ31ρσ2−δ1δ2+δ1δ31ω2βγ2σ2,θ73=ρδ31γ2,and we obtain the following results.Theorem 7.
System (34) always has one tumor-free and ECs-free equilibrium F0, which is unstable.Theorem 8.
System (34) has one tumor-free equilibrium F1 when ρ>δ1δ2/σ2 and θ<θ71, which is LAS if the inequality θ<θ70 holds.Theorem 9.
System (34) has one tumor-dominant equilibrium F2 when ω2<ω2∗, which is LAS if the inequality θ>θ72 holds.Theorem 10.
System (34) has a unique coexistence equilibrium F∗ when ω2<ω2∗,θ0<θ<θ2 or ω2>ω2∗,θ0<θ<θ3 hold. Further F∗ is LAS, when condition (55) is satisfied.


## 3.1.3. Combined Immunotherapy Model

Whenσ1>0andσ3>0, (36) can be simplified as follows:(58)αx1−βx−xy=0,ρyz−θyu−δ1y+σ1=0,σ2+ω2xz−δ2z=0,γ1y+γ2z−δ31u=0,where δ31=δ3δ4+ξσ3/δ4.The analysis of (58) is similar to that of (37), and the equilibria of (34) can be obtained as follows:(59)N0=x0,y0,z0,u0,v0=0,y0,σ2δ2,γ1δ2y0+γ2σ2δ2δ31,σ3δ4,where y0=−B+B2−4AC/2A with A=θγ1δ2>0,B=δ31δ1δ2−ρσ2+θγ2σ2,andC=−σ1δ2δ31<0.(60)N∗=x∗,y∗,z∗,u∗,v∗=x∗,α1−βx∗,σ2δ2−ω2x∗,αγ11−βx∗δ2−ω2x∗+γ2σ2δ31δ2−ω2x∗,σ3δ4,where x∗∈0,μ and satisfies the equation c3x3+c2x2+c1x+c0=0 with coefficients(61)c3=θα2β2γ1ω2>0,c2=−θα2βγ1ω2−αβθαγ1βδ2+ω2+δ1δ31ω2<0,c1=θα2αβγ1δ2+αγ1ω2+βγ2σ2+αδ31δ1ω2−βρσ2−δ1δ2−σ1ω2δ31,c0=σ1δ2δ31+αδ31ρσ2−δ1δ2−θααγ1δ2+γ2σ2.For system (34) at N0=x0,y0,z0,u0,v0, Jacobian matrix is expressed as follows:(62)JN0=α−y000000−σ1y0ρy0−θy00ω2z00−δ2000γ1γ2−δ31ξu00000−δ4.The corresponding characteristic equation is(63)JN0−λI=λ+δ4JE0−λI=0.Thus, the stability analysis ofN0 is similar to that of E0. Set θ8=αδ31ρσ2−δ1δ2+σ1δ2δ31/ααγ1δ2+σ2γ2, and we can get the following result.Theorem 11.
System (34) always has one tumor-free equilibrium N0, which is LAS if the inequality θ<θ8 holds.
For system (34) at N∗=x∗,y∗,z∗,u∗,v∗, Jacobian matrix is given as follows:(64)JN∗=−αβx∗−x∗0000−σ1y∗ρy∗−θy∗0ω2z∗0ω2x∗−δ2000γ1γ2−δ31ξu∗0000−δ4.
The corresponding characteristic equation is(65)JN∗−λI=λ+δ4JE∗−λI=0.
Thus, the stability analysis ofN∗ is similar to that of E∗. Set θ9=ρδ31/γ2 and substitute δ31 for δ3 in (47) to get new conditions:(66)h4:B4=αβx∗σ2z∗θγ1y∗+δ31σ1y∗+ρδ31−θγ2ω2x∗y∗z∗>0,h5:B1B2−B3=k1θ+k2>0,h6:B1B2B3−B1B4−B32=B3k1θ+k2−B12B4=k3θ2+k4θ+k5>0,where(67)k1=σ1+δ31y∗γ1,k2=δ31+σ1y∗+σ2z∗α2β2x∗2+αβx∗+σ2z∗σ1y∗+δ31+αβx∗σ2z∗+δ31σ1y∗σ1y∗+δ31−ρω2x∗y∗z∗,k3=k1αβx∗+σ2z∗γ1y∗,k4=k2γ1y∗+k1σ1δ31y∗αβx∗+σ2z∗+k1αβx∗σ2z∗δ31+σ1y∗−αβx∗+δ31+σ1y∗+σ2z∗2γ1y∗+ρ+γ2ω2x∗y∗z∗,k5=k2αβx∗+σ2z∗σ1δ31y∗−αβx∗σ2z∗k2δ31+σ1y∗−αβx∗+δ31+σ1y∗+σ2z∗2σ1δ31y∗+k2−δ31ρω2x∗y∗z∗.
Therefore, we can obtain the following result.Theorem 12.
System (34) has a coexistence equilibrium N∗ when ω2>βδ2,σ1≥ααργ1+δ1γ2/γ2,θ9<θ<θ8 or ω2>βδ2, σ1≤ααργ1+δ1γ2/γ2,θ8<θ<θ9 or ω2<βδ2,θ>θ8 hold. Further N∗ is LAS when condition (66) is satisfied.


## 3.2. Numerical Simulations

### 3.2.1. ACI Model (σ1>0,σ3=0)

We conduct numerical simulations of the ACI model. We chooseω2=0.0004 and the other parameter values are given in (22). We find the stability region of E∗=x∗,y∗,z∗,u∗, as shown in Figure 5(a). E∗ is stable in region I and unstable in region II.Case (a): we choose a pointG1=σ1,θ=0.4,0.5 in the region I, then system (34) has one coexistence equilibrium:(68)E∗=11.560,1.598,7.543,6.993.The eigenvalues of characteristic equation of (47) are(69)−0.269,−0.041,−0.1388+0.334iand−0.1388−0.334i,soE∗ is stable, as shown in Figure 6(b).Case (b): we choose a pointG2=σ1,θ=0.4,0.55 in the region II; then, system (34) has one coexistence equilibrium:(70)E∗=69.131,1.409,13.895,11.962.The eigenvalues of characteristic equation of (47) are(71)−0.830,−0.021,0.0326+0.6039iand0.0326−0.6039i,soE∗ is unstable, as shown in Figures 6(c) and 6(d).Figure 5
(a) Stability region ofE∗ and Hopf bifurcation curve in θ−σ1 parameter plane. E∗ is stable in region I and unstable in region II. (b) ACI curve. E0 is stable in region III and unstable in region IV. (c) The bifurcation diagram of x with respect to θ for fixed σ1=0.4. Hopf bifurcation point appears at θACI∗=0.5354. (d) The bifurcation diagram of x with respect to θ for different σ1=0,0.1,0.3,0.5,0.7,0.9 (the blue curves represent the stable steady states of system (34), while the red curves represent the unstable steady states).
(a)(b)(c)(d)Figure 6
Chooseσ1=0.4, and the other parameter values are given in (22). (a) When θ=0.48, E00,1.782,6.909,6.596 is LAS. (b) When θ=0.50, E∗11.560,1.598,7.543,6.993 is LAS. (c) When θ=0.55, E∗69.131,1.409,13.895,11.962 is unstable. (d) When θ=0.55, the 3D phase portrait depicts TCs, ECs, and Tregs.
(a)(b)(c)(d)According to the stability condition of tumor-free equilibriumE0 in Theorem 5, ACI curve can be obtained as(72)θ=αδ3ρσ2−δ1δ2+σ1δ2δ3ααγ1δ2+γ2σ2=0.09391σ1+0.452.We find the stability region ofE0=0,y0,z0,u0 shown in Figure 5(b). E0 is stable in the region III and unstable in the region IV. With the increase of σ1, region III gradually increases, while region IV gradually decreases.Case (c): we choose a pointG3=σ1,θ=0.4,0.48 in the region III; then, system (34) has a tumor-free equilibrium:(73)E0=0,1.782,6.909,6.596.The eigenvalues of Jacobian matrix ofE0 are(74)−0.055,−0.146,−0.237+0.358iand−0.237−0.358i,soE0 is stable, as shown in Figure 6(a).Case (d): we choose a pointG4=σ1,θ=0.4,0.53 in the region IV; then, system (34) has a tumor-free equilibrium:(75)E0=0,1.141,6.909,6.212.The eigenvalues of Jacobian matrix ofE0 are(76)0.494,−0.055,−0.3+0.297iand−0.3−0.297i,so E0 is unstable.We choose ACI parameterσ1=0.4 to study the relationship between the number of TCs x and the parameter θ (see Figure 5(c)). By some calculations, we have θ3=ρδ3/γ2=0.6,θ5=αδ3ρσ2−δ1δ2+σ1δ2δ3/ααγ1δ2+γ2σ2=0.489. And Hopf bifurcation point appears at θACI∗=0.5354. We can obtain the following result.Proposition 5.
When0<θ<θ5, E0 exists and is LAS, and E∗ is nonexistent. When θ5<θ<θACI∗, E0 exists and is unstable, and E∗ exists and is LAS. When θACI∗<θ<θ3, E0 exists and is unstable, and E∗ exists and is unstable. When θ>θ3, E0 exists and is unstable, and E∗ is nonexistent.
Next, we choose different ACI parametersσ1=0,0.1,0.3,0.5,0.7,0.9, to study the relationship between the number of TCs x and the parameter θ (see Figure 5(d)).
By comparing the curves in Figure5(d), we see that, as σ1 increases gradually, the stable region of tumor-free equilibrium E0 of system (34) gradually increases (the intersection point of the curve and the x-coordinate gradually moves to the right) and the stability region of equilibrium E∗ in system (34) gradually increases (the blue curve gradually moves upward). This indicates that increasing the injection volume of adoptive T cells can not only delay the inhibitory effect of Tregs on tumor immune response but also help the immune system to remove more TCs. It also helps the immune system to control more TCs, keeping them at a stable state.


### 3.2.2. MAI Model (σ1=0,σ3>0)

We conduct numerical simulations of the MAI model. We chooseω2=0.0004,ξ=0.05,andδ4=0.25 and the other parameter values are given in (22). We find the stability region of F∗=x∗,y∗,z∗,u∗, as shown in Figure 7(a). F∗ is stable in region I and unstable in region II.Case (a): we choose a pointM1=σ3,θ=0.4,0.6 in region I; then, system (34) has one coexistence equilibrium:(77)F∗=2.413,1.628,7.032,5.002,1.6.The eigenvalues of Jacobian matrix ofF∗ are(78)−0.25,−0.146+0.336i,−0.146−0.336i,−0.0499+0.033iand−0.0499−0.033i,soF∗ is stable, as shown in Figure 8(b).Case (b): we choose a pointM2=σ3,θ=0.4,0.635 in region II; then, system (34) has one coexistence equilibrium:(79)F∗=27.345,1.546,8.624,5.93,1.6.The eigenvalues of Jacobian matrix ofF∗ are(80)−0.4334,−0.25,−0.0644,0.0172+0.429iand0.0172−0.429i,soF∗ is unstable, as shown in Figures 8(c) and 8(d).According to the stability condition of tumor-free equilibriumF1 in Theorem 7, MAI curve can be obtained as(81)θ=δ3δ4+ξσ3ρσ2−δ1δ2/αγ1δ2+σ2γ2δ4=0.3616σ3+0.452.We find the stability region ofF1, as shown in Figure 7(b).F0 is stable in region III and unstable in region IV. With the increase of σ3, region III gradually increases, while region IV gradually decreases.Case (c): we choose a pointM3=σ3,θ=0.4,0.5 in region III; then, system (34) has a tumor-free equilibrium:(82)F1=0,2.556,6.909,5.349,1.6.The eigenvalues of Jacobian matrix ofF1 are(83)−0.92,−0.25,−0.055,−0.164+0.428iand−0.164−0.428i,soF1 is stable, as shown in Figure 8(a).Case (d): we choose a pointM4=σ3,θ=0.4,0.6 in region IV; then, system (34) has a tumor-free equilibrium:(84)F1=0,1.575,6.909,4.903,1.6.Figure 7
(a) Stability region ofF∗ and Hopf bifurcation curve in θ−σ3 parameter plane. F∗ is stable in region I and unstable in region II. (b) MAI curve. F1 is stable in region III and unstable in region IV. (c) The bifurcation diagram of x with respect to θ for fixed σ3=0.4. Hopf bifurcation point appears at θMAI∗=0.627. (d) The bifurcation diagrams of x with respect to θ for different σ3=0,0.1,0.3,0.5,0.7,0.9 (the blue curves represent the stable steady states of system (34), while the red curves represent the unstable steady states).
(a)(b)(c)(d)Figure 8
Chooseσ3=0.4, and the other parameter values are given in (22). (a) When θ=0.5, F10,2.556,6.909,5.349,1.6 is LAS. (b) When θ=0.6, F∗2.413,1.628,7.032,5.002,1.6 is LAS. (c) When θ=0.635, F∗27.345,1.546,8.624,5.93,1.6 is unstable. (d) When θ=0.635, the 3D phase portrait depicts TCs, ECs, and Treg.
(a)(b)(c)(d)The eigenvalues of Jacobian matrix ofF1 are(85)0.061,−0.25,−0.055,−0.164+0.338iand−0.164−0.338i,so F1 is unstable.We choose MAI parameterσ3=0.4 to study the relationship between the number of TCs x and the parameter θ (see Figure 7(c)). By some calculations, we have θ70=δ31ρσ2−δ1δ2/αγ1δ2+γ2σ2=0.597, θ71=δ31ρσ2−δ1δ2/γ2σ2=0.703, and θ73=ρδ31/γ2=0.792. And Hopf bifurcation point appears at θMAI∗=0.627. We can obtain the following result.Proposition 6.
When0<θ<θ70, F0 exists and is unstable, F1 exists and is LAS, and F∗ exists and is LAS. When θ70<θ<θMAI∗, F0 exists and is unstable, F1 exists and is unstable, and F∗ exists and is unstable. When θMAI∗<θ<θ71, F0 exists and is unstable, F1 exists and is unstable, and F∗ exists and is unstable. When θ71<θ<θ73, F0 exists and is unstable, F1 is nonexistent, and F∗ exists and is unstable. When θ>θ3, F0 exists and is unstable, F1 is nonexistent, and F∗ is nonexistent.
Next, we choose different MAI parametersσ3=0,0.1,0.3,0.5,0.7,0.9, to study the relationship between the number of TCs x and the parameter θ (see Figure 7(d)).
By comparing the curves in Figure7(d), we can find that, with the gradual increase of σ3, the stability region of the tumour-free equilibrium F1 of system (34) gradually increases (the intersection point of curves and x-coordinate gradually moves to the right). This shows that the increase of antibody injection quantity can help to slow down Tregs inhibition of tumor immune responses. The stability region of equilibrium F∗ of system (34) gradually increases (the blue curve gradually moves upward), which means that increasing the amount of antibody injected can help the immune system to control more TCs. By comparing Figures 5 and 7, we find that the effect of MAI is better than that of AIC in delaying the inhibitory effect of Tregs on tumor immune response (at the same injection dose, the intersection point of curves and x-coordinate in Figure 7 moves to the right more widely than that in Figure 5). AIC is more effective than MAI in controlling TCs (the blue curve in Figure 7 moves up more than the blue curves in Figure 5 at the same injection dose).


### 3.2.3. Combined Immunotherapy Model (σ1>0,σ3>0)

We conduct numerical simulations of the combined immunotherapy model. We chooseω2=0.0004,ξ=0.05,andδ4=0.25, and the other parameter values are given in (22). We find the stability region of N∗=x∗,y∗,z∗,u∗,v∗, as shown in Figure 9(a). When σ3 increases gradually, the stability region of N∗ also increases gradually.Figure 9
(a) Choose differentσ3=0.1,0.5,0.8. Stability region of N∗ and Hopf bifurcation curve in θ−σ1 parameter plane. N∗ is stable to the left of the curve and unstable to the right of the curve. (b) Combined immunotherapy surface. N0 is stable below the surface and unstable above the surface. (c) When σ1=0.4andσ3=0.4, the bifurcation diagram of x with respect to θ for system (34). Hopf bifurcation point appears at θCI∗=0.727. (d) Choose different σ1,σ3=0.4,0,0,0.4,0.4,0.4,0.7,0.7. The bifurcation diagram of x with respect to θ (the stable steady state is represented by the blue curve, and the unstable one corresponds to the red curve).
(a)(b)(c)(d)Figure 10
(a) Whenσ1=σ3=0.5,θ=0.4, N00,8.3,6.909,7.505,2 is LAS. (b) When σ1=0.4,σ3=0.5,θ=0.7, N∗11.56,1.59,7.54,4.99,2 is LAS. (c) When σ1=0.4,σ3=0.5,θ=0.78, N∗77.887,1.381,15.936,9.698,2 is unstable. (d) When σ1=0.4,σ3=0.5,θ=0.78, the 3D phase portrait depicts TCs, ECs, and Tregs.
(a)(b)(c)(d)Case (a): we choose a pointL1=θ,σ1=0.7,0.4 on the left of the curve σ3=0.5; then, system (34) has one coexistence equilibrium:(86)N∗=11.56,1.59,7.54,4.99,2.The eigenvalues of Jacobian matrix of (64) are(87)−0.06,−0.176,−0.25,−0.22−0.37iand−0.22+0.37i,soN∗ is stable, as shown in Figure 10(b).Case (b): we choose a pointL2=θ,σ1=0.78,0.4 on the right of the curve σ3=0.5; then, system (34) has one coexistence equilibrium:(88)N∗=77.887,1.381,15.936,9.698,2.The eigenvalues of Jacobian matrix of (64) are(89)−0.906238,−0.0248,−0.25,0.00641−0.66iand0.00641+0.66i,soN∗ is unstable, as shown in Figures 10(c) and 10(d).According to the stability condition of tumor-free equilibriumN0 in Theorem 11, combined immunotherapy surface can be obtained as(90)θ=αδ31ρσ2−δ1δ2+σ1δ2δ31ααγ1δ2+σ2γ2=0.055σ1+0.26471.366σ3+1.707.Therefore, the stability region of tumor-free equilibriumN0 of system (34) is obtained (see Figure 9(b)).Case (c): we choose a pointL3=σ1,σ2,θ=0.5,0.5,0.4 below surface (90); then, system (34) has a tumor-free equilibrium:(91)N0=0,8.3,6.909,7.505,2.The eigenvalues of Jacobian matrix of (64) are(92)−6.664,−0.055,−0.25,−0.204+0.69iand−0.204−0.69i,soN0 is stable, as shown in Figure 10(a).Case (d): we choose a pointL4=σ1,σ2,θ=0.5,0.5,0.8 above surface (90); then, system (34) has a tumor-free equilibrium:(93)N0=0,0.886,6.909,4.32,2.The eigenvalues of Jacobian matrix of (64) are(94)0.75,−0.055,−0.25,−0.453+0.311iand−0.453−0.311i,so N0 is unstable.We choose combined immunotherapy parameterσ1=0.4andσ3=0.4 to study the relationship between the number of TCs x and the parameter θ (see Figure 9(c)). By some calculations, we have θ8=αδ31ρσ2−δ1δ2+σ1δ2δ31/ααγ1δ2+σ2γ2=0.646,andθ9=ρδ31/γ2=0.791. And Hopf bifurcation point appears at θCI∗=0.727. We can obtain the following result.Proposition 7.
When0<θ<θ8, N0 exists and is LAS, and N∗ is nonexistent. When θ8<θ<θCI∗, N0 exists and is unstable, and N∗ exists and is LAS. When θCI∗<θ<θ9, N0 exists and is unstable, and N∗ exists and is unstable. When θ>θ9, N0 exists and is unstable, and N∗ is nonexistent.
For combined immunotherapy, by comparing the curves in Figure9(d), we know that, with the increase of σ1 and σ3, the stability region of the tumour-free equilibrium of N0 of system (34) gradually increases (the intersection point of curves and x-coordinate gradually moves to the right), and the stability region of equilibrium N∗ of system (34) gradually increases (the blue curve gradually moves upward). By comparing Figure 9 with Figures 5 and 7, it can be seen that combined immunotherapy has better effects on delaying the inhibitory effect of Tregs on tumor immune response and helps the immune system to control more TCs than ACI or MAI.


## 3.2.1. ACI Model (σ1>0,σ3=0)

We conduct numerical simulations of the ACI model. We chooseω2=0.0004 and the other parameter values are given in (22). We find the stability region of E∗=x∗,y∗,z∗,u∗, as shown in Figure 5(a). E∗ is stable in region I and unstable in region II.Case (a): we choose a pointG1=σ1,θ=0.4,0.5 in the region I, then system (34) has one coexistence equilibrium:(68)E∗=11.560,1.598,7.543,6.993.The eigenvalues of characteristic equation of (47) are(69)−0.269,−0.041,−0.1388+0.334iand−0.1388−0.334i,soE∗ is stable, as shown in Figure 6(b).Case (b): we choose a pointG2=σ1,θ=0.4,0.55 in the region II; then, system (34) has one coexistence equilibrium:(70)E∗=69.131,1.409,13.895,11.962.The eigenvalues of characteristic equation of (47) are(71)−0.830,−0.021,0.0326+0.6039iand0.0326−0.6039i,soE∗ is unstable, as shown in Figures 6(c) and 6(d).Figure 5
(a) Stability region ofE∗ and Hopf bifurcation curve in θ−σ1 parameter plane. E∗ is stable in region I and unstable in region II. (b) ACI curve. E0 is stable in region III and unstable in region IV. (c) The bifurcation diagram of x with respect to θ for fixed σ1=0.4. Hopf bifurcation point appears at θACI∗=0.5354. (d) The bifurcation diagram of x with respect to θ for different σ1=0,0.1,0.3,0.5,0.7,0.9 (the blue curves represent the stable steady states of system (34), while the red curves represent the unstable steady states).
(a)(b)(c)(d)Figure 6
Chooseσ1=0.4, and the other parameter values are given in (22). (a) When θ=0.48, E00,1.782,6.909,6.596 is LAS. (b) When θ=0.50, E∗11.560,1.598,7.543,6.993 is LAS. (c) When θ=0.55, E∗69.131,1.409,13.895,11.962 is unstable. (d) When θ=0.55, the 3D phase portrait depicts TCs, ECs, and Tregs.
(a)(b)(c)(d)According to the stability condition of tumor-free equilibriumE0 in Theorem 5, ACI curve can be obtained as(72)θ=αδ3ρσ2−δ1δ2+σ1δ2δ3ααγ1δ2+γ2σ2=0.09391σ1+0.452.We find the stability region ofE0=0,y0,z0,u0 shown in Figure 5(b). E0 is stable in the region III and unstable in the region IV. With the increase of σ1, region III gradually increases, while region IV gradually decreases.Case (c): we choose a pointG3=σ1,θ=0.4,0.48 in the region III; then, system (34) has a tumor-free equilibrium:(73)E0=0,1.782,6.909,6.596.The eigenvalues of Jacobian matrix ofE0 are(74)−0.055,−0.146,−0.237+0.358iand−0.237−0.358i,soE0 is stable, as shown in Figure 6(a).Case (d): we choose a pointG4=σ1,θ=0.4,0.53 in the region IV; then, system (34) has a tumor-free equilibrium:(75)E0=0,1.141,6.909,6.212.The eigenvalues of Jacobian matrix ofE0 are(76)0.494,−0.055,−0.3+0.297iand−0.3−0.297i,so E0 is unstable.We choose ACI parameterσ1=0.4 to study the relationship between the number of TCs x and the parameter θ (see Figure 5(c)). By some calculations, we have θ3=ρδ3/γ2=0.6,θ5=αδ3ρσ2−δ1δ2+σ1δ2δ3/ααγ1δ2+γ2σ2=0.489. And Hopf bifurcation point appears at θACI∗=0.5354. We can obtain the following result.Proposition 5.
When0<θ<θ5, E0 exists and is LAS, and E∗ is nonexistent. When θ5<θ<θACI∗, E0 exists and is unstable, and E∗ exists and is LAS. When θACI∗<θ<θ3, E0 exists and is unstable, and E∗ exists and is unstable. When θ>θ3, E0 exists and is unstable, and E∗ is nonexistent.
Next, we choose different ACI parametersσ1=0,0.1,0.3,0.5,0.7,0.9, to study the relationship between the number of TCs x and the parameter θ (see Figure 5(d)).
By comparing the curves in Figure5(d), we see that, as σ1 increases gradually, the stable region of tumor-free equilibrium E0 of system (34) gradually increases (the intersection point of the curve and the x-coordinate gradually moves to the right) and the stability region of equilibrium E∗ in system (34) gradually increases (the blue curve gradually moves upward). This indicates that increasing the injection volume of adoptive T cells can not only delay the inhibitory effect of Tregs on tumor immune response but also help the immune system to remove more TCs. It also helps the immune system to control more TCs, keeping them at a stable state.


## 3.2.2. MAI Model (σ1=0,σ3>0)

We conduct numerical simulations of the MAI model. We chooseω2=0.0004,ξ=0.05,andδ4=0.25 and the other parameter values are given in (22). We find the stability region of F∗=x∗,y∗,z∗,u∗, as shown in Figure 7(a). F∗ is stable in region I and unstable in region II.Case (a): we choose a pointM1=σ3,θ=0.4,0.6 in region I; then, system (34) has one coexistence equilibrium:(77)F∗=2.413,1.628,7.032,5.002,1.6.The eigenvalues of Jacobian matrix ofF∗ are(78)−0.25,−0.146+0.336i,−0.146−0.336i,−0.0499+0.033iand−0.0499−0.033i,soF∗ is stable, as shown in Figure 8(b).Case (b): we choose a pointM2=σ3,θ=0.4,0.635 in region II; then, system (34) has one coexistence equilibrium:(79)F∗=27.345,1.546,8.624,5.93,1.6.The eigenvalues of Jacobian matrix ofF∗ are(80)−0.4334,−0.25,−0.0644,0.0172+0.429iand0.0172−0.429i,soF∗ is unstable, as shown in Figures 8(c) and 8(d).According to the stability condition of tumor-free equilibriumF1 in Theorem 7, MAI curve can be obtained as(81)θ=δ3δ4+ξσ3ρσ2−δ1δ2/αγ1δ2+σ2γ2δ4=0.3616σ3+0.452.We find the stability region ofF1, as shown in Figure 7(b).F0 is stable in region III and unstable in region IV. With the increase of σ3, region III gradually increases, while region IV gradually decreases.Case (c): we choose a pointM3=σ3,θ=0.4,0.5 in region III; then, system (34) has a tumor-free equilibrium:(82)F1=0,2.556,6.909,5.349,1.6.The eigenvalues of Jacobian matrix ofF1 are(83)−0.92,−0.25,−0.055,−0.164+0.428iand−0.164−0.428i,soF1 is stable, as shown in Figure 8(a).Case (d): we choose a pointM4=σ3,θ=0.4,0.6 in region IV; then, system (34) has a tumor-free equilibrium:(84)F1=0,1.575,6.909,4.903,1.6.Figure 7
(a) Stability region ofF∗ and Hopf bifurcation curve in θ−σ3 parameter plane. F∗ is stable in region I and unstable in region II. (b) MAI curve. F1 is stable in region III and unstable in region IV. (c) The bifurcation diagram of x with respect to θ for fixed σ3=0.4. Hopf bifurcation point appears at θMAI∗=0.627. (d) The bifurcation diagrams of x with respect to θ for different σ3=0,0.1,0.3,0.5,0.7,0.9 (the blue curves represent the stable steady states of system (34), while the red curves represent the unstable steady states).
(a)(b)(c)(d)Figure 8
Chooseσ3=0.4, and the other parameter values are given in (22). (a) When θ=0.5, F10,2.556,6.909,5.349,1.6 is LAS. (b) When θ=0.6, F∗2.413,1.628,7.032,5.002,1.6 is LAS. (c) When θ=0.635, F∗27.345,1.546,8.624,5.93,1.6 is unstable. (d) When θ=0.635, the 3D phase portrait depicts TCs, ECs, and Treg.
(a)(b)(c)(d)The eigenvalues of Jacobian matrix ofF1 are(85)0.061,−0.25,−0.055,−0.164+0.338iand−0.164−0.338i,so F1 is unstable.We choose MAI parameterσ3=0.4 to study the relationship between the number of TCs x and the parameter θ (see Figure 7(c)). By some calculations, we have θ70=δ31ρσ2−δ1δ2/αγ1δ2+γ2σ2=0.597, θ71=δ31ρσ2−δ1δ2/γ2σ2=0.703, and θ73=ρδ31/γ2=0.792. And Hopf bifurcation point appears at θMAI∗=0.627. We can obtain the following result.Proposition 6.
When0<θ<θ70, F0 exists and is unstable, F1 exists and is LAS, and F∗ exists and is LAS. When θ70<θ<θMAI∗, F0 exists and is unstable, F1 exists and is unstable, and F∗ exists and is unstable. When θMAI∗<θ<θ71, F0 exists and is unstable, F1 exists and is unstable, and F∗ exists and is unstable. When θ71<θ<θ73, F0 exists and is unstable, F1 is nonexistent, and F∗ exists and is unstable. When θ>θ3, F0 exists and is unstable, F1 is nonexistent, and F∗ is nonexistent.
Next, we choose different MAI parametersσ3=0,0.1,0.3,0.5,0.7,0.9, to study the relationship between the number of TCs x and the parameter θ (see Figure 7(d)).
By comparing the curves in Figure7(d), we can find that, with the gradual increase of σ3, the stability region of the tumour-free equilibrium F1 of system (34) gradually increases (the intersection point of curves and x-coordinate gradually moves to the right). This shows that the increase of antibody injection quantity can help to slow down Tregs inhibition of tumor immune responses. The stability region of equilibrium F∗ of system (34) gradually increases (the blue curve gradually moves upward), which means that increasing the amount of antibody injected can help the immune system to control more TCs. By comparing Figures 5 and 7, we find that the effect of MAI is better than that of AIC in delaying the inhibitory effect of Tregs on tumor immune response (at the same injection dose, the intersection point of curves and x-coordinate in Figure 7 moves to the right more widely than that in Figure 5). AIC is more effective than MAI in controlling TCs (the blue curve in Figure 7 moves up more than the blue curves in Figure 5 at the same injection dose).


## 3.2.3. Combined Immunotherapy Model (σ1>0,σ3>0)

We conduct numerical simulations of the combined immunotherapy model. We chooseω2=0.0004,ξ=0.05,andδ4=0.25, and the other parameter values are given in (22). We find the stability region of N∗=x∗,y∗,z∗,u∗,v∗, as shown in Figure 9(a). When σ3 increases gradually, the stability region of N∗ also increases gradually.Figure 9
(a) Choose differentσ3=0.1,0.5,0.8. Stability region of N∗ and Hopf bifurcation curve in θ−σ1 parameter plane. N∗ is stable to the left of the curve and unstable to the right of the curve. (b) Combined immunotherapy surface. N0 is stable below the surface and unstable above the surface. (c) When σ1=0.4andσ3=0.4, the bifurcation diagram of x with respect to θ for system (34). Hopf bifurcation point appears at θCI∗=0.727. (d) Choose different σ1,σ3=0.4,0,0,0.4,0.4,0.4,0.7,0.7. The bifurcation diagram of x with respect to θ (the stable steady state is represented by the blue curve, and the unstable one corresponds to the red curve).
(a)(b)(c)(d)Figure 10
(a) Whenσ1=σ3=0.5,θ=0.4, N00,8.3,6.909,7.505,2 is LAS. (b) When σ1=0.4,σ3=0.5,θ=0.7, N∗11.56,1.59,7.54,4.99,2 is LAS. (c) When σ1=0.4,σ3=0.5,θ=0.78, N∗77.887,1.381,15.936,9.698,2 is unstable. (d) When σ1=0.4,σ3=0.5,θ=0.78, the 3D phase portrait depicts TCs, ECs, and Tregs.
(a)(b)(c)(d)Case (a): we choose a pointL1=θ,σ1=0.7,0.4 on the left of the curve σ3=0.5; then, system (34) has one coexistence equilibrium:(86)N∗=11.56,1.59,7.54,4.99,2.The eigenvalues of Jacobian matrix of (64) are(87)−0.06,−0.176,−0.25,−0.22−0.37iand−0.22+0.37i,soN∗ is stable, as shown in Figure 10(b).Case (b): we choose a pointL2=θ,σ1=0.78,0.4 on the right of the curve σ3=0.5; then, system (34) has one coexistence equilibrium:(88)N∗=77.887,1.381,15.936,9.698,2.The eigenvalues of Jacobian matrix of (64) are(89)−0.906238,−0.0248,−0.25,0.00641−0.66iand0.00641+0.66i,soN∗ is unstable, as shown in Figures 10(c) and 10(d).According to the stability condition of tumor-free equilibriumN0 in Theorem 11, combined immunotherapy surface can be obtained as(90)θ=αδ31ρσ2−δ1δ2+σ1δ2δ31ααγ1δ2+σ2γ2=0.055σ1+0.26471.366σ3+1.707.Therefore, the stability region of tumor-free equilibriumN0 of system (34) is obtained (see Figure 9(b)).Case (c): we choose a pointL3=σ1,σ2,θ=0.5,0.5,0.4 below surface (90); then, system (34) has a tumor-free equilibrium:(91)N0=0,8.3,6.909,7.505,2.The eigenvalues of Jacobian matrix of (64) are(92)−6.664,−0.055,−0.25,−0.204+0.69iand−0.204−0.69i,soN0 is stable, as shown in Figure 10(a).Case (d): we choose a pointL4=σ1,σ2,θ=0.5,0.5,0.8 above surface (90); then, system (34) has a tumor-free equilibrium:(93)N0=0,0.886,6.909,4.32,2.The eigenvalues of Jacobian matrix of (64) are(94)0.75,−0.055,−0.25,−0.453+0.311iand−0.453−0.311i,so N0 is unstable.We choose combined immunotherapy parameterσ1=0.4andσ3=0.4 to study the relationship between the number of TCs x and the parameter θ (see Figure 9(c)). By some calculations, we have θ8=αδ31ρσ2−δ1δ2+σ1δ2δ31/ααγ1δ2+σ2γ2=0.646,andθ9=ρδ31/γ2=0.791. And Hopf bifurcation point appears at θCI∗=0.727. We can obtain the following result.Proposition 7.
When0<θ<θ8, N0 exists and is LAS, and N∗ is nonexistent. When θ8<θ<θCI∗, N0 exists and is unstable, and N∗ exists and is LAS. When θCI∗<θ<θ9, N0 exists and is unstable, and N∗ exists and is unstable. When θ>θ9, N0 exists and is unstable, and N∗ is nonexistent.
For combined immunotherapy, by comparing the curves in Figure9(d), we know that, with the increase of σ1 and σ3, the stability region of the tumour-free equilibrium of N0 of system (34) gradually increases (the intersection point of curves and x-coordinate gradually moves to the right), and the stability region of equilibrium N∗ of system (34) gradually increases (the blue curve gradually moves upward). By comparing Figure 9 with Figures 5 and 7, it can be seen that combined immunotherapy has better effects on delaying the inhibitory effect of Tregs on tumor immune response and helps the immune system to control more TCs than ACI or MAI.


## 4. Discussion and Conclusion

Tregs-mediated tumor immune escape is one of the core mechanisms of tumor immune regulation. And Tregs have been found to mediate tumor evasion and immune escape in many different solid tumors [46]. The study on Tregs has a very high research value and application prospect in the immunotherapy of tumors. If the activity of Tregs is controlled or blocked during the tumor immune response, or a barrier is set to prevent Tregs from migrating into the tumor microenvironment, then the effect of tumor immunotherapy can be improved [47].First, we developed a mathematical model to study the inhibitory role of Tregs in the tumor immune system. For the lower recognition of tumor antigens by the immune system, the stronger the inhibition effect of Tregs on ECs, TCs can easily escape the control of the immune system (see Figure4(a)). When the immune system is highly sensitive to tumor antigens, the immune system activates ECs; the stronger the inhibition effect of Tregs on ECs, the more complicated interactions between TCs and immune cells (see Figure 4(b)).Second, we incorporated the previous mathematical model with three types of immunotherapy to obtain ACI model, MAI model, and combined both ACI and MAI model. Through the theoretical analysis and numerical simulations, we found that ACI can control more TCs, but have no obvious effects on reducing the inhibitory effect of Tregs on ECs (see Figure5). MAI can effectively reduce the inhibitory effect of Tregs on ECs, but cannot control more TCs (see Figure 7). However, combination immunotherapy with ACI and MAI is more effective than single immunotherapy. It can not only significantly reduce the inhibitory effect of Tregs on ECs but also help the immune system to kill TCs to the maximum extent (see Figure 9). Therefore, we recommend the use of combined immunotherapy in the treatment of tumors. Besides, clinical trials are needed to further evaluate the safety and efficacy of combined immunotherapy.This paper focused on the general process of tumor immune response with negative feedback. Using the mathematical model, it is possible to simulate the state of tumors in the immune system at different inhibition states. The results of the study can contribute to the understanding of tumor immunity; at the same time, it also provides new ideas for the treatment of tumors. However, due to the complexity and heterogeneity of tumor microenvironment, there is still a certain gap between the mathematical model and the description of the real interactions between the tumor and immune system. Therefore, specific tumor microenvironment and heterogeneous tumoral populations should be considered in practical application to make the model more realistic [34]. In addition to the bilinear incidence model considered in this paper, nonbilinear model with saturation incidence should be employed in the further study [35]. Besides, other important factors such as immune activation delays [36–38, 48], stochastic effects [39, 40, 49], and impulsive perturbations [41, 42, 50] can be considered in the modelling of the tumor immune system. The tumor immune response dynamics in vivo is very complex and not well understood primarily because the measurements of the necessary parameters are difficult in vivo [43].



---

*Source: 4834165-2020-08-12.xml*