Modeling tumor transport and growth with poroelastic biopolymer networks

Abstract

The mechanical properties of the extracellular matrix (ECM) regulate tumor growth and invasion in the tumor microenvironment. Models of biopolymer networks have been used to investigate the impact of elasticity and viscoelasticity of the ECM on tumor behavior. Under tumor compression, these networks also show poroelastic behavior that is governed by the resistance to water flow through their pores. This work investigates the hypothesis that stress-dependent transport properties of biopolymer networks regulate tumor growth. Here, alginate hydrogels are used as a model ECM system with tunable ionic and hybrid ionic/covalent crosslinking. Hydrogel stiffness, viscoelasticity, and stress relaxation behavior were characterized using stepwise axial compression. Among these properties, we find poroelastic fluid outflow dominates ECM stress relaxation, as the measured water flux was significantly affected under compression. Continuum mechanics-based modeling was developed to formulate and calculate the chemical potential gradients of water (solvent) in the hydrogels under compression. This framework was extended into an advection–diffusion framework to quantify growth factor (solute) distribution under varying strengths of stress and diffusion indexed by the relative strength of convective to diffusive transport, characterized by the Péclet number. An agent-based computational simulation showed that the Péclet numbers based on our experimental timescales strongly influenced tumor growth over longer, more physiologic timescales. Together, these results highlight the important role of water flux and transport in three-dimensional biopolymer networks.

---

1 Introduction

The extracellular matrix (ECM) plays a critical role in regulating cell behavior, including proliferation, migration, and differentiation, through both biochemical and mechanical cues. Mechanical cues of ECM arise from the structure and properties of biopolymer networks and their interstitial fluid flow. Stiffness, viscoelasticity, and plasticity significantly impact cancer cell mechanobiology, as cancer cells sense and respond to variations in mechanical forces determined by these mechanical properties.^1–6 Interstitial fluid flow within the network's pores regulates cancer cell migration and transport of nutrients, including metabolites, proteins, and cytokines, which have downstream effects on cancer cell behaviors.^7–9 Under deformation, both the biopolymer network and interstitial fluid flow contribute to the ECM's stress relaxation properties. Biopolymer networks relax stress via structural reorganization (viscoelasticity) and fluid efflux from the porous network (poroelasticity). Viscoelastic properties correspond to shorter timescales (0.1–100 s), while poroelastic relaxation manifest at longer timescales (>1000 s).¹⁰ Oscillatory shear rheology and stress relaxation with small deformation (≤10% strain) and short timescales (0.1–100 s) are used to characterize viscoelasticity of biopolymer networks.^2,11,12 However, the compressive poro-mechanical behavior of biopolymer networks under large uniaxial or biaxial deformations—on the order of 30–40%, which is comparable to strains reported in tumor spheroid experiments—and over longer timescales (1000–100 000 s) remains understudied.^13,14 Such strain levels are physiologically relevant, as in vitro studies have shown that tumor spheroid growth in agarose gels can generate strains of approximately 30–40% and stresses exceeding 10 kPa, while measurements in glioma cases report strains in surrounding brain tissue of up to about 24%.¹⁵

Tumor proliferation pushes radially outward on the surrounding tumor microenvironment (TME), leading to self-compression due to physical confinement.^5,16 Compressive stress experienced by solid tumors impacts tumor growth dynamics, morphology, and therapy resistance, ultimately shaping disease progression. In responses to compressive stress, tumor cells intrinsically undergo glycolytic reprogramming, enter quiescence that can impair drug efficacy, and activate mechanosensitive YAP/TAZ, PI3K–AKT, and Wnt–β-catenin pathways, collectively enhancing survival and invasiveness.^17–20 Solid stresses during compression can cause collapse of blood vessels and affect lymphatic drainage.^16,21 Compression during tumor growth also generates a pressure gradient and outward convective flux.²² Solid-phase compression and interstitial fluid flow are intrinsically coupled by the poroelastic nature of the ECM.^23,24 Under large-scale, slow deformation, poroelastic models show interstitial fluid flow that influences the availability and distribution of important solutes like oxygen, growth factors, nutrients, and cytokines mediating tumor-immune interaction.^18,25 Thus, investigating the poroelastic properties of tumor ECM may provide new strategies to mitigate compression-driven resistance, as well as improve predictive models of growth factor transport and availability in the tumor microenvironment (TME).²⁰

We hypothesize that poroelastic properties govern tumor behavior at relevant scalings of tumor growth compared to the faster and shorter length-scales of viscoelastic properties.²⁶ We developed a poroelastic experimental, theoretical, and computational framework for characterizing compression-induced fluid flux in advection–diffusion simulations to quantify how deformation-driven fluid flow modulates growth factor distribution. Alginate provides a natural polysaccharide-based network model system that inherently demonstrates compressive strain-stiffening properties,^27,28 which is a key property of native ECM.^29,30 Alginate's nanoporous structure mimics the restrictive transport environment of the TME, which limits interstitial fluid flow,²⁴ compared with collagen-based systems that lack the ability to control growth factor transport. Alginate has been studied widely for its ability to control growth factor delivery for regenerative and cancer therapy applications.^31,32 In this work, we used stepwise large (up to 40% strain) compression of an alginate hydrogel to examine ECM deformation and outward fluid flow of tunable biopolymer networks. A customized rheological setup and finite element analysis (FEA) were utilized to evaluate the volumetric flux and chemical potential gradient of fluid within the polymer network, providing a basis for fluid dynamics within the matrix for solute transport phenomena and tumor tissue agent-based modeling. Agent-based modeling (ABM) was used to simulate the spatiotemporal dynamics of individual entities, or agents, interacting within a shared environment. In the context of solid tumor growth, each agent represents a cancer cell, which resides within a domain representing the TME. Cancer cells proliferate based on user-defined growth rates.^33,34 In this study, the ABM provided direct quantification of growth factor distribution-mediated differences in tumor cell growth. We used the computed growth factor concentration profiles under compression to simulate tumor growth dynamics using an agent-based modeling platform. Overall, our results offer key insights into the role of compressive stress-induced interstitial flow in ECM remodeling during tumor progression.

2 Results and discussion

2.1 Fluid flow dominates in biopolymer networks with large and slow deformation

Stress relaxation behavior reveals the mechanisms of ECM relaxation through its two components: the biopolymer network and the interstitial fluid. The biopolymer network relaxes stress via structural reorganization (viscoelasticity), while the interstitial fluid contributes through fluid transport and efflux from the porous network (poroelasticity).²⁶ Stepwise compression, involving a series of 10% compressions, 2500 s of stress relaxation, and time sweep after complete relaxation, was used to measure relaxation time constants, storage moduli, and tan(δ) (viscoelasticity). Relaxation time constants were obtained by fitting stress–relaxation curves to an exponential decay model based on the generalized Maxwell model, while storage moduli and tan(δ) were directly measured from time sweep rheology (Fig. 1A). Ionic-crosslinked alginate (via calcium ions) and covalent-crosslinked alginate (via Nb–Tz click chemistry) were tested under stepwise compression to assess their relaxation time constants. Unconfined compression is used in both our experiments and the corresponding FEA simulations. This choice reflects the fact that, in many tumors, the surrounding host tissue is relatively compliant compared with the tumor itself. Numerous studies have shown that solid tumors often exhibit higher stiffness than the adjacent normal tissue due to extracellular matrix remodeling and collagen deposition.^35,36 After fast compression (strain rate of 0.25% s⁻¹) from 0 to 10% strain, the relaxation time constant of covalent alginate was approximately three times higher than that of ionic alginate, indicating faster relaxation in ionic alginate. Under larger strain compressions, covalent alginate exhibited distinct trends and significantly different relaxation time constants compared to ionic alginate (Fig. 1B). These differences likely arise from the nature of crosslinks, as chemical hydrogels with rigid covalent bonds are less likely to break and rearrange and thus relax slower than physical hydrogels.²⁶ However, after slow compressions (strain rate of 0.025% s⁻¹), both covalent alginate and ionic alginate showed similar relaxation time constants and trends (Fig. 1C).

Fig. 1 Fluid flow dominates stress relaxation behavior during large and slow deformations. (A) Stepwise compression and stress relaxations with time sweep, which quantifies relaxation time constants, storage moduli, and tan(δ). (B) and (C) Relaxation time constants at various strains for ionic and covalent crosslinked alginate (Alg) after fast, 0.25% s⁻¹ strain rate (B) and slow, 0.025% s⁻¹ strain rate (C) compressions. (D) Compression on gel disc generates the outward convective flux, which is measured by a customized rheometer setup via recording volume change over time. (E) Ionic alginate's flux and pore pressure difference change during stress relaxation under 10–40% compressive strain. (F) Permeability of ionic alginate as compressive strain increases.

