Analysis and Mean-Field Limit of a Hybrid PDE-ABM Modeling Angiogenesis-Regulated Resistance Evolution

1. Introduction

Cancer remains a leading cause of global mortality, accounting for 16.8% of all deaths and 22.8% of noncommunicable diseases-related deaths [1]. Despite advances in targeted therapy, immunotherapy, and radiotherapy, therapeutic resistance persists as the primary barrier to achieving durable cures, often resulting in relapse or refractory disease [2]. Here, we focus on how resistance mechanisms—particularly those driven by tumor angiogenesis—undermine current treatments and perpetuate poor prognosis. Resistance arises from both intrinsic tumor heterogeneity (e.g., spatial and temporal genetic diversity [3]) and adaptive responses to therapy, with the tumor microenvironment (TME) playing a pivotal role. Among these factors, dysregulated angiogenesis emerges as a critical facilitator of resistance. Tumors hijack pro-angiogenic pathways (e.g., vascular endothelial growth factor/vascular endothelial growth factor receptor (VEGF/VEGFR), hypoxia-inducible factor 1$\alpha$ (HIF-1$\alpha$)) to sustain growth, while suppressing anti-angiogenic signals (e.g., thrombospondin-1 (TSP-1)) [4]. This "angiogenic switch" not only fuels tumor progression but also establishes a resilient vascular niche that evades therapy. Although anti-angiogenic agents (e.g., bevacizumab) initially improve progression-free survival by blocking VEGF, their long-term efficacy is limited by rapid resistance [5,6]. Critically, tumors activate compensatory pathways (e.g., fibroblast growth factor 2 (FGF-2), Interleukin-8 (IL-8), and Angiopoietin-2 (ANGPT-2)) to restore angiogenesis [7], while stromal components (e.g., tumor-associated macrophages (TAMs), myeloid-derived suppressor cells (MDSCs)) further sustain vascular remodeling via cytokine secretion [8]. These findings reveal a paradox: angiogenesis inhibition may inadvertently select for more aggressive, adaptive phenotypes. To overcome this challenge, we argue that disrupting the dynamic crosstalk between tumor cells and the angiogenic TME is essential. Recent studies suggest that combinatorial targeting of both VEGF and compensatory pathways (e.g., FGF-2/IL-8) may delay resistance [9]. However, a systematic framework to predict tumor adaptive responses and optimize combination regimens is still lacking. In this study, we propose a computational hybrid model to identify vulnerabilities in angiogenic networks and investigate the bidirectional coupling between angiogenesis and tumor resistance evolution. To address this need, we first contextualize existing modeling paradigms that inform such integrative frameworks.

Continuum models efficiently capture bulk tumor growth and transient dynamics. The Keller-Segel framework describes cell migration via partial differential equations (PDEs), where endothelial cells follow chemical gradients [10]. Such models assume continuous cell densities, overlooking individual cell heterogeneity. The Fisher-KPP equation [11,12] links proliferation and motility to model tumor spread as reaction-diffusion waves, but cannot resolve rare stochastic events like resistance mutations. Discrete models track individual cell behaviors and microenvironment heterogeneity. Anderson et al. [13] pioneered agent-based models (ABMs) to simulate capillary sprouting through stochastic tip-cell migration rules, but lack coupling to continuum metabolite fields (e.g., oxygen gradients) and microenvironmental feedback (e.g., hypoxia-driven angiogenesis). Despite these advances, both approaches face inherent limitations: Continuum models lack single-cell resolution for stochastic events (e.g., mutations, phenotype switching) and spatiotemporal heterogeneity, while discrete models are computationally prohibitive for large-scale dynamics (e.g., diffusion, fluid flow). This dichotomy motivates hybrid frameworks that transcend scale constraints.

Hybrid discrete-continuum (HDC) models integrate agent-based rules with reaction-diffusion equations to reflect cancer’s multiscale nature [14,15,16]. Some hybrid models address drug resistance evolution. Gevertz et al. [17] coupled continuous drug and oxygen fields with agent-based tumor dynamics to study preexisting and acquired resistance, but assumed stationary vessels, uniform resistance acquisition capacity, and continuous drug administration. Picco et al. [18] modeled environmentally mediated drug resistance (EMDR) using drug and signaling fields, generating EMDR by increasing vascular density but neglecting vascular remodeling. Other hybrid models coupled tumor growth with angiogenesis [19,20]. Macklin et al. [21] developed a multiscale framework integrating extracellular matrix, angiogenesis, and tumor progression. Despite these innovations, a critical gap persists: Most models simulate vascular remodeling or resistance in isolation, omitting spatiotemporal feedback where vessel adaptation alters resistance trajectories. One exception is Wang et al. [22], which proposed a hybrid model coupling hypoxia-induced angiogenesis with resistance selection, but left unresolved three limitations: (i) inadequate mathematical analysis of the hybrid system, insufficient numerical scheme analysis, and absence of derived mean-field limits for agent dynamics. Our work directly addresses these gaps while extending their framework.

We advance the hybrid paradigm in [22] through three innovations: First, we introduced a spatially resolved model unifying stochastic phenotypic switching via a Gillespie algorithm to capture drug-induced resistance transitions. Second, we incorporate metabolite-regulated angiogenesis where vessel agents respond to dynamic hypoxia (o) and tumor angiogenic factor (TAF) (c) gradients. Third, we establish bidirectional coupling enabling emergent feedback between hypoxia, angiogenesis, and resistance evolution. Consequently, this framework yields three contributions: rigorous mathematical analysis of the coupled PDE-ABM system and numerical analysis of PDE discretization schemes; integration of exact stochastic resistance (Gillespie ABM) with physics-driven angiogenesis (PDE-coupled ABM); and derivation of analytically tractable continuum mean-field limits from discrete rules. The mean-field limit integrates the hybrid system into one PDE system, linking stochastic microdynamics with deterministic macrodynamics. It motivates innovative hybrid numerical schemes that selectively retain the full ABM structure in critical regions while replacing it with efficient PDE surrogates elsewhere. By unifying PDE analysis, stochastic analysis, and numerical analysis, our framework moves beyond descriptive modeling to establish rigorous mathematical foundations for hybrid PDE-ABM systems.

The paper is structured as follows: Section 2 details the hybrid formulation. Section 3 links our model to classical mathematical biology frameworks and outlines contributions. Section 4 provides mathematical analysis of the coupled system. Section 5 presents numerical analysis, including consistency, stability, convergence, and conservation properties. Section 6 concludes.

2. Material & Methods

Building on the Introduction’s biological foundations, we formalize a hybrid PDE-ABM with bidirectional coupling between continuum fields (TAF, drug, oxygen) and discrete agents (tumor and vessel cells). Reaction-diffusion equations govern TAF (c), drug (d), and oxygen (o) dynamics. Endothelial cell (n) migration follows increasing TAF concentrations via the chemotaxis equation:

$$\begin{array}{c}\hfill {\partial }_{t}n=D_n\mathrm{\Delta }n-\nabla ·\left({\displaystyle \frac{{\chi }_0}{1+\alpha c}}n\nabla c\right).\end{array}$$

Discrete tumor and vessel agents interact with the TME through stochastic rule-based ABM dynamics. We denote tumor, normoxic, hypoxic, tip, and vessel agent sets at time t as ${\mathrm{\Lambda }}_t$, ${\mathrm{\Lambda }}_t^n$, ${\mathrm{\Lambda }}_t^h$, $T_t$, and $V_t$, respectively. Local oxygen levels partition the tumor population into normoxic (${\mathrm{\Lambda }}_t^n$) and hypoxic (${\mathrm{\Lambda }}_t^h$) subpopulations. Full model details appear in [22]. We present the non-dimensionalized system for $(n,c,d,o)$; the original dimensional system and non-dimensionalization procedure derive from Appendix A.

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\partial }_{t}n=D_n\mathrm{\Delta }n-\nabla ·\left({\displaystyle \frac{{\chi }_0}{1+\alpha c}}n\nabla c\right),\hfill \\ {\partial }_{t}c=D_c\mathrm{\Delta }c-{\xi }_{c}c+\eta \sum _{a\in {\mathrm{\Lambda }}_t^h}{\chi }_a-\lambda c\sum _{v\in V_t}{\chi }_v,\hfill \\ {\partial }_{t}d=D_d\mathrm{\Delta }d-{\xi }_{d}d-{\rho }_{d}d\sum _{a\in {\mathrm{\Lambda }}_t}{\chi }_a+S_d\sum _{v\in V_t}{\chi }_v,\hfill \\ {\partial }_{t}o=D_o\mathrm{\Delta }o-{\xi }_{o}o-{\rho }_o\sum _{a\in {\mathrm{\Lambda }}_t}{\chi }_a+S_o\left(1-o\right)\sum _{v\in V_t}{\chi }_v.\hfill \end{array}\right.\end{array}$$

(2.1)

Here, $D_n$, $D_c$, $D_d$, $D_o$ are diffusion coefficients; ${\chi }_0$ is the chemotactic sensitivity; $\alpha$ is the saturation parameter; ${\xi }_c$, ${\xi }_d$, ${\xi }_o$ are decay rates; ${\rho }_d$, ${\rho }_o$ are cellular uptake rates; $S_d$, $S_o$ are supply rates from blood vessels; and $\eta$, $\lambda$ are production/uptake rates of TAF. ${\chi }_a$, ${\chi }_v$ are characteristic functions that equal 1 within disks of cellular radius $R_c$:

