PaperlayerView exact revision on arXiv

Source-linked AI summary

Computational Oncology of Chemotaxis-Driven Tumour--Immune Spatial Patterning and Stability

Zonghao Liu, Jiguang Yu, Lei Su, Louis Shuo Wang, Yang Du, Jingfeng Liu

arXiv:2607.03813v2math.APcs.NAmath.NA

TL;DR

Spatial tumour–immune dynamics require models that connect local interactions with tissue-scale patterns while respecting biological constraints. This paper analyzes and simulates a reaction–diffusion–chemotaxis model, finding that chemotaxis can destabilize homogeneous coexistence without changing the tumour-invasion threshold.

  • Problem

    Tumour progression is spatially structured, and computational models must preserve nonnegative densities and no-flux boundaries while producing interpretable stability criteria.

    Source: §1.1 · pp. 1–2

  • Method

    The paper analyzes a nondimensional tumour–immune–chemokine reaction–diffusion–chemotaxis system and tests its predictions with a conservative finite-volume scheme and numerical diagnostics.

    Source: §8 · p. 37

  • Results

    Chemotaxis leaves the tumour-invasion threshold unchanged but can trigger finite-wavelength oscillatory instability above a critical sensitivity; simulations reproduce analytical thresholds and equilibria.

    Source: §8 · p. 37, §7.8 · p. 37

  • Takeaways & Limitations

    The framework connects immune-control thresholds and chemotaxis-driven spatial patterning, with simulations mapping the stability boundary across three parameter planes.

    Source: §8 · p. 37

    The a priori estimates do not establish full global L∞-boundedness of the three-component chemotaxis system in arbitrary spatial dimensions.

    Source: §4.2 · p. 17

Abstract

from arXiv · show

We develop a reaction--diffusion--chemotaxis model for spatial tumour--immune--chemokine dynamics that couples logistic tumour growth, immune-mediated killing, chemokine-dependent immune recruitment, chemotactic migration, and signal production. For the nondimensional system, we establish local classical solvability, nonnegativity, a uniform tumour-density bound, and global mass estimates for the immune and chemokine components. The tumour-free equilibrium is stable precisely when the baseline immune-control index satisfies \(\sigma_0/\delta>1\), whereas positive homogeneous coexistence is characterized by a scalar nonlinear equation. Linearization in the Neumann Laplacian eigenbasis yields a mode-dependent cubic dispersion relation, showing that chemotaxis does not alter the tumour-invasion threshold but can destabilize homogeneous coexistence through a finite-wavelength oscillatory instability above a critical sensitivity \(\xi_c\). A conservative finite-volume discretization with upwind chemotactic fluxes and implicit backward differentiation formula time integration is used to test these predictions. Numerical experiments recover the analytical equilibria and growth rates, identify the dominant unstable mode, reproduce the transition to spatial heterogeneity, and quantify the effects of immune recruitment, decay, and diffusion on the stability boundary. Grid-refinement, mass-balance, residual, and nonnegativity diagnostics support the computational reliability of the results.

1. Introduction

The introduction motivates spatial tumour modelling as a way to connect local cell–signal interactions with tissue-scale patterns, then presents a tumour–immune–chemokine feedback model. It situates this framework among mathematical and computational oncology approaches and outlines the paper’s analysis and simulations.

1.1. Motivation: computational tumour modelling

  • Spatial tumour models connect local cell–signal interactions to tissue-scale patterns, while nonnegativity, no-flux boundaries, and interpretable stability criteria help distinguish biological behaviour from numerical artifacts.The paper frames these requirements as biologically meaningful, not merely technical.Source: §1.1 · pp. 1–2

1.2. Biological setting: tumour–immune–chemokine interactions

  • Tumour cells grow logistically and produce chemokines; chemokines recruit and guide immune effectors, which locally kill tumour cells.Tumour–immune contacts provide an additional chemokine source.Source: §1.2 · pp. 2, 4, Figure 1 · p. 3

