A Unified Geometric Framework for Developmental Analysis of Spatial Transcriptomic Data

Mary Chriselda Antony Oliver; Kaitlyn Hohmeier; Tuyen Tran; Alejandra Castillo; Caroline Moosmüller; Shiying Li

Source-linked AI summary

A Unified Geometric Framework for Developmental Analysis of Spatial Transcriptomic Data

Mary Chriselda Antony Oliver, Kaitlyn Hohmeier, Tuyen Tran, Alejandra Castillo, Caroline Moosmüller, Shiying Li

arXiv:2608.15306v1Aug 15, 2026stat.MLcs.LGmath.MG

TL;DR

Spatial transcriptomics makes developmental dynamics difficult to reconstruct from unaligned samples, so this paper uses GW geometry and curvature to model network evolution, recovering observed developmental trends.

Problem
Methods must better represent, compare, and integrate temporal and spatial dependencies in biological data across scales, conditions, and sample types.
Method
The framework represents spatial transcriptomic stages as measure hypernetworks and uses GW geometry, geodesics, and Ollivier–Ricci curvature to analyze network evolution.
Results
Curvature distributions along GW geodesics recover developmental trends observed in real data and can distinguish networks across spatial locations and time points.
Takeaways & Limitations
Agreement among empirical patterns, GW interpolations, curvature signatures, and COOT embeddings indicates that OT geometry captures biologically meaningful information across scales.
Takeaways & Limitations
Gromov–Wasserstein optimization remains computationally challenging for large-scale transcriptomic datasets.

Abstract

from arXiv · show

High-throughput single-cell and spatial transcriptomic technologies provide high-resolution snapshots of heterogeneous cellular states, but their destructive nature prevents repeated measurements of the same cells over time. Consequently, temporal and spatial dynamics must be inferred from independently sampled, unaligned cell populations, making it challenging to reconstruct developmental trajectories. Optimal transport (OT) offers a geometric framework for aligning cell populations and inferring developmental trajectories, but many existing approaches focus on modeling the evolution of distributions of cells in gene expression space rather than the relational structure encoded by gene expression networks. To address this limitation, we introduce a geometric framework for analyzing the spatiotemporal evolution of gene expression networks through embeddings in Gromov--Wasserstein (GW) space. By representing each developmental stage as a graph combining gene expression and spatial proximity, our approach enables comparisons of network structure across time, continuous interpolation between developmental stages via GW geodesics, and quantification of network-level changes using Ollivier-Ricci curvature. We evaluate our framework on a spatiotemporal transcriptomic \textit{Drosophila} dataset and show that GW geodesic interpolations reproduce main trends in curvature dynamics observed in empirical gene expression networks. Agreement with higher-order Co-Optimal Transport (COOT) distances, which jointly represent spatial and temporal information, further validates the framework and suggests that hypernetwork representations successfully record salient biological changes across time. In general, our approach provides a unified geometric approach to study dynamically evolving biological networks.

1 Introduction

Spatial transcriptomics provides rich molecular snapshots across space and time, but destructive sampling, correlations, and alignment requirements complicate reconstruction of developmental dynamics. The paper addresses these challenges with a unified geometric framework combining graph curvature, Gromov–Wasserstein geometry, and co-optimal transport.

Challenges: Inferring developmental dynamics remains difficult because measurements are destructive, spatially or temporally correlated, and require alignment across individuals or biological scales.These limitations hinder full use of temporal information and understanding of cell types and interactions.
Motivation: Spatial transcriptomic measurements reveal molecular and biological changes across space and time that static measurements cannot capture.They support tissue-development reconstruction, disease-elimination analyses, and characterization of structural variation such as within tumors.
Geometric approach: The framework uses optimal transport, especially Gromov–Wasserstein distance, to compare spatial-temporal biological systems while respecting intrinsic geometry and differing network sizes.GW is suited to gene co-expression networks across developmental stages because graphs can have different numbers of nodes.
Main contributions: Edgewise Ollivier-Ricci curvature distributions provide discriminative signatures for comparing and classifying gene expression networks across temporal stages.The framework evaluates developmental changes using stage-wise OR curvatures in a Drosophila spatial transcriptomics dataset.
Main contributions: GW geodesic interpolation models developmental progression between initial and terminal stages, while hypernetworks and COOT capture multilayer spatial-temporal organization across Drosophila development.Together, these tools support developmental-stage characterization, trajectory analysis, and higher-order network comparison.

