Nutrient Availability and Immune Cytotoxic Strength Determine Tumor Control Versus Escape: A Spatiotemporal Reaction–Diffusion Analysis

1. Introduction

The failure of the immune system to eliminate a growing tumor — a phenomenon termed immune escape — is one of the central unsolved problems in cancer biology and a primary obstacle to durable immunotherapy responses [14,22]. In solid tumors, immune control is not determined solely by the intrinsic cytotoxic capacity of effector cells, but is profoundly shaped by the spatial and metabolic architecture of the tumor microenvironment (TME). Nutrient depletion, in particular, drives the formation of hypoxic and necrotic cores that physically exclude immune infiltration and metabolically exhaust cytotoxic T lymphocytes before they can reach the tumor interior [10,34]. Understanding the quantitative conditions — in terms of immune killing rate, recruitment efficiency, and nutrient consumption — under which immune pressure is sufficient to contain tumor growth is therefore of direct relevance to the design of immunotherapy regimens and the identification of patients likely to respond to immune checkpoint blockade.

Mathematical modeling has emerged as a powerful tool for investigating these threshold phenomena, because the nonlinear coupling between tumor proliferation, immune dynamics, and nutrient transport makes it difficult to isolate causal mechanisms experimentally [5,15]. A model, as defined by Ayeni [4], is an intellectual construct which represents reality and can be manipulated to predict the consequences of future actions. Spatially resolved reaction–diffusion models are particularly informative because they capture the spatial heterogeneity of solid tumors: the necrotic core, the proliferating rim, and the peripheral immune infiltrate, all of which are clinically observable via imaging [19,38]. The spatial tumor–immune modeling framework builds on Mambili-Mamboundou et al. [26], who studied the interaction of tumour infiltrating cytotoxic lymphocytes with solid tumours using coupled reaction–diffusion equations, and demonstrated the existence of a stable tumour dormant state through stability analysis and numerical simulation.

The present study makes three interconnected contributions. First, we develop and rigorously analyze a three-component spatiotemporal model coupling tumor cell density, cytotoxic immune cell density, and a diffusible nutrient field, establishing positivity, uniform boundedness, and global well-posedness via analytic semigroup theory [17,32]. Second, equilibrium and stability analysis identifies two qualitatively distinct long-term states — a chronically coexisting tumor–immune state and a tumor-dominated escape state — and characterizes the parameter thresholds separating them in terms of biologically interpretable dimensionless groups: the immune cytotoxic strength $\chi$, the tumor-driven immune recruitment ratio $\eta /\mu$, and the nutrient consumption rate $\theta$. Third, numerical simulations under two representative parameter regimes illustrate these contrasting biological outcomes with spatial resolution, including the formation of a necrotic core under nutrient depletion consistent with avascular solid tumor morphology. To integrate these simulations reliably, the multiscale stiffness arising from widely separated reaction and diffusion timescales is handled by a new class of high-order A-stable implicit integrators (SSHBBDF) constructed via off-grid collocation, whose fourth-order convergence is rigorously verified. The capacity to distinguish immune-controlled tumor suppression from sustained tumor escape has direct implications for therapeutic decision-making, where computational predictions of disease trajectories increasingly inform treatment planning and immunotherapy protocol design.

The remainder of this paper is organized as follows. Section 2 presents the dimensional model, nondimensionalization, and governing equations. Section 3 establishes the analytical properties of the system. Section 4 introduces the SSHBBDF framework and derives specific schemes. Section 5 provides the stability analysis. Section 6 describes the numerical implementation. Section 7 presents results and biological interpretation. Section 8 concludes.

2. Mathematical Formulation

2.1. Dimensional Model

We consider a one-dimensional tissue segment of physical length $L_a$, representing a vascularized region in which tumor cells, immune cells, and a metabolic nutrient evolve and interact. The spatiotemporal dynamics of tumor cell density $T_a(x_a,t_a)$, immune cell density $I_a(x_a,t_a)$, and nutrient concentration $N_a(x_a,t_a)$ are governed by the dimensional reaction–diffusion system

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial T_a}{\partial t_a}}& =D_{T,a}{\displaystyle \frac{{\partial }^2{T}_a}{\partial x_a^2}}+p_a{T}_a\left(1-{\displaystyle \frac{T_a}{K_a}}\right)-{\gamma }_a{T}_a{I}_a+{\alpha }_a{T}_a{N}_a,\hfill \end{array}$$

(1)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial I_a}{\partial t_a}}& =D_{I,a}{\displaystyle \frac{{\partial }^2{I}_a}{\partial x_a^2}}+{\eta }_a{T}_a-{\mu }_a{I}_a,\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill {\displaystyle \frac{\partial N_a}{\partial t_a}}& =D_{N,a}{\displaystyle \frac{{\partial }^2{N}_a}{\partial x_a^2}}-{\delta }_a{T}_a{N}_a,\hfill \end{array}$$

(3)

