Reduced-order modeling of solute transport within physiologically realistic solid tumor microenvironment

TL;DR

This paper presents a new computational model to simulate how plasma moves through the complex environment of solid tumors, incorporating realistic biological features.

Contribution

The study introduces a calibrated reduced-order model that integrates electrohydrodynamic forces and multiphase simulations for tumor plasma transport.

Findings

EHD increases inlet plasma intensity at the fenestra compared to non-EHD models.

Plasma perfusion in tumors shows two-stage kinetics with an initial advection-dominated phase.

The RAD model accurately reproduces plasma propagation after calibration.

Abstract

Solid tumors are characterized by densely packed extracellular matrices and limited vascularization, creating significant resistance to both diffusive and convective transport. Tumor growth depends on complex flow–structure interactions across multiple scales, while vascular abnormalities and enhanced permeability elevate interstitial pressure in the tumor microenvironment.

In this study, we developed an integrated computational framework with a theoretical modeling framework that couples three phase, viscous-laminar, transient simulations of glycocalyx-patched tumor vessel-resolving plasma, red blood cells (RBCs), and white blood cells (WBCs) and tracking their volume fractions with a calibrated reverse advection-diffusion (RAD) model for intratumoral plasma transport. The reduced-order tumor microenvironment model incorporates electrohydrodynamic (EHD) force at the tumor vessel wall via glycocalyx patches on the luminal surface.

At the fenestra, EHD increases inlet plasma intensity relative to a non-EHD framework across all 15 numerical models (means: 0.576 non-EHD vs. 0.722 EHD; gain
25.34%
). Numerical simulations of plasma perfusion in both the tumor ECM domain and a microfluidic benchmark exhibit two-stage kinetics, with an initial advection-dominated regime. The RAD model reproduces this behavior and, after temporal calibration, matches the observed propagation.

By using fully resolved, EHD-inclusive multiphase CFD simulations to calibrate a reduced-order RAD model parameterized by measurable geometric features, we bridge the gap between classical Darcy–Starling perfusion models and fully resolved CFD. The resulting framework provides a tractable mechanism-grounded tool for quantifying plasma progression in dense solid tumors and for establishing the baseline transport capacity of the tumor extracellular matrix, independent of solute-specific biochemical properties.

Figures

Physiologically inspired schematic and simulation geometry of a tumor vessel coupled to the surrounding tumor ECM domain: (A) Endothelial capillary schematic showing the endothelial glycocalyx (red patches), red blood cells (RBCs), white blood cells (WBCs), and a plasma-only fenestra that connects the vessel to the tumor extracellular matrix (ECM). This schematic motivates the reduced-order model; in panel (C), red lines indicate glycocalyx-patched regions where electrohydrodynamic (EHD) forcing is applied, whereas black lines represent non-patched endothelial surfaces with no applied EHD forcing. (B) Representative histology of a tumor vessel and a scanned cross-section of a dense pancreatic tumor implanted in mice; this morphology informed the simulated vessel. (C) Two dimensional biomimetic domain of the tumor microenvironment, where the reduced-order tumor vessel is connected to the tumor stroma via a fenestra. For numerical simulations, the domain in (C) is partitioned into two segments. The modeling strategy restricts the entry of suspended particulates at the fenestra by selecting a fenestra diameter smaller than typical particle sizes. (D–F) The two segments of the test geometry with their corresponding meshes: (D) a portion of the tumor vessel showing the circular wall depression that coincides with the fenestra location; (E) the ECM region with the fenestra (tumor-transport inlet); and (F) a zoomed-in view of the fenestra region with finer mesh elements. Plasma pressure and velocity fields extracted at the fenestra center (panel D) from the vessel-scale simulations define the inlet boundary conditions for the ECM transport domain (panel E), providing a one-way coupling from the intravascular to the extravascular space. Additional details of the ECM geometry are provided in Section 2.1.2.

Reduced-order tumor vessel and ECM geometries used in the present numerical framework (A) Reduced-order tumor vessel with inlet (right) and outlet (left). The x marks the admissible interval for the fenestra center used in the parametric sweep; varying x changes the model index (Models 1–15 per Equation 1). (B) Zoomed view from the actual reduced-order CAD model showing alternating surface segments red line shows glycocalyx-patched regions where localized electrohydrodynamic (EHD) forcing is applied on the luminal surface of the tumor vessel; while blue line shows non-patched endothelial wall regions with no EHD forcing on the luminal surface of the tumor vessel (implementation details in Section 2.2.1). (C) Representative H&E histology image of a tumor slice illustrating irregular microvascular morphology; this image is shown for physiological context only and is not used to define the simulation geometry. (D) Planar (2D) surrogate of the tumor ECM domain used for transport simulations; populated with circular fiber bundles to achieve the prescribed packing fraction. The ECM domain is coupled to the vessel only through the fenestra inlet boundary conditions derived from vessel-scale simulations.

Microfluidic visualization of dye-driven perfusion through a 3D-printed tumor ECM (A) Experimental setup for the microfluidic test. Red–dyed deionized water was used to mimic plasma perfusion through a tumor ECM domain 3D-printed in clear photopolymer resin (C, green box). (B) Harvard Apparatus syringe pump delivering a constant 2 mL/min flow. (D,E) Time-lapse percolation within the ECM at 20 s and 120 s, respectively, showing propagation of the dyed front inside the ECM; the green frame indicates the zoomed field of view from (C) during flow. Scale bar 1 cm.

EHD flow yields enhanced, more persistent plasma penetration at the fenestra compared with non-EHD flow (A,B) Plasma (primary phase) contour inside the tumor vessel with EHD (A) and without EHD (B). (C,D) Zoomed fenestra ROI from panels A and B (same simulation Model 5). The dashed rectangle delineates the identical physical area used for pixel counting and threshold analysis. (E) Fenestra-ROI intensity distributions show a clear rightward shift under EHD with minimal overlap, indicating stronger local plasma penetration; mean intensities 0.650 (EHD) and 0.525 (non-EHD). (F) Threshold (survival) curves summarize robustness across cutoffs the EHD curve dominates across thresholds, and the area-between-curves (ABC) equals 0.125 (Equations 21-22), demonstrating more persistent high-intensity pixels only with EHD.

Structural sensitivity of EHD-induced plasma extravasation across vessel geometries. EHD flow consistently yields higher mean plasma penetration than non-EHD across all 15 models: (A) For models 1–15 (fenestra center translated from near inlet to near outlet), grouped bars compare Non–EHD and EHD mean calibrated intensity for each model; higher values indicate greater plasma penetration. EHD exceeds Non–EHD in every model. (B) Per–model means (dots) with grand mean and 95% bootstrap CI (square with whiskers) for each condition. Across models, mean Non–EHD =0.576 and mean EHD =0.722 (mean Δ=0.146 , 95% CI [0.124, 0.166] ; Wilcoxon signed–rank p=6.1×10−5 ).

Spatiotemporal spreading of plasma in the tumor ECM shows an advection-to-diffusion transition in wetted area growth: (A–D) Contours of plasma volume fraction at (A) 0.01 s, (B) 0.05 s, (C) 0.10 s, and (D) 0.20 s. White circles denote impermeable fibers; blue is unfilled interstitium. The 0–1 color bar maps local plasma fraction. The wetted region expands from a thin inlet plume to a broad front, illustrating a qualitative shift from advection-driven advance toward diffusion-controlled spreading. (E) Wetted area versus time on log–log axes. The red dashed line marks the transition time t*=1.665×10−1 sec , separating the advection-dominant region (left) from the diffusion dominant region (right).

CFD and experiments both show advection-dominated growth regimes in ECM wetting, with different transition times: (A–C) Snapshots show front advance at t=0.01, 0.10, 0.20,sec respectively. (D) log–log growth of area A(t) with an advection-dominated regime transitioning to diffusion at tCFD*=0.2277 sec . (E) Early diffusion-dominated growth (A∼t1/2) transitioning to advection (A∼t) at texp*=4.1 sec .

CFD-informed reverse advection-diffusion RAD model showing solute plume evolution through a reduced-order fiber-dense tumor ECM: (A) The simplified rectangular computational domain (10×5 μm2) used for the reduced-order RAD model, containing randomly distributed circular ECM fiber bundles occupying 30% of the rectangular domain area. All fibers have equal radius and act as impermeable structures, mimicking dense tumor stroma and restricting solute motion. (B) The high-resolution unstructured triangular finite-element mesh generated for the theoretical analysis. Whereas (C) provides a zoom view of the mesh in the fenestra region, highlighting refined discretization around fiber boundaries to resolve steep concentration gradients. (D–H) The predicted spatiotemporal evolution of plasma-borne solute under the advection-diffusion transport model at t=0.1 sec, 0.5 sec, 1.0 sec, 2.0 sec, and 3.0 sec , respectively, where C~ represent the normalized concentration field.

Comparison of CFD and RAD model for plasma transport transition: (A–D) show the CFD-predicted plasma propagation at t=0.1,1,2 , and 3 sec , respectively. Plasma first penetrates swiftly along preferential pathways between fibers, then progressively spreads laterally and fills the matrix, reflecting a transition from advection-dominated entry to diffusion-driven equilibration. (E) Shows the occupied interstitial area A(t) as a function of time for the CFD solution, demonstrating a steep early-time growth regime followed by a slowed diffusion-controlled regime as the domain fills. (F) Shows the corresponding theoretical RAD model prediction for A(t) , showing the same trend as the CFD solution.

