







Cells sense and respond to their environment through signalling pathways, and the dynamics of these pathways shape cell fate even within genetically identical populations. Two largely separate computational traditions describe this behaviour: mechanistic differential-equation models and representation-learning methods. Mechanistic models encode pathway topology and kinetics but cannot easily represent variation arising outside the modelled pathway. Representation learning, instead, maps genome-wide measurements onto low-dimensional manifolds but offers no mechanistic account of how the resulting cell states execute their functions. Reconciling these views, explaining signalling heterogeneity in a manner that is at once data-driven and mechanistically interpretable, has remained difficult. Here we introduce deep mechanistic models (DMMs), which couple semi-supervised representation learning to an ordinary-differential-equation model of EGFR/MAPK signalling, trained end-to-end so that the learnt representation and mechanistic parametrisation inform each other. Applying DMMs to multiplexed signalling data from 63 breast cancer cell lines, we show that the models generalise to held-out cell lines and attribute most heterogeneity to pathway-extrinsic factors, namely baseline ERBB2 activation and a Ca²⁺/p38 signalling axis, rather than to variation in core MAPK components. Where the models fail, the discrepancies pinpoint rare signalling-altering mutations and recurrent programmes, including a putative AMPK– BRAF MEK-inhibitor-resistance axis and a cytoskeletal programme. We further find that mechanistic integration of EGFR receptor levels reshapes the learnt representation, rendering a molecular and a systems-level account of the same data equivalent in predictive power. DMMs thus offer a general framework for fusing mechanism with learning, naturally extendable to further modalities such as imaging, and simultaneously turn model failure into a systematic route to discover and evaluate candidate biology.
Gene regulatory networks: from correlative models to causal explanations
Gene regulatory networks (GRNs) explain how the genome controls cellular behaviour and tissue morphogenesis, serving to connect molecular mechanism to functional output. Single-cell technologies now provide descriptions of these networks with unprecedented detail, but this advance has also revealed gene regulatory systems that are too complex for our existing conceptual frameworks. GRNs, which should provide mechanistic explanations, are increasingly reduced to statistical correlations — ‘hairballs’ that fail to capture molecular causation. Here, we explore why this dilemma exists and propose a path forward. We argue that methods in ‘representation learning’ can be used to model GRNs, without needing to capture every molecular detail. For this framework, we advocate three linked principles: models must be inherently mechanistic, with structures grounded in cellular and evolutionary biology; molecular principles and constraints must be used to reduce the solution space for learning GRN models; and more sophisticated forms of experimental perturbation and synthetic biological engineering are needed to train models and test predictions. By reimagining GRNs through these principles, we can bridge the gap from data abundance to new conceptual understanding.

VAPOR: A variational autoencoder with transport operators to disentangle cellular gene expression dynamics of co-occurring biological processes in time and space
Abstract Single-cell and spatial transcriptomics enable the analysis of cellular states and dynamics in gene expression, revealing how diverse biological processes relate to these states over time and space. To study these dynamics, trajectory inference methods order cells along computationally inferred paths to reconstruct gradual transitions in cell states. However, by encouraging smooth and continuous trajectories, these approaches tend to conflate co-occurring processes-such as proliferation, maturation, and spatial organization-that are jointly reflected in gene expression, potentially overlooking process-specific gene expression dynamics. To address this, we developed VAPOR, which integrates a variational autoencoder with transport operators to model and disentangle cellular gene expression dynamics for potentially co-occurring biological processes. VAPOR inputs single-cell (or spatial) gene expression data into a variational autoencoder (VAE) to learn the latent states of cells and then models their latent dynamics as an ordinary differential equation. The latent dynamics are further decomposed into process-specific components parameterized by transport operators (TOs) and their corresponding process weights. Each TO defines a process-specific dynamics, and its weight for each cell quantifies the process's contribution to the cell dynamics. After assessment by simulation studies, we applied VAPOR with benchmarking to real data, including time-course scRNA-seq from postconceptual human brain development, spatial transcriptomics of the mouse hippocampus, and cross-species scRNA-seq spanning human and macaque first-trimester forebrain development. In these applications, VAPOR has identified a variety of temporal and spatial co-occurring processes, such as cell cycle, gliogenesis, neurogenesis, and neuronal migration, along with associated dynamic genes, including those species-specific to human and macaque development. VAPOR is available as an open-source tool for general-purpose use. ### Competing Interest Statement The authors have declared no competing interest.