To experimentally measure fluid flow, a customized rheometer setup was used to obtain volume changes during compression (Fig. 1D) of the alginate hydrogel, which is used to calculate volumetric flux of the fluid, q.³⁷ As compressive strain increases, fluid flux becomes increasingly negative, indicating fluid is leaving the matrix at a faster speed. The measured flux of 0–0.6 µm s⁻¹ (Fig. 1E) matches in situ interstitial fluid speed in neoplastic tissue (0.5–0.6 µm s⁻¹),²⁴ suggesting our alginate hydrogel and testing system is a reasonable model for ECM in tumor microenvironment. Darcy's law, , quantifies volumetric flux, q, of the fluid (with viscosity, η) moving through a matrix with permeability, k, under a pore pressure gradient, Δμ, over distance L. Using Darcy's law, ionic alginate's permeability for each stress relaxation was determined by linearly fitting flux versus pore pressure difference (Fig. 1E). The permeability of ionic alginate decreased with increasing compressive strain (Fig. 1F), consistent with network densification. As the alginate disc is compressed and fluid is expelled, the polymer network becomes denser, making fluid permeation more difficult. Oscillatory shear rheology of the ionic gel was performed after the normal stress had fully relaxed (before each 10% compression step). Shear storage moduli and tan(δ) of the ionic alginate slightly increased, which was consistent with densification of the matrix (Fig. S1). These results together demonstrate the poroelastic behavior of alginate hydrogel networks.

2.2 Water flux varies non-linearly as a function of compressive strain

To better understand compression-induced fluid flow dynamics, finite element analysis was conducted to compute the chemical potential (pore pressure) of the fluid within the network. The finite element model (FEM) was first constructed and validated using experimental stress–strain curves at three different strain rates (Fig. S2), which showed quantitative agreement between the model and experiments. The chemical potential, μ, is defined as μ = p − π, where p is the local hydrostatic pressure and π is the osmotic pressure. Since the hydrogel disc is confined between two impermeable plates, fluid can only exit through the lateral surface of the gel cylinder during compression. At the boundary (r_edge), fluid escapes easily due to the absence of network confinement in the direction of flow. In contrast, fluid at the center (r₀), remains trapped due to confinement by the polymer network in the radial direction (Fig. 2A). As a result, the chemical potential varies with radial distance, forming a gradient, , under compression. Volumetric flux at the boundary, q(r_edge), can be calculated from Darcy's law, , with the viscosity of fluid, η, and permeability, K (Fig. 2B).

Fig. 2 Nonlinear outward convective fluid flow under compression. (A) Schematic of a high chemical potential fluid generated at the center of the hydrogel disc. (B) Computational result: chemical potential gradient under compression. (C) Computational result: chemical potential builds up as strain increases. (D) Computational result: change in chemical potential gradient at gel boundary, , under increasing strain. (E) Computational result: rate of change in chemical potential gradient, , at gel boundary under increasing strains. (F) Experimental result: stage 1: volumetric flux and chemical potential gradient get more negative faster, meaning fluid speed leaving the gel increases at a higher rate. (G) Experimental result: stage 2: volumetric flux and chemical potential gradient decrease more slowly, meaning fluid speed leaving the gel increases at a lower rate.

As strain increases, fluid chemical potential builds up (Fig. 2C). Due to the sudden generation of high chemical potential at the gel center from zero chemical potential at 0% strain, the chemical potential gradient at the boundary showed a fast decrease at the beginning of compression (stage 1) and followed by a slow and steady decrease at stage 2, showing a nonlinear water efflux under compression (Fig. 2D and E). A more negative gradient indicates a greater chemical potential difference, leading to fluid leaving the hydrogel at a higher speed. The finite element simulation aligned with experimental findings, revealing two stages of fluid flux behavior during compression. In Fig. 2F and G, fluid started with a fast decrease in flux at stage 1, then a more gentle and steady flux decrease at stage 2. Interestingly, two stages were also observed in stress–strain curves and correspond to two stages in flux change: when stage 1 transitioned to stage 2 for flux, fast-increasing pressure at stage 1 transitioned to the slow-increasing pressure at stage 2 as well (Fig. 2F and G).

In stage 1, due to the sudden generation of high-pressure fluid at the center and its inability to immediately leave the gel, a large increase in chemical potential gradient is induced. As compressive strain increases with time, the high-pressure fluid in the center reaches the gel boundary and leaves without any network confinement, releasing the large difference in chemical potential gradient at stage 1 and entering stage 2. A transition region between stages 1 and 2 was observed (Fig. 2F, S3, SI). This phenomenon may result from the high-pressure fluid hitting the gel boundary and thus needing time to equilibrate and transition into stage 2. However, since stage 1 and transition stage only occur at the start of compressions and tumor growth is a long timescale process, we focused solely on the steady flux and chemical potential gradient change of stage 2 for the growth factor concentration profile and tumor growth simulation.

2.3 Péclet number governs growth factor transport and distribution in the alginate gel geometry

Next, we modeled the effects of fluid flow on growth factor transport. Mitogenic growth factors are usually produced locally within the tumor microenvironment by tumor and stromal cells. Tumor cells acquire the ability to produce mitogenic growth factors for self-sufficiency, which is one of the two major hallmarks of cancer, along with the ability to sustain angiogenesis.^38–41 Compression-induced stress influences both interstitial fluid flux and ECM permeability. Several studies have analyzed the impact of compression and permeability changes on solute convection and diffusion transport in TME, with the two processes showing different sensitivities to permeability changes based on the characteristics of the biopolymer networks and the solute properties.^9,42–44 Since solute distribution depends on the relative strengths of convective and diffusive transport, we characterize the balance with the Péclet number in the advection–diffusion solute transport partial differential equation (PDE) modeling framework.

We hypothesized that compression-dependent changes in transport properties of the biopolymer networks impact the longer timescale behavior of tumor growth due to limiting growth factor convective transport. We simulated growth factor transport in the alginate gel subjected to 30% compression over 5 h in “real-time” as opposed to simulating in a relative time space. Using the baseline run as a reference, we systematically varied the Péclet number to assess how changes in advective versus diffusive transport influence growth factor concentration fields. Fig. 3 shows the ranges of Pe number corresponding to transport in our alginate gel domain of characteristic length L = 4 mm. The diffusion coefficient (D) of a typical growth factor, e.g. vascular endothelial growth factor or VEGF, (size 35–50 kDa) is set at 100 µm² s⁻¹.⁴⁵ The temporal chemical potential gradient profiles obtained (Fig. S4, SI) from simulating 30% compression over a 5-hour period in the continuum finite element modeling framework, result in a maximum velocity (v) of ∼5.5 µm s⁻¹ (Fig. S6, SI). Given that permeability variations remain within a narrow range (Fig. S5, SI), we assume the diffusion coefficient to be spatially and temporally invariant for the purposes of this model. This gives us for the selected D–v combination. Assuming that permeability changes are minimal and the non-dimensionalized velocity profiles scale with compression under small strain rates in a fixed geometry, we utilize the same spatiotemporal profiles and scale them for simulating lower Pe cases at 0.2208, 2.208, 5, and 22.08. Note that if large compressive strains alter permeability significantly, the profile shape must be recomputed rather than merely scaled.

Fig. 3 Physiological ranges of interstitial fluid velocity (v) and solute diffusivity (D) drawn from tumor and hydrogel studies were mapped to Pe values. We thus solved for three representative transport regimes: Pe = 0.2208, 2.208 (diffusion-biased), 5, 22.08 (mixed), and 220.8 (advection-dominated) to represent the spectrum of in vivo and in vitro TME transport behavior. The underlying v and D ranges and their primary literature sources are summarized in the methods section.

We modeled the concentration dynamics of a prototypical growth factor under two scenarios: (1) growth factor saturated domain – where the tumor microenvironment is initially saturated with growth factors; initial concentration (non-dimensionalized) of the growth factor set to 1 and the boundary condition across all domain edges is set to 0 (Dirichlet boundary condition); (2) growth factor depleted domain – where the tumor microenvironment initially lacks growth factors, with growth factors supplied from the surrounding regions (e.g. blood vessels), diffusing into the tumor microenvironment; initial concentration of growth factor set to 0 and Dirichlet boundary conditions set to 1. The PDE setup, non-dimensionalization, and solution methodology have been outlined in the methods section. In the results section, we only show results from the growth factor saturated domain scenario. The growth factor depleted domain simulations and results have been discussed in the SI (S8 and S9).

In the growth factor saturated domain simulation, two transport processes operate simultaneously in radially outward directions: (i) passive diffusion driven by the emerging concentration gradient carries the growth factor radially outwards, and (ii) a radially outward interstitial fluid flux (generated by compression) carries the growth factor with it (Fig. 4A). At Pe = 0.2208, diffusion dominates, and the core concentration reduces to 15% of the maximum attainable concentration C₀ in 6 h, with gradients penetrating deep into the domain. In the mixed-transport regimes (Pe = 2.208 and 5.0), advective flux begins to counterbalance diffusion, such that the central concentration remains near 0.4C₀ and 0.7C₀, respectively, at 360 min. At Pe = 22.08, the profiles are almost flat in the core at ≈0.95C₀, with depletion confined to a narrow boundary layer. Finally, under advection-dominated conditions (Pe = 220.8), the concentration remains essentially uniform (C* ≈ 1) throughout the entire domain over the full 360 min. Together, our numerical solutions reveal that the Pe number governs both solute retention/availability and spatial distribution in the domain.