Temporal evolution of plasma penetration inside tumor ECM domain at different flow times: Panels show snapshots of plasma perfusion inside the tumor domain at (A) t = 0.1666 s, (B) 0.1800 s, (C) 0.1900 s, and (D) 0.2000 s within the diffusion-dominated region. Red denotes the high-fraction bulk, green + yellow denote the graded rim, and blue denotes the unperfused region. Panels show a persistent graded rim leading a slowly advancing bulk with a gradual recession of blue—indicates gradient-limited (diffusive) penetration through a tortuous, fiber dense matrix, consistent with subunit area-growth scaling.

Scaling reconciliation of wetted-area growth between experiment and CFD: Wetted-area dynamics align between experiment and CFD after time and thickness scaling. Red dashed: experiment vs. raw CFD (before scaling). Blue solid: experiment vs. scaled CFD (after scaling). Post-scaling, the offset is nearly constant and the trends coincide, indicating matched advection-dominated growth.

Histology-to-geometry reconstruction pipeline for future specimen-specific vessel–ECM modeling. (A) H&E histology of an orthotopic human pancreatic tumor slice illustrating irregular microvasculature; this specimen informed the vessel geometry. (B) Imaging-based vessel extracted from the tumor slice and used to define the curved vessel–ECM interface. (C) Tumor ECM domain (3D view) populated with circular fiber bundles to achieve the prescribed packing fraction; the curved red vessel indicates the histology-informed vessel–ECM boundary used to locate the fenestra.

Tables

TABLE 1. Mesh-refinement study: domain-averaged pressure and velocity. Differences between the 1.00 μm mesh and the finer 0.76 μm mesh are negligible, indicating spatial convergence.

Element size ( μ m)	Total elements	Pavg (Pa)	Vavg (m/s)
5.0	8k	2193.67	3.1×10−5
3.0	19k	2180.66	4.3×10−5
1.00	77k	2180.26	3.43×10−5
0.76	80k	2180.48	3.44×10−5

TABLE 2. Color-occupancy ratios for Figures 6A–D. Ratios are computed from pixel counts with blue normalized to 1; white pore pixels are excluded. Increasing R:B and (G+Y):B indicate deeper penetration and a thickening front.

Panel	Time (sec)	R:B	(G+Y):B	R:(G+Y):B
(A)	0.01	0.0020	0.0006	0.0020:0.0006:1
(B)	0.05	0.0740	0.0096	0.0740:0.0096:1
(C)	0.10	0.4386	0.0396	0.4386:0.0396:1
(D)	0.20	6.8276	0.0930	6.8276:0.0930:1

TABLE 3. Segment-wise fits on (log10t,log10A) for log10A=m log10t+b .

Segment	t –range (sec)	m	b	R2
Advection region	[1.60×10−3,1.664×10−1]	2.21	5.23	0.99
Diffusion dominant region	[1.666×10−1,2.00×10−1]	0.41	3.75	0.98

TABLE 4. Color–composition ratios in the diffusion–dominant region in Figure 10.

Panel	Time (sec)	R:B	(G+Y):B	(G+Y):R
A	0.1666	4.63	0.161	0.0348
B	0.1800	5.63	0.147	0.0261
C	0.1900	6.188	0.137	0.0221
D	0.2000	6.75	0.130	0.0193

Full text

Introduction

1

Solute transport in solid tumors is fundamentally governed by the structural and mechanical properties of the surrounding tissue (Gullino, 1980; Sefidgar et al., 2014). Irregular vascular architecture, dense stromal packing, and elevated interstitial resistance disrupt flow pathways and reduce transport efficiency (Vaupel et al., 1990; Jain et al., 2014; Koumoutsakos et al., 2013; Stylianopoulos et al., 2018; Jiang et al., 2020). These conditions create spatial variability in pressure-driven and diffusive mechanisms, producing highly non-uniform intratumoral distributions. Despite the enhanced permeability and retention effect often reported for nanoparticle delivery (Maeda, 2012), perfusion in tumors can be much lower than in surrounding healthy tissues due to vascular constriction and abnormal leakiness into the extracellular space (Stylianopoulos and Jain, 2013).

Solid tumors, which constitute the majority of cases, present as noncystic masses categorized as benign or malignant (Gavhane et al., 2011; Rahman et al., 2023) and exhibit high resistance to diffusion (Jain, 1996; Jain, 1998) with limited blood and lymphatic drainage that impedes advection (Chim and Mikos, 2018). Accurate, mechanics-based quantification of spatiotemporal solute distributions is therefore essential for clinical decision-making (d’Esposito et al., 2018), as flow–structure interactions and vascular leakiness elevate interstitial pressure and vary across patient-specific tumor topography (Soltani and Chen, 2011; Zhang et al., 2013). Beyond hydrodynamics, electrohydrodynamic (EHD) introduces electrical body forces that can bias plasma transport. In vivo, the endothelial glycocalyx presents spatially patched and heterogeneous surface charge distributions that can generate localized EHD responses at the vessel wall (Weinbaum et al., 2007; Huilgo and Phan-Thien, 1997; Sumets et al., 2018; Teodoro et al., 2024). These forces act directly on the plasma phase and can therefore influence transport behavior in regions where plasma extravasation is initiated, such as near fenestral openings. In confined tumor microvascular geometries, localized electric body forces may become comparable to hydrodynamic stresses at the vessel wall, motivating the need for explicit vessel-scale modeling of EHD effects. Resolving these interactions is essential to capture the mechanistic “gatekeeping” behavior of the vessel wall, which lies outside the scope of purely pressure-driven descriptions of transvascular transport. Blood flow is multiphase and non-Newtonian: plasma with suspended RBCs and WBCs forms a heterogeneous mixture whose collective effects give rise to shear-dependent rheology and dispersion behavior (Merrill et al., 1963; Baskurt and Meiselman, 2003). Two phase hemodynamics models (Jung et al., 2006a) improve over single-phase assumptions yet remain incomplete; three-phase formulations (Akash et al., 2023) better capture particulate interactions and motivate coupling geometry and biology in a predictive framework. Prior multiphase–EHD studies (Rahmat, 2017) provide valuable numerical background but typically do not incorporate explicit three-phase blood or patched glycocalyx charge (Kang et al., 2018). In parallel, classical mean-transport models (Darcy’s law (Baish et al., 2011); Starling’s equation (Sugihara-Seki and Fu, 2005)) remain useful at the tissue scale but overlook key features of tumor hemodynamics, including non-Newtonian, particulate rheology, and geometry-specific microstructure. At the tissue scale, Darcy–Starling models provide an effective mean description of transport but treat the plasma flux across the vessel wall as a prescribed quantity. In contrast, the present study interrogates the vessel-wall-scale physics that govern the initiation of plasma entry into the interstitium, including localized electrohydrodynamic forcing and multiphase interactions near fenestral openings. Likewise, single-phase computational fluid dynamics (CFD), which treats blood as a homogeneous continuum, cannot resolve the interactive dynamics and spatial distribution of suspended particulate matter (Womersley, 1957; Attinger, 1964; Srivastava, 2007). Recent computational studies have explored tumor transport using porous-media, discrete–continuum, and hybrid multiscale formulations, providing valuable insight into pressure-driven perfusion and solute redistribution in complex vascularized tissues (Vavourakis et al., 2018; Sweeney et al., 2024; Wirthl et al., 2020). Complementary review articles have further emphasized the importance of multiscale coupling, structural heterogeneity, and transport bottlenecks in cancer modeling (Mahesh et al., 2022; Budhwani et al., 2022; Lorenzo et al., 2024; Hadjicharalambous et al., 2021). However, these approaches generally prescribe vessel-wall transport phenomenologically and do not explicitly resolve electrohydrodynamic effects arising from heterogeneous glycocalyx charge at fenestral interfaces. The present study addresses this gap by directly coupling vessel-resolved, EHD-inclusive multiphase simulations to downstream ECM transport models.

Building on these advances and remaining gaps, we develop a three-phase simulation framework that explicitly resolves blood mechanics and glycocalyx-mediated EHD inside tumor vessels to mechanistically quantify—and, under controlled conditions, enhance plasma ingress into the tumor extracellular matrix (ECM) in a biologically faithful manner. Building on these high-fidelity simulations, we introduce a CFD-informed reverse advection-diffusion (RAD) model: an advection–diffusion framework whose boundary data and closures are calibrated from the resolved vessel physics and whose coefficients are expressed in terms of measurable geometric descriptors such as stromal packing fraction (Sefidgar et al., 2014; Sefidgar et al., 2015; Soltani and Chen, 2011). Here, “reverse” denotes a transport regime in which the net progression of plasma-borne species within the ECM region is biased in the opposite direction to the expected forward perfusion direction predicted by pressure-driven Darcy flow, due to geometric obstruction and elevated interstitial resistance. The fundamental objective of this study is to quantify the mechanistic barriers governing plasma ingress from a microvessel into the tumor extracellular matrix (ECM). While plasma transport ultimately enables the delivery of nutrients and other carried species in vivo, the present work focuses specifically on the fluid-mechanical behavior of the carrier phase itself. We emphasize plasma because its perfusion kinetics represent a key transport-limiting step in dense, poorly vascularized tumor tissue. Within this framework, plasma volume fraction computed in CFD, dye-front advancement observed in microfluidic experiments, and solute fields obtained from reduced-order transport modeling are treated as complementary, scale-dependent descriptors of the same underlying physical process: the spatiotemporal wetting and saturation of the tumor interstitium. By unifying these metrics, the model is designed to interrogate how localized vessel-scale mechanisms, such as electrohydrodynamic (EHD) forcing, bias the initiation and depth of interstitial plasma transport. Our aim is to retain the mechanistic fidelity of vessel-resolved three-phase CFD—including localized EHD at patched glycocalyx segments–while bridging the gap between classical Darcy–Starling formulations and modern fluid dynamics by embedding these vessel-scale effects into a reduced-order transport model. This motivates the following questions:

Q1. As a marker of enhanced biological realism, how would the inclusion of electrohydrodynamic (EHD) effects change the plasma extravasation across the vessel wall at the fenestra, relative to a non-EHD baseline?

Q2. As a parametric tool, could we set up a simplified, reduced-order reverse advection-diffusion (RAD), informed by vessel-scale CFD-derived inlet conditions, to reproduce the spatiotemporal progression of plasma transport in the ECM without explicitly resolving multiphase or electrohydrodynamic effects?

To address these questions, we first perform three phase CFD of tumor-vessel flow with EHD localized to the vessel domain to quantify changes in plasma-perfusion metrics relative to a non-EHD baseline; we then develop and calibrate the RAD model to express penetration and advection diffusion transitions explicitly in terms of packing fraction and geometry, thereby bridging the gap between classical mean-transport descriptions and fully resolved CFD. Finally, we benchmark the model’s key predictions (front advance, wetted-area scaling, and transition time) against complementary microfluidic validation.

Extension of prior work and current scope

1.1

The present study builds directly upon our previous multiphase modeling framework for tumor microvessels and plasma transport in the extracellular matrix (ECM) (Akash et al., 2023). In that earlier work, blood flow within the tumor vessel was modeled as a three-phase system (plasma–RBC–WBC), and plasma transport in the ECM was described using a two-phase formulation, focusing on fenestra-mediated plasma leakage under purely hydrodynamic conditions. The current study retains this multiphase foundation but introduces several new physical, methodological, and validation components that substantially extend the scope and insight of the framework. Electrohydrodynamic (EHD) physics: The primary physical extension is the inclusion of electrohydrodynamic (EHD) forcing at the vessel wall through glycocalyx-mediated surface charge heterogeneity. This mechanism was not considered in the prior work and enables direct investigation of how localized electrokinetic forces bias plasma ingress at the fenestra. Incorporating EHD reveals a consistent enhancement in plasma perfusion intensity relative to non-EHD conditions, with a mean gain of 25.34% across all cases. Systematic spatial sensitivity analysis. Beyond introducing EHD, the present study expands the analysis from a limited set of geometries to an ensemble of $[eqn]$ vessel models, generated by translating the fenestra center axially along the vessel over an effective length of $[eqn]$ m. This systematic parametric sweep allows us to quantify how fenestra location modulates EHD-driven plasma transport and to assess the robustness of EHD effects across spatial variability, rather than relying on a single representative configuration. Experimental validation. The computational framework is further strengthened through direct benchmarking against complementary microfluidic experiments using 3D-printed biomimetic ECM structures. These experiments validate plasma front advancement and perfusion trends predicted by the simulations, providing experimental grounding that was not present in the prior computational-only study. Reverse advection-diffusion (RAD) theoretical framework integration: In addition to the high-fidelity multiphase simulations, the present study introduces a reduced-order reverse advection–diffusion (RAD) model to describe plasma progression within the ECM. The RAD model is not a multiphase formulation and does not solve the Navier–Stokes equations. Instead, it provides a theoretical, continuum description of plasma transport based on an advection–diffusion equation posed on a simplified ECM geometry. The RAD model serves as a computationally efficient surrogate that reproduces the spatiotemporal progression of plasma transport observed in CFD, while deliberately abstracting away detailed fluid mechanics in favor of mechanistic interpretability and parametric exploration.

Methods

2

Domain reconstruction

2.1

Development of reduced-order tumor vessel

2.1.1

To interrogate electrohydrodynamic (EHD)-driven changes in plasma perfusion, we employed a biomimetic tumor microvessel coupled to an extracellular matrix (ECM) domain (Figure 1C). The reduced-order vessel geometry follows the configuration in our prior work and its validated in silico capillary model (Akash et al., 2023; Akash et al., 2024; Akash et al., 2021): a 500 $[eqn]$ m-long two dimensional planar channel of height 150 $[eqn]$ m with a single fenestra, modeled as a rounded wall depression, which serves as the sole plasma-ECM exchange interface. The vessel wall contains a micron-scale rounded endothelial depression that defines the geometric access region for plasma extravasation. This depression (Figure 1D) is $[eqn]$ deep and $[eqn]$ in diameter, sized to permit plasma entry while excluding larger formed elements such as red blood cells (RBCs, 7– $[eqn]$ ) (Nithiarasu, 2022). It is designed to fully contain the sub-micron fenestral opening that mediates plasma transport into the extracellular matrix. The effective fenestra (transport aperture) has a characteristic streamwise length of approximately $[eqn]$ m (Figure 1E). The larger depression ensures that the entire fenestral region is geometrically resolved within the vessel wall and that particulate exclusion (e.g., RBCs) is enforced upstream of the transport aperture. A single fenestra is modeled as a representative plasma–ECM exchange site, under the assumption that fenestrae along tumor capillaries act as spatially localized and largely independent transport units at the scales considered. This representation enables isolation of local plasma entry mechanics without introducing additional geometric repetition, while preserving the governing transport physics at the vessel–ECM interface. Endothelial-scale context, including endothelial cells, glycocalyx, RBCs, and white blood cells (WBCs), is illustrated in Figure 1A, and the near-fenestra mesh is shown in Figure 1D. To represent heterogeneous endothelial surface properties that mediate EHD, we applied wall-attached glycocalyx patches (Figure 2B) with fixed streamwise extent of $[eqn]$ Dull and Hahn (2022); Weinbaum et al. (2021). Patches were distributed in a randomized, nonperiodic pattern (some adjacent, others separated by a single gap) and held fixed across all cases for reproducibility. Numerically, patched regions were implemented via a user-defined forcing term localized to wall-adjacent control volumes, adding an electric body force consistent with the prescribed patched surface charge (details in Section 2.2.1). In the schematic and layout, red wall segments denote patched, glycocalyx-bearing regions where EHD forcing is applied (Figures 1A,C), whereas black wall segments in Figure 1C indicate non-patched wall with no EHD forcing.

Physiologically inspired schematic and simulation geometry of a tumor vessel coupled to the surrounding tumor ECM domain: (A) Endothelial capillary schematic showing the endothelial glycocalyx (red patches), red blood cells (RBCs), white blood cells (WBCs), and a plasma-only fenestra that connects the vessel to the tumor extracellular matrix (ECM). This schematic motivates the reduced-order model; in panel (C), red lines indicate glycocalyx-patched regions where electrohydrodynamic (EHD) forcing is applied, whereas black lines represent non-patched endothelial surfaces with no applied EHD forcing. (B) Representative histology of a tumor vessel and a scanned cross-section of a dense pancreatic tumor implanted in mice; this morphology informed the simulated vessel. (C) Two dimensional biomimetic domain of the tumor microenvironment, where the reduced-order tumor vessel is connected to the tumor stroma via a fenestra. For numerical simulations, the domain in (C) is partitioned into two segments. The modeling strategy restricts the entry of suspended particulates at the fenestra by selecting a fenestra diameter smaller than typical particle sizes. (D–F) The two segments of the test geometry with their corresponding meshes: (D) a portion of the tumor vessel showing the circular wall depression that coincides with the fenestra location; (E) the ECM region with the fenestra (tumor-transport inlet); and (F) a zoomed-in view of the fenestra region with finer mesh elements. Plasma pressure and velocity fields extracted at the fenestra center (panel D) from the vessel-scale simulations define the inlet boundary conditions for the ECM transport domain (panel E), providing a one-way coupling from the intravascular to the extravascular space. Additional details of the ECM geometry are provided in Section 2.1.2.

Sensitivity of perfusion to the location of plasma entry was probed by generating $[eqn]$ geometries in which the fenestra center was shifted axially along the vessel while keeping fenestra shape and the patch pattern unchanged (Figure 2A). With $[eqn]$ clearance zones at inlet and outlet, admissible centers span $[eqn]$ . We define

$$[eqn]$$

Reduced-order tumor vessel and ECM geometries used in the present numerical framework (A) Reduced-order tumor vessel with inlet (right) and outlet (left). The x marks the admissible interval for the fenestra center used in the parametric sweep; varying x changes the model index (Models 1–15 per Equation 1). (B) Zoomed view from the actual reduced-order CAD model showing alternating surface segments red line shows glycocalyx-patched regions where localized electrohydrodynamic (EHD) forcing is applied on the luminal surface of the tumor vessel; while blue line shows non-patched endothelial wall regions with no EHD forcing on the luminal surface of the tumor vessel (implementation details in Section 2.2.1). (C) Representative H&E histology image of a tumor slice illustrating irregular microvascular morphology; this image is shown for physiological context only and is not used to define the simulation geometry. (D) Planar (2D) surrogate of the tumor ECM domain used for transport simulations; populated with circular fiber bundles to achieve the prescribed packing fraction. The ECM domain is coupled to the vessel only through the fenestra inlet boundary conditions derived from vessel-scale simulations.

which places one non-overlapping fenestra per configuration over the effective length $[eqn]$ . The first admissible center at $[eqn]$ is Model 1, and the sequence proceeds in steps of $[eqn]$ so that Model 15 corresponds to $[eqn]$ .

