Physics-Informed Neural Networks for Osteosarcoma Tumor-Immune Dynamics

1  Introduction

Osteosarcoma constitutes the predominant malignant bone tumor in adolescent populations, with an annual global incidence of approximately 4–5 cases per million [1]. Epidemiological analysis of the SEER database spanning 1975–2017 indicates peak incidence among individuals aged 10–24 years, with significantly higher rates in males [2]. The disease presents substantial therapeutic challenges due to the heterogeneous tumor microenvironment, where malignant cells engage in dynamic interactions with immune components and extracellular matrix, generating behaviors resistant to simple characterization. Existing mathematical frameworks based on ordinary differential equations fail to capture spatial variability and mechanical influences on tumor growth and dissemination [3,4]. This deficiency has motivated the development of multiphysics models employing partial differential equations to provide more realistic representations of tumor-immune interactions and long-term evolutionary trajectories.

Physics-Informed Neural Networks (PINNs) have emerged as a methodology addressing limitations of both traditional machine learning and conventional numerical methods [5,6]. Whereas data-driven approaches require substantial experimental datasets often unavailable in biological contexts, PINNs incorporate physical laws directly into the learning objective [7]. In the context of the present coupled tumor-immune-mechanical problem, three specific advantages motivate the PINN choice over conventional finite element method (FEM) or finite difference method (FDM) solvers. First, the mesh-free character circumvents the generation and adaptive refinement of meshes over the circular domain, which would otherwise be required to capture localized gradients near the tumor core. Second, automatic differentiation delivers exact spatial and temporal derivatives entering the residual functionals, avoiding discretization errors associated with stiff terms such as the rapid IFN-γ decay. Third, once trained, the network provides a continuous surrogate Ω×𝒯→R8 that can be queried at arbitrary points without additional solver invocations, which simplifies downstream sensitivity and inverse studies. Recent work has extended PINN formulations to coupled mechano-biological problems, including lymph-node remodeling under inflammation and homeostasis [8] and solid-mechanics problems with complex geometries [9], further supporting the applicability of this class of methods to multiphysics biomedical models.

A complementary motivation arises from the biology of macrophage polarization. In four-variable tumor-immune models, macrophage activation is typically expressed as a direct response to cancer cell density [5]. In reality, polarization toward the anti-tumoral M1 phenotype is primarily cytokine-mediated, with IFN-γ acting as the dominant activating signal [8,10]. The present work introduces IFN-γ as an explicit fifth variable so that T-cells activated by tumor cells produce IFN-γ, which in turn polarizes macrophages and directly inhibits tumor proliferation. This single modification makes the activation cascade fully explicit.

The aims of this paper are therefore three-fold: (i) to formulate an extended five-variable tumor-immune-mechanical model that embeds the IFN-γ-mediated activation cascade within a Biot poroelastic framework; (ii) to develop a dedicated PINN training protocol capable of handling the temporal stiffness induced by the rapid IFN-γ decay through curriculum learning and adaptive loss weighting; and (iii) to systematically investigate the influence of macrophage spatial distribution on tumor dynamics across four initial configurations representative of distinct immune-infiltration patterns.

The main findings can be summarized as follows. The proposed framework achieves stable convergence across all four scenarios, with total loss values in the range [0.056,0.158] and mechanical residuals below 3.2×10−5. The boundary-concentrated macrophage configuration yields the lowest biological loss and the strongest early tumor suppression, while the cluster configuration exhibits the highest residual due to asymmetric spatial gradients. Across all scenarios, the model predicts initial immune-mediated control followed by progressive macrophage depletion consistent with exhaustion patterns observed in solid tumors, with pressure-driven advection contributing cumulatively to tumor displacement.

The rest of the paper is organized as follows. Section 2 presents the coupled biological-mechanical model. Section 3 describes the neural network architecture and training methodology. Section 4 reports computational results, computational performance, and biological interpretation. Section 5 closes the paper with conclusions and limitations.

2  Mathematical Model

This section presents the complete mathematical formulation governing osteosarcoma tumor-immune-mechanical dynamics. The biological system is based on the reduced model proposed in [4], extended here with the IFN-γ variable and coupled to a poroelastic mechanical component. The full model is first presented in compact form, then discussed with biological and mechanical interpretation.

Model Formulation

The tumor microenvironment occupies a two-dimensional circular domain Ω⊂R2 with radius R=0.01 m, representing a tissue cross-section:

Ω={(x,y)∈R2:x2+y2≤R2}.

The boundary is denoted ∂Ω={(x,y)∈R2:x2+y2=R2}. The time evolution is considered over the interval 𝒯=[0,T] with T=500 days. All field variables are functions of the spatiotemporal coordinates (x,y,t)∈Ω×𝒯. Throughout the manuscript, the domain is assumed fixed in time: this is a common simplifying assumption for studies focusing on mass-exchange and transport dynamics over physiologically short time scales, and its limitations are discussed in Section 5.

The model comprises five biological fields Ψ=(C,m,Tc,N,Iγ) and three mechanical fields Φ=(ux,uy,p), defined on Ω×𝒯. Cancer cell density C(x,y,t) [cells/m2], macrophage density M(x,y,t) [cells/m2], cytotoxic T lymphocyte density Tc(x,y,t) [cells/m2], necrotic cell density N(x,y,t) [cells/m2], and interferon-gamma concentration Iγ(x,y,t) [dimensionless, normalized to a reference physiological concentration] constitute the biological state. The displacement field u(x,y,t)=(ux,uy)T [m] and interstitial pressure p(x,y,t) [Pa] characterize the mechanical state.

The biological variables satisfy

C,m,Tc,N,Iγ:Ω×𝒯→R+,(1)

