Wentao Gao
College of Mathematics and Statistics, Northwest Normal University, Lanzhou, China.
DOI: 10.4236/jamp.2026.142042   PDF    HTML   XML   113 Downloads   368 Views   Citations

Abstract

In this paper, we discuss a class of dynamic models of the interaction between tumors and the immune system with antigen delay and Michaelis-Menten type inhibition terms. The Michaelis-Menten type function $\frac{\beta E\left(t\right)T\left(t\right)}{a+T\left(t\right)}$ and $\frac{\alpha E\left(t\right)T\left(t\right)}{a+T\left(t\right)}$ is used to describe the immune response of effector cells interacting with tumor cells, the linear antigen stimulation term $cT$ is used to describe the linear recruitment effect of tumor antigens on effector cells, and the stimulation delay of tumor antigens in the immune system is introduced. The dynamic behavior of the model is studied through qualitative analysis and numerical simulation. Saddle-node bifurcation may occur both in the case with and without time delay. Contrary to the case without time delay, stimulation delay may lead to some complex dynamic behaviors and biological phenomena. In the presence of time delay, the existence condition of Hopf bifurcation at the equilibrium point is obtained. Further discussion shows that the model may exhibit bistability under certain conditions, that is, the growth and development state of the tumor depends on its initial state. Finally, numerical simulation is used to verify the accuracy of the relevant theoretical results, and the corresponding biological significance is briefly discussed.

Keywords

Tumor Immunity, Antigen, Michaelis-Menten Type, Equilibrium Point, Time Delay, Stability, Hopf Bifurcation

Share and Cite:

1. Introduction

Tumors are benign or malignant abnormally growing neoplastic tissues with no physiological function in the human body. Malignant tumors (i.e., cancers) are usually caused by the uncontrolled rapid proliferation of cells, which affects the quality of life of patients and threatens their health. The immune system is one of the systems in the human body that defends against the invasion of foreign pathogenic microorganisms, and it has the functions of immune surveillance, defense and regulation. Studies have found that the surface antigens of tumor cells may be the same as or different from those of normal cells, so immune responses may occur. On the other hand, the immune system is a very complex network composed of immune organs, immune cells and immune active substances. The immune system can promote and inhibit the development of tumor cells. Today, we recognize that the immune system plays a dual role in cancer: it can exert anti-tumor effects by destroying cells or inhibiting their growth, and can also promote tumor progression by selecting cancer cells more adaptable to the immune-active host environment or creating conditions conducive to tumor growth in the tumor microenvironment. Due to the complexity of the interaction between tumors and the immune system, the relevant mechanisms are not fully understood, and there are still some issues to be further explored [1] [2]. Mathematical modeling and analysis have played a significant role in this regard [3]-[7].

In the process of tumor induction and growth, the complex interaction between tumor cells and effector cells is determined by three key factors: the malignant potential of the tumor, the antigenicity of the tumor, and the immune response of the host [8] [9]. The malignant potential of a tumor refers to its ability to metastasize, escape and destroy the immune system; the antigenicity of a tumor is defined as the initial size of the effector cell population that can be stimulated after the introduction of antigens, which is related to the tumor and varies significantly among different patients and cancer types. A larger value indicates that the antigens presented by tumor cells are more easily recognized, and a smaller value indicates weaker antigenicity; the immune response reflects the inhibitory effect of the host’s immune system on tumor growth, that is, the proliferation process of effector cells in the presence of tumor cells, and this proliferation process depends on the antigenicity of the tumor [10]-[12]. The classic tumor-immune interaction model was first proposed by Kuznetsov et al. [11] and Kirschner and Panetta [12]. In reference [11], Kuznetsov et al. considered the interaction between tumor cells and effector cells and established the following model:

$\{\begin{array}{l}\frac{\text{d}E}{\text{d}t}=s+\frac{cET}{\epsilon +T}-\beta ET-\mu E,\hfill \\ \frac{\text{d}T}{\text{d}t}=rT\left(1-\frac{T}{K}\right)-\alpha ET.\hfill \end{array}$ (1.1)

where $E=E\left(t\right)$ and $T=T\left(t\right)$ represent the number of effector cells and tumor cells at time $t$ , respectively. $s$ is the normal flow rate of effector cells from external sources to the tumor site (independent of the presence of tumor cells); the saturation term $\frac{cET}{\epsilon +T}$ (with saturation effect) describes the process of tumor cells producing effector cells through antigen stimulation; $\beta$ is the rate at which tumor cells kill or inactivate effector cells; $\mu$ is the natural mortality rate of effector cells. It is assumed that in the absence of effector cells, tumor cells follow the logistic growth model with an intrinsic growth rate $r$ and a carrying capacity $K$ ; effector cells kill tumor cells at a rate $\alpha$ through the law of mass action.

In addition to the tumor-free equilibrium point, the model can have at most three positive equilibrium points (tumor-present equilibrium points). Phase diagram analysis shows that the model exhibits various phenomena, including immune stimulation of tumor growth, tumor escape, and formation of tumor dormant state. However, this model does not study complex dynamic behaviors such as periodic solutions and Hopf bifurcation. Kirschner and Panetta [12] established another theoretical model to explore the effect of effector cells on tumor growth and regression by studying the role of the cytokine interleukin-2 (IL-2) in a single tumor site. It is assumed that IL-2 is mainly produced by effector cells and can induce the generation of corresponding effector cells. The model is as follows:

$\{\begin{array}{l}\frac{\text{d}E}{\text{d}t}=s_1+cT-\frac{{\beta }_{1}EI}{{\epsilon }_1+I}-{\mu }_{1}E,\hfill \\ \frac{\text{d}T}{\text{d}t}=rT\left(1-\frac{T}{K}\right)-\frac{\alpha ET}{\epsilon +T},\hfill \\ \frac{\text{d}I}{\text{d}t}=s_2+\frac{{\beta }_{2}EI}{{\epsilon }_2+I}-{\mu }_{2}I,\hfill \end{array}$ (1.2)

where the biological meanings of $E\left(t\right)$ and $T\left(t\right)$ are consistent with those in model (1.1), $I=I\left(t\right)$ represents the concentration of the cytokine interleukin-2 (IL-2) at time $t$ , and the specific biological meanings of the parameters refer to reference [12]. Kirschner and Panetta regarded parameters $s_1$ and $s_2$ as therapeutic terms, mainly to explore the effect of adoptive cell immunotherapy (ACI). This model can also have at most three tumor-present equilibrium points, and Hopf bifurcation may occur to produce periodic solutions within a specific parameter range. However, this study only conducted numerical analysis for some appropriate parameter values and did not obtain general theoretical results. In addition, unlike model (1.1), model (1.2) is a three-dimensional system, which is difficult to intuitively show the interaction effect between tumors and the immune system and related clinical phenomena. Since then, combining the modeling ideas and assumptions of references [11] [12], researchers have introduced time delays, random terms and detailed immune mechanisms, proposed more mathematical models of tumor-immune systems, and discussed the dynamic properties of these models to explain other clinical phenomena.

Delisi, Adam [13] [14] and others studied ordinary differential equation models of the interaction between tumor cells and effector cells. Both models used Michaelis-Menten type inhibition functions to represent the inhibitory effect of effector cells on tumor cells. The research results showed that within a certain range, the growth of effector cells will increase the survival rate of tumor cells. In addition, they also gave the threshold condition for tumor growth to be uncontrolled by the immune system and become malignant tumors. On the basis of references [11] [13], Kirschner et al. [12] assumed that the stimulation term of tumors on effector cells is a linear term and first proposed a mathematical model of the interaction between tumor cells, immune effector cells and IL-2. The study found that the antigenicity of tumors has a very significant impact on the dynamic properties of the model and explained the reasons for tumor recurrence. On the basis of reference [12], Yang [15] and others added pulsed immunotherapy to establish a mathematical model related to tumor-immunity and pulsed immunotherapy. The study found that the initial density of effector cells, the ratio of effector cells to tumor cells, and the cycle of immunotherapy are crucial for cancer treatment. Zhang et al. used the Michaelis-Menten form to represent the inhibitory effect of tumors on effector cells in reference [16], conducted dynamic analysis on it, and further considered the impact of the inhibition rate coefficient of tumor cells on effector cells on the dynamic behavior of the model. Galach replaced the saturated form stimulation growth rate of effector cells in model (1.1) with the bilinear form $\theta ET$ in reference [17]. This model only includes two variables of tumors and effector cells, and conducts local dynamic analysis on it. Li et al. added an antigen term in reference [18], proposed and studied a simple model of tumor-immune interaction. Through qualitative and quantitative analysis, the model has complex dynamic behaviors. The models of all assume that tumor growth follows a logistic model in the absence of immune system effects, which results in the tumor cell count eventually being bounded. The obtained results can explain some biological phenomena.

As we all know, time delay plays an important role in describing the interaction between tumors and the immune system. Delays may be due to the time lag caused by tumor cell proliferation [19] [20], the growth process of effector cells stimulated by tumor cells [21]-[24], In reference [22], the local stability of the model is obtained by analyzing the characteristic equations of the model at the corresponding equilibria, the sufficient conditions on the global stability are found by applying the Fluctuation Lemma and constructing the different convergent sequences. The obtained results show that, compared to the results for the model without time delay, the time delay of tumor action can affect the stability of tumor equilibrium of the model as the stimulation effect of the tumor cells is strong enough, while the delay is harmless for the stability of tumor equilibrium under the neutralization of tumor cells. For the appropriate neutralization of tumor cells on effector cells, the bistability of the tumor free equilibrium and the stronger tumor equilibrium can appear. In the case of stimulation of tumor cells, the sufficiently large time delay can lead to the appearance of a stable periodic solution by Hopf bifurcation. The differentiation of effector cells, and the neutralization process of effector cells by tumor cells [24]. When there is only the neutralization delay, the model has a uniform upper bound while when there is only the stimulation delay, the bound varies with the delay. The paper [25] presents three quantities with clear biological significance to determine the asymptotic states of tumor progression, while also analyzing the differences in asymptotic states under two ways of describing anti-tumor immunity. the model exhibits rich dynamical behaviors including super-critical and sub-critical Bogdanov-Takens bifurcations (consisting of Hopf bifurcation, saddle-node bifurcation, and homoclinic bifurcation) and saddle-node bifurcation of nonconstant periodic solutions (leading to the appearance of two periodic orbits) as the parameters vary; A large number of studies have been conducted on dynamic models of tumor-immune systems with delays in the literature [26]-[28]. In particular, the research on models described by two-dimensional delay differential equations has obtained rich theoretical results, including the oscillation of solutions, the existence of periodic solutions, Hopf bifurcation, chaos, etc.

In this paper, our purpose is to explore the effect of tumor antigen stimulation delay by qualitatively studying a delayed tumor-immune system model containing antigens and Michaelis-Menten type inhibition functions. We focus on analyzing the effect of tumor cells on the immune system (including stimulation and inhibition), and illustrate the existence of Hopf bifurcation and saddle-node bifurcation. In addition, numerical simulation reveals the impact of antigen delay on the asymptotic state of tumor development and proves the complexity of the dependence of the asymptotic state on initial conditions for different delay values.

The rest of the paper is organized as follows. In the next section, we will establish a two-dimensional model of the interaction between tumors and the immune system. Then, in Section 3, we obtain the non-negativity and boundedness of the solutions of the model. In Section 4, the existence of equilibrium points of the model is obtained. We study the global behavior of the model without delay in Section 5. Next is the analysis of the model with delay in Section 6, where numerical simulation shows the complexity of dynamic behavior. Finally, the paper briefly summarizes and discusses the impact of delay.

2. Model Construction

As mentioned in the previous section, regarding the interaction process between tumors and the immune system, it is sufficient to consider two variables: tumor cells (T) and effector cells (E). The following is the form of the two-dimensional model of the tumor-immune system [11] [18] [22] [24] [25] [29]:

$\{\begin{array}{l}\frac{\text{d}E}{\text{d}t}=s+{\Phi }_1-{\Phi }_2-\mu E,\hfill \\ \frac{\text{d}T}{\text{d}t}=rT\left(1-\frac{T}{K}\right)-{\Phi }_3,\hfill \end{array}$ (2.1)

where it is assumed that in the absence of interaction, the growth of effector cells and tumor cells follows the equations $E^{\prime }=s-\mu E$ and $T^{\prime }=r_{1}T\left(1-\frac{T}{K}\right)$ , respectively. Here, ${\Phi }_1$ represents the effector cells recruited due to the stimulation of tumor antigens, ${\Phi }_2$ describes the anti-immunity of tumors, and ${\Phi }_3$ reflects the killing and destruction of tumors by the immune system. The expressions of ${\Phi }_i\left(i=1,2,3\right)$ are listed in Table 1, where $t$ and ${\tau }_i\left(i=1,2\right)$ are the corresponding delays, and the biological meanings of other parameters can be found in the above references.

Table 1. Expressions of the interaction between tumor cells and effector cells.

${\Phi }_1$	${\Phi }_2$	${\Phi }_3$	References
$\frac{\sigma E\left(t\right)T\left(t\right)}{1+\epsilon T\left(t\right)}$	$\beta E\left(t\right)T\left(t\right)$	$\alpha E\left(t\right)T\left(t\right)$	[11]
$\sigma T\left(t\right)$	$\beta E\left(t\right)T\left(t\right)$	$\frac{\alpha E\left(t\right)T\left(t\right)}{1+\epsilon T\left(t\right)}$	[18]
$\sigma E\left(t\right)T\left(t-\tau \right)$	$\beta E\left(t\right)T\left(t-\tau \right)$	$\alpha E\left(t\right)T\left(t\right)$	[22]
$\sigma E\left(t\right)T\left(t-{\tau }_1\right)$	$\beta E\left(t\right)T\left(t-{\tau }_2\right)$	$\alpha E\left(t\right)T\left(t\right)$	[24]
$\sigma T\left(t-\tau \right)$	$\beta E\left(t\right)T\left(t\right)$	$\alpha E\left(t\right)T\left(t\right)$	[25]
$\frac{\sigma E\left(t-\tau \right)T\left(t-\tau \right)}{1+\epsilon T\left(t-\tau \right)}$	$\beta E\left(t\right)T\left(t\right)$	$\alpha E\left(t\right)T\left(t\right)$	[29]

Compared with existing tumor immunology models (such as those in Table 1), which mainly focus on instantaneous forms of immune suppression or single nonlinear stimulation terms, there are the following gaps: they do not simultaneously integrate the dual-module effects of antigen-delayed stimulation and saturation inhibition;

In this paper, we will use the Michaelis-Menten type inhibition function to describe the immune response of the interaction between effector cells and tumor cells. The increase in effector cells caused by tumor antigenicity is proportional to the concentration of tumor cells. Considering

${\Phi }_1=cT\left(t-{\tau }_1\right),{\Phi }_2=\frac{\beta E\left(t\right)T\left(t\right)}{a+T\left(t\right)},{\Phi }_3=\frac{\alpha E\left(t\right)T\left(t\right)}{a+T\left(t\right)}$

in model (2.1), where ${\tau }_1$ is the stimulation delay of tumor antigens, we get the following model:

$\{\begin{array}{l}\frac{\text{d}E\left(t\right)}{\text{d}t}=s+cT\left(t-{\tau }_1\right)-\frac{\beta E\left(t\right)T\left(t\right)}{a+T\left(t\right)}-\mu E\left(t\right),\hfill \\ \frac{\text{d}T\left(t\right)}{\text{d}t}=rT\left(t\right)\left(1-\frac{T\left(t\right)}{K}\right)-\frac{\alpha E\left(t\right)T\left(t\right)}{a+T\left(t\right)}.\hfill \end{array}$ (2.2)

In [23], the authors focused on the impact of tumor anti-immunity (i.e., $\beta$ ) and immune anti-tumor (i.e., $\alpha$ ) on the model, and discussed the impact of delay on dynamic behavior. Model (2.2) also allows at most two tumor-present equilibrium points, and there are no periodic solutions in the absence of delay. The addition of stimulation delay can also lead to the occurrence of some bifurcations, such as Hopf bifurcation and saddle-node bifurcation. In this paper, for model (2.2), we focus on the impact of tumors on the immune system (reflected by parameters $\eta$ and $\delta$ ) and antigen delay $T$ . Through theoretical analysis, the conditions determining the stability of the tumor-present equilibrium are clearly expressed by the relationship between $\eta$ and $\delta$ . Numerical simulation shows the complexity of the dependence of tumor development state on initial conditions.

For the convenience of mathematical analysis, we perform dimensionless transformation on Equation (2.2):

Let

$\overline{t}=\mu t,\overline{x}=\frac{\mu }{s}E,\overline{y}=\frac{T}{K},$

and denote

$\eta =\frac{cK}{s},\delta =\frac{\beta K}{\mu },q=\frac{\alpha s}{{\mu }^2},\rho =\frac{r}{\mu },h=\frac{a}{K},\tau =\mu {\tau }_1.$

At the same time, we still write $\overline{t}$ as $t$ , then model (2.2) becomes:

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=1+\eta y\left(t-\tau \right)-\frac{\delta xy}{h+y}-x,\hfill \\ \frac{\text{d}y}{\text{d}t}=\rho y\left(1-y\right)-\frac{qxy}{h+y}.\hfill \end{array}$ (2.3)

3. Non-Negativity and Boundedness of Model Solutions

Model (2.3) must be studied under the following initial conditions. Defined in space:

$\phi =\left({\phi }_1,{\phi }_2\right),C_{+}=\left\{\phi \in C\left(\left[-\tau ,0\right],{ℝ}_{+}^2\right):x\left(\theta \right)={\phi }_1\left(\theta \right),y\left(\theta \right)={\phi }_2\left(\theta \right)\right\},$ (3.1)

where ${\phi }_i\left(\theta \right)\ge 0$ , $\theta \in \left[-\tau ,0\right]$ , ${\phi }_i>0$ , $i=1,2$ .$\phi =\left({\phi }_1,{\phi }_2\right)\in C\left(\left[-\tau ,0\right],{ℝ}_{+}^2\right)$ , where $C$ denotes the Banach space of continuous functions from $\left[-\tau ,0\right]\left(\tau >0\right)$ to ${ℝ}_{+}^2$ with the norm:

$\Vert \phi \Vert =\text{max}\left\{\underset{-\tau \le \theta \le 0}{\text{sup}}\left|{\phi }_1\left(\theta \right)\right|,\underset{-\tau \le \theta \le 0}{\text{sup}}\left|{\phi }_2\left(\theta \right)\right|\right\},\phi \in C_{+}.$

${ℝ}_{+}^2=\left\{\left(x,y\right)|x\ge 0,y\ge 0\right\}.$

Theorem 1. Any solution of model (2.3) under initial condition (3.1) is defined on $\left[0,+\infty \right)$ and is positive for all $t>0$ .

Proof. Let $\left(x\left(t\right),y\left(t\right)\right)$ be any solution of model (2.3) under initial condition (3.1).

From the first equation of model (2.3), we have:

${\frac{\text{d}x}{\text{d}t}|}_{x=0}=1+\eta y\left(t-\tau \right)>0.$

The solution is:

$x\left(t\right)=x\left(0\right){\text{e}}^{-{\displaystyle {\int }_0^t\left(1+\frac{\delta y\left(s\right)}{h+y\left(s\right)}\right)\text{d}s}}+{\displaystyle {\int }_0^t\left(1+\eta y\left(s-\tau \right)\right){\text{e}}^{{\displaystyle {\int }_s^t\left(1+\frac{\delta y\left(\xi \right)}{h+y\left(\xi \right)}\right)\text{d}\xi }}\text{d}s}>0.$

Since $x\left(t\right)>0$ , the solution of the equation $\frac{\text{d}y}{\text{d}t}=\rho y\left(1-y\right)-\frac{qxy}{h+y}$ is:

$y\left(t\right)=y\left(0\right){\text{e}}^{{\displaystyle {\int }_0^t\left(\rho \left(1-s\right)-\frac{qs}{h+s}\right)\text{d}s}}>0.$

because $y\left(0\right)>0$ and the exponential function is always positive.

Theorem 2. There exists a constant $M>0$ such that all solutions $\left(x\left(t\right),y\left(t\right)\right)$ of model (2.3) satisfy:

$\underset{t\to \infty }{\lim \sup }x\left(t\right)\le M,\underset{t\to \infty }{\lim \sup }y\left(t\right)\le M.$

Proof. From the equation of $y\left(t\right)$ :

$\frac{\text{d}y}{\text{d}t}=\rho y\left(1-y\right)-\frac{qxy}{h+y}\le \rho y\left(1-y\right).$

Thus, we have:

$y\left(t\right)\le \frac{y\left(0\right)}{y\left(0\right)+\left(1-y\left(0\right)\right){\text{e}}^{-\rho t}},$

so,

$\underset{t\to +\infty }{\lim \sup }y\left(t\right)\le 1.$

Since $y\left(t\right)\le 1$ , for any $\epsilon >0$ , there exists $t_1>0$ such that $y\left(t\right)<1+\epsilon$ when $t>t_1$ . Then, when $t>t_1$ , from the first equation of system (2.3), we have:

$\frac{\text{d}x}{\text{d}t}=1+\eta y\left(t-\tau \right)-\frac{qxy}{h+y}\le 1+\eta y\left(t-\tau \right)\le 1+\eta \left(1+\epsilon \right),$

Therefore, $\underset{t\to +\infty }{\lim \sup }x\left(t\right)\le 1+\eta \left(1+\epsilon \right)$ . Due to the arbitrariness of $\epsilon$ , we get $\underset{t\to +\infty }{\lim \sup }x\left(t\right)\le 1+\eta$ . Thus, the system (2.3) has a positive invariant set:

$D=\left\{\left(x,y\right)\right|x\in \left(0,1+\eta \left],y\in \right[0,1\right]\},$

that is, the system (2.3) is bounded. The dynamic properties of system (2.3) will be studied on $D$ below.

4. Existence of Equilibrium Points

When $\tau =0$ , i.e., without considering stimulation delay, the DDE system (2.3) becomes the following ODE system:

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=1+\eta y-\frac{\delta xy}{h+y}-x,\hfill \\ \frac{\text{d}y}{\text{d}t}=\rho y\left(1-y\right)-\frac{qxy}{h+y}.\hfill \end{array}$ (4.1)

Setting $\frac{\text{d}x}{\text{d}t}=0$ and $\frac{\text{d}y}{\text{d}t}=0$ in system (4.1), we get:

$\{\begin{array}{l}1+\eta y-\frac{\delta xy}{h+y}-x=0,\hfill \\ \rho y\left(1-y\right)-\frac{qxy}{h+y}=0.\hfill \end{array}$ (4.2)

On the one hand, the system (4.1) has a tumor-free equilibrium point $E_0\left(1,0\right)$ . On the other hand, when $y\ne 0$ , from the second equation of (4.2), we have:

$x=\frac{\rho \left(1-y\right)\left(h+y\right)}{q}.$ (4.3)

Substituting Equation (4.3) into the first equation of (4.2), we get:

$f\left(y\right)=\rho \left(\delta +1\right)y^2+\left[q\eta -\delta \rho -\rho +\rho h\right]y+\left(q-\rho h\right)=0.$ (4.4)

Let $\theta =q\eta -\delta \rho -\rho +\rho h$ , then Equation (4.4) becomes:

$f\left(y\right)=\rho \left(\delta +1\right)y^2+\theta y+\left(q-\rho h\right)=0.$

Since $f\left(y\right)>0$ when $y\ge 1$ , the positive solution of $f\left(y\right)=0$ only exists in the interval (0, 1). For the quadratic function $f\left(y\right)$ , $f\left(0\right)=q-\rho h$ , $f\left(1\right)=q\left(1+\eta \right)>0$ , and $f^{\prime }\left(y\right)=2\rho \left(\delta +1\right)y+\theta$ .

Let the discriminant of equation $f\left(y\right)=0$ be $\text{Δ}={\theta }^2-4\rho \left(\delta +1\right)\left(q-\rho h\right)$ , then:

1) When $q<\rho h$ , since $f\left(0\right)=q-\rho h<0$ , $f\left(y\right)=0$ has a unique positive solution in (0, 1):

$y_1=\frac{-\theta +\sqrt{\text{Δ}}}{2\rho \left(\delta +1\right)}.$

2) When $q=\rho h$ and $\eta <\frac{\delta +1-h}{h}$ , $f\left(y\right)=0$ has only one positive solution in (0, 1):

$y_2=\frac{\delta +1-h\eta -h}{\delta +1}.$

3) When $q>\rho h$ , $\eta <\frac{\left(\delta +1-h\right)\rho }{q}$ and $\text{Δ}>0$ , $f\left(y\right)=0$ has two positive solutions in (0, 1):

$y_{3,4}=\frac{-\theta \mp \sqrt{\text{Δ}}}{2\rho \left(\delta +1\right)}.$

4) When $q>\rho h$ , $\eta <\frac{\left(\delta +1-h\right)\rho }{q}$ and $\text{Δ}=0$ , $f\left(y\right)=0$ has a double positive solution in (0, 1):

$y_5=\frac{-\theta }{2\rho \left(\delta +1\right)}.$

Furthermore, substituting $y_j\left(j=1,\cdots ,5\right)$ into Equation (4.2) accordingly, we get:

$x_1=\frac{\left(\rho \delta +\rho -q\eta -\rho h+\sqrt{\text{Δ}}\right)\left(2\rho h\delta +\delta \rho +\rho -q\eta +\rho h+\sqrt{\text{Δ}}\right)}{4q\rho {\left(\delta +1\right)}^2},$

$x_2=\frac{\rho \left(2+h+h\eta \right)}{q\left(\delta +1\right)},$

$x_{3,4}=\frac{\left(\rho \delta +\rho +q\eta +\rho h\pm \sqrt{\text{Δ}}\right)\left(2\rho h\delta +\delta \rho +\rho -q\eta +\rho h\pm \sqrt{\text{Δ}}\right)}{4q\rho {\left(\delta +1\right)}^2},$

$x_5=\frac{\left(\rho \delta +\rho +q\eta +\rho h\right)\left(2\rho h\delta +\delta \rho +\rho -q\eta +\rho h\right)}{4q\rho {\left(\delta +1\right)}^2}.$

Thus, the system (4.1) has positive equilibrium points $E_j\left(x_j,y_j\right)$ (called tumor-present equilibrium points) under corresponding conditions. For $E_j$ , according to the relationship between $y_j$ and $f\left(y\right)$ , it is easy to know that the condition:

$f^{\prime }\left(y_i\right)>0\left(i=1,2,4\right),f^{\prime }\left(y_3\right)<0,f^{\prime }\left(y_5\right)=0.$

Note the condition:

$q>\rho h,\eta <\frac{\left(\delta +1-h\right)\rho }{q},\text{Δ}\ge 0.$

Equivalent to condition:

$\delta >h-1,q>\rho h,\eta \le \frac{\left(\delta +1-h\right)\rho -2\sqrt{\rho \left(\delta +1\right)\left(q-\rho h\right)}}{q}.$

Summarizing the above discussions, model (2.3) and (4.1) always have a tumor-free equilibrium point $E_0\left(1,0\right)$ . Furthermore, the existence of their tumor-present equilibrium points is given by the following conclusion:

Theorem 3. The existence of tumor-present equilibrium points of system (2.3) and (4.1) is as follows:

1) When $q<\rho h$ , system (4.1) has a unique tumor-present equilibrium point $E_1\left(x_1,y_1\right)$ ;