$$\begin{array}{c}\hfill {\chi }_a(x,t)=\left\{\begin{array}{cc}1\hfill & if∥x-a^{\left(X,Y\right)}\left(t\right)∥\le R_c,\hfill \\ 0\hfill & otherwise,\hfill \end{array}\right.and{\chi }_v(x,t)=\left\{\begin{array}{cc}1\hfill & if∥x-v^{\left(X,Y\right)}\left(t\right)∥\le R_c,\hfill \\ 0\hfill & otherwise.\hfill \end{array}\right.\end{array}$$

Here, $a^{(X,Y)},v^{(X,Y)}$ denote the position center of a and v. We choose a closed domain $U\subset {\mathbb{R}}^2$ with $C^{\infty }$ boundary $\partial U$ and impose homogeneous Neumann conditions:

$$\begin{array}{c}\hfill \overrightarrow{n}{·\nabla u|}_{\partial U}=0,\end{array}$$

where $\overrightarrow{n}$ is the outward unit normal at $\partial U$.

Tumor and vessel agents evolve on a lattice ${\widehat{U}}_{\mathrm{ABM}}$ discretizing the ABM domain $U_{\mathrm{ABM}}={[0,1]}^2$ into $N^2$ squares of side length $\mathrm{\Delta }x=1/N$. Tip cell motion uses the Von Neumann neighborhood structure (four orthogonal neighbors), while agent-microenvironment interactions use the Moore structure (eight neighbors including diagonals).

We construct a PDE domain $U_{\mathrm{PDE}}\subset {\mathbb{R}}^2$ with $C^{\infty }$ boundary approximating $U_{\mathrm{ABM}}$. Appendix B details the mapping $U_{\mathrm{ABM}}\mapsto U_{\mathrm{PDE}}$ and the numerical derivation of the discrete lattice ${\widehat{U}}_{\mathrm{ABM}}$ in the ABM. Briefly, define the open $\delta$-neighborhood $U_{\delta }$ of $U_{\mathrm{ABM}}$ with $\delta <\frac{1}{6}\mathrm{\Delta }x$. Convolve its characteristic ${\chi }_{U_{\delta }}$ with a compactly supported $C^{\infty }$ mollifier ${\psi }_{\delta }$, and set

$$\begin{array}{c}\hfill U_{\mathrm{PDE}}=\left\{(x,y)\in {\mathbb{R}}^2:({\chi }_{U_{\delta }}*{\psi }_{\delta })(x,y)>1-ϵ\right\},\end{array}$$

for small $ϵ>0$ ensuring $U_{\mathrm{PDE}}\subset U_{\delta /6}$. This yields a domain $\partial U_{\mathrm{PDE}}$ whose $C^{\infty }$ boundary $\partial U_{\mathrm{PDE}}$ is proximal to $\partial U_{\mathrm{ABM}}$. We obtain the discrete approximation ${\widehat{U}}_{\mathrm{ABM}}$ by sampling $\mathrm{\Psi }=({\chi }_{U_{\delta }}*{\psi }_{\delta })$ on the ABM lattice, thresholding at level $1-ϵ$, and identifying boundary agents as those with Moore neighbors outside $U_{\mathrm{PDE}}$. Algorithm A1 summarizes the construction.

Each tumor cell $a\in {\mathrm{\Lambda }}_t$ has state:

$$\begin{array}{c}\hfill a=\left\{a^{(X,Y)},a^o,a^d,a^{\mathrm{dam}},a^{\mathrm{death}},a^{\mathrm{age}},a^{\mathrm{mat}}\right\},\end{array}$$

denoting position, local oxygen, accumulated drug, DNA damage, death threshold, age, and maturation time. Each tip cell b has state:

$$\begin{array}{c}\hfill b=\{b^{(X,Y)}\left(t\right),b^{\mathrm{age}}\left(t\right)\},\end{array}$$

for position and age since the last branching. Agent states update via Markov processes under the time discretization $0=t_0<t_1<t_2<\dots <t_M=T$ of $[0,T]$:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}x\left(t_{k+1}\right)={\mathcal{F}}_{x,k}(c(·,t_k),d(·,t_k),o(·,t_k),\mathcal{M}\left(t_k\right)),\hfill \\ x\left(t\right)=x\left(t_{k+1}\right),t\in [t_k,t_{k+1}),\hfill \end{array}\right.\end{array}$$

(2.2)

where ${\mathcal{F}}_{x,k}$ is the Markov transition operator [22]. Algorithm A2 summarizes the transition operator ${\mathcal{F}}_{x,k}$ for each mark x. for tumor and endothelial tip cells. Local cell density includes tumor-vessel interactions in our current model:

$$\begin{array}{c}\hfill F(x,t)\sum _{a\in {\mathrm{\Lambda }}_t}{\chi }_{B_{R_c}}\left(x-a^{(X,Y)}\left(t\right)\right)+\sum _{v\in V_t}{\chi }_{B_{R_c}}\left(x-v^{(X,Y)}\left(t\right)\right)\end{array}$$

(2.3)

We decouple phenotype switching from cell division. Unlike classical mutation models linking resistance acquisition to division events, stochastic $S\rightleftharpoons R$ transitions occur through continuous-time Markov processes independent of division cycles. Mutation elevates the death threshold, DNA repair rate, and oxygen consumption rate but reduces the proliferation rate. Phenotypic traits

$$\begin{array}{c}\hfill x\in \{\mathrm{death}\mathrm{threshold},\mathrm{repair}\mathrm{rate},\mathrm{oxygen}\mathrm{consumption}\mathrm{rate},\mathrm{proliferation}\mathrm{rate}\}\end{array}$$

take sensitive ($x_S$) and resistant ($x_R$) values:

$$\begin{array}{c}\hfill x_S:\mathrm{drug}-\mathrm{sensitive}\mathrm{trait}\mathrm{value},\\ \hfill x_R:\mathrm{drug}-\mathrm{resistant}\mathrm{trait}\mathrm{value}.\end{array}$$

Cells stochastically switch between these states according to the reaction channels

$$\begin{array}{cc}\hfill & x_S\stackrel{k_{SR}}{⟶}{x}_R\left(\mathrm{induction}\mathrm{of}\mathrm{resistance}\right),\hfill \\ \hfill & x_R\stackrel{k_{RS}}{⟶}{x}_S\left(\mathrm{reversion}\mathrm{to}\mathrm{sensitivity}\right),\hfill \end{array}$$

where baseline rates $k_1$ and $k_2$, drug sensitivity coefficients ${\alpha }_d$ and ${\beta }_d$, sensitive-to-resistant mutation rate $k_{SR}$ and resistant-to-sensitive mutation rate $k_{RS}$

$$\begin{array}{c}\hfill k_{SR}(x,t)=k_1+{\alpha }_{d}d(x,t),k_{RS}(x,t)=max\{0,k_2-{\beta }_{d}d(x,t)\}\end{array}$$

govern phenotypic switching. Each tumor cell $a_i\in {\mathrm{\Lambda }}_t$ bears a phenotypic mark ${\sigma }_i\in \{S,R\}$, updated independently of the cell cycle via the Gillespie algorithm [23]. Waiting times follow an exponential distribution with mean $1/A_0$ ($A_0$: total propensity):

$$\begin{array}{c}\hfill \mathrm{\Delta }t={\displaystyle \frac{1}{A_0}}log\left({\displaystyle \frac{1}{r_1}}\right),r_1\sim \mathrm{Uniform}(0,1).\end{array}$$

Algorithm 1 Stochastic phenotype switching via the Gillespie Direct Method

Require: Initial tumor agents ${\mathrm{\Lambda }}_0$, positions $a_i^{(X,Y)}$, phenotypes ${\sigma }_i\in \{S,R\}$, mutation rates $k_{SR}(x,t)$, $k_{RS}(x,t)$, and ABM time horizon $T_{\mathrm{ABM}}$. 1: Set $t\leftarrow 0$. 2: while$t<T_{\mathrm{ABM}}$do 3:     for each agent $a_i\in {\mathrm{\Lambda }}_t$ do 4:         Compute local mutation rates: $k_{SR}^i\leftarrow k_{SR}(a_i^{(X,Y)}\left(t\right),t)$, $k_{RS}^i\leftarrow k_{RS}(a_i^{(X,Y)}\left(t\right),t)$. 5:         Set propensity $A_i\leftarrow k_{SR}^i$, if ${\sigma }_i=S$; $A_i\leftarrow k_{RS}^i$ if ${\sigma }_i=R$. 6:     end for 7:     Compute total propensity $A_0\leftarrow {\sum }_{i=1}^{|{\mathrm{\Lambda }}_t|}{A}_i$. 8:     if $A_0=0$ then 9:         break ▹ No further switches possible 10:     end if 11:     Draw $r_1,r_2\stackrel{\mathrm{i}.\mathrm{i}.\mathrm{d}.}{\sim }\mathrm{Uniform}(0,1)$. 12:     Compute exponential waiting time $\mathrm{\Delta }t\leftarrow {\displaystyle \frac{1}{A_0}}log\left({\displaystyle \frac{1}{r_1}}\right)$. 13:     Find minimal j satisfying ${\sum }_{i=1}^j{A}_i\ge r_2{A}_0$. 14:     Flip ${\sigma }_j$: ${\sigma }_j\leftarrow R$ if ${\sigma }_j=S$; ${\sigma }_j\leftarrow S$ if ${\sigma }_j=R$. 15:     Update $t\leftarrow t+\mathrm{\Delta }t$. 16:     Update: (i) local drug $d(a_i^{(X,Y)}\left(t\right),t)$, (ii) mutation rates $k_{SR}^i$, $k_{RS}^i$, (iii) ${\mathrm{\Lambda }}_t$ and $\left\{{\sigma }_i\right\}$, (iv) phenotype-associated traits: death threshold, repair rate, oxygen consumption rate, proliferation rate. 17: end while

Each switching event selects the reacting agent by propensities, flips its phenotype ($S\leftrightarrow R$), and updates time-dependent quantities. Algorithm 1 encapsulates phenotype switching procedures. This Gillespie framework permits resistance emergence independent of cell division, capturing stochastic gene expression fluctuations. It thus recovers resistance in rare subpopulations that deterministic ordinary differential equation (ODE)/PDE models average out, and naturally quantifies treatment schedule effects on resistance probability and timing. Table 1 summarizes all parameters in our hybrid PDE-ABM, categorized by continuum PDEs, discrete ABM rules, and numerical implementation.

Our model uniquely captures bidirectional feedback: how hypoxia induces angiogenesis, which alters drug delivery and selects for resistance. The theoretical framework justifies how small changes in local oxygen or vessel density can produce nonlinear shifts in resistance trajectories. These feedback mechanisms align with clinical observations where anti-angiogenic therapies may paradoxically increase resistance risk. Thus, the model could be used to explore: (i) critical thresholds in vessel remodeling that shift tumor response from chemosensitive to chemoresistant, and (ii) timing effects, e.g., whether anti-angiogenic pre-treatment amplifies or delays resistant onset.

3. Connection to Classical Mathematical Biology Frameworks and Contributions

Having formalized our multiscale framework, we contextualize its biological relevance through its mathematical lineage. This section demonstrates that limiting cases of our hybrid PDE-ABM system reduce rigorously to three foundational frameworks: Keller-Segel chemotaxis, Fisher-KPP invasion waves, and Galton-Watson branching processes. These connections are mathematically exact. For a small chemokine concentration c such that $1+\alpha c\approx 1$, the endothelial cell equation simplifies to the classical Keller-Segel form:

$$\begin{array}{c}\hfill {\partial }_{t}n=D_n\mathrm{\Delta }n-{\chi }_0\nabla ·(n\nabla c).\end{array}$$

Our model generalizes Keller-Segel by incorporating nonlinear chemotactic saturation and ABM coupling.

While endothelial migration follows chemotactic principles, tumor proliferation exhibits distinct spatiotemporal dynamics. This divergence aligns with the Fisher-KPP paradigm, governed by the equation:

$$\begin{array}{c}\hfill {\partial }_t\rho =D\mathrm{\Delta }\rho +\alpha \rho \left(1-{\displaystyle \frac{\rho }{K}}\right),\end{array}$$

where $\rho (x,t)$ denotes density, D diffusion, $\alpha$ growth rate, and K carrying capacity. Under normoxic conditions ($o\approx o_{max}$) and negligible drug ($d\to 0$), tumor density equation (Equation (4.5) reduces to

$$\begin{array}{c}\hfill {\partial }_t\rho =D\mathrm{\Delta }\rho +{\alpha }_n\rho \left(1-{\displaystyle \frac{\rho }{F_{max}}}\right).\end{array}$$

This predicts traveling wavefronts with minimal speed $2\sqrt{D{\alpha }_n}$. Crucially, our full model extends this framework by incorporating microenvironmental heterogeneity (oxygen gradients, drug distribution), which distorts wavefronts into clinically observed asymmetric invasion patterns.

Complementing continuum-scale dynamics, cellular-scale proliferation exhibits stochastic foundations. Absent spatial coupling and microenvironmental feedback (e.g., well-mixed conditions), tumor lineage expansion follows a Galton-Watson branching process. A normoxic cell divides with probability $p_d$ or dies with probability $p_{death}$ in time interval $\mathrm{\Delta }t$. yielding the offspring distribution:

$$\begin{array}{c}\hfill \mathbb{P}(\mathrm{offspring}=2)=p_d,\mathbb{P}(\mathrm{offspring}=0)=p_{death},\mathbb{P}(\mathrm{offspring}=1)=1-p_d-p_{death}.\end{array}$$

For two phenotypes (drug-sensitive S, drug-resistant R), this extends to a two-type Galton-Watson process with expected offspring matrix:

$$\begin{array}{c}\hfill M=\left(\begin{array}{cc}\mathbb{E}\left[X^{(S\to S)}\right]& \mathbb{E}\left[X^{(S\to R)}\right]\\ \mathbb{E}\left[X^{(R\to S)}\right]& \mathbb{E}\left[X^{(R\to R)}\right]\end{array}\right).\end{array}$$

Crucially, classical branching models restrict phenotypic switching to division events. In contrast, our framework implements phenotype switching via a continuous-time Markov process simulated using the Gillespie algorithm (Algorithm 1), enabling asynchronous phenotype evolution independent of the cell cycle. Figure 1 illustrates interactions between branching proliferation and stochastic phenotype switching.

Our mathematically rigorous, biologically grounded HDC framework extends classical theories through multiscale feedback, stochastic switching, and microenvironmental coupling. It advances hybrid modeling by integrating agent-based stochastic phenotypic switching, continuum reaction-diffusion metabolite fields, and bidirectional vascular-resistance coupling. We provide analytical guarantees for model well-posedness, numerical robustness, and mean-field limits. Our contributions distinguish themselves from existing oncology hybrid models (e.g., [39,40,41]) in four key aspects:

(i)

Mean-Field Limits: Unlike direct continuum PDE-ABM coupling or variance reduction methods [39,40], we prove convergence from the stochastic agent-based description to a deterministic PDE limit using chemical master equations (CMEs) and moment closures. This yields a unified continuum system (4.6), linking stochastic microdynamics to deterministic macrodynamics. It motivates innovative hybrid numerical schemes retaining ABM structures in critical regions while employing efficient PDE surrogates elsewhere.

(ii)

Well-Posedness Analysis: While most stochastic-hybrid models prioritize computation over analytical soundness [40,41], we establish complete well-posedness for our hybrid PDE-ABM system. This includes the existence, uniqueness, and regularity of weak solutions (Theorems 1-3), plus analysis of numerical schemes: consistency (Section 5.2), conditional stability (Theorem 6), convergence (Theorem 7), nonnegativity (Theorem 8) and mass conservation (Theorem 9).

(iii)

Stochastic-Deterministic Feedback: Unlike [39], who model virus spread in lung tissue, our hybrid system couples cell phenotype switching, angiogenic fields, and vascular remodeling. This reveals nonlinear resistance dynamics driven by environment-phenotype feedback, underpinned by rigorous mathematics and mean-field approximations.

(iv)

Clinical interpretability: While stochastic interacting particle fields (SIPFs) [41] or variance reduction techniques [40] emphasize algorithmic efficiency, our model explicitly links parameters (e.g., mutation rates, oxygen thresholds, drug uptake) to clinical outcomes (e.g., tumor front speed, resistance emergence). This offers a framework for optimizing treatment timing and dosage.

Thus, our work uniquely combines mathematical rigor (well-posedness, mean-field limits) with biologically interpretable hybrid modeling, delivering methodological novelty and translational potential beyond current stochastic-hybrid cancer models.

4. Mathematical Analysis of the Coupled PDE-ABM System & Results

Biological motivation supports the hybrid model design, but rigorous mathematical analysis remains essential to establish its validity and predictive reliability. This section systematically investigates the coupled PDE-ABM framework through three analytical pillars: (i) well-posedness of the PDE-ABM system (Section 4.1), guaranteeing solution existence, uniqueness, and regularity; (ii) master equation and mean-field analysis (Section 4.2), connecting discrete stochastic dynamics to deterministic continuum-scale approximations; and (iii) validity regimes for moment closure approximations, linking microscopic randomness with macroscopic behavior. This triad ensures mathematical robustness for therapeutic applications, including chemotherapy optimization and anti-angiogenic interventions.

4.1. Well-Posedness

We establish existence, uniqueness, and regularity of weak solutions for the coupled PDE-ABM system governing the continuous fields $\mathbf{u}={(c,d,o)}^{\top }$. Without agent interactions, the PDEs reduce to a classical linear parabolic system with guaranteed well-posedness. We formalize this foundation using the Lions–Gelfand triple framework [Theorem 10.9 [42]; see [43] for detailed proofs.

Lemma 1

(J.-L. Lions). Let $V\subset H\subset V^{*}$ be a Gelfand triple, with V and H Hilbert spaces and V densely and continuously embedded into H. The scalar product and norm for H are $(·,·)$ and $|·|$, and the norm for V is $\parallel ·\parallel$. Let $T>0$. Suppose for almost every $t\in [0,T]$, $a(t;u,v):V\times V\to \mathbb{R}$ is a time-dependent bilinear form satisfying:

(i)

For all $u,v\in V$, the function $t\mapsto a(t;u,v)$ is measurable.

(ii)

$\left|a\right(t;u,v\left)\right|\le M\parallel u\parallel \parallel v\parallel$ for almost every $t\in [0,T]$ and all $u,v\in V$.

(iii)

$a(t;v,v)\ge {\alpha \parallel v\parallel }^2-C{\left|v\right|}^2$ for almost every $t\in [0,T]$ and all $v\in V$.

where $\alpha >0$, M, and C are constants. Then for $f\in L^2(0,T;V^{*})$ and $u_0\in H$, there exists a unique function

$$\begin{array}{c}\hfill u\in L^2(0,T;V)\cap C([0,T];H),{\partial }_{t}u\in L^2(0,T;V^{*})\end{array}$$

satisfying the weak formulation:

$$\begin{array}{c}\hfill 〈{\partial }_{t}u\left(t\right),v〉+a(t;u\left(t\right),v)=〈f\left(t\right),v〉\mathrm{for}a.e.t\in (0,T),\forall v\in V,\end{array}$$

with initial condition $u\left(0\right)=u_0$.

We now apply this result to the coupled reaction-diffusion equations on the finite time interval $t\in [0,T]$:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\partial }_t\mathbf{u}=\mathcal{D}\mathrm{\Delta }\mathbf{u}+\mathcal{G}(x,t,\mathbf{u}),x\in U,t\in (0,T],\hfill \\ \mathbf{u}(x,0)={\mathbf{u}}_0\left(x\right)=(c_0\left(x\right),d_0\left(x\right),o_0\left(x\right)),x\in U,\hfill \\ {\partial }_n\mathbf{u}=0,x\in U,t\in (0,T],\hfill \end{array}\right.\end{array}$$

(4.1)

where $\mathcal{D}=\mathrm{diag}(D_c,D_d,D_o)$ and $\mathcal{G}={(G_1,G_2,G_3)}^{\top }$ incorporates agent-dependent reaction terms:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{G}_1(x,t,\mathbf{u})=-{\xi }_{c}c+\eta \sum _{a\in {\mathrm{\Lambda }}_t^h}{\chi }_a-\lambda c\sum _{v\in V_t}{\chi }_v,\hfill \\ G_2(x,t,\mathbf{u})=-{\xi }_{d}d-{\rho }_{d}d\sum _{a\in {\mathrm{\Lambda }}_t}{\chi }_a+S_d\sum _{v\in V_t}{\chi }_v,\hfill \\ G_3(x,t,\mathbf{u})=-{\xi }_{o}o-{\rho }_o\sum _{a\in {\mathrm{\Lambda }}_t}{\chi }_a+S_o\left(1-o\right)\sum _{v\in V_t}{\chi }_v.\hfill \end{array}\right.\end{array}$$

We omit the analysis of the endothelial density equation:

$$\begin{array}{c}\hfill {\partial }_{t}n=D_n\mathrm{\Delta }n-\nabla ·(\chi \left(c\right)n\nabla c).\end{array}$$

Classical two-dimensional Keller–Segel systems exhibit chemotactic blow-up, that is, unbounded density in finite time [44,45], but admit global solutions under small-mass conditions. In particular, consider the following parabolic-elliptic Keller-Segel system:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\partial }_{t}n-\kappa \mathrm{\Delta }n+\nabla ·(\chi n\nabla c)=0,\hfill & x\in U\subset {\mathbb{R}}^2,t>0,\hfill \\ -\mathrm{\Delta }c=n-\langle n\rangle ,\hfill & x\in U,t>0,\hfill \\ n(x,0)=n_0\left(x\right),\hfill & x\in U,\hfill \\ \overrightarrow{n}·\nabla n=\overrightarrow{n}·\nabla c=0,\hfill & x\in \partial U,t\ge 0,\hfill \end{array}\right.\end{array}$$

(4.2)

where $\overrightarrow{n}$ is the outward unit normal to $\partial U$ and $\langle n\rangle ={\displaystyle \frac{1}{\left|U\right|}}{\int }_{U}n(x,t)\mathrm{d}x$ denotes the spatial average of n. Assume the initial mass satisfies:

$$\begin{array}{c}\hfill {\int }_U{n}_0\left(x\right)\mathrm{d}x<{\displaystyle \frac{C_{\mathrm{opt}}^{*}\kappa }{\chi }},C_{\mathrm{opt}}^{*}=\left\{\begin{array}{cc}8\pi ,\hfill & \mathrm{if}U={\mathbb{R}}^2,\hfill \\ 4\pi ,\hfill & \mathrm{if}U\mathrm{is}\mathrm{a}\mathrm{bounded},\mathrm{connected}{C}^2\mathrm{domain}.\hfill \end{array}\right.\end{array}$$

Then, a global-in-time weak solution exists under the following conditions:

(i)

Bounded domain case: If U is a bounded, connected $C^2$ domain and $n_0\in L^{\infty }\left(U\right)$,

(ii)

Whole space case: If $U={\mathbb{R}}^2$, with ${(1+|x|}^2)n_0\in L^1\left({\mathbb{R}}^2\right)$ and $n_{0}log{n}_0\in L^1\left({\mathbb{R}}^2\right)$.

This follows from [44,46].1

Discrete agents represent the endothelial population in our hybrid formulation, replacing continuous density and further regularizing dynamics. Each tip occupies at most one lattice site. Vessel elongation occurs with a fixed period $\tau$, and branching requires the tip cell age to exceed the threshold $\psi$ ($b^{\mathrm{age}}>\psi$). Anastomosis stabilizes the system: tip encounters cause emerging, reducing the tip count. These rules ensure a finite agent population $|T_t|$ and prevent uncontrolled mass concentration, consistent with lattice-based chemotaxis ABMs [13,28,47,48]. We thus focus our analysis on the $(c,d,o)$ subsystem. To manage PDE-ABM coupling, we introduce:

Assumption A1

(PDE-ABM Coupling).

(i)

The spatial domain $U\subset {\mathbb{R}}^2$ is bounded with $C^{\infty }$ boundary $\partial U$.

(ii)

The initial data ${\mathbf{u}}_0=(c_0,d_0,o_0)$ is nonnegative and belongs to ${[L^2\left(U\right)\cap L^{\infty }\left(U\right)]}^3$.

(iii)

For all $t\in [0,T]$ and lattice sites $x\in {\widehat{U}}_{\mathrm{ABM}}$:

$$\begin{array}{c}\hfill F(x,t)\sum _{a\in {\mathrm{\Lambda }}_t}{\chi }_{B_{R_c}}(x-a^{(X,Y)}\left(t\right))+\sum _{v\in V_t}{\chi }_{B_{R_c}}(x-v^{(X,Y)}\left(t\right))\le F_{max}.\end{array}$$

Assumption A1 (iii) imposes a uniform local density bound, not global agent boundedness, as tumor proliferation may increase $|{\mathrm{\Lambda }}_t|+|V_t|$ over time. This aligns with the biological crowding constraint $F_{max}$ in ABM rules: lattice sites accommodate finitely many agents within radius $R_c$, and cell division halts at local density $F_{max}$. Consequently, global agent growth does not compromise the local boundedness of PDE source terms, which is essential for coupled system well-posedness.

We define the weak formulation by seeking $\mathbf{u}={(c,d,o)}^{\top }$ satisfying $u\in L^2(0,T;H^1\left(U\right))$ and ${\partial }_{t}u\in L^2(0,T;H^1{\left(U\right)}^{\prime })$ for each $u\in \{c,d,o\}$, with:

$$\begin{array}{c}\hfill \langle {\partial }_{t}u,\varphi \rangle +D_u{(\nabla u,\nabla \varphi )}_{L^2\left(U\right)}={(G_u,\varphi )}_{L^2\left(U\right)},\forall \varphi \in H^1{\left(U\right)}^3,a.e.,t\in (0,T),\end{array}$$

and initial conditions $u(x,0)=u_0\left(x\right)$. Here ${(·,·)}_{L^2\left(U\right)}$ denotes the $L^2$ inner product, and $\langle ·,·\rangle$ is the $H^1\left(U\right)$ duality pairing.

Theorem 1

(Existence and Uniqueness of Weak Solutions). Assumption A1 guarantees unique weak solutions for the coupled PDE-ABM system (4.1) with:

$$\begin{array}{c}\hfill u\in L^2(0,T;H^1\left(U\right))\cap C([0,T];L^2\left(U\right)),{\partial }_{t}u\in L^2(0,T;H^1{\left(U\right)}^{\prime }),u\in \{c,d,o\}.\end{array}$$

Proof. 

Apply Lemma 1 to each component via the Gelfand triple $H^1\left(U\right)\subset L^2\left(U\right)\subset H^1{\left(U\right)}^{\prime }$ and coercive, bounded bilinear forms $a_u(t;v,w)=D_u(\nabla v,\nabla w)$ for $u\in \{c,d,o\}$.    □

We rigorously couple ABM and PDE dynamics through time discretization $0=t_0<t_1<t_2<\dots <t_n=T$ and operator splitting:

Theorem 2

(Well-posedness of the Coupled PDE-ABM System). Under Assumption A1 on $[0,T]$ with time steps $t_k$:

(i)

The PDE subsystem admits a unique weak solution with Theorem 1 regularity.

(ii)

The ABM subsystem defines càdlàg Markov process almost surely for agent marks

$$\begin{array}{c}\hfill \mathcal{M}=\{a^{(X,Y)},a^o,a^d,a^{\mathrm{dam}},a^{\mathrm{death}},a^{\mathrm{age}},a^{\mathrm{mat}},b^{(X,Y)},b^{\mathrm{age}}\},\end{array}$$

Proof. 

Use mathematical induction over intervals $I_k=[t_{k-1},t_k]$. For $k=1$, Theorem 1 establishes existence, uniqueness, and regularity on $I_1=[0,t_1]$ with initial data $(c_0,d_0,o_0)$ and fixed agent $({\mathrm{\Lambda }}_{t_0},V_{t_0})$. Assume the result holds on $I_k=[t_{k-1},t_k]$. For $I_{k+1}=[t_k,t_{k+1}]$:

(i)

Solve the parabolic subsystem on $I_{k+1}^r=(t_k,t_{k+1}]$:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\partial }_t\mathbf{u}=\mathcal{D}\mathrm{\Delta }\mathbf{u}+\mathcal{G}(x,t_k,\mathbf{u}),\hfill & x\in U,t\in I_{k+1}^r,\hfill \\ \mathbf{u}(x,t_k)=(c(x,t_k),d(x,t_k),o(x,t_k)),\hfill & x\in U,\hfill \\ {\partial }_n\mathbf{u}=0,\hfill & x\in \partial U,t\in I_{k+1}^r,\hfill \end{array}\right.\end{array}$$

(4.3)

which has unique weak solutions by Theorem 1.

(ii)

Evolve ABM marks on $I_{k+1}^l=[t_k,t_{k+1})$ via:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}x\left(t_{k+1}\right)={\mathcal{F}}_{x,k+1}(c(·,t_k),d(·,t_k),o(·,t_k),\mathcal{M}\left(t_k\right)),\hfill \\ x\left(t\right)=x\left(t_{k+1}\right),t\in I_{k+1}^l,\hfill \end{array}\right.\end{array}$$