1.3. Related work and modelling context

  • The model adds a continuum account of chemokine-directed immune chemotaxis to tumour–immune models and spatial approaches, while using a conservative finite-volume scheme with positivity and residual diagnostics.The stated contributions also include analytical thresholds for tumour-free control, coexistence, and finite-wavelength instability.Source: §1.3 · pp. 4–5

1.4. Outline

  • The paper formulates and nondimensionalizes the model, develops its numerical scheme and mathematical analysis, identifies equilibria and stability, then reports experiments verifying theoretical results..Source: §1.4 · p. 5

2. Methods (Model Formulation)

The methods formulate tumour, immune, and chemokine dynamics as a reaction–diffusion–chemotaxis system on a tissue domain. The formulation specifies no-flux boundaries and nonnegative initial data, then nondimensionalizes the system for analysis and simulation.

2.1. Dimensional tumour–immune–chemokine model

  • The dimensional model represents tumour, immune-effector, and chemokine fields on a bounded smooth tissue domain in one, two, or three dimensions.Tumour growth is logistic, immune cells kill tumour cells and move up chemokine gradients, and tumour cells and tumour–immune interactions produce chemokine.Source: §2.1 · pp. 5–6

2.2. Boundary and initial conditions

  • No-flux boundaries describe a closed tissue region, and initial tumour, immune, and chemokine fields are nonnegative and sufficiently smooth when classical solutions are considered.The stability analysis imposes componentwise Neumann conditions.Source: §2.2 · p. 8

2.3. Nondimensionalization

  • The dimensional system is nondimensionalized using a reference tissue length, yielding tumour, immune, and chemokine variables scaled by carrying capacity, tumour-killing relative to growth, and recruitment half-saturation.The resulting system makes explicit how recruitment and chemotaxis may stabilize homogeneous coexistence or destabilize it into spatial patterns.Source: §2.3 · pp. 8–9
  • Chemokine diffusion, production by tumour and tumour–immune interactions, and degradation are represented in the third nondimensional equation.The nondimensional system also includes diffusion coefficients, chemotactic sensitivity, immune recruitment and loss, and chemokine production and degradation parameters.Source: §2.3 · pp. 8–9

2.2. Boundary and initial conditions

  • Table 1 gives representative dimensional parameter orders for solid-tumour tissue, using values to set scales rather than to model a specific tumour.The no-flux formulation supports use of the Neumann Laplacian eigenbasis in stability and dispersion calculations.Source: §2.2 · p. 8, Table 1 · p. 9

2.3. Nondimensionalization

  • Table 2 lists the nondimensional parameters used in system (2.4).The nondimensional no-flux conditions are retained after rescaling.Source: §2.3 · p. 9, Table 2 · p. 10
  • The nondimensional base set encodes a three-order motility separation, with d3 = 1, d2 = 0.1, and d1 = 0.01; fast chemokine turnover is deliberately not imposed to retain analytical transparency.The retained turnover rescales but does not remove the chemotaxis-driven finite-wavelength instability.Source: §2.3 · p. 10

3. Methods (Numerical Scheme)

The scheme uses conservative finite-volume fluxes with no-flux boundaries, alongside mass-balance, positivity, and residual diagnostics. Upwind chemotactic fluxes and implicit BDF stepping address directed immune migration and stiff dynamics; the analysis distinguishes explicit-step positivity guarantees from monitored BDF runs.

3.1. Finite-volume discretization and no-flux treatment

  • Centered diffusion fluxes and upwind chemotactic immune fluxes form a conservative finite-volume scheme with zero boundary fluxes.The upwind choice is first-order for chemotaxis but improves robustness and positivity where immune gradients are steep.Source: §3.1 · p. 11

3.2. Mass-balance diagnostics

  • Discrete mass-balance identities follow from interior flux cancellation and provide checks on temporal accuracy and boundary-flux implementation.The mass discrepancies compare the two sides of the discrete balance laws.Source: §3.2 · p. 13

