Physics-Informed Neural Networks for Biological $2\mathrm{D}{+}t$ Reaction-Diffusion Systems

Abstract: Physics-informed neural networks (PINNs) provide a powerful framework for learning governing equations of dynamical systems from data. Biologically-informed neural networks (BINNs) are a variant of PINNs that preserve the known differential operator structure (e.g., reaction-diffusion) while learning constitutive terms via trainable neural subnetworks, enforced through soft residual penalties. Existing BINN studies are limited to $1\mathrm{D}{+}t$ reaction-diffusion systems and focus on forward prediction, using the governing partial differential equation as a regulariser rather than an explicit identification target. Here, we extend BINNs to $2\mathrm{D}{+}t$ systems within a PINN framework that combines data preprocessing, BINN-based equation learning, and symbolic regression post-processing for closed-form equation discovery. We demonstrate the framework's real-world applicability by learning the governing equations of lung cancer cell population dynamics from time-lapse microscopy data, recovering $2\mathrm{D}{+}t$ reaction-diffusion models from experimental observations. The proposed framework is readily applicable to other spatio-temporal systems, providing a practical and interpretable tool for fast analytic equation discovery from data.

• The paper introduces an end-to-end pipeline that combines biologically informed neural networks with PySR symbolic regression to learn closed-form density-dependent diffusion and growth functions from noisy 2D+t microscopy data.
• Applied to three PC9 lung cancer replicates, the method consistently recovered increasing diffusion and decreasing growth, with exponential diffusion forms and mostly linear growth laws, while forward simulations matched observed cell counts well for two replicates.
• The study shows that early stopping is essential for reliable equation discovery because longer training can reduce aggregate validation loss while overfitting the PDE and producing unrealistically low diffusion estimates; each full run takes about 10 minutes on a CPU laptop.

Overview

This paper extends biologically-informed neural networks (BINNs) to $2\mathrm{D}{+}t$ reaction–diffusion systems and embeds them in a complete equation-learning pipeline that combines data preprocessing, BINN training, and symbolic regression (SR) post-processing. Prior BINN work was restricted to $1\mathrm{D}{+}t$ systems and used the governing PDE primarily as a regulariser for forward prediction rather than as an explicit identification target. The authors address both gaps: they generalise the architecture to two spatial dimensions plus time, and they add a PySR-based post-processing stage that converts learned neural surrogates into closed-form analytic expressions. The framework is demonstrated on time-lapse microscopy data of PC9 lung cancer cell populations, recovering interpretable density-dependent diffusion and growth functions from three experimental replicates.

Methodological framework

The governing model is a nonlinear reaction–diffusion equation,

$$\frac{\partial u}{\partial t} = \nabla \cdot \left[D(u)\,\nabla u\right] + G(u)\,u,$$

where $D(u)$ is a density-dependent diffusivity and $G(u)$ a per-capita growth rate; the equation-learning task is to recover these constitutive functions from observations of $u(\mathbf{x},t)$ alone.

The BINN comprises three multilayer perceptrons trained jointly: $\mathrm{NN}_u$ maps $(\mathbf{x},t)$ to predicted density $\hat u$ (three hidden layers of 64 SiLU units, softplus output); $\mathrm{NN}_D$ and $1\mathrm{D}{+}t$0 map $1\mathrm{D}{+}t$1 to diffusion and growth predictions (narrow 4-unit hidden layers, softplus and linear outputs respectively). Training minimises a weighted sum of an ordinary least-squares data loss and a PDE residual loss evaluated at randomly sampled collocation points, with derivatives computed by automatic differentiation. Notably, no monotonicity or bounding constraints are imposed on $1\mathrm{D}{+}t$2 or $1\mathrm{D}{+}t$3, so any structure recovered in these functions emerges purely from the data.

Several design choices target robustness in low-data, high-noise regimes. Data are split into fixed random training/validation partitions (five 80/20 splits per replicate form a small ensemble), early stopping halts training when validation loss fails to improve by at least 5% within a patience window, and loss weights are normalised so that $1\mathrm{D}{+}t$4. The biological loss term is set to zero because the experimental setting provides little prior information restricting the functional forms of $1\mathrm{D}{+}t$5 and $1\mathrm{D}{+}t$6.

The SR stage uses PySR's evolutionary search over expressions built from $1\mathrm{D}{+}t$7 and arithmetic operators, with expression complexity capped at ten nodes. Ensemble predictions are constructed only over well-supported densities — the intersection of central 90% intervals across TV splits — weighted by replicate-specific density distributions. SR is run ten times with different seeds; candidate expressions are reduced to coefficient-free templates, the most frequent template is selected (ties broken toward simplicity), and coefficients come from the lowest-error matching candidate.

Application to lung cancer microscopy data

The pipeline is applied to three replicates of PC9 cell cultures imaged every 4 hours for up to 2.5 days. Raw cell coordinates are binned at approximately 0.1 mm resolution — chosen to balance the spatial scale of cell–environment sensing against the need for enough cells per bin to yield a smooth density field — producing tensors of size $1\mathrm{D}{+}t$8 with $1\mathrm{D}{+}t$9 across replicates.