Fig. 4 Influence of Péclet number on growth-factor distribution. (A) At t = 0 min, the concentration field is uniform across the 8 mm × 8 mm domain and boundary conditions are set to 0. Black arrows indicate the imposed outward convective flux, while red arrows show the direction of diffusion. (B) Two-dimensional concentration contours at t = 0 min (left) and t = 360 min (right) illustrate how the initially uniform field evolves under a moderate flow regime (here, Pe = 5.0). By 360 min, a nearly square-shaped plateau of high concentration persists in the interior, bounded by steep boundary layers. (C) Radial profiles of at six time points (0, 60, 120, 180, 240, 360 min) for five Péclet numbers (0.2208, 2.208, 5.0, 22.08, 220.8). Low Pe leads to significant interior depletion by diffusion, whereas raising Pe progressively preserves interior growth-factor levels and confines depletion to narrow boundary regions.

2.4 Growth factor concentration gradients modulate tumor proliferation dynamics

Next, we quantified the difference in tumor growth propensities due to the spatial heterogeneity in access to growth factor levels by simulating tumor growth in an ABM. In our model, the primary source of spatial heterogeneity is manifested through a gradient in growth factor (solute) concentration due to the chemical potential gradient of water (solvent) established by poroelastic effects. Other sources of heterogeneity such as non-uniformity in the matrix or spatial effects in the stroma are not considered in our base model and can be incorporated in future extensions of our ABM. The Péclet number strongly influences the local availability of growth factors to tumor cells. In low-Pe regimes, diffusion dominates, leading to systemic depletion of growth factors and distinct spatial heterogeneity, whereas in high-Pe regimes, advective transport preserves elevated, more uniform growth factor concentrations throughout the domain. Growth factor concentration levels were used as input to the tumor ABM in PhysiCell, systemically varying from 0 to 1. This framework was used to generate a family of tumor growth curves as a function of growth factor levels. Fig. 5 summarizes the observed tumor growth dynamics over 7 days in a small 0.4 mm × 0.4 mm tumor microenvironment subsection under different growth factor availability. When the environment is replete (1) with growth factors or moderately supplied (0.64, 0.32), the cell population grows exponentially. Intermediate growth factor availability (0.16) yields slower growth, which is approximately linear. Below a threshold, net proliferation ceases: the 0.08 case remains nearly steady, while lower growth factor levels show a progressive decline in growth.

Fig. 5 Temporal profiles of tumor load (cancer cell population count) – family of curves generated from ABM simulations under varying growth factor levels. The listed growth factor levels are used as initial and boundary conditions for the PhysiCell ABM. Shaded envelopes indicate inter-replicate variability, which widens with higher growth-factor availability due to amplified stochastic proliferation. Overall, tumor load scales monotonically with growth factor availability, exhibiting net growth above 0.08 and net regression below it.

To compare tumor growth under different Péclet number conditions and spatial regions, we also conducted the ABM simulations using growth factor profiles in two distinct regions in our computational grid (region A at r = 0 mm and region B at r = 3.6 mm), under different Péclet number conditions. The advection–diffusion PDE is solved on an 8 mm × 8 mm domain, while an ABM probes two 0.4 × 0.4 mm grids (regions A and B) (Fig. 6A). The growth factor availability at regions A and B under different Pe values has been summarized in a heatmap (Fig. 6B). The observed tumor load (total number of cancer cells) after a 7-day tumor growth ABM simulation under the corresponding growth factor availabilities has been summarized in a bar plot (Fig. 6C). The cell counts remain low when Pe = 0.2208 and increase two to three-fold at higher Pe values. Under diffusion-dominated conditions (Pe = 0.2208, 2.208), mean tumor load at t = 7 days for region A was significantly higher than region B (p < 0.0001). Thus, region A tumors grew significantly more than region B tumors. However, at higher Péclet numbers, no significant difference in mean tumor load was observed between regions A and B.

Fig. 6 Coupling growth factor transport and tumor growth results across scales (A) setup of the computational study – PDE domain on the 8 mm × 8 mm domain and tumor growth ABM domain on 0.4 mm × 0.4 mm regions of the domain. (B) Heatmap – summarizing the region A and region B growth factor (GF) concentrations (at t = 6 h). Spatial heterogeneity in growth factor availability between regions A and B reduces as Pe increases. Under high Pe conditions, growth factor levels are both elevated and more uniform across the domain. (C) Tumor load (total number of cancer cells) at t = 7 days/10 080 min at regions A and B, under different Pe. Tumor load bars show the μ ± σ from ten independent ABM stochastic runs per condition. Welch's two-sample t-test was applied to each region A–B pair; a single asterisk indicates a significant difference with p < 0.01 and three asterisks indicate p < 0.0001.

3 Discussion

By measuring fluxes under stepwise compression, we demonstrated that mechanical relaxation from fluid flow dominates in biopolymer networks due to a decrease in permeability with each increment of compression. Previous studies showed interstitial fluid flow regulates cell migration⁸ and solute transport,⁹ while viscoelasticity regulates cell invasion² and proliferation via short length-scale mechanotransduction.^1,3 Our work suggests that fluid flow (poroelasticity) and solute transport dominate cell behavior in larger deformations (>10%) of biopolymer networks.

We identified a non-linear behavior of fluid flux during compression with a customized rheometer and finite element modeling. Both experiments and simulation show that fluid experiences a fast-changing stage 1, and then a slower and steadily changing stage 2 (Fig. 2D and G), which is interestingly reflected in the corresponding two stages observed in stress–strain curves (Fig. 2F and G). In stage 1, due to the sudden generation of high-pressure fluid at the center and its inability to immediately leave the gel, a fast and large change in chemical potential gradient is induced. As strain increases, the high-pressure fluid in the center reaches the gel boundary and leaves without any network confinement, releasing the large difference in chemical potential and entering stage 2. Interestingly, a transition region (Fig. 2F, S3, SI) may result from the high-pressure fluid impacting the gel boundary, which requires time to equilibrate to stage 2. However, since stage 1 and the transition stage only occur at the start of compressions and tumor growth is a long timescale process, we focused solely on the steady flux and chemical potential gradient change at stage 2 for the following concentration profile and tumor growth simulation. A previous study identified the non-linear stress–strain relationship to fluid flow based on experiments and a theoretical model.⁴⁶ Our work expanded on this finding to show a non-linear flux behavior during compression by experiments and FEM simulation.

In this study, we coupled poroelastic mechanics to an advection–diffusion framework that explicitly computes compression-induced interstitial fluid flow and its effect on growth factor (solute) gradients. We studied this system under different Péclet number conditions and observed the spatial variations in the solute distribution. In the growth factor saturated environment, higher Péclet number conditions help retain more growth factors and promote more spatially uniform distribution of growth factors as compared to low Péclet number conditions (Fig. 4C). Our model is consistent with previous findings that growth-induced solid stress and interstitial fluid pressure can influence tumor growth by suppressing proliferation and collapsing vessels,^19,47,48 and that these forces drive vascular remodeling and heterogeneity in hypoxia.^49–51

We employed an ABM to simulate and compare tumor growth under different Péclet number conditions. Our results show that Péclet number significantly impacts tumor progression. In a growth factor-saturated tumor microenvironment, high Péclet number (higher compression) promotes tumor growth, while low Péclet number results in slower progression or even regression of tumors (Fig. 6C). Under high Péclet number conditions, growth factors are more evenly distributed, resulting in similar local tumor growth rates across different regions of the tumor microenvironment (as shown by Fig. 6C high Pe cases, where no significant difference in mean tumor load was observed between regions A and B). In contrast, under low Péclet number conditions, steep solute concentration gradients emerge, resulting in heterogeneous growth factor availability and regionally variable tumor growth (as shown by Fig. 6C low Pe cases, where mean tumor load was found to be significantly different between regions A and B). Liedekerke et al. developed an ABM for tumor cell growth by calibrating model parameters to mouse colon carcinoma compression experiments, using a strain-dependent Hill function for cell cycle progression.⁵² However, these previous studies do not account for compression-induced mitogen transport. Our results showed transport properties alone can contribute to region-specific differences in local tumor growth due to spatial differences in solute availability within the TME. We expect compression-induced changes in transport may also affect immune responses through the lymphatic drainage of tumor antigens and cells, as well as immune infiltration and activation through the availability of chemokines and cytokines. Our results are focused on isolating poroelastic transport from other biological variables to clarify its role in substrate distribution and tumor growth by using poroelastic-driven advection as a standalone regulator of Péclet number and consequently, growth factor distribution. Future work will validate these findings using tumor spheroids encapsulated in the tunable biopolymer networks under intermittent compression.