3.3. BDF time stepping with positivity diagnostics

  • Implicit BDF time stepping addresses stiffness, while Newton failures or positivity violations trigger time-step reduction.The method does not itself guarantee positivity; reported runs monitor component minima, which remain nonnegative up to solver tolerance.Source: §3.3 · p. 13

3.4. Residual diagnostics

  • Post-hoc residuals, alongside positivity and mass-balance checks, assess whether computed patterns are consistent with the PDE system rather than discretization artifacts.The residuals use discrete time-derivative, Laplacian, and divergence operators.Source: §3.4 · p. 14

3.5. Consistency, positivity, and convergence status of the scheme

  • The explicit forward-Euler scheme preserves nonnegativity under a computable step-size restriction, but these properties do not prove convergence for the nonlinear cross-diffusion system.Observed grid-refinement order 1.92 is numerical self-convergence evidence in a smooth diffusion-dominated regime; formal spatial order is generally first order.Source: §3.5 · pp. 14–15

3.6. Extension to two-dimensional rectangular domains

  • In two dimensions, applying the finite-volume construction by coordinate direction preserves the conservative structure and uses the same positivity, mass-balance, and residual diagnostics.Diffusive face fluxes use centered differences and chemotactic immune fluxes use upwinding.Source: §3.6 · pp. 15–16

4. Results (Mathematical Well-posedness)

The analysis establishes local classical solvability and nonnegativity for the nondimensional system, then derives a uniform tumour-density bound and global-in-time mass control for immune and chemokine components. These estimates do not establish full global L∞ boundedness in arbitrary spatial dimension.

4.1. Classical solvability and positivity

  • For compatible nonnegative initial data on a bounded smooth domain, the system has a unique nonnegative classical solution up to a maximal existence time.The result relies on normal parabolicity and quasi-positivity of the reaction terms under no-flux boundary conditions.Source: §4.1 · p. 16

4.2. A priori estimates

  • The tumour density remains between zero and max{1, ‖u0‖L∞(Ω)} throughout the classical solution's interval of existence.Source: §4.2 · p. 16
  • Immune and chemokine components have global-in-time mass control, but these estimates do not imply full global L∞ boundedness in arbitrary dimension.Further boundedness requires additional assumptions or estimates; later results are stated for classical solutions on their existence interval or globally defined bounded solutions.Source: §4.2 · p. 17

5. Results (Equilibria)

The tumour-free equilibrium is linearly stable above the baseline immune-control threshold, and chemotaxis does not change that invasion threshold. Positive homogeneous coexistence is reduced to a scalar equation, with an interior equilibrium guaranteed under stated parameter conditions.

5.1. Tumour-free equilibrium

  • The tumour-free equilibrium is linearly stable when σ0>δ, unstable when σ0<δ, and non-hyperbolic at equality; chemotaxis does not alter this threshold.The dominant stability sign is determined by 1−σ0/δ.Source: §5.1 · p. 18

5.2. Coexistence equilibria

  • Every positive homogeneous coexistence equilibrium has 0<u*<1 and is determined by a scalar equation F(u*)=0.The immune and chemokine equilibrium expressions reduce the coexistence condition to this single equation.Source: §5.2 · p. 19
  • If σ0<δ and F(1)>0, at least one interior coexistence equilibrium exists; F(1)>0 holds when σ0>0 or when σ1>0 and α>0.Source: §5.2 · p. 19
  • When σ0=σ1=0, or σ0=0 with σ1>0 and α=0, the endpoint u=1 is an immune-free tumour-dominant state, not positive coexistence.Source: §5.2 · p. 19

6. Results (Linear Stability Analysis)

Linearization in the Neumann eigenbasis shows how chemotaxis can destabilize finite spatial modes while leaving the homogeneous mode unaffected. Under homogeneous-mode stability, a finite critical sensitivity is attained, and its first crossing may be stationary or oscillatory.

