Juanlan Wang
College of Mathematics and Statistics, Northwest Normal University, Lanzhou, China.
DOI: 10.4236/jamp.2026.147132   PDF    HTML   XML   11 Downloads   90 Views   Citations

Abstract

Tumor cells interact with the immune system via multiple nonlinear feedbacks, which can trigger various clinically relevant dynamical outcomes including tumor clearance, immune escape, long-term tumor dormancy, and recurrent periodic oscillations. This study investigates a two-stage T-lymphocyte-tumor model. First, the positivity and boundedness of solutions with nonnegative initial values are rigorously established. Next, the existence conditions and local stability criteria for the tumor-free equilibrium and positive equilibria are derived. Then, the Sotomayor theorem is applied to analyze the transcritical bifurcation at the boundary equilibrium. The Poincare-Andronov-Hopf bifurcation theory is further used to obtain the conditions for a Hopf bifurcation at a positive equilibrium. Moreover, normal form theory and the center manifold theorem are employed to characterize the direction of the Hopf bifurcation and the stability of the resulting periodic solutions. Finally, MATLAB simulations verify the stability of the equilibria and the Hopf bifurcation.

Keywords

Tumor-Immune Model, Equilibrium, Stability, Routh-Hurwitz Criterion, Transcritical Bifurcation, Hopf Bifurcation

Share and Cite:

1. Introduction

Tumor initiation, progression, and recurrence are not caused only by the autonomous growth of a single cell population. Instead, they arise from the long-term coupled evolution of tumor cells, effector immune cells, cytokines, stromal components, and metabolic resources. The immune system can recognize and eliminate abnormal cells. However, under persistent selection pressure, it may also promote the enrichment of weakly immunogenic clones and facilitate immune escape. This dual role is commonly summarized as the “elimination-equilibrium-escape” process of cancer immunoediting [1]. As the tumor microenvironment, tumor-promoting inflammation, and immune escape have been incorporated into the expanded framework of cancer hallmarks, quantitatively revealing the nonlinear feedback between tumor growth and immune regulation has become an important problem in mathematical oncology and tumor immunology [2].

Mathematical models can transform mechanisms such as antigen stimulation, immune-cell recruitment and maturation, effector-cell killing, tumor resource competition, and therapeutic intervention into tractable dynamical systems. Compared with statistical analyses that describe correlations at only a single time point, dynamical models can study long-term behaviors such as thresholds, steady states, bistability, dormancy, recurrence, and periodic oscillations within a unified framework. Parameter sensitivity and bifurcation analyses can also identify the key mechanisms that determine transitions between system states [3] [4]. In recent years, mathematical models have also been increasingly integrated with clinical data, virtual patients, and quantitative systems pharmacology. This integration provides a useful methodological basis for evaluating immunotherapy strategies and making individualized predictions [5].

Research on tumor-immune dynamics can be traced back to the immune-surveillance model developed by DeLisi and Rescigno. That model represented the interaction between lymphocytes and tumor cells as a predator-prey relationship. It showed that a low-dimensional nonlinear system can already generate distinct outcomes, including tumor clearance, persistent coexistence, and immune escape [6]. Later, Rescigno and DeLisi introduced a two-stage structure of immature and mature lymphocytes. This structure clearly separated immune-cell recruitment, maturation, and killing effects across different time scales [7]. Kuznetsov et al. combined experimental data, parameter estimation, and global bifurcation analysis. Their study revealed threshold effects, sneaking-through behavior, tumor dormancy, and recurrence-like oscillations in an immunogenic tumor model [8]. Kirschner and Panetta further discussed adoptive immunotherapy and long-term recurrence within a tumor-cell-effector-cell-IL-2 framework [9]. de Pillis et al. constructed and validated a detailed cell-mediated immune-response model by coupling NK cells, CD8⁺ T cells, and tumor cells [10].

Building on these studies, the research focus gradually expanded from the existence of equilibria to complex attractors and bifurcation mechanisms. Eftimie et al. systematically reviewed the structure, scales, and main dynamical conclusions of non-spatial tumor-immune models. They emphasized that low-dimensional ordinary differential equation models remain especially valuable for threshold identification and clear mechanistic interpretation [3]. Liu, Ruan, and Zhu proved the existence of stable periodic oscillations in a two-stage tumor-immune model. They interpreted these oscillations as recurrence-like behavior in which tumor burden and immune response repeatedly rise and fall [11]. Pang et al. discussed steady and oscillatory dynamics from the perspective of antitumor immune responses [12]. Li et al. further analyzed how intrinsic tumor growth, antigen stimulation, and immune-killing parameters affect equilibrium stability and Hopf bifurcation in a two-stage lymphocyte framework [13].

Current tumor-immune modeling is moving toward more detailed cellular states, stronger data constraints, and greater clinical translatability. For example, Lai et al. divided cytotoxic T cells according to their degree of exhaustion and revealed how T-cell functional decline affects tumor clearance, equilibrium, and escape outcomes [14]. Yao et al. combined an ordinary differential equation model of the tumor immune microenvironment with deep reinforcement learning to optimize patient-specific combination therapy with immune checkpoint inhibitors [15]. Arulraj et al. emphasized the important role of multi-omics data in the parameterization, calibration, and validation of quantitative systems pharmacology models [16]. These studies have promoted a shift from mechanistic explanation toward prediction and decision-making. However, they have also introduced challenges such as high parameter dimensionality, limited identifiability, and less transparent analytical structure.

Therefore, high-dimensional and data-driven models cannot replace clear and well-structured low-dimensional dynamical analyses. For models that include immune recruitment, cell maturation, saturating killing, and Logistic tumor growth, a complete theoretical study should establish the biological feasibility of solutions and the stability of equilibria. It should also identify the state-exchange mechanism at boundary equilibria and verify the transversality of Hopf bifurcation. In addition, center manifold and normal form theories should be used to determine the direction of periodic-solution branches and their orbital stability [17] [18]. In particular, cross-validation between the first Lyapunov coefficient and numerical critical values can prevent the oscillation type from being judged only from time series or phase portraits. This procedure improves the reproducibility and theoretical reliability of the model conclusions.

Based on the two-stage lymphocyte tumor-immune model proposed by Li et al. [13], this study retains the transformation of immature T lymphocytes into mature T lymphocytes, Logistic tumor growth, and the antigen-stimulation term. The killing effect of mature T lymphocytes on tumor cells is rewritten as a Holling type-II saturating functional response to describe the biologically realistic saturation of the killing efficiency per effector cell under a high tumor burden. The model is given as follows:

$\{\begin{array}{l}\frac{\text{d}{L}_1}{\text{d}t}=b-\left({\delta }_0+\mu \right)L_1+\epsilon T,\\ \frac{\text{d}{L}_2}{\text{d}t}=\mu L_1-\delta L_2,\\ \frac{\text{d}T}{\text{d}t}=rT\left(1-\frac{T}{K}\right)-s\frac{L_{2}T}{a+T}.\end{array}$ (1.1)

where b is the production rate of immature T lymphocytes in the absence of tumor cells; ${\delta }_0$ and $\delta$ are the natural death-rate coefficients of immature and mature T lymphocytes, respectively; $\mu$ is the rate coefficient for the transformation of immature T lymphocytes into mature T lymphocytes; $\epsilon$ is the recruitment-rate coefficient of immature T lymphocytes induced by tumor-antigen stimulation; r is the maximum tumor-cell growth-rate coefficient; K is the carrying capacity of the tumor microenvironment, namely the maximum number of tumor cells that the microenvironment can support; and s is the killing-rate coefficient of mature T lymphocytes against tumor cells. The parameter $a$ is the half-saturation constant, and all parameters are positive.

For convenient analysis, we introduce the following nondimensional transformation for system (1.1):

$x=l_1{L}_1,y=l_2{L}_2,z=l_{0}T,\tau =t_{0}t$

where

$l_1=\frac{{\delta }_0+\mu }{b},l_2=\frac{{\left({\delta }_0+\mu \right)}^2}{\mu b},l_0=\frac{1}{K},t_0={\delta }_0+\mu$

and the nondimensional parameters are

$\alpha =\frac{\epsilon K}{b},\beta =\frac{\delta }{{\delta }_0+\mu },\gamma =\frac{r}{{\delta }_0+\mu },m=\frac{s\mu b}{{\left({\delta }_0+\mu \right)}^3},\eta =\frac{a}{K}$

For convenience, we continue to use the notation

$\tau =t$

Then system (1.1) can be rewritten as

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=1-x+\alpha z,\\ \frac{\text{d}y}{\text{d}t}=x-\beta y,\\ \frac{\text{d}z}{\text{d}t}=\gamma z\left(1-z\right)-m\frac{yz}{\eta +z}.\end{array}$ (1.2)

This study mainly investigates the qualitative properties and bifurcation structure of system (1.2). The main contributions are summarized as follows. First, the positivity and boundedness of solutions with nonnegative initial values are established, and the existence conditions and local stability criteria for the semi-trivial and positive equilibria are systematically derived. Second, the Sotomayor theorem is used to analyze the transcritical bifurcation at the boundary equilibrium and to clarify the dynamical mechanism of equilibrium exchange near the tumor-invasion threshold. Third, the nondimensional parameter $\gamma$ , which corresponds to the maximum tumor-cell growth rate, is selected as the bifurcation parameter. The pure-imaginary-root condition and the transversality criterion for Hopf bifurcation at a positive equilibrium are obtained. Center manifold and normal form theories are then used to construct the first Lyapunov coefficient and determine the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions. Finally, MATLAB simulations, equilibrium calculations, and eigenvalue verification are used to cross-check the theoretical results. While retaining a clear and interpretable low-dimensional structure, the analysis connects three important dynamical states: the tumor-invasion threshold, stable coexistence, and recurrence-like periodic oscillation.

2. Positivity and Boundedness of the Model Solutions

Theorem 2.1: If the initial conditions satisfy $x\left(0\right)>0$ , $y\left(0\right)>0$ , $z\left(0\right)>0$ , then for every $t>0$ , we have $x\left(t\right)>0$ , $y\left(t\right)>0$ , $z\left(t\right)>0$ .

Proof: The vector field is continuous on the boundary of the region $x>0$ , $y>0$ , $z>0$ and satisfies the Lipschitz condition. By the existence and uniqueness theorem, the solution $\left(x,y,z\right)$ of system (1.2) exists uniquely on $\left(0,\tau \right)$ . Moreover,

$\begin{array}{l}x\left(t\right)=x\left(0\right){\text{e}}^{-t}+{\text{e}}^{-t}\left[{\displaystyle {\int }_0^t\left[1+\alpha z\left(s\right)\right]{\text{e}}^s\text{d}s}\right],\\ y\left(t\right)=y\left(0\right){\text{e}}^{-\beta t}+{\text{e}}^{-\beta t}\left[{\displaystyle {\int }_0^{t}x\left(s\right){\text{e}}^{\beta s}\text{d}s}\right],\\ z\left(t\right)=z\left(0\right)\text{exp}\left[{\displaystyle {\int }_0^t\left[\gamma \left(1-z\left(s\right)\right)-m\frac{y\left(s\right)}{\eta +z\left(s\right)}\right]\text{d}s}\right]\end{array}$

Since $x\left(0\right)>0$ , $y\left(0\right)>0$ , $z\left(0\right)>0$ , it follows that $x\left(t\right)>0$ , $y\left(t\right)>0$ , $z\left(t\right)>0$ for every $t>0$ .

Theorem 2.2: For every nonnegative initial value, the corresponding solution of system (1.2) is ultimately bounded. More precisely, the rectangular set

$\Omega =\left\{\left(x,y,z\right)\in {ℝ}_{+}^3:0\le z\le 2,0\le x\le 2\left(1+2\alpha \right),0\le y\le \frac{4\left(1+2\alpha \right)}{\beta }\right\}$

is an absorbing region: every solution with a nonnegative initial value enters Ω after a finite time and remains in a bounded subset of ${ℝ}_{+}^3$ .

Proof: From the third equation of system (1.2),

$\frac{\text{d}z}{\text{d}t}=\gamma z\left(1-z\right)-\frac{myz}{\eta +z}\le \gamma z\left(1-z\right).$ (2.1)

Let $u^{\prime }=\gamma u\left(1-u\right)$ with $u\left(0\right)=z\left(0\right)$ . The scalar comparison principle gives $0\le z\left(t\right)\le u\left(t\right)$ , and $u\left(t\right)\to 1$ . Hence there exists $t_1>0$ such that $0\le z\left(t\right)\le 2$ for all $t\ge t_1$ .

For $t\ge t_1$ , the first equation satisfies

$\frac{\text{d}x}{\text{d}t}\le 1+2\alpha -x.$ (2.2)

Comparison with $v^{\prime }=1+2\alpha -v$ shows that there exists $t_2\ge t_1$ such that $0\le x\left(t\right)\le 2\left(1+2\alpha \right)$ for all $t\ge t_2$ . Consequently, for $t\ge t_2$ ,

$\frac{\text{d}y}{\text{d}t}\le 2\left(1+2\alpha \right)-\beta y.$ (2.3)

A final comparison with the corresponding linear equation yields a time $t_3\ge t_2$ such that $0\le y\left(t\right)\le 4\left(1+2\alpha \right)/\beta$ for all $t\ge t_3$ . Thus, each nonnegative solution eventually enters the explicitly defined absorbing region Ω, proving ultimate boundedness.

The positivity and boundedness results show that the system cannot produce negative cell populations or finite-time blow-up. Therefore, equilibrium and bifurcation analyses can be carried out within the biologically feasible region.

3. Existence and Stability of Equilibria

3.1. Existence of Equilibria

The equilibria of system (1.2) satisfy

$\{\begin{array}{l}1-x+\alpha z=0,\\ x-\beta y=0,\\ \gamma z\left(1-z\right)-m\frac{yz}{\eta +z}=0.\end{array}$ (3.1)

The boundary equation $z=0$ gives the semi-trivial equilibrium

$E_0=\left(1,\frac{1}{\beta },0\right).$

For a positive equilibrium, $z^{\text{*}}>0$ , the first two equations give

$x^{\text{*}}=1+\alpha z^{\text{*}},y^{\text{*}}=\frac{1+\alpha z^{\text{*}}}{\beta }.$

After division of the third equation by $z^{\text{*}}>0$ and substitution of these expressions, one obtains

$\gamma \left(1-z^{\text{*}}\right)-\frac{m\left(1+\alpha z^{\text{*}}\right)}{\beta \left(\eta +z^{\text{*}}\right)}=0.$

Multiplying by $\beta \left(\eta +z^{\text{*}}\right)$ and collecting powers of $z^{\text{*}}$ yields the scalar quadratic

$f\left(z\right)=a_1{z}^2+a_{2}z+a_3=0,$ (3.2)

where

$a_1=\gamma \beta ,$

$a_2=m\alpha +\gamma \beta \left(\eta -1\right),$

$a_3=m-\gamma \beta \eta .$

when $\text{Δ}\ge 0$ , the solutions of Equation (3.2) are

$z_1=\frac{-a_2+\sqrt{\text{Δ}}}{2a_1},$

$z_2=\frac{-a_2-\sqrt{\text{Δ}}}{2a_1}.$

where

$\begin{array}{l}\text{Δ}=a_2^2-4a_1{a}_3={\left[m\alpha +\gamma \beta \left(\eta -1\right)\right]}^2-4\gamma \beta \left(m-\gamma \beta \eta \right)\hfill \end{array}$

A root is biologically admissible only when $0<z<1$ . For every such root, $x=1+\alpha z>0$ and $y=\left(1+\alpha z\right)/\beta >0$ , so it determines a feasible positive equilibrium. The upper bound $z<1$ also follows directly from the equilibrium identity $\gamma \left(1-z\right)=m\left(1+\alpha z\right)/\left[\beta \left(\eta +z\right)\right]>0$ . Therefore, every positive equilibrium satisfies $0<z<1$ .

Theorem 3.1. System (1.2) has the semi-trivial equilibrium $E_0=\left(1,\frac{1}{\beta },0\right)$ . For a positive solution $E^{\text{*}}=\left(x^{\text{*}},y^{\text{*}},z^{\text{*}}\right)$ , the following conclusions hold:

1) If any one of the following conditions holds, model (1.2) has a unique positive equilibrium $E^{\text{*}}=\left(x^{\text{*}},y^{\text{*}},z^{\text{*}}\right)$ .

(H₁) $m<\gamma \beta \eta$ .

(H₂) $m=\gamma \beta \eta$ , $\eta <\frac{1}{1+\alpha }$ .

(H₃) $m>\gamma \beta \eta$ , $m\alpha <\gamma \beta \left(1-\eta \right)$ , and ${\left[m\alpha +\gamma \beta \left(\eta -1\right)\right]}^2-4\gamma \beta \left(m-\gamma \beta \eta \right)=0$ .

2) If