All geometries were constructed and meshed in ANSYS Workbench 2023 R1 using DesignModeler (DM) for CAD and the integrated mesher (ANSYS Inc., Canonsburg, Pennsylvania). Because the computational domain is two-dimensional (Figures 1C–F), unstructured triangular meshes were employed, with local refinement around the fenestra and patched-wall regions to resolve steep pressure/velocity gradients while keeping the far field economical.

Development of the tumor extracellular matrix (ECM) domain

2.1.2

We constructed the adjacent ECM region and populated it to achieve a prescribed in-plane packing fraction. The ECM geometry was designed to match a representative solid packing fraction of 0.30, supported by quantitative color-decomposition analysis of tumor cross-sectional histology in our prior work (local range 0.23–0.41) (Akash et al., 2023) and by independent MRI measurements reporting absolute tumor interstitial fluid volume fractions of approximately 27.5%, consistent with a comparably dense solid ECM environment (Kim et al., 2004). In the ECM model, this packing fraction was implemented by representing stromal crowding by $[eqn]$ (Figure 2D) circular bundles placed within the ECM such that

$$[eqn]$$

with a target mean bundle diameter $[eqn]$ (range 3.0– $[eqn]$ ). Here $[eqn]$ is the total fiber cross-sectional area and $[eqn]$ is the ECM domain area (Equation 2). Fixing $[eqn]$ sets the porosity $[eqn]$ , which enters the transport closures for the effective diffusivity $[eqn]$ (via tortuosity) and the hydraulic permeability $[eqn]$ , thereby limiting both diffusive penetration and pressure-driven perfusion. The imposed fiber field was checked against orientation statistics from the same segmentation workflow and against reported stromal-alignment trends in tumor tissues (Eekhoff and Lake, 2020). Coupling to the vessel is applied along the extracted boundary locus at the fenestra by prescribing the plasma-field inlet (and associated flux) from the vessel simulation; this enforces a smooth interface without explicitly meshing the vessel solid within the ECM grid. The two dimensional ECM domain (Figure 2D) was meshed with approximately $[eqn]$ unstructured triangular elements (see also Figure 1E), with local refinement near the vessel–ECM interface. Near the boundary locus corresponding to the fenestra, a characteristic length scale $[eqn]$ was used in the near-wall mass-transfer closure and scaled locally with the extracted wall geometry. In the simulations, the multiphase transport calculation in the vessel yields the plasma volume fraction at the artificial circular depression (fenestra), which is then used as the inlet condition for the subsequent ECM computation (Figure 1E), where intratumoral plasma transport is tracked.

Simulation setup

2.2

Simulation method for blood flow inside the tumor vessel

2.2.1

We model the in silico transport of blood within tumor vessels as an unsteady, viscous–laminar, two dimensional (planar) flow. Numerical solutions are obtained in ANSYS Inc., 2023 R1 using its Eulerian multiphase model (no custom inter-particle equations). Plasma is the primary (continuous) phase; red blood cells (RBCs) and white blood cells (WBCs) are treated as two dispersed phases. The formulation follows the classical multiphase CFD framework (Gidaspow, 1994; Anderson and Jackson, 1967) and is appropriate for in vivo hematocrits of 30%–55% (Jung et al., 2006a). Closures (mixture shear-thinning, interphase drag, agglomeration shape factor) follow (Akash et al., 2023). The governing conservation equations for the tumor-vessel simulations, together with the associated electrohydrodynamic and rheological closure relations, are given in Equations 3–11 and are solved numerically.

Let $[eqn]$ index plasma $[eqn]$ , RBCs $[eqn]$ , and WBCs $[eqn]$ , with volume fractions $[eqn]$ and velocities $[eqn]$ in 2D, and let $[eqn]$ . Interpenetrating phases satisfy

$$[eqn]$$

and, under phasewise incompressibility,

$$[eqn]$$

Phase momenta evolve as

$$[eqn]$$

with a common mixture pressure $[eqn]$ , viscous stress $[eqn]$ , and interphase exchange $[eqn]$ satisfying $[eqn]$ . Phase densities used here are $[eqn]$ , $[eqn]$ , $[eqn]$ (Nithiarasu, 2022; Schwartz and Conley, 2020). The indicator $[eqn]$ activates electrohydrodynamic (EHD) forcing only in the tumor–vessel subdomain.

The electroquasistatic potential satisfies

$$[eqn]$$

where $[eqn]$ is the permittivity $[eqn]$ , $[eqn]$ is the vacuum permittivity, $[eqn]$ is the case-specific relative permittivity, and $[eqn]$ is the free charge density $[eqn]$ . The Kelvin (electrohydrodynamic) body force $[eqn]$ added to Equation 5 is

$$[eqn]$$

as in continuum EHD theory (Landau et al., 2013; Saville, 1997; Castellanos et al., 2003). In this study $[eqn]$ is treated as spatially uniform unless otherwise noted, so $[eqn]$ . The Fluent UDF applies the phase-aggregated force components (where $[eqn]$ is the electric potential)

$$[eqn]$$

where $[eqn]$ are phase charge densities and $[eqn]$ ; $[eqn]$ is the zeta potential applied on EHD-active glycocalyx patches. Electrical properties used in Equations 6–8 are: $[eqn]$ (Squires and Bazant, 2006); $[eqn]$ ; phase conductivities $[eqn]$ (Pauly and Schwan, 1966), $[eqn]$ (Zhbanov and Yang, 2015), $[eqn]$ (Yang et al., 1999); and a case-specific $[eqn]$ stated with each EHD case.

The non-Newtonian shear-thinning behavior of blood is introduced in this numerical framework as a phenomenological mixture-level closure calibrated to experimental rheology, rather than as an emergent consequence of resolved red blood cell deformation or aggregation. While RBCs and WBCs are modeled as dispersed Eulerian phases to capture phase segregation and transport heterogeneity, the present formulation does not resolve cell-scale mechanics required for self-consistent rheological emergence. Specifically, because the Eulerian dispersed phases are treated as interpenetrating continua without explicit cell-membrane deformation or aggregation dynamics, phase kinematics are decoupled from the constitutive mixture-level stress response, thereby necessitating the use of an empirical rheological closure. The blood’s shear-thinning behavior is represented by (Jung and Hassanein, 2008; Jung et al., 2006b). With hematocrit $[eqn]$ (RBC volume fraction) and scalar shear rate $[eqn]$ ,

$$[eqn]$$

where $[eqn]$ is the mixture dynamic viscosity (Pa $[eqn]$ s), $[eqn]$ (Nithiarasu, 2022), and we take $[eqn]$ (Schwartz and Conley, 2020). The effective RBC contribution is captured through the hematocrit-dependent mixture law; no separate constant $[eqn]$ is prescribed. The coefficients $[eqn]$ and $[eqn]$ are dimensionless Carreau–Yasuda–type fits to $[eqn]$ :

$$[eqn]$$

and, for $[eqn]$ ,

$$[eqn]$$

where $[eqn]$ is the (dimensionless) Yasuda exponent. Interphase drag follows Schiller–Naumann (Naumann and Schiller, 1935); RBC agglomeration enters via the dynamic shape factor

$$[eqn]$$

while WBCs use $[eqn]$ . Characteristic cell sizes used in closures are $[eqn]$ and $[eqn]$ (Nithiarasu, 2022; Schwartz and Conley, 2020). The transport equations inside tumor vessel were discretized using the finite-volume method implemented in ANSYS FLUENT. An implicit Eulerian multiphase formulation was employed, with pressure–velocity coupling handled via the phase-coupled SIMPLE algorithm. Spatial discretization used least-squares cell-based gradient reconstruction, second-order schemes for pressure, and second-order upwind discretization for convective terms. Temporal integration was performed using a second-order implicit transient formulation with a constant time step $[eqn]$ s. Across the ensemble of 15 tumor-vessel simulations, the Courant–Friedrichs–Lewy (CFL) number remained below unity, ensuring numerical stability and temporal convergence throughout the simulation duration.

Boundary conditions for the tumor-vessel calculation are: pressure inlet $[eqn]$ and pressure outlet $[eqn]$ (Zhao et al., 2007). At the outlet, the backflow volume fraction is set to RBC $[eqn]$ and WBC $[eqn]$ to prevent unphysical WBC re-entry while remaining consistent with the dispersed-phase specification. Electric potentials $[eqn]$ are applied only on tumor–vessel boundaries where EHD is active; elsewhere $[eqn]$ . Time integration employed a constant time step $[eqn]$ ; each vessel case was advanced for 2000 iterations, matching the ECM simulation, with convergence monitored via mass and momentum residuals and outlet flux balances. Computations were run in parallel on four Xeon cores at 3.1 GHz. From each of the fifteen vessel geometries, plasma pressure and velocity were sampled at the fenestra-center location (wall depression); the arithmetic averages of these fields define the inlet data subsequently used in the ECM (2D) plasma-perfusion simulation.

Mesh refinement study and numerical stability

2.2.1.1