Density learning. The trained $$\frac{\partial u}{\partial t} = \nabla \cdot \left[D(u)\,\nabla u\right] + G(u)\,u,$$0 smooths stochastic binning noise while preserving population-scale spatial structure, including heterogeneous initial conditions with clustered high-density regions and near-empty areas that fill in over time through combined diffusion and growth. Variability across TV splits is acknowledged as expected given the stochasticity of the data.

Diffusion identification. All learned diffusion functions increase monotonically with density despite no imposed monotonicity constraint, which the authors interpret as evidence of shared underlying biology. Magnitudes are on the order of $$\frac{\partial u}{\partial t} = \nabla \cdot \left[D(u)\,\nabla u\right] + G(u)\,u,$$1 at low densities, rising roughly twofold (replicate 2) to tenfold (replicates 1 and 3) at high densities; this discrepancy is attributed to differences in initial conditions, since replicate 2 reaches high plate coverage earlier, limiting observable high-diffusion regimes. Symbolic regression recovers exponential forms: replicates 1 and 2 share the template $$\frac{\partial u}{\partial t} = \nabla \cdot \left[D(u)\,\nabla u\right] + G(u)\,u,$$2 (with $$\frac{\partial u}{\partial t} = \nabla \cdot \left[D(u)\,\nabla u\right] + G(u)\,u,$$3), while replicate 3 yields $$\frac{\partial u}{\partial t} = \nabla \cdot \left[D(u)\,\nabla u\right] + G(u)\,u,$$4.

Growth identification. Learned growth functions decrease monotonically across all replicates, again without constraints, and show stronger cross-replicate and within-replicate agreement than diffusion. The recovered symbolic forms are predominantly linear: $$\frac{\partial u}{\partial t} = \nabla \cdot \left[D(u)\,\nabla u\right] + G(u)\,u,$$5 for replicates 2 and 3, and $$\frac{\partial u}{\partial t} = \nabla \cdot \left[D(u)\,\nabla u\right] + G(u)\,u,$$6 for replicate 1. The authors conclude that growth is more readily identifiable than diffusion in this setting and plays the dominant role in the learned dynamics.

Forward validation. Total cell-count trajectories computed from (i) the density MLP, (ii) forward solves using the MLP surrogates $$\frac{\partial u}{\partial t} = \nabla \cdot \left[D(u)\,\nabla u\right] + G(u)\,u,$$7, and (iii) forward solves using the symbolic pair $$\frac{\partial u}{\partial t} = \nabla \cdot \left[D(u)\,\nabla u\right] + G(u)\,u,$$8 closely reproduce observed counts for replicates 1 and 2. Replicate 3 is less accurate at later times: its count trajectory has greater curvature, and a sharp increase between the third and penultimate time points is not captured, suggesting additional dynamics beyond the assumed reaction–diffusion structure. This is a substantive limitation on the generality of the recovered equations for that replicate.

Computational cost and convergence behaviour

A full pipeline run for one replicate — five TV splits plus ten SR runs — takes approximately 10 minutes on a standard laptop (Apple M3 Pro, CPU-only), making the framework practical for routine EQL workflows.

An ablation over early-stopping patience reveals an important failure mode. Increasing patience from 500 to 1000 to 2000 epochs raises median runtime from roughly 141 s to 292 s to 578 s, while validation loss improves only marginally (~3% per doubling). More critically, although total validation loss continues decreasing with longer training, the PDE residual reaches a minimum and then increases — signalling a transition from physically meaningful solutions to overfitting. This manifests concretely in the learned diffusion function, which is progressively driven toward unrealistically low values at later epochs, whereas growth remains comparatively stable. The practical implication is that aggressive convergence can degrade equation discovery even when aggregate loss metrics suggest improvement; $$\frac{\partial u}{\partial t} = \nabla \cdot \left[D(u)\,\nabla u\right] + G(u)\,u,$$9 is selected as preferable on both cost and fidelity grounds.

Limitations and open questions

The authors identify several limitations directly. Diffusion was identified less robustly than growth, partly because growth dominates the observed dynamics and partly because high plate coverage limits spatial redistribution, restricting the data's information content about diffusive transport. The choice of bin size involves a trade-off between biological scale and statistical smoothing that is resolved heuristically here. Uncertainty quantification is handled only through ensembling over TV splits and repeated SR runs rather than principled probabilistic treatment. Replicate 3's unexplained late-time dynamics indicate that the fixed reaction–diffusion ansatz may be incomplete for some experimental conditions. Open questions include which experimental designs enhance diffusion identifiability, how to provide more rigorous uncertainty estimates, and whether the framework extends to more complex biological dynamics beyond the current ansatz.

Conclusion

This work delivers a computationally lightweight, end-to-end framework for discovering closed-form $D(u)$0 reaction–diffusion equations from noisy biological imaging data, combining BINN-based learning of constitutive terms with template-selection symbolic regression. Applied to PC9 lung cancer cell populations, it recovers monotonically increasing, exponentially structured diffusion and monotonically decreasing, approximately linear growth functions consistent across replicates, without imposing those structures as constraints. The demonstration that meaningful equation discovery is feasible under substantial individual-level stochasticity — at a cost of minutes on commodity hardware — is the paper's principal contribution, tempered by the acknowledged difficulty of identifying diffusion when growth dominates and by evidence that the reaction–diffusion ansatz does not fully explain all replicates.