(H₄) $m>\gamma \beta \eta$ , $m\alpha <\gamma \beta \left(1-\eta \right)$ , and ${\left[m\alpha +\gamma \beta \left(\eta -1\right)\right]}^2-4\gamma \beta \left(m-\gamma \beta \eta \right)>0$ .

then Equation (3.2) has two equilibria:

$E_1=\left(1+\alpha z_1,\frac{1+\alpha z_1}{\beta },z_1\right),E_2=\left(1+\alpha z_2,\frac{1+\alpha z_2}{\beta },z_2\right)$

3.2. Local Stability of Equilibria

Let

$\begin{array}{l}{f}_1\left(x,y,z\right)=1-x+\alpha z,\\ f_2\left(x,y,z\right)=x-\beta y,\\ f_3\left(x,y,z\right)=\gamma z\left(1-z\right)-m\frac{yz}{\eta +z}.\end{array}$ (3.3)

The Jacobian matrix of the model at the equilibrium $E^{\text{*}}$ is

$J\left(x,y,z\right)=\left(\begin{array}{ccc}-1& 0& \alpha \\ 1& -\beta & 0\\ 0& -\frac{mz}{\eta +z}& \gamma \left(1-2z\right)-\frac{m\eta y}{{\left(\eta +z\right)}^2}\end{array}\right).$

Here, $s=-\frac{mz}{\eta +z}$ , $s_0=\gamma \left(1-2z\right)-\frac{m\eta y}{{\left(\eta +z\right)}^2}$ .

The stability of an equilibrium is determined below by calculating the eigenvalues of the Jacobian matrix.

Theorem 3.2: If $m>\eta \gamma \beta$ holds, then the semi-trivial equilibrium $E_0$ is locally asymptotically stable. If $m<\gamma \beta \eta$ , then $E_0$ is unstable.

Proof: The Jacobian matrix of the model at $E_0$ is

$J\left(E_0\right)=\left(\begin{array}{ccc}-1& 0& \alpha \\ 1& -\beta & 0\\ 0& 0& \gamma -\frac{m}{\beta \eta }\end{array}\right).$

The characteristic polynomial of $J\left(E_0\right)$ is $f\left(\lambda \right)=\left(\lambda +1\right)\left(\lambda +\beta \right)\left(\lambda -\left(\gamma -\frac{m}{\beta \eta }\right)\right)=0$ .

The eigenvalues of $J\left(E_0\right)$ are

${\lambda }_1=-1,{\lambda }_2=-\beta ,{\lambda }_3=\gamma -\frac{m}{\beta \eta }$

If ${\lambda }_3<0$ , i.e., $m>\gamma \beta \eta$ , then $f\left(\lambda \right)$ has three negative real roots, and $E_0$ is locally asymptotically stable. If ${\lambda }_3>0$ , i.e., $m<\gamma \beta \eta$ , then at least one root is positive, and therefore $E_0$ is unstable.

Biologically, $m$ represents the threshold at which the basal mature immune level can suppress tumor invasion. When immune-cell killing is sufficiently strong, a small tumor perturbation cannot invade. When immune killing is insufficient, the tumor can enter the system and establish coexistence or oscillation with immune cells.

Theorem 3.3: Suppose that (H₁) holds. So that system (1.2) has a unique positive equilibrium $E^{\text{*}}$ . Then $E^{\text{*}}$ is locally asymptotically stable if and only if

(H₅) $b_1>0,b_2>0,b_3>0,b_1{b}_2-b_3>0.$

Proof: The Jacobian matrix of the system (1.2) at the equilibrium $E^{\text{*}}$ is

$J\left(E^{\text{*}}\right)=\left(\begin{array}{ccc}-1& 0& \alpha \\ 1& -\beta & 0\\ 0& -a_{32}& a_{33}\end{array}\right).$

where

$a_{32}=m\frac{z^{\text{*}}}{\eta +z^{\text{*}}},a_{33}=\gamma \left(1-2z^{\text{*}}\right)-\frac{m\eta y^{\text{*}}}{{\left(\eta +z^{\text{*}}\right)}^2}.$

the characteristic polynomial of the Jacobian at $E^{\text{*}}$ is

${\lambda }^3+b_1{\lambda }^2+b_2\lambda +b_3=0$ (3.4)