How AI Can Advance Mathematical Biology: Opportunities, Challenges, and Future Directions
Mathematical biology has long relied on mechanistic models, including ordinary and partial differential equations, stochastic systems, and agent-based models, to study biological processes across scales. These approaches remain central because they provide structure, interpretability, and biological insight. However, modern biological and biomedical data, including longitudinal clinical records, medical imaging, and multi-omics measurements, are often high-dimensional, noisy, heterogeneous, and incomplete. These features make model calibration, simulation, and uncertainty quantification increasingly difficult. As a result, artificial intelligence (AI) and machine learning (ML) are playing a growing role in mathematical biology, not as replacements for mechanistic modeling, but as complementary tools that can support prediction, hybrid modeling, inference, and control. In this review, we organize these uses of AI from a mathematical biology perspective and examine how data-driven and mechanistic approaches interact across different modeling tasks. The selected examples span biomedical, epidemiological, ecological, evolutionary, biochemical, and population-level systems, while the review is intended to be representative rather than exhaustive. We emphasize recurring challenges that strongly affect biological credibility and practical usefulness, including interpretability, identifiability, generalization, and uncertainty quantification. Our goal is to clarify the roles AI can play in mathematical biology and to highlight the opportunities and limitations that arise when flexible learning methods are integrated with biologically grounded modeling.

Benchmarking gene expression reconstruction from single-cell latent representations
Single-cell transcriptomics is typically modeled in low-dimensional latent representations that improve the signal-to-noise ratio of the data. Such representations underpin data integration, cell state discovery, and perturbation prediction, with applications ranging from large-scale organ atlases to latent trajectory modeling. Recent virtual cell approaches further leverage these representations to predict cellular responses as distributional shifts in latent space. Each of these applications ultimately requires faithful gene expression reconstruction from latent spaces for biological interpretation, enabling gene-level analysis of predicted perturbed or batch-corrected cells. Yet representation choice is typically treated as an implementation detail rather than a primary modeling decision, with no systematic evaluation of how well latent representations support gene expression reconstruction. Here, we introduce ReconEval, a benchmark for evaluating gene expression reconstruction from single-cell latent spaces. We benchmark two classes of latent representations: end-to-end trained models such as PCA, autoencoders, and variational autoencoders, and pretrained single-cell foundation model embeddings coupled to newly trained decoders. Reconstruction is evaluated both directly and after latent-space perturbation prediction. Across perturbational and observational datasets totaling over 100 million cells, our metric suite quantifies statistical fidelity; biological signal preservation, including differential expression, coexpression, cell-cycle structure, cytokine response and pathway activity; and perturbation-specific effects. We find that autoencoders achieve the strongest stand-alone reconstruction at low dimensionality, while variational regularization does not improve generalization in reconstruction. Frozen foundation model embeddings retain recoverable gene-level information, with reconstruction quality depending strongly on decoder architecture and pretraining objective. In latent perturbation modeling, high-dimensional PCA matches foundation model embeddings, while low-dimensional AE embeddings are optimal for flow-based generative models. Overall, reconstruction depends critically on the interplay between representation and downstream model, and simpler representations can outperform complex alternatives given appropriate capacity. Our benchmark establishes reconstruction as a critical evaluation axis for single-cell foundation models. We envision it improving the biological interpretability of latent-space modeling, a prerequisite for future virtual cell models to be validated by domain experts and grounded in biology. ### Competing Interest Statement F.J.T. consults for Immunai, CytoReason, Valinor Industries, Bioturing and Phylo Inc., and has ownership interest in RN.AI Therapeutics, Dermagnostix, and Cellarity. The remaining authors declare no competing interests.