for $x_a\in (0,L_a)$ and $t_a>0$. All dimensional variables and parameters are listed in Table 1. Homogeneous Neumann boundary conditions,

$${\displaystyle \frac{\partial T_a}{\partial x_a}}={\displaystyle \frac{\partial I_a}{\partial x_a}}={\displaystyle \frac{\partial N_a}{\partial x_a}}=0\mathrm{at}{x}_a=0,L_a,t_a>0,$$

(4)

model a biologically isolated tissue segment in which neither cells nor nutrients flux across the boundaries. The system is closed by nonnegative initial profiles $T_a(x_a,0)=T_{a,0}\left(x_a\right)$, $I_a(x_a,0)=I_{a,0}\left(x_a\right)$, $N_a(x_a,0)=N_{a,0}\left(x_a\right)$.

2.2. Nondimensionalization

Following the classical nondimensionalization strategy used in tumor reaction-diffusion models [26,38], we introduce the scalings

$$x={\displaystyle \frac{x_a}{L_a}},t=p_a{t}_a,T={\displaystyle \frac{T_a}{K_a}},I={\displaystyle \frac{I_a}{I_0}},N={\displaystyle \frac{N_a}{N_0}},$$

(5)

where $I_0$ and $N_0$ are reference immune and nutrient density scales. The resulting dimensionless parameters are

$$D_T={\displaystyle \frac{D_{T,a}}{p_a{L}_a^2}},D_I={\displaystyle \frac{D_{I,a}}{p_a{L}_a^2}},D_N={\displaystyle \frac{D_{N,a}}{p_a{L}_a^2}},$$

(6)

$$\chi ={\displaystyle \frac{{\gamma }_a{I}_0}{p_a}},\alpha ={\displaystyle \frac{{\alpha }_a{N}_0}{p_a}},\eta ={\displaystyle \frac{{\eta }_a{K}_a}{p_a{I}_0}},\mu ={\displaystyle \frac{{\mu }_a}{p_a}},\theta ={\displaystyle \frac{{\delta }_a{K}_a}{p_a{N}_0}}.$$

(7)

The dimensionless parameters and their biological interpretations are summarized in Table 2.

2.3. The Governing System

The nondimensional Tumor–Immune–Nutrient system on $x\in [0,1]$, $t>0$ is

$$\begin{array}{cc}\hfill T_t& =D_T{T}_{xx}+T(1-T)-\chi TI+\alpha TN,\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill I_t& =D_I{I}_{xx}+\eta T-\mu I,\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill N_t& =D_N{N}_{xx}-\theta TN,\hfill \end{array}$$

(10)

with no-flux boundary conditions $T_x=I_x=N_x=0$ at $x=0,1$ and nonnegative initial data $T(x,0)=T_0\left(x\right)$, $I(x,0)=I_0\left(x\right)$, $N(x,0)=N_0\left(x\right)$.

The biological mechanisms encoded in (8)–(10) are as follows. Tumor proliferation follows a logistic model modulated by immune-mediated killing at rate $\chi$ and nutrient-driven enhancement at rate $\alpha$. Immune dynamics are governed by tumor-stimulated recruitment at rate $\eta$ and apoptotic decay at rate $\mu$. The nutrient field diffuses rapidly relative to cellular motion, reflecting the physical reality that $D_N\gg D_T$ for small metabolic substrates, and is consumed by tumor metabolic activity at rate $\theta$. This creates a starvation feedback loop in which high tumor density locally depletes nutrients, subsequently limiting further proliferation and creating spatial heterogeneity.

3. Analytical Properties of the System

3.1. Positivity and Invariant Region

A biologically meaningful model must guarantee that all state variables remain nonnegative and uniformly bounded for all finite time. Both properties follow from standard parabolic comparison arguments.

Theorem 1

(Positivity). Let $T_0,I_0,N_0\ge 0$ be continuous on $[0,1]$. Under homogeneous Neumann boundary conditions, the solution satisfies $T(x,t)\ge 0$, $I(x,t)\ge 0$, $N(x,t)\ge 0$ for all $x\in [0,1]$, $t>0$.

Proof. 

At $T=0$ the reaction term in (8) vanishes, and diffusion with no-flux boundaries cannot create negativity. At $I=0$, equation (9) gives $I_t=\eta T\ge 0$. At $N=0$, equation (10) gives $N_t=D_N{N}_{xx}\ge 0$ at interior spatial minima. Hence ${\mathbb{R}}_{+}^3$ is positively invariant. □

Theorem 2

(Uniform Boundedness). There exist constants $T_{max},I_{max},N_{max}>0$ depending only on initial data and parameters such that $0\le T\le T_{max}$, $0\le I\le I_{max}$, $0\le N\le N_{max}$ for all $x\in [0,1]$, $t>0$.

Proof. 