and are governed by the following system of coupled reaction-diffusion-advection equations:

∂C∂t−∇⋅(DC∇C)+∇⋅(Cvf)=λCC(1−CC0)−δCTcTcC+λCMMC−δCIγIγC−δCC,(2)

∂M∂t−∇⋅(DM∇M)+∇⋅(Mvf)=λMIγIγM−δMm,(3)

∂Tc∂t−∇⋅(DTc∇Tc)+∇⋅(Tcvf)=λTcCCTc−δTcTc,(4)

∂N∂t=αNCδCTcTcC−δNN,(5)

∂Iγ∂t=λIγTcTc−δIγIγ.(6)

Eq. (2) governs cancer cell evolution, incorporating proliferation following logistic growth with carrying capacity C0, killing by cytotoxic T-cells at rate δCTc, promotion by macrophages at rate λCM, direct inhibition by IFN-γ at rate δCIγ, and natural death at rate δC. The diffusion coefficient DC quantifies random cell motility, while the advection term ∇⋅(Cvf) couples cancer cell transport to the interstitial fluid velocity vf.

Eq. (3) describes macrophage dynamics, where activation at rate λMIγ occurs through IFN-γ-mediated polarization toward the M1 phenotype. In the present formulation, the total macrophage density M is treated as a net pool dominated by the anti-tumoral M1 phenotype; M2 polarization and phenotype switching are not resolved explicitly, and the rate λMIγ therefore represents the net activation balance toward the M1 compartment. This simplification is consistent with the reduced formulation of [4] and with the cluster 1 cohort of the TARGET data, characterized by a moderately polarized M1-dominant infiltrate. Macrophage death proceeds at rate δM. The implications of this closed-pool assumption (in particular the absence of circulatory monocyte recruitment) are discussed in Section 5.

Eq. (4) governs cytotoxic T-cell evolution, with activation by cancer cells at rate λTcC and natural death at rate δTc. Eq. (5) models necrotic cell accumulation from T-cell-mediated cancer cell killing, with conversion factor αNC and removal rate δN. Necrotic cells are assumed immobile and to occupy the same reference specific volume as the cancer cells they replace; a refined treatment accounting for edema- or fragmentation-induced volume differences would require distinct specific volumes in Eq. (11) and is left to future work.

Finally, Eq. (6) governs IFN-γ dynamics with production by activated cytotoxic T-cells at rate λIγTc and natural decay at rate δIγ. The rapid decay rate δIγ=33.27 day−1 reflects the short half-life of this cytokine in vivo, introducing temporal stiffness into the system.

Remark 1: Eq. (6) models IFN-γ without spatial diffusion or advection. This simplification is justified by the separation of scales: the characteristic diffusion length of IFN-γ over its half-life (t1/2≈0.02 days) is ℓ=DIγt1/2∼𝒪(10−4) m for typical cytokine diffusivities DIγ∼10−10 m2/s [11], which is small compared to the domain radius R=0.01 m. Given the rapid decay, IFN-γ concentration at any point is effectively determined by local T-cell production.

The interstitial fluid velocity vf appearing in the advection terms of Eqs. (2)–(4) couples the biological system to the mechanical component through Darcy’s law:

vf=−κ∇p,(7)

where κ=5.0×10−14 m2/(Pa⋅s) denotes the hydraulic conductivity and p represents the interstitial fluid pressure. Parameter values for the biological component are presented in Table 1 and derive from fitting to experimental gene expression data from the TARGET cohort, corresponding to cluster 1 patients characterized by moderate immune infiltration and intermediate prognosis [3].

Tissue mechanical behavior follows Biot’s poroelastic theory [12], modeling the interaction between solid-phase deformation and fluid flow through porous media. Let u:Ω×𝒯→R2 denote the displacement field and p:Ω×𝒯→R the fluid pressure.

The governing equations comprise quasi-static mechanical equilibrium:

−∇⋅σ=0,(8)

and fluid pressure evolution:

∂∂t(c0p+α∇⋅u)−∇⋅(κ∇p)=0,(9)

where c0 denotes the specific storage coefficient, α the Biot–Willis coefficient, and κ the hydraulic conductivity already introduced in Eq. (7).

The total stress tensor σ accounts for elastic deformation, fluid pressure, and volumetric changes from cellular growth and death:

σ=2με+λ¯(∇⋅u)I−αpI−KV∗I,(10)

where ε=12(∇u+(∇u)T) denotes the strain tensor, μ and λ¯ are Lamé parameters (the overbar is introduced to distinguish the Lamé constant from the biological rate constants λ(⋅)), I is the identity tensor, K the bulk modulus, and V∗ the normalized volume change. The Lamé parameters relate to Young’s modulus E and Poisson’s ratio ν through:

μ=E2(1+ν),λ¯=Eν(1+ν)(1−2ν).

The volume change term V∗ couples mechanical and biological components, defined as:

V∗=AC(C−C0)+AN(N−N0),(11)

where AC=1/C0 and AN=1/N0 represent specific volumes associated with cancer and necrotic cells, respectively. Here the reference parameter C0 plays the dual role of logistic carrying capacity in Eq. (2) and of reference concentration for the volumetric coupling, while N0 denotes the reference necrotic concentration.

Substituting Eqs. (10) into (8) yields:

−μΔu−(μ+λ¯)∇(∇⋅u)+α∇p+K∇V∗=0,(12)

where Δu=(∇2ux,∇2uy)T denotes the componentwise (vector) Laplacian; throughout the rest of the manuscript ∇2 without vector decoration refers to the scalar Laplacian applied to scalar fields. Mechanical parameter values are reported in Table 2.

The problem requires specification of initial conditions at t=0:

C(x,0)=0.02C0,M(x,0)=M0(x),Tc(x,0)=Tc0≈1.90×107 cells/m2,N(x,0)=0.16N0,Iγ(x,0)=0.868,u(x,0)=0,p(x,0)=p0=0 Pa,(13)

where M0(x) varies with spatial position according to scenario-specific distributions, described in Section 4. The value Iγ(x,0)=0.868 represents a physiological baseline normalized to a reference cytokine concentration; it is significant but sub-saturated, since the reference equilibrium level corresponds to Iγ=1. This initial condition is imposed uniformly in space.

Homogeneous Neumann boundary conditions enforce no-flux for the biological fields:

∂C∂n=∂M∂n=∂Tc∂n=∂Iγ∂n=0on ∂Ω×𝒯,(14)

where n denotes the outward unit normal. Mechanical boundary conditions impose traction-free surfaces and continuity with surrounding tissue pressure:

σ⋅n=0,p=p0on ∂Ω×𝒯.

3  Physics-Informed Neural Network Approach

The PINN framework employs two separate neural networks coupled through the governing equations: a biological network 𝒩bio:R3→R5 mapping spatiotemporal coordinates to biological variables, and a mechanical network 𝒩mech:R3→R3 mapping coordinates to pressure and displacement components. The biological network produces:

𝒩bio(x,y,t;θbio)=[C(x,y,t),M(x,y,t),Tc(x,y,t),N(x,y,t),Iγ(x,y,t)]T,(15)

where θbio denotes trainable parameters (weights and biases). The mechanical network produces:

𝒩mech(x,y,t;θmech)=[p(x,y,t),ux(x,y,t),uy(x,y,t)]T(16)

with trainable parameters θmech.

The biological network architecture in Eq. (15) comprises four hidden layers with neuron counts [50, 100, 100, 50], implementing a bottleneck structure. This configuration follows from empirical observations showing that moderate depth (three to five layers) with symmetric expansion-contraction patterns effectively captures multiscale features in biological PDE systems [13,14]. The mechanical network in Eq. (16) employs three hidden layers with counts [30, 60, 30], reflecting the reduced representational capacity required for mechanical fields exhibiting greater regularity compared to biological dynamics [9]. Hidden layers utilize hyperbolic tangent activation functions, h(l)=tanh⁡(W(l)h(l−1)+b(l)), where h(l) denotes the output of layer l, W(l) the weight matrix, and b(l) the bias vector.

The biological network output layer applies the softplus activation to ensure strict positivity of the predicted concentrations:

𝒩bio(x,t;θbio)=softplus⁡(W(L)h(L−1)+b(L)),(17)

with softplus⁡(z)=ln⁡(1+ez) and L the total number of layers. Because softplus⁡(z)>0 for every z∈R, the output layer does not admit negative values, which prevents the non-physical oscillations that an identity or ReLU output would permit when populations approach extinction. However, the softplus output can become arbitrarily small without ever reaching exactly zero; the practical implication of this asymptotic behavior on the predicted tumor extinction in Scenario 4 is discussed in Section 4 and in the Limitations of Section 5. Network parameters are initialized using the Glorot normal distribution:

Wij(l)∼𝒩(0,2nin+nout),bi(l)∼𝒩(0,0.1),

where nin and nout denote input and output dimensions of layer l.

3.1 Loss Function Construction

Training proceeds by minimizing a composite loss function enforcing PDE residuals, initial conditions, and boundary conditions:

ℒtotal=ℒPDE+λICℒIC+λBCℒBC,(18)

where λIC and λBC balance component contributions. The PDE loss decomposes into biological and mechanical parts:

ℒPDE=ℒbio+ℒmech.(19)

The biological loss aggregates squared residuals of Eqs. (2)–(6):

ℒbio=1Nf∑j=1Nf[wC|RC|2+wM|RM|2+wTc|RTc|2+wN|RN|2+wIγ|RIγ|2](xj,tj),(20)

where {(xj,tj)}j=1Nf represent collocation points, w(⋅) are component weights, and R(⋅) denote PDE residuals. The cancer cell residual from Eq. (2) takes the form:

RC=∂C^∂t−DC∇2C^+∇⋅(C^v^f)−λCC^(1−C^C0)+δCTcT^cC^−λCMM^C^+δCIγI^γC^+δCC^,

where hat notation indicates network predictions and v^f=−κ∇p^ derives from Eq. (7) using the mechanical network output. The IFN-γ residual from Eq. (6) yields:

RIγ=∂I^γ∂t−λIγTcT^c+δIγI^γ.(21)

Residuals for remaining biological variables follow analogously from Eqs. (3)–(5). The mechanical loss combines displacement and pressure equation residuals:

ℒmech=1Nf∑j=1Nf[wu|Ru|2+wp|Rp|2](xj,tj).

The displacement residual from Eq. (12) reads:

Ru=−μΔu^−(μ+λ¯)∇(∇⋅u^)+α∇p^+K∇V^∗,

where V^∗=AC(C^−C0)+AN(N^−N0) from Eq. (11) couples to biological network outputs. The pressure residual from Eq. (9) becomes:

Rp=∂∂t(c0p^+α∇⋅u^)−∇⋅(κ∇p^).

Automatic differentiation computes all spatial and temporal derivatives appearing in residual expressions, ensuring exact gradient evaluation without introducing discretization errors.

The initial and boundary condition losses are

ℒIC=ℒICbio+ℒICmech,ℒICbio=1N0∑j=1N0[|C^j−0.02C0|2+|M^j−M0(xj)|2+|T^c,j−Tc0|2+|N^j−0.16N0|2+|I^γ,j−0.868|2],ℒICmech=1N0∑j=1N0[|u^j|2+|p^j−p0|2],

evaluated on {xj}j=1N0 sampled in Ω at t=0, and