4 Conclusions

We have demonstrated that poroelastic fluid flow dominates alginate biopolymer network behavior under larger compression greater than 10%. We showed experimentally and computationally that interstitial fluid transport contributes to stress relaxation. We observed a nonlinear two-stage fluid flux behavior under compression governed by evolving chemical potential gradients and permeability that mirrors corresponding changes in stress–strain behavior. Advection–diffusion simulations of growth factor transport showed tumor growth was strongly influenced by the Péclet number. In growth factor-saturated conditions, high Péclet number promoted proliferation, while slower growth or regression was observed with low Péclet number. Further, spatial heterogeneity in tumor proliferation was differentially controlled by Péclet number, with higher tumor load accumulating in regions with higher growth factor availability. Together, these results suggest that compressive stress exerted on biopolymer networks can impact tumor growth by controlling fluid-driven growth factor transport.

5 Materials and methods

5.1 Overall experimental and computational workflow

Flux and permeability were measured by compressing an alginate gel in a custom setup to model compression of the surrounding ECM during expansive tumor growth (Fig. 7). Continuum-mechanics modeling and finite-element analysis were then used to map different compression levels to outward water flux, based on the temporal profiles of the chemical potential gradient of water. The hydrogel was modeled as a hyperelastic solid network (using a Flory–Rehner free energy density function) saturated with an incompressible fluid obeying Darcy's law. This model explicitly accounts for large-deformation non-linearities, where permeability (K) evolves as a function of the local solid volume fraction (J) according to the Carman–Kozeny relation. The FE analysis outputs the transient chemical potential gradients (Δμ) and fluid velocities (v) generated by compressive strain. The fluid velocity fields derived from the FEA were mapped onto a 2D grid to solve the advection–diffusion equation (ADE) for a representative growth factor. We utilized the Péclet number (Pe = vL/D) to characterize the transport regime, systematically varying Pe to represent physiological ranges found in the tumor microenvironment (TME). This step quantifies how compression-driven convection alters the spatial concentration profile (C) of growth factors relative to passive diffusion. The spatially heterogeneous growth factor concentration profiles calculated in the PDE step were used as static environmental inputs for an ABM utilizing the PhysiCell platform. Tumor cells were modeled as discrete agents within a 0.4 mm × 0.4 mm subdomain of the larger transport grid. Cellular proliferation was stochastically determined by local growth factor availability using a Hill function saturation curve. This framework utilizes a one-way sequential coupling (mechanics → transport → biology). While tumor growth in vivo generates solid stress that further deforms the ECM (closing the loop), this study isolates the specific contribution of externally applied or accumulated compressive stress on transport phenotypes to determine how the resulting fluid flux regulates growth factor availability and subsequent proliferation dynamics. Below, we describe each method in detail:

Fig. 7 Schematic of the multiscale theoretical and computational framework of tumor growth in poroelastic biopolymer networks.

5.2 Alginate functionalization

Ultra-pure sodium alginate with very low viscosity (ProNova UP VLVG, NovaMatrix) was used, which we will be referred as VLVG in the remainder of the text. VLVG has an approximate molecular weight of 32 kDa and a guluronate-to-mannuronic monomer ratio of ≥1.5. VLVG is covalently modified by conjugating either 5-norbornene-2-methylamine (Nb, TCI America) or (4-(6-methyl-1,2,4,5-tetrazin-3-yl)phenyl)methanamine hydrochloride (Tz, Karebay Biochem), generating Nb-VLVG and Tz-VLVG, respectively.^11,28 A theoretical degree of substitution (DS) of 25% was targeted. Initially, VLVG was dissolved at a 1% w/v concentration in a 2-(N-morpholino)ethanesulfonic acid (MES) buffer with 0.3 M NaCl at pH 6.5, adjusted using NaOH. Next, N-(3-dimethylaminopropyl)-N′-ethylcarbodiimide hydrochloride (EDC) and N-hydroxysuccinimide (NHS) were added at a 2.5-fold molar excess relative to alginate carboxylic groups and continuously stirred at 4 °C. Nb or Tz was then added dropwise at a concentration of 4.7 mM per gram of VLVG. EDC, NHS, Nb, and Tz were divided into four aliquots, which were introduced at four intervals, each spanning two hours. The reaction continued for 18 hours at 4 °C with stirring. After completion, the solution underwent four rounds of ultracentrifugation at 14 000 rpm for 15 minutes each, followed by filtration through a 0.22 µm membrane. The product was then dialyzed using a 3.5 kDa molecular weight cutoff membrane (Spectra Por 3, Spectrum Labs) against a stepwise NaCl gradient of 150, 100, 50, and 0 mM in Milli-Q water over four days, with the buffer replaced twice per concentration each day. The final product was lyophilized and stored at −80 °C.

5.3 Fabrication of alginate hydrogels

VLVG alginate was prepared as previously described.^11,28 A buffered salt solution was prepared by adding 20 mM 50× N-2-hydroxyethylpiperazine-N-2-ethane sulfonic acid (HEPES, Gibco) into 1× Hank's balanced salt solution (HBSS, Gibco) and adjusting the pH to 7.5 with 10 M NaOH. Stock solutions of unmodified alginate (VLVG) and amine-functionalized alginates (Nb-VLVG and Tz-VLVG) were prepared at a concentration of 5 wt% in HBSS/HEPES buffer. Precipitated calcium carbonate (CaCO₃; Specialty Minerals) was dispersed in Milli-Q water at a 10 wt% concentration and ultrasonicated at 50% amplitude for 15 seconds to form a slurry. Freshly prepared D-(+)-gluconic acid δ-lactone (GDL; Sigma-Aldrich) was dissolved at 0.4 g mL⁻¹ in HBSS/HEPES buffer immediately before use. For purely ionically crosslinked hydrogels, HBSS–HEPES buffer, CaCO₃ with a concentration of 30 mM, VLVG with a concentration of 1.5 wt%, and GDL solution (CaCO₃ : GDL molar ratio of 0.3) were sequentially added and fully stirred on ice. For click-modified alginates (Nb-VLVG and Tz-VLVG), the HBSS/HEPES buffer, CaCO₃, Nb-VLVG stock solutions, and GDL solution with the same concentrations as unmodified VLVG alginate, were added sequentially with thorough mixing before introducing the Tz-VLVG stock solution. The click-modified alginate was made with an Nb-VLVG : Tz-VLVG ratio of 1 : 1, where covalent crosslinks were established via Nb–Tz click chemistry and ionic crosslinking was achieved with Ca²⁺.

5.4 Modulus and relaxation characterization with stepwise compression

Hydrogel solution prepared as described above was loaded in a customized acrylic mold and polymerized for 5 hours at room temperature to form a hydrogel disc (D × H = 8 mm × 2 mm). The hydrogel disc was preconditioned for tests by soaking in Milli-Q water for 24 hours. Gel disc was transferred to a stress-controlled rheometer (DHR-30, TA instrument) and confined between two impermeable plates (top: 20 mm steel geometry; bottom: immersion cup) by lowering the top plate in contact with the hydrogel surface. Then, Milli-Q water was added to submerge the hydrogel disc and conditioned until reaching a stable axial force. Stepwise compression was composed of 10% compressions followed by 2500 s of stress relaxations, ending after 8 steps and reaching 80% strain total compression. Elastic modulus was determined by performing a linear fit on the stress–strain curve for each 10% compression. Relaxation time constant for stress relaxation was characterized by fitting the stress–relaxation curve to an exponential decay function based on the generalized Maxwell model.

5.5 Flux and permeability measurement with customized compression setup

The preconditioned hydrogel disc was transferred to a stress-controlled rheometer (DHR-20, TA Instruments) and positioned between two impermeable plates (top: 20 mm steel geometry; bottom: glass sheet placed above a camera). Once the geometry fully contacted the top surface of the gel, a Moticam S20 high-resolution camera was used to capture the stepwise compression described above, ending with 40% strain compression. The recorded video was processed through ImageJ to quantify the cross-sectional area of hydrogel disc over time. Assuming V(t) = πr(t)²h(t), where πr(t)² is the cross-sectional area of the gel disc and h(t) is the gel thickness at time t, hydrogel volume, V(t), was calculated with πr² obtained from ImageJ analysis and h(t) from the rheometer. The volumetric flux, q(t), is a measure of flow rate per unit area and is mathematically expressed as with unit m s⁻¹. Volumetric flow rate, , denotes the volumetric flow rate in unit m³ s⁻¹. The fluid flow area, A(t) = 2πr(t)·h(t), denotes the effective area available for fluid outflow, where we assume gel with cylindrical geometry based on experimental observations. Permeability is an intrinsic parameter and characterizes porous material's ability to allow fluid to flow through, which can be characterized by Darcy constant, K, with unit m². K is calculated from Darcy's law with viscosity, η, of fluid traveling over distance L, where a chemical potential gradient Δμ will cause fluid flow. A Python script was used to synchronize and process both the stress data from the rheometer and the video-flux data, generating a curve of flux vs. normal stress (Fig. 1E). Considering the equation for the Darcy permeability, a slope, , is obtained to calculate permeability K by fitting a linear function to the curve. By taking L ≈ r(t), we can approximate the chemical potential gradient as , where σ is the average normal Cauchy stress on the top plate. We did so by following^53,54 who give a simple expression for fluid flow through compressed open cell foams in terms of the measured compressive stress and dimensions of the foam. We fitted the last linear part of the curve at the low flux region for each compression-relaxation step, because at the later phase of relaxation step the chemical potential gradient in the radial direction can be approximated as linear in agreement with our assumption above.³⁷