Since $N_t-D_N{N}_{xx}=-\theta TN\le 0$, the maximum principle gives $N\le N_{max}:={max}_x{N}_0\left(x\right)$. Dropping the nonpositive term $-\chi TI$ in (8) and using $N\le N_{max}$ yields a logistic upper bound $T\le T_{max}:=1+\alpha N_{max}$. Finally, with $T\le T_{max}$, equation (9) gives $I\le I_{max}:=\eta T_{max}/\mu$. The compact rectangular region $\mathcal{E}=[0,T_{max}]\times [0,I_{max}]\times [0,N_{max}]$ is therefore positively invariant. □

3.2. Well-Posedness of the Reaction–Diffusion System

For completeness, we summarize the analytical properties ensuring that the tumor–immune–nutrient system is mathematically well-defined. Full proofs follow standard results in analytic semigroup theory and are omitted for brevity.

Linear Operator.

Let $U={(T,I,N)}^{\top }$, $X={\left(L^2(0,1)\right)}^3$, and define the diffusion operator

$$A=\mathrm{diag}(D_T\Delta ,D_I\Delta ,D_N\Delta ),$$

with homogeneous Neumann boundary conditions. Each scalar Laplacian with Neumann boundary conditions is self-adjoint, dissipative, and has compact resolvent. Hence A generates an analytic, contractive $C_0$-semigroup on X [17,32].

Nonlinear Reaction Term.

The nonlinear mapping

$$F\left(U\right)=\left(\begin{array}{c}T(1-T)-\chi TI+\alpha TN\\ \eta T-\mu I\\ -\theta TN\end{array}\right)$$

is locally Lipschitz on X. On the positively invariant region established in Section 3.1, F is globally Lipschitz since each component is polynomial of degree at most two and all partial derivatives are uniformly bounded on $\mathcal{E}$.

Global Well-Posedness.

The semilinear evolution equation

$$U_t=AU+F\left(U\right),U\left(0\right)=U_0\ge 0,$$

follows the abstract framework developed in Egwurube [16] and extended in Orakwelu et al. [29], where existence, asymptotic stability, and contraction of solutions of evolution equations of the form $du/dt=Au+f(u,t)$ on Banach spaces were established using strongly continuous semigroups and resolvent estimates. Global well-posedness follows from Pazy [32] and Engel and Nagel [17]. The tumor–immune–nutrient PDE system is globally well-posed, positive, and uniformly bounded. These properties guarantee that the model is mathematically consistent and suitable for high-order implicit time integration methods such as the SSHBBDF schemes developed in this work.

3.3. Equilibrium Analysis

Setting all time derivatives to zero and neglecting diffusion yields two biologically distinct steady states from the factorization $-\theta T^{*}{N}^{*}=0$.

The tumor-free equilibrium $E_0=(0,0,N_0)$ arises when $T^{*}=0$, giving $I^{*}=0$ and leaving nutrients at homeostatic level $N_0$. This represents healthy tissue in which the tumor has been eliminated or failed to establish.

The coexistence equilibrium $E_1$ arises when $T^{*}\ne 0$, forcing $N^{*}=0$ due to complete metabolic exhaustion, and yields

$$T^{*}={\displaystyle \frac{\mu }{\mu +\chi \eta }},I^{*}={\displaystyle \frac{\eta }{\mu }}{T}^{*},N^{*}=0,$$

(11)

representing a chronic disease state in which the tumor persists at a size limited by immune pressure despite complete nutrient depletion.

3.4. Local Stability

The Jacobian of the reaction system at $(T^{*},I^{*},N^{*})$ is

$$J=\left(\begin{array}{ccc}1-2T^{*}-\chi I^{*}+\alpha N^{*}& -\chi T^{*}& \alpha T^{*}\\ \eta & -\mu & 0\\ -\theta N^{*}& 0& -\theta T^{*}\end{array}\right).$$

(12)

Theorem 3

(Instability of $E_0$).The tumor-free equilibrium $E_0=(0,0,N_0)$ is unstable for any $N_0>0$.

Proof. 

At $E_0$ the Jacobian is upper triangular with diagonal entries ${\lambda }_1=1+\alpha N_0>0$, ${\lambda }_2=-\mu <0$, ${\lambda }_3=0$. The positive eigenvalue ${\lambda }_1$ implies that any small tumor inoculation will grow initially rather than be cleared. □

Theorem 4

(Stability of $E_1$).The coexistence equilibrium $E_1$ is locally asymptotically stable.

Proof. 

At $E_1$ with $N^{*}=0$, the eigenvalue ${\lambda }_3=-\theta T^{*}<0$ decouples the nutrient direction. Using the equilibrium relation $1-T^{*}-\chi I^{*}=0$, the remaining $2\times 2$ tumor–immune block has trace $-T^{*}-\mu <0$ and determinant $\mu T^{*}+\chi \eta T^{*}>0$. By the Routh–Hurwitz criterion, all eigenvalues of $J\left(E_1\right)$ have negative real parts, confirming local asymptotic stability. □

4. The SSHBBDF Framework