where ${\mathcal{F}}_{x,k+1}$ is the Markov transition operator. Piecewise constant agent trajectories with right-continuous and left limits ensure the càdlàg Markov property.

Induction extends the result to $[0,T]$.    □

Remark 1

(Justification of Stepwise Decoupling). Theorem 2 solves the PDE subsystem on each interval $I_{k+1}^r=(t_k,t_{k+1}]$ with fixed agent configuration $({\mathrm{\Lambda }}_{t_k},V_{t_k})$. This approach directly implements the hybrid model’s operator–splitting scheme: PDE fields evolve continuously via (4.3) on $(t_k,t_{k+1}]$. While agent updates occur discretely at $t_k$. Consequently, PDE source terms $\mathcal{G}(x,t_k,\mathbf{u})$ are piecewise constant in time, preserving well-posedness as a linear parabolic problem on each subinterval. The ABM update (2.2) defines a Markov map at mesh points.

Since the alternating-update rule defines the hybrid model itself, no splitting error accumulates. The piecewise-constant ABM forcing is an exact model component. The induction argument formalizes this sequential construction, guaranteeing existence and uniqueness on each $I_k$ and hence on $[0,T]$.

We examine the decoupled system to highlight the ABM’s role. The resulting linear parabolic PDEs with Neumann boundary conditions admit a classical solution via spectral decomposition of the Neumann Laplacian (proof in [49]):