$\begin{array}{l}{b}_1=\frac{2\gamma {\left(z^{\text{*}}\right)}^2+\left(1+\beta +\eta \gamma -\gamma \right)z^{\text{*}}+\eta \left(1+\beta \right)}{\eta +z^{\text{*}}},\\ b_2=\frac{2\gamma \left(1+\beta \right){\left(z^{\text{*}}\right)}^2+\left(\beta -\gamma +\eta \gamma +\eta \gamma \beta \right)z^{\text{*}}+\beta \eta }{\eta +z^{\text{*}}},\\ b_3=\frac{2\gamma \beta {\left(z^{\text{*}}\right)}^2+\left(\alpha m+\beta \eta \gamma -\beta \gamma \right)z^{\text{*}}}{\eta +z^{\text{*}}}.\end{array}$

For a monic cubic polynomial, the Routh-Hurwitz criterion is necessary and sufficient. Therefore, all eigenvalues have negative real parts if and only if the four strict inequalities in (H₅) hold, proving the local asymptotic stability statement.

If the strict inequalities in (H₅) fail and $p$ has no root on the imaginary axis, then at least one eigenvalue has a positive real part and $E^{\text{*}}$ is unstable. The boundary cases must be separated from this genuinely unstable case. If $b_3=0$ , then $\lambda =0$ is an eigenvalue. If

$b_1>0,b_2>0,b_3>0,b_1{b}_2=b_3,$

then

$p\left(\lambda \right)=\left(\lambda +b_1\right)\left({\lambda }^2+b_2\right),$

so the eigenvalues are $-b_1$ and $\pm i\sqrt{b_2}$ . This is the Hopf boundary rather than a genuinely unstable case.

Several key parameters can be varied over suitable ranges to analyze their effects on equilibrium stability and system dynamics. Consider $\gamma$ as an example. As $\gamma$ gradually increases, the stability of the semi-trivial equilibrium $E_0$ may change. A previously stable $E_0$ may lose stability, while the stability region of the positive equilibrium $E^{\text{*}}$ may shrink or expand. This occurs because $\gamma$ is the maximum tumor-cell growth rate, and its variation disturbs the overall immune balance. Similarly, variation in $m$ , which represents the killing rate of mature T lymphocytes against tumor cells, may also change equilibrium stability and affect the occurrence and threshold of Hopf bifurcation. This sensitivity analysis provides a clearer understanding of the importance of each parameter and the response of the system to changes in different biological processes.

4. Bifurcation Analysis of the System

4.1. Transcritical Bifurcation

Select $m$ as the bifurcation parameter. If $m=m_c=\gamma \beta \eta$ and $\eta \ne \frac{1}{1+\alpha }$ , then system (1.2) undergoes a transcritical bifurcation near the equilibrium $E_0=\left(1,\frac{1}{\beta },0\right)$ .

Proof: Let

$F\left(x,y,z,m\right)=\left(\begin{array}{c}{f}_1\left(x,y,z\right)\\ f_2\left(x,y,z\right)\\ f_3\left(x,y,z,m\right)\end{array}\right).$

The Jacobian matrix of the system at $E_0$ is

$J\left(E_0\right)=\left(\begin{array}{ccc}-1& 0& \alpha \\ 1& -\beta & 0\\ 0& 0& \gamma -\frac{m}{\beta \eta }\end{array}\right).$

When $m=m_c=\gamma \beta \eta$ , we have $\gamma -\frac{m}{\beta \eta }=0$ .

Therefore,

$J\left(E_0;m_c\right)=\left(\begin{array}{ccc}-1& 0& \alpha \\ 1& -\beta & 0\\ 0& 0& 0\end{array}\right).$

The eigenvalues are ${\lambda }_1=-1,{\lambda }_2=-\beta ,{\lambda }_3=0$ .

For the zero eigenvalue, the corresponding eigenvectors of $J\left(E_0;m_c\right)$ and $J^{\text{T}}\left(E_0;m_c\right)$ are

$V=\left(\begin{array}{c}\alpha \\ \alpha /\beta \\ 1\end{array}\right),W=\left(\begin{array}{c}0\\ 0\\ 1\end{array}\right),W^{\text{T}}V=1.$

We now verify the nondegeneracy conditions for the transcritical bifurcation in the Sotomayor theorem.

$F_m\left(E_0,m_c\right)={\left(\begin{array}{l}\frac{\partial f_1}{\partial m}\hfill \\ \frac{\partial f_2}{\partial m}\hfill \\ \frac{\partial f_3}{\partial m}\hfill \end{array}\right)|}_{\left(E_0,m_c\right)}=\left(\begin{array}{l}0\hfill \\ 0\hfill \\ 0\hfill \end{array}\right).$

$D{F}_m\left(E_0,m_c\right)V={\left(\begin{array}{ccc}\frac{\partial f_{1m}}{\partial x}& \frac{\partial f_{1m}}{\partial y}& \frac{\partial f_{1m}}{\partial z}\\ \frac{\partial f_{2m}}{\partial x}& \frac{\partial f_{2m}}{\partial y}& \frac{\partial f_{2m}}{\partial z}\\ \frac{\partial f_{3m}}{\partial x}& \frac{\partial f_{3m}}{\partial y}& \frac{\partial f_{3m}}{\partial z}\end{array}\right)\left(\begin{array}{l}{v}_1\hfill \\ v_2\hfill \\ v_3\hfill \end{array}\right)|}_{\left(E_0,m_c\right)}=\left(\begin{array}{c}0\\ 0\\ -\frac{1}{\beta \eta }\end{array}\right).$

$\begin{array}{l}{D}^2{F}_m\left(E_0,m_c\right)\left(V,V\right)\\ ={\left(\begin{array}{c}\frac{{\partial }^2{f}_1}{\partial x^2}{v}_1^2+\frac{{\partial }^2{f}_1}{\partial y^2}{v}_2^2+\frac{{\partial }^2{f}_1}{\partial z^2}{v}_3^2+2\frac{{\partial }^2{f}_1}{\partial x\partial y}{v}_1{v}_2+2\frac{{\partial }^2{f}_1}{\partial x\partial z}{v}_1{v}_3+2\frac{{\partial }^2{f}_1}{\partial y\partial z}{v}_2{v}_3\\ \frac{{\partial }^2{f}_2}{\partial x^2}{v}_1^2+\frac{{\partial }^2{f}_2}{\partial y^2}{v}_2^2+\frac{{\partial }^2{f}_2}{\partial z^2}{v}_3^2+2\frac{{\partial }^2{f}_2}{\partial x\partial y}{v}_1{v}_2+2\frac{{\partial }^2{f}_2}{\partial x\partial z}{v}_1{v}_3+2\frac{{\partial }^2{f}_2}{\partial y\partial z}{v}_2{v}_3\\ \frac{{\partial }^2{f}_3}{\partial x^2}{v}_1^2+\frac{{\partial }^2{f}_3}{\partial y^2}{v}_2^2+\frac{{\partial }^2{f}_3}{\partial z^2}{v}_3^2+2\frac{{\partial }^2{f}_3}{\partial x\partial y}{v}_1{v}_2+2\frac{{\partial }^2{f}_3}{\partial x\partial z}{v}_1{v}_3+2\frac{{\partial }^2{f}_3}{\partial y\partial z}{v}_2{v}_3\end{array}\right)|}_{\left(E_0,m_c\right)}\\ =\left(\begin{array}{c}0\\ 0\\ \frac{2\gamma }{\eta }\left[1-\eta \left(1+\alpha \right)\right]\end{array}\right).\end{array}$