A systematic mesh-refinement study was performed (see Table 1) for the tumor-vessel domain to verify spatial convergence of the multiphase CFD simulations and to ensure that the reported plasma perfusion metrics are independent of discretization resolution. The refinement analysis focused on domain-averaged pressure and velocity, which directly influence plasma transport and extravasation near the fenestral interface. Element size was progressively reduced while maintaining identical boundary conditions and solver settings. Between the $[eqn]$ and $[eqn]$ meshes, the change in domain-averaged pressure is $[eqn]$ (approximately $[eqn]$ relative to $[eqn]$ ), and the change in average velocity magnitude is $[eqn]$ (approximately $[eqn]$ relative to $[eqn]$ ). These sub-percent variations indicate spatial convergence of the numerical solution. Given the marginal accuracy gain relative to the increased computational cost associated with further refinement, the $[eqn]$ mesh (77k elements) was selected for all 15 vessel geometries to ensure consistent resolution across the parametric ensemble. Temporal discretization employed a fully implicit second-order scheme with a constant time step of $[eqn]$ sec. The resulting Courant–Friedrichs–Lewy (CFL) number remained below unity across all vessel realizations, ensuring numerical stability and temporal convergence.

TABLE 1: Mesh-refinement study: domain-averaged pressure and velocity. Differences between the 1.00 μm mesh and the finer 0.76 μm mesh are negligible, indicating spatial convergence.

Simulation method for plasma perfusion inside the ECM tumor region

2.2.2

Transport within the tumor extracellular matrix (ECM) domain was modeled as an unsteady, viscous–laminar, two-dimensional two-phase flow using the Eulerian multiphase formulation in ANSYS Inc., 2023 R1. A shared pressure field is solved while phasewise continuity and momentum equations advance the primary (plasma) and secondary (air) phases. Interphase coupling is represented through drag exchange terms with coefficients evaluated from the local Reynolds number, consistent with the Eulerian framework (Afolabi and Lee, 2014). Electrohydrodynamic forcing is not applied in the ECM calculation (EHD is restricted to the vessel domain), so the momentum balances here contain only hydrodynamic and interphase exchange terms. The ECM simulations consider plasma as the invading liquid phase and air as a displaced secondary phase, enabling explicit tracking of the advancing plasma front. For each phase $[eqn]$ (plasma $[eqn]$ and air $[eqn]$ ), volume fractions $[eqn]$ satisfy the constraint

$$[eqn]$$

and evolve according to the phase continuity equations

$$[eqn]$$

The phase momentum equations are solved in mixture form,

$$[eqn]$$

where $[eqn]$ is the mixture velocity, $[eqn]$ is the volume-fraction-weighted mixture density, and ** τ ** denotes the viscous stress tensor (Equations 12-14). Interphase momentum exchange is captured implicitly through the Eulerian multiphase closure implemented in FLUENT (ANSYS Inc., 2023; Ishii and Hibiki, 2010).

The ECM transport equations were discretized using the finite-volume method implemented in ANSYS FLUENT. An implicit Eulerian multiphase formulation was employed, with pressure–velocity coupling handled via the phase-coupled SIMPLE algorithm. Spatial discretization used least-squares cell-based gradient reconstruction, second-order schemes for pressure, and second-order upwind discretization for convective terms. Temporal integration was performed using a second-order implicit transient formulation with a constant time step $[eqn]$ sec. Throughout the ECM simulations, the Courant–Friedrichs–Lewy (CFL) number remained below unity, indicating stable time stepping for the advective transport. Each ECM simulation was advanced for 2000 time steps, matching the duration of the corresponding tumor-vessel simulations. Convergence at each time step was assessed by monitoring residual reduction for mass and momentum equations, as well as global flux balances at the ECM outlets. All computations were executed in parallel on four Xeon cores at 3.1 GHz with a typical wall-clock time of 3–4 h.

Geometric coupling between tumor vessel and ECM domains: the fenestra geometry is explicitly resolved only within the tumor vessel domain, where the curved wall depression serves to localize plasma extrusion while restricting the entry of suspended cellular phases based on size. In contrast, the ECM domain models plasma-only transport and therefore does not require explicit replication of the fenestral curvature. Accordingly, the vessel-to-ECM interface is treated as a transport inlet as a flat sampling plane, rather than a direct geometric continuation of the vessel wall. This planar inlet representation enables consistent mapping of time-resolved pressure and velocity fields from the vessel simulations into the ECM domain while avoiding curvature-induced numerical bias in the downstream transport calculations. This domain partitioning strategy follows our prior validated modeling approach and preserves the relevant transport physics at the vessel–ECM interface without introducing additional geometric constraints.

Boundary conditions of numerical simulation for plasma perfusion inside the ECM tumor region: The ECM inlet condition is constructed from the tumor vessel simulations in a time-resolved manner. For each of the $[eqn]$ vessel cases, the plasma pressure $[eqn]$ and velocity $[eqn]$ were sampled at the fenestra-center location at every time step $[eqn]$ , with $[eqn]$ t = $[eqn]$ . At each $[eqn]$ , the ensemble-averaged inlet fields supplied to the ECM domain are defined as

$$[eqn]$$

thereby preserving the transient vessel-scale dynamics in the downstream ECM transport simulations. To create a clean plasma-advance front in the interstitial spaces, the secondary phase was taken as air. The ECM domain was initialized as air-filled (air volume fraction $[eqn]$ , plasma volume fraction $[eqn]$ ), and plasma displaced air over time under the inlet fields (Equation 15). The air-filled initialization of the ECM domain is adopted as a controlled numerical reference state to enable clear identification of plasma entry and displacement within the ECM region. This formulation facilitates robust front tracking and consistent quantification of plasma penetration dynamics, without aiming to reproduce the fully saturated interstitial fluid environment in vivo. The tumor slice geometry has three outlets; each outlet was set to a fixed pressure of $[eqn]$ , following prior studies (Wu et al., 2013) and consistent with elevated interstitial pressures reported in tumors (Sven and Josipa, 2007; Stylianopoulos et al., 2018; Heldin et al., 2004). No-slip conditions are imposed on solid ECM boundaries, and surface tension effects are neglected due to the low capillary number regime considered. Material properties matched those used elsewhere in this study: plasma density and viscosity $[eqn]$ and $[eqn]$ (Nithiarasu, 2022); air density and viscosity $[eqn]$ and $[eqn]$ (Porterfield and Kruse, 1995). These settings, together with the inlet fields $[eqn]$ and outlet pressures, define a closed initial–boundary-value problem for plasma permeation through the ECM microstructure, from which we compute the space–time evolution of plasma volume fraction, interstitial velocities, and cumulative perfusion metrics.

Experimental setup

2.3

To replicate the physical features of a tumorous microenvironment, we fabricated a biomimetic pillar array using a Form $[eqn]$ stereolithography (SLA) resin printer (Formlabs Inc., United States). The structure (Figure 3C) was printed (3D version of Figure 2D) in clear photopolymer resin, achieving a spatial resolution of $[eqn]$ m. The design matched the packing fraction of 0.30 observed within tumors. After printing, the structure was rinsed in isopropyl alcohol and post-cured under UV light. The flow through the pillar network was driven using red-dyed deionized water as the working fluid. The fluid was supplied at a constant volumetric flow rate using a Harvard Apparatus syringe pump (Figures 3A,B), enabling precise control of perfusion conditions within the range of $[eqn]$ (Figures 3D,E are for $[eqn]$ ). This range ensured laminar flow within the micro-pillared structure and preserved quasi-two dimensional dynamics consistent with the assumptions of the numerical models. The transparent resin structure was illuminated from below using a diffuse white light source, while the evolving liquid front was imaged from above with a digital camera. The recorded image sequences were post-processed using image analysis software to extract the time evolution of the wetted area. Binarized frames were used to calculate the area of liquid penetration as a function of time, from which the spreading dynamics were quantified.

Microfluidic visualization of dye-driven perfusion through a 3D-printed tumor ECM (A) Experimental setup for the microfluidic test. Red–dyed deionized water was used to mimic plasma perfusion through a tumor ECM domain 3D-printed in clear photopolymer resin (C, green box). (B) Harvard Apparatus syringe pump delivering a constant 2 mL/min flow. (D,E) Time-lapse percolation within the ECM at 20 s and 120 s, respectively, showing propagation of the dyed front inside the ECM; the green frame indicates the zoomed field of view from (C) during flow. Scale bar 1 cm.

CFD-informed reverse advection-diffusion (RAD) model

2.4

In classical porous-media descriptions, such as Darcy–Starling formulations, advection is assumed to proceed in the direction of decreasing pressure from the vessel into the tissue, while diffusion proceeds down concentration gradients. However, in dense tumor extracellular matrices (ECM), elevated interstitial fluid pressure (IFP), strong geometric confinement by fiber bundles, and limited pore connectivity substantially suppress forward pressure-driven convection. As a consequence, the effective transport of plasma-borne species can oppose the direction predicted by pressure-driven models, despite the constitutive flux laws remaining unchanged (Baish et al., 2011; Sugihara-Seki and Fu, 2005). In this work, we introduce a reverse advection–diffusion (RAD) model to describe this regime. Importantly, the term “reverse” does not imply a reversal of diffusion or advection physics, nor negative diffusivity or flow inversion within the ECM. Instead, it denotes a transport regime in which the net progression of plasma-borne species within the ECM region is biased in the opposite direction to the expected forward perfusion direction predicted by pressure-driven Darcy flow, due to geometric obstruction and elevated interstitial resistance. The RAD model is therefore a conceptual and mathematical framework for capturing transport opposition relative to classical perfusion assumptions, rather than a direct solver for interstitial biofluid flow. The RAD model is formulated as a reduced-order theoretical transport surrogate, distinct from the multiphase ECM simulations described in Section 2.2.2. No momentum or pressure equations are solved in the ECM domain. Instead, solute transport is modeled directly through an advection–diffusion equation with prescribed advection strength, allowing isolation of transport kinetics without resolving multiphase or interstitial flow dynamics. In this sense, the RAD model is not a multiphase ECM model and is not dynamically coupled to the plasma–air displacement simulations of Section 2.2.2. Rather, it provides a simplified representation of plasma progression once entry into the ECM has occurred. The governing transport equation for solute concentration $[eqn]$ is