2 Background

The framework represents spatiotemporal gene-expression data as metric-measure graphs and hypernetworks, enabling structural comparisons in Gromov–Wasserstein space. It uses GW geodesics and barycenters to interpolate between developmental stages while incorporating spatial and higher-order relationships.

Graph and metric-measure representations: Gene co-expression data are modeled as finite graphs or metric measure spaces, with genes as nodes and thresholded pairwise associations as edges.Associations commonly derive from correlations of gene expression across samples; weak associations can be excluded by thresholding.
Gromov–Wasserstein geometry: The p-Gromov-Wasserstein distance compares metric measure spaces up to isomorphism, and the p = 2 case provides a geodesic structure for network interpolation.GW geodesics were theoretically characterized, while the implementation focuses on p = 2.
Measure hypernetworks: Measure hypernetworks extend the representation to spatial-temporal data and support comparisons across multiple structural layers, including rectangular relationships.Unlike square kernels indexed by one space, rectangular kernels relate distinct spaces such as cells and genes.
Measure hypernetworks: At a fixed developmental stage, a spatial transcriptomic dataset is encoded as a hypernetwork whose cell–gene function records each gene’s expression at each spatial location.The two underlying spaces are spatial locations or cells and genes, each equipped with a probability measure.

3 Methodology and analysis · 3.1 Curvature as geometric signatures

The framework represents each developmental stage as a weighted gene co-expression network, analyzes curvature in graph and metric-measure settings, and approximates interpolated spaces using GW barycenters. The graph approximation separates topology construction from edge-weight assignment based on unweighted and weighted shortest-path information.

3.1 Curvature as geometric signatures: Each developmental stage is represented by an undirected weighted gene co-expression network whose nodes are selected genes and whose edges encode retained nonnegative co-expression strengths.
3.1 Curvature as geometric signatures: In the graph setting, edge-level Ollivier-Ricci curvature is computed on retained co-expression edges using the graph metric.
3.1 Curvature as geometric signatures: For trajectory analysis via GW geodesics, each observed graph is treated as a finite metric-measure space with graph shortest-path distances and a uniform probability measure on nodes.Interpolated spaces are finite metric-measure spaces rather than observed gene co-expression networks, so their local measures are constructed from metric structure.
3.1 Curvature as geometric signatures: Interpolated spaces are converted into graph approximations whose common-gene nodes use GW barycenter topology and weights derived from the corresponding weighted GW barycenter.The resulting graph-based construction supports curvature analysis of the interpolated spaces.
3.1.1 Graph approximations of interpolated mm-spaces: The graph approximation uses two steps: GW barycenter distances from unweighted shortest-path matrices determine connectivity, while weighted shortest-path barycenters assign edge weights.
3.1.1 Graph approximations of interpolated mm-spaces: Two vertices are connected when their unweighted barycenter distance is sufficiently close to one, with the prescribed tolerance δ∈(0, 1) controlling this criterion.
3.1.1 Graph approximations of interpolated mm-spaces: For disconnected graphs, infinite unweighted distances are replaced by a large finite γ≫0, while missing interactions in weighted matrices receive value zero before applying the two-step construction.

3.2 Inference with curvature and Gromov-Wasserstein geodesics

The framework evaluates whether GW geodesics capture developmental progression by matching real-graph curvature distributions to interpolated points and measuring recovered temporal ordering. It also classifies spatial slices against stage-specific curvature distributions using nearest-neighbor discrepancies.