Block methods for differential equations were introduced by Milne [27] and developed by many researchers [11,18,37,40]. To overcome the Dahlquist stability barrier, hybrid block methods incorporating at least one non-integer function evaluation were proposed by Butcher [7], Gear [20], and Gragg and Stetter [21]. The Backward Differentiation Formula, the most widely used class of multistep methods for stiff ODEs, was introduced by Curtiss and Hirschfelder [12]. The multistep collocation approach underlying the SSHBBDF derivation follows the continuous interpolant framework established by Adee et al. [1] and Onumanyi et al. [31], who developed hybrid block formulas of order four using non-integer off-grid points to improve starting accuracy, demonstrating that the placement of intra-step evaluation points has a decisive influence on both order and global error. Within the past decade, researchers have clearly demonstrated the dependence of the accuracy of block hybrid methods on the number and placement of grid points in the derivation process [2,28,30,33,35].

We present a class of self-starting SSHBBDF of the form

$${\displaystyle \sum _{j=0}^k}{\alpha }_j{y}_{n+j}+{\displaystyle \sum _{\rho =1}^m}{\alpha }_{{\nu }_{\rho }}{y}_{n+{\nu }_{\rho }}=h{\beta }_k{f}_{n+k},$$

(13)

where $k\in {\mathbb{Z}}^{+}$, $|{\alpha }_k|=1$, and ${\nu }_{\rho }\in (0,k)$ are $m=2n+1$, $n=1,2,3,\dots$, equally spaced non-integer off-grid points defined by

$$x_{n+{\nu }_{\rho }}=x_n+{\displaystyle \frac{\rho }{m+1}}h,\rho =1,2,\dots ,m.$$

(14)

Definition 1.

A Single Step Hybrid Block BDF is the coupled system