Theorem 3

(Well-Posedness of Sole PDE System). The system

$$\begin{array}{c}\hfill {\partial }_{t}u=D_u\mathrm{\Delta }u-{\xi }_{u}u\mathrm{in}U{\partial }_n\mathbf{u}=0\mathrm{on}\partial Uu(x,0)=u_0\left(x\right),\end{array}$$

(4.4)

admits a unique classical solution:

$$\begin{array}{c}\hfill u(x,t)=\sum _{k=1}^{\infty }{a}_k^u{e}^{-(D_u{\lambda }_k+{\xi }_u)t}{\omega }_k,\end{array}$$

where the Neumann Laplacian eigenbasis $\left\{{\omega }_k\right\}$ satisfies:

$$\begin{array}{c}\hfill -\mathrm{\Delta }{\omega }_k={\lambda }_k{\omega }_k\mathrm{in}U{\partial }_n{\omega }_k=0\mathrm{on}\partial U,{\omega }_k\in H^1\left(U\right)\cap C^{\infty }\left(\overline{U}\right),\end{array}$$

with eigenvalues $0={\lambda }_1<{\lambda }_2\le {\lambda }_3\le \dots \to +\infty$. Coefficients are $a_k^u={(u_0,{\omega }_k)}_{L^2\left(U\right)}$. Each $u\in \{c,d,o\}$ is spatially smooth for $t>0$ and temporally analytic for $x\in U$.

The spectral result confirms that without ABM interactions, the PDE subsystem maintains full regularity and decouples completely. In the coupled system, agent density constraints ensure well-posedness despite nonlinear agent-driven reaction terms. Thus, our hybrid framework is both analytically well-posed and biologically interpretable, providing a solid foundation for clinical applications.

Weak solution existence establishes a rigorous mathematical foundation for hybrid models integrating microscopic interactions with macroscopic behaviors. The proof of well-posedness for the coupled PDE-ABM system guarantees that the model behaves predictably under biologically relevant conditions, e.g., bounded agent density due to cell crowding and realistic spatial constraints. This is critical for ensuring that simulated phenomena like angiogenesis, hypoxia-driven resistance, and spatial tumor heterogeneity reflect biological processes rather than numerical artifacts. From a clinical perspective, this mathematical well-posedness lays the foundation for future optimization of treatment regimens (e.g., dosing schedules) within a consistent modeling framework.

4.2. Master Equation Formulation and Mean-Field Limit

Having established deterministic well-posedness, we now analyze stochastic aspects of agent dynamics. Integrating stochastic agent-based dynamics within a continuum PDE framework requires principled linkages between microscopic randomness and macroscopic observables. The CME provides this foundation by capturing the probabilistic evolution of the agent state space. This section derives a lattice-based CME for the hybrid PDE-ABM system, computes first-order moment equations, and establishes conditions justifying a deterministic mean-field approximation. By unifying the hybrid PDE-ABM into a continuous description (4.6), our model connects microscale interactions with macroscale dynamics, replacing the computationally demanding ABM structure in dense tumor regions with an efficient PDE solver.

We analyze a coarse-grained mesoscopic representation of the ABM component. Consider a regular lattice ${\widehat{U}}_{\mathrm{ABM}}=\{x_1,x_2,\dots ,x_{N^2}\}$ with spacing $\mathrm{\Delta }x={\displaystyle \frac{1}{N}}>0$ over the spatial domain. The local agent density function $M_i\left(t\right)$ denotes the tumor agent count at site $x_i\in {\widehat{U}}_{\mathrm{ABM}}$ and time t.

Let $\overrightarrow{M}\left(t\right)(M_1\left(t\right),\dots ,M_{N^2}\left(t\right))$ represent the discrete ABM configuration at time t. A continuous-time Markov process on ${\mathbb{Z}}_{+}^{N^2}$ governs the stochastic evolution of $\overrightarrow{M}\left(t\right)$, with state-dependent transition rates determined by agent rules, including proliferation, death, motility, and phenotypic switching.

Let $P(\overrightarrow{M},t)\mathbb{P}[\mathrm{ABM}\mathrm{state}=\overrightarrow{M}\mathrm{at}\mathrm{time}t]$ denote the probability mass function over all ABM configurations. The CME governing the evolution of $P(\overrightarrow{M},t)$ is:

$$\begin{array}{c}\hfill {\displaystyle \frac{d}{dt}}P(\overrightarrow{M},t)=\sum _{{\overrightarrow{M}}^{\prime }\ne \overrightarrow{M}}\left[T\left(\overrightarrow{M}\right|{\overrightarrow{M}}^{\prime })P({\overrightarrow{M}}^{\prime },t)-T\left({\overrightarrow{M}}^{\prime }\right|\overrightarrow{M})P(\overrightarrow{M},t)\right],\end{array}$$

where $T\left({\overrightarrow{M}}^{\prime }\right|\overrightarrow{M})$ denotes the transition rate from configuration $\overrightarrow{M}$ to ${\overrightarrow{M}}^{\prime }$. Each transition corresponds to elementary stochastic events with specific forms:

(i)

Brownian Movement: each cell jumps to neighboring site $x_j\in \mathcal{N}\left(i\right)$ at rate

$$\begin{array}{c}\hfill T_{ij}\left(M_i\right)={\displaystyle \frac{{\epsilon }^2}{2\mathrm{\Delta }{x}^2}}{M}_i,\end{array}$$

where $\mathcal{N}\left(i\right)$ denotes the Von Neumann neighborhood (four orthogonal neighbors) of $x_i$. The coefficient ${\displaystyle \frac{{\epsilon }^2}{2}}$ originates from the Fokker-Planck equation ${\partial }_t\rho ={\displaystyle \frac{{\epsilon }^2}{2}}\rho$ for tumor agents undergoing Brownian motion: $d{a}^{(X,Y)}\left(t\right)=\epsilon d{W}_t^a$ with ${\left\{W_t^a\right\}}_{a\in {\mathrm{\Lambda }}_t}$ independent standard Brownian motions.

(ii)

Regulated Proliferation: Using the Heaviside function H to identify normoxic states, only normoxic cells proliferate at rate

$$\begin{array}{c}\hfill {\lambda }_i\left(M_i\right)={\alpha }_{n}H(o_i-o_{\mathrm{hyp}}){\left(1-{\displaystyle \frac{M_i}{M_{max}(x_i,\mathrm{\Delta }x)}}\right)}^{+}{M}_i.\end{array}$$

(iii)

Drug-Induced Death: Apoptosis occurs at constant rate ${\kappa }_d$: ${\mu }_i\left(M_i\right)={\kappa }_d{d}_i{M}_i$.

(iv)

Mutation: mutation preserves cell counts per voxel, and thus does not appear in the CME.

Define the expected local agent density:

$$\begin{array}{c}\hfill {\overline{m}}_i\left(t\right)=\mathbb{E}\left(M_i\left(t\right)\right)=\sum _{\overrightarrow{M}}{M}_i\left(t\right)\mathbb{P}(\overrightarrow{M},t).\end{array}$$

The cádlág property $M_i\left(t\right)$ yields:

$$\begin{array}{c}\hfill {\displaystyle \frac{d}{dt}}{\overline{m}}_i\left(t\right)=\sum _{\overrightarrow{M}}{M}_i\left(t\right){\displaystyle \frac{d}{dt}}\mathbb{P}(\overrightarrow{M},t).\end{array}$$

The CME governs the change rate due to proliferation, death, and movement:

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\overline{m}}_i}{dt}}& =\mathbb{E}\left[{\lambda }_i\left(M_i\right)\right]-\mathbb{E}\left[{\mu }_i\left(M_i\right)\right]\hfill \\ \hfill & +\sum _{x_j\in \mathcal{N}\left(i\right)}\left({\displaystyle \frac{{\epsilon }^2}{2\mathrm{\Delta }{x}^2}}\mathbb{E}\left[M_j\right]-{\displaystyle \frac{{\epsilon }^2}{2\mathrm{\Delta }{x}^2}}\mathbb{E}\left[M_i\right]\right).\hfill \end{array}$$

This equation is unclosed because first-order moments depend on higher-order correlations through nonlinear functions. A tractable macroscopic approximation requires closure. The mean-field approximation assumes statistical independence between agent states, providing the first-moment closure:

$$\begin{array}{c}\hfill \mathbb{E}\left[f\left(M_i\right)\right]\approx f\left(\mathbb{E}\left[M_i\right]\right)=f\left({\overline{m}}_i\right).\end{array}$$

Applying this closure yields:

$$\begin{array}{cc}\hfill {\displaystyle \frac{d{\overline{m}}_i}{dt}}& ={\displaystyle \frac{{\epsilon }^2}{2\mathrm{\Delta }{x}^2}}\sum _{x_j\in \mathcal{N}\left(i\right)}({\overline{m}}_j-{\overline{m}}_i)\hfill \\ \hfill & +\left[{\alpha }_{n}H(o_i-o_{\mathrm{hyp}}){\left(1-{\displaystyle \frac{{\overline{m}}_i}{M_{max}(x_i,\mathrm{\Delta }x)}}\right)}^{+}-{\kappa }_d{d}_i\right]{\overline{m}}_i.\hfill \end{array}$$

Define the density $\rho (x_i,t){\displaystyle \frac{{\overline{m}}_i\left(t\right)}{\mathrm{\Delta }{x}^2}}$ gives ${\overline{m}}_i\left(t\right)=\rho (x_i,t)\mathrm{\Delta }{x}^2$, producing the CME for $\rho$:

$$\begin{array}{cc}\hfill {\displaystyle \frac{d\rho }{dt}}(x_i,t)& ={\displaystyle \frac{{\epsilon }^2}{2\mathrm{\Delta }{x}^2}}\sum _{x_j\in \mathcal{N}\left(i\right)}(\rho (x_j,t)-\rho (x_i,t))\hfill \\ \hfill & +\left[{\alpha }_{n}H(o_i-o_{\mathrm{hyp}}){\left(1-{\displaystyle \frac{\rho (x_i,t)}{M_{max}(x_i,\mathrm{\Delta }x)/\mathrm{\Delta }{x}^2}}\right)}^{+}-{\kappa }_d{d}_i\right]\rho (x_i,t)\hfill \end{array}$$

The deterministic approximation holds under the law of large numbers when local agent populations are large and suppress fluctuations. The continuum limit yields the macroscopic density equation:

$$\begin{array}{c}\hfill {\partial }_t\rho ={\displaystyle \frac{{\epsilon }^2}{2}}\mathrm{\Delta }\rho +\left[{\alpha }_{n}H(o-o_{\mathrm{hyp}}){\left(1-{\displaystyle \frac{\rho }{{\rho }_{max}}}\right)}^{+}-{\kappa }_{d}d\right]\rho ,\end{array}$$

where ${\rho }_{max}=\underset{\mathrm{\Delta }x\to 0}{lim}{\displaystyle \frac{M_{max}(x_i,\mathrm{\Delta }x)}{\mathrm{\Delta }{x}^2}}$. In the continuum limit, the local cell count is

$$\begin{array}{c}\hfill F(x,t)={\int }_{B_{R_c}\left(x\right)}[{\rho }_{\mathrm{tumor}}(y,t)+{\rho }_{\mathrm{vessel}}(y,t)]dy,\end{array}$$

where ${\rho }_{\mathrm{tumor}}$ and ${\rho }_{\mathrm{vessel}}$ denote tumor and vessel agent spatial densities. The uniform bound $F(x,t)\le F_{max}$ implies the tumor cell density constraint:

$$\begin{array}{c}\hfill {\rho }_{\mathrm{tumor}}^{max}={\displaystyle \frac{F_{max}}{|B_{R_c}|}}-{\displaystyle \frac{1}{|B_{R_c}|}}{\int }_{B_{R_c}\left(x\right)}{\rho }_{\mathrm{vessel}}(y,t)dy\le {\displaystyle \frac{F_{max}}{|B_{R_c}|}},\end{array}$$

where $|B_{R_c}|=\pi R_c^2$ begin the disk area. The maximum local cell count follows as

$$\begin{array}{c}\hfill M_{max}(x_i,\mathrm{\Delta }x)={\rho }_{\mathrm{tumor}}^{max}\mathrm{\Delta }{x}^2={\displaystyle \frac{F_{max}}{\pi R_c^2}}\mathrm{\Delta }{x}^2,\end{array}$$

thus ${\rho }_{max}=\underset{\mathrm{\Delta }x\to 0}{lim}{\displaystyle \frac{M_{max}(x_i,\mathrm{\Delta }x)}{\mathrm{\Delta }{x}^2}}={\displaystyle \frac{F_{max}}{\pi R_c^2}}$.

Therefore, applying first-order moment closure $\mathbb{E}\left[f\left(M_i\right)\right]\approx f\left(\mathbb{E}\left[M_i\right]\right)$ yields a deterministic approximation of the tumor density dynamics. This closure is valid under large-population and weak-correlation regimes.

Theorem 4

(Mean–Field Limit for the Tumor Cell Density). Let ${\left\{M_i\left(t\right)\right\}}_{i=1}^{N^2}$ denote tumor cell agent counts evolving under the stochastic rules in Section 2 on a lattice ${\widehat{U}}_{\mathrm{ABM}}$ with spacing $\mathrm{\Delta }x>0$. Assume:

(i)

Initial agent positions and states are i.i.d. with macroscopic density ${\rho }_0\in L^1\left(U\right)\cap L^{\infty }\left(U\right)$.

(ii)

The maximum local occupancy $M_{max}(x_i,\mathrm{\Delta }x)$ satisfies

$$\begin{array}{c}\hfill M_{max}(x_i,\mathrm{\Delta }x)={\displaystyle \frac{F_{max}}{\pi R_c^2}}\mathrm{\Delta }{x}^2,\end{array}$$

where $F_{max}$ is the maximal local cell–vessel occupancy value and $R_c$ the cell radius.

Then, as $\mathrm{\Delta }x\to 0$ and $N_t\to \infty$, the tumor cell density satisfies:

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}{\partial }_t\rho \hfill & ={\displaystyle \frac{{\epsilon }^2}{2}}\mathrm{\Delta }\rho +\left[{\alpha }_{n}H(o-o_{\mathrm{hyp}}){\left(1-{\displaystyle \frac{\rho }{{\rho }_{max}}}\right)}^{+}-{\kappa }_{d}d\right]\rho ,\hfill \\ \rho (0,x)\hfill & ={\rho }_0\left(x\right).\hfill \end{array}\right.\end{array}$$

(4.5)

Here ${\rho }_{max}=F_{max}/\left(\pi R_c^2\right)$ is the macroscopic carrying capacity, H the Heaviside function, ${(·)}^{+}$ the positive part, and d the drug concentration. Diffusion arises from cellular Brownian motion; the reaction term models normoxic proliferation and drug-induced apoptosis.

Theorem 4 provides the leading–order deterministic approximation via first–order moment closure. Agent-based simulations are computationally intensive for tumor populations. The mean-field limit unifies the hybrid model into a computationally efficient continuum description:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\partial }_t\rho ={\displaystyle \frac{{\epsilon }^2}{2}}\mathrm{\Delta }\rho +\left[{\alpha }_{n}H(o-o_{\mathrm{hyp}}){\left(1-{\displaystyle \frac{\rho }{{\rho }_{max}}}\right)}^{+}-{\kappa }_{d}d\right]\rho ,\hfill \\ {\partial }_{t}n=D_n\mathrm{\Delta }n-\nabla ·\left({\displaystyle \frac{{\chi }_0}{1+\alpha c}}n\nabla c\right),\hfill \\ {\partial }_{t}c=D_c\mathrm{\Delta }c-{\xi }_{c}c+\eta \rho H(o-o_{\mathrm{hyp}})-\lambda cn,\hfill \\ {\partial }_{t}d=D_d\mathrm{\Delta }d-{\xi }_{d}d-{\rho }_{d}d\rho +S_{d}n,\hfill \\ {\partial }_{t}o=D_o\mathrm{\Delta }o-{\xi }_{o}o-{\rho }_o\rho +S_o\left(1-o\right)n.\hfill \end{array}\right.\end{array}$$

(4.6)

The mean-field equation (4.6) allows one to interpret complex, stochastic agent-based behavior through simpler, deterministic PDEs. This bridge from microscale (single-cell events like stochastic phenotype switching) to macroscale (tumor cell density evolution) enables: (i) analytical estimation of tumor front velocity, aiding prognosis (e.g., time to vascular invasion or metastatic risk), (ii) sensitivity to environmental factors (e.g., hypoxia, drug gradients) to be encoded directly in continuous equations, which may inform spatially-resolved treatment planning, and (iii) approximate solutions in real time, supporting integration into clinical decision-support tools, where full ABM simulations may be computationally prohibitive. On the other hand, mean-field approximations neglect fluctuations, so we clarify valid parameter regimes and biological relevance.

Remark 2

(Validity Regime of the Moment Closure). Equation (4.5) uses first–order moment closure $\mathbb{E}\left[f\left(M_i\right)\right]\approx f\left(\mathbb{E}\left[M_i\right]\right)$, neglecting second-order moments $\mathrm{Var}\left(M_i\right)$ and higher–order cumulants. This requires:

(i)

Large local populations:$M_i\gg 1$ per site $x_i$, making demographic noise $O\left(M_i^{-1/2}\right)$ under central limit scalings;

(ii)

Weak correlations:Negligible spatial correlations between neighboring voxels at scale $\mathrm{\Delta }x$, ensured by fast mixing (large ε) or weakly interactions.

Fluctuations then contribute only $O\left(M_i^{-1}\right)$ corrections, so (4.5) captures population dynamics to leading order.

Biologically, these hold in dense tumor regions where local interactions average stochasticity. They fail in:

(i)

Early growth stages, invasion fronts, or metastatic niches with small $M_i$ and non-negligible extinction;

(ii)

Strongly heterogeneous microenvironments where correlations and clustering invalidate the weak–fluctuation assumption.

Solutions include second-order moment closure (e.g., pair approximations) or hybrid models retaining ABM structures in critical regions while using PDE surrogates elsewhere. These are active multiscale modeling research areas.

5. Numerical Analysis & Results

We rigorously evaluate the numerical properties of the hybrid PDE-ABM scheme. First, we detail the discretization framework for both continuous and agent-based components. We then analyze consistency, stability, and convergence, verifying critical biological features including nonnegativity and mass conservation. Finally, we validate theoretical convergence through grid refinement tests.

5.1. Hybrid Time-Stepping Framework

We discretize a square computational domain of area $2.5\times {10}^{-5}{\mathrm{m}}^2$ using a uniform Cartesian grid with $\mathrm{\Delta }x=\mathrm{\Delta }y=5\times {10}^{-5}$ m, yielding a $100\times 100$ mesh. We evolve TAF, drug, and oxygen fields via a semi-implicit alternating direction implicit (ADI) scheme, updating reaction terms explicitly. Endothelial cell density advances using an explicit forward Euler method (5.1), subject to the CFL condition from Theorem 6.

We select the endothelial cell PDE time step $\mathrm{\Delta }t$ to ensure stability:

$$\begin{array}{c}\hfill \mathrm{\Delta }t\le min\left({\displaystyle \frac{\mathrm{\Delta }{x}^2}{4D_n}},{\displaystyle \frac{\mathrm{\Delta }x}{{2\parallel \chi \left(c\right)\nabla c\parallel }_{\infty }}}\right)=0.0050,\end{array}$$

where $\parallel \chi \left(c\right)\nabla c\parallel \approx 0.9922$ denotes the chemotactic flux supremum norm for our parameters. Choosing $\mathrm{\Delta }t=0.005$ ($\mathrm{\Delta }t=2.88\times {10}^2$ s) with $\mathrm{\Delta }x=0.01$ provides a conservative stability margin.

We update the ABM component with a larger time step $\mathrm{\Delta }{t}^{\prime }=0.1$ ($\mathrm{\Delta }{t}^{\prime }=5.76\times {10}^3$s), which is synchronized with the PDE component at each $\mathrm{\Delta }{t}^{\prime }$ interval. The hybrid time-stepping balances computational efficiency with stability for stiff PDE-ABM couplings. All simulations enforce homogeneous Neumann boundary conditions.

Table 2 summarizes the numerical parameters and scheme configurations.

We detail the numerical scheme for endothelial cell chemotaxis. The discretization employs second-order central differences in space and forward Euler integration in time [22]. The chemotaxis-diffusion equation yields deterministic and probabilistic formulations. The deterministic update is:

$$\begin{array}{c}\hfill {\displaystyle \frac{n_{i,j}^{k+1}-n_{i,j}^k}{\mathrm{\Delta }t}}=\left({\displaystyle \frac{D_n}{\mathrm{\Delta }{x}^2}}{\delta }_x^2{n}_{i,j}^k+{\displaystyle \frac{D_n}{\mathrm{\Delta }{y}^2}}{\delta }_y^2{n}_{i,j}^k-{\displaystyle \frac{1}{\mathrm{\Delta }x}}{\delta }_x{F}_{i,j}-{\displaystyle \frac{1}{\mathrm{\Delta }y}}{\delta }_y{G}_{i,j})\right)\end{array}$$

(5.1)

where F and G encode chemotactic drift computed at cell edges via linear interpolation.

For ABM compatibility, we express the update probabilistically:

$$\begin{array}{c}\hfill n_{i,j}^{k+1}=n_{i,j}^k{P}_0+n_{i+1,j}^k{P}_1+n_{i-1,j}^k{P}_2+n_{i,j+1}^k{P}_3+n_{i,j-1}^k{P}_4,\end{array}$$

(5.2)