5.6 Continuum modeling and finite element analysis

5.6.1 Kinematics. We focus here on the poro-elastic modeling of the alginate gel compression experiments. The gel is modeled as a continuum. Each material point in this continuum consists of a mixture of solid and fluid molecules and is considered homogeneous and isotropic in the initial undeformed configuration. Upon (large) deformation, the solid network can become anisotropic with different tangent moduli in different directions, governed by the local deformation field. Letting x(X,t) and X denote the current and initial positions of a specific material point, respectively, we can define the deformation gradient tensor and its determinant as:

The initial (reference) configuration is taken to be the swollen state that the alginate network forms with water at equilibrium. We assume that the alginate material and the interstitial liquid are both incompressible and that volume change of the gel is caused purely by the movement of liquid. Let C(X,t) denote the volume fraction of liquid in the reference state, then:

C(X,t) = J − ϕ^ref_s, (2)

where ϕ^ref_s denotes the volume fraction of the polymer in the reference configuration. Now, if we consider a volume V with a normal N(X) and surface S(X) in the reference configuration, then mass balance in the absence of any liquid source can be written as:where Q denotes the liquid flux per unit reference area and ρ denotes the density of the liquid. Utilizing the divergence theorem, we can rewrite this equation in localized form as:

If we consider the loading to be quasi-static, which means that the load is applied slowly, we can write our momentum balance equation expressed in the reference configuration as:

Div P = 0, (5)

where P is the first Piola–Kirchhoff stress tensor. This type of chemo-elastic formulation has been adopted in previous studies^55,56 and it is equivalent to the theory of porous media (TPM), including under large deformations. TPM explicitly accounts for the relative motion between solid and fluid by formulating separate equilibrium and mass balance equations for the solid and fluid constituents, which allow a material point in the Eulerian space to be occupied by both solid and fluid particles. It also incorporates body forces exerted on both the solid and fluid phases by the relative motion. It has been shown by Stracuzzi et al.⁵⁷ that for incompressible solid and fluid constituents, the more detailed TPM framework is mathematically equivalent to the chemo-elastic formulation employed here. Therefore, the relative motion between the solid and fluid has been implicitly incorporated in the FEM implementation.

5.6.2 Constitutive laws. Experimental results indicate that the viscous effects are negligible and less important for long time scales and large compressive strains. Consequently, dissipation is caused only by the flow of liquid through pores and we consider the network to be elastic. Based on Flory–Rehner theory,⁵⁸ the total Helmholtz free energy per reference volume of the polymer network and liquid mixture can be written as a summation of the free energy density due to the deformation of the network Ψ_net(F), and the free energy density of interaction of liquid and solid Ψ_int(C):

Ψ(F,C) = Ψ_net(F) + Ψ_int(C). (6)

Utilizing Legendre transformation, we define a new free energy density as:

[capital Psi, Greek, circumflex] = Ψ(F,C) − μC. (7)

Here is the chemical potential of the interstitial liquid.⁵⁶ Based on standard continuum mechanical derivations that use the dissipation inequality, it can be shown that:

The free energy Ψ_net(F) and Ψ_int(C) used in this work is based on:^57,59,60

Notice that C in eqn (9) is replaced by J according to eqn (2). Here G₁ and G₂ are two constants. G₁ is the small strain shear modulus of the dry polymer, while G₂ is related to the shear modulus at large strain. A two-term polynomial hyper-elastic model is used to characterize the strain-stiffening behaviors of the polymer under compression. A ln J term is included in the expression for the free energy density following standard texts in nonlinear elasticity (e.g., Ogden,⁶¹ Holzapfel⁶²) and biophysics (e.g., Boal⁶³). Ī₁ is the modified invariant of the left Cauchy–Green tensor B. The isochoric part of I₁ is chosen to account for the incompressibility of both liquid and solid phases. π₀ is the initial osmotic pressure, and the value of α₁ is chosen such that the reference state is stress free. β₁ is a power law coefficient, which is slightly larger than 1 in the current formulation. Utilizing eqn (9), the Cauchy stress can be calculated as:

In order to set the value of α₁, we require σ(I,0) = 0 at t = 0, so that:

G₁ and G₂ are fitted to experimental stress–strain data from compression experiments as shown in the SI.

5.6.3 Liquid flux. In this work, Darcy's law is used for the constitutive behavior of the liquid flux. Further, we assume the material is isotropic so that this relation can be written as:where η is the viscosity of the liquid and K is the permeability. K is a decreasing function of the current solid phase volume fraction ϕ_s which evolves with local deformation. We use the Carman–Kozeny expression for the permeability:⁶⁴where c₀ is a constant less than 1 and σ is the specific surface area. Here we assume that the alginate chains form tubes with circular cross section, which results in c₀ = 0.5 and . Note that the current solid phase volume fraction is related to the reference solid volume fraction by , which depends on the local deformation. Therefore, the permeability is not a constant in this formulation. Since it depends explicitly on J, which varies in space and time under non-uniform deformation, the permeability is itself a field variable and is thus a function of position and time, K = K(x,t).

To solve the full boundary value problem we use (a) impermeable boundary conditions on the top and bottom plates and permeable boundary conditions on the lateral surfaces, (b) displacement boundary condition on the top surface, fixed boundary condition on the bottom surface, and traction free boundary conditions on the lateral surface. The model is implemented in commercial finite element software ABAQUS Standard in the SOILS module with constitutive laws input through a user subroutine. C3D8P brick (Taylor-Hood) elements are used for the discretization. The parameters used for finite element calculations are presented in the SI.

5.7 Growth factor transport

5.7.1 Advection–diffusion PDE and Péclet number. Compression of the poroelastic gel induces a range of outward fluid fluxes under varying levels of stress and strain rates. This flux is characterized using the chemical potential gradient and permeability profiles obtained from the constitutive finite element modeling. FEM framework was used to simulate 30% compression of the gel over a 5-hour period, to represent the growth factor transport dynamics at a time scale comparable to the cellular growth time scale (∼a few hours). Chemical potential as a function of radial position was recorded every 10 minutes during the first hour and then at hourly intervals for the remainder of the simulation (i.e., 2, 3, 4, and 5 hours). A logistic function was fitted to the chemical potential gradient data to characterize the radial profiles (eqn (14b)). Λ represents the maximum attainable value of the chemical potential gradient in the domain (unit Pa), k is the steepness coefficient (unit m⁻¹) and r₀ is the inflection point (unit m). Permeability was also extracted as a function of time from the FEM simulations, computed using eqn (13).

The advection term u(r,t) associated with outward fluid flux was estimated using eqn (14a) and (14b).

The fluid flux induced by compression influences the concentration profiles of diffusible substrates (e.g., nutrients, growth factors, metabolites) within the tumor microenvironment. To quantify this effect, we analyze growth factor concentration profiles within an advection–diffusion PDE framework, incorporating fluid flux as a governing transport mechanism in addition to diffusion.

The advection–diffusion equation was solved for an 8 mm × 8 mm two-dimensional square grid, designed to mimic the compression experiments on the alginate gel in silico. This grid represents a tumor microenvironment assumed to be isotropic and symmetric. In this PDE setup, we are not accounting for growth factor depletion due to consumption by tumor cells. We also do not account for heterogeneous perfusion, transvascular exchange or extravascular binding in the 2D tumor microenvironment. These terms can be accounted for using an additional reaction term in the PDE. Baxter and Jain provided an extensive theoretical framework to model advection–diffusion-reaction PDEs for macromolecules solute transport in the tumor microenvironment.^23,65–67 This can be adapted for growth factors.

The advection–diffusion equation in two-dimensional Cartesian coordinates describes the transport of a scalar concentration field, C(x,y,t), under the combined effects of advection and diffusion (eqn (15)). The spatial coordinates x and y represent length dimensions, while t denotes time. The velocity components in these directions are u_x and u_y, respectively. The diffusion coefficients are assumed constant throughout the domain, with D_x = D_y = D.