6.1. Mode-wise stability matrix

  • Linearized no-flux conditions reduce to componentwise Neumann conditions, permitting mode-by-mode analysis with a cubic characteristic polynomial.The Neumann Laplacian eigenbasis diagonalizes the linearized fluxes and yields the modal stability matrix.Source: §6.1 · p. 20, §6.2 · p. 20
  • Chemotactic sensitivity enters the modal matrix through a term proportional to xi v* mu_k, which vanishes at mu_0=0 and can destabilize finite modes.The chemotactic amplification grows with the spatial eigenvalue, allowing a band of finite wavelengths to destabilize while the uniform mode remains unaffected.Source: §6.1 · p. 20, §6.2 · p. 22

6.2. Dispersion relation and critical sensitivity

  • If coexistence is stable to the homogeneous mode, a finite critical sensitivity xi_c>0 is attained at some finite spatial mode k*>=1.For each nonzero mode, a2 and a3 decrease affinely with sensitivity, and the mode growth rate eventually becomes positive.Source: §6.2 · p. 21
  • The first stability loss is stationary when a3(k*;xi_c)=0, or oscillatory when a1a2-a3=0 with a3>0.At a Hopf crossing, the marginal frequency satisfies omega*^2=a2.Source: §6.2 · p. 21
  • For the computational parameter set, the first unstable mode is k*=5 at xi_c=12.479, with Hopf frequency omega*=0.989 and a3=4.681.The Hopf condition a1a2-a3 is −1.5 × 10^-11 at onset.Source: §6.2 · p. 21
  • At 1.6 xi_c, the fastest-growing mode is k=6, distinct from the first unstable mode k*=5 at onset.The critical onset mode and the dominant mode above onset need not coincide.Source: §6.2 · p. 21

7. Results (Numerical Experiments)

The experiments test the model’s analytical thresholds, equilibria, instability predictions, and numerical scheme. They recover the predicted transition from homogeneous coexistence to chemotaxis-driven patterns and report diagnostics supporting numerical reliability.

7.1. Experiment 1: tumour-free threshold

  • The measured tumour rates match predictions: decay is −0.6677 versus −0.6667 predicted, while invasion is 0.6635 versus 0.6667 predicted, then saturates at u*=0.531.The tumour-free threshold is σ0=δ; chemotaxis does not alter the linear tumour growth rate.Source: §7.1 · pp. 22–23, Figure 4 · p. 24

7.2. Experiment 2: coexistence equilibrium verification

  • The numerical solution converges monotonically to the analytical coexistence equilibrium, with total L2 error 4.4 × 10^-11 and final fields matching the analytical values.Componentwise errors are 1.1 × 10^-11 for u, 4.7 × 10^-12 for v, and 7.7 × 10^-12 for w.Source: §7.2 · p. 25, Figure 5 · p. 26

7.3. Experiment 3: dispersion relation and unstable modes

  • At ξ=1.6ξc, intermediate modes become unstable; the dominant mode k*=6 has predicted growth 0.267 and measured growth 0.273, with active-band correlation 0.980.Homogeneous and weakly chemotactic cases remain stable across modes; strongly damped modes reach the nonlinear noise floor.Source: §7.3 · p. 27, Figure 6 · p. 28

7.4. Experiment 4: chemotaxis-driven pattern formation

  • Across ξc, simulations transition from homogeneous fields to onset modulation and then sustained heterogeneous patterns, with immune aggregation at chemokine maxima and local tumour depletion.At 0.7ξc the amplitude decays to about 10^-9; near ξc it remains about 5 × 10^-3; at 1.8ξc the tumour and immune amplitudes reach about 0.69 and 6.9.Source: §7.4 · p. 29, Figure 7 · p. 30

7.5. Experiment 5: two-parameter sensitivity maps

  • Stronger chemokine-induced recruitment lowers the critical sensitivity, whereas greater immune decay or diffusion raises it and shrinks the unstable region.The stability boundary is mapped in the (ξ, σ1), (ξ, δ), and (ξ, d2) planes.Source: §7.5 · p. 31, Figure 8 · p. 32