ℒBC=ℒBCbio+ℒBCmech,ℒBCbio=1Nb∑j=1Nb[|∂C^∂n|2+|∂M^∂n|2+|∂T^c∂n|2+|∂I^γ∂n|2](xj,tj),ℒBCmech=1Nb∑j=1Nb[|σ^(xj,tj)⋅n(xj)|2+|p^(xj,tj)−p0|2],

with {(xj,tj)}j=1Nb sampling ∂Ω×𝒯.

3.2 Computational Implementation

The spatial domain Ω is the circular region with radius R=0.01 m representing initial tumor extent, while the temporal domain spans 𝒯=[0,500] days capturing long-term evolution. Collocation points totaling Nf=20,000 sample the spatiotemporal domain through probabilistic distributions ensuring uniform coverage with enhanced density in regions of interest. For each point, the angular coordinate θ∼𝒰(0,2π) is uniform, the radial coordinate r=RU with U∼𝒰(0,1) ensures uniform areal density (here r is the sampling radius, not the domain radius), and the temporal coordinate t=500⋅B with B∼Beta(1.5,3.0) emphasizes early evolution phases where dynamics proceed most rapidly. Cartesian coordinates follow from x=rcos⁡θ, y=rsin⁡θ.

Initial condition enforcement employs N0=2000 points uniformly distributed spatially at t=0. Cancer cells follow Gaussian distributions centered in the domain. Macrophage distributions vary by scenario as detailed in Section 4. T-cells exhibit nearly uniform concentrations with slight elevation near tumor regions. Necrotic cells and IFN-γ initialize to the baseline values specified in Eq. (13). Boundary condition enforcement samples Nb=2000 points along ∂Ω with uniform angular spacing and random temporal sampling across 𝒯.

3.3 Training Strategy

Network parameters are optimized using the Adam algorithm [15] with adaptive learning rate scheduling and gradient clipping. Exponentially decaying first and second moment estimates of the gradient use β1=0.9, β2=0.999, ε =10−8. The learning rate η follows polynomial decay from ηmax=10−3 to ηmin=10−5. Curriculum learning gradually increases the PDE loss weight:

λPDE(e)=min(1.0, 0.1+0.9⋅eNwarmup)

with Nwarmup=1000 epochs, allowing the network to first learn initial and boundary conditions before enforcing full PDE constraints. Gradient clipping at threshold 1.0 prevents instabilities from large gradient magnitudes. Early stopping terminates training after 150 consecutive non-improving epochs, restoring parameters corresponding to the minimum observed loss.

Component loss weights are set to λIC=10.0, λBC=5.0, wC=1.0, wM=1.0, wTc=1.0, wN=0.5, wIγ=0.1. The elevated initial condition weight reflects the well-known importance of accurate initialization for temporal PDE integration [13], while the reduced wIγ compensates for the fast-scale contribution of the stiff IFN-γ residual to the loss landscape. These values were determined starting from uniform weights and progressively reducing wIγ until stable convergence was obtained in a preliminary Scenario 1 run; the same weights were then used without further tuning across all scenarios.

The full procedure is reported in Algorithm 1.

4  Results

Computational experiments were executed on a workstation equipped with an AMD Ryzen 7 Pro 4750U CPU (8 cores, 16 threads, base clock 1.7 GHz) and 16 GB RAM, running a Linux distribution. Implementation utilized Python 3.10 with TensorFlow 2.15 (CPU backend, graph-compiled automatic differentiation; no GPU acceleration was used), NumPy 1.24, and Matplotlib 3.8.

Four scenarios were investigated, distinguished by initial macrophage spatial distributions (Fig. 1): Scenario 1 employed uniform distribution with M0=5.0×106 cells/m2 throughout the domain. Scenario 2 concentrated macrophages in a concentric ring at radius r=0.005 m with Gaussian profile width σ=0.002 m. Scenario 3 positioned two Gaussian clusters at diagonal corners (−0.004,−0.004) and (0.004,0.004) m. Scenario 4 accumulated macrophages within 0.002 m of domain boundaries following exponential decay inward. Initial distributions for the remaining biological variables remained identical across scenarios (Fig. 2). Throughout this section we adopt the convention that executed runs are described in the past tense while steady-pattern physical interpretations use the present tense.

Figure 1: Initial macrophage distributions (cells/m2) defining the four simulation scenarios.

Figure 2: Overview of the initial conditions for the cellular and cytokine distributions.

Training employed batch size 1000 with Nf=20,000 collocation points, N0=2000 initial condition points, and Nb=2000 boundary points. Adam parameters: β1=0.9, β2=0.999, ε =10−8, with learning rate decaying polynomially from 10−3 to 10−5. Curriculum learning ramped the PDE weight from 0.1 to 1.0 over 1000 epochs. Early stopping patience: 150 epochs. Maximum training was capped at 2000 epochs.

4.1 Convergence Behavior

Table 3 presents final loss values and training characteristics. Scenario 1 converged after 897 epochs with ℒtotal=0.0616, ℒbio=0.0089, ℒmech=6.0×10−6, requiring 7.2 h. Initial loss ℒtotal(0)=155.04 decreased rapidly during curriculum learning, reaching ℒtotal(200)=0.0809 by epoch 200, with subsequent monotonic descent: ℒtotal(400)=0.0697, ℒtotal(600)=0.0630, and early stopping at epoch 897.