To nondimensionalize the equation, characteristic scales were introduced based on the advection time scale: characteristic length L, velocity U, concentration C₀, and advection time . The resulting non-dimensional variables are defined as , , , and , with velocity components and .

The dimensionless Péclet number, the ratio between rates of advection to diffusion, given by Pe = UL/D, characterizes the relative influence of advection to diffusion on the characteristic length scale L. The convection–diffusion equation was non-dimensionalized using the advection time scale(eqn (16)).

The characteristic length L is defined as 4 mm, consistent with the radial dimensions of the experimental setup used to study poroelastic effects under different conditions. The characteristic velocity U corresponds to the maximum radial velocity observed from the FEM simulations. The characteristic concentration C₀ is set to a nominal growth factor concentration of 10⁻⁷ M, reflecting typical concentration of growth factors in tumor microenvironment.

Diffusion coefficient of a growth factor (or other macromolecular solute in a tumor tissue or engineered gels) presented in literature typically ranges from 10⁻⁸ to 10⁻¹² m² s⁻¹, with variations influenced by the size, charge, porosity, and interaction with extracellular matrix components. (EGF in rat brain extra cellular space:⁶⁸ 10⁻¹¹ m² s⁻¹, VEGF in type 1 collagen hydrogel:⁴⁵ 10⁻¹⁰ m² s⁻¹, IGF in fibrin gels:⁶⁹ 10⁻⁸ m² s⁻¹). Interstitial fluid velocities in normal tissues are in the range 0.1–1 µm s⁻¹, whereas velocities in tumor tissues have been noted to be as high as 55 µm s⁻¹ in some cases.⁷⁰

5.7.2 Numerical scheme for solving the PDE. The computational domain was defined over x*, y* ∈ [0, 1] with symmetrical Dirichlet boundary conditions for growth factor concentration on all 4 edges of the grid. The initial condition was set as a uniform concentration of C* across the entire domain.

The advection–diffusion equation was discretized using the method of lines (MOL), with spatial derivatives approximated via finite-difference method on a uniform grid with spacing Δx* = Δy* = 0.05. The advection term was discretized using an upwind scheme (order 1) for numerical stability. The advection term is updated at discrete time intervals throughout the PDE numerical scheme to account for the change in chemical potential gradients and mobility with time. This update is done every 10 min for the first hour and then, every hour, until the end of the simulation. Adjusting the frequency of updates to the advection term influences both the accuracy of the numerical solution and the computational efficiency. By updating more frequently during the initial phases of the simulation, we capture initially faster transient behaviors with greater precision.

Time integration was performed using the implicit second-order TRBDF2 (trapezoidal backward differentiation formula 2) method, which is suitable for stiff problems.

The numerical solution was implemented in Julia 1.11.3, utilizing ModelingToolkit.jl for symbolic modeling and equation generation, MOLFiniteDifference.jl for spatial discretization via the method of lines approach, and DifferentialEquations.jl for numerical solving of the advection–diffusion equation.

5.8 Agent-based modeling for tumor growth under different growth factor availability

To investigate the impact of compression-induced growth factor gradients on tumor progression, an ABM was developed using PhysiCell (Version 1.14.0),⁷¹ an open-source, physics-based cell simulator for multicellular systems. PhysiCell incorporates BioFVM,⁷² a parallelized diffusive transport solver for 2D and 3D biological simulations, to simulate growth factor transport dynamics within the tumor microenvironment. Growth factor concentrations ranging from 0 to 10⁻⁷ M were used to parameterize the tumor microenvironment by specifying the initial and boundary conditions for substrate fields within the PhysiCell framework.

All ABM simulations were executed on the Bridges-2⁷³ cluster at the Pittsburgh Supercomputing Center. To capture stochastic variability, ten independent ABM simulations were run per condition. Results are summarized as mean ± standard deviation. Between-condition differences in the means were evaluated with Welch's two-sample t-test (α = 0.05, two-tailed), which does not assume equal variances. Snapshots of the simulation were saved as scalable vector graphics (SVGs) at specified time intervals for visualization.

The ABM was simulated on a 0.4 mm × 0.4 mm square domain (a small sub-section of the 8 × 8 mm alginate gel) with an initial tumor cell population of 200 cells. The simulation was run for 7 days (10 080 minutes), with tumor cells proliferating in response to local growth factor concentrations. Each cell sampled its local microenvironment and determined its growth rate based on a sigmoidal relationship between growth factor concentration and cellular proliferation (S7, SI).

Growth factor-receptor interactions typically exhibit saturable binding kinetics. This is depicted by a sigmoidal curve, where low concentrations result in minimal proliferation due to insufficient receptor activation, while increasing concentrations enhance proliferation rates up to a critical threshold. Beyond this point, receptor saturation limits further increase in proliferation despite higher extracellular growth factor levels.^74–76 In our ABM, growth factor level – cellular growth relationship was defined by a saturation growth rate of 10⁻⁴, a Hill coefficient of 2, and a half-max concentration of 10⁻⁹ M. These values were selected as an approximation for a generic growth factor but can be adjusted depending on the specific growth factor or available experimental data. Additional parameters used in the ABM were adapted from previously published PhysiCell simulations⁷¹ and are listed in the supplementary section (Table S1).

Previous computational studies have built agent-based and hybrid models to study the impact of vascular structure and function on emergent cell population behaviors.⁷⁷ These frameworks have also been used to understand how angiogenesis drives tumor progression and to model how vascular features influence drug delivery.^78–80 Since our ABM is designed as a computational platform to emulate cancer cell growth within an alginate gel geometry, we do not incorporate certain physiologically relevant factors, such as ECM heterogeneity, vascular heterogeneity, and nutrient competition, into the current model. Only heterogeneity in growth factor gradients generated by compression-driven interstitial flow is represented. Individual cancer cells are assumed to be phenotypically identical, and proliferation is governed solely by local growth factor concentration. Direct effects of solid stress, hydrostatic pressure, or matrix stiffness on cell cycle progression are not modeled. These simplifications help to isolate the effect of mechanical compression on growth factor transport and demonstrate its impact on tumor proliferation.

Future extensions of these models could couple multiscale mechanotransduction models, ECM remodeling, nutrient dynamics, and vascular heterogeneity to capture the joint influence of stress, stiffness, and biochemical cues on tumor growth.

Author contributions

Authors ZL, YR, SK, PKP, RR, PAJ, and KHV contributed to the conceptualization, data curation, writing, and review and editing of this manuscript. ZL, YR, SK, PM, and JK contributed to the investigation, data analysis, and visualization. All authors reviewed and approved final version of the manuscript.

Conflicts of interest

The authors have no conflicts to disclose.

Data availability

The data that support the findings will be available in Dryad repository at DOI: https://doi.org/10.5061/dryad.2z34tmq1h after publication.

Supplementary information (SI): supplementary figures (S1–S9) and Tables S1 and S2. See DOI: https://doi.org/10.1039/d6sm00060f.

Acknowledgements

This study was supported by the REgeneration and ReSToration of Functions in ORal HEalth (RESTORE) Prize funded by the Center for Innovation & Precision Dentistry (CiPD) and the Department of Oral and Maxillofacial Surgery's Schoenleber Fund, at the University of Pennsylvania. KHV received support from the National Institutes of Health for this work, NIH/NIDCR R00DE030084 and L70DE035343. This work was partially supported by NSF through the University of Pennsylvania Materials Research Science and Engineering Center (MRSEC) (DMR-2309043) through a seed grant to KHV and PKP. PKP and YR acknowledge NSF grant DMR 2212162. PKP acknowledges partial support from NIH grant R01 HL148227. RR and SK acknowledge funding by the National Institutes of Health Grants GM136259 (also PAJ), EB017753 and CA250044. We acknowledge the use of the Bridges-2 system, a resource of the Pittsburgh Supercomputing Center (PSC), part of the ACCESS program under grant MCB200101. The graphical abstract was edited using NotebookLM and Gemini Pro.

References