Geodesic-based temporal classification: Three discrepancy methods compare curvature distributions: scalar mean or median, raw 1-Wasserstein distance, and symmetrized KL divergence.The Wasserstein distance captures location, spread, and shape, whereas standardized KDEs make the KL divergence focus on shape differences and symmetrization removes reference-distribution dependence.
Geodesic-based temporal classification: Each real graph is assigned to the closest of 40 points along the discretized GW geodesic by minimizing a curvature-distribution discrepancy.The predicted geodesic index is the argmin over discrepancy values computed between the real graph and every interpolation point.
Geodesic-based temporal classification: Temporal-order recovery is evaluated by Kendall’s τ and Spearman’s ρ between the true order of five graphs and their predicted order.Values near +1 indicate that a discrepancy measure recovers the correct developmental ordering.
Slice-to-stage curvature classification: A second task assigns each of 50 spatial-slice observations to one of five developmental stages by nearest-neighbor comparison with stage-specific reference curvature distributions.Each stage contributes 10 observations, and uniform random assignment has expected accuracy 0.2.

3.3 Trend analysis with OR curvature

The section evaluates whether GW trajectories capture developmental changes in average Ollivier-Ricci curvature. It compares endpoint-geodesic predictions with observed stages and assesses directional consistency against a piecewise path through all five stages.

Stagewise predictive accuracy: Stagewise evaluation compares observed weighted graphs across five developmental stages with predictions from the GW geodesic connecting the initial embryonic and terminal larval stages.Developmental pseudo-times are derived from reported developmental times.
Stagewise predictive accuracy: The comparison uses average edgewise OR-curvature distributions, measuring relative l1 error and the R2 score of GW-geodesic-based predictions.The curvature predictions are evaluated against the corresponding observed curvature distributions.
Trajectory-level trend consistency: Trajectory-level analysis compares average curvatures along the endpoint geodesic from G1 to G5 with a piecewise GW path concatenating interstage geodesics through all five stages.The interstage paths G1 → G2, G2 → G3, G3 → G4, and G4 → G5 are rescaled using reported developmental times.
Trajectory-level trend consistency: Directional agreement between curvature trends is quantified using sign alignment, based on the signs of successive trend differences.The relaxed score introduces 𝜖 to account for small oscillations, with 𝜖 chosen as a percentage of maximum variation.

3.4 Trend analysis with dynamic OR curvature

This section compares edge-wise distance-scaled dynamic Ollivier–Ricci curvature distributions across five developmental-stage networks over diffusion scales, using Wasserstein discrepancies to quantify stage differences. The analysis summarizes distributional changes across scales while treating comparisons as descriptive rather than independent biological replication.

Method: Edge-wise distance-scaled dynamic OR curvature distributions are compared across five observed stage-specific gene co-expression networks over diffusion parameters τ.The analysis uses 19 τ values equally spaced on the log10(τ) scale from 0.01 to 10, comparing empirical distributions because retained edge sets differ.
Method: Distance scaling preserves sensitivity to network changes encoded in graph-weight scale, unlike original curvature under uniform rescaling.Uniform rescaling leaves original curvature unchanged but changes distance-scaled curvature.
Method: Each stage-scale distribution is summarized by its median, interquartile range, and ECDF, while pairwise differences use 1-Wasserstein discrepancies before and after median-centering and interquartile-range scaling.The centered/scaled comparison measures differences remaining after removing location and scale effects.
Across diffusion scales: Wasserstein discrepancies are averaged over log10(τ) from τmin = 0.01 to τmax = 10 using the trapezoidal rule.For consecutive stages, the analysis also records which transition has the largest discrepancy at each sampled τ.
Interpretation: The comparisons are descriptive and do not treat edge-level curvature values as independent biological replicates.The procedure is summarized in Algorithm 3 and applied to the Drosophila data in Section 4.4.