We next verify the transversality conditions.

$\begin{array}{l}{W}^{\text{T}}{F}_m\left(E_0,m_c\right)=0,\\ W^{\text{T}}\left[D{F}_m\left(E_0,m_c\right)V\right]=-\frac{1}{\beta \eta }\ne 0,\\ W^{\text{T}}\left[D^2{F}_m\left(E_0,m_c\right)\left(V,V\right)\right]=\frac{2\gamma }{\eta }\left[1-\eta \left(1+\alpha \right)\right]\ne 0.\end{array}$

Therefore, by the Sotomayor theorem, system (1.2) undergoes a transcritical bifurcation at the boundary equilibrium $E_0$ .

4.2. Hopf Bifurcation

When the bifurcation value changes, the stability of the model changes abruptly, and a limit cycle “emerges” around an equilibrium. When an equilibrium changes from stable to unstable, or from unstable to stable, the topological structure of the system solutions changes in a small neighborhood of the bifurcation value, and a periodic solution is generated. We therefore use the Poincare-Andronov-Hopf bifurcation theorem to discuss the Hopf bifurcation of system (1.2).

Lemma 4.1: Let $\Omega \subset {ℝ}^3$ be an open set containing $O\left(x_1,x_2,x_3\right)$ , and let $S\subset ℝ$ be an open set containing 0. Suppose that $f:\Omega \times S\to R^3$ is analytic and that $f\left(0,\rho \right)=0$ for every $\rho \in S$ . Assume that the variational matrix $Df\left(0,\rho \right)$ has one real eigenvalue $\gamma \left(\rho \right)$ and a pair of complex-conjugate eigenvalues $\alpha \left(\rho \right)\pm i\beta \left(\rho \right)$ . At $\rho =0$ , suppose that

$\gamma \left(0\right)<0$ , $\alpha \left(0\right)=0$ , $\beta \left(0\right)>0$ ,

and that the eigenvalues cross the imaginary axis with a nonzero speed, namely,

$\frac{\text{d}\alpha \left(\rho \right)}{\text{d}\rho }\ne 0$ .

Then the differential system

$X=f\left(X,\rho \right)$

undergoes a Hopf bifurcation at the equilibrium O when $\rho =0$ .

Definition 4.1 Assume that (H₁) holds, and define the Routh-Hurwitz discriminant function

$H\left(\gamma \right)=b_1\left(\gamma \right)b_2\left(\gamma \right)-b_3\left(\gamma \right).$ (4.1)

If there exists $\gamma ={\gamma }_H$ such that

$b_1\left({\gamma }_H\right)>0,b_2\left({\gamma }_H\right)>0,b_3\left({\gamma }_H\right)>0,H\left({\gamma }_H\right)=0,$

From the characteristic equation

${\lambda }^3+b_1{\lambda }^2+b_2\lambda +b_3=0$ (4.2)

Substituting the purely imaginary root $\lambda =\xi +i\omega$ into the characteristic equation gives

${\left(\xi +i\omega \right)}^3+b_1{\left(\xi +i\omega \right)}^2+b_2\left(\xi +i\omega \right)+b_3=0$ (4.3)

Choose $\gamma ={\gamma }_c$ such that $\alpha \left({\gamma }_c\right)=0$ . Then the characteristic equation becomes

${\left(i\omega \right)}^3+b_1{\left(i\omega \right)}^2+b_2\left(i\omega \right)+b_3=0$ (4.4)

Separating the real and imaginary parts gives