with transition probabilities:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{P}_0=\begin{array}{ccc}\hfill & 1-4{\displaystyle \frac{D_n\mathrm{\Delta }t}{\mathrm{\Delta }{x}^2}}\hfill & \hfill -{\displaystyle \frac{\mathrm{\Delta }t{\chi }_0}{4\mathrm{\Delta }{x}^2}}\left({\displaystyle \frac{1}{1+\alpha c_{i+1,j}}}+{\displaystyle \frac{1}{1+\alpha c_{i,j}}}\right)(c_{i+1,j}-c_{i,j})\\ \hfill & \hfill & \hfill +{\displaystyle \frac{\mathrm{\Delta }t{\chi }_0}{4\mathrm{\Delta }{x}^2}}\left({\displaystyle \frac{1}{1+\alpha c_{i-1,j}}}+{\displaystyle \frac{1}{1+\alpha c_{i,j}}}\right)(c_{i,j}-c_{i-1,j})\\ \hfill & \hfill & \hfill -{\displaystyle \frac{\mathrm{\Delta }t{\chi }_0}{4\mathrm{\Delta }{x}^2}}\left({\displaystyle \frac{1}{1+\alpha c_{i,j+1}}}+{\displaystyle \frac{1}{1+\alpha c_{i,j}}}\right)(c_{i,j+1}-c_{i,j})\\ \hfill & \hfill & \hfill +{\displaystyle \frac{\mathrm{\Delta }t{\chi }_0}{4\mathrm{\Delta }{x}^2}}\left({\displaystyle \frac{1}{1+\alpha c_{i,j-1}}}+{\displaystyle \frac{1}{1+\alpha c_{i,j}}}\right)(c_{i,j}-c_{i,j-1})\end{array}\hfill \\ P_1=\begin{array}{ccc}\hfill & {\displaystyle \frac{D_n\mathrm{\Delta }t}{\mathrm{\Delta }{x}^2}}\hfill & \hfill -{\displaystyle \frac{\mathrm{\Delta }t{\chi }_0}{4\mathrm{\Delta }{x}^2}}\left({\displaystyle \frac{1}{1+\alpha c_{i+1,j}}}+{\displaystyle \frac{1}{1+\alpha c_{i,j}}}\right)(c_{i+1,j}-c_{i,j})\end{array}\hfill \\ P_2=\begin{array}{ccc}\hfill & {\displaystyle \frac{D_n\mathrm{\Delta }t}{\mathrm{\Delta }{x}^2}}\hfill & \hfill +{\displaystyle \frac{\mathrm{\Delta }t{\chi }_0}{4\mathrm{\Delta }{x}^2}}\left({\displaystyle \frac{1}{1+\alpha c_{i-1,j}}}+{\displaystyle \frac{1}{1+\alpha c_{i,j}}}\right)(c_{i,j}-c_{i-1,j})\end{array}\hfill \\ P_3=\begin{array}{ccc}\hfill & {\displaystyle \frac{D_n\mathrm{\Delta }t}{\mathrm{\Delta }{x}^2}}\hfill & \hfill -{\displaystyle \frac{\mathrm{\Delta }t{\chi }_0}{4\mathrm{\Delta }{x}^2}}\left({\displaystyle \frac{1}{1+\alpha c_{i,j+1}}}+{\displaystyle \frac{1}{1+\alpha c_{i,j}}}\right)(c_{i,j+1}-c_{i,j})\end{array}\hfill \\ P_4=\begin{array}{ccc}\hfill & {\displaystyle \frac{D_n\mathrm{\Delta }t}{\mathrm{\Delta }{x}^2}}\hfill & \hfill +{\displaystyle \frac{\mathrm{\Delta }t{\chi }_0}{4\mathrm{\Delta }{x}^2}}\left({\displaystyle \frac{1}{1+\alpha c_{i,j-1}}}+{\displaystyle \frac{1}{1+\alpha c_{i,j}}}\right)(c_{i,j}-c_{i,j-1})\end{array}\hfill \end{array}\right.\end{array}$$

(5.3)

The implementation uses five intervals

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill R_0& =[0,P_0],\hfill \\ \hfill R_j& =\left(\sum _{i=0}^{j-1}{P}_i,\sum _{i=0}^j{P}_i\right]\mathrm{for}j=1,\dots ,4.\hfill \end{array}\end{array}$$

(5.4)

We normalize $P_i$ such that $P_0+\dots +P_4=1$ (Assumption A2). Drawing $r\sim \mathrm{Uniform}[0,1]$, cells move according to the interval $R_j$ containing r. Sampling from $\mathrm{Uniform}\left[0,\sum _{i=0}^4{P}_i\right]$ with directional selection via (5.4) is equivalent to the sampling from $\mathrm{Uniform}[0,1]$ post-normalization.

Assumption A2

(Unity Summation). After computing $P_i$ from (5.3), we normalize:

$$\begin{array}{c}\hfill P_0+P_1+P_2+P_3+P_4=1.\end{array}$$

We implement the ADI method for TAF, drug, and oxygen equations (2.1). For a general reaction diffusion system ${\partial }_{t}u=D\mathrm{\Delta }u+f(u,x,y,t)$, the ADI scheme is:

$$\begin{array}{cc}\hfill \left(1-{\displaystyle \frac{D\mathrm{\Delta }t}{2\mathrm{\Delta }{x}^2}}{\delta }_x^2\right)U^{k+1/2}& =\left(1+{\displaystyle \frac{D\mathrm{\Delta }t}{2\mathrm{\Delta }{y}^2}}{\delta }_y^2\right)U^k+{\displaystyle \frac{\mathrm{\Delta }t}{2}}f\left(U^k\right),\hfill \\ \hfill \left(1-{\displaystyle \frac{D\mathrm{\Delta }t}{2\mathrm{\Delta }{y}^2}}{\delta }_y^2\right)U^{k+1}& =\left(1+{\displaystyle \frac{D\mathrm{\Delta }t}{2\mathrm{\Delta }{x}^2}}{\delta }_x^2\right)U^{k+1/2}+{\displaystyle \frac{\mathrm{\Delta }t}{2}}f\left(U^k\right),\hfill \end{array}$$

where $U^{k+1/2}\approx u\left((k+1/2)\mathrm{\Delta }t\right)$ and ${\delta }_x^2$, ${\delta }_y^2$ denote second-order central differences.

5.2. Consistency Verification

Section 5.1 established the discretization framework. We now examine approximation fidelity through a consistency analysis, quantifying local truncation errors. The endothelial chemotaxis scheme (5.1) and the ADI scheme yields:

$$\begin{array}{c}\hfill {\tau }_{\mathrm{endo}}=\mathcal{O}(\mathrm{\Delta }t+h^2),{\tau }_{\mathrm{ADI}}=\mathcal{O}(\mathrm{\Delta }{t}^2+h^2),\mathrm{for}\mathrm{\Delta }x=\mathrm{\Delta }y=h.\end{array}$$

Appendix D details derivations. The endothelial scheme has $\mathcal{O}(\mathrm{\Delta }t+h^2)$ truncation error, indicating first-order temporal and second-order spatial accuracy. The ADI scheme achieves $\mathcal{O}(\mathrm{\Delta }{t}^2+h^2)$, reflecting second-order convergence. These distinct characteristics motivate adaptive time-stepping for endothelial components while permitting larger ADI steps.

For grid discretization $V={\left\{V_{i,j}\right\}}_{i,j=0}^N$, define the discrete p-norm ($1\le p<\infty$):

$$\begin{array}{c}\hfill {\parallel V\parallel }_p{\left(\sum _{i,j=0}^N{V}_{i,j}^p\mathrm{\Delta }{x}^2\right)}^{1/p}.\end{array}$$

We establish iterative ${\ell }_2$ norm control:

Theorem 5

(Discrete Energy Bound). Under one-sided homogeneous Neumann boundary approximation (5.15) on $U={[0,1]}^2$, the scheme satisfies:

$$\begin{array}{c}\hfill \mathrm{\Delta }{x}^2\parallel n^k{\parallel }_2\le {\parallel n^{k+1}\parallel }_2\le \sqrt{\parallel n^k{\parallel }_2^2+\mathcal{O}\left(\mathrm{\Delta }t\right)}\end{array}$$

Proof. 

For any $1\le p<\infty$,

$$\begin{array}{c}\hfill {\parallel V\parallel }_p\le {\parallel V\parallel }_{\infty }.\end{array}$$

(5.5)

Cauchy-Schwartz inequality and $(N+1)\mathrm{\Delta }x=1$ yield:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \parallel n^k{\parallel }_1=\sum _{i,j=0}^N{n}_{i,j}\mathrm{\Delta }{x}^2& \le {\left(\sum _{i,j=0}^N\mathrm{\Delta }{x}^2\right)}^{1/2}{\left(\sum _{i,j=0}^N{\left(n_{i,j}^k\right)}^2\mathrm{\Delta }{x}^2\right)}^{1/2}\hfill \\ \hfill & ={\left(\sum _{i,j=0}^N{\left(n_{i,j}^k\right)}^2\mathrm{\Delta }{x}^2\right)}^{1/2}={\parallel n^k\parallel }_2.\hfill \end{array}\end{array}$$

(5.6)

Equations (5.5) and (5.6), and mass conservation (5.16) imply:

$$\begin{array}{c}\hfill \mathrm{\Delta }{x}^2\parallel n^k{\parallel }_2\le \mathrm{\Delta }{x}^2\parallel n^k{\parallel }_{\infty }\le \parallel n^k{\parallel }_1=\parallel n^{k+1}{\parallel }_1\le {\parallel n^{k+1}\parallel }_2\end{array}$$

Given $P_i\ge 0$ and $\sum _{i=0}^4{P}_i=1$ (Assumption A2), Jensen’s inequality applied to Equation (5.2) gives:

$$\begin{array}{c}\hfill {\left(n_{i,j}^{k+1}\right)}^2\le {\left(n_{i,j}^k\right)}^2{P}_0+{\left(n_{i+1,j}^k\right)}^2{P}_1+{\left(n_{i-1,j}^k\right)}^2{P}_2+\left(n_{i,j+1}^2\right)P_3+\left(n_{i,j-1}^2\right)P_4.\end{array}$$

Thus:

$$\begin{array}{c}\hfill \parallel n^{k+1}{\parallel }_2^2\le \sum _{i,j=0}^N\left({\left(n_{i,j}^k\right)}^2{P}_0+{\left(n_{i+1,j}^k\right)}^2{P}_1+{\left(n_{i-1,j}^k\right)}^2{P}_2+{\left(n_{i,j+1}^k\right)}^2{P}_3+{\left(n_{i,j-1}^k\right)}^2{P}_4\right).\end{array}$$

The difference expands as:

$$\begin{array}{cc}\hfill \parallel n^k{\parallel }_2^2-& \sum _{i,j=0}^N\left({\left(n_{i,j}^k\right)}^2{P}_0+{\left(n_{i+1,j}^k\right)}^2{P}_1+{\left(n_{i-1,j}^k\right)}^2{P}_2+\left(n_{i,j+1}^2\right)P_3+\left(n_{i,j-1}^2\right)P_4\right)\hfill \\ \hfill =& \sum _{i,j=0}^N\left({\left(n_{i,j}^k\right)}^2-{\left(n_{i+1,j}^k\right)}^2\right)P_1+\sum _{i,j=0}^N\left({\left(n_{i,j}^k\right)}^2-{\left(n_{i-1,j}^k\right)}^2\right)P_2\hfill \\ \hfill +& \sum _{i,j=0}^N\left({\left(n_{i,j}^k\right)}^2-{\left(n_{i,j+1}^k\right)}^2\right)P_3+\sum _{i,j=0}^N\left({\left(n_{i,j}^k\right)}^2-{\left(n_{i,j-1}^k\right)}^2\right)P_4\hfill \\ \hfill \stackrel{\mathrm{telescoping}}{\underset{\mathrm{summation}}{=}}& \sum _{j=0}^N\left({\left(n_{0,j}^k\right)}^2-{\left(n_{N,j}^k\right)}^2\right)(P_1-P_2)+\sum _{i=0}^N\left({\left(n_{i,0}^k\right)}^2-{\left(n_{i,N}^k\right)}^2\right)(P_3-P_4),\hfill \end{array}$$

using one-sided Neumann approximations (5.15). The difference $P_1-P_2$ satisfies:

$$\begin{array}{c}\hfill P_1-P_2=-{\displaystyle \frac{\mathrm{\Delta }t}{4\mathrm{\Delta }{x}^2}}\left(\chi \left(c_{i+1,j}\right)(c_{i+1,j}-c_{i,j})+\chi \left(c_{i-1,j}\right)(c_{i,j}-c_{i-1,j})+\chi \left(c_{i,j}\right)(c_{i+1,j}-c_{i,j})\right).\end{array}$$