1. K. H. Vining and D. J. Mooney, Mechanical forces direct stem cell behaviour in development and regeneration, Nat. Rev. Mol. Cell Biol., 2017, 18, 728–742.
2. A. Elosegui-Artola, A. Gupta, A. J. Najibi, B. R. Seo, R. Garry, C. M. Tringides, I. de Lázaro, M. Darnell, W. Gu, Q. Zhou, D. A. Weitz, L. Mahadevan and D. J. Mooney, Matrix viscoelasticity controls spatiotemporal tissue organization, Nat. Mater., 2023, 22, 117–127.
3. H.-P. Lee, F. Alisafaei, K. Adebawale, J. Chang, V. B. Shenoy and O. Chaudhuri, The nuclear piston activates mechanosensitive ion channels to generate cell migration paths in confining microenvironments, Sci. Adv., 2021, 7, eabd4058.
4. K. M. Wisdom, K. Adebowale, J. Chang, J. Y. Lee, S. Nam, R. Desai, N. S. Rossen, M. Rafat, R. B. West, L. Hodgson and O. Chaudhuri, Matrix mechanical plasticity regulates cancer cell migration through confining microenvironments, Nat. Commun., 2018, 9, 4144.
5. S. Nam, V. K. Gupta, H.-P. Lee, J. Y. Lee, K. M. Wisdom, S. Varma, E. M. Flaum, C. Davis, R. B. West and O. Chaudhuri, Cell cycle progression in confining microenvironments is regulated by a growth-responsive TRPV4-PI3K/Akt-p27Kip1 signaling axis, Sci. Adv., 2019, 5, eaaw6171.
6. K. Parihar, J. Nukpezah, D. V. Iwamoto, K. Cruz, F. J. Byfield, L. Chin, M. E. Murray, M. G. Mendez, A. S. van Oosten, A. Herrmann, E. E. Charrier, P. A. Galie, M. Donlick, T. Lee, P. A. Janmey and R. Radhakrishnan, Tissue-dependent mechanosensing by cells derived from human tumors, npj Biol. Phys. Mech., 2025, 2, 19.
7. K. M. Stroka, H. Jiang, S.-H. Chen, Z. Tong, D. Wirtz, S. X. Sun and K. Konstantopoulos, Water Permeation Drives Tumor Cell Migration in Confined Microenvironments, Cell, 2014, 157, 611–623.
8. W. J. Polacheck, J. L. Charest and R. D. Kamm, Interstitial flow influences direction of tumor cell migration through competing mechanisms, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 11115–11120.
9. P. A. Netti, D. A. Berk, M. A. Swartz, A. J. Grodzinsky and R. K. Jain, Role of Extracellular Matrix Assembly in Interstitial Transport in Solid Tumors1, Cancer Res., 2000, 60, 2497–2503.
10. Q.-M. Wang, A. C. Mohan, M. L. Oyen and X.-H. Zhao, Separating viscoelasticity and poroelasticity of gels with different length and time scales, Acta Mech. Sin., 2014, 30, 20–27.
11. K. H. Vining, A. Stafford and D. J. Mooney, Sequential modes of crosslinking tune viscoelasticity of cell-instructive hydrogels, Biomaterials, 2019, 188, 187–197.
12. O. Chaudhuri, L. Gu, D. Klumpers, M. Darnell, S. A. Bencherif, J. C. Weaver, N. Huebsch, H.-P. Lee, E. Lippens, G. N. Duda and D. J. Mooney, Hydrogels with tunable stress relaxation regulate stem cell fate and activity, Nat. Mater., 2016, 15, 326–334.
13. A. C. Shieh, Biomechanical forces shape the tumor microenvironment, Ann. Biomed. Eng., 2011, 39(5), 1379–1389.
14. G. Helmlinger, P. A. Netti, H. C. Lichtenbeld, R. J. Melder and R. K. Jain, Solid stress inhibits the growth of multicellular tumor spheroids, Nat. Biotechnol., 1997, 15(8), 778–783.
15. T. Selbekk, R. Brekken, O. Solheim, S. Lydersen, T. A. Hernes and G. Unsgaard, Tissue motion and strain in the human brain assessed by intraoperative ultrasound in glioma patients, Ultrasound Med. Biol., 2010, 36(1), 2–10.
16. R. K. Jain, J. D. Martin and T. Stylianopoulos, The Role of Mechanical Forces in Tumor Growth and Therapy, Annu. Rev. Biomed. Eng., 2014, 16, 321–346.
17. Q. Liu, Q. Luo, Y. Ju and G. Song, Role of the mechanical microenvironment in cancer development and progression, Cancer Biol. Med., 2020, 17, 282–292.
18. T. Stylianopoulos, The Solid Mechanics of Cancer and Strategies for Improved Therapy, J. Biomech. Eng., 2017, 139, 021004.
19. F. Mpekris, S. Angeli, A. P. Pirentis and T. Stylianopoulos, Stress-mediated progression of solid tumors: effect of mechanical stress on tissue oxygenation, cancer cell proliferation, and drug delivery, Biomech. Model. Mechanobiol., 2015, 14, 1391–1402.
20. J. A. Linke, L. L. Munn and R. K. Jain, Compressive stresses in cancer: characterization and implications for tumour progression and treatment, Nat. Rev. Cancer, 2024, 24, 768–791.
21. T. Stylianopoulos, J. D. Martin, V. P. Chauhan, S. R. Jain, B. Diop-Frimpong, N. Bardeesy, B. L. Smith, C. R. Ferrone, F. J. Hornicek, Y. Boucher, L. L. Munn and R. K. Jain, Causes, consequences, and remedies for growth-induced solid stress in murine and human tumors, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 15101–15108.
22. M. Soltani and P. Chen, Numerical Modeling of Fluid Flow in Solid Tumors, PLoS One, 2011, 6, e20344.
23. L. T. Baxter and R. K. Jain, Transport of fluid and macromolecules in tumors. I. Role of interstitial pressure and convection, Microvasc. Res., 1989, 37, 77–104.
24. S. R. Chary and R. K. Jain, Direct measurement of interstitial convection and diffusion of albumin in normal and neoplastic tissues by fluorescence photobleaching, Proc. Natl. Acad. Sci. U. S. A., 1989, 86, 5385–5389.
25. R. K. Jain, R. T. Tong and L. L. Munn, Effect of vascular normalization by antiangiogenic therapy on interstitial hypertension, peritumor edema, and lymphatic metastasis: insights from a mathematical model, Cancer Res., 2007, 67, 2729–2735.
26. M. R. Islam and M. L. Oyen, Load-Relaxation Characteristics of Chemical and Physical Hydrogels as Soft Tissue Mimics, Exp. Mech., 2021, 61, 939–949.
27. T. Distler, E. Schaller, P. Steinmann, A. R. Boccaccini and S. Budday, Alginate-based hydrogels show the same complex mechanical behavior as brain tissue, J. Mech. Behav. Biomed. Mater., 2020, 111, 103979.
28. K. Zhang, Z. Li, Y.-C. Chen, I.-C. Yoon, A. Graham and K. Vining, Shear and compressive stiffening of dual-cross-linked alginate hydrogels with tunable viscoelasticity, ACS Appl. Bio Mater., 2025, 8(5), 3899–3908 . PMID: 40274262.
29. M. Perepelyuk, L. Chin, X. Cao, A. V. Oosten, V. B. Shenoy, P. A. Janmey and R. G. Wells, Normal and Fibrotic Rat Livers Demonstrate Shear Strain Softening and Compression Stiffening: A Model for Soft Tissue Mechanics, PLoS One, 2016, 11, e0146588.
30. A. S. G. van Oosten, X. Chen, L. Chin, K. Cruz, A. E. Patteson, K. Pogoda, V. B. Shenoy and P. A. Janmey, Emergence of tissue-like mechanics from fibrous networks confined by close-packed cells, Nature, 2019, 573, 96–101.
31. J. A. Rowley, G. Madlambayan and D. J. Mooney, Alginate hydrogels as synthetic extracellular matrix materials, Biomaterials, 1999, 20(1), 45–53.
32. T. P. Richardson, M. C. Peters, A. B. Ennett and D. J. Mooney, Polymeric system for dual growth factor delivery, Nat. Biotechnol., 2001, 19, 1029–1034.
33. J. West, M. Robertson-Tessi and A. R. A. Anderson, Agent-based methods facilitate integrative science in cancer, Trends Cell Biol., 2023, 33, 300–311.
34. S. Kemkar, M. Tao, A. Ghosh, G. Stamatakos, N. Graf, K. Poorey, U. Balakrishnan, N. Trask and R. Radhakrishnan, Towards verifiable cancer digital twins: tissue level modeling protocol for precision medicine, Front. Physiol., 2024, 15, 1473125.
35. R. S. Stowers, S. C. Allen, K. Sanchez, C. L. Davis, N. D. Ebelt, C. Van Den Berg and L. J. Suggs, Extracellular matrix stiffening induces a malignant phenotypic transition in breast epithelial cells, Cell. Mol. Bioeng., 2017, 10(1), 114–123.
36. J. Huang, L. Zhang, D. Wan, L. Zhou, S. Zheng, S. Lin and Y. Qiao, Extracellular matrix and its therapeutic potential for cancer treatment, Signal Transduction Targeted Ther., 2021, 6(1), 153.
37. P. Mollenkopf, J. A. Kochanowski, Y. Ren, K. H. Vining, P. A. Janmey and P. K. Purohit, Poroelasticity and permeability of fibrous polymer networks under compression, Soft Matter, 2025, 21, 2400–2412.
38. D. Hanahan and R. A. Weinberg, The Hallmarks of Cancer, Cell, 100, 57–70, 2000.
39. X. Zhang, D. Nie and S. Chakrabarty, Growth factors in tumor microenvironment, Front. Biosci., 2010, 15, 151–165.
40. E. Witsch, M. Sela and Y. Yarden, Roles for Growth Factors in Cancer Progression, Physiology, 2010, 25, 85–101.
41. X. Jiang, J. Wang, X. Deng, F. Xiong, S. Zhang, Z. Gong, X. Li, K. Cao, H. Deng, Y. He, Q. Liao, B. Xiang, M. Zhou, C. Guo, Z. Zeng, G. Li, X. Li and W. Xiong, The role of microenvironment in tumor angiogenesis, J. Exp. Clin. Cancer Res., 2020, 39, 204.
42. R. C. Evans and T. M. Quinn, Dynamic Compression Augments Interstitial Transport of a Glucose-Like Solute in Articular Cartilage, Biophys. J., 2006, 91, 1541–1547.
43. J. Ferruzzi, M. Sun, A. Gkousioudi, A. Pilvar, D. Roblyer, Y. Zhang and M. H. Zaman, Compressive Remodeling Alters Fluid Transport Properties of Collagen Networks – Implications for Tumor Growth, Sci. Rep., 2019, 9, 17151.
44. S. Ramanujan, A. Pluen, T. D. McKee, E. B. Brown, Y. Boucher and R. K. Jain, Diffusion and convection in collagen gels: implications for transport in the tumor interstitium, Biophys. J., 2002, 83, 1650–1660.
45. Y. Shin, S. Han, J. S. Jeon, K. Yamamoto, I. K. Zervantonakis, R. Sudo, R. D. Kamm and S. Chung, Microfluidic assay for simultaneous culture of multiple cell types on surfaces or within hydrogels, Nat. Protoc., 2012, 7, 1247–1259.
46. M. T. J. J. M. Punter, B. E. Vos, B. M. Mulder and G. H. Koenderink, Poroelasticity of (bio)polymer networks during compression: theory and experiment, Soft Matter, 2020, 16(5), 1298–1305.
47. T. Stylianopoulos, J. D. Martin, M. Snuderl, F. Mpekris, S. R. Jain and R. K. Jain, Co-evolution of solid stress and interstitial fluid pressure in tumors during progression: Implications for vascular collapse, Cancer Res., 2013, 73, 3833–3841.
48. H. Yan, D. Ramirez-Guerrero, J. Lowengrub and M. Wu, Stress generation, relaxation and size control in confined tumor growth, PLoS Comput. Biol., 2021, 17, e1009701.
49. M. Welter and H. Rieger, Interstitial fluid flow and drug delivery in vascularized tumors: A computational model, PLoS One, 2013, 8, 1–23.
50. H. Rieger and M. Welter, Integrative models of vascular remodeling during tumor growth, Wiley Interdiscip. Rev.:Syst. Biol. Med., 2015, 7(3), 113–129.
51. M. Welter, T. Fredrich, H. Rinneberg and H. Rieger, Computational model for tumor oxygenation applied to clinical data on breast tumor hemoglobin concentrations suggests vascular dilatation and compression, PLoS One, 2016, 11, 1–42.
52. P. V. Liedekerke, J. Neitsch, T. Johann, K. Alessandri, P. Nassoy and D. Drasdo, Quantitative cell-based model predicts mechanical stress response of growing tumor spheroids over various growth conditions and cell lines, PLoS Comput. Biol., 2019, 15, e1006273.
53. L. J. Gibson, Cellular Solids, MRS Bull., 2003, 28, 270–274.
54. M. F. Ashby and L. J. Gibson, Cellular solids: structure and properties, Press Syndicate of the University of Cambridge, Cambridge, UK, 1997, pp. 175–231.
55. S. A. Chester and L. Anand, A coupled theory of fluid permeation and large deformations for elastomeric materials, J. Mech. Phys. Solids, 2010, 58(11), 1879–1906.
56. W. Hong, Z. Liu and Z. Suo, Inhomogeneous swelling of a gel in equilibrium with a solvent and mechanical load, Int. J. Solids Struct., 2009, 46(17), 3282–3289.
57. A. Stracuzzi, E. Mazza and A. E. Ehret, Chemomechanical models for soft tissues based on the reconciliation of porous media and swelling polymer theories, Z. Angew. Math. Mech., 2018, 98(12), 2135–2154.
58. P. J. Flory and J. Rehner Jr, Statistical mechanics of cross-linked polymer networks ii. swelling, J. Chem. Phys., 1943, 11(11), 521–526.
59. K. Garyfallogiannis, P. K. Purohit and J. L. Bassani, Cracks in tensile-contracting and tensile-dilating poroelastic materials, Int. J. Solids Struct., 2024, 286, 112563.
60. J. E. Bischoff, E. M. Arruda and K. Grosh, A new constitutive model for the compressibility of elastomers at finite deformations, Rubber Chem. Technol., 2001, 74(4), 541–559.
61. R. W. Ogden, Non-linear elastic deformations, Courier Corporation, 1997.
62. G. A. Holzapfel, Nonlinear solid mechanics: a continuum approach for engineering science, 2002.
63. D. H. Boal, Mechanics of the Cell, Cambridge University Press, 2012.
64. R. Schulz, N. Ray, S. Zech, A. Rupp and P. Knabner, Beyond kozeny-carman: predicting the permeability in porous media, Transp. Porous Media, 2019, 130, 487–512.
65. L. T. Baxter and R. K. Jain, Transport of fluid and macromolecules in tumors. II. Role of heterogeneous perfusion and lymphatics, Microvasc. Res., 1990, 40, 246–263.
66. L. T. Baxter and R. K. Jain, Transport of fluid and macromolecules in tumors. III. Role of binding and metabolism, Microvasc. Res., 1991, 41, 5–23.
67. L. T. Baxter and R. K. Jain, Transport of fluid and macromolecules in tumors. IV. A microscopic model of the perivascular distribution, Microvasc. Res., 1991, 41, 252–272.
68. R. G. Thorne, S. Hrabetová and C. Nicholson, Diffusion of epidermal growth factor in rat brain extracellular space measured by integrative optical imaging, J. Neurophysiol., 2004, 92, 3471–3481.
69. J. V. Nauman, P. G. Campbell, F. Lanni and J. L. Anderson, Diffusion of Insulin-Like Growth Factor-I and Ribonuclease through Fibrin Gels, Biophys. J., 2007, 92, 4444–4450.
70. G. Baronzio, G. Parmar and M. Baronzio, Overview of Methods for Overcoming Hindrance to Drug Delivery to Tumors, with Special Attention to Tumor Interstitial Fluid, Front. Oncol., 2015, 5, 165.
71. A. Ghaffarizadeh, R. Heiland, S. H. Friedman, S. M. Mumenthaler and P. Macklin, PhysiCell: An open source physics-based cell simulator for 3-D multicellular systems, PLoS Comput. Biol., 2018, 14, e1005991.
72. A. Ghaffarizadeh, S. H. Friedman and P. Macklin, BioFVM: an efficient, parallelized diffusive transport solver for 3-D biological simulations, Bioinformatics, 2016, 32, 1256–1258.
73. Bridges-2, 2020.
74. C.-H. Chen and S. C. Chen, Cell growth factor activity: New quantitative method in cell culture assay, Exp. Cell Res., 1981, 136, 43–51.
75. T. Domagala, K. Nicky, S. Fiona, J. Robert N., F. Louis, G. Detlef, L. Irit, S. Joseph, S. William, H. Geoffrey J., B. Antony W. and E. C. Nice, Stoichiometry, Kinetic and Binding Analysis of the Interaction between Epidermal Growth Factor (EGF) and the Extracellular Domain of the EGF Receptor, Growth Factors, 2000, 18, 11–29,  DOI:10.3109/08977190009003231.
76. D. J. Knauer, H. S. Wiley and D. D. Cunningham, Relationship between epidermal growth factor receptor occupancy and mitogenic response. Quantitative analysis using a steady state model system, J. Biol. Chem., 1984, 259, 5623–5631.
77. J. S. Yu and N. Bagheri, Modular microenvironment components reproduce vascular dynamics de novo in a multi-scale agent-based model, Cell Syst., 2021, 12(8), 795–809.e9.
78. T. Duswald, E. A. B. F. Lima, J. T. Oden and B. Wohlmuth, Bridging scales: A hybrid model to simulate vascular tumor growth and treatment response, Comput. Meth. Appl. Mech. Eng., 2024, 418, 116566.
79. C. M. Phillips, E. A. B. F. Lima, R. T. Woodall, A. Brock and T. E. Yankeelov, A hybrid model of tumor growth and angiogenesis: In silico experiments, PLoS One, 2020, 15(4), e0231137.
80. J. M. Benítez, L. García-Mozos, A. Santos, F. Montáns and L. Saucedo Mora, A simple agent-based model to simulate 3d tumor-induced angiogenesis considering the evolution of the hypoxic conditions of the cells, Eng. Comput., 2022, 38(5), 4115–4133.

---

Footnote

† These authors contributed equally to this work.

---