2) When $q=\rho h$ , $\delta >h-1$ and $\eta <\frac{\delta +1-h}{h}$ , system (4.1) has a unique tumor-present equilibrium point $E_2\left(x_2,y_2\right)$ ;

3) When $q>\rho h$ , $\delta >h-1$ and $\eta <\frac{\left(\delta +1-h\right)\rho -2\sqrt{\rho \left(\delta +1\right)\left(q-\rho h\right)}}{q}$ , system (4.1) has two different tumor-present equilibrium points $E_3\left(x_3,y_3\right)$ and $E_4\left(x_4,y_4\right)$ ;

4) When $q>\rho h$ , $\delta >h-1$ and $\eta =\frac{\left(\delta +1-h\right)\rho -2\sqrt{\rho \left(\delta +1\right)\left(q-\rho h\right)}}{q}$ , system (4.1) has a unique tumor-present equilibrium point $E_5\left(x_5,y_5\right)$ .

According to the expression of function $f\left(y\right)$ and Theorem 3, for the case $\delta \le h-1$ , system (4.1) has a unique tumor-present equilibrium point $E_1$ when $q<\rho h$ , and no tumor-present equilibrium points when $q\ge \rho h$ ; for the case $\delta >h-1$ , the existence conditions of equilibrium points of system (4.1) are relatively complex. To intuitively show these conditions, using the function