Assuming bounded first-order partial derivatives of c, we have $|P_1-P_2|=\mathcal{O}\left({\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}}\right)$. Similarly, $|P_3-P_4|=\mathcal{O}\left({\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{\Delta }x}}\right)$. The extra terms thus bound by $4\parallel n^k{\parallel }_{\infty }^2\mathcal{O}\left(\mathrm{\Delta }t\right)\le 4{\parallel n^0\parallel }_{\infty }^2\mathcal{O}\left(\mathrm{\Delta }t\right)=\mathcal{O}\left(\mathrm{\Delta }t\right)$.    □

We simulate truncation errors by excluding tumor and vessel agents (initialized to zero). This simplifies the TAF, drug, and oxygen PDEs to:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{\partial }_{t}c=D_c\mathrm{\Delta }c-{\xi }_{c}c,\hfill \\ {\partial }_{t}d=D_d\mathrm{\Delta }d-{\xi }_{d}d,\hfill \\ {\partial }_{t}o=D_o\mathrm{\Delta }o-{\xi }_{o}o.\hfill \end{array}\right.\end{array}$$

Homogeneous Neumann conditions apply on $U={[0,1]}^2$. We use the ghost point method to handle Neumann conditions, e.g., at $x=0$, we have $U_{-1,j}^k=U_{1,j}^k$ for $U\in \{c,d,o\}$.2 This treatment introduces $\mathcal{O}\left(\mathrm{\Delta }{x}^2\right)$ error, in alignment with $\mathcal{O}(\mathrm{\Delta }{t}^2+\mathrm{\Delta }t\mathrm{\Delta }{x}^2)$ truncation error of the ADI scheme. We use the following exact solutions as initial conditions:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}c(x,y,t)=e^{-(4{\pi }^2{D}_c+{\xi }_c)t}cos\left(2\pi x\right)cos\left(2\pi y\right),\hfill \\ d(x,y,t)=e^{-(16{\pi }^2{D}_d+{\xi }_d)t}cos\left(4\pi x\right)cos\left(4\pi y\right),\hfill \\ o(x,y,t)=4e^{-\left(4{\pi }^2{D}_{o}t{\xi }_o\right)t}cos\left(2\pi x\right)cos\left(2\pi y\right).\hfill \end{array}\right.\end{array}$$

We compute convergence order via:

$$\begin{array}{c}\hfill {\mathrm{od}}_u\left(h\right){\displaystyle \frac{log\left(\frac{{\mathrm{LTE}}_u\left(h\right)}{{\mathrm{LTE}}_u(h/2)}\right)}{log\frac{h}{h/2}}},\end{array}$$

where ${\mathrm{LTE}}_u\left(h\right)$ is the local truncation error for $u\in \{c,d,o\}$ at resolution h. Figure 2 demonstrates convergence. Table 3 lists orders for $h\in \{0.02,0.01,0.005,0.0025\}$ with $\mathrm{\Delta }t={10}^{-7}$. Observed orders approach 2, confirming ADI’s second-order spatial accuracy.

In summary, the endothelial discretization exhibits $\mathcal{O}(\mathrm{\Delta }t+h^2)$ truncation error (first-order time, second-order space). The ADI scheme $\mathcal{O}(\mathrm{\Delta }{t}^2+h^2)$ (second-order time and space). These results align with theory and substantiate the convergence results in subsequent sections.

5.3. Stability Enforcement

While consistency guarantees alignment between the discretization and governing equations, it does not ensure that bounded or biologically meaningful solutions over time. This subsection analyzes the numerical stability of the explicitly discretized endothelial cell scheme, which exhibits conditional stability. We invoke the Lax–Richtmyer equivalence theorem to confirm convergence.

Theorem 6

(Conditional Stability). The ADI scheme maintains unconditional stability for pure linear diffusion. The finite difference scheme for endothelial chemotaxis (Equation 5.1) achieves conditional stability under the time step constraint:

$$\begin{array}{c}\hfill \mathrm{\Delta }t\le min\left\{{\displaystyle \frac{\mathrm{\Delta }{x}^2}{4D_n}},{\displaystyle \frac{\mathrm{\Delta }x}{{2\parallel \chi \left(c\right)\nabla c\parallel }_{\infty }}}\right\}.\end{array}$$

(5.7)

Subject to the additional constraints:

$$\begin{array}{c}\hfill \mathrm{\Delta }x\le 2{\displaystyle \frac{D_n}{{\parallel \chi \left(c\right)\nabla c\parallel }_{\infty }}},\mathrm{\Delta }t\le {\displaystyle \frac{1}{{\displaystyle \frac{4D_n}{\mathrm{\Delta }{x}^2}}+2{\displaystyle \frac{{\parallel \chi \left(c\right)\nabla c\parallel }_{\infty }}{\mathrm{\Delta }x}}}},\end{array}$$

(5.8)

all motility probabilities $P_0,P_1,P_2,P_3,P_4$ remain non-negative.

Proof. 

Define the auxiliary expressions:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}{A}_{i+{\displaystyle \frac{1}{2}},j}{\displaystyle \frac{\mathrm{\Delta }t{\chi }_0}{4\mathrm{\Delta }{x}^2}}\left({\displaystyle \frac{1}{1+\alpha c_{i+1,j}}}+{\displaystyle \frac{1}{1+\alpha c_{i,j}}}\right)(c_{i+1,j}-c_{i,j})\hfill \\ A_{i-{\displaystyle \frac{1}{2}},j}{\displaystyle \frac{\mathrm{\Delta }t{\chi }_0}{4\mathrm{\Delta }{x}^2}}\left({\displaystyle \frac{1}{1+\alpha c_{i-1,j}}}+{\displaystyle \frac{1}{1+\alpha c_{i,j}}}\right)(c_{i,j}-c_{i-1,j})\hfill \\ B_{i,j+{\displaystyle \frac{1}{2}}}{\displaystyle \frac{\mathrm{\Delta }t{\chi }_0}{4\mathrm{\Delta }{x}^2}}\left({\displaystyle \frac{1}{1+\alpha c_{i,j+1}}}+{\displaystyle \frac{1}{1+\alpha c_{i,j}}}\right)(c_{i,j+1}-c_{i,j})\hfill \\ B_{i,j-{\displaystyle \frac{1}{2}}}{\displaystyle \frac{\mathrm{\Delta }t{\chi }_0}{4\mathrm{\Delta }{x}^2}}\left({\displaystyle \frac{1}{1+\alpha c_{i,j-1}}}+{\displaystyle \frac{1}{1+\alpha c_{i,j}}}\right)(c_{i,j}-c_{i,j-1})\hfill \end{array}\right.\end{array}$$

It follows that

$$\begin{array}{c}\hfill max\left\{\left|A_{i+{\displaystyle \frac{1}{2}},j}\right|,\left|A_{i-{\displaystyle \frac{1}{2}},j}\right|,\left|B_{i,j+{\displaystyle \frac{1}{2}}}\right|,\left|B_{i,j-{\displaystyle \frac{1}{2}}}\right|\right\}\le {\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }x}}{\parallel \chi \left(c\right)\nabla c\parallel }_{\infty }\end{array}$$

(5.9)

The requirement

$$\begin{array}{c}\hfill {\displaystyle \frac{D_n\mathrm{\Delta }t}{\mathrm{\Delta }{x}^2}}\ge {\displaystyle \frac{\mathrm{\Delta }t}{2\mathrm{\Delta }x}}{\parallel \chi \left(c\right)\nabla c\parallel }_{\infty }\Rightarrow D_n\ge {\displaystyle \frac{\mathrm{\Delta }x}{2}}{\parallel \chi \left(c\right)\nabla c\parallel }_{\infty },\end{array}$$

(5.10)

ensures $P_1,P_2,P_3,P_4\ge 0$. Express $P_0$ from Equation (5.3) as:

$$\begin{array}{c}\hfill P_0=1-4{\displaystyle \frac{D_n\mathrm{\Delta }t}{\mathrm{\Delta }{x}^2}}-A_{i+{\displaystyle \frac{1}{2}},j}+A_{i-{\displaystyle \frac{1}{2}},j}-B_{i,j+{\displaystyle \frac{1}{2}}}+B_{i,j-{\displaystyle \frac{1}{2}}}\end{array}$$

Equation (5.9) yields:

$$\begin{array}{c}\hfill P_0\ge 1-4{\displaystyle \frac{D_n\mathrm{\Delta }t}{\mathrm{\Delta }{x}^2}}-{\displaystyle \frac{2\mathrm{\Delta }t}{\mathrm{\Delta }x}}{\parallel \chi \left(c\right)\nabla c\parallel }_{\infty }.\end{array}$$

The $\mathrm{\Delta }t$ constraint becomes:

$$\begin{array}{c}\hfill \mathrm{\Delta }t\le {\displaystyle \frac{1}{{\displaystyle \frac{4D_n}{\mathrm{\Delta }{x}^2}}+2{\displaystyle \frac{{\parallel \chi \left(c\right)\nabla c\parallel }_{\infty }}{\mathrm{\Delta }x}}}}.\end{array}$$

(5.11)

Equation (5.10) bounds $\mathrm{\Delta }x$:

$$\begin{array}{c}\hfill \mathrm{\Delta }x\le 2{\displaystyle \frac{D_n}{{\parallel \chi \left(c\right)\nabla c\parallel }_{\infty }}}.\end{array}$$

Under (5.11), this necessarily implies:

$$\begin{array}{c}\hfill \mathrm{\Delta }t\le {\displaystyle \frac{1}{{\displaystyle \frac{4D_n}{\mathrm{\Delta }{x}^2}}+2{\displaystyle \frac{{\parallel \chi \left(c\right)\nabla c\parallel }_{\infty }}{\mathrm{\Delta }x}}}}\le {\displaystyle \frac{1}{4{\displaystyle \frac{D_n}{{\left(2{\displaystyle \frac{D_n}{{\parallel \chi \left(c\right)\nabla c\parallel }_{\infty }}}\right)}^2}}+2{\displaystyle \frac{{\parallel \chi \left(c\right)\nabla c\parallel }_{\infty }}{\left(2{\displaystyle \frac{D_n}{{\parallel \chi \left(c\right)\nabla c\parallel }_{\infty }}}\right)}}}}={\displaystyle \frac{D_n}{{2\parallel \chi \left(c\right)\nabla c\parallel }_{\infty }^2}}.\end{array}$$

This constitutes the chemotaxis CFL condition, where ${\displaystyle \frac{\mathrm{\Delta }{x}^2}{4D_n}}$ originates from diffusion-dominated constraints in two dimensions. With nonnegative $P_i$ and $\sum _{i=0}^4{P}_i=1$ (Assumption A2), update (5.2) forms a convex combination. Thus:

$$\begin{array}{c}\hfill \parallel n^{k+1}{\parallel }_{\infty }\le \parallel n^k{\parallel }_{\infty }\le {\parallel n^0\parallel }_{\infty }\end{array}$$

Boundedness of the maximum norm over all time k establishes stability under the CFL and grid conditions.    □

Figure 3 visualizes the stability criteria, showing the admissible $(\mathrm{\Delta }x,\mathrm{\Delta }t)$ region. This empirically validates Theorem 6’s CFL-type condition ensuring nonnegative $P_i$. The theorem further guarantees that the update forms a convex combination, preserving endothelial density nonnegativity. This maintains biological realism and prevents unphysical numerical artifacts like negative concentrations.

Consistency and conditional stability imply convergence via the Lax-Richtmyer equivalence theorem for linear initial value problems.

Theorem 7

(Convergence). Let $n(x,y,t)$ be the exact endothelial chemotaxis solution and $n_{i,j}^k$ its numerical approximation. Under Theorem 6’s time step constraints and sufficient regularity of n and c, the scheme converges with:

$$\begin{array}{c}\hfill \parallel n_{i,j}^k-n(x_i,y_j,t_k)\parallel =\mathcal{O}(\mathrm{\Delta }t+\mathrm{\Delta }{x}^2).\end{array}$$

5.4. Biological Constraints Preservation

In addition to consistency, stability, and convergence, biological realism requires that the numerical scheme preserves the nonnegativity intrinsic to the modeled system. We now analyze whether the discretization maintains nonnegativity for endothelial cell densities. The explicit update scheme takes the form of a convex combination of neighboring cell densities:

$$\begin{array}{c}\hfill n_{i,j}^{k+1}=\sum _{\ell =0}^4{P}_{\ell }{n}_{\mathrm{neighbor}(\ell )}^k,\end{array}$$

where $P_{\ell }\ge 0$ and ${\sum }_{\ell =0}^4{P}_{\ell }=1$ (Assumption A2), provided the CFL condition holds. This structure guarantees:

$$\begin{array}{c}\hfill n_{i,j}^{k+1}\ge 0\mathrm{whenever}{n}_{i,j}^k\ge 0,\forall i,j.\end{array}$$

Theorem 8

(Nonnegativity). The finite difference update preserves nonnegativity of $n_{i,j}^k$ at each time step, assuming nonnegative initial data and satisfaction of the CFL constraints (5.7) and (5.8).

The second biological constraint is the conservation law, such as mass conservation or population invariance. Numerical schemes must respect this constraint, particularly in the absence of sources, sinks, or boundary fluxes. We therefore verify whether the total endothelial cell count $\sum _{i,j}{n}_{i,j}^k$ remains constant under the proposed numerical scheme. Summing the probabilistic motility update (5.1) over all grid points $(i,j)$, with $\mathrm{\Delta }x=\mathrm{\Delta }y$, yield:

$$\begin{array}{c}\hfill \sum _{i,j}{n}_{i,j}^{k+1}=\sum _{i,j}{n}_{i,j}^k+\mathrm{\Delta }t\sum _{i,j}\left({\displaystyle \frac{D_n}{\mathrm{\Delta }{x}^2}}{\delta }^2{n}_{i,j}^k-{\displaystyle \frac{1}{\mathrm{\Delta }x}}(\delta F_{i,j}+\delta G_{i,j})\right)\end{array}$$

(5.12)

We expand the second right-hand side term as:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \sum _{i,j}{\delta }^2{n}_{i,j}^k& =\sum _{i,j=0}^N\left(n_{i-1,j}^k+n_{i+1,j}^k+n_{i,j-1}^k+n_{i,j+1}^k-4n_{i,j}^k\right)\hfill \\ \hfill & =\sum _{j=0}^N\left(n_{-1,j}-n_{0,j}+n_{N+1,j}-n_{N,j}\right)+\sum _{i=0}^N(n_{i,-1}-n_{i,0}+n_{i,N+1}-n_{i,N})\hfill \end{array}\end{array}$$

(5.13)

The telescoping summations simplify to:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \sum _{i,j=0}^N(\delta F_{i,j}+\delta G_{i,j})& =\sum _{j=0}^N(F_{N+{\displaystyle \frac{1}{2}},j}-F_{-{\displaystyle \frac{1}{2}},j})+\sum _{j=0}^N(G_{i,j+N+{\displaystyle \frac{1}{2}}}-G_{i,-{\displaystyle \frac{1}{2}}})\hfill \\ \hfill & =\chi {\left(c\right)}_{N+{\displaystyle \frac{1}{2}},j}{n}_{N+{\displaystyle \frac{1}{2}},j}{\displaystyle \frac{c_{N+1,j}-c_{N,j}}{\mathrm{\Delta }x}}-\chi {\left(c\right)}_{-{\displaystyle \frac{1}{2}},j}{n}_{-{\displaystyle \frac{1}{2}},j}{\displaystyle \frac{c_{0,j}-c_{-1,j}}{\mathrm{\Delta }x}}\hfill \\ \hfill & +\chi {\left(c\right)}_{i,N+{\displaystyle \frac{1}{2}}}{n}_{i,N+{\displaystyle \frac{1}{2}}}{\displaystyle \frac{c_{i,N+1}-c_{i,N}}{\mathrm{\Delta }x}}-\chi {\left(c\right)}_{i,-{\displaystyle \frac{1}{2}}}{n}_{i,-{\displaystyle \frac{1}{2}}}{\displaystyle \frac{c_{i,0}-c_{i,-1}}{\mathrm{\Delta }x}}\hfill \end{array}\end{array}$$

(5.14)

We impose homogeneous Neumann boundary conditions using one-sided approximations:

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial \varphi }{\partial n}}=0\mathrm{on}\partial U\end{array}$$

for $\varphi \in \{n,c,d,o\}$. At $x=0$, this yields:

$$\begin{array}{c}\hfill {\displaystyle \frac{U_{0,j}^k-U_{-1,j}^k}{\mathrm{\Delta }x}}=0\Rightarrow U_{-1,j}^k=U_{0,j}^k,\end{array}$$

(5.15)

Analogous conditions hold at other boundaries: $x=1$ implies $U_{N+1,j}^k=U_{N,j}^k$, $y=0$ implies $U_{i,-1}^k=U_{i,0}^k$, and $y=1$ implies $U_{i,N+1}^k=U_{i,N}^k$. Combining (5.12)–(5.14) with these boundary conditions ensures exact mass conservation for pure diffusion or chemotaxis processes:

$$\begin{array}{c}\hfill \sum _{i,j}{n}_{i,j}^{k+1}=\sum _{i,j}{n}_{i,j}^k.\end{array}$$

(5.16)

Theorem 9

(Discrete Mass Conservation for Cell Density Update). Let $n_{i,j}^k$ denote the discrete approximation of endothelial cell density at time step k and grid point $(i,j)$. Under the probability finite difference scheme (5.2) and one-sided Neumann approximation (5.15), the total mass is conserved:

$$\begin{array}{c}\hfill \sum _{i,j}{n}_{i,j}^{k+1}=\sum _{i,j}{n}_{i,j}^k.\end{array}$$

The Neumann boundary conditions reflect biologically plausible zero-flux constraints at tissue boundaries. This mass conservation property ensures that endothelial cell mass remains invariant up to machine precision in the absence of external sources or sinks, maintaining physical and biological realism.

Having established the analytical and qualitative properties of the numerical scheme, we conduct a spatial grid refinement study to empirically validate the consistency and stability results derived in Section 5.2-Section 5.4. This test simulates the system in the absence of tumor and vessel agents.

We evaluate grid spacings $\mathrm{\Delta }x\in \{0.02,0.01,0.005\}$ with fixed time steps $\mathrm{\Delta }t=0.005$ for PDE updates and $\mathrm{\Delta }{t}^{\prime }=0.1$ for ABM decisions (disabled here). The simulation runs until $T=30$, and convergence is measured using the step-change metric at each time step k:

$$\begin{array}{c}\hfill E_h^k\underset{i,j}{max}\left(|c_{i,j}^k\left(h\right)-c_{i,j}^{k-1}\left(h\right)|+|o_{i,j}^k\left(h\right)-o_{i,j}^{k-1}\left(h\right)|\right),\end{array}$$

where c and o represent TAF and oxygen concentrations, respectively. We choose

$$\begin{array}{c}\hfill c(x,y,0)=1-\sqrt{{(x-0.5)}^2+{(y-0.5)}^2},d(x,y,0)=0,o(x,y,0)=0.5.\end{array}$$

as initial conditions. No drug treatment or cellular dynamics are applied. Figure 4 shows that ${log}_{10}{E}_h^k$ decays consistently across grid levels. Beyond $t=5$, the decay becomes nearly linear on a logarithmic scale, with least-squares regression revealing closely matched slopes for all resolutions. This confirms the scheme’s temporal stability and spatial consistency. Table 4 reports numerical values at $t=30$ and estimated convergence slopes.

In summary, the hybrid numerical framework, combining explicit probabilistic updates for endothelial dynamics with ADI-based diffusion solvers, demonstrates mathematical consistency, conditional stability, and convergence. Rigorous verification of non-negativity, mass conservation, and discrete energy control ensures biologically faithful and numerically reliable simulations of tumor growth and vascular remodeling.

6. Discussion

This study introduced a hybrid PDE-ABM framework (HDC model) that rigorously couples stochastic cell-level dynamics with tissue-scale continuum processes. Our approach formalizes a hybrid methodology with precise mathematical guarantees for formulation, numerical discretization, and closed-form mean field representation. We analyze the well-posedness of the hybrid model and provide comprehensive numerical evidence for its consistency, stability, convergence, and conservation properties. The model integrates stochastic phenotypic switching, metabolite-regulated angiogenesis, and hypoxia-resistance feedback, thereby addressing key limitations in previous hybrid tumor models. Earlier studies often relied on heuristic simulations without mathematical and numerical analysis or treated angiogenesis and resistance in isolation [17,22]. Our framework extends the work by Wang et al. [22] by embedding an exact Gillespie-based mutation and phenotype switching algorithm with drug-modulated mutation rates, deriving a continuum mean-field limit that connects discrete rules to continuum equations, and establishing mathematical well-posedness and numerical robustness for the hybrid scheme. The mean-field limit integrates the hybrid model into a unified PDE system (4.6), which motivates hybrid numerical schemes that retain full ABM structure in critical regions and replace it with efficient PDE surrogates elsewhere. These methodological advances enable rigorous and robust multi-scale tumor modeling not achieved in previous descriptive or simulation-based studies.

Beyond its methodological contributions, the HDC model yields several biological and clinical insights. First, model-predicted resistance timelines can inform adaptive therapy protocols designed to suppress or delay resistance [50,51]. Second, simulations of oxygen gradients can guide hypoxia-targeted interventions and optimize the application of hypoxia-activated prodrugs, which may enhance the efficacy of subsequent chemotherapy or radiotherapy [52,53,54,55]. Third, our framework also motivates experimental validation and calibration, such as in vitro tracking of endothelial tip cell dynamics and resistance evolution. Designed for integration with experimental and clinical data, the model supports a translational modeling pipeline. For example, oxygen partial pressure measurement or PET imaging with 18F-FMISO can calibrate the simulated oxygen distribution $o(x,t)$, while microvessel density measurements from CD31-stained immunohistochemistry can validate the vascular network $V_t$ generated by ABM [52,56]. Time-resolved single-cell RNA-seq data under treatment provide drug-induced mutation rates $k_{SR}\left(d\right)$ by linking transcriptional states to phenotypic transitions. Longitudinal MRI during therapy offers radiomic profiles (e.g., tumor volume changes) for comparison with simulated resistance trajectories [57,58]. By supporting multi-scale calibration and validation, the model enables in silico exploration of therapy optimization strategies.

Several simplifying assumptions in the current framework merit consideration. First, simulations and analysis occur in a two-dimensional domain, which reduces the complexity of three-dimensional vascular geometry and cell-TME interactions compared to in vivo tumors. Extending the framework to three dimensions would capture the full complexity of tumor vasculature but requires greater computational resources. Second, the model treats TAF and oxygen dynamics with simplified production, diffusion, and uptake mechanisms. Actual angiogenesis involves multiple signaling pathways such as VEGF and Dll4/Notch [59]), and diverse cell phenotypes including stalk cells [60] and pericytes [61]), which are not explicitly represented. Including these features may enhance biological fidelity. Third, the immune system remains unmodeled, even though immune infiltration, suppression, and therapy-induced immunogenicity are recognized as critical regulators of tumor evolution and resistance. This omission limits the model’s applicability to immunotherapy scenarios. Finally, the model does not account for pharmacokinetic or pharmacodynamic effects, nor does it simulate treatment-induced vascular normalization or off-target toxicity. Incorporating these elements would improve clinical relevance. Despite these simplifications, the current framework offers a tractable and mathematically robust platform for multiscale tumor dynamics, serving as a foundation for future biological extensions and clinical applications.

Future work can extend this foundation in several directions. Extending to three-dimensional geometries will allow investigation of how vascular topology and spatial drug gradients influence resistance dynamics. Incorporating additional microenvironmental features, such as interstitial fluid flow and matrix density, can improve predictions of drug distribution and cell motility. Beyond single-agent therapy, simulating sequential and combination therapies, particularly integrating immunotherapy and anti-angiogenic strategies, could identify schedules that exploit transient vessel normalization to suppress resistance. Linking model predictions with patient-specific imaging, such as DCE-MRI, would support personalized therapy strategies. Finally, the mean-field derivation also enables reduced ODE/PDE systems that approximate complex spatial-stochastic dynamics and facilitate analytical control strategies.

In summary, the proposed HDC model establishes a rigorous and extensible framework for coupling discrete stochastic cell behavior with continuum tissue-scale dynamics. By uniting mathematical rigor, multiscale biological realism, and translational capacity, this work advances the field from descriptive hybrid models to predictive, clinically actionable in silico platforms for adaptive cancer therapy.