Scenario 2 demonstrated stable convergence from initial ℒtotal(0)=51.32, progressing through ℒtotal(200)=0.1949 and ℒtotal(400)=0.1750, reaching final ℒtotal=0.1577 at epoch 631 (5.1 h). The concentric ring geometry produced spatial gradients raising the total loss relative to Scenario 1, while the biological loss remained at the same ℒbio=0.0089. Scenario 3 showed more fluctuating behavior from the highest initial loss ℒtotal(0)=194.88, with persistent oscillations between 0.101–0.108 through epochs 400–700 before stabilizing at ℒtotal=0.1067 (epoch 728, 5.9 h) and the highest biological loss (ℒbio=0.0127), consistent with asymmetric residuals in the cluster transition regions.

Scenario 4 achieved optimal performance despite starting from ℒtotal(0)=148.81, running the full 2000 epochs (14.3 h) and reaching ℒtotal=0.0557. Biological loss reached ℒbio=0.0024, the lowest across all scenarios by 73%, while ℒmech=1.2×10−6 remained minimal. The boundary configuration appears to create spatial gradients naturally compatible with no-flux boundary conditions, reducing constraint conflicts. We note that Scenario 4 never triggered early stopping, suggesting that the current patience threshold may be conservative for high-gradient configurations and that further training could potentially improve Scenarios 2–3 as well.

Mechanical loss components remained consistently low across scenarios (ℒmech<3.2×10−5), confirming excellent satisfaction of the poroelastic constraints. The elevation of biological loss relative to mechanical components reflects increased system complexity arising from the stiff IFN-γ kinetics.

4.2 Computational Performance

Per-epoch training time averaged 26±2 s on the hardware described above, yielding 7.2 h for the shortest run (Scenario 1, 897 epochs) and 14.3 h for the longest (Scenario 4, 2000 epochs). After training, inference at an arbitrary (x,y,t) is executed in roughly 10−3 s, so that post-training evaluation over a 200×200×200 visualization grid (8×106 points) requires on the order of 2 min. A direct, like-for-like FEM benchmark of the fully coupled poroelastic-biological system was not performed in the present work, and we prefer to report the PINN costs in isolation rather than comparing against figures we did not personally measure. We note, however, that adaptive mesh-based solvers for coupled tumor-immune reaction-diffusion-advection problems of comparable size typically report forward-simulation times of the order of minutes to tens of minutes on standard desktop hardware, and that stiff cytokine decay terms such as the one arising from δIγ are known to impose tight time-step restrictions on explicit schemes [8]. For this problem class, the PINN is therefore unlikely to outperform FEM as a one-shot forward solver; its competitive advantage lies rather in (i) continuous inference over Ω×𝒯 without remeshing, (ii) natural extension to inverse problems by making physical parameters trainable, and (iii) robustness against stiffness-induced adaptive-stepping instabilities. A rigorous FEM validation study is explicitly identified as future work in Section 5. The scaling of cost and accuracy with the collocation density was probed with two auxiliary runs on Scenario 1 with Nf∈{5000,10,000} (keeping N0,Nb unchanged). The final biological loss increased from ℒbio=0.0089 at Nf=20,000 to ℒbio≈0.013 at Nf=10,000 and ℒbio≈0.021 at Nf=5000, while training time scaled approximately linearly. The density Nf=20,000 was therefore retained as a reasonable accuracy-cost compromise for the present architecture and domain size. Regarding design-choice sensitivity, the bottleneck architecture [50,100,100,50] was compared on Scenario 1 during the development phase against a uniform-width [100,100,100,100] and a shallower [50,100,50]; the bottleneck reached the lowest biological loss by approximately 15% at comparable training cost, motivating its selection. A preliminary local sensitivity on the two parameters governing the novel IFN-γ pathway, δIγ and λMIγ, perturbed one-at-a-time by ±20% on Scenario 1, preserved the qualitative dynamics but shifted the timing of the cancer peak by 15–20 days; a systematic sensitivity analysis over the full parameter set is deferred to future work via inverse-PINN techniques.

4.3 Spatiotemporal Evolution of Cancer Cells

Fig. 3 displays cancer cell distributions at representative time points across scenarios. For readability, color scales are normalized per sub-panel because the predicted cancer concentrations span more than twelve orders of magnitude across scenarios; a global logarithmic scale was tested and preserved the qualitative comparison but compressed intra-scenario features. The values quoted below are the raw field values and should be read together with the panel-specific color legends.

Figure 3: Spatiotemporal evolution of the cancer cell fields across the four different macrophage distribution scenarios.

Scenario 1 (uniform macrophages) showed an initial Gaussian concentration centered at the origin with peak C(0,0,0)=2.69×108 cells/m2 expanding through combined diffusion and advection. By t=100 days, concentration increased to Cmax(100)=8.94×109 cells/m2 (≈0.67C0), representing 33-fold growth reflecting early proliferation dominance. The spatial distribution shifted toward (x,y)=(0.003,0.002) m, driven by pressure gradient-induced fluid flow.

At t=300 days, peak concentration reached Cmax(300)=9.58×109 cells/m2 (0.71C0), approaching carrying capacity with growth slowing per the logistic term λCC(1−C/C0) in Eq. (2). The distribution homogenized while maintaining rightward bias, with tumor diameter expanding from initial 0.004 m to approximately 0.012 m. By t=500 days, concentration declined to Cmax(500)=3.21×109 cells/m2 (0.24C0), indicating late immune control. The spatial pattern revealed concentration shifting toward (−0.003,−0.003) m with fragmented distribution suggesting heterogeneous immune infiltration.

Scenario 2 (concentric ring) created spatial heterogeneity in immune pressure. Central tumor growth proceeded rapidly to Cmax(100)=9.41×109 cells/m2 by t=100 days, but expansion toward the ring region encountered elevated macrophage-mediated killing. At t=300 days, asymmetric development emerged with higher density in the upper-right quadrant (C≈9.71×109 cells/m2) and suppressed growth in the lower-left (C≈7.2×109 cells/m2) where macrophages remained elevated. By t=500 days, the peak Cmax(500)=9.48×109 cells/m2 persisted in the upper-right, while central regions showed lower density.