$\{\begin{array}{l}-b_1{\omega }^2+b_3=0,\\ -{\omega }^3+b_2\omega =0.\end{array}$ (4.5)

From the second equation in (4.5), we obtain

$\omega ={\omega }_0=\sqrt{b_2}$ , $b_3=b_1{b}_2$

${\lambda }_{1,2}=\pm i\sqrt{b_2}\ne 0$

Since $b_2>0$ , ${\omega }_0>0$ is a positive real number.

The sum of the three roots is ${\lambda }_1+{\lambda }_2+{\lambda }_3=i\omega +\left(-i\omega \right)+{\lambda }_3=-b_1$ .

Therefore, ${\lambda }_3=-b_1<0$ .

This shows that the system has a pair of purely imaginary conjugate eigenvalues and one negative real eigenvalue at the positive equilibrium.

Theorem 4.1: Select $\gamma$ as the bifurcation parameter. Suppose that system (1.2) has a unique positive equilibrium $E^{\text{*}}$ when $m<\gamma \beta \eta$ . If there exists $\gamma ={\gamma }_H$ satisfying

$H\left({\gamma }_H\right)=b_1\left({\gamma }_H\right)b_2\left({\gamma }_H\right)-b_3\left({\gamma }_H\right)=0.$

and

${\frac{\text{d}\left(b_1\left(\gamma \right)b_2\left(\gamma \right)-b_3\left(\gamma \right)\right)}{\text{d}\gamma }|}_{\gamma ={\gamma }_H}\ne 0$ (4.6)

then system (1.2) undergoes a Hopf bifurcation at $E^{\text{*}}\left({\gamma }_H\right)$ .

Proof: When $\gamma ={\gamma }_H$ , the linearized matrix has a pair of purely imaginary conjugate eigenvalues $\pm i{\omega }_H$ , where

${\omega }_H=\sqrt{b_2\left({\gamma }_H\right)}>0,$

The other eigenvalue is $-b_1\left({\gamma }_H\right)$ . Let

$\lambda \left(\gamma \right)=\xi \left(\gamma \right)+i\omega \left(\gamma \right)$

be the eigenvalue branch passing through $i{\omega }_H$ . Write the characteristic polynomial as

$P\left(\lambda ,\gamma \right)={\lambda }^3+b_1\left(\gamma \right){\lambda }^2+b_2\left(\gamma \right)\lambda +b_3\left(\gamma \right).$

Differentiating $P\left(\lambda \left(\gamma \right),\gamma \right)=0$ with respect to $\gamma$ gives

$\left(3{\lambda }^2+2b_1\lambda +b_2\right)\frac{\text{d}\lambda }{\text{d}\gamma }+\frac{\text{d}{b}_1}{\text{d}\gamma }{\lambda }^2+\frac{\text{d}{b}_2}{\text{d}\gamma }\lambda +\frac{\text{d}{b}_3}{\text{d}\gamma }=0$

Thus,

$\frac{\text{d}\lambda }{\text{d}\gamma }=-\frac{P_{\gamma }\left(\lambda ,\gamma \right)}{P_{\lambda }\left(\lambda ,\gamma \right)}=-\frac{\frac{\text{d}{b}_1}{\text{d}\gamma }{\lambda }^2+\frac{\text{d}{b}_2}{\text{d}\gamma }\lambda +\frac{\text{d}{b}_3}{\text{d}\gamma }}{3{\lambda }^2+2b_1\lambda +b_2}$

At $\gamma ={\gamma }_H$ and $\lambda =i{\omega }_H$ , and $b_3\left({\gamma }_H\right)=b_1\left({\gamma }_H\right)b_2\left({\gamma }_H\right)$ , we have

$P_{\lambda }\left(i{\omega }_H,{\lambda }_H\right)=-2b_2+2b_1{\omega }_{H}i$

and

$P_{\gamma }\left(i{\omega }_H,{\gamma }_H\right)=-\frac{\text{d}{b}_1}{\text{d}\gamma }{b}_2+\frac{\text{d}{b}_3}{\text{d}\gamma }+i{\omega }_H\frac{\text{d}{b}_2}{\text{d}\gamma }$

Taking the real part further gives

$a_H:={\frac{\text{d}}{\text{d}\gamma }\mathrm{Re}\lambda \left(\gamma \right)|}_{\gamma ={\gamma }_H}=-\frac{H^{\prime }\left({\gamma }_H\right)}{2\left[b_1{\left({\gamma }_H\right)}^2+b_2\left({\gamma }_H\right)\right]}.\ne 0$ (4.7)

Since the denominator is positive, $H^{\prime }\left({\gamma }_H\right)\ne 0$ is equivalent to $a_H\ne 0$ . condition (4.7) implies that a pair of complex-conjugate eigenvalues crosses the imaginary axis with a nonzero speed. Thus, the transversality condition holds. By the Poincare-Andronov-Hopf bifurcation theorem, the system undergoes a Hopf bifurcation at $E^{\text{*}}\left({\gamma }_H\right)$ . The proof is complete.

5. Direction of the Hopf Bifurcation and Stability of Periodic Solutions

To determine the direction of the Hopf bifurcation and the stability of the bifurcating periodic solutions more precisely, we use center manifold and normal form theories to calculate the first Lyapunov coefficient. We first translate the positive equilibrium of system (1.2) to the origin.

Let

$x_1=x-x^{\text{*}},y_2=y-y^{\text{*}},z_3=z-z^{\text{*}},$

For convenience, we continue to use $x,y,z$ in place of $x_1,y_1,z_1$ . Then system (5.1) becomes

$\{\begin{array}{l}\frac{\text{d}x}{\text{d}t}=1-\left(x+x^{\text{*}}\right)+\alpha \left(z+z^{\text{*}}\right),\\ \frac{\text{d}y}{\text{d}t}=\left(x+x^{\text{*}}\right)-\beta \left(y+y^{\text{*}}\right),\\ \frac{\text{d}z}{\text{d}t}=\gamma \left(z+z^{\text{*}}\right)\left(1-z-z^{\text{*}}\right)-m\frac{\left(y+y^{\text{*}}\right)\left(z+z^{\text{*}}\right)}{\eta +z+z^{\text{*}}}.\end{array}$ (5.1)

Thus, the equilibrium $E^{\text{*}}=\left(x^{\text{*}},y^{\text{*}},z^{\text{*}}\right)$ of system (1.2) is shifted to the origin $\left(0,0,0\right)$ , and the system becomes

$\frac{\text{d}U}{\text{d}t}=J_{H}U+F\left(U\right)+\frac{1}{2}B\left(U,U\right)+\frac{1}{6}C\left(U,U,U\right)+\cdots$ (5.2)

where

$U={\left(x,y,z\right)}^{\text{T}},F\left(U\right)={\left(F_1,F_2,F_3\right)}^{\text{T}}$ ,

the components of $F\left(U\right)$ are

$F_1=0,$

$F_2=0,$

$F_3=-\gamma z^2-\frac{m\eta }{{\left(\eta +z^{\text{*}}\right)}^2}yz+\frac{my\eta }{{\left(\eta +z^{\text{*}}\right)}^3}{z}^2+\cdots$

When $\gamma ={\gamma }_H$ , the eigenvalues at $E^{\text{*}}=\left(x^{\text{*}},y^{\text{*}},z^{\text{*}}\right)$ are ${\lambda }_{1,2}=\pm i{\omega }_H\left({\omega }_H=\sqrt{b_2}\right)$ and ${\lambda }_3=-b_1$ . Let ${\lambda }_3\ne 0$ . Choose the right and left eigenvectors $p,q$ , and use the Hermitian inner product $\langle p,q\rangle ={\overline{p}}^{\text{T}}q$, such that

$J_{H}q=i{\omega }_{H}q,J_H^{\text{T}}p=-i{\omega }_{H}p,\langle p,q\rangle =1,$

$q^{\left(0\right)}=\left(\begin{array}{c}\alpha \left(\beta +i{\omega }_H\right)\\ \alpha \\ \left(1+i{\omega }_H\right)\left(\beta +i{\omega }_H\right)\end{array}\right).$

This choice automatically satisfies the first two rows of the eigenvector equations. The third row follows from the characteristic equation

$\text{det}\left(i{\omega }_{H}I-J_H\right)=0$

After normalization, $\langle p,q\rangle =1$ .

At $\left(y^{\text{*}},z^{\text{*}};{\gamma }_H\right)$ , the function $f_3\left(y,z;\gamma \right)=\gamma z\left(1-z\right)-\frac{myz}{\eta +z}$ satisfies

$\begin{array}{l}{f}_{3,yz}=-\frac{m\eta }{{\left(\eta +z^{\text{*}}\right)}^2},f_{3,zz}=-2{\gamma }_H+\frac{2my\eta }{{\left(\eta +z^{\text{*}}\right)}^3},\\ f_{3,yzz}=\frac{2m\eta }{{\left(\eta +z^{\text{*}}\right)}^3},f_{3,zzz}=-\frac{6my\eta }{{\left(\eta +z^{\text{*}}\right)}^4},\end{array}$

Therefore, for arbitrary $x={\left(x_1,x_2,x_3\right)}^{\text{T}}$ , $y={\left(y_1,y_2,y_3\right)}^{\text{T}}$ , and $z={\left(z_1,z_2,z_3\right)}^{\text{T}}$ , we have

$B\left(u,v\right)=\left(\begin{array}{c}0\\ 0\\ f_{3,yz}\left(u_2{v}_3+u_3{v}_2\right)+f_{3,zz}{u}_3{v}_3\end{array}\right)$ (5.3)

$C\left(u,v,w\right)=\left(\begin{array}{c}0\\ 0\\ f_{3,yzz}\left(u_2{v}_3{w}_3+u_3{v}_2{w}_3+u_3{v}_3{w}_2\right)+f_{3,zzz}{u}_3{v}_3{w}_3\end{array}\right)$ (5.4)

The formula for calculating the first Lyapunov coefficient is as follows:

$\begin{array}{c}{l}_1\left(0\right)=\frac{1}{2{\omega }_H}\mathrm{Re}\{\langle p,C\left(q,q,\overline{q}\right)\rangle -2\langle p,B\left(q,J_H^{-1}B\left(q,\overline{q}\right)\right)\rangle \\ +\langle p,B\left(\overline{q},{\left(2i{\omega }_{H}I-J_H\right)}^{-1}B\left(q,q\right)\right)\rangle \}.\end{array}$(5.5)

All vectors and matrices in (5.5) are uniquely determined by $E^{\text{*}}$ , ${\gamma }_H$ , and the model parameters. Therefore, they can be evaluated directly by symbolic substitution or numerical calculation.

Let $\lambda \left(\gamma \right)$ denote the eigenvalue branch near ${\gamma }_H$ satisfying

$\lambda \left({\gamma }_H\right)=i{\omega }_H$

and denote the transversality derivative by

$a_H={\frac{\text{d}}{\text{d}\gamma }\mathrm{Re}\lambda \left(\gamma \right)|}_{\gamma ={\gamma }_H}\ne 0,$

Define

${\mu }_2=-\frac{l_1\left(0\right)}{a_H},{\beta }_2=2l_1\left(0\right).$ (5.6)

From the normal form of the Hopf bifurcation, the amplitude of the bifurcating periodic solution satisfies

$r^2=\frac{\gamma -{\gamma }_H}{{\mu }_2}+o\left(\left|\gamma -{\gamma }_H\right|\right).$ (5.7)

Therefore, if ${\mu }_2>0$ , the periodic-solution branch appears on the side $\gamma >{\gamma }_H$ .

If ${\mu }_2<0$ , the periodic-solution branch appears on the side $\gamma <{\gamma }_H$ . The stability of the periodic solution is determined by the sign of ${\beta }_2$ .

Theorem 5.1 Assume that (H₁) and condition (H₅) hold. Then system (1.2) undergoes a Hopf bifurcation near $\gamma ={\gamma }_H$ , and the following conclusions hold:

1) If $l_1\left(0\right)<0$ , the bifurcating periodic solution is orbitally asymptotically stable, and the Hopf bifurcation is supercritical.

2) If $l_1\left(0\right)>0$ , the bifurcating periodic solution is unstable, and the Hopf bifurcation is subcritical.

