Learning developmental dynamics from lineage context

Gokul Gowri; Nic Fishman; William N. Colgan; Minming Wang; Yusuke Kijima; Sachin A. Kumar; Aleksandr Zimin; Luke W. Koblan

q-bio.QMOct 29, 2026

Learning developmental dynamics from lineage context

Authors: Gokul Gowri, Nic Fishman, William N. Colgan, Minming Wang, Yusuke Kijima, Sachin A. Kumar, Aleksandr Zimin, Luke W. Koblan,

Organizations: Department of Electrical Engineering and Computer Science, Massachusetts Institute of Technology, Cambridge, MA, 02139, USA · Whitehead Institute for Biomedical Research, Cambridge, MA, 02142, USA · Howard Hughes Medical Institute, Chevy Chase, MD, 20815, USA · David H. Koch Center for Integrative Cancer Research, Massachusetts Institute of Technology, Cambridge, MA, 02139, USA · Department of Biology, Massachusetts Institute of Technology, Cambridge, MA, 02139, USA · Department of Statistics, Harvard University, Cambridge, MA, USA · Department of Biomedical Informatics, Harvard Medical School, Boston, MA, USA · Department of Genetics, Yale Stem Cell Center, Yale School of Medicine, New Haven, CT, 06519, USA · Max Planck Institute of Molecular Genetics, Max Planck Society, Berlin, 14195, Germany · Department of Pediatric Oncology, Dana-Farber Cancer Institute and Division of Hematology/Oncology, Boston Children’s Hospital, Harvard Medical School, Boston, MA, USA · Broad Institute of MIT and Harvard, Cambridge, MA, USA · Kempner Institute for the Study of Natural and Artificial Intelligence at Harvard University, Allston, MA, USA · Harvard Data Science Initiative, Cambridge, MA, USA · Department of Mathematics, Massachusetts Institute of Technology, Cambridge, MA, 02139, USA

Abstract