$$[eqn]$$

where $[eqn]$ is the effective diffusivity, and $[eqn]$ is the imposed advection field; in practice we take $[eqn]$ and advection is imposed solely in the direction normal to the fenestra through $[eqn]$ . This imposed advection does not represent resolved convection within the ECM bulk; rather, it serves as a phenomenological representation of the transport-driving influence of plasma entry at the vessel–ECM interface. The RAD model is termed “CFD-informed” as the magnitude of $[eqn]$ is derived directly from vessel-scale multiphase CFD simulations. Specifically, fifteen independent tumor-vessel CFD simulations are performed with varying fenestra locations. From each simulation, the magnitude of the plasma entry velocity at the fenestra is extracted. A single characteristic advection scale is then defined as the arithmetic ensemble average of these velocity magnitudes across all vessel geometries, and this value is used to prescribe $[eqn]$ in Equation 16. No temporal averaging is performed, and no spatial velocity field is transferred into the ECM domain. Thus, CFD-informed here denotes parameterization of transport strength using CFD-derived entry conditions, not volumetric or transient coupling between CFD and ECM models. The RAD model is posed on a reduced-order rectangular ECM domain of size $[eqn]$ with a $[eqn]$ fenestra through which plasma enters. A total of 40 circular fiber inclusions are placed with minimum spacing (no overlap) to maintain a packing fraction $[eqn]$ (Andersson et al., 2008; Nier et al., 2016; Spatarelu et al., 2019; Yeasin, 2025), yielding an average fiber radius $[eqn]$ . The geometry (Figure 8A) is reconstructed in Wolfram Mathematica and meshed (Figures 8 B,C) via ToElementMesh into an unstructured triangular grid ( $[eqn]$ elements), providing sufficient resolution near the fenestra and fiber boundaries. We initialize with $[eqn]$ and apply a constant inlet concentration of $[eqn]$ —twelve orders of magnitude above the physiological protein concentration in plasma ( $[eqn]$ (Leeman et al., 2018))–at the fenestra. To match CFD timescales in this reduced domain, A constant inlet concentration is prescribed at the fenestra, while zero-flux (Neumann) conditions are enforced on all fiber surfaces and lateral boundaries to prevent artificial back-diffusion. The RAD model is designed to reproduce the spatiotemporal progression trends of plasma transport observed in CFD and microfluidic experiments, rather than absolute concentration levels. Accordingly, inlet concentration magnitude and diffusivity ( $[eqn]$ (Muraoka et al., 2008)) are treated as effective parameters calibrated to match CFD-resolved transport timescales.

Finally, to verify consistency between the theoretical RAD framework and full numerical transport simulations under identical boundary forcing, a two-phase Eulerian CFD simulation was performed on the same reduced geometry using ANSYS Inc, 2023 R1. The domain was constructed in ANSYS Inc, 2023 R1 and discretized with an unstructured triangular mesh (approximately $[eqn]$ elements) with local refinement near the fenestra and along the patched-wall segments. To accommodate the concentration inlet prescribed by the theoretical model, only the outlet boundary condition was modified for this validation case; all other settings (solver options, material properties, numerics) were kept identical to Section 2.2.2. Geometry and boundary placements were identical to those used for the RAD domain. This validation exercise does not constitute coupling between models, but rather demonstrates that the simplified RAD formulation reproduces key transport trends when driven by the same fenestral entry conditions.

Results

3

Comparison of plasma perfusion from vessel to ECM with and without EHD

3.1

We quantify how the electrohydrodynamic (EHD) force in the tumor vessel modifies the plasma delivered through the fenestra into the ECM, relative to an otherwise identical non-EHD baseline. Using the representative geometry (Model 5), we select the same region of interest (ROI)–the fenestra–in both simulations (Figures 4A,B). Pixel colors in the ROI are mapped to a dimensionless plasma intensity $[eqn]$ producing two paired empirical distributions over the identical ROI: $[eqn]$ and $[eqn]$ . The histogram comparison in Figure 4E shows a decisive rightward shift with EHD. Denote the sample means by $[eqn]$ and $[eqn]$ , and standard deviations by $[eqn]$ and $[eqn]$ and let $[eqn]$ be the pooled standard deviation. The standardized mean difference (Cohen’s $[eqn]$ ) (Cohen, 2013) is

$$[eqn]$$

EHD flow yields enhanced, more persistent plasma penetration at the fenestra compared with non-EHD flow (A,B) Plasma (primary phase) contour inside the tumor vessel with EHD (A) and without EHD (B). (C,D) Zoomed fenestra ROI from panels A and B (same simulation Model 5). The dashed rectangle delineates the identical physical area used for pixel counting and threshold analysis. (E) Fenestra-ROI intensity distributions show a clear rightward shift under EHD with minimal overlap, indicating stronger local plasma penetration; mean intensities 0.650 (EHD) and 0.525 (non-EHD). (F) Threshold (survival) curves summarize robustness across cutoffs the EHD curve dominates across thresholds, and the area-between-curves (ABC) equals 0.125 (Equations 21-22), demonstrating more persistent high-intensity pixels only with EHD.

For Model 5, $[eqn]$ increases from 0.525 (non–EHD) to 0.650 (EHD), yielding $[eqn]$ using Equation 17 (“very large” separation). A threshold-free dominance metric is the probability of superiority (AUC),

$$[eqn]$$

from which Cliff’s delta follows as

$$[eqn]$$

Here, $[eqn]$ (from Equation 18) and $[eqn]$ (from Equation 19) indicate minimal overlap between conditions. A robust location shift is summarized by the Hodges–Lehmann median (Equation 20) difference,

$$[eqn]$$

To assess areal coverage at clinically relevant intensities, we examine the exceedance (survival) curves (Figure 4F) over a threshold $[eqn]$ :

$$[eqn]$$

Figure 4F shows $[eqn]$ (black) lies strictly above $[eqn]$ (orange) across the inspected range, i.e., larger high-intensity area for any chosen cutoff. The area between curves (ABC) provides a cutoff-agnostic summary:

$$[eqn]$$

At a stringent threshold of $[eqn]$ , the high-intensity fraction is $[eqn]$ (EHD) versus $[eqn]$ (non-EHD), clearly showing a strongly perfused area under EHD. In Figure 4, panels E–F jointly indicate that EHD (via patched, glycocalyx-mediated wall forcing) increases both how much of the fenestra area attains high plasma fraction (larger exceedance across $[eqn]$ ) and how strongly it is filled (right-shifted intensity distribution). These inlet-side gains are exactly the mechanism that yields the faster and broader percolation observed downstream in the ECM domain (Section 3.2).

Structural sensitivity: across all 15 models, Figures 5A,B demonstrates that the mean calibrated intensity under EHD is consistently higher than under Non–EHD. Figure 5A shows this pattern model–by–model, while panel B summarizes the across–model effect: grand means of 0.722 (EHD) versus 0.576 (Non–EHD), with a mean improvement $[eqn]$ (95% CI $[eqn]$ ; Wilcoxon signed–rank $[eqn]$ ). Taken together, these results indicate that the EHD-induced enhancement of plasma extravasation is robust to changes in fenestra location and vessel geometry, rather than being specific to a single configuration. Our objective here is to establish the condition effect (EHD vs. Non–EHD), not to optimize fenestra placement; accordingly, we do not rank or compare the EHD models against one another by location.

Structural sensitivity of EHD-induced plasma extravasation across vessel geometries. EHD flow consistently yields higher mean plasma penetration than non-EHD across all 15 models: (A) For models 1–15 (fenestra center translated from near inlet to near outlet), grouped bars compare Non–EHD and EHD mean calibrated intensity for each model; higher values indicate greater plasma penetration. EHD exceeds Non–EHD in every model. (B) Per–model means (dots) with grand mean and 95% bootstrap CI (square with whiskers) for each condition. Across models, mean Non–EHD =0.576 and mean EHD =0.722 (mean Δ=0.146 , 95% CI [0.124, 0.166] ; Wilcoxon signed–rank p=6.1×10−5 ).

Plasma percolation inside the tumor ECM domain

3.2

With the tumor-domain inlet fields $[eqn]$ constructed from the $[eqn]$ vessel simulations (Section 2.2.2; Equation 15), the flow advances the air–plasma interface through the tumor ECM domain. The spatiotemporal evolution exhibits a coherent percolation pattern: an early thin, inlet-aligned jet penetrates along preferential low-resistance pathways; intermediate times show lateral broadening and progressive saturation of interstitial voids; and late times approach a nearly uniform, high–plasma-fraction interior behind the front. This progression is visualized in Figures 6A–D via contours of local plasma volume fraction, where red denotes the advancing high-fraction front, yellow/green the graded rim, and blue the not-yet-perfused interstitium. To quantify front advance and rim thickening, we compute an image-based color occupancy ratio for each frame by normalizing the blue (unfilled) pixels to $[eqn]$ and reporting the triplet $[eqn]$ after excluding white pores from the counts. A rise in $[eqn]$ signals deeper filling, and a rise in $[eqn]$ signals a thickening front. As summarized in Table 2, $[eqn]$ grows monotonically from $[eqn]$ at $[eqn]$ sec to $[eqn]$ at $[eqn]$ sec, with a concurrent increase in $[eqn]$ . Consistent with the visual maps, early times are advection-dominant (small $[eqn]$ , thin inlet-aligned plume); by mid-times the rim thickens as cross-stream spreading increases; at late time most accessible voids are filled (blue nearly exhausted) and the field behind the front approaches homogeneous red.