This result is consistent with the stable limit cycle observed in the numerical simulations. It shows that when the negative feedback among tumor-cell growth, antigen stimulation, T-cell maturation, and immune killing has a sufficiently strong phase difference, the system can shift from stable coexistence to stable periodic oscillation. Biologically, this stable limit cycle represents a recurrence-like process in which tumor burden and immune-cell levels repeatedly rise and fall.

6. Numerical Simulations

MATLAB is used to perform numerical simulations of the system. To make the numerical results directly comparable with the preceding theoretical analysis, the nondimensional parameter $\gamma$ , corresponding to the maximum tumor-cell growth rate, is selected as the main bifurcation parameter. The remaining parameters are fixed as

$\alpha =5.5,\beta =1.5,m=0.75,\eta =0.25,$

This set is used as an illustrative nondimensional case study rather than as a patient-calibrated parameter set. This set is used because it preserves a biologically feasible positive equilibrium and produces a tumor-invasion threshold together with two Hopf crossings, allowing stable coexistence and recurrence-like periodic oscillations to be compared within a single sweep of $\gamma$ . Quantitative calibration and validation against experimental or clinical data are left for future work.

The initial conditions are

$\left(x\left(0\right),y\left(0\right),z\left(0\right)\right)=\left(1.9,0.8,0.16\right).$

An adaptive Runge-Kutta method is used to solve the differential equations. A sufficiently long transient is removed when plotting the bifurcation and amplitude diagrams. The positive equilibrium is obtained by numerically solving the algebraic equations, while local stability is verified jointly by the eigenvalues of the Jacobian matrix and the Routh-Hurwitz discriminant function. For this parameter set, the transcritical threshold at the boundary equilibrium and the two Hopf critical values are

${\gamma }_T=2.000000,{\gamma }_{H1}=2.631082,{\gamma }_{H2}=\mathrm{3.474665.}$

The eigenvalues and normal-form quantities at the two Hopf critical points are listed in Table 1. Numerical differentiation of the critical eigenvalue branch gives $a_{H1}=0.256683>0$ and $a_{H2}=-0.165608<0$ . Evaluation of (5.5) gives $l_1\left({\gamma }_{H1}\right)=-0.533949$ and $l_1\left({\gamma }_{H2}\right)=-0.056076$ . Because both first Lyapunov coefficients are negative, both Hopf points are supercritical under the standard convention, and the emerging periodic solutions are orbitally asymptotically stable. The branch-side coefficients are ${\mu }_{2,1}=2.080187>0$ and ${\mu }_{2,2}=-0.338606<0$ . Therefore, the stable periodic branch lies on the side $\gamma >{\gamma }_{H1}$ at the first crossing and on the side $\gamma <{\gamma }_{H2}$ at the second crossing. These two orientations delimit the stable oscillatory interval observed numerically around $\gamma =3.00$ .

Table 1. Numerical information at the Hopf critical points.

Critical point	$\gamma$	Positive equilibrium $E^{\text{*}}=\left(x^{\text{*}},y^{\text{*}},z^{\text{*}}\right)$	${\omega }_H$	$a_H$	$l_1\left(0\right)$	Type and cycle stability
${\gamma }_{H1}$	2.631082	(1.760760, 1.173840, 0.138320)	0.625239	0.256683	−0.533949	Supercritical; stable
${\gamma }_{H2}$	3.474665	(2.681178, 1.787452, 0.305669)	0.915099	−0.165608	−0.056076	Supercritical; stable

6.1. Phase Portraits and Time Evolution for Different Values of $\gamma$

We first examine how the system trajectories change for several representative values of $\gamma$ . Figure 1 shows the phase portrait and time-evolution curves for $\gamma =2.40$ . This parameter value lies to the left of the first Hopf critical value. The positive equilibrium is locally asymptotically stable, and the trajectory converges to a stable coexistence state after a brief damped oscillation. Biologically, this means that the tumor burden and immune-cell levels eventually reach a relatively stable dynamic balance.

When $\gamma =3.00$ , the parameter lies between the two Hopf critical values. The positive equilibrium loses local stability, and the system trajectory forms a stable closed orbit near the positive equilibrium. Figure 2 shows a clear limit cycle in the phase portrait, while the time-evolution curves display sustained periodic oscillations. This indicates a phase-lagged feedback among tumor growth, antigen stimulation, T-cell maturation, and immune killing.

When $\gamma =4.00$ , the parameter passes the second Hopf critical value, and the positive equilibrium becomes locally stable again. Figure 3 shows that although the system still undergoes an initial oscillatory adjustment, the trajectory eventually converges to a new positive equilibrium. Compared with the case $\gamma =2.40$ , the equilibrium tumor burden is higher. This result shows that a stronger intrinsic tumor-growth ability raises the final coexistence level of the system.

Figure 1. Phase portrait and time evolution of the system at $\gamma =2.40$ .

Figure 2. Stable periodic oscillation of the system at $\gamma =3.00$ .

Figure 3. Stable positive equilibrium of the system at $\gamma =4.00$ .

6.2. Comparison of Phase Portraits and Evolution Curves under Parameter Variation

To compare the effect of $\gamma$ on the system dynamics more clearly, Figure 4 presents the x-z phase portraits for $\gamma =2.40$ , $\gamma =3.00$ , and $\gamma =4.00$ side by side. The trajectories converge to equilibria when $\gamma =2.40$ and $\gamma =4.00$ , whereas the trajectory remains on a closed curve when $\gamma =3.00$ . This comparison shows that stable periodic oscillations occur only inside the Hopf interval.

Figure 5 further compares the time evolution of the mature T-lymphocyte variable $y\left(\tau \right)$ and the tumor variable $z\left(\tau \right)$ . For $\gamma =3.00$ , both curves maintain sustained oscillations, reflecting the long-term mutual restriction between immune killing and tumor growth. For $\gamma =2.40$ and $\gamma =4.00$ , the oscillations gradually decay, and the system eventually enters a stable coexistence state.