3.5 Quantifying developmental changes via COOT distance

The section represents each developmental stage as a gene–cell interaction matrix and compares stages using COOT distances that jointly align spatial and molecular structures. Classical multidimensional scaling embeds the resulting pairwise developmental-stage dissimilarities in one dimension while preserving their induced geometry.

Stage representation: Each developmental stage is represented by a rectangular gene–cell interaction matrix whose entries encode gene expression levels at spatial locations.The matrix H_s has dimensions n_c × n_g, with n_c spatial locations and n_g genes.
COOT distance: COOT defines developmental-stage dissimilarity by jointly optimizing couplings over cell and gene spaces using block coordinate descent.This produces the discrete hypernetwork p-distance D_ss′ = COOT_p(H_s, H_s′).
COOT distance: The resulting transport-based discrepancy captures differences in joint spatial and molecular relational organization rather than entrywise differences.Pairwise dissimilarities form a symmetric matrix D ∈ R^(S×S).
Geometric embedding: Classical multidimensional scaling produces one-dimensional embeddings of developmental stages that preserve the COOT-induced geometry.The embeddings are denoted y_1, …, y_S ∈ R.

4 Application to Drosophila data

The Drosophila analysis constructs stage-specific gene co-expression graphs from five spatial transcriptomic developmental stages and evaluates developmental trajectories using GW geodesics, curvature, and COOT distances. These analyses recover stage ordering, classify slices, and identify L3 as structurally divergent from earlier stages.

Data and preprocessing: The dataset spans five Drosophila developmental stages with gene expression measurements and spatial coordinates, enabling joint analysis of transcriptional variation and spatial organization.Each tissue section has a unique slice ID providing a consistent spatial reference frame.
Data and preprocessing: The analysis selects the union of stage-specific top 50 highly variable genes, computes pairwise Pearson correlations, and retains correlations with |a_i_j| ≥ 0.2 to form weighted sparse graphs.The graph retains the top 5% of correlations by absolute value, with μ=0.2 fixed throughout.
Trajectory prediction: GW-geodesic curvature distributions recover the developmental trend, with KL divergence achieving Kendall τ=0.9016 and Spearman ρ=0.9528.The mean-curvature discrepancy measure performs worse than the Wasserstein and KL divergence methods.
Stage classification: The Wasserstein method achieves perfect classification accuracy for assigning all spatial slices to their correct developmental stages.Per-slice discrepancy heatmaps visualize both classifications and their relative closeness to reference stages.
Trajectory validation: GW-geodesic predictions reduce square error by 37.5% relative to the constant across-stage mean predictor.The endpoint geodesic also shows directional consistency with the piecewise trajectory, with S_ε=0.82 versus chance baseline bC_ε=0.73 when ε=0.008.
Network-level developmental changes: COOT distances remain approximately 0.006–0.010 among E14–16, E16–18, L1, and L2 but increase to approximately 0.024–0.028 for comparisons involving L3.At diffusion scales τ=1 through τ=10, every raw 1-Wasserstein comparison involving L3 exceeds comparisons among the earlier stages, identifying L3 as structurally divergent.

5 Discussion

The study presents GW geometry and Ollivier–Ricci curvature as a framework for analyzing developmental changes in spatial transcriptomic co-expression networks. It reports biologically informative multiscale network comparisons while identifying inference and computational limitations and proposing extensions.