Dynamic Landscape Analysis of cell fate decisions provides predictive models of neural development from single-cell data
Building a mechanistic understanding of cell fate decisions remains a fundamental goal of developmental biology, with implications for stem cell therapies, regenerative medicine and understanding disease mechanisms. Single-cell transcriptomics provides a detailed picture of the cellular states observed during these decisions, but building dynamic and predictive models from these data remains a challenge. Here, we present dynamic landscape analysis (DLA), an integrative framework that applies dynamical systems theory to identify stable cell states, map transition pathways, and generate a predictive cell fate decision landscape from single-cell data. Applying this framework to vertebrate neural tube development revealed that progenitor specification by Sonic Hedgehog (Shh) can be captured in a landscape with an unexpected topology in which initially divergent lineages converge to the same fate through multiple distinct routes. The model accurately predicted cellular responses and cell fate allocation for unseen dynamic signalling regimes. Cross-species validation using human embryonic organoid data demonstrated conservation of this decision-making architecture. By modelling the dynamic responses that drive cell fate decisions, the DLA framework provides a quantitative and generative framework for extracting mechanistic insights from high-dimensional single-cell data.
Quantifying cell-state densities in single-cell phenotypic landscapes using Mellon
Cell-state density characterizes the distribution of cells along phenotypic landscapes and is crucial for unraveling the mechanisms that drive diverse biological processes. Here, we present Mellon, an algorithm for estimation of cell-state densities from high-dimensional representations of single-cell data. We demonstrate Mellon’s efficacy by dissecting the density landscape of differentiating systems, revealing a consistent pattern of high-density regions corresponding to major cell types intertwined with low-density, rare transitory states. We present evidence implicating enhancer priming and the activation of master regulators in emergence of these transitory states. Mellon offers the flexibility to perform temporal interpolation of time-series data, providing a detailed view of cell-state dynamics during developmental processes. Mellon facilitates density estimation across various single-cell data modalities, scaling linearly with the number of cells. Our work underscores the importance of cell-state density in understanding the differentiation processes, and the potential of Mellon to provide insights into mechanisms guiding biological trajectories.

Temporal tissue dynamics from a spatial snapshot
Physiological and pathological processes such as inflammation and cancer emerge from interactions between cells over time1. However, methods to follow cell populations over time within the native context of a human tissue are lacking because a biopsy offers only a single snapshot. Here we present one-shot tissue dynamics reconstruction (OSDR), an approach to estimate a dynamical model of cell populations based on a single tissue sample. OSDR uses spatial proteomics to learn how the composition of cellular neighbourhoods influences division rate, providing a dynamical model of cell population change over time. We apply OSDR to human breast cancer data2–4, and reconstruct two fixed points of fibroblasts and macrophage interactions5,6. These fixed points correspond to hot and cold fibrosis7, in agreement with co-culture experiments that measured these dynamics directly8. We then use OSDR to discover a pulse-generating excitable circuit of T and B cells in the tumour microenvironment, suggesting temporal flares of anticancer immune responses. Finally, we study longitudinal biopsies from a triple-negative breast cancer clinical trial3, in which OSDR predicts the collapse of the tumour cell population in responders but not in non-responders, based on early-treatment biopsies. OSDR can be applied to a wide range of spatial proteomics assays to enable analysis of tissue dynamics based on patient biopsies.

ResolVI: addressing noise and bias in spatial transcriptomics
Technologies for estimating RNA expression at high throughput, in intact tissue slices and with high spatial resolution (spatial transcriptomics) shed new light on how cells communicate and tissues function. A fundamental step in analyzing data generated by subcellular resolution spatial transcriptomics technologies is quantification, namely, segmenting the plane into regions, each approximating a cell, and then collating the molecules inside each region to estimate the cellular expression profile. Despite many advances in this area, a persistent problem is that of the incorrect assignment of molecules to cells, which limits many current applications to the level of a priori-defined cell subsets and complicates the discovery of novel cell states. Here we develop resolVI, a model that operates downstream of any segmentation algorithm to generate a probabilistic representation, correcting for the misassignment of molecules, as well as for batch effects and other nuisance factors. We demonstrate that resolVI improves our ability to distinguish between cell states, identify subtle expression changes in space and perform integrated analysis across datasets. ResolVI is available as open source software within scvi-tools.