Figure 4. Comparison of the x-z phase portraits for different values of $\gamma$ .

Figure 5. Comparison of the evolution curves of $y\left(\tau \right)$ and $z\left(\tau \right)$ for different values of $\gamma$ .

In summary, the numerical simulations agree with the preceding theoretical analysis. When $\gamma$ lies outside the Hopf interval, the system approaches a stable positive equilibrium. When $\gamma$ lies between the two Hopf critical values, the system produces stable periodic oscillations. These oscillations provide a natural dynamical explanation for recurrence-like behavior in which tumor burden and immune-cell numbers repeatedly increase and decrease.

7. Conclusion

This study presents a detailed dynamical analysis of a two-stage tumor-immune model with T-lymphocyte action. The positivity and boundedness of solutions, the existence of equilibria, and their local stability are systematically discussed. By constructing the Jacobian matrix and applying the Routh-Hurwitz criterion, sufficient conditions for the local asymptotic stability of the positive equilibrium are obtained. The Sotomayor theorem is then used to analyze the transcritical bifurcation at the boundary equilibrium. In addition, Poincare-Andronov-Hopf bifurcation theory is applied to determine the conditions for Hopf bifurcation at a positive equilibrium. The system is further examined from the perspectives of Hopf bifurcation direction and periodic-solution stability, and the first Lyapunov coefficient is used to determine the stability of the bifurcating periodic orbit. The numerical simulations are consistent with the theoretical analysis. When the tumor-growth parameter $\gamma$ lies outside the Hopf interval, the system eventually approaches a stable positive equilibrium. When $\gamma$ lies between the two Hopf critical values, the system produces stable periodic oscillations. The phase portraits and time-evolution curves reveal a clear nonlinear coupling among tumor growth, antigen stimulation, T-cell maturation, and immune killing. This feedback mechanism can explain periodic fluctuations in tumor burden and immune-cell levels.

Biologically, a stable positive equilibrium represents long-term coexistence between tumor cells and immune cells. A stable limit cycle induced by a Hopf bifurcation represents an alternating oscillatory process between tumor recurrence and immune suppression. When immune killing is relatively strong, a small tumor perturbation cannot continue to expand. When the intrinsic tumor-growth ability increases or immune-killing efficiency becomes relatively weaker, the system may shift from stable coexistence to periodic oscillation, or even approach a stable state with a higher tumor burden. Therefore, immunotherapy should consider not only the instantaneous killing intensity, but also the duration and replenishment of effector cells, as well as the balance between immune-killing efficiency and tumor-growth ability.

Overall, the model retains a low dimension and strong analytical tractability. It provides a clear dynamical explanation of stable coexistence, recurrence-like oscillations, and immune control in tumor-immune interactions. Future studies may extend the present framework by introducing NK cells, dendritic cells, macrophages, regulatory T cells, cytokines, immune checkpoints, CAR-T therapy, treatment pulses, time delays, random perturbations, or spatial diffusion. Parameter estimation and model validation based on experimental data may further improve the ability of the model to represent realistic tumor-immune processes.

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.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Schreiber, R.D., Old, L.J. and Smyth, M.J. (2011) Cancer Immunoediting: Integrating Immunity’s Roles in Cancer Suppression and Promotion. Science, 331, 1565-1570.[CrossRef] [PubMed]
[2]	Hanahan, D. (2022) Hallmarks of Cancer: New Dimensions. Cancer Discovery, 12, 31-46.[CrossRef] [PubMed]
[3]	Eftimie, R., Bramson, J.L. and Earn, D.J.D. (2011) Interactions between the Immune System and Cancer: A Brief Review of Non-Spatial Mathematical Models. Bulletin of Mathematical Biology, 73, 2-32.
[4]	Murray, J.D. (2002) Mathematical Biology I: An Introduction. 3rd Edition, Springer.
[5]	Butner, J.D., Dogra, P., Chung, C., Pasqualini, R., Arap, W., Lowengrub, J., et al. (2022) Mathematical Modeling of Cancer Immunotherapy for Personalized Clinical Translation. Nature Computational Science, 2, 785-796.[CrossRef] [PubMed]
[6]	Delisi, C. and Rescigno, A. (1977) Immune surveillance and Neoplasia—I. A Minimal Mathematical Model. Bulletin of Mathematical Biology, 39, 201-221.[CrossRef]
[7]	Rescigno, A. and Delisi, C. (1977) Immune Surveillance and Neoplasia—II. A Two-Stage Mathematical Model. Bulletin of Mathematical Biology, 39, 487-497.[CrossRef]
[8]	Kuznetsov, V.A., Makalkin, I.A., Taylor, M.A. and Perelson, A.S. (1994) Nonlinear Dynamics of Immunogenic Tumors: Parameter Estimation and Global Bifurcation Analysis. Bulletin of Mathematical Biology, 56, 295-321.[CrossRef] [PubMed]
[9]	Kirschner, D. and Panetta, J.C. (1998) Modeling Immunotherapy of the Tumor-Immune Interaction. Journal of Mathematical Biology, 37, 235-252.[CrossRef] [PubMed]
[10]	de Pillis, L.G., Radunskaya, A.E. and Wiseman, C.L. (2005) A Validated Mathematical Model of Cell-Mediated Immune Response to Tumor Growth. Cancer Research, 65, 7950-7958.[CrossRef] [PubMed]
[11]	Liu, D., Ruan, S. and Zhu, D. (2012) Stable Periodic Oscillations in a Two-Stage Cancer Model of Tumor and Immune System Interactions. Mathematical Biosciences and Engineering, 9, 347-368.
[12]	Pang, L., Liu, S., Zhang, X. and Tian, T. (2020) Mathematical Modeling and Dynamic Analysis of Anti-Tumor Immune Response. Journal of Applied Mathematics and Computing, 62, 473-488.[CrossRef]
[13]	Li, J., Chen, Y., Guo, J., Wu, H., Xi, X. and Zhang, D. (2025) Dynamical Analysis of a Simple Tumor-Immune Model with Two-Stage Lymphocytes. Mathematical Methods in the Applied Sciences, 48, 10016-10027.[CrossRef]
[14]	Lai, N., Farman, A. and Byrne, H.M. (2025) The Impact of T-Cell Exhaustion Dynamics on Tumour-Immune Interactions and Tumour Growth. Bulletin of Mathematical Biology, 87, Article No. 61.[CrossRef] [PubMed]
[15]	Yao, Y., Chen, Y.F. and Zhang, Q. (2024) Optimized Patient-Specific Immune Checkpoint Inhibitor Therapies for Cancer Treatment Based on Tumor Immune Microenvironment Modeling. Briefings in Bioinformatics, 25, bbae547.[CrossRef] [PubMed]
[16]	Arulraj, T., Wang, H., Ippolito, A., et al. (2024) Leveraging Multi-Omics Data to Empower Quantitative Systems Pharmacology in Immuno-Oncology. Briefings in Bioinformatics, 25, bbae131.
[17]	Kuznetsov, Y.A. (2004) Elements of Applied Bifurcation Theory. 3rd Edition, Springer.
[18]	Perko, L. (2001) Differential Equations and Dynamical Systems. 3rd Edition, Springer.