Scenario 3 (two clusters) generated pronounced asymmetry. Early growth (t=100 days) showed selective suppression near cluster locations: concentration reached only C≈5.2×109 cells/m2 at (−0.004,−0.004) m vs. C≈9.1×109 cells/m2 in the opposite low-macrophage region. By t=300 days, preferential expansion toward the upper-left quadrant achieved C≈9.83×109 cells/m2. Late evolution (t=500 days) showed overall decrease to Cmax(500)=7.26×109 cells/m2 but maintained a spatial gradient.

Scenario 4 (boundary concentration) exhibited qualitatively different dynamics. The centralized tumor, expanding through a macrophage-poor interior, reached only Cmax(100)=2.84×10−4 cells/m2 at t=100 days. Because the model is continuous and the domain area is πR2≈3.14×10−4 m2, such a concentration corresponds to fractional cell numbers per unit area with no direct biological counterpart. We therefore interpret the Scenario 4 prediction at t∈[50,150] days as effective extinction of the tumor field, not as a literal cell count. This extreme suppression is consistent with the advective mechanism: pressure gradients drive fluid flow per Eq. (7), advecting cancer cells via ∇⋅(Cvf) toward boundaries where macrophage concentrations peak, producing enhanced tumor-immune contact and elevated killing rates. The subsequent increase to Cmax(300)=1.07×104 and Cmax(500)=1.08×104 cells/m2 should be read as continuous-field regrowth from vanishingly small residuals preserved by the softplus activation (which never outputs exactly zero), i.e., model-level recurrence from sub-threshold residual disease, not biological re-seeding from an exogenous source. Imposing a numerical extinction threshold Cmin below which C is clamped to zero would eliminate this artifact and is identified as a future refinement in Section 5.

4.4 Temporal Dynamics of Immune Components

Fig. 4 presents domain-averaged temporal profiles for all biological variables and pressure across scenarios. All axis quantities carry the units reported in the caption: cell densities in cells/m2, IFN-γ dimensionless, pressure in Pa, time in days.

Figure 4: Temporal dynamics of domain-averaged cell densities, cytokine concentration, and interstitial fluid pressure over 500 days.

Cancer cell evolution (panel A) shows a universal early growth phase reaching peaks between t=100–150 days, followed by scenario-dependent trajectories. Scenario 1 declined gradually from ⟨C⟩(100)≈6.2×109 cells/m2 to ⟨C⟩(500)≈2.1×109 cells/m2. Scenarios 2–3 maintained elevated concentrations: ⟨C⟩(500)≈7.8×109 cells/m2 (Scenario 2), ⟨C⟩(500)≈5.9×109 cells/m2 (Scenario 3), indicating partial immune escape. Scenario 4 exhibited a dramatically suppressed trajectory with ⟨C⟩(100)≈1.8×10−4 cells/m2, remaining below 104 cells/m2 throughout evolution.

Macrophage temporal profiles (panel B) revealed scenario-dependent depletion patterns. Scenario 1 maintained relatively stable ⟨M⟩≈5.1×107 cells/m2 through t=200 days, then declined to ⟨M⟩(500)≈3.8×107 cells/m2 (25% reduction). Scenarios 2–3 showed more pronounced depletion: Scenario 2 decreased from 5.0×107 to ⟨M⟩(500)=2.1×107 cells/m2 (58% reduction). Scenario 4 demonstrated the most severe depletion, dropping from initial domain-average ⟨M⟩(0)≈5.0×107 cells/m2 to ⟨M⟩(500)=1.4×107 cells/m2 (72% reduction), with boundary regions declining to M(500)=8.2×106 cells/m2 (95% local reduction). This depletion reflects the death rate δM=0.03 day−1 operating continuously while IFN-γ-mediated activation λMIγIγM diminishes with declining cytokine levels. We emphasize that in the present closed-system formulation the only counterbalance to macrophage death is IFN-γ-mediated activation. In vivo, circulatory monocyte recruitment would partially replenish the tissue macrophage pool; the depletion values predicted here should therefore be regarded as upper bounds relative to the physiological scenario, and the predicted long-term quasi-equilibria would shift toward higher macrophage densities once a recruitment source term is introduced. This limitation is revisited in Section 5.

T-cell dynamics (panel C) exhibited biphasic behavior. Initial Tc(0)=1.90×107 cells/m2 increased during early cancer expansion as the activation term λTcCCTc dominated, reaching peaks between t=80–120 days: Tcmax=2.78×107 (Scenario 1), 2.91×107 (Scenario 2), 2.64×107 (Scenario 3), 2.52×107 cells/m2 (Scenario 4). Subsequent decline occurred as cancer concentrations decreased. By t=500 days: Tc(500)=1.82,1.76,1.69,1.58×107 cells/m2, within 10%–20% of baseline.

Necrotic cell accumulation (panel D) reflected T-cell-mediated killing history. Starting from N(0)=1.07×108 cells/m2, concentrations initially decreased to minima around t=50 days as the removal rate δNN exceeded production, then peaked at t=150–180 days: Nmax=1.32,1.41,1.28,1.18×108 cells/m2 for Scenarios 1–4, stabilizing at N(500)=1.14,1.19,1.08,0.97×108 cells/m2.

IFN-γ temporal evolution (panel E) exhibited characteristic biphasic profiles. From Iγ(0)=0.868, concentrations rose rapidly during T-cell expansion to t=90–110 day peaks: Iγmax=0.951,0.968,0.932,0.913. Peak timing correlated closely with T-cell maxima, confirming the production term λIγTcTc as the primary driver. By t=500 days: Iγ(500)=0.821,0.796,0.763,0.713. The rapid decay rate δIγ=33.27 day−1 maintains cytokine levels within 10%–15% of the equilibrium value Iγeq=(λIγTc/δIγ)Tc for the instantaneous T-cell concentration.