An atlas-scale generative model for unified representation learning of bulk RNA-seq data
Public bulk RNA-seq repositories contain hundreds of thousands of samples, creating opportunities for large-scale representation learning, but integration across studies remains challenging because of heterogeneous annotations, experimental protocols, and technical variation. While pre-trained foundation models are now widely available for single-cell RNA-seq, comparable resources for bulk RNA-seq remain scarce, motivating a model that learns a unified, tissue-aware representation directly from bulk data. We trained a supervised variational autoencoder (VAE) on a compendium of 118,263 bulk RNA-seq samples that we assembled from TCGA, GTEx, and ARCHS4 and mapped to 42 tissue categories. The model classifies tissue of origin at 94.9% balanced accuracy (weighted F1 96.2%) and compresses 16,115 genes into a 121-dimensional latent space. Tissue identity is the primary organizing axis of the latent space, while source effects remain secondary. To assess the impact of data volume, we constructed training sets at three different scales (38K, 75K, and 118K samples). Our results demonstrated that reconstruction fidelity improved incrementally with each expansion of the dataset, but with diminishing returns. We validated the model on an independent cohort of 734 paediatric tumour samples from TARGET, achieving 84.6% agreement with the expected tissue of origin. The trained model and code are available at GitHub ([https://github.com/BIMSBbioinfo/flexynesis\_tissue\_vae_manuscript][1]) with an interactive web application. ### Competing Interest Statement The authors have declared no competing interest. [1]: https://github.com/BIMSBbioinfo/flexynesis_tissue_vae_manuscript

Spatial Transcriptomics As Rasterized Image Tensors (STARIT) characterizes cell states with subcellular molecular heterogeneity
Imaging-based spatially resolved transcriptomics (imSRT) technologies provide high-throughput molecular-resolution spatial characterization of genes within cells. Conventional analysis methods to identify cell-types and states in imSRT data rely on gene count matrices derived from tallying the number of mRNA molecules detected for each gene per segmented cell, thereby overlooking subcellular heterogeneity that can be useful in defining cell states. To take advantage of the molecular-resolution information in imSRT data and potentially identify cell-states based on subcellular heterogeneity, we developed STARIT (Spatial Transcriptomics As Rasterized Image Tensors). STARIT converts transcripts within segmented cells in imSRT data into an image-based tensor representation that can be combined with deep learning computer vision models for downstream analysis. Using simulated and real imSRT data, we demonstrate that STARIT distinguishes transcriptionally distinct cell-types and further separates cell states based on subcellular transcript localization, which conventional gene count analysis fails to capture. By providing a standardized framework to encode subcellular molecular information in imSRT data, STARIT will enable deeper insights into subcellular heterogeneity and enhance the identification and characterization of cell-types and states that are overlooked by gene count representations.
Selecting synthetic data for successful simulation-based transfer learning in dynamical biological systems
Accurate prediction of the temporal dynamics of biological systems is crucial for informing timely and effective interventions, e.g., in ecological or epidemiological contexts, or for treatment adjustments in therapy. While machine learning has proven its capabilities in generalizing the underlying non-linear dynamics of such systems, unlocking its predictive power is often restrained by the limited availability of large, curated datasets. To supplement real-world data, informing machine learning by transfer learning with synthetic data derived from simulations using ordinary differential equations has emerged as a promising solution. However, the success of this approach highly depends on the designed characteristics of the synthetic data.

Comparing phenotypic manifolds with Kompot: Detecting differential abundance and gene expression at single-cell resolution
Single-cell studies are frequently designed to compare across conditions such as health and disease. However, existing computational approaches typically rely on grouping cells into discrete populations before making comparisons, which can limit resolution for detecting state-dependent changes. Here, we introduce Kompot, a statistical framework for comparative analysis of multi-condition single-cell data. Kompot quantifies both differential abundance, capturing how cells redistribute across the phenotypic space, and differential expression, identifying condition-specific transcriptional changes that may be localized, heterogeneous, or oppositely regulated across states. By modeling cell density and gene expression as continuous functions over a shared cell-state representation, Kompot enables single-cell–resolution inference with principled uncertainty estimates, without requiring predefined clusters or cell types. Applying Kompot to aging murine bone marrow, we identified a continuum of shifts in hematopoietic stem cell and mature cell states, transcriptional remodeling of monocytes independent of compositional changes, and divergent regulation of oxidative stress response genes across cell types. We demonstrate the utility of Kompot in disease settings by identifying cell-state and gene expression changes associated with improved efficacy of combinatorial immunotherapy in melanoma. Additionally, Kompot enables multi-sample comparative analysis by accounting for sample-to-sample heterogeneity. By capturing both global and cell-state–specific effects of perturbation, the Kompot framework is broadly applicable to dissecting condition-specific effects in complex single-cell landscapes. ### Competing Interest Statement The authors have declared no competing interest. National Institutes of Health, R35GM147125, R01CA292932, T32GM136534, S10OD028685 The Mark Foundation for Cancer Research, https://ror.org/00v7th354, Endeavor Award Edward P. Evans Foundation, https://ror.org/03h22gm35, Discovery Research Grant Brotman Baty Institute, https://ror.org/03jxvbk42, Pilot Award

A unified network systems approach uncovers a core program underlying T follicular helper cell differentiation
Characterizing multi-scale processes underlying immune-state progression is critical for defining their function. This becomes pertinent for functionally diverse and plastic immune cells, such as T follicular helper (Tfh) cells. Here, we adopt a multi-scale network-systems approach that incorporates both regulatory and protein-protein interactions. This approach integrates diverse data types, captures regulation across levels of immune system organization, and recapitulates known Tfh differentiation drivers. Further, we present CoreNet, a core Tfh gene set that is conserved between humans and mice, across tissue types and disease contexts, and is consistent across data modalities. Using CoreNet, we implicate NR3C1 and interleukin (IL)-12 in the regulation of Tfh differentiation. Notably, IL-12 is permissive for differentiation of Tfh precursors but blocks differentiation into germinal center Tfh cells. Overall, this work elucidates networks with unexplored roles governing Tfh differentiation across species and tissues, while providing a generalizable framework. CoreNet is accessible through an interactive web server: https://pitt-csi.shinyapps.io/tfhcorenet/.
GEMOT: Towards Mechanistic World Models for Biology
Scientific discovery seeks mechanisms that explain observations and predict beyond the measurements that produced them. Whereas large language models (LLMs) encode knowledge implicitly, mechanistic world models organise it as a parsimonious set of explicit, modular, reusable mechanisms whose predictions can be scored against data.
An Integrated Computational Model of the Bone Microenvironment in Bone-Metastatic Prostate Cancer
Abstract Bone metastasis will impact most men with advanced prostate cancer. The vicious cycle of bone degradation and formation driven by metastatic prostate cells in bone yields factors that drive cancer growth. Mechanistic insights into this vicious cycle have suggested new therapeutic opportunities, but complex temporal and cellular interactions in the bone microenvironment make drug development challenging. We have integrated biologic and computational approaches to generate a hybrid cellular automata model of normal bone matrix homeostasis and the prostate cancer-bone microenvironment. The model accurately reproduces the basic multicellular unit bone coupling process, such that introduction of a single prostate cancer cell yields a vicious cycle similar in cellular composition and pathophysiology to models of prostate-to-bone metastasis. Notably, the model revealed distinct phases of osteolytic and osteogenic activity, a critical role for mesenchymal stromal cells in osteogenesis, and temporal changes in cellular composition. To evaluate the robustness of the model, we assessed the effect of established bisphosphonate and anti-RANKL therapies on bone metastases. At approximately 100% efficacy, bisphosphonates inhibited cancer progression while, in contrast with clinical observations in humans, anti-RANKL therapy fully eradicated metastases. Reducing anti-RANKL yielded clinically similar results, suggesting that better targeting or dosing could improve patient survival. Our work establishes a computational model that can be tailored for rapid assessment of experimental therapies and delivery of precision medicine to patients with prostate cancer with bone metastases. Cancer Res; 74(9); 2391–401. ©2014 AACR.
Modeling Osteocyte Network Formation: Healthy and Cancerous Environments
Advanced cancers, such as prostate and breast cancers, commonly metastasize to bone. In the bone matrix, dendritic osteocytes form a spatial network allowing communication between osteocytes and the osteoblasts located on the bone surface. This communication network facilitates coordinated bone remodelling. In the presence of a cancerous microenvironment, the morphology of this network changes. Commonly osteocytes appear to be either overdifferentiated (i.e., there are more dendrites than healthy bone) or underdeveloped(i.e., dendrites do not fully form). In addition to structural changes, histological sections from metastatic breast cancer xenografted mice show that number of osteocytes per unit area is different between healthy bone and cancerous bone. We present a stochastic agent-based model for bone formation incorporating osteoblasts and osteocytes that allows us to probe both network structure and density of osteocytes in bone. Our model both allows for the simulation of our spatial network model and analysis of mean-field equations in the form of integro-partial differential equations. We considered variations of our model to study specific physiological hypotheses related to osteoblast differentiation; for example predicting how changing biological parameters, such as rates of bone secretion, rates of cancer formation and rates of osteoblast differentiation can allow for qualitatively different network morphologies. We then used our model to explore how commonly applied therapies such as bisphosphonates (e.g. zoledronic acid) impact osteocyte network formation.