7.6. Experiment 6: grid convergence and residual certification

  • Grid refinement gives observed order 1.92; component minima stay positive and PDE residuals remain small and decay in time.The L2 differences decrease from 5.0 × 10^-5 to 3.5 × 10^-6; maximum residuals are about 2 × 10^-3 for Ru, 1 × 10^-2 for Rv, and 6 × 10^-4 for Rw.Source: §7.6 · p. 33, Figure 9 · p. 34

7.7. Experiment 7: two-dimensional pattern and immune-infiltration foci

  • In two dimensions, immune cells form sharp foci at chemokine maxima, where tumour density is depleted, while intervening regions remain immune-poor and tumour-rich.At ξ=1.4ξc, the pattern remains strictly positive in every component.Source: §7.7 · p. 35, Figure 10 · p. 36

7.8. Summary of verification

  • Across six experiments, simulations reproduce analytical thresholds and equilibria within solver tolerance and confirm chemotaxis-driven patterns, with positivity, mass-balance, and residual diagnostics supporting numerical reliability.Table 4 collects quantitative agreement between theory and computation.Source: Table 4 · p. 37, §7.8 · p. 37

8. Conclusion

The paper analyzes a tumour–immune–chemokine reaction–diffusion–chemotaxis model and establishes analytical bounds and stability criteria. Its analysis and numerical experiments identify a finite-wavelength oscillatory instability above a critical chemotactic sensitivity, while full global boundedness remains unresolved.

Conclusion

  • The model couples tumour growth, immune killing, chemokine-driven immune recruitment, and chemotactic migration; analysis establishes local solvability, nonnegativity, a tumour bound, and global immune and chemokine mass estimates.The estimates do not establish global L∞-boundedness for the full three-component system.Source: §8 · p. 37
  • Tumour-free stability is determined by σ0/δ, while linearization about coexistence yields a finite-wavelength instability above critical sensitivity ξc, oscillatory for the parameters studied.Chemotactic sensitivity enters the mode-wise stability matrix through the wavenumber-amplified coupling ξv*μk.Source: §8 · p. 37
  • A conservative finite-volume scheme and its diagnostics support numerical verification of the thresholds, equilibria, dispersion relation, pattern transition, stability maps, and computed patterns.The scheme uses upwind chemotactic fluxes and discrete no-flux boundaries while preserving discrete positivity and mass balance.Source: §8 · p. 37
  • Global boundedness for the full system remains open; proposed directions include additional structural conditions, weakly nonlinear onset analysis, higher-dimensional pattern classification, and treatment terms.The stated possible additional structures include small chemotactic sensitivity, stronger logistic damping, signal-dependent desensitization, or dimension-dependent regularity estimates.Source: §8 · p. 37

A. Proofs of Technical Estimates

The appendix records technical arguments for local solvability, positivity, and a priori estimates, rather than proving global boundedness of the full system.

A. Proofs of Technical Estimates

  • The technical estimates concern the nondimensional system on a bounded smooth domain with no-flux boundary conditions.The initial data are assumed nonnegative.Source: A. Proofs of Technical Estimates · p. 37
  • The appendix’s stated purpose is to support the equilibrium, stability, and numerical verification results with estimates sufficient for those analyses.It explicitly does not seek a full global boundedness theorem for the three-component chemotaxis system.Source: A. Proofs of Technical Estimates · p. 37
  • Nonnegative initial data are an assumption for the technical arguments collected in the appendix.Source: A. Proofs of Technical Estimates · p. 37

A.1. Normal parabolicity, local solvability, and positivity

The appendix establishes normal parabolicity and local solvability, then derives positivity and a priori bounds using comparison and integral estimates.

A.1. Normal parabolicity, local solvability, and positivity

  • The principal diffusion symbol has eigenvalues bounded below by d*|ζ|², giving normal ellipticity with uniform sectoriality constants when |v| is bounded.The no-flux conditions are compatible with the principal part, and standard quasilinear parabolic theory yields a unique maximal classical solution for sufficiently regular compatible data.Source: §A.1 · p. 38
  • The negative-part estimates and continuation argument establish that nonnegative initial data remain nonnegative throughout the classical solution’s interval of existence.Source: §A.1 · pp. 38–39