Spatiotemporal spreading of plasma in the tumor ECM shows an advection-to-diffusion transition in wetted area growth: (A–D) Contours of plasma volume fraction at (A) 0.01 s, (B) 0.05 s, (C) 0.10 s, and (D) 0.20 s. White circles denote impermeable fibers; blue is unfilled interstitium. The 0–1 color bar maps local plasma fraction. The wetted region expands from a thin inlet plume to a broad front, illustrating a qualitative shift from advection-driven advance toward diffusion-controlled spreading. (E) Wetted area versus time on log–log axes. The red dashed line marks the transition time t=1.665×10−1 sec , separating the advection-dominant region (left) from the diffusion dominant region (right).*

TABLE 2: Color-occupancy ratios for Figures 6A–D. Ratios are computed from pixel counts with blue normalized to 1; white pore pixels are excluded. Increasing R:B and (G+Y):B indicate deeper penetration and a thickening front.

A complementary, domain–integrated metric is the wetted area starting with $[eqn]$ (Figure 6E) of non-blue pixels. On log–log axes we analyze $[eqn]$ so that a power law $[eqn]$ appears as a straight line:

$$[eqn]$$

with exponent $[eqn]$ . A two–segment least–squares fit in $[eqn]$ identifies a break point $[eqn]$ that minimizes the total residual sum of squares. Using this procedure we obtain $[eqn]$ , which partitions the record into advection region: t $[eqn]$ sec, and diffusion dominant region: t $[eqn]$ sec. Within any $[eqn]$ , ordinary least squares gives (see Table 3)

$$[eqn]$$

TABLE 3: Segment-wise fits on (log10t,log10A) for log10A=m log10t+b .

To summarize the transition at $[eqn]$ , we report the one–sided model predictions, $[eqn]$ . The advection side extrapolation is therefore higher at the breakpoint (offset $[eqn]$ relative to the diffusion-side value; $[eqn]$ relative to the mean of the two), which is typical for an unconstrained two–line fit and indicates an abrupt change in kinetics rather than a perfectly continuous crossover. The effective exponent drops from super–linear to subunit from Equations 23-24 ( $[eqn]$ to 0.41), consistent with a shift from advection–driven growth to a diffusion–controlled late regime; this interpretation agrees with the RGB composition (persistent graded rim).

Experiment–CFD validation of wetted-area growth

3.3

Top view images (Figures 3D,E) from the experiments were segmented to obtain wetted pixel counts over time and converted to physical area, yielding $[eqn]$ in Figure 7E. The experiment exhibits advection dominated growth overall, with an early shoulder consistent with $[eqn]$ before a clear linear trend $[eqn]$ becomes established; we identify the knee at $[eqn]$ . For a like-for-like benchmark, we ran a separate 2D CFD configuration (distinct from Section 2.2.2) using a pressure-driven inlet to match the experimental $[eqn]$ ; plasma and air properties were unchanged, and the front evolution appears in Figures 7A–C. Applying the same area-extraction to the simulation yields the same qualitative behavior; advection-dominated growth with a brief early shoulder—and a much earlier knee at $[eqn]$ (Figure 7D). Thus, both experiment and CFD are advection-dominated; the primary discrepancy is the transition timescale ( $[eqn]$ versus $[eqn]$ ), which we analyze in Section 4.4 in terms of inlet start-up/wetting effects and geometric idealization.

CFD and experiments both show advection-dominated growth regimes in ECM wetting, with different transition times: (A–C) Snapshots show front advance at t=0.01, 0.10, 0.20,sec respectively. (D) log–log growth of area A(t) with an advection-dominated regime transitioning to diffusion at tCFD=0.2277 sec . (E) Early diffusion-dominated growth (A∼t1/2) transitioning to advection (A∼t) at texp*=4.1 sec .*

Comparison of theoretical model predictions with corresponding numerical simulation outcomes

3.4

The reverse advection–diffusion (RAD) model reproduces the same two-stage transport seen in the like-for-like CFD of the simplified rectangular domain: an inlet-aligned, advection-driven entry plume followed by diffusion-dominated homogenization as the interior fills (theory snapshots in Figures 8D–H; CFD snapshots in Figures 9A–D). To quantify kinetics, we track the occupied interstitial area $[eqn]$ for both models on log–log axes (Figures 9E,F). In each case the slope transitions from a steep, advection-driven regime to a shallow, diffusion-controlled regime, with breakpoints at $[eqn]$ (Figure 9E) and $[eqn]$ (Figure 9F).

CFD-informed reverse advection-diffusion RAD model showing solute plume evolution through a reduced-order fiber-dense tumor ECM: (A) The simplified rectangular computational domain (10×5 μm2) used for the reduced-order RAD model, containing randomly distributed circular ECM fiber bundles occupying 30% of the rectangular domain area. All fibers have equal radius and act as impermeable structures, mimicking dense tumor stroma and restricting solute motion. (B) The high-resolution unstructured triangular finite-element mesh generated for the theoretical analysis. Whereas (C) provides a zoom view of the mesh in the fenestra region, highlighting refined discretization around fiber boundaries to resolve steep concentration gradients. (D–H) The predicted spatiotemporal evolution of plasma-borne solute under the advection-diffusion transport model at t=0.1 sec, 0.5 sec, 1.0 sec, 2.0 sec, and 3.0 sec , respectively, where C~ represent the normalized concentration field.

Comparison of CFD and RAD model for plasma transport transition: (A–D) show the CFD-predicted plasma propagation at t=0.1,1,2 , and 3 sec , respectively. Plasma first penetrates swiftly along preferential pathways between fibers, then progressively spreads laterally and fills the matrix, reflecting a transition from advection-dominated entry to diffusion-driven equilibration. (E) Shows the occupied interstitial area A(t) as a function of time for the CFD solution, demonstrating a steep early-time growth regime followed by a slowed diffusion-controlled regime as the domain fills. (F) Shows the corresponding theoretical RAD model prediction for A(t) , showing the same trend as the CFD solution.

We summarize the post-transition timescale offset by the simple ratio

$$[eqn]$$

indicating that the theoretical curve evolves about 1.295 times more slowly than the CFD result on this domain. If we rescale the theoretical time axis by this factor,

$$[eqn]$$

it will bring the $[eqn]$ trajectories into close agreement in both the early advection regime and the late diffusion-dominated regime (Equations 25–26). Thus, the reduced-order RAD model captures the same mechanistic shift observed in CFD and, after a single temporal calibration, serves as an accurate and computationally efficient surrogate for fully resolved flow simulations on this geometry.

Discussion

4

Diffusion dominant late regime with subunit area scaling

4.1

In the late regime ( $[eqn]$ sec)), wetted–area growth is governed by diffusion even though the log–log exponent is subunit $[eqn]$ . A persistent graded rim (green + yellow) leading the bulk (red), quantified by slowly varying color ratios (Table 4), is consistent with diffusion–limited penetration through a tortuous, fiber–dense matrix. In diffusion–limited flow, the characteristic penetration length obeys $[eqn]$ , so observables that scale with length (e.g., a band–like front across a roughly fixed width) grow as $[eqn]$ , not linearly (Crank, 1979). Consequently, the area can grow sublinearly $[eqn]$ and still be diffusion–controlled. In ideal band–like geometries one expects $[eqn]$ $[eqn]$ ; values below 0.5 are expected in complex matrices where tortuosity, dead–ends, or partial channelization slow the effective front (subdiffusion) (Sahimi, 2011; Metzler and Klafter, 2000; Berkowitz et al., 2000; Edery et al., 2014).

Temporal evolution of plasma penetration inside tumor ECM domain at different flow times: Panels show snapshots of plasma perfusion inside the tumor domain at (A) t = 0.1666 s, (B) 0.1800 s, (C) 0.1900 s, and (D) 0.2000 s within the diffusion-dominated region. Red denotes the high-fraction bulk, green + yellow denote the graded rim, and blue denotes the unperfused region. Panels show a persistent graded rim leading a slowly advancing bulk with a gradual recession of blue—indicates gradient-limited (diffusive) penetration through a tortuous, fiber dense matrix, consistent with subunit area-growth scaling.

Independent image–based composition metrics corroborate a diffusion–controlled tail. Post $[eqn]$ snapshots show a persistent graded rim (green + yellow) leading the high–fraction bulk (red) while the unperfused region (blue) depletes gradually. Quantitatively, the $[eqn]$ increases only slowly from 4.63 to 6.75 between t = 0.1666 and 0.2000 sec, and the $[eqn]$ remains appreciable (0.161–0.130). Crucially, the rim–to–bulk ratio $[eqn]$ stays non-vanishing and declines only mildly $[eqn]$ , indicating a broad, gradient–limited front rather than a sharp advective discontinuity at high Péclet number.

Modeling assumptions for blood rheology

4.2

In the present framework, red and white blood cells are represented as dispersed phases to capture phase fraction distributions, and electrohydrodynamic (EHD) effects near fenestrae. However, the microscale mechanisms responsible for blood’s non-Newtonian behavior—such as RBC deformation, aggregation, and membrane dynamics—are not explicitly resolved. Accordingly, the Carreau–Yasuda viscosity law is introduced phenomenologically as a constitutive closure to represent the net rheological response observed experimentally, rather than as a property emerging self-consistently from the multiphase dynamics. This separation of transport representation and constitutive rheology avoids double-counting and is consistent with established continuum blood flow modeling practices.