A central goal of biology is to predict the dynamics of multicellular systems. Lineage context (the states of a cell's sisters, cousins, etc.) could inform these predictions, but how best to leverage it remains unclear. Here we introduce lineage inpainting, a generative modeling approach that uses lineage context with tree-structured transformers and stochastic interpolants to model cell state dynamics across lineage trees. Applied to a lineage-traced transcriptional atlas of mouse embryogenesis, it recovers dynamics not resolved by gene expression alone, including multiple transcriptional convergence events where distinct lineages are drawn to a common state. Despite being fit on unperturbed development, lineage inpainting can generate hypotheses about perturbation responses by simulating counterfactuals. These simulations agree in many cases with observed developmental phenotypes, including a single-cell measurement of the Tbx5 knockout response. Our results show that lineage information provides a rich source of context for learning multicellular dynamics.

Figures & tables

Figure 1: Reconstructing developmental dynamics using lineage inpainting. \captionmore

Figure 2: Lineage inpainting leverages tree topology and lineage context to reconstruct synthetic developmental processes. \captionmore

Figure 3: Summarizing cell state transitions inferred by lineage inpainting. \captionmore

Figure 4: Lineage inpainting captures transcriptional convergence dynamics. \captionmore

Figure 5: Lineage inpainting suggests developmental effects of genetic perturbations. \captionmore

\widehat{\mathcal{X}}_{i}=I_{\theta}(\mathcal{T}_{i},A_{i},x_{i,A_{i}}),\hskip 18.49988pti\in\mathcal{I}\cup\mathcal{I}^{\prime}.

Algorithm 1 Imputation alignment loss

Cell type	Hox genes
Hindbrain (r2–r4)	Hoxa2 , Hoxb2 , Hoxb1
Cranial motor neuron	Hoxb1
Cranial neural crest	Hoxa2
Pharyngeal arch mesenchyme	Hoxa2 , Hoxb2 , Hoxa3 , Hoxd3 , Hoxd4
Hindbrain (r5–r8)	Hoxa3 , Hoxd3 , Hoxb4 , Hoxd4
Spinal cord (cervical)	Hoxb4 , Hoxd4

Table 1: Cell types with stable Hox gene expression between E9.0 and E9.5

Figure 6: Additional results on synthetic and semi-synthetic benchmarks. (a) Fate accuracy of reconstructions across coherent and crossing dynamics from Fig. 2 , decomposed by transition. (b) Type accuracy of model predictions from random cell type motif benchmark of Fig. 2 , including additional baseline methods. WaddingtonOT and moslin are non-generative methods, here we evaluate using the discrete coupling. (c) Same as (b), but for the cell-type context motifs. (d) Comparison of fate accuracy metrics with distributional metrics, on the random cell type motif benchmarks. MMD refers to maximum mean discrepancy. Colors match panel (c).

Figure 7: Inspecting the learned attention heads. (a) Learned patristic distance bias scale ( $\gamma$ ), across each head in the sampler of a lineage inpainting model fit to Colgan et al. [36] . One dot per each of four heads across four layers. (b) Fraction of attention paid to observed (input condition) nodes, across interpolant times, across each head. Hue matches layers from panel (a). (c) For an example input subtree, the normalized attention matrix of attention head 2 from layer 2, at different interpolant times. Columns represent key nodes, and have ticks hued by observed vs. masked.

Figure 8: Model evaluation on held-out lineage trees. The metrics of Fig. 3 b , computed over the 11 clone trees that were held out. Lower is better for every metric. Filled circles are models trained on the remaining 45 trees, while open circles are the models from the main text. Points show the mean, and bars show the 95% interval over subtrees (clonal), stage pairs (marginal), stage transitions (linkage) or (cell type, Hox gene) pairs (Hox-code).

Figure 9: Comparing lineage inpainting models with Qiu et al. [9] manually curated transitions. (a) Agreement with Qiu et al. [9] transition map, at the aggregate level of edges between cell types. Match denotes an edge in both graphs, including identity edges between the same cell type label at different times, reverse denotes an edge in Qiu et al. [9] which lineage inpainting infers the opposite direction, boundary denotes an edge between cell types which are labeled as spatially adjacent in Qiu et al. [9] , or nodes which share an ancestor in Qiu et al. [9] , unmapped edges involve a cell type with no corresponding annotation in Qiu et al. [9] . (b) Same as (a), but now measured at the single cell transition level rather than the aggregate cell type level. (c) - (f) List of top 10 most frequent transitions from each non-match category.

Figure 10: Additional results on transcriptional convergence analysis. (a) Table of top 10 and bottom 10 cell types as ranked by convergence score. Convergence and OT discrepancy refer to normalized scores from Fig. 4 b . Leaves refers to number of observed cells of each type from which scores are computed. (b) Same as Fig. 4 b , except where scores have not been normalized by local nearest-neighbor distances. (c) Mean expression (per cell type) of each of the top four genes per transcriptional program. LP refers to lateral plate direction, paraxial refers to paraxial mesoderm direction, post. refers to posterior direction, ant. refers to anterior direction. EC = endothelial cell, SCI = spinal cord interneuron, SMN = somatic motor neuron. (d) Fraction of divergent genes considered erased, per cell type, across different divergence (E8.5) and memory (E9.5) cutoffs.

Figure 11: Additional results on perturbation effect comparisons with MGD annotations. (a) Across all MGD annotated (gene, cell type) pairs, the predicted depletion of the cell type under gene knockout against the mean expression of the gene within the cell type. (b) The top six transcription factors with low cell type mean expression, and high predicted depletion. Chosen by difference between rank across 34 cell types with corresponding MGD annotations. (c) Representative trees with inferred expression along edges for cell type gene pairs with low expression and high predicted depletion.

Figure 12: Additional results on Tbx5 mutant profiling. (a) UMAP of cells, prior to quality filtering, based on SNP profile, hued by inferred embryo identity. (b) Number of cells, counts per cell, and genes per cell, across embryos. UMI = unique molecular identifier. (c) Correlation of measured nearest-neighbor distortion (per reference E9.5 cell) and Tbx5 expression. (d) Comparison of model predicted distortion agreement and reference cell Tbx5 expression correlation, across embryos. (e) Energy distance of Tbx5 crispant cells to E9.0 reference and E9.5 reference cells of corresponding coarse cell type label. (f) UMAP of neuroectoderm cells from E9.5 reference and Tbx5 crispant. Neurons hued, all other cells grey. (g) Same as (f), hued by dataset rather than cell type. (h) Fraction of neuroectoderm labeled as a neuron, across reference and crispant embryos.

Figure 13: The higher-order components of the alignment loss tighten the lineage structure. The alternating VP-SDE model trained with the state marginal alone, with the path added, with the common-ancestor triple added, and with all three loss components together.

Figure 14: The additive baseline improves the coupling metrics. Using the full grid of models we compare the plain tree transformer and the tree transformer using the additive baseline in Sec. 4.1.2 . The additive baseline improves essentially all metrics across all models.

Figure 15: Comparing the imputer and sampler models on ancestral inference metrics. Using the full grid of models we compute all our metrics for both the imputer and the sampler, plotting them against each other. Points below the diagonal demonstrate better performance as a sampler than as an imputer. We use markers to indicate the three different sampler models we test: the one-step model, the CFM, and VP-SDE.

Figure 16: Iterative samplers sharpen the distributional metrics over the one-step model. Each row is an imputer and a pairing; open circles show the one-step model, and colored circles show flow matching (CFM) and variance-preserving diffusion (VP-SDE). Left is better in every panel.

Figure 17: Comparing model performance against different training setups and gradient flows. Using the models from the grid we can compare three methods to leverage the imputer in the sampling model. We compare training on the raw imputed trees, the nearest neighbors, and a training procedure that alternates between the imputer and the nearest neighbors. The raw and nearest-neighbor settings are shown with and without gradients flowing from the sampler to the imputer; the alternating setting retains this gradient.

Figure 18: Iterative models with the additive baseline infer similar cell-type transition graphs. We compute the cell-type transition graph as in the main text figures for each model. We then convert these graphs into a vector, where each dimension of the vector represents the edge weight of an edge in the underlying graphs. We compute the cosine similarity between these vectors to represent how similar the induced cell-type transitions are across models. Missing edges have zero weight, and vectors use the union of time-resolved edges. The color scale transforms cosine similarity $c$ as $-\log_{10}\max(1-c,10^{-6})$ .

Figure 19: A path table can leave the whole-tree law partially identified. Panel a contrasts the independent fork law $P$ with the Bernstein triple $Q$ , and panel b shows their masses on the eight ordered assignments. Panel c maps those masses to the product moment representation. The distributions are the same except for the triple-product moment. Panel d normalizes this single squared difference: the whole-fork comparison is $1$ , while the parent–child comparison is $0$ because every pair marginal agrees.

Figure 20: Fixed-size feature tables trade sampling variance for partial identification. (a) Purple marks show the path table, and teal marks show the three common-ancestor tables; marginal tables are redundant with these joint tables. (b–d) At $\rho=0$ , $P_{0}$ and $Q_{0}$ agree on the selected tables although their four-grandchild laws differ; panel d shows the resulting compatible set. (e) Purple polygons project the compatible sets onto two omitted product-moment combinations, with $P_{\rho}$ in black and the numerical worst case $Q_{\rho}$ in orange. (f–h) For $(P_{\rho},Q_{\rho})$ , the heatmaps show squared bias for the whole-tree energy, feature sampling variance, and their sum for the V-statistic at a fixed number $n$ of independent trees per population, averaging child orderings exactly. Colors are relative to whole-tree MSE, and the black contour marks equal MSE.

depth	nodes	orderings	placements	whole-tree coordinates	feature coordinates
2	7	8	22	56	66
3	15	128	134	1,920	402
4	31	32,768	646	$1.0\times 10^{6}$	1,938
5	63	$2.1\times 10^{9}$	2,790	$1.4\times 10^{11}$	8,370
6	127	$9.2\times 10^{18}$	11,430	$1.2\times 10^{21}$	34,290

Table 2: Exhaustive ordering enumeration requires exponentially more coordinates than feature-placement enumeration. For a complete binary tree of depth

d

with

N

nodes, the whole-tree statistic averages over the

2^{2^{d}-1}

orderings of each tree,

N

coordinates each, while the feature representation lists every ancestor–middle–descendant triple at three distinct generations and every common-ancestor triple with distinct ordered descendants, three coordinates per placement. Marginal coordinates are omitted from this representation-size comparison.

Component	Dimensions
State embedding (scVI)
Input genes	$3{,}956$
Encoder	$2$ hidden layers, width $512$
Latent width	$50$
Decoder	$2$ hidden layers, width $512$
Inputs and embeddings

Table 3: Architecture of the state embedding and lineage models.

May 27, 2026

Applications of temporal graph learning for predicting the dynamics of biological systems

Biological foundation models have shown strong performance in single-cell representation learning by applying transformer architectures directly to gene-expression matrices. However, these approaches predominantly operate in static settings and do not explicitly model the temporal evolution of developmental programs in the cell. Modeling such dynamics is important for understanding how cellular states progressively emerge, differentiate, and reorganize during development or disease progression. In this work-in-progress paper, we investigate an alternative temporal graph-based perspective in which cellular states are represented through pseudotime-resolved gene regulatory networks and modeled as evolving graph structures over persistent gene identities. Starting from single-cell transcriptomic data, we infer pseudotime trajectories, discretize cells into developmental snapshots, reconstruct one gene regulatory network per snapshot, and apply temporal graph neural networks to forecast biological states. We evaluate this framework on two publicly available mouse developmental datasets, erythroid gastrulation and pancreatic endocrinogenesis, considering three complementary tasks: gene-expression forecasting, link prediction, and out-degree centrality prediction. Our results show that graph-based models outperform well-known foundation-model such as scGPT and scFoundation, suggesting that explicitly modeling evolving regulatory structure provides useful information beyond static pretrained representations. For link prediction and centrality forecasting, temporal graph learning captures non-trivial regulatory dynamics and enables the identification of temporally important gene hubs. Overall, our findings support temporal graph learning as a promising direction for modeling dynamic biological systems and as a complementary paradigm to current foundation model approaches in single-cell biology.

Manuel Dileo, Andrea Sottoriva

Center for Computational Biology, Human Technopole Foundation, Milan, Italy

Temporal GNNs Temporal Link Prediction Single-Cell Transcriptomics

Ask PDF

Sep 27, 2026

KoopCell: Koopman-Based Generative Model for Learning Single-Cell Dynamics from Distribution Snapshots

Learning population dynamics from temporally sparse, unpaired distribution snapshots is a fundamental challenge in developmental biology. Recent approaches based on neural differential equations and flow matching can interpolate between observed population snapshots, but may struggle to extrapolate beyond the training horizon and often lack an explicit mechanism for modeling developmental branching. We propose KoopCell, a unified generative framework based on Koopman-Mori-Zwanzig theory that jointly learns representations and predictive linear latent dynamics. Theoretically, using the weak continuity equation, we derive a closed-form least-squares estimator for the Koopman generator from distribution snapshots and establish convergence guarantees under suitable assumptions. To model branching dynamics, we further develop KoopCell-M, which incorporates non-Markovian memory into the latent Koopman dynamics through a Markovian embedding. Experiments on synthetic systems and three scRNA-seq datasets demonstrate the ability of our framework to recover Koopman spectra, model branching through memory, and scale to predicting high-dimensional gene expression distributions, achieving state-of-the-art performance among the evaluated methods.

Wanfeng Lu, Yutong Zhang, Keyi Zhou +3

School of Mathematical Sciences, SCMS, SCAM, and CCSB, Fudan University, China · Research Institute of Intelligent Complex Systems, Fudan University, China · State Key Laboratory of Medical Neurobiology and MOE Frontiers Center for Brain Science, Institutes of Brain Science, Fudan University, China +1

Koopman Operator Learning Latent Dynamics Modeling Generative Modeling Single-Cell Transcriptomics

Ask PDF

Apr 27, 2026

Learning biophysical models of gene regulation with probability flow matching

Cellular differentiation is governed by gene regulatory networks, the high-dimensional stochastic biochemical systems that determine the transcriptional landscape and mediate cellular responses to signals and perturbations. Although single-cell RNA sequencing provides quantitative snapshots of the transcriptome, current methods for inferring gene-regulatory dynamics often lack mechanistic interpretability and fail to generalize to unseen conditions. Here we introduce Probability Flow Matching (PFM), a scalable framework for learning biophysically consistent stochastic processes directly from time-resolved single-cell measurements. Applying PFM to three hematopoiesis datasets, we show that models with similar interpolation accuracy can encode fundamentally different dynamics, with only biophysically consistent formulations accurately capturing mechanisms of lineage transitions, fate specification, and gene perturbation responses. We further demonstrate that PFM accommodates unbalanced populations, enabling simultaneous inference of cellular proliferation and death dynamics. Together, these results establish PFM as a flexible, scalable framework for integrating mechanistic modeling with single-cell omics.

Suryanarayana Maddu, Victor Chardès, Michael J. Shelley

Center for Computational Biology, Flatiron Institute, New York, NY, USA, 10010 · Molecular & Cellular Biology Department, Harvard University, Cambridge, MA 02138 · Courant Institute of Mathematical Sciences, New York University, New York, NY, USA, 10012

Flow Matching Single-Cell Perturbation Modeling Gene Regulatory Network Inference Single-Cell Transcriptomics

Ask PDF

Abstract

Explore similar work

Applications of temporal graph learning for predicting the dynamics of biological systems

KoopCell: Koopman-Based Generative Model for Learning Single-Cell Dynamics from Distribution Snapshots

Learning biophysical models of gene regulation with probability flow matching