Contributions: The framework represents developmental gene-expression patterns as spatially and temporally varying co-expression networks analyzed with GW geometry and Ollivier–Ricci curvature.It was developed for a spatial transcriptomic Drosophila melanogaster dataset.
Main conclusions: Ollivier–Ricci curvature distributions distinguish gene-expression networks across spatial locations and developmental time points, providing a scale-dependent view of developmental differences.This supports prior work on recovering biological information from gene-expression networks while incorporating temporal and spatial variation.
Limitations: Correlation-thresholded co-expression networks are sensitive to correlation thresholds, gene selection, and treatment of negative correlations, while GW optimization remains computationally challenging for large datasets.The study calls for alternative network-inference strategies and biologically informed interaction networks.
Future directions: Future work could add biologically informed edge weights, extend Ollivier–Ricci curvature to hypergraphs, and use unbalanced GW or COOT variants to accommodate changing cell populations and gene regulation.Proposed edge sources include regulatory interactions, ligand–receptor communication, and multi-omics data.
Broader significance: GW-based geometric methods extend transcriptomic analysis beyond conventional distribution alignment by connecting empirical patterns, GW interpolations, curvature signatures, and COOT embeddings.Their agreement indicates that OT geometry encodes biologically meaningful information across multiple scales for developmental and cellular processes.

A Supplementary Figures and Tables

The supplementary figures examine gene-correlation structure, classification rankings and discrepancies, and curvature behavior across developmental stages. They show largely stable gene relationships except for marked changes at L3 and provide curvature visualizations under specified thresholds and scaling parameters.

Gene correlations: Gene-correlation heatmaps across five Drosophila stages show visually similar correlation values and gene-relationship structure, except for large changes at L3.The stages are E14-16, E16-18, L1, L2, and L3, using common genes across all stages.
Classification results: The supplementary results include rank-ordering results for the mean method.
Classification results: Discrepancy plots compare mean and median classification methods, marking correct classifications with stars and incorrect classifications with x symbols.Black boxes identify the corresponding correctly or incorrectly classified class.
Curvature analyses: Average Ollivier-Ricci curvatures are plotted for geodesics between different fruit-fly developmental stages using a gene-network threshold of μ=0.2.
Curvature analyses: Histograms visualize distance-scaled dynamic OR curvature values κ_ds for L1, L2, and L3 at τ=0.1 and τ=1.Kernel density estimates are included for visualization.

B Supplementary Details · B.1 Chance baseline for relaxed sign alignment

The section defines consecutive increments for two scalar trends and derives a relaxed sign-alignment chance baseline under independent, symmetric sign assumptions. The derivation identifies disagreement as opposite signs among increments whose magnitudes meet the threshold 𝜖, and notes a corresponding empirical baseline.

B.1 Chance baseline for relaxed sign alignment: The increments are defined as Δf_i = f(t_i+1) − f(t_i) and Δg_i = g(t_i+1) − g(t_i) for i = 1, . . . , n−1.
B.1 Chance baseline for relaxed sign alignment: Lemma 1 assumes independent random increments from two scalar trends, with signs equally likely positive or negative conditional on magnitude at least 𝜖.
B.1 Chance baseline for relaxed sign alignment: Taking the complement of the relaxed-disagreement probability yields the chance baseline for the relaxed sign-alignment score.
B.1 Chance baseline for relaxed sign alignment: Under the relaxed criterion, disagreement occurs only when both increments have magnitude at least 𝜖 and opposite signs.
B.1 Chance baseline for relaxed sign alignment: Independence is used to obtain the probability of relaxed disagreement from the joint magnitude and sign conditions.
B.1 Chance baseline for relaxed sign alignment: Conditional on both increments meeting the magnitude threshold, independent equally likely signs make opposite signs occur with probability 1/2.
B.1 Chance baseline for relaxed sign alignment: The section also defines a corresponding empirical relaxed chance baseline.

B.2 Additional Data Preprocessing

For classification, the preprocessing uses a bootstrapped approach that generates multiple correlation matrices per developmental stage from randomized HGV gene instances across slices. It also uses a union of the same genes across slices for each sample.

B.2 Additional Data Preprocessing: The preprocessing generates 10 correlation matrices for each developmental stage using different random instances of HGV genes across slices.This bootstrapped construction is motivated by prior work cited as.
B.2 Additional Data Preprocessing: For each sample, the method considers the union of the same genes across different slices, analogous to Section 4.1.