$$\left\{\begin{array}{c}{\displaystyle {\displaystyle \sum _{j=0}^k}{\alpha }_j{y}_{n+j}+{\displaystyle \sum _{\rho =1}^m}{\alpha }_{{\nu }_{\rho }}{y}_{n+{\nu }_{\rho }}=h{\beta }_k{f}_{n+k},}\hfill \\ {\displaystyle {\displaystyle \sum _{j=0}^{k-1}}{\overline{\alpha }}_j{y}_{n+j}+{\displaystyle \sum _{\rho =1}^m}{\overline{\alpha }}_{{\nu }_{\rho }}{y}_{n+{\nu }_{\rho }}=h{\overline{\beta }}_k{f}_{n+k}+h{\overline{\beta }}_{{\nu }_{\widehat{\rho }}}{f}_{n+{\nu }_{\widehat{\rho }}},\widehat{\rho }=1,\dots ,m,}\hfill \end{array}\right.$$

(15)

where all coefficients are determined by the collocation conditions of Section 4.2.

Conjecture 5

(Generalized Intra-Step Optimization Principle). Let $k\ge 1$ and $m=2n+1$, $n\in \mathbb{N}$. Define ${\mathcal{V}}_m=\{\nu :0<{\nu }_1<\dots <{\nu }_m<k\}$. Then there exists a nonempty subset ${\mathcal{V}}_m^{*}\subseteq {\mathcal{V}}_m$ such that each $\nu \in {\mathcal{V}}_m^{*}$ yields a hybrid scheme whose leading local truncation error is strictly smaller than that of any other configuration in ${\mathcal{V}}_m$, while its stability region satisfies $S_m\left(\nu \right)\supseteq S_m^{\mathrm{BDF}}$.

4.1. Derivation of the Continuous Approximation

We seek a polynomial $Y\left(t\right)={\sum }_{j=0}^{p+q-1}{b}_j{x}^j$, where $p=m+1$ interpolation points and $q=k$ collocation points determine the degree $p+q-1$. Coefficients $b_j$ are obtained by imposing interpolation conditions at the $m+1$ grid and off-grid values,

$${\displaystyle \sum _{j=0}^{m+1}}{b}_j{x}_{n+{\nu }_{\rho }}^j=y_{n+{\nu }_{\rho }},\rho =0,1,\dots ,m,$$

(16)

and a collocation condition at $x_{n+k}$,

$${\displaystyle \sum _{j=0}^{m+1}}j{b}_j{x}_{n+k}^{j-1}=f_{n+k}.$$

(17)

This system $AB=C$ is solved by matrix inversion. Evaluating $Y\left(x\right)$ at $x_{n+k}$ yields the main formula (13), and evaluating $Y^{\prime }\left(x\right)$ at the off-grid points yields the supplementary formulas.

4.2. Specific Schemes

SSHBBDF with $m=3$.

Setting $k=1$, $m=3$, and off-grid points $\{1/4,1/2,3/4\}$, the derivation yields the main formula

$$-y_{n+1}-\frac{36}{25}{y}_{n+{\displaystyle \frac{1}{2}}}+\frac{16}{25}{y}_{n+{\displaystyle \frac{1}{4}}}+\frac{48}{25}{y}_{n+{\displaystyle \frac{3}{4}}}-\frac{3}{25}{y}_n=-\frac{3}{25}h{f}_{n+1},$$

(18)

together with supplementary formulas at each off-grid point:

$$\begin{array}{cc}\hfill \frac{26}{25}{y}_n-\frac{138}{25}{y}_{n+{\displaystyle \frac{1}{2}}}+\frac{78}{25}{y}_{n+{\displaystyle \frac{1}{4}}}+\frac{34}{25}{y}_{n+{\displaystyle \frac{3}{4}}}& =\frac{1}{25}h{f}_{n+1}-h{f}_{n+{\displaystyle \frac{1}{4}}},\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill -\frac{12}{25}{y}_{n+{\displaystyle \frac{1}{2}}}+\frac{72}{25}{y}_{n+{\displaystyle \frac{1}{4}}}-\frac{152}{75}{y}_{n+{\displaystyle \frac{3}{4}}}-\frac{28}{75}{y}_n& =-\frac{1}{25}h{f}_{n+1}-h{f}_{n+{\displaystyle \frac{1}{2}}},\hfill \end{array}$$

(20)

$$\begin{array}{cc}\hfill \frac{34}{75}{y}_n+\frac{186}{25}{y}_{n+{\displaystyle \frac{1}{2}}}-\frac{66}{25}{y}_{n+{\displaystyle \frac{1}{4}}}-\frac{394}{75}{y}_{n+{\displaystyle \frac{3}{4}}}& =\frac{3}{25}h{f}_{n+1}-h{f}_{n+{\displaystyle \frac{3}{4}}}.\hfill \end{array}$$

(21)

SSHBBDF with $m=5$.

Setting $k=1$, $m=5$, and off-grid points $\{1/6,1/3,1/2,2/3,5/6\}$, the analogous collocation derivation yields a $6\times 6$ linear system whose solution gives the coefficient matrices $A_5^{\left(1\right)}$, $A_5^{\left(0\right)}$, and $B_5^{\left(1\right)}$ recorded in Section 5.

5. Analysis of the SSHBBDF Class

5.1. Block Representation

The SSHBBDF update for each block is written as

$$A^{\left(1\right)}{Y}_{\omega }=A^{\left(0\right)}{Y}_{\omega -1}+h{B}^{\left(1\right)}{F}_{\omega }+h{B}^{\left(0\right)}{F}_{\omega -1},\omega =1,2,\dots ,$$

(22)

where $Y_{\omega }$ and $F_{\omega }$ are block vectors of solution and function values, and $A^{\left(1\right)},A^{\left(0\right)},B^{\left(1\right)},B^{\left(0\right)}$ are $(m+1)\times (m+1)$ matrices. For $m=3$ these are

$$A^{\left(1\right)}=\left(\begin{array}{cccc}-\frac{16}{25}& \frac{36}{25}& -\frac{48}{25}& 1\\ -\frac{78}{25}& \frac{138}{25}& -\frac{34}{25}& 0\\ -\frac{72}{25}& \frac{12}{25}& \frac{152}{75}& 0\\ \frac{66}{25}& -\frac{186}{25}& \frac{394}{75}& 0\end{array}\right),A^{\left(0\right)}=\left(\begin{array}{cccc}0& 0& 0& \frac{3}{25}\\ 0& 0& 0& -\frac{26}{25}\\ 0& 0& 0& \frac{28}{75}\\ 0& 0& 0& -\frac{34}{75}\end{array}\right),$$

(23)

$$B^{\left(1\right)}=\left(\begin{array}{cccc}0& 0& 0& \frac{3}{25}\\ 1& 0& 0& -\frac{1}{25}\\ 0& 1& 0& \frac{1}{25}\\ 0& 0& 1& -\frac{3}{25}\end{array}\right).$$

(24)

For $m=5$, the $6\times 6$ matrices $A_5^{\left(1\right)}$, $A_5^{\left(0\right)}$, and $B_5^{\left(1\right)}$ follow directly from the collocation derivation; their entries are given explicitly in the supplementary material.

5.2. Local Truncation Error

The local truncation error operator associated with the SSHBBDF is

$$L[y\left(x\right),h]={\displaystyle \sum _{j=0}^k}{\alpha }_{j}y(x+jh)+{\displaystyle \sum _{\rho =1}^m}{\alpha }_{{\nu }_{\rho }}y(x+{\nu }_{\rho }h)-h{\beta }_k{y}^{\prime }(x+kh)-h{\displaystyle \sum _{\widehat{\rho }=1}^m}{\beta }_{{\nu }_{\widehat{\rho }}}{y}^{\prime }(x+{\nu }_{\widehat{\rho }}h).$$

(25)

Expanding in Taylor series about x gives $L[y,h]={\sum }_{q=0}^{\infty }{c}_q{h}^q{y}^{\left(q\right)}\left(x\right)$. The method has order p if $c_0=\dots =c_p=0$ and $c_{p+1}\ne 0$ [24]. Our calculations demonstrate that the SSHBBDF class achieves high order with reassuringly small error constants. For the SSHBBDF with $m=3$ and $m=5$ off-grid points, symbolic computation yields the following principal local truncation error terms for each formula in the block:

$$m\left(3\right)=\left\{\begin{array}{c}{\displaystyle \frac{3h^5{y}^{\left(5\right)}\left(x\right)}{32000}}+O\left(h^6\right)\hfill \\ \\ -{\displaystyle \frac{29h^5{y}^{\left(5\right)}\left(x\right)}{128000}}+O\left(h^6\right)\hfill \\ \\ {\displaystyle \frac{31h^5{y}^{\left(5\right)}\left(x\right)}{192000}}+O\left(h^6\right)\hfill \\ \\ -{\displaystyle \frac{37h^5{y}^{\left(5\right)}\left(x\right)}{128000}}+O\left(h^6\right),\hfill \end{array}\right.$$

$$m\left(5\right)=\left\{\begin{array}{c}{\displaystyle \frac{5h^7{y}^{\left(7\right)}\left(x\right)}{24004512}}+O\left(h^8\right)\hfill \\ \\ -{\displaystyle \frac{53h^7{y}^{\left(7\right)}\left(x\right)}{96018048}}+O\left(h^8\right)\hfill \\ \\ {\displaystyle \frac{h^7{y}^{\left(7\right)}\left(x\right)}{4445280}}+O\left(h^8\right)\hfill \\ \\ -{\displaystyle \frac{167h^7{y}^{\left(7\right)}\left(x\right)}{960180480}}+O\left(h^8\right)\hfill \\ \\ {\displaystyle \frac{59h^7{y}^{\left(7\right)}\left(x\right)}{240045120}}+O\left(h^8\right)\hfill \\ \\ -{\displaystyle \frac{23h^7{y}^{\left(7\right)}\left(x\right)}{32006016}}+O\left(h^8\right),\hfill \end{array}\right.$$

$$m\left(3\right):\mathcal{O}\left(h^5\right),m\left(5\right):\mathcal{O}\left(h^7\right).$$

(26)

Since the leading term of the local truncation error is $c_{p+1}{h}^{p+1}{y}^{(p+1)}$ for a method of order p [24], the relations in (26) confirm that the SSHBBDF with $m=3$ is of 4 and the SSHBBDF with $m=5$ is of 6.

5.3. Zero Stability

As $h\to 0$, the block system (22) reduces to $A^{\left(1\right)}{Y}_{\omega }=A^{\left(0\right)}{Y}_{\omega -1}$. The first characteristic polynomial of the class is

$$\chi \left(R\right)=det\left[R{A}^{\left(1\right)}-A^{\left(0\right)}\right]=\xi R^m(1+R),\xi \in \mathbb{Q}.$$

(27)

All roots satisfy $|R_j|\le 1$ and those with $|R_j|=1$ have multiplicity one, so the SSHBBDF class is zero stable [30].

5.4. Linear Stability

Applying the SSHBBDF to the test equation $y^{\prime }=\lambda y$, $\lambda <0$, yields $Y_{\omega }=D\left(z\right)Y_{\omega -1}$ where $D\left(z\right)={(A^{\left(1\right)}-z{B}^{\left(1\right)})}^{-1}{A}^{\left(0\right)}$ and $z=\lambda h$. The stability function $R\left(z\right)$ is the dominant eigenvalue of $D\left(z\right)$.

For $m=3$ the stability function is

$$R\left(z\right)=-{\displaystyle \frac{3z^3+44z^2+288z+768}{3z^4-25z^3+140z^2-480z+768}}.$$

(28)

For $m=5$ the stability function is

$$R\left(z\right)=-{\displaystyle \frac{2\left(5z^5+137z^4+2025z^3+18360z^2+97200z+233280\right)}{10z^6-147z^5+1624z^4-13230z^3+75600z^2-272160z+466560}}.$$

(29)

The stability regions defined by $\left|R\right(z\left)\right|\le 1$ are plotted in Figure 1. The SSHBBDF with $m=3$ is A-stable, as the stability region contains the entire left half-plane. The scheme with $m=5$ is $L_0$-stable [9], satisfying ${max}_{z\le 0}\left|R\left(z\right)\right|\le 1$ and ${lim}_{z\to -\infty }R\left(z\right)=0$.

6. Numerical Implementation

6.1. Method of Lines Spatial Discretization

We discretize the spatial domain $[0,1]$ using a uniform grid of $M+1$ nodes $x_i=i\Delta x$, $\Delta x=1/M$, $i=0,\dots ,M$. The Laplacian is approximated by standard second-order central differences,

$${\displaystyle \frac{{\partial }^{2}U}{\partial x^2}}{|}_{x_i}\approx {\displaystyle \frac{U_{i+1}\left(t\right)-2U_i\left(t\right)+U_{i-1}\left(t\right)}{{(\Delta x)}^2}},$$

(30)

with ghost points enforcing the Neumann conditions: $U_{-1}=U_1$ and $U_{M+1}=U_{M-1}$, giving modified boundary stencils

$$\Delta U{|}_{x_0}\approx {\displaystyle \frac{2U_1-2U_0}{{(\Delta x)}^2}},\Delta U{|}_{x_M}\approx {\displaystyle \frac{2U_{M-1}-2U_M}{{(\Delta x)}^2}}.$$

(31)

Substituting into (8)–(10) yields the semi-discrete system of $3(M+1)$ coupled ODEs, written compactly as

$${\displaystyle \frac{dY}{dt}}=\mathbf{A}Y+F\left(Y\right),Y\left(0\right)=Y_0,$$

(32)

where $Y\left(t\right)={[T_0,\dots ,T_M,I_0,\dots ,I_M,N_0,\dots ,N_M]}^{\top }\in {\mathbb{R}}^{3(M+1)}$, $\mathbf{A}$ is the sparse block-diagonal diffusion matrix, and $F\left(Y\right)$ contains the nonlinear reaction terms.

6.2. Time Integration

The SSHBBDF schemes are applied to integrate (32). Modified Newton iteration is employed at each time step to handle the nonlinear reaction terms, with the diffusion component treated by Gaussian elimination with partial pivoting. The block structure of the SSHBBDF allows simultaneous evaluation at all off-grid points within each step, enabling natural parallelization and reducing the number of function evaluations relative to sequential approaches.

7. Numerical Results

7.1. Parameter Values

Table 3 and Table 4 present two parameter sets representing biologically distinct scenarios. Tumor proliferation rates and immune decay rates follow de Pillis et al. [13], who validated a tumor-immune model against mouse and human chromium release assay data. Diffusion coefficients for tumor cells are consistent with melanoma cell migration measurements reported by Stepien et al. [39], and the nutrient diffusion coefficient is set substantially higher than cellular diffusivities to reflect the rapid equilibration of small metabolic substrates relative to cell-scale motion [19].

7.2. Convergence Verification

Since this is the first numerical study of this exact three-component Tumor–Immune–Nutrient system, no reference solution from the prior literature is available. Temporal convergence is therefore assessed through a four-level Richardson step-doubling procedure using successively halved time-step sizes $\Delta t\in \{0.01,0.005,0.0025,0.00125\}$. All solutions are evaluated at the shared output point $T_{\mathrm{final}}=1$, the unique grid point common to all four meshes, which eliminates interpolation error entirely. The off-grid stage values at $t+h/4$, $t+h/2$, $t+3h/4$ are internal block quantities and are excluded from the comparison.

The Experimental Order of Convergence (EOC) is estimated via

$$\mathrm{EOC}={log}_2\left({\displaystyle \frac{∥Y^{\Delta t}-Y^{\Delta t/2}∥}{∥Y^{\Delta t/2}-Y^{\Delta t/4}∥}}\right),$$

(33)

where the norm $\parallel ·\parallel$ is evaluated over all $N=60$ spatial nodes at $T_{\mathrm{final}}$. Three norms are computed in parallel to confirm robustness of the estimate:

$${\parallel \mathbf{d}\parallel }_{\infty }=\underset{1\le j\le 60}{max}|d_j{|,\parallel \mathbf{d}\parallel }_2={\left({\displaystyle \sum _{j=1}^{60}}{\left|d_j\right|}^2\right)}^{1/2}{,\parallel \mathbf{d}\parallel }_1={\displaystyle \sum _{j=1}^{60}}\left|d_j\right|,$$

(34)

where $d_j=|Y_j^{\Delta t}-Y_j^{\Delta t/2}|$ is the pointwise difference at node j. The four step sizes yield three successive difference vectors

$${\mathbf{d}}_{12}:\Delta t=0.01\to 0.005,{\mathbf{d}}_{23}:\Delta t=0.005\to 0.0025,{\mathbf{d}}_{34}:\Delta t=0.0025\to 0.00125,$$

and two independent EOC estimates,

$${\mathrm{EOC}}_1={log}_2{\displaystyle \frac{\parallel {\mathbf{d}}_{12}\parallel }{\parallel {\mathbf{d}}_{23}\parallel }},{\mathrm{EOC}}_2={log}_2{\displaystyle \frac{\parallel {\mathbf{d}}_{23}\parallel }{\parallel {\mathbf{d}}_{34}\parallel }}.$$

Table 5 records the global EOC values under all three norms. All six estimates lie in the narrow band $[3.98,4.00]$, in excellent agreement with the theoretical order $p=4$ of the SSHBBDF scheme with $m=3$, as established analytically in Section 7.2. The values of ${\mathrm{EOC}}_2$ are marginally closer to $4.0$ than those of ${\mathrm{EOC}}_1$, consistent with asymptotic convergence: as $\Delta t$ decreases, higher-order error terms become negligible and the EOC converges monotonically toward the true order $p=4$. The minor deviation from $4.0$ at coarser refinement levels is a preasymptotic effect fully expected for finite step sizes.

Figure 2 shows the three pointwise difference vectors $|{\mathbf{d}}_{12}|$, $|{\mathbf{d}}_{23}|$, and $|{\mathbf{d}}_{34}|$ across all 60 spatial nodes. The successive geometric reduction in magnitude as $\Delta t$ is halved is visually pronounced: each curve lies approximately $2^4=16$ times below the previous one, which is the hallmark signature of a 4th-order scheme and provides direct graphical confirmation of the values in Table 5. The y-axis upper bound of $1.62\times {10}^{-12}$ at the finest refinement level reflects the very high absolute accuracy achieved.

Figure 3 overlays the pointwise ${\mathrm{EOC}}_1$ and ${\mathrm{EOC}}_2$ at each spatial node $j=1,\dots ,60$. Both curves cluster tightly around the theoretical value $p=4$ across the entire spatial domain, demonstrating that convergence is spatially uniform and not merely an average effect. The ${\mathrm{EOC}}_2$ curve sits marginally closer to the reference line, confirming the expected asymptotic behaviour. Isolated spikes near $j=18$ and $j=55$–56 occur where the solution value is very close to zero, making the ratio of successive differences sensitive to floating-point rounding; they carry no physical significance and do not indicate a local loss of accuracy.

7.3. Spatiotemporal Dynamics

Figure 4 and Figure 5 present three-dimensional spatiotemporal plots of tumor density, nutrient concentration, and immune density under the two parameter regimes.

Under parameter set 1, tumor density decreases in regions of high immune activity and nutrient availability, indicating effective immunological suppression. The immune population concentrates near the tumor boundary, and nutrient levels remain sufficient to sustain ongoing immune function. The system evolves toward the stable coexistence equilibrium $E_1$, in which tumor size is limited by immune pressure rather than by metabolic exhaustion.

Under parameter set 2, sustained tumor expansion overwhelms the immune response. Rapid nutrient depletion creates a nutrient-starved core in the spatial interior, consistent with the necrotic cores observed in solid avascular tumors [19,38]. The immune density remains low and spatially peripheral, unable to penetrate the metabolically depleted region. This regime demonstrates therapeutic failure and corresponds dynamically to growth toward a tumor-dominated state far from $E_1$.

8. Conclusion

This study investigated the conditions under which cytotoxic immune responses can suppress solid tumor growth versus being overwhelmed by aggressive expansion and metabolic competition. The central biological finding is that the outcome — immune-mediated control or tumor escape — is determined by a combination of three dimensionless parameters: the immune cytotoxic strength $\chi$, the tumor-driven immune recruitment ratio $\eta /\mu$, and the nutrient consumption rate $\theta$. When $\chi$ and $\eta /\mu$ are sufficiently large relative to $\theta$, the system evolves toward the chronic coexistence equilibrium $E_1$, in which immune pressure limits tumor size despite complete nutrient depletion at steady state. When these conditions are not met, rapid nutrient exhaustion creates a metabolically hostile interior that excludes immune infiltration, generating a spatially peripheral immune response and a necrotic core consistent with avascular solid tumor morphology observed clinically.

These results have direct implications for immunotherapy design. The instability of the tumor-free equilibrium $E_0$ — which holds for any positive nutrient level — establishes mathematically that passive immune surveillance alone is insufficient to eliminate an established tumor, consistent with clinical experience. The stability of $E_1$ suggests that partial immune control, rather than eradication, may be the achievable target in many solid tumor contexts, and that interventions which boost $\chi$ (e.g., checkpoint inhibitors) or $\eta$ (e.g., adoptive cell transfer) may shift a patient from the escape regime into the controllable coexistence regime. The spatial simulations further suggest that nutrient competition — and therefore anti-angiogenic or metabolic co-therapies — may critically affect whether immune cells can physically access the tumor interior.

The mathematical framework supporting these conclusions is rigorous: positivity, uniform boundedness, global well-posedness, and local asymptotic stability of $E_1$ are all established analytically. Numerical simulations are performed using a fourth-order A-stable SSHBBDF implicit integrator, whose convergence is verified to fourth order across all spatial nodes and three norms ($\mathrm{EOC}\in [3.98,4.00]$), ensuring that all reported dynamics reflect the true behavior of the biological model rather than numerical artifacts.

Future directions include extension to two- and three-dimensional tumor geometries, incorporation of angiogenesis and drug pharmacokinetics, parameter estimation from clinical imaging data, and application to patient-specific treatment planning for immunotherapy protocols. Computational frameworks of the kind developed here represent a step toward evidence-based tools for predicting individual tumor–immune trajectories and identifying the combination of biological conditions most amenable to immune-mediated tumor control.

As computational models of tumor-immune dynamics grow in clinical relevance, the development of reliable, high-order numerical frameworks of the kind presented here represents a step toward more robust evidence-based tools for cancer treatment policy.