A.2. Tumour-density bound

  • Comparison with a logistic ordinary differential equation gives 0 ≤ u(t,x) ≤ max{1, ||u0||L∞(Ω)} throughout the solution’s existence interval.The tumour equation is bounded above by the logistic diffusion equation because v ≥ 0.Source: §A.2 · p. 39

A.3. Immune mass estimate; A.4. Chemokine mass estimate

  • Integrating the immune and chemokine equations over Ω and using no-flux conditions yields scalar differential inequalities for their total masses.The chemokine estimate uses the tumour-density bound.Source: §A.3 · pp. 39–40, §A.4 · p. 40

A.5. Why these estimates do not imply full global boundedness

The available mass estimates do not establish full global boundedness for the three-component chemotaxis system. The appendix explains the limitation through chemotactic aggregation and gives the linearization setup used for the stability calculation.

A.5. Why these estimates do not imply full global boundedness

  • The estimates do not establish global L∞-boundedness for the full three-component chemotaxis system in arbitrary dimension.The stated results are local classical well-posedness, positivity, tumour-density boundedness, and mass estimates; stronger a priori control would be needed for a full global boundedness theorem.Source: §A.5 · p. 41
  • Chemotactic aggregation can turn chemokine gradients into immune-density concentration, which L1 bounds on immune cells and chemokine do not control.The passage identifies possible additional structure, including small chemotactic sensitivity, strong logistic damping, or stronger control of ∇w.Source: §A.5 · p. 41

B.1. Linearization about coexistence

  • The coexistence stability calculation begins by linearizing the system about a spatially homogeneous equilibrium and retaining first-order terms.The linearized chemokine equation includes diffusion, tumour and immune perturbations, and chemokine decay.Source: B. Details of the Routh–Hurwitz… · p. 41, §B.1 · p. 41
  • At the homogeneous equilibrium, the chemotactic term linearizes to a positive immune–chemokine coupling after projection onto a Neumann eigenmode.The sign of this contribution matters for the mode-wise stability calculation.Source: §B.1 · p. 41

B.2. Projection onto Neumann modes

Projection onto Neumann modes makes the chemotactic contribution vanish for the homogeneous mode while allowing it to affect finite spatial modes. The resulting cubic stability conditions show how chemotactic sensitivity can contribute to instability, although other modal crossings may set the onset.

B.2. Projection onto Neumann modes

  • Chemotaxis has no effect through this term on the homogeneous mode, but can destabilize finite spatial modes with μ_k > 0.The chemotactic contribution is proportional to μ_k, which is zero for the homogeneous mode and positive for nonzero modes.Source: §B.2 · p. 42

B.3. Characteristic polynomial

  • The modal characteristic polynomial is cubic, and each mode is linearly stable exactly when all its roots lie in the open left half-plane.The polynomial follows from the mode matrix and its expansion.Source: §B.3 · p. 42, §B.4 · p. 43

B.4. Routh–Hurwitz conditions

  • Routh–Hurwitz stability requires a1(k), a2(k), and a3(k) to be positive, with a1(k)a2(k) > a3(k).Since A_k, B_k, and C_k are positive, a1(k) is always positive; the remaining conditions encode competing damping, recruitment, and chemotactic amplification.Source: §B.4 · p. 43

B.5. Chemotactic contribution to the constant term

  • When Λ_k > 0, increasing chemotactic sensitivity decreases a3(k; ξ), so this condition can provide a chemotactic destabilization threshold.The stated threshold is useful only when Λ_k > 0 and A_k t − p s > 0; if Λ_k ≤ 0, the mode is already unstable at ξ = 0.Source: §B.5 · p. 43

B.6. Alternative instability through the determinant condition

  • Instability can also begin through a2(k) = 0 or a1(k)a2(k) − a3(k) = 0; for the Section 7 parameter set, the onset is the latter Hopf-type crossing.The explicit threshold from a3(k) is a sufficient destabilization criterion, not the only instability mechanism.Source: §B.6 · p. 44

Related papers

Loading 2607.03813v2…