Pressure evolution (panel F) coupled directly to tumor growth through the volume change term V∗ in Eqs. (11) and (12). Maxima pmax=4.21,4.38,3.87,0.42 Pa occurred at t=100–120 days, with relaxation to p(500)=1.82,1.96,1.64,0.18 Pa. The lower pressure in Scenario 4 reflects much smaller cancer concentrations and correspondingly reduced volumetric forcing. These pressure magnitudes lie in the range of interstitial fluid pressure reported by Netti et al. [11] for solid tumors at early stages.

4.5 Biological Interpretation

The extended PINN model captures essential immune-tumor feedback mechanisms operating across multiple timescales. Early evolution exhibits immune-dominated dynamics: T-cell activation by cancer drives IFN-γ production, which simultaneously inhibits tumor proliferation and enhances macrophage M1 polarization. This creates negative feedback: immune activation increases IFN-γ, which suppresses tumor growth and activates macrophages, enhancing killing and reducing cancer burden. The loop operates most effectively when spatial overlap between tumor and immune cells remains high, explaining why Scenarios 2–4 show stronger early suppression than the uniform Scenario 1.

Intermediate evolution transitions toward tumor escape. As concentrations approach carrying capacity C0, proliferation slows through the logistic term regardless of immune pressure. Simultaneously, sustained activity depletes macrophage populations, reducing killing capacity. T-cell decline, as tumor suppression decreases activation stimulus, further lowers IFN-γ production, creating positive feedback enabling immune evasion.

Late evolution reaches scenario-dependent quasi-equilibria. Scenario 1 maintains moderate cancer burden with depleted but functional immune response. Scenario 2 sustains high tumor load, indicating partial escape. Scenario 3 shows intermediate behavior with persistent spatial heterogeneity. Scenario 4 achieves the strongest suppression through early intense immune engagement, though severe macrophage depletion suggests transient rather than stable control.

A practical consequence concerns the therapeutic window, understood as the time interval during which the tumor field remains below a clinically relevant threshold. In the boundary configuration, the model predicts ⟨C⟩<104 cells/m2 for the entire simulated interval of 500 days, which however coincides with a 72% macrophage loss, signaling that this state is unlikely to be stable beyond the present simulation window without external immune support (e.g., monocyte recruitment or checkpoint-inhibitor therapy). A precise quantification of this window against a patient-specific viability or Allee threshold will require calibration against longitudinal imaging data and is left to future work.

IFN-γ kinetics prove crucial for these dynamics. The rapid decay rate prevents accumulation, requiring continuous T-cell production to maintain physiological levels. This creates responsiveness: activation immediately elevates IFN-γ, rapidly modulating tumor inhibition and macrophage polarization; conversely, T-cell decline produces swift IFN-γ reduction within days, allowing tumor recovery. This responsiveness contrasts with slower cellular timescales, positioning IFN-γ as a high-frequency signaling molecule mediating immune-tumor communication.

Mechanical coupling through pressure-driven advection influences tumor topology. The boundary configuration of Scenario 4 generates inward fluid flow opposing outward cancer diffusion, enhancing early immune contact. Advective velocities, while small compared to cell migration, produce cumulative transport effects over hundreds of days comparable to diffusive displacement. We emphasize that this interpretation suggests rather than demonstrates a causal role for advective transport, since the present study does not include a mechanical-only ablation run.

5  Conclusions and Limitations

An extended Physics-Informed Neural Network framework was presented for modeling osteosarcoma tumor-immune-mechanical dynamics, incorporating interferon-gamma as a fifth biological variable coupled with poroelastic tissue mechanics. Four macrophage distribution scenarios were investigated. The framework achieved stable convergence in all cases, with total loss values between 0.0557 and 0.1577. Mechanical residuals remained below 3.2×10−5 throughout. The boundary-concentrated configuration reached the lowest biological loss, while the cluster scenario exhibited the highest due to asymmetric gradient complexity. Temporal stiffness from the IFN-γ decay rate was managed through curriculum learning and adaptive component weighting without compromising biological fidelity.

The predicted dynamics exhibited initial immune-mediated tumor suppression followed by progressive macrophage depletion across all scenarios, consistent with observed patterns of immune exhaustion in solid tumors. Pressure-driven advection influenced tumor spatial evolution, and boundary macrophage accumulation appeared to enhance early immune-tumor contact through inward fluid transport.

Limitations

Several limitations of the present work should be acknowledged. First, the predicted fields have not been benchmarked against an independent reference solution (FEM or analytical). The low loss values indicate satisfaction of the governing equations at the collocation points, but not pointwise accuracy of the predicted fields; a quantitative FEM verification is identified as a priority follow-up. Second, all simulations used a single parameter set (cluster 1 of TARGET) without a formal sensitivity study; local sensitivities with respect to λCM, δCIγ, αNC, δIγ, and λMIγ are planned via inverse-PINN techniques. Third, the model assumes a fixed spatial domain over 500 days, neglecting large-strain tumor growth and associated mesh motion; in the parameter regime studied, predicted displacements remained small enough (<10−4 m) to support this approximation, but at larger strains or longer horizons a moving-domain or ALE formulation would be needed. Fourth, macrophage dynamics are closed-system: circulatory monocyte recruitment is not modeled, so the predicted depletion magnitudes are upper bounds. Fifth, M1/M2 phenotypes are not resolved separately; λMIγ represents a net M1-dominant activation. Sixth, the softplus output layer prevents true extinction of any biological variable, which in Scenario 4 produced continuous-field concentrations below any physiologically meaningful value; a future refinement will introduce a numerical extinction threshold Cmin and treat the system as piecewise below this threshold. Seventh, no uncertainty quantification was performed: reporting mean ± standard deviation of key metrics over multiple seeds of the collocation point sampling and weight initialization is a natural next step.