$\eta =\frac{\left(\delta +1-h\right)\rho -2\sqrt{\rho \left(\delta +1\right)\left(q-\rho h\right)}}{q}$ of $q$ , the corresponding equilibrium point existence regions are shown in Figure 1 on the plane $\left(q,\eta \right)$ , where:

Figure 1. Existence of tumor-present equilibrium points.

$D_0=\left\{\left(q,\eta \right)|q>\rho h,\eta >\frac{\left(\delta +1-h\right)\rho -2\sqrt{\rho \left(\delta +1\right)\left(q-\rho h\right)}}{q}\right\},$

$D_1=\left\{\left(q,\eta \right)|0<q<\rho h\right\}\cup \left\{\left(q,\eta \right)|q=\rho h,\eta <\frac{\delta +1-h}{h}\right\},$

$D_2=\left\{\left(q,\eta \right)|q>\rho h,\eta <\frac{\left(\delta +1-h\right)\rho -2\sqrt{\rho \left(\delta +1\right)\left(q-\rho h\right)}}{q}\right\},$

$D_3=\left\{\left(q,\eta \right)|q>\rho h,\eta =\frac{\left(\delta +1-h\right)\rho -2\sqrt{\rho \left(\delta +1\right)\left(q-\rho h\right)}}{q}\right\}.$

From Theorem 3, we know: when $\left(q,\eta \right)\in D_0$ , system (4.1) has no tumor-present equilibrium points; when $\left(q,\eta \right)\in D_1$ , system (4.1) has one tumor-present equilibrium point $E_1$ or $E_2$ ; when $\left(q,\eta \right)\in D_2$ , system (4.1) has two tumor-present equilibrium points $E_3$ and $E_4$ ; when $\left(q,\eta \right)\in D_3$ , system (4.1) has a unique tumor-present equilibrium point $E_5$ .

At the same time, it is easy to see from Figure 1 that the existence of tumor-present equilibrium points of system (4.1) depends on parameter $\eta$ (antigenicity). When $\eta \ge \frac{\delta +1-h}{h}$ , with the increase of $q$ , the tumor-present equilibrium points of system (4.1) change from none to one. When $\delta >h-1$ and $\eta <\frac{\delta +1-h}{h}$ , with the increase of $q$ , the tumor-present equilibrium points of system (4.1) change from none to two, then to one. At this time, system (4.1) undergoes saddle-node bifurcation when changing from the tumor-free equilibrium point to two tumor-present equilibrium points.

5. Global Stability Analysis of ODE Model (4.1)

We start with the local dynamic behavior of system (4.1). In a two-dimensional planar system, there are two types of locally asymptotically stable equilibrium points: foci and nodes. Different types lead to different ways of the system’s trajectories converging to the equilibrium points. In this section, we first discuss the stability of the tumor-free equilibrium point $E_0\left(1,0\right)$ and the tumor-present equilibrium point $E_i\left(x_i,y_i\right)$ of system (4.1), and finally analyze the type of stable equilibrium points.

5.1. Stability of the Tumor-Free Equilibrium Point

First, the Jacobian matrix of system (4.1) at the tumor-free equilibrium point $E_0$ is:

$J\left(E_0\right)=\left(\begin{array}{cc}-1& \eta -\frac{\delta }{h}\\ 0& \rho -\frac{q}{h}\end{array}\right).$

Its two eigenvalues are ${\lambda }_1=-1$ and ${\lambda }_2=\rho -\frac{q}{h}$ . Thus, when $\rho <\frac{q}{h}$ , $E_0$ is locally asymptotically stable; when $\rho >\frac{q}{h}$ , $E_0$ is unstable.

When $\rho =\frac{q}{h}$ , the two eigenvalues of the Jacobian matrix of system (4.1) at point $E_0$ are ${\lambda }_1=-1$ and ${\lambda }_2=0$ , so $E_0$ is a high-order equilibrium point. In this case, to study the behavior of $E_0$ on the invariant set $D$ , we consider the following two cases:

Case 1: $\eta -\frac{\delta }{h}\ne 0$

First, perform a translation transformation on system (4.1): let $u=x-1$ , $v=y$ , then system (4.1) becomes:

$\{\begin{array}{l}\frac{\text{d}u}{\text{d}t}=\eta v-u-\frac{\delta uv}{h+v}-\frac{\delta v}{h+v},\hfill \\ \frac{\text{d}v}{\text{d}t}=\frac{q}{h}v-\frac{q}{h}{v}^2-\frac{quv}{h+v}-\frac{qv}{h+v}.\hfill \end{array}$ (5.1)

translating the equilibrium point $E_0$ to the origin. Then, transform the linear part of system (5.1) into Jordan canonical form. From the eigenvalues −1 and 0, we get two eigenvectors ${\left(1,0\right)}^{\text{T}}$ and ${\left(\eta -\frac{\delta }{h},1\right)}^{\text{T}}$ of the linear part of system (5.1). Thus, perform an invertible linear transformation on system (5.1):

$\left(\begin{array}{l}u\hfill \\ v\hfill \end{array}\right)=\left(\begin{array}{cc}1& \eta -\frac{\delta }{h}\\ 0& 1\end{array}\right)\left(\begin{array}{c}z\\ w\end{array}\right),$ (5.2)

then system (5.1) becomes:

$\{\begin{array}{l}\frac{\text{d}z}{\text{d}t}=-z+G\left(z,w\right),\hfill \\ \frac{\text{d}w}{\text{d}t}=H\left(z,w\right).\hfill \end{array}$ (5.3)

where

$\begin{array}{c}G\left(z,w\right)=-\frac{\delta }{h}\left(z+\left(\eta -\frac{\delta }{h}\right)w\right)w+\frac{\delta }{h^2}\left(z+\left(\eta -\frac{\delta }{h}\right)w+1\right)w^2\\ -\frac{\delta }{h^3}\left(z+\left(\eta -\frac{\delta }{h}\right)w+1\right)w^3,\end{array}$

$H\left(z,w\right)=-\frac{q}{h}{w}^2-\frac{q}{h}\left(z+\left(\eta -\frac{\delta }{h}\right)w\right)w+\frac{q}{h^2}\left(z+\left(\eta -\frac{\delta }{h}\right)w+1\right)w^2-\frac{q}{h^3}{w}^3.$

By the Center Manifold Theorem [30], assume the local center manifold of system (5.3) at the origin is:

$z=h\left(w\right)=a{w}^2+b{w}^3+o\left(w^3\right),h\left(0\right)=0,h^{\prime }\left(0\right)=0.$ (5.4)

Using the invariance of the center manifold, we solve:

$a=-\frac{\delta }{h}\left(\eta -\frac{\delta }{h}-\frac{1}{h}\right),$

$b=\frac{2q\delta }{h^2}\left(1+\eta -\frac{1}{h}-\frac{\delta }{h}\right)\left(\eta -\frac{\delta }{h}-\frac{1}{h}\right)+\frac{{\delta }^2}{h^2}\left(\eta -\frac{\delta }{h}-\frac{1}{h}\right)+\frac{\delta }{h^2}\left(\eta -\frac{\delta }{h}\right)-\frac{\delta }{h^3}.$

Thus, the center manifold of system (4.1) at the origin is:

$\begin{array}{l}z=h\left(w\right)=-\frac{\delta }{h}\left(\eta -\frac{\delta }{h}-\frac{1}{h}\right)w^2+\frac{\delta }{h^2}[2q\left(1+\eta -\frac{1}{h}-\frac{\delta }{h}\right)\left(\eta -\frac{\delta }{h}-\frac{1}{h}\right)\\ +\delta \left(\eta -\frac{\delta }{h}-\frac{1}{h}\right)+\left(\eta -\frac{\delta }{h}\right)-\frac{1}{h}]w^3+o\left(w^3\right).\end{array}$ (5.5)

Substitute Equation (5.5) into the second equation of system (5.3):

$\begin{array}{c}\frac{\text{d}w}{\text{d}t}=A{w}^2+B{w}^3+o\left(w^3\right)\\ =-\frac{q}{h}\left[1+\eta -\frac{\delta }{h}-\frac{1}{h}\right]w^2+\frac{q}{h}\left[\frac{\delta }{h}\left(\eta -\frac{\delta }{h}-\frac{1}{h}\right)+\frac{1}{h}\left(\eta -\frac{\delta }{h}\right)-\frac{1}{h^2}\right]w^3\\ +o\left(w^3\right).\end{array}$ (5.6)

where

$A=-\frac{q}{h}\left(1+\eta -\frac{\delta }{h}-\frac{1}{h}\right),$

$B=\frac{q}{h}\left[\frac{\delta }{h}\left(\eta -\frac{\delta }{h}-\frac{1}{h}\right)+\frac{1}{h}\left(\eta -\frac{\delta }{h}\right)-\frac{1}{h^2}\right].$

If $A>0$ , i.e., $\eta <\frac{\delta +1-h}{h}$ , then the tumor-free equilibrium point $E_0\left(1,0\right)$ is unstable;

If $A<0$ , i.e., $\eta >\frac{\delta +1-h}{h}$ , then the tumor-free equilibrium point $E_0\left(1,0\right)$ is locally asymptotically stable;

If $A=0$ , i.e., $\eta =\frac{\delta +1-h}{h}$ , then the stability of the tumor-free equilibrium point $E_0\left(1,0\right)$ needs further discussion.

When $\eta =\frac{\delta +1-h}{h}$ holds, substituting into Equation (5.6), the stability of $E_0$ is determined by the higher-order terms.

$\begin{array}{c}\frac{\text{d}w}{\text{d}t}=\frac{q}{h}\left[\frac{\delta }{h}\left(\eta -\frac{\delta }{h}-\frac{1}{h}\right)+\frac{1}{h}\left(\eta -\frac{\delta }{h}\right)-\frac{1}{h^2}\right]\\ =\frac{q}{h^2}\left(-\delta -1\right).\end{array}$ (5.7)

According to the Local Center Manifold Theorem [30], the tumor-free equilibrium point $E_0\left(1,0\right)$ is locally asymptotically stable at this time.

Case 2: $\eta -\frac{\delta }{h}=0$

When $\eta -\frac{\delta }{h}=0$ , i.e., $\eta =\frac{\delta }{h}$ , substituting into system (4.1), we get:

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=1+\frac{\delta }{h}y-\frac{\delta xy}{h+y}-x,\hfill \\ \frac{\text{d}y}{\text{d}t}=\rho y\left(1-y\right)-\frac{qxy}{h+y}.\hfill \end{array}$ (5.8)

Substituting into Equation (5.6), we get:

$\frac{\text{d}w}{\text{d}t}=-\frac{q}{h}\left(1-\frac{1}{h}\right)w^2+\frac{q}{h}\left[\frac{\delta }{h}\left(\eta -\frac{\delta }{h}-\frac{1}{h}\right)+\frac{1}{h}\left(\eta -\frac{\delta }{h}\right)-\frac{1}{h^2}\right]w^3+o\left(w^3\right).$ (5.9)

If $h<1$ , the tumor-free equilibrium point $E_0\left(1,0\right)$ is unstable;

If $h>1$ , the tumor-free equilibrium point $E_0\left(1,0\right)$ is locally asymptotically stable;

If $h=1$ , the stability of the tumor-free equilibrium point $E_0\left(1,0\right)$ is determined by the higher-order terms.

$\frac{\text{d}w}{\text{d}t}=\frac{q}{h}\left[\frac{\delta }{h}\left(\eta -\frac{\delta }{h}-\frac{1}{h}\right)+\frac{1}{h}\left(\eta -\frac{\delta }{h}\right)-\frac{1}{h^2}\right]w^3=q\left(-\delta -1\right)w^3+o\left(w^3\right).$ (5.10)

According to the Local Center Manifold Theorem [30], the tumor-free equilibrium point $E_0\left(1,0\right)$ is locally asymptotically stable at this time.

In summary, the stability of the tumor-free equilibrium point is given by the following Theorem 4:

Theorem 4. (1) The tumor-free equilibrium point $E_0\left(1,0\right)$ of system (4.1) is locally asymptotically stable if any of the following conditions is satisfied:

a) $\rho <\frac{q}{h}$ ;

b) $\rho =\frac{q}{h}$ , $\eta \ne \frac{\delta }{h}$ , $\eta \ge \frac{\delta +1-h}{h}$ ;

c) $\rho =\frac{q}{h}$ , $\eta =\frac{\delta }{h}$ , $h\ge 1$ .

2) The tumor-free equilibrium point $E_0\left(1,0\right)$ of system (4.1) is unstable if any of the following conditions is satisfied:

a) $\rho >\frac{q}{h}$ ;

b) $\rho =\frac{q}{h}$ , $\eta \ne \frac{\delta }{h}$ , $\eta <\frac{\delta +1-h}{h}$ ;

c) $\rho =\frac{q}{h}$ , $\eta =\frac{\delta }{h}$ , $h<1$ .

5.2. Stability of Tumor-Present Equilibrium Points

For tumor-present equilibrium point $E_i\left(x_i,y_i\right),i=1,2,3,4,5$ , discussed in Section 4, the Jacobian matrix of system (4.1) at $E_i\left(x_i,y_i\right),i=1,2,3,4,5$ , is:

$J\left(E_i\right)=\left(\begin{array}{cc}-1-\frac{\delta y_i}{h+y_i}& \eta -\frac{\delta h{x}_i}{{\left(h+y_i\right)}^2}\\ -\frac{q{y}_i}{h+y_i}& \rho \left(1-2y_i\right)-\frac{qh{x}_i}{{\left(h+y_i\right)}^2}\end{array}\right).$

From Equation (4.2), we know $x_i=\frac{\rho \left(1-y_i\right)\left(h+y_i\right)}{q}$ . Its trace and determinant are:

$\begin{array}{c}trJ\left(E_i\right)=-1-\frac{\delta y_i}{h+y_i}+\rho \left(1-2y_i\right)-\frac{qh{x}_i}{{\left(h+y_i\right)}^2}\\ =-1-\frac{\delta y_i}{h+y_i}+\frac{\rho y_i\left(1-h-2y_i\right)}{h+y_i}\end{array}$

$\begin{array}{c}detJ\left(E_i\right)=\left(-1-\frac{\delta y_i}{h+y_i}\right)\left[\rho \left(1-2y_i\right)-\frac{qh{x}_i}{{\left(h+y_i\right)}^2}\right]+\frac{q{y}_i}{h+y_i}\left(\eta -\frac{\delta h{x}_i}{{\left(h+y_i\right)}^2}\right)\\ =\frac{y_i{f}^{\prime }\left(y_i\right)}{h+y_i}\end{array}$ (5.11)

From the geometric properties of $f\left(y\right)$ , we have

$f^{\prime }\left(y_1\right)>0,f^{\prime }\left(y_2\right)>0,f^{\prime }\left(y_3\right)<0,f^{\prime }\left(y_4\right)>0.$

From the above Equation (5.11), it can be obtained:

$detJ\left(E_1\right)>0,detJ\left(E_2\right)>0,detJ\left(E_3\right)<0,detJ\left(E_4\right)>0.$

From the above Equation (5.11), $trJ\left(E_i\right)$ is negative if and only if one of the following conditions holds:

1) $h\ge 1$ ;

2) $h<1$ , and $y\ge \frac{1-h}{2}$ ;

3) $h<1$ , and $y<\frac{1-h}{2},\rho >\frac{h+y+\delta y}{y\left(1-h-2y\right)}$ .

From the above analysis, regarding the stability of the tumor-present equilibrium points of model (4.1), we have the following conclusion:

Theorem 5. 1) System (4.1) has a positive equilibrium point $E_3\left(x_3,y_3\right)$ which is a saddle point and thus unstable;

2) System (4.1) has positive equilibrium points $E_1\left(x_1,y_1\right)$ , $E_2\left(x_2,y_2\right)$ , $E_4\left(x_4,y_4\right)$ . When condition (1), (2) or (3) is satisfied, system (4.1) is locally asymptotically stable at the tumor-present equilibrium point $E_1\left(x_1,y_1\right)$ , $E_2\left(x_2,y_2\right)$ , $E_4\left(x_4,y_4\right)$ ;

3) When $detJ\left(E_5\right)=0$ , system (4.1) has a high-order tumor-present equilibrium point, and its stability needs further discussion.

At this equilibrium point, the eigenvalues of the Jacobian matrix are ${\lambda }_1=a_{11}+a_{22}$ and ${\lambda }_2=0$ , and the corresponding eigenvectors are ${\left(1,{\xi }_1\right)}^{\text{T}}$ and ${\left({\xi }_2,1\right)}^{\text{T}}$ , where ${\xi }_1=\frac{a_{22}}{a_{12}}$ , ${\xi }_2=-\frac{a_{12}}{a_{11}}$ . Next, we study the dynamic properties at this high-order tumor-present equilibrium point $E_5\left(x_5,y_5\right)$ .

To study the dynamic behavior of system (4.1) near the tumor-present equilibrium point $E_i\left(x_i,y_i\right)$ , perform a translation transformation on system (4.1): let $u_1=x-x^{\text{*}}$ , $v_1=y-y^{\text{*}}$ and substitute into system (4.1), we get:

$\{\begin{array}{l}\frac{\text{d}{u}_1}{\text{d}t}=a_{11}{u}_1+a_{12}{v}_1+M\left(u_1,v_1\right),\hfill \\ \frac{\text{d}{v}_1}{\text{d}t}=a_{21}{u}_1+a_{22}{v}_1+N\left(u_1,v_1\right).\hfill \end{array}$ (5.12)

where

$a_{11}=-1-\frac{\delta y^{\text{*}}}{h+y^{\text{*}}},a_{12}=\eta -\frac{\delta h{x}^{\text{*}}}{{\left(h+y^{\text{*}}\right)}^2},a_{21}=-\frac{q{y}^{\text{*}}}{h+y^{\text{*}}},a_{22}=\rho \left(1-2y^{\text{*}}\right)-\frac{qh{x}^{\text{*}}}{{\left(h+y^{\text{*}}\right)}^2}.$

$M\left(u_1,v_1\right)=-\frac{\delta }{h+y^{\text{*}}}{u}_1{v}_1+\frac{\delta x^{\text{*}}}{{\left(h+y^{\text{*}}\right)}^2}{v}_1^2-\frac{\delta }{{\left(h+y^{\text{*}}\right)}^2}{u}_1{v}_1^2+\frac{2\delta x^{\text{*}}}{{\left(h+y^{\text{*}}\right)}^3}{v}_1^3,$

$N\left(u_1,v_1\right)=-\rho v_1^2-\frac{q}{h+y^{\text{*}}}{u}_1{v}_1+\frac{q{x}^{\text{*}}}{{\left(h+y^{\text{*}}\right)}^2}{v}_1^2-\frac{q}{{\left(h+y^{\text{*}}\right)}^2}{u}_1{v}_1^2+\frac{2q{x}^{\text{*}}}{{\left(h+y^{\text{*}}\right)}^3}{v}_1^3.$

Perform a linear transformation on system (5.12):

$\left(\begin{array}{l}{u}_1\hfill \\ v_1\hfill \end{array}\right)=\left(\begin{array}{ll}1\hfill & {\xi }_2\hfill \\ {\xi }_1\hfill & 1\hfill \end{array}\right)\left(\begin{array}{l}{z}_1\hfill \\ w_1\hfill \end{array}\right),$

so:

$\{\begin{array}{l}\frac{\text{d}{z}_1}{\text{d}t}={\lambda }_1{z}_1+\overline{M}\left(z_1,w_1\right),\hfill \\ \frac{\text{d}{w}_1}{\text{d}t}=\overline{N}\left(z_1,w_1\right).\hfill \end{array}$(5.13)

where

$\overline{M}\left(z_1,w_1\right)=\frac{1}{1-{\xi }_1{\xi }_2}\left(M\left(z_1+{\xi }_2{w}_1,{\xi }_1{z}_1+w_1\right)-{\xi }_{2}N\left(z_1+{\xi }_2{w}_1,{\xi }_1{z}_1+w_1\right)\right),$

$\overline{N}\left(z_1,w_1\right)=\frac{1}{1-{\xi }_1{\xi }_2}\left(N\left(z_1+{\xi }_2{w}_1,{\xi }_1{z}_1+w_1\right)-{\xi }_{1}M\left(z_1+{\xi }_2{w}_1,{\xi }_1{z}_1+w_1\right)\right).$

Similar to Section 5.1, the local center manifold of system (5.13) at the origin is found to be:

$\overline{z}=h\left(w\right)=\overline{a}{w}^2+o\left(w^2\right).$ (5.14)

where

$\overline{a}=-\frac{1}{{\lambda }_1\left(1-{\xi }_1{\xi }_2\right)}\left[\rho {\xi }_2+{\xi }_2\left(\frac{q{\xi }_2-\delta }{h+y^{\text{*}}}\right)+\frac{\delta x^{\text{*}}-q{\xi }_2{x}^{\text{*}}}{{\left(h+y^{\text{*}}\right)}^2}\right].$

Substitute Equation (5.14) into the second equation of (5.13), we have:

$\begin{array}{c}\frac{\text{d}w}{\text{d}t}=\frac{1}{1-{\xi }_1{\xi }_2}[\left(-\rho -\frac{q}{h+y^{\text{*}}}{\xi }_2+\frac{q{x}^{\text{*}}}{{\left(h+y^{\text{*}}\right)}^2}\right)\\ -{\xi }_1\left(-\frac{\delta }{h+y^{\text{*}}}{\xi }_2+\frac{\delta x^{\text{*}}}{{\left(h+y^{\text{*}}\right)}^2}\right)]w^2+o\left(w^2\right)\\ \triangleq Q{w}^2+o\left(w^2\right).\end{array}$

(H1) $Q<0$ .

Theorem 6. When conditions $trJ\left(E_5\right)<0$ and $detJ\left(E_5\right)=0$ are satisfied, system (4.1) has a high-order tumor-present equilibrium point $E_5\left(x_5,y_5\right)$ . When condition (H1) is satisfied, the tumor-present equilibrium point is locally asymptotically stable.

Finally, we analyze the type of the stable tumor-present equilibrium point $E^{\text{*}}\left(x^{\text{*}},y^{\text{*}}\right)$ . Theorems 5 and 6 show that there exist locally asymptotically stable tumor-present equilibrium points. For the tumor-present equilibrium point $E_i\left(x_i,y_i\right)$ of system (4.1), when $\kappa =trJ{\left(E_i\right)}^2-4detJ\left(E_i\right)<0$ , $E_i$ is a focus; when $\kappa \ge 0$ , $E_i$ is a node.

Choose $B\left(x,y\right)=\frac{1}{xy}$ as the Dulac function, and denote the functions on the right-hand side of system (4.1) as $M\left(x,y\right)$ and $N\left(x,y\right)$ , respectively. Then:

$\frac{\partial \left(BM\right)}{\partial x}+\frac{\partial \left(BN\right)}{\partial y}=-\frac{1}{x^{2}y}\left(1+\eta y-\frac{\delta xy}{h+y}\right)-\frac{1}{x{y}^2}\left(\rho y\left(1-y\right)-\frac{qxy}{h+y}\right)<0,$

For all $x>0$ , $y>0$ . Thus, according to the Bendixson-Dulac Theorem, system (4.1) has no closed trajectories, i.e., no periodic solutions.

Remark. The tumor-free equilibrium point $E_0$ of system (4.1) is globally asymptotically stable on the region $D$ if it is locally asymptotically stable and there are no other stable equilibrium points.

6. Dynamic Analysis of DDE Model (2.3)

In this section, we first study the local dynamic behavior of system (2.3). Then we discuss the stability change of the equilibrium point with respect to delay, and verify the impact of delay through numerical simulation, and provide some biological explanations.

The tumor-free equilibrium point of model (2.3) is $P_0\left(1,0\right)$ , and the possible positive equilibrium points are denoted as $P_k\left(E_k,T_k\right)\left(k=1,2,3,4,5\right)$ . This section mainly discusses the local stability of the tumor-free equilibrium point and positive equilibrium points and the Hopf bifurcation of model (2.3).

6.1. Local Stability of $P_0$

In model (2.3), let $x=E-E^{\prime }$ , $y=T-T^{\prime }$ , then the linear system of model (2.3) at any equilibrium point $P^{\prime }\left(E^{\prime },T^{\prime }\right)$ is:

$\{\begin{array}{l}\frac{\text{d}x\left(t\right)}{\text{d}t}=-\left(1+\frac{\delta T^{\prime }}{h+T^{\prime }}\right)x-\frac{\delta h{E}^{\prime }}{{\left(h+T^{\prime }\right)}^2}y\left(t\right)+\eta y\left(t-\tau \right),\hfill \\ \frac{\text{d}y\left(t\right)}{\text{d}t}=-\frac{q{T}^{\prime }}{h+T^{\prime }}x\left(t\right)+\left(\rho \left(1-2T^{\prime }\right)-\frac{qh{E}^{\prime }}{{\left(h+T^{\prime }\right)}^2}\right)y\left(t\right).\hfill \end{array}$ (6.1)

Theorem 7. For all $\tau \ge 0$ , if $\rho <\frac{q}{h}$ , then the tumor-free equilibrium point $P_0$ is locally asymptotically stable; if $\rho >\frac{q}{h}$ , then $P_0$ is unstable.

Proof. At the tumor-free equilibrium point $P_0\left(1,0\right)$ , the characteristic equation of Equation (6.1) is $\left(\lambda +1\right)\left(\lambda -\left(\rho -\frac{q}{h}\right)\right)=0$ , and the eigenvalues are −1 and $\rho -\frac{q}{h}$ , which are independent of $\tau$ . Thus, for any $\tau \ge 0$ , when $\rho <\frac{q}{h}$ , $P_0$ is locally asymptotically stable; when $\rho >\frac{q}{h}$ , $P_0$ is unstable.

6.2. Local Stability and Hopf Bifurcation of $P_k$

In Equation (6.1), let $x\left(t\right)=C_1{\text{e}}^{\lambda t}$ , $y\left(t\right)=C_2{\text{e}}^{\lambda t}$ (where $C_1$ and $C_2$ are non-negative constants), then the characteristic equation of Equation (6.1) at the positive equilibrium point $P_k\left(E_k,T_k\right)\left(k=1,2,3,4,5\right)$ is:

${\lambda }^2+a_1\lambda +a_2+a_3{\text{e}}^{-\lambda \tau }=0,$ (6.2)

where

$\begin{array}{l}{a}_1=1+\frac{\delta y^{\text{*}}}{h+y^{\text{*}}}-\frac{\rho y^{\text{*}}\left(1-h-2y^{\text{*}}\right)}{h+y^{\text{*}}},\\ a_2=\left(-1-\frac{\delta y^{\text{*}}}{h+y^{\text{*}}}\right)\left[\rho \left(1-2y^{\text{*}}\right)-\frac{qh{x}^{\text{*}}}{{\left(h+y^{\text{*}}\right)}^2}\right]+\frac{q{y}^{\text{*}}}{h+y^{\text{*}}}\left(\eta -\frac{\delta h{x}^{\text{*}}}{{\left(h+y^{\text{*}}\right)}^2}\right),\\ a_3=-\frac{q\eta y^{\text{*}}}{h+y^{\text{*}}}.\end{array}$

Theorem 8. 1) When $\tau =0$ , the positive equilibrium point of model (2.3) is stable if and only if $a_1>0$ and $a_2+a_3>0$ ;

2) There exists ${\tau }_0^{\text{*}}=\frac{1}{{\omega }_0}\text{arccos}\left(\frac{{\omega }_0^2-a_2}{a_3}\right)$ such that when $0<\tau <{\tau }_0^{\text{*}}$ , if

$a_1^2-2a_2<0$ and $a_1^4-4a_1^2{a}_2+4a_2^2>0$ , then $P_k\left(k=1,2,3,4,5\right)$ is stable; when $\tau >{\tau }_0^{\text{*}}$ , $P_k\left(k=1,2,3,4,5\right)$ is unstable, where $\pm i{\omega }_0$ is a pair of pure imaginary roots of Equation (6.2) and

${\omega }_0^2=\frac{1}{2}\left[-\left(a_1^2-2a_2\right)\pm \sqrt{{\left(a_1^2-2a_2\right)}^2-4\left(a_2^2-a_3^2\right)}\right];$

3) When $\tau ={\tau }_0^{\text{*}}$ , model (2.3) undergoes Hopf bifurcation at $\tau ={\tau }_0^{\text{*}}$ .

Proof. 1) When $\tau =0$ , the equilibrium point of model (2.3) is locally asymptotically stable if and only if all roots of the equation $H\left(\lambda \right)={\lambda }^2+a_1\lambda +a_2+a_3=0$ have negative real parts. According to the Routh-Hurwitz criterion [30], all roots of $H\left(\lambda \right)=0$ have negative real parts if and only if $a_1>0$ and $a_2+a_3>0$ .

2) Let $\lambda =i\omega \left(\omega >0\right)$ be a root of Equation (6.2). Substitute it into Equation (6.2) and separate the real and imaginary parts:

$\{\begin{array}{l}{a}_3\cos \left(\omega \tau \right)={\omega }^2-a_2,\hfill \\ a_3\sin \left(\omega \tau \right)=a_1\omega ,\hfill \end{array}$ (6.3)

Square and add both sides of the two equations in (6.3), and let $\xi ={\omega }^2$ , then we have:

${\xi }^2+\left(a_1^2-2a_2\right)\xi +\left(a_2^2-a_3^2\right)=0.$ (6.4)

If:

$a_1^2-2a_2<0\text{and}{a}_2^2-a_3^2<0,$ (6.5)

or:

$a_1^2-2a_2>0\text{and}{\left(a_1^2-2a_2\right)}^2-4\left(a_2^2-a_3^2\right)=a_1^4-4a_1^2{a}_2+4a_3^2\ge 0,$ (6.6)

holds, then Equation (6.4) has at least one positive root. Also, $a_3<0$ , and from (1), we assume $a_2+a_3>0$ , so $a_2>-a_3>0$ , thus $a_2^2-a_3^2>0$ , which contradicts (6.5). Therefore, Equation (6.4) has at least one positive real root ${\xi }_0$ if and only if (6.6) holds,

${\xi }_0={\omega }_0^2=\frac{1}{2}\left[-\left(a_1^2-2a_2\right)\pm \sqrt{{\left(a_1^2-2a_2\right)}^2-4\left(a_2^2-a_3^2\right)}\right],$

Then $H\left(\lambda \right)=0$ has a pair of pure imaginary roots $\pm i{\omega }_0$ . and ${\omega }_0=\sqrt{{\xi }_0}$ , From Equation (6.4), we get ${\tau }_0^{\text{*}}=\frac{1}{{\omega }_0}\text{arccos}\left(\frac{{\omega }_0^2-a_2}{a_3}\right)$ . Thus, when $\tau \in \left[0,{\tau }_0^{\text{*}}\right)$ , the equilibrium point $P_k$ is stable; when $\tau >{\tau }_0^{\text{*}}$ , $P_k$ is unstable.

3) From (2), Equation (6.2) has a pair of pure imaginary roots $\pm i{\omega }_0$ . Let $\lambda \left(\tau \right)=\delta \left(\tau \right)+i\omega \left(\tau \right)$ be the root of (6.2) under the conditions $\delta \left({\tau }_0^{\text{*}}\right)=0$ and $\omega \left({\tau }_0^{\text{*}}\right)={\omega }_0$ . Differentiate both sides of Equation (6.2) with respect to $\tau$ :

$\left(2\lambda +a_1-a_3\tau {\text{e}}^{-\lambda \tau }\right)\frac{\text{d}\lambda }{\text{d}\tau }=a_3\lambda {\text{e}}^{-\lambda \tau }.$

so:

${\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}=\frac{2\lambda +a_1-a_3\tau {\text{e}}^{-\lambda \tau }}{a_3\lambda {\text{e}}^{-\lambda \tau }}=-\frac{2\lambda +a_3}{\lambda \left({\lambda }^2+a_1\lambda +a_3\right)}-\frac{\tau }{\lambda }.$

Let $H\left({\xi }_0\right)={\xi }_0^2+\left(a_1^2-2a_2\right){\xi }_0+a_2^2-a_3^2$ .

Then:

$\begin{array}{c}\text{sign}\left\{{{\left(\frac{\text{d}\mathrm{Re}\left(\lambda \right)}{\text{d}\tau }\right)}^{-1}|}_{\tau ={\tau }_0^{\text{*}}}\right\}=\text{sign}\left\{{{\left(\frac{\text{d}\mathrm{Re}\left(\lambda \right)}{\text{d}\tau }\right)}^{-1}|}_{\lambda =i{\omega }_0}\right\}\\ =\text{sign}\left\{\mathrm{Re}\left({{\left(\frac{\text{d}\lambda }{\text{d}\tau }\right)}^{-1}|}_{\lambda =i{\omega }_0}\right)\right\}\\ =\text{sign}\left\{\frac{2{\omega }_0^2+a_1^2-2a_2}{{\left({\omega }_0^2-a_2\right)}^2+a_1^2{\omega }_0^2}\right\}\\ =\text{sign}\left\{\frac{H^{\prime }\left({\omega }_0^2\right)}{{\left({\omega }_0^2-a_2\right)}^2+a_1^2{\omega }_0^2}\right\}.\end{array}$

Thus, when $H^{\prime }\left({\omega }_0^2\right)=2{\omega }_0^2+a_1^2-2a_2>0$ , i.e., the transversality condition ${\frac{\text{d}\mathrm{Re}\left(\lambda \right)}{\text{d}\tau }|}_{\tau ={\tau }_0^{\text{*}}}={\frac{\text{d}\delta \left(\tau \right)}{\text{d}\tau }|}_{\tau ={\tau }_0^{\text{*}}}>0$ holds, so model (2.3) undergoes Hopf bifurcation at $\tau ={\tau }_0^{\text{*}}$ .

6.3. Numerical Simulation

In this section, Matlab is used for numerical simulation analysis of the conclusions.

For system (2.3), select parameters $\eta =0.5$ , $\delta =1.5$ , $q=0.8$ , $\rho =0.4$ , $h=1.0$ and different $\tau$ values for numerical simulation. As shown in Figure 2, all solutions converge to the point $E_0\left(1,0\right)$ , which indicates that the tumor-free equilibrium point $E_0\left(1,0\right)$ is locally asymptotically stable, i.e., tumor cells will be eliminated by the immune system.

(a)

(b)

(c)

Figure 2. Temporal dynamics and phase portrait of System (2.3) at the tumor-free equilibrium.

Select parameters $\eta =1.1$ , $\delta =1.8$ , $q=0.9$ , $\rho =0.9$ , $h=1.2$ . As shown in Figure 3, there exists a tumor-present equilibrium point $E^{\text{*}}\left(x^{\text{*}},y^{\text{*}}\right)=\left(0.99,0.37\right)$ . The stability of system (2.3) at the tumor-present equilibrium point $E^{\text{*}}\left(x^{\text{*}},y^{\text{*}}\right)$ varies with different $\tau$ values.

(a)

(b)

Figure 3. Stability of System (2.3) at the endemic equilibrium $E^{\text{*}}\left(x^{\text{*}},y^{\text{*}}\right)$ .

Select parameters $\eta =0.65$ , $\delta =2.3$ , $q=1.3$ , $\rho =0.9$ , $h=1.4$ for numerical simulation. As shown in Figure 4, at this time $q<\rho h$ , so the tumor-free equilibrium point $E_0\left(1,0\right)$ is locally asymptotically stable, and $trJ\left(E_2\right)<0$ , $detJ\left(E_2\right)>0$ , i.e., conditions are satisfied. This indicates that the tumor-present equilibrium point $E_2=\left(0.8668,0.2336\right)$ is locally asymptotically stable, i.e., the system exhibits bistability at this time.

Figure 4. Bistability of System (2.3).

7. Conclusion

This paper discusses a class of models of the interaction between tumors and the immune system with antigen delay and Michaelis-Menten type inhibition terms. For the convenience of analysis, the proposed model is first subjected to dimensionless transformation to simplify the model. On the basis of obtaining the existence conditions of the tumor-free equilibrium point and tumor-present equilibrium point of the model, the local dynamic behavior of the tumor-free equilibrium point is analyzed by using the Center Manifold Theorem, and the existence of periodic solutions of the system is excluded by constructing a Dulac function, so as to obtain the global dynamic behavior of the model. The analysis results show that time delay has an important impact on the stability of the positive equilibrium point. The asymptotic stability or instability of the positive equilibrium point depends on the size of the time delay $T$ . There exists a critical value ${\tau }_0^{\text{*}}$ such that when $\tau <{\tau }_0^{\text{*}}$ , the positive equilibrium point is stable; when $\tau >{\tau }_0^{\text{*}}$ , the positive equilibrium point is unstable, and the system may produce Hopf bifurcation when the positive equilibrium point is unstable. The introduction of antigenicity can cause saddle-node bifurcation of the model and the phenomenon that the tumor-present equilibrium point and the tumor-free equilibrium point are stable at the same time. The occurrence of this bistability means that the growth and development outcome of the tumor depends on their initial state (as shown in Figure 4). This bistability phenomenon indicates that tumor growth is related to the initial state: when the initial values of both tumor cells and effector cells are very small, or the initial value of tumor cells is very small and the initial value of effector cells is very large, the tumor cells will eventually disappear; when the initial value of tumor cells is very large and the initial value of effector cells is very small, or the initial values of both tumor cells and effector cells are very large, the tumor tends to the positive equilibrium point. This bistability creates an “opportunity window” in treatments such as immunotherapy or chemotherapy, shifting the system from a “tumor equilibrium point” to a “tumor-free equilibrium point”; this conclusion provides theoretical support for the clinical strategy of “early intervention and cell reinfusion”. In addition, this paper further gives the judgment conditions for whether the stable tumor-present equilibrium point is a focus or a node, providing a theoretical basis for clinically judging the nature of the final development stage of the tumor. According to the results obtained, the dynamic behavior of the proposed model is complex to a certain extent, which reflects the interaction mechanism between tumor cells and immune cells to a certain extent, but it still cannot fully show some more complex dynamic phenomena (such as the existence of B-T bifurcation). This needs to be further explored for the established model.

Acknowledgements

Sincere thanks to the members of JAMP for their professional performance, and special thanks to managing editor Hellen XU for a rare attitude of high quality.