Modeling assumptions and interpretation of plasma transport inside ECM region

4.3

Within the current framework, plasma transport within the tumor extracellular matrix (ECM) is modeled using a two-phase transport formulation to provide a well-defined reference state for quantifying plasma entry and penetration from the vessel wall. Although the in vivo ECM is continuously saturated with interstitial fluid, the air phase used here should be interpreted as a numerical placeholder representing the resident interstitial volume rather than a physical gas phase. In this sense, the model captures liquid–liquid replacement behavior, with plasma displacing an initially present fluid, but does so through an immiscible front-tracking formulation rather than a concentration-based mixing model. This approach enables clear identification of transport initiation and penetration depth, avoids numerical ambiguity associated with weak concentration gradients, and facilitates consistent comparison across vessel realizations and inlet conditions. Consequently, the resulting plasma volume fraction and penetration metrics should be interpreted as relative indicators of ECM perfusion efficiency and transport trends, and are generalizable in a comparative sense rather than as direct predictions of absolute in vivo interstitial solute concentrations. Future extensions of the framework may incorporate fully saturated liquid–liquid mixing formulations to recover detailed concentration fields while retaining the multiscale vessel–ECM coupling strategy developed in this work.

On the experiment–CFD timescale discrepancy and thickness scaling

4.4

We first quantify the temporal discrepancy by comparing the post-transition analysis. The corresponding timescale ratio is

$$[eqn]$$

This result indicates that the experimental timescale is approximately 509 times (from Equation 27) longer than the computational timescale, emphasizing the profound effect of physical constraints present in the 3D experimental setup but absent in the 2D numerical model.

Under the constant–flow (pump–driven) conditions of the experiment, mass conservation in a thin layer of thickness $[eqn]$ gives $[eqn]$ , so the characteristic time obeys $[eqn]$ at leading order. Interpreting the 2D CFD as a per–depth surrogate of effective thickness $[eqn]$ , the transition–time ratio therefore scales as

$$[eqn]$$

With $[eqn]$ = $[eqn]$ and the measured ratio $[eqn]$ , we infer from Equation 28

$$[eqn]$$

which explains the $[eqn]$ gap directly from the 3D–to–2D thickness difference. To reconcile the curves graphically (Figure 11), we rescale CFD time and area (from Equations 28–29) by

$$[eqn]$$ $$[eqn]$$

and plot the residual as (using Equations 30–31)

$$[eqn]$$

Scaling reconciliation of wetted-area growth between experiment and CFD: Wetted-area dynamics align between experiment and CFD after time and thickness scaling. Red dashed: experiment vs. raw CFD (before scaling). Blue solid: experiment vs. scaled CFD (after scaling). Post-scaling, the offset is nearly constant and the trends coincide, indicating matched advection-dominated growth.

After an initial transient, the residual flattens (blue curve), indicating consistent advection–dominated growth once thickness and time are reconciled; any remaining offset is attributable to early–time device effects (pump start–up/compliance, wetting/pinning and capillary entry on the resin surface, transient air evacuation) that are absent from the idealized 2D model.

On the use of a large inlet concentration in the RAD model

4.5

The reverse advection-diffusion (RAD) model matches CFD on both the advective and diffusion regimes, indicating that key transport physics are retained without full Navier–Stokes resolution; hence RAD is well-suited for rapid sweeps and design use. To theoretically model the diffusive component for the perfusion dynamics within the extracellular domain, we have imposed a finitely large value ( $[eqn]$ in model units) for plasma concentration at time $[eqn]$ as a pragmatic approximation (Crank, 1979) mimicking a persistent, high concentration fenestral ‘point’ source. A sufficiently large (but finite) value for initial inlet concentration creates a steep concentration gradient that effectively drives diffusion while maintaining the stability and robustness of the analytical solver. This approach aligns with standard practice in transport modeling, where high boundary concentrations are employed as Dirichlet boundary conditions to represent sustained sources, thus allowing for a realistic yet computationally feasible launch point for the complex intra-tumoral transport process without sacrificing physical relevance.

The present surrogate intentionally simplifies geometry and driving fields to isolate advection–diffusion obstruction physics: an idealized rectangle with uniformly sized, non-contacting circular fibers and a spatially simple inlet advection profile. Future work will incorporate histology-derived fiber statistics (size distribution, alignment, and contacts) and spatially varying $[eqn]$ , informed by localized CFD or microfluidic measurements, thereby extending RAD from rapid screening toward tissue-specific predictions while preserving computational tractability.

Specimen-specific tumor vessel-ECM geometry reconstruction: future model extensions

4.6

While the present study employs idealized vessel–ECM geometries to enable controlled mechanistic interrogation of plasma transport and electrohydrodynamic effects, we have also developed a histology-to-geometry reconstruction pipeline intended for future specimen-specific modeling. Representative H&E-stained tumor sections (Figures 1B, 2C) were obtained from pancreatic tumor tissue, which was collected, fixed in formaldehyde for 24 h, paraffin embedded, and sectioned at 5 $[eqn]$ m. Sections were deparaffinized in Histo® Clear and ethanol, rehydrated, stained with hematoxylin (5 min) and eosin (1 min), dehydrated in 95%–100% ethanol, cleared in xylene, mounted, and imaged using a Zeiss Axio Observer Z1 microscope, following the same histological preparation and imaging protocol described in our prior work (Akash et al., 2023). An in-house image-analysis framework was applied to these histological sections to detect the luminal periphery using intensity-based edge detection combined with morphological cleanup. The extracted edge contours were converted into ordered coordinate sets along the vessel wall and imported into a CAD environment, where they were aligned to a common reference frame. A spline-based interpolation procedure was then used to generate a continuous vessel-wall representation with $[eqn]$ continuity, reducing pixelation artifacts inherent to raster-based images while preserving geometric fidelity. This pipeline enables reconstruction of curved vessel–ECM interfaces directly from histological data (Figures 12A,B) and supports the generation of adjacent ECM domains populated with fiber bundles at prescribed packing fractions (Figure 12C). Importantly, this histology-informed reconstruction capability is presented here as a methodological extension and is not utilized in the numerical simulations reported in the present study. The simulations intentionally rely on idealized geometries to isolate electrohydrodynamic and multiphase transport mechanisms, thereby avoiding confounding anatomical variability. Integration of specimen-specific vessel morphologies derived from this pipeline into the multiscale transport framework represents a natural and ongoing direction for future investigation.

Histology-to-geometry reconstruction pipeline for future specimen-specific vessel–ECM modeling. (A) H&E histology of an orthotopic human pancreatic tumor slice illustrating irregular microvasculature; this specimen informed the vessel geometry. (B) Imaging-based vessel extracted from the tumor slice and used to define the curved vessel–ECM interface. (C) Tumor ECM domain (3D view) populated with circular fiber bundles to achieve the prescribed packing fraction; the curved red vessel indicates the histology-informed vessel–ECM boundary used to locate the fenestra.

Scope of the study: fluid-mechanical foundations vs. pharmacological modeling

4.7

It is important to clarify the scope and interpretation of the present study. This work does not aim to predict drug delivery or molecular pharmacokinetics within tumor tissue. Instead, it focuses on the fundamental fluid-mechanical mechanisms governing plasma perfusion and transport within the tumor extracellular matrix (ECM). Plasma is treated as the carrier phase whose ingress, penetration, and redistribution establish the baseline transport capacity of the tumor extracellular matrix (ECM), independent of any specific solute payload. Accordingly, plasma volume fraction, penetration depth, and front advancement are interpreted as mechanistic indicators of ECM perfusion efficiency rather than as direct surrogates for drug concentration or therapeutic efficacy. By isolating plasma transport physics, the framework enables systematic interrogation of how vessel-scale flow conditions and localized electrohydrodynamic effects bias interstitial transport pathways, providing mechanistic insight that can later inform, but does not itself constitute, drug-delivery modeling.

The main takeaways

5

We integrate three-phase, glycocalyx-patched (EHD) CFD with a histology-informed ECM domain to quantify how endothelial EHD alters plasma delivery at a fenestra and its downstream percolation. We then use a CFD-informed reverse advection–diffusion (RAD) model as a corroborative, geometry-aware surrogate for rapid exploration. The key findings are that across all 15 fenestra positions, adding EHD increases the inlet-side mean calibrated plasma intensity from 0.576 (Non-EHD) to 0.722 (EHD), representing a $[eqn]$ gain. Downstream, both the ECM simulations and a microfluidic benchmark exhibit the same two-stage transport (advection-dominated entry followed by diffusion-driven spread), and the CFD-informed RAD reproduces these kinetics successfully.

In future work, we will extend the ECM-domain simulations to fully three-dimensional geometries, capturing transverse leakage paths and thickness-dependent effects that are absent in 2D. The reverse advection–diffusion model will be augmented with explicit electrohydrodynamic forcing (glycocalyx-patch charge and field–flow coupling) to incorporate EHD into the closures for effective advection and dispersion. Experimentally, we will develop refined organ-on-chip platforms with tunable wall charge and co-cultured endothelium–stroma to control EHD conditions in vitro, enabling direct quantification of perfusion trends under EHD versus Non-EHD and more rigorous validation of the modeling framework. Together, these advances will broaden biological realism and improve the predictive scope of EHD-aware tumor perfusion modeling.