Future extensions include three-dimensional domain formulation, patient-specific parameter calibration via inverse PINN approaches [7], systematic validation against FEM reference solutions, and incorporation of macrophage M1/M2 phenotype switching.

Acknowledgement: Pasquale De Luca and Livia Marcellino are members of the Gruppo Nazionale Calcolo Scientifico–Istituto Nazionale di Alta Matematica (GNCS-INdAM).

Funding Statement: This work was partially supported by INdAM-GNCS project 2026 “Metodi Numerici e Machine Learning per Sistemi Dinamici Complessi: Identificazione, Approssimazione”, ref. no. CUP E53C25002010001.

Author Contributions: The authors confirm contribution to the paper as follows: study conception and design: Pasquale De Luca; data collection: Pasquale De Luca and Livia Marcellino; analysis and interpretation of results: Pasquale De Luca and Livia Marcellino. All authors reviewed and approved the final version of the manuscript.

Availability of Data and Materials: The source code implementing the PINN framework described in this paper is available upon reasonable request from the corresponding author.

Ethics Approval: Not applicable.

Conflicts of Interest: The authors declare no conflicts of interest.

References

1. Mirabello L, Troisi RJ, Savage SA. International osteosarcoma incidence patterns in children and adolescents, middle ages and elderly persons. Int J Cancer. 2009;125(1):229–34. doi:10.1002/ijc.24320. [Google Scholar] [PubMed] [CrossRef]

2. Cole S, Gianferante DM, Zhu B, Mirabello L. Osteosarcoma: a surveillance, epidemiology, and end results program-based analysis from 1975 to 2017. Cancer. 2022;128(11):2107–18. doi:10.1002/cncr.34163. [Google Scholar] [PubMed] [CrossRef]

3. Le T, Su SY, Kirshtein A, Shahriyari L. Data-driven mathematical model of osteosarcoma. Cancers. 2021;13(10):2367. doi:10.3390/cancers13102367. [Google Scholar] [PubMed] [CrossRef]

4. Le T, Su SY, Shahriyari L. Investigating optimal chemotherapy options for osteosarcoma patients through a mathematical model. Cells. 2021;10(8):2009. doi:10.3390/cells10082009. [Google Scholar] [PubMed] [CrossRef]

5. Di Vicino A, De Luca P, Marcellino L. First experiences on exploiting physics-informed neural networks for approximating solutions of a biological model. In: Paszynski M, Barnard AS, Zhang YJ, editors. Computational science—ICCS 2025 Workshops. Cham, Switzerland: Springer; 2025. Vol. 15910, p. 18–26. doi:10.1007/978-3-031-97567-7_2. [Google Scholar] [CrossRef]

6. Valentino C, Pagano G, Conte D, Paternoster B, Colace F, Casillo M. Step-by-step time discrete physics-informed neural networks with application to a sustainability PDE model. Math Comput Simul. 2025;230:541–58. doi:10.1016/j.matcom.2024.11.023. [Google Scholar] [CrossRef]

7. Raissi M, Perdikaris P, Karniadakis GE. Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. J Comput Phys. 2019;378:686–707. doi:10.1016/j.jcp.2018.10.045. [Google Scholar] [CrossRef]

8. Wang MY, Li B, Feng XQ, Gao H. A mechano-immunological framework for lymph node remodeling during inflammation and homeostasis. J Mech Phys Solids. 2026;206(Part A):106347. doi:10.1016/j.jmps.2025.106347. [Google Scholar] [CrossRef]

9. Haghighat E, Raissi M, Moure A, Gomez H, Juanes R. A physics-informed deep learning framework for inversion and surrogate modeling in solid mechanics. Comput Methods Appl Mech Eng. 2021;379:113741. doi:10.1016/j.cma.2021.113741. [Google Scholar] [CrossRef]

10. Frieboes HB, Jin F, Chuang YL, Wise SM, Lowengrub JS, Cristini V. Three-dimensional multispecies nonlinear tumor growth—II: tumor invasion and angiogenesis. J Theor Biol. 2010;264(4):1254–78. doi:10.1016/j.jtbi.2010.02.036. [Google Scholar] [PubMed] [CrossRef]

11. Netti PA, Baxter LT, Boucher Y, Skalak R, Jain RK. Time-dependent behavior of interstitial fluid pressure in solid tumors: implications for drug delivery. Cancer Res. 1995;55(22):5451–8. [Google Scholar] [PubMed]

12. Biot MA. General theory of three-dimensional consolidation. J Appl Phys. 1941;12(2):155–64. doi:10.1063/1.1712886. [Google Scholar] [CrossRef]

13. Yazdani A, Lu L, Raissi M, Karniadakis GE. Systems biology informed deep learning for inferring parameters and hidden dynamics. PLoS Comput Biol. 2020;16(11):e1007575. doi:10.1371/journal.pcbi.1007575. [Google Scholar] [PubMed] [CrossRef]

14. Lagergren JH, Nardini JT, Lavigne GM, Rutter EM, Flores KB. Biologically-informed neural networks guide mechanistic modeling from sparse experimental data. PLoS Comput Biol. 2020;16(12):e1008462. doi:10.1371/journal.pcbi.1008462. [Google Scholar] [PubMed] [CrossRef]

15. Kingma DP, Ba J. Adam: a method for stochastic optimization. arXiv:1412.6980. 2014. [Google Scholar]