Benchmarking zero-shot single-cell foundation model embeddings for cellular dynamics reconstruction

22 Reconstructing cellular trajectories from time-resolved single-cell transcriptomics is fundamental to 23 understanding processes from embryonic development to cancer progression. While single -cell 24 foundation models (scFMs) promise universal biological representations through large -scale 25 pretraining, their capacity to capture the non-linear dynamics governing cell-fate decisions remains 26 uncharacterized. Here we systematically benchmark multiple scFMs across challenging biomedical 27 scenarios involving branching lineages and continuous state transitions . By c oupling zero -shot 28 scFM embeddings with dynamic optimal transport, we evaluated their performance against a 29 traditional highly variable gene (HVG) baseline in backtracking progenitor states, interpolating 30 transition intermediates, and extrapolating future fates. We find that zero -shot scFM embeddings 31 underperform the HVG baseline across diverse biological systems, particularly in recovering the 32 distributional complexity of unobserved cells. Mechanistic analysis reveals that current scFM 33 architectures tend to over -compress subtle temporal signals, causing an artificial "linearization" of 34 branched biological structures that may obscure critical divergence points in disease progression . 35 Our findings suggest that while scFMs provide unified cell -state views, the HVG baseline remains 36 more robust for trajectory inference, identifying a fundamental "temporal -compression" bottleneck 37 that must be addressed to develop next-generation, dynamics-aware foundation models. 38 39

40 Understanding how cells change state over time in response to development, differentiation, or 41 external perturbations is a central problem in biology 1. High -throughput single -cell profiling has 42 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint enabled the characterization of molecular states across vast cellular populations, providing detailed 43 snapshots of transcriptional heterogeneity1,3–6. However, these assays are destructive, the same cell 44 cannot be observed as it changes over time or in response to perturbations. Consequently, time -45 series single-cell experiments yield a sequence of unaligned snapshots rather than direct longitudinal 46 measurements. Reconstructing continuous cellular dynamics from such snapshots therefore 47 requires computational strategies that (i) embed sparse, noisy, high -dimensional single -cell 48 transcriptomes into informative low -dimensional representations, and (ii) estimate transportation 49 maps across time to approximate population-level flows.7,8 50 Traditionally, cell embeddings for trajectory analysis are obtained by selecting highly variable genes 51 and applying dimensional -reduction methods such as principal component analysis, after which 52 lineage structure and cellular dynamics are inferred in the resulting low -dimensional space 9. In 53 parallel, single-cell foundation models like Geneformer2, scGPT10, Genecompass11, are pretrained 54 on large collections of scRNA -seq profiles to learn general -purpose representations of cells and 55 genes, which are then used in a zero -shot manner for downstream analyses. These models have 56 been reported to transfer to diverse tasks, including cell clustering, annotation, batch correction, and 57 gene-level applications such as regulatory network inference 2,10–18. Because they are exposed to a 58 wide range of tissues, conditions, and species during pretraining, their embeddings are expected to 59 capture shared structure in gene-expression space, to be robust to technical variation, and to support 60 data-efficient analysis in new settings15. If such properties extend to temporal processes, single-cell 61 foundation models might provide more informative representations for trajectory inference than 62 conventional HVG-based embeddings, by better preserving subtle transitional states and stabilizing 63 estimates of population -level flows across time. Nonetheless, existing benchmark studies have 64 primarily focused on conventional, largely static single -cell analysis tasks , such as cell clustering, 65 annotation, and batch correction, and have shown that, in zero-shot settings, single-cell foundation 66 models often provide limited or inconsistent gains over simple HVG-based baselines19–21. In contrast, 67 the performance of these models on explicitly dynamical tasks has been much less explored. Recent 68 benchmarks on perturbation prediction suggest that foundation models do not necessarily 69 outperform traditional approaches in this setting 22–25, while their utility for reconstructing cellular 70 dynamics and trajectories from snapshot time-series data has not yet been systematically evaluated. 71 It therefore remains unclear whether embeddings produced by such foundation models actually 72 confer an advantage over an HVG -based baseline for reconstructing cellular dynamics from 73 snapshot time-series data. 74 Cell dynamic processes can be modeled computationally using trajectory inference methods, which 75 order cells along a trajectory based on similarities in their expression patterns 6. Many widely used 76 approaches like graph-based methods (e.g. Monocle 26) reconstruct branched manifolds from 77 neighborhood graphs, and RNA velocity based methods (e.g. scVelo27–30) infer local transcriptional 78 dynamics for individual cells based on unspliced and spliced mRNA. However, these approaches do 79 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint not explicitly model the coupling between successive sampling times. In contrast, optimal 80 transport(OT) based approaches (e.g. Waddington-OT31) are explicitly formulated for time -series 81 designs, in which distinct cell populations are sampled at successive time points 8,32–34. Given 82 embeddings of cells at two or more times, optimal transport estimates couplings that describe how 83 the population at an earlier time is probabilistically mapped into the population at a later time, thereby 84 providing a principled framework for reconstructing population-level flows from destructive snapshot 85 measurements8,31,35–38. In single -cell biology, OT has been used to align distributions across 86 developmental stages or treatment time courses and to recover trajectories and fate biases in 87 systems such as reprogramming, differentiation, cancer development and response to perturbation38. 88 Extensions based on unbalanced OT 39–42 allow total mass to vary between time points and have 89 been introduced to better accommodate proliferation, death, and population expansion or contraction 90 dynamically37. However, single -cell transcriptomes are high -dimensional and sparse, and directly 91 fitting OT in raw expression space is statistically unstable and computationally demanding 8. 92 Consequently, OT-based analyses are typically performed on low-dimensional embeddings obtained 93 by selecting HVG and applying dimensionality reduction methods, and the inferred dynamics depend 94 critically on this embedding choice 43. This dependence motivates a systematic benchmark of 95 whether embeddings learned by single -cell foundation models provide advantages over simple 96 HVG-based baselines for OT-based reconstruction of cellular dynamics from time-series data. 97 Although promising zero-shot transfer has been reported for general single -cell analysis tasks, the 98 suitability of scFM embeddings for trajectory inference and dynamical reconstruction remains 99 unresolved. To address this gap, a systematic benchmark was conducted to test whether zero-shot 100 embeddings from five published scFMs improve trajectory inference relative to a baseline 101 constructed from highly variable genes projected with principal component analysis (HVG -PCA) 102 across a diverse range of complex biological systems, including hematopoietic lineage branching, 103 embryoid body development, the epithelial-to-mesenchymal transition (EMT) and the directed in vitro 104 differentiation of stem cells into pancreatic βcells. We evaluated three canonical scenarios that 105 address fundamental inquiries in developmental and disease biology: backtracking seeks precursor 106 or progenitor states that precede observed samples 44; interpolation aims to recover transient 107 intermediates that connect adjacent sampled time points 45; and extrapolation anticipates 108 downstream states beyond the last observation under continued progression of the process 46. And 109 performance was quantified with complementary metrics: (i) distribution recovery; (ii) pseudotime 110 correlation; (iii) local velocity coherence47–49. This benchmarking design isolates embedding choice 111 from downstream modeling, enabling a direct assessment of whether foundation model embeddings 112 preserve the intricate non -linearities of cellular dynamics or introduce structural biases that could 113 misinform the design of engineered biological systems. 114 115 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint 116

117 Overview of benchmark framework. 118 To assess whether zero -shot embeddings from single -cell foundation models (scFMs ) provide 119 advantages for trajectory inference, a benchmark was designed that isolates representation quality 120 from downstream modeling. Given the sparsity and dimensionality of single -cell expression data, 121 performing trajectory inference directly in the original gene expression space is both computationally 122 challenging and sensitive to noise. Consequently, in practice, cellular dynamics are typically inferred 123 after projecting cells into a lower-dimensional representation, which makes the choice of embedding 124 a central determinant of downstream trajectory reconstruction 7. As illustrated in Fig. 1a , time -125 stamped snapshot expression profiles were first mapped by each scFM, as well as an HVG -PCA 126 baseline, into low dimensional cell embeddings. Trajectory inference was then performed directly in 127 this embedding space to reconstruct cellular dynamics, thereby separating the effect of 128 representation learning from that of downstream dynamical inference. 129 Our benchmark spans multiple published time-series snapshot datasets (Methods), covering a range 130 of experimental systems and scales, from approximately 3,000 to 49,000 cells (Fig. 1b, a complete 131 list of datasets is provided in Supplementary Table 1), encompassing fundamental processes such 132 as differentiation, development, pathological transitions and cellular reprogramming. Each dataset 133 is processed by six embedding pipelines: five foundation model embedders (Geneformer, 134 Genecompass, scGPT, UCE, and scFoundation ) and a n HVG-based baseline. On top of these 135 representations, we apply four trajectory inference methods that are closely related through optimal 136 transport formulations and their entropy-regularized variants, including Dynamical Optimal Transport 137 (DOT), Unbalanced Dynamical Optimal Transport (UOT), Dynamical Schrödinger Bridge 50,51, and 138 Regularized Unbalanced Optimal Transport (RUOT) 37,52–54 , yielding a comprehensive set of 139 “embedding × inference” combinations for systematic comparison. 140 To evaluate temporal generalization, we partition each dataset by sampling time points and construct 141 three complementary tasks (Fig. 1b, right). In backtracking, models are fit on later time points and 142 used to reconstruct initial states. In interpolation, intermediate time points are held out and predicted 143 from other time points. In extrapolation, models are fit on early time points and evaluated on later, 144 unseen time points (Methods). This design mirrors practical scenarios with incomplete temporal 145 coverage and tests whether embeddings and inference methods can recover cellular dynamics 146 beyond the observed time window. 147 Because different embedders and trajectory inference methods generate representations in distinct 148 coordinate systems, we align both the cell embeddings and inferred trajectories into a shared latent 149 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint space prior to evaluation (Fig. 1c)55. This alignment normalizes for coordinate differences, ensuring 150 that performance comparisons reflect recovered dynamical structure rather than representational 151 arbitrariness. Within this shared space, performance is quantified using three complementary 152 metrics that capture the distributional, directional, and ordinal aspects of cellular dynamics (Fig. 1d, 153 Methods). Specifically, we employ the Wasserstein-1 distance (Earth Mover’s Distance, EMD) to 154 quantify discrepancies between predicted and observed cell state distributions at held-out time points 155 (distributional); a velocity coherence score to assess the local consistency of inferred velocity vectors 156 (directional); and a pseudotime correlation, defined as the Spearman correlation between inferred 157 pseudotime and the reference chronological order (ordinal). Together, these metrics provide an 158 integrated assessment of each pipeline’s ability to reconstruct cell state distributions, flow directions, 159 and temporal progression. 160 161 162 163 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint 164 Fig. 1 | Benchmark workflow and evaluation pipeline for reconstruction of cellular dynamics in 165 foundation model embedding space. a, Overview of the workflow. Single-cell transcriptomic snapshot data 166 are embedded by large pretrained models to obtain cell representations. Trajectory inference is subsequently 167 performed in the embedding space to reconstruct cellular dynamics. b, Detailed benchmark design. Snapshot 168 data are embedded by five foundation models and one baseline method (i.e. PCA embedding with highly 169 variable genes) to generate time -resolved cell embeddings. Dynamic optimal transport ( DOT), Unbalanced 170 Dynamical Optimal Transport (UOT), Dynamical Schrödinger Bridge, and Regularized Unbalanced Optimal 171 Transport (RUOT), four dynamic inference methods are then applied. Data are partitioned into training and 172 test sets based on sampling time points to simulate three dynamical tasks: backtracking, interpolation, and 173 extrapolation. c, For fair comparison, embeddings and inferred trajectories from different models are aligned 174 into a unified latent space, in which all performance metrics are computed. d, Illustration of the three 175 evaluation metrics used in this study. Wasserstein -1 distance (W1) quantifies the distributional divergence 176 between predicted and observed cell state distributions at held-out time points. Velocity coherence measures 177 local directional consistency of inferred dynamics by computing the agreement of velocity vectors among 178 neighbouring cells: for a representative cell (highlighted in orange), velocity vectors of its K nearest neighbours 179 (circled) are compared using cosine similarity. Pseudotime correlation is defined as the Spearman rank 180 correlation between the inferred pseudotime ordering and the reference pseudotime. 181 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint HVG outperform zero -shot foundation model embeddings across the 182 majority of tasks and metrics. 183 To evaluate whether zero-shot embeddings from single-cell foundation models provide advantages 184 for optimal transport (OT) based trajectory reconstruction beyond a simple HVG–PCA baseline, we 185 applied the benchmark framework described above to six embedders in combination with four 186 trajectory inference methods across multiple time -series single -cell datasets. In the aligned 187 evaluation space, HVG –PCA embeddings yield the strongest distributional recovery across tasks 188 and datasets. As summarized by the Wasserstein-1 (EMD) distance (lower is better; Fig. 2a), HVG–189 PCA attains the lowest discrepancies between predicted and observed held -out cell -state 190 distributions in backtracking, interpolation, and extrapolation settings. For each dataset –task 191 combination, points represent the mean performance across the four inference methods, with 192 whiskers indicating the corresponding range. Among the foundation models, Geneformer and scGPT 193 are the most competitive, but both remain inferior to the HVG baseline, whereas scFoundation , 194 based on direct expression value encoding (Methods), shows the weakest performance. This 195 suggests that traditional HVG representations better preserve the multi -modal heterogeneity 196 inherent in complex transitions, whereas foundation model embeddings may suffer from 197 representation collapse, particularly in the challenging regimes of backtracking and extrapolation 198 where the model must infer progenitor states or future disease states beyond the observed window. 199 Temporal ordering accuracy, quantified by the Spearman correlation with reported pseudotime, 200 exhibits a similar overall pattern. Because pseudotime lacks an absolute ground truth, correlations 201 are computed against the reference pseudotime provided in the original data publications after 202 alignment in the shared latent space (Methods). As shown in Fig. 2b , the HVG –PCA baseline 203 achieves the highest pseudotime correlation on the EMT dataset, reaching a Spearman’s ρ of 0.892 204 under the interpolation setting. In datasets with more complex temporal and branching structures 205 (EBdata and HSPC), inspection of the DPT56 pseudotime annotations reveals weaker concordance 206 with experimental sampling time ( Supplementary Fig. 1). Specifically, pseudotime values are not 207 uniformly distributed across sampling time points and do not exhibit a clear monotonic increase with 208 progression. In these settings, pseudotime correlation is therefore treated as a secondary metric and 209 interpreted alongside distributional and velocity-based measures. 210 Directional consistency further supports this overall trend. Velocity coherence quantifies the local 211 smoothness of inferred dynamics by measuring agreement among velocity vectors within cell 212 neighborhoods (Methods). As shown in Fig. 2c, velocity coherence (higher is better) is reported 213 across tasks for each dataset and trajectory inference method; points correspond to individual 214 inference methods, while the red line indicates the median value across datasets and embeddings. 215 Overall, HVG–PCA achieves the highest velocity coherence across the majority of datasets and task 216 regimes, indicating that OT -based trajectory inference tends to recover smoother and more self -217 consistent dynamical flows in HVG–PCA space. Notably, in the more challenging backtracking and 218 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint extrapolation settings, certain foundation model embeddings , particularly GeneCompass and 219 scFoundation exhibit comparable performance to HVG –PCA and, in some datasets, achieve the 220 highest velocity coherence. These results suggest that while HVG –PCA provides the most robust 221 substrate for directional consistency overall, specific foundation models can be competitive under 222 particular dynamical regimes. 223 To provide qualitative intuition for these quantitative differences, we visualize inferred trajectories on 224 the EMT dataset in Fig. 2d . In a two -dimensional PCA projection of the embedding space, 225 predictions at held-out time points substantially overlap the observed cell populations for HVG–PCA 226 across all task settings, effectively reconstructing unobserved cell state in EMT process. In contrast, 227 trajectories inferred in the scFoundation embedding space show limited overlap and fail to 228 reconstruct the withheld distributions to a comparable extent. Taken together, these results indicate 229 that, under zero -shot settings, HVG –PCA embeddings provide a more reliable substrate for 230 trajectory inference than current single -cell foundation model embeddings, with the largest 231 performance gaps observed in distributional accuracy and local directional coherence, and with 232 backtracking and extrapolation representing the most challenging regimes. 233 234 235 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint 236 Fig. 2 | Zero-shot embeddings underperform the highly variable gene (HVG) baseline across three 237 dynamical reconstruction tasks. a, Distribution of Wasserstein-1 (W₁) distances between predicted and 238 ground-truth cell distributions, evaluated on all datasets for back tracking, interpolation, and extrapolation. 239 Lower W₁ values indicate more accurate recovery of the true distributions. b, Spearman correlation between 240 inferred pseudotime and reference pseudotime provided by the original dataset publications, computed 241 across all tasks and datasets. Higher correlation coefficients reflect stronger agreement with temporal 242 trajectories. c, Local velocity coherence across datasets and tasks. Scores range from 0 to 1, with higher 243 values denoting greater consistency among velocity vectors within local cell neighborhoods. d, Trajectory 244 reconstructions on the EMT dataset for the best -performing embedding (HVG) and the worst -performing 245 embedding (scFoundation). Six subplots are organized into two rows —top row: HVG; bottom row: 246 scFoundation—and three columns corresponding, from left to right, to back tracking, interpolation, and 247 extrapolation. In each subplot, true cell s are shown as circles, predicted cells as “×” markers, and inferred 248 velocity vectors are overlaid as streamlines. 249 250 251 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint Sensitivity analyses confirm robustness to alignment strategy, reference 252 space, and latent dimensionality. 253 Because our benchmark integrates heterogeneous embeddings and reconstructs cellular dynamics 254 in a reduced latent space, it is important to assess whether the observed performance differences 255 depend on specific alignment choices or modeling hyperparameters. We therefore conducted a 256 series of sensitivity analyses to evaluate the robustness of our conclusions with respect to 257 embedding alignment, reference space selection, and latent dimensionality. 258 We first examined the role of embedding alignment, which is required to place representations 259 produced by different models which are often differing in scale, orientation, and dimensionality into 260 a common coordinate system prior to evaluation. Using the EMT extrapolation task as an illustrative 261 example (Fig. 3a), alignment of HVG, Geneformer, and scFoundation embeddings to a consensus 262 latent space substantially reduces inter -embedding differences in scale and orientation, enabling 263 direct comparison of inferred trajectories and predicted cell state distributions. Aligned embeddings 264 for all datasets are visualized in Supplementary Fig. 2–6. In contrast, analyses performed in the 265 original, unaligned embedding spaces are confounded by pronounced scale mismatches and 266 rotations, which impair comparability and can bias downstream metrics. 267 Having established the necessity of alignment, we next examined whether the benchmark 268

depend on the specific alignment strategy or choice of reference space. Comparing 269 aligned and unaligned evaluations reveals that the absence of alignment introduces modest shifts in 270 Wasserstein-1 distances and can alter the relative ranking of embedders in some cases ( Fig. 3b), 271 underscoring the importance of embedding alignment for fair comparison. Nevertheless, when 272 alignment is applied, varying the alignment procedure yields stable relative rankings of embedders 273 by distributional accuracy, as measured by the Wasserstein-1 distance. To further assess potential 274 reference-space bias, we repeated the interpolation task while alternately treating each embedding 275 as the reference space for alignment (Fig. 3c). Although absolute W1 values shift modestly, which 276 often favoring the embedding chosen as the reference , the overall performance ordering across 277 embedders remains largely unchanged. Consistent trends are observed for pseudotime correlation 278 and velocity coherence under different alignment and reference space choices (Supplementary Fig. 279 7), indicating that the robustness of our conclusions is not specific to a single evaluation metric. 280 Based on these sensitivity analyses, all main results are reported using Generalized Procrustes 281 Analysis (GPA) alignment to a consensus latent space (Methods)57. 282 We additionally assessed sensitivity to a key optimal transport hyperparameter: the dimensionality 283 of the latent space used to learn cellular dynamics. Across latent dimensions of 5, 10, 20, and 50, 284 pseudotime correlations preserve the overall ranking of embedders ( Fig. 3d), although the optimal 285 dimensionality varies across embeddings. For example, in the Genecompass embedding space, a 286 higher-dimensional latent representation substantially outperforms lower -dimensional alternatives, 287 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint indicating that appropriate dimensionality selection can materially improve performance within a 288 given embedding. 289 Velocity coherence exhibits a similar pattern of robustness at the cross -embedding level (Fig. 3e), 290 with relative performance largely preserved across latent dimensionalities. Within individual 291 embeddings, lower dimensional latent spaces generally yield higher local velocity coherence, 292 consistent with the recovery of smoother and more locally self -consistent flow fields on a more 293 compact manifold. Sensitivity analyses across all datasets are provided in Supplementary Fig. 8. 294 Taken together, these sensitivity analyses demonstrate that our central findings i.e. the superior 295 performance of HVG-based embeddings over zero-shot foundation model embeddings across most 296 tasks and evaluation metrics, are robust to the choice of alignment strategy, reference space, and 297 latent dimensionality. This robustness supports the use of the benchmark as a general and reliable 298 framework for evaluating embedding representations in dynamic single cell settings. 299 300 301 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint 302 Fig. 3 | Sensitivity analyses with respect to alignment, reference space, and embedding dimensionality. 303 a, Visualization of the EMT dataset illustrating evaluations with and without alignment. For HVG, Geneformer, 304 scFoundation, and an integrated “merge” representation (combining model -specific embeddings), observed 305 cell states and predicted states are displayed as scatter plots. The aligned condition is shown in the top row 306 and the unaligned condition in the bottom row; columns correspond to the four representations. Observed and 307 predicted states are distinguished by marker shape . b, Wasserstein-1 distance evaluated across the three 308 dynamical tasks under alternative alignment settings. “Unalign” denotes evaluation in each model’s native 309 embedding space without alignment. “ Aligned” denotes orthogonal Procrustes alignment to a reference 310 embedding, corresponding to an orthonormal rotation and reflection that minimizes squared distances between 311 embeddings. c, Effect of the choice of reference space on Wasserstein-1 distance. Results are reported after 312 mapping embeddings into different reference spaces; “consensus” denotes a shared latent space obtained by 313 aggregating embeddings aligned across all models. d, Heatmaps showing the effect of embedding 314 dimensionality on pseudotime correlation across the three dynamical tasks. The horizontal axis indicates the 315 number of retained principal components. Color encodes the Spearman correlation between inferred 316 pseudotime and the reference pseudotime, with blue indicating lower and red indicating higher correlation 317 (values approaching 1 indicate stronger agreement). e, Heatmaps showing the effect of embedding 318 dimensionality on local velocity coherence across the three tasks. The horizontal axis indicates the number of 319 retained principal components. Color encodes coherence, with blue indicating lower and red indicating higher 320 values (values approaching 1 indicate stronger local directional consistency of velocity vectors). 321 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint Zero-shot foundation model embeddings attenuate temporal and 322 branching structure relevant for trajectory inference . 323 To investigate why zero -shot foundation model embeddings exhibit systematic limitations on 324 dynamics tasks, we examined how they represent temporal structure. We defined a time variance 325 ratio (TVR), the fraction of total variance attributable to time (Methods), to quantify the extent to 326 which time-scale differences are preserved. As illustrated by a linear discriminant analysis on the 327 EMT dataset ( Fig. 4a)58, foundation model embeddings exhibit markedly reduced TVR relative to 328 the HVG baseline, indicating diminished temporal separability. Pairwise energy distances between 329 time points, visualized as a heatmap, further support this observation: time points exhibit increased 330 proximity and reduced discriminative resolution in foundation model spaces than in the HVG 331 representation. Among the foundation models, scFoundation retains comparatively greater temporal 332 separation, although still below HVG and without translating into consistent improvements across 333 downstream dynamics metrics. 334 We next examined the relationship between temporal variance preservation and dynamics 335 performance. Across four datasets and three dynamical scenarios (Fig. 4b), higher TVR is generally 336 associated with lower Wasserstein -1 distances, indicating improved recovery of held -out cell-state 337 distributions when time-scale differences are preserved. However, this negative association is not 338 universal. For example, in the EMT dataset ( Supplementary Fig. 9), Geneformer, scGPT, and 339 GeneCompass display very low TVR, with time points nearly co -localized in the embedding space, 340 yet achieve low W1 distances. In this regime, predicting held -out distributions becomes 341 comparatively trivial because observed and target states are already close in the embedding. Similar 342 analyses relating TVR to pseudotime correlation and local velocity coherence in multiple datasets 343 (Supplementary Fig. 10, 11 ) further support that compression of time -structured variation is 344 generally accompanied by degraded dynamical reconstruction. 345 Motivated by this pattern, we asked whether foundation model embeddings might be attenuating 346 temporal signal in a manner analogous to batch -effect correction. In the EMT dataset, explicitly 347 applying Harmony 59 batch correction similarly reduces between time variance and equalizes W1 348 across embeddings, but at the cost of lower pseudotime correlation and velocity coherence 349 (Supplementary Fig. 12). This controlled perturbation supports the view that removing time -350 structured variation, whether by design or implicitly in foundation model embeddings, tends to impair 351 dynamical reconstruction. 352 We then asked whether a similar “compression” occurs along branching, fate related axes. In the 353 Veres human pancreatic differentiation dataset, HVG embeddings preserve two clear branches from 354 NKX6-1⁺ progenitors toward SC -β and SC -EC fates ( Fig. 4c , left), and extrapolated RUOT 355 trajectories follow these reference branches in a fate-specific manner (Fig. 4c, middle). By contrast, 356 in the Geneformer embedding, SC -β and SC -EC populations are less clearly separated, and 357 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint extrapolated flows are less tightly aligned with the two terminal branches (Fig. 4c, right). Consistently, 358 cosine similarity analyses show increased similarity between SC -β and SC -EC states in multiple 359 foundation model embedding spaces ( Fig. 4d ). An embedding incorporating protein sequence 360 information during pretraining (UCE) exhibits comparatively improved separation between these two 361 endocrine populations, but does not yield consistently better dynamical performance. A similar 362 pattern is observed in mouse hematopoiesis ( Fig. 4e, f ), in this bifurcating system, large -model 363 embeddings increase the similarity between Neutrophil and Monocyte populations relative to HVG, 364 and RUOT-based velocity streamlines become less branch-aligned in the large-model space than in 365 the HVG space, indicating that branch -specific biological variation is systematically compressed in 366 zero-shot foundation model embeddings. 367 Overall, temporal and branching analyses together indicate that zero -shot foundation model 368 embeddings systematically compress both time -scale differences and branch -specific biological 369 variation. Such compression of biologically meaningful structure, likely arising from an 370 over-correction of variation treated as nuisance (e.g. batch -like effects), offers a plausible 371 explanation for the consistently weaker performance of downstream dynamical methods in 372 large-model embedding spaces relative to HVG-based embeddings. 373 374 375 376 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint 377 Fig. 4 | Batch correction like temporal compression and impaired branch resolution in large -model 378 embeddings and their association with diminished performance in dynamical reconstruction. a, EMT 379 dataset: for each embedding model, an LDA visualization colored by sampling time and a heatmap of pairwise 380 energy distances between time points are shown; panel titles report the Temporal Variance Ratio (TVR, higher 381 indicates greater temporal separation). Reduced TVR together with attenuated between–timepoint distances 382 indicates temporal compression. b, Association between temporal separation and reconstruction error across 383 four datasets and three tasks. Each scatter plot relates TVR to the Wasserstein -1 distance (W1); points are 384 colored by embedding model. A negative association is observed, indicating that lower temporal separation is 385 linked to higher distributional error. c, Human pancreatic differentiation (Veres dataset). Left: reference 386 trajectory depicting differentiation from NKX6 -1⁺ progenitors toward SC -β and SC -EC fates. Middle: HVG 387 embedding visualized with UMAP and extrapolation velocity streamlines inferred by RUOT (regularized 388 unbalanced optimal transport). Right: Geneformer embedding visualized with UMAP and RUOT extrapolation 389 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint streamlines. In the HVG embedding, SC-β and SC-EC branches remain well separated and extrapolated flows 390 follow the reference branches; in the Geneformer embedding, these fates are partly merged and extrapolated 391 flows become less branch -aligned. d, Pairwise similarity between SC -β cells and SC -EC cells across 392 embeddings from different single -cell foundation models (scFMs), quantified by cosine similarity. e, Mouse 393 HSPC. From left to right: reference trajectory; HVG embedding with UMAP visualization and RUOT 394 extrapolation velocity streamlines; Geneformer embedding with UMAP visualization and RUOT extrapolation 395 velocity streamlines. f, Pairwise similarity between Neutrophil and Monocyte populations across scFM 396 embeddings. 397 398 399 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint Global benchmarking reveals consistent advantages of HVG -based 400 embeddings across methods and datasets 401 To synthesize the benchmark results across all evaluated settings, we next summarized dynamical 402 performance for all combinations of embedding spaces, trajectory inference methods, tasks, and 403 datasets using an integrated heatmap representation (Fig. 5). This global view enables simultaneous 404 assessment of consistency and variability across the full benchmarking landscape. 405 Across the majority of settings, the HVG -based baseline exhibits superior performance relative to 406 zero-shot single-cell foundation model embeddings, most prominently in distributional accuracy as 407 measured by the EMD. This advantage is observed consistently across datasets, task regimes, and 408 optimal-transport–based inference methods, indicating that the superiority of HVG -based 409 representations is not driven by a particular modeling choice or experimental system. Among 410 foundation models, Geneformer generally achieves stronger distributional recovery, whereas 411 scFoundation tends to yield higher local velocity coherence in some settings; however, no single 412 foundation model embedding consistently dominates across all metrics or tasks. 413 In Fig. 5, we summarize global benchmark performance using rank-based heatmaps rather than raw 414 metric values. For each dataset, task, and trajectory inference method, embeddings are ranked 415 according to distributional accuracy (Wasserstein -1 distance) and local directional consistency 416 (velocity coherence), and these ranks are then aggregated to provide a comparative overview across 417 settings. This rank-based representation emphasizes relative performance trends and reduces the 418 influence of dataset specific scale differences across metrics. Heatmaps of the corresponding raw 419 metric values for all settings are provided in the Supplementary Note3, Fig13, Fig14. 420 Taken together, this global analysis confirms that the performance advantages of HVG -based 421 embeddings over zero -shot foundation model embeddings are robust across inference methods, 422 datasets, and task settings. At the same time, the heterogeneity revealed by the heatmap highlights 423 that current foundation models exhibit metric - and task-specific strengths, suggesting that future 424 improvements will require embedding strategies that better preserve temporal and biological-specific 425 structure rather than relying on zero-shot representations alone. 426 427 428 429 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint 430 Fig. 5 | Global rank-based summary of embedding performance across datasets and tasks. Heatmap 431 showing the relative ranking of embedding methods across all datasets, three dynamical tasks (backcasting, 432 interpolation, and extrapolation), and optimal-transport–based trajectory inference settings. For each dataset–433 task–inference combination, embeddings are ranked based on distributional accuracy (EMD) and local velocity 434 coherence, with higher ranks indicating better relative performance. For interpolation tasks involving multiple 435 intermediate test time points, ranks are aggregated across these settings. The top two performing embeddings 436 for each dataset –task combination are highlighted. Heatmaps of the corresponding raw metric values are 437 provided in the Supplementary Information (Supplementary Note 3). 438 439 440 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint

441 In this study, a systematic benchmark was conducted to compare single -cell foundation models 442 (scFMs) with a highly variable gene (HVG) baseline for cellular dynamics reconstruction across three 443 tasks: backtracking, interpolation, and extrapolation. All methods were evaluated in a shared aligned 444 embedding space using complementary metrics that capture different aspects of dynamical 445 reconstruction: (i) distributional recovery, quantified by the Wasserstein -1 distance (Earth Mover’s 446 Distance), (ii) global agreement with reference pseudotime, quantified by Spearman correlation, and 447 (iii) local velocity coherence, which measures neighborhood -level consistency of inferred velocity 448 vectors. 449 Based on their training paradigm and architecture, scFMs pretrained with transformer models on 450 large single-cell corpora would be expected to capture gene –gene dependencies and higher order 451 structure through self -attention, thereby encoding information relevant to cellular dynamics and 452 potentially outperforming a simple HVG baseline. However, contrary to this expectation, we found 453 that the HVG-based embedding consistently outperformed scFM embeddings in most settings and 454 across metrics. This pattern was robust to the choice of alignment procedure and latent 455 dimensionality. Among scFMs, Geneformer and scGPT achieved comparatively stronger 456 performance on the Wasserstein -1 metric, whereas scFoundation yielded higher local velocity 457 coherence. Pseudotime correlation appeared to be more sensitive to data topology: in datasets with 458 bifurcating trajectories, several scFM embeddings tended to linearize branches, thereby obscuring 459 biologically meaningful branching structure while artificially inflating pseudotime correlations. Beyond 460 intrinsic data geometry, our analyses also indicated that zero -shot scFM embeddings can over 461 correct batch effects, reducing temporal separability and degrading performance on dynamical tasks. 462 Taken together, these results show that, in their current zero -shot form, scFM embeddings are not 463 superior to a straightforward HVG-based baseline for reconstructing cellular dynamics, and may in 464 fact suppress temporal variation by treating it as batch-like noise. 465 A possible explanation for these observations lies in the inductive biases introduced by current self-466 supervised training objectives and model architectures. Most scFMs are trained with reconstruction 467 or masked token prediction losses that encourage invariance to technical perturbations and 468 emphasize stable co-expression structure. Such objectives are well aligned with tasks that require 469 robust cell identity representations, but they inherently prioritize features that are frequent and 470 persistent across cells60. By contrast, transient transcriptional programs that govern short lived state 471 transitions are less consistently sampled and therefore contribute less to the self-supervised signal. 472 As a result, the learned embeddings are dominated by static, identity related variation (general 473 signals), whereas temporally restricted dynamics (specific signals) are under -represented. This 474 interpretation is consistent with recent benchmarking studies reporting that zero -shot embeddings 475 from scFMs often underperform HVG -based baselines on tasks that depend on fine -grained or 476 context specific variation, including perturbation responses prediction25,60. Together, these 477 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint observations suggest that current scFMs preferentially capture general, context -agnostic structure 478 at the expense of process-specific temporal signals. 479 Architectural and encoding choices likely reinforce this bias. Rank - or bin-based preprocessing of 480 gene expression compresses dynamic range and down-weights gradual expression changes, while 481 pooling gene token embeddings into a single cell level representation further compresses local 482 structure on the expression manifold 18. These design choices stabilize training and enhance 483 invariance to nuisance variation, but they also smooth out fine -scale curvature and reduce 484 separability between time points and branches. In developmental and perturbation trajectories, 485 where relevant signals are often concentrated in a subset of genes and confined to specific temporal 486 windows, such smoothing can collapse biologically distinct intermediate states into a more 487 homogeneous representation22,24,61. This offers a geometric explanation for the reduced temporal 488 variance and branch blurring observed in zero-shot scFM embeddings. 489 Taken together, our findings imply that general purpose scFMs , as currently trained, may be 490 intrinsically better suited to static tasks such as cell type classification and batch integration than to 491 reconstructing dynamical processes. This functional divergence is rooted in the model’s architecture 492 and pretraining objectives, which we illustrate conceptually in Fig. 6. Current scFMs map diverse 493 biological contexts into a relatively uniform embedding space with a collapsed geometry. By 494 prioritizing identity -related stability, these models treat most temporal variation as noise to be 495 removed, thereby obscuring short -lived transitional programs and directly degrading downstream 496 performance on dynamics. Closing this gap will likely require making dynamical structure an explicit 497 target of representation learning rather than a purely downstream objective. In particular, models will 498 need to preserve temporal differences between states and maintain the separation of distinct 499 developmental or perturbation trajectories. Overall, d eveloping scFMs that explicitly balance 500 robustness to technical variation with retention of biologically meaningful temporal structure 501 represents a key avenue for future work. 502 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint 503 Fig. 6 | Conceptual understanding of temporal information loss in scFM embeddings: Inductive biases 504 in current single-cell foundation models (scFMs) prioritize stable cell identity over transient dynamical 505 signals. Input (Left): Time-resolved biological data with structured trajectories and distinct branching lineages 506 (t0 – t3). Mechanism (Center): Self-supervised pretraining introduces an inductive bias (represented by the 507 balance) that favors persistent “General Signals” (e.g., cell type) while down weighting transient, process 508 specific “Specific Signals”. Embedding (Right): The resulting latent space exhibits collapsed geometry, where 509 temporal variation is suppressed and distinct cellular states overlap. Impact (Bottom): Loss of fine -grained 510 temporal information degrades downstream dynamical tasks, manifested as Noisy Velocity vectors, 511 Ambiguous Trajectories, and Weak Transport between time points. 512 513 514 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint

515 Benchmark framework overview 516 The benchmark framework comprises three core steps: (i) scRNA -seq data are reduced to a 517 low-dimensional cell space using different embedding methods; (ii) trajectory -inference algorithms 518 are applied in this space to model cellular dynamics; and (iii) predicted cell states and trajectories 519 are aligned to a consensus reference space, where three evaluation metrics are computed. 520 521 Foundation models and baseline 522 For all single -cell foundation models (scFMs), we first obtain cell embeddings from the model 523 encoders. In this study, we use the [CLS] token embedding from the encoder as the representation 524 of each cell, construct a time-series of cell embeddings with the original hidden dimension, and then 525 apply PCA to reduce them to a low-dimensional space suitable for trajectory inference. 526 Geneformer. Geneformer V2 is a foundation al transformer model pretrained on ~104 million non -527 cancer human single-cell transcriptomes. It encodes each cell using a rank-based representation in 528 which genes are ordered by relative expression and scaled using population -level statistics , 529 emphasizing cell-state-informative genes while reducing the influence of ubiquitous housekeeping 530 genes. For each cell, the top 4,096 ranked genes are selected as the input sequence. The encoded 531 sequences pass through multiple Transformer encoder layers and are trained using a masked-gene 532 prediction objective, enabling fully self-supervised learning on unlabeled data. In this study, we use 533 the Geneformer-12L-4096i-104M model, which is pretrained on 104 million human single cells and 534 takes 4,096 ranked genes as input, encoded by a 12-layer Transformer encoder. 535 Genecompass. Genecompass is a knowledge -informed cross -species foundation model self -536 supervised pre-trained on an extensive dataset scCompass-126M, which comprises over 120 million 537 high-quality human and mouse single -cell transcriptomes. For each cell, the model takes the top 538 2048 ranked genes after normalization and ranking of gene expression values, and incorporates four 539 types of prior biological knowledge (including GRN, promoter sequence, gene family and co -540 expressions) encoded within a unified embedding space. It employs a 12 -layer Transformer 541 architecture and adopts a masked language modeling objective which randomly masks 15% of gene 542 inputs in each cell. 543 UCE. UCE is a cross-species single-cell foundation model, which employs a 33 -layer Transformer 544 encoder with 650 million parameters. It is pre-trained on over 36 million cells spanning eight species 545 in a completely self-supervised way without any data annotations. Gene embeddings are initialized 546 using the pretrained protein language model ESM -262, incorporating protein biological knowledge, 547 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint genes are sampled based on expression probability and sorted by chromosomal order; the 548 pretraining task uses masked binary prediction, training the model to determine whether a gene is 549 expressed in a cell. This enables UCE to establish a unified biological latent space capable of 550 representing any cell, regardless of tissue or species. In this study, we use the 33-layer UCE model 551 to obtain cell embeddings. 552 scFoundation. scFoundation is a large -scale pretrained model with 100M parameters, 553 which models 19,264 genes and is pre -trained on over 50 million scRNA -seq data. It adopts an 554 asymmetric Transformer-like encoder-decoder architecture xTrimoGene, whose embedding module 555 converts continuous gene expression values into high -dimensional learnable vectors rather than 556 using discretized values, and employs an asymmetric encoder-decoder structure to handle the high 557 sparsity of single-cell data, enabling efficient learning of all gene relationships. 558 scGPT. scGPT is a single-cell foundation model for multiple single-cell omics data analysis applying 559 large-scale pretrained Transformer. The model takes 1,200 HVGs as input and introduces a cell -560 dependent value binning technique to discretely encode gene expression levels. It employs a special 561 attention masking mechanism to adapt the Transformer to non -sequential omics data while 562 incorporating condition tokens to encode metadata, enabling it to simultaneously learn context -563 aware embeddings for both cells and genes. The pre-trained scGPT can be optimized to achieve 564 superior performance across diverse downstream applications through fine-tuning. Here we apply 565 the version whole-human model pretrained on 33 million normal human cells. 566 HVG baseline. As a traditional baseline, raw transcript counts are library-size normalized to 1 × 105 567 and log1p-transformed using Scanpy v1.11.3 [ref. 63]. The top 2,000 highly variable genes are then 568 selected, and PCA is applied to obtain a low-dimensional space for trajectory inference. 569 570 Trajectory inference methods 571 To build a comprehensive framework for benchmarking scFM embeddings from cellular dynamics, 572 we apply four OT -based methods to reconstruct cell -state transitions from time -series scRNA-seq 573 datasets. The core idea of optimal transport relies on finding a mapping that minimizes a specific 574 cost functional when transporting a probability mass from an initial distribution to a terminal state. In 575 the context of temporally resolved snapshots data, this problem can be interpreted as an evolution 576 of cell density 𝜌(𝒙, 𝑡) governed by a continuity equation and energy minimization principle. This 577 section provides a brief review of four typical OT-based trajectory inference methods. 578 579 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint The four methods considered are: Dynamical Optimal Transport (DOT), Unbalanced Dynamical 580 Optimal Transport (UOT), Dynamical Schrödinger Bridge, and Regularized Unbalanced Optimal 581 Transport (RUOT). In practice, all four approaches were implemented using the Python package 582 DeepROUT53, with method -specific behavior controlled primarily through two parameters: sigma, 583 which modulates the strength of entropy or Tikhonov regularization, and use_mass, which specifies 584 whether to enforce mass conservation (balanced OT) or allow for unbalanced transport. Adjusting 585 these parameters allows us to flexibly switch between standard, regularized, and unbalanced OT 586 formulations while maintaining a unified computational framework across all embedding spaces and 587 tasks. 588 Dynamical Optimal Transport (DOT) 589 Dynamical optimal transport (DOT) is a landmark work which models cell development process as 590 a deterministic transport characterized by 𝑑𝑿𝒕 = 𝒗(𝑿𝒕, 𝑡)𝑑𝑡 [ref. 64]. Mathematically, this framework 591 involves finding the velocity field 𝒗(𝒙, 𝑡) and cell density 𝜌(𝒙, 𝑡) that minimize the total transport cost: 592 593 i𝑛𝑓 "#(𝒙,'),𝒗(𝒙,')* 1 1 1 2 ℝ! - . ‖𝒗(𝒙, 𝑡)‖/ /𝜌(𝒙, 𝑡)𝑑𝒙𝑑𝑡 (1) 594 with 𝒗(𝒙, 𝑡) and 𝜌(𝒙, 𝑡) subject to the continuity equation constraint and boundary conditions: 595 596 𝜕𝜌(𝒙, 𝑡) 𝜕𝑡 + ∇08𝒗(𝒙, 𝑡)𝜌(𝒙, 𝑡)9 = 0, 𝜌(∙ ,0) = 𝜌., 𝜌(∙, 𝑇) = 𝜌- (2) 597 598 However, the strict conversation assumption in equation (2) fails to consider mass changing which 599 is popular in biological systems due to cell proliferation and apoptosis. 600 Unbalanced Dynamical Optimal Transport (UOT) 601 To address the limitation of mass conservation, unbalanced dynamical optimal transport (UOT) 602 introduces a growth term 𝑔(𝒙, 𝑡) that represents cell growth or death rate, and rewrite equation (2) 603 as: 604 605 𝜕𝜌(𝒙, 𝑡) 𝜕𝑡 + ∇08𝒗(𝒙, 𝑡)𝜌(𝒙, 𝑡)9 = 𝑔(𝒙, 𝑡)𝜌(𝒙, 𝑡), 𝜌(∙ ,0) = 𝜌., 𝜌(∙, 𝑇) = 𝜌- (3) 606 607 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint The Wasserstein–Fisher–Rao (WFR) distance is often used to quantify the overall shipping cost with 608 respect to both kinetic and growth energy: 609 610 i𝑛𝑓 "#(𝒙,'),𝒗(𝒙,'),1(𝒙,')* 1 1 ?1 2 ‖𝒗(𝒙, 𝑡)‖/ / + 𝛼|𝑔(𝒙, 𝑡)|/ /B 𝜌(𝒙, 𝑡)𝑑𝒙𝑑𝑡 ℝ! - . (4) 611 612 where 𝛼 is a hyperparameter to balance the penalty of growth against transport. While UOT 613 successfully models the unbalanced cellular dynamics, it remains a deterministic transport model 614 that neglects the intrinsic stochasticity in gene expression and cell differentiation. 615 Dynamical Schrödinger Bridge 616 In this case, biological processes is modelled as SDE dynamics characterized by 𝑑𝑿𝒕 = 𝒗(𝑿𝒕, 𝑡)𝑑𝑡 +617 𝝈(𝑿𝒕, 𝑡)𝑑𝑾𝒕, where 𝑾𝒕 is a standard Brownian motion and 𝝈(𝒙𝒕, 𝑡) denotes the random fluctuations 618 in the system. The Schrödinger Bridge problem seeks to identify the most likely evolutionary path 619 between an initial distribution 𝜌. and a terminal distribution 𝜌-, which can be formulated as solving 620 min 2"𝑿3#",2$ 𝑿3#$ 𝒟458𝜇[.,-] 𝑿 J𝜇[.,-] 𝒀 9 with the reference process defined as 𝑑𝒀𝒕 = 𝝈(𝒀𝒕, 𝑡)𝑑𝑾𝒕. This problem 621 can be equivalently transformed into a dynamical form: 622 i𝑛𝑓 "#(𝒙,'),𝒗(𝒙,')* 1 1 L1 2 𝒗𝑻(𝒙, 𝑡)𝒂;𝟏(𝒙, 𝑡)𝒗(𝒙, 𝑡)N 𝜌(𝒙, 𝑡)𝑑𝒙𝑑𝑡 ℝ! - . (5) 623 624 where 𝒂(𝒙, 𝑡) = 𝝈(𝒙, 𝑡)𝝈𝑻(𝒙, 𝑡) and all pairs 8𝜌(𝒙, 𝑡), 𝒗(𝒙, 𝑡)9 satisfy following constraints: 625 626 ∂ρ(𝐱, t) ∂t + ∇=8𝐯(𝐱, t)ρ(𝐱, t)9 = 1 2 ∆𝐱 8𝐚(𝐱, t), ρ(𝐱, t)9, ρ(∙ ,0) = ρ., ρ(∙, T) = ρ? (6) 627 628 This formulation interprets trajectory inference as a stochastic optimal control problem rather than 629 purely deterministic transport, but doesn’t incorporate unbalanced effects. 630 Regularized Unbalanced Optimal Transport (RUOT) 631 Integrating the stochasticity in Schrödinger Bridge problem with cell proliferation effects in UOT, 632 Regularized Unbalanced Optimal Transport (RUOT) assumes particles are governed by the SDE 633 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint model 𝑑𝒙𝒕 = 𝒗(𝒙𝒕, 𝑡)𝑑𝑡 + 𝜎(𝑡)𝑰𝑑𝒘𝒕 . RUOT seeks to minimize a generalized shipping cost 634 accounting for velocity, diffusion, and growth penalties simultaneously: 635 636 i𝑛𝑓 "#(𝒙,'),𝒗(𝒙,'),1(𝒙,')* 1 1 L1 2 ‖𝒗(𝒙, 𝑡)‖/ / + 𝛼Ψ8𝑔(𝒙, 𝑡)9N 𝜌(𝒙, 𝑡)𝑑𝒙𝑑𝑡 ℝ! - . (7) 637 638 where Ψ(∙) represents the growth penalty function, and all pairs 8𝜌(𝒙, 𝑡), 𝒗(𝒙, 𝑡), 𝑔(𝒙, 𝑡)9 are 639 constrained by the unnormalized Fokker-Planck equation: 640 641 𝜕𝜌(𝒙, 𝑡) 𝜕𝑡 + ∇08𝒗(𝒙, 𝑡)𝜌(𝒙, 𝑡)9 = 1 2 𝜎/(𝑡)∆𝒙 𝜌(𝒙, 𝑡) + 𝑔(𝒙, 𝑡)𝜌(𝒙, 𝑡), 𝜌(∙ ,0) = 𝜌., 𝜌(∙, 𝑇) = 𝜌- (8) 642 with vanishing boundary condition: lim|0|→B 𝜌(𝒙, 𝑡) = 0. 643 644 Dynamical Reconstruction Tasks 645 To systematically evaluate the generalization capability of embeddings, we formulated three 646 dynamical reconstruction scenarios: i nterpolation, extrapolation, and backtracking, based on the 647 partition of sampling time points 𝒯 = {𝑡., 𝑡C, … , 𝑡4}. For each task, we hold out a specific time point 648 𝑡'DE' and train the trajectory inference models (DOT, UOT, Schrödinger Bridge, RUOT) on the 649 remaining subset 𝒯𝓉𝓇𝒶𝒾𝓃, subsequently inferring the cell state distribution at held-out time points. 650 Interpolation. This task evaluates the ability to recover transient intermediate states within the 651 observed temporal window. We hold out an intermediate time point 𝑡'DE' = 𝑡K (where 0 < 𝑘 < 𝐾) 652 and train the model using the remaining time points 𝒯𝓉𝓇𝒶𝒾𝓃 = 𝒯 ∖ {𝑡K} . The reconstruction is 653 performed by integrating the learned dynamics starting from the initial distribution 𝜌⋅,'". The cell state 654 𝑋' evolves according to the forward stochastic differential equation (SDE): 655 656 𝑑𝑋' = 𝑣(𝑋', 𝑡)𝑑𝑡 + 𝜎𝑑𝑊', 𝑡 ∈ [𝑡., 𝑡K] (9) 657 658 659 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint subject to the initial condition 𝑋'" ∼ 𝜌⋅,'" . Here 𝑊' denotes a standard Brownian motion . For 660 deterministic transport, we set 𝜎 to 0. For unbalanced settings (UOT, RUOT), the weight of each cell 661 is simultaneously modulated by the growth rate 𝑔(𝑥, 𝑡): 662 𝑑𝑚' 𝑑𝑡 = 𝑔(𝑋', 𝑡)𝑚' (10) 663 664 The predicted distribution 𝜌⋅,'%r is obtained by pushing forward the population through the time interval 665 [𝑡., 𝑡K] and comparing it against the held-out ground truth. 666 Extrapolation. In this scenario, we assess the capacity to predict future developmental states 667 beyond the training horizon. We hold out the last time point 𝑡'DE' = 𝑡4 and train on the preceding 668 sequence 𝒯𝓉𝓇𝒶𝒾𝓃 = {𝑡., … , 𝑡4;C}. Similar to the interpolation task, the inference is conducted by 669 simulating the process starting from 𝑡. up to the target time 𝑡4 . The dynamics follow the same 670 forward SDE and growth equations defined above. 671 Backtracking. Backtracking tests the capability to reconstruct progenitor states. The model is 672 trained on later time points 𝒯𝓉𝓇𝒶𝒾𝓃 = {𝑡C, … , 𝑡4} to predict the unobserved initial state at 𝑡'DE' = 𝑡.. 673 Unlike forward tasks, this requires solving the time -reversal of the diffusion process. The reverse 674 dynamics 𝑋Msss, where 𝜏 = 𝑇 − 𝑡 denotes reversed time, are governed by the backward SDE: 675 676 𝑑𝑋Msss = 𝑣̅(𝑋Msss, 𝜏)𝑑𝜏 + 𝜎𝑑𝑊Mssss （11） 677 678 Here, the effective reverse drift 𝑣̅ is determined by the forward velocity 𝑣 and the score function of 679 the density field: 680 681 𝑣̅(X, 𝜏) = −𝑣(𝑋, 𝑇 − 𝜏) + 𝜎/∇ log 𝜌(𝑋, 𝑇 − 𝜏) （12） 682 683 For deterministic transport (where 𝜎 → 0 ), the score term vanishes, and the reverse velocity 684 simplifies to 𝑣̅ = −𝑣. For unbalanced transport, the growth dynamics are similarly inverted. The mass 685 evolution follows: 686 687 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint 𝑑𝑚M 𝑑𝜏 = −𝑔(𝑋Msss, 𝑇 − 𝜏)𝑚M （13） 688 689 We simulate this coupled reverse process starting from the observed population at 𝑡C to estimate the 690 progenitor distribution 𝜌⋅,'"r . 691 692 Alignment 693 Alignment was used to map heterogeneous cell embeddings into a common space so that 694 downstream metrics could be computed on comparable scales. Each aligner was fitted using cells 695 from training time points and then applied to held-out time point. Both the reconstructed trajectories 696 and the observed cell-expression embeddings were transformed into this common reference space 697 before metric computation. 698 Procrustes alignment. The Procrustes problem seeks a rotation and a single global scale that best 699 align X to Y in the least -squares sense while approximately preserving relative geometry. The 700 optimal transformation is obtained from the singular value decomposition of the cross -covariance 701 between X and Y, with the rotation given by the orthogonal factor and the uniform scale derived from 702 the trace of the aligned matrices 65. This constraint limits distortion and promotes stable, near -703 isometric alignment when dimensionalities are comparable. 704 Generalized Procrustes Analysis . To reduce dependence on an arbitrary reference and avoid 705 favoring any single embedding, Generalized Procrustes Analysis (GPA) was employed to produce 706 a consensus space. GPA iteratively (i) aligns each embedding to the current consensus using 707 Procrustes and (ii) updates the consensus as the average of the aligned configurations; iterations 708 continue until changes fall below a tolerance 66. The consensus was initialized as the mean 709 configuration across embedding spaces. The resulting per -embedding transforms, estimated on 710 training time points, were then fixed and applied to test data. In this manner, all embeddings were 711 compared in a reference-free, consensus-defined space. 712 713 Evaluation Metrics 714 All metrics were computed in a common reference latent space to ensure comparability across 715

and tasks. Under the alignment setting, both observed cells at the held-out time point and 716 model-generated predictions were first projected into this shared space. For each dataset and task, 717 the optimal transport (OT) trajectory model defined a time -conditioned mapping 𝑓(𝑡) from an initial 718 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint snapshot at 𝑡. to later states. Predictions at the held -out time 𝑡∗ were obtained by evaluating 𝑓(𝑡∗) 719 and by interpolating 100 intermediate time points between 𝑡. and 𝑡∗ (fixed a priori). Unless stated 720 otherwise, the neighborhood size was K = 5 for nearest-neighbor operations. 721 Distributional recovery (Wasserstein-1 Distance). Distributional agreement between predicted 722 and observed cells at the held-out time point was quantified using the Wasserstein-1 distance (Earth 723 Mover’s Distance, EMD). Let Ρ = {𝑥O}(𝑗 = 1 … 𝑛) and 𝑄 = {𝑦O}(𝑗 = 1 … 𝑚) denote the aligned 724 embeddings of predicted and observed cells, treated as empirical measures with uniform weights. 725 With Euclidean ground cost 𝐶(𝑥, 𝑦) = ‖𝑥 − 𝑦‖/ the distance was computed as 726 727 𝑊1(𝑃, 𝑄) = 𝑖𝑛𝑓P∈R(S,T)𝐸(0,U)∼P[‖𝑥 − 𝑦‖/] (14) 728 Pseudotime correlation. Global temporal ordering accuracy was quantified using the Spearman 729 rank correlation between predicted and reference pseudotime values. Because pseudotime does 730 not admit a universal ground truth, dataset -specific strategies were used to define the reference 731 ordering. For datasets in which pseudotime annotations were provided by the original studies, these 732 values were used directly as reference pseudotime. For datasets lacking published pseudotime 733 labels (e.g., the Mouse HSPC dataset), reference pseudotime was estimated using diffusion 734 pseudotime (DPT), which infers a one -dimensional progression by integrating diffusion distances 735 along the cell -state manifold56. After OT fitting, 100 evenly spaced intermediate time points were 736 interpolated between 𝑡. and 𝑡∗. For each observed cell, the predicted pseudotime was defined as 737 the mean interpolation time of its K = 5 n nearest predicted neighbors among all interpolated cell 738 states in the aligned latent space. Spearman’s ρ was then computed across observed cells between 739 the predicted pseudotime and the reference pseudotime. 740 Local velocity coherence. Local consistency of directional change was assessed via the coherence 741 of predicted velocities in the aligned evaluation space. After fitting the OT model, a velocity vector 742 𝒗(𝑿𝒕, 𝑡) was obtained for each cell, representing the predicted instantaneous direction of change at 743 state 𝑿𝒕. All generated cells at held-out time points were mapped, together with their velocity vectors, 744 into the aligned space. For each generated cell, the K = 5 nearest generated neighbors were 745 identified based on their positions in this space, and local velocity coherence was defined as the 746 mean cosine similarity between the velocity of the index cell and the velocities of its neighbors, with 747 velocity vectors normalized to unit length prior to similarity computation unless stated otherwise 49. 748 Higher values indicate greater agreement in the predicted directions. 749 750 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint Time variance ratio 751 To quantify how much temporal information is encoded in each embedding space, we computed the 752 time variance ratio (TVR) for each embedding. 753 For a given embedding, let 𝑋 ∈ ℝW×Y be the cell -by-dimension matrix, with cells assigned to time 754 groups 𝑘 = 1. . . 𝐾. We treat time as a categorical grouping variable and decompose the total 755 variance into “between-time” and “within-time” components. We define TVR as: 756 757 𝑇𝑉𝑅 = 𝑡𝑖𝑚𝑒 𝑣𝑎𝑟𝑖𝑎𝑛𝑐𝑒 𝑡𝑜𝑡𝑎𝑙 𝑣𝑎𝑟𝑖𝑎𝑛𝑐𝑒 = ∑ 𝑛K‖𝜇K − 𝜇‖/K ∑ ‖𝑥Z − 𝜇‖/Z (15) 758 759 where xZ is the embedding of cell 𝑖, 𝜇 is the overall mean embedding, 𝜇K is the mean embedding of 760 time group 𝑘, and 𝑛K is the number of cells in group 𝑘. A higher TVR indicates that a larger fraction 761 of the embedding variance is explained by time. 762 763 Data collection and preprocessing 764 We assembled five published time -series single -cell RNA -sequencing (scRNA -seq) datasets for 765 benchmarking ( Supplementary Table 1), encompassing diverse biological systems, temporal 766 resolutions, and dataset sizes (~3,000 –49,000 cells). The datasets include: EMT (GEO: 767 GSE147405)67, Mouse HSPC (Weinreb et al.) 68, Veres (GEO: GSE114412) 69, embryoid bodies 768 (EBdata) (Moon et al.) 70, and HSPC (GEO: GSE226824) 71. Raw count matrices and associated 769 metadata were retrieved to satisfy the input requirements of single-cell foundation models (scFMs). 770 To minimize processing -dependent artifacts in the learned representations, only minimal quality 771 control (QC) was applied: cells with fewer than 200 detected genes and genes expressed in fewer 772 than 3 cells were excluded. No normalization, scaling, batch correction, ambient RNA removal, or 773 imputation was applied prior to embedding. For non-human datasets (Mouse HSPC), human–mouse 774 one-to-one ortholog mapping was performed to enable compatibility with human -trained scFMs. 775 Time-point annotations and lineage labels were retained as provided in the original publications. 776 EMT: This dataset captures the epithelial –mesenchymal transition (EMT) induced by transforming 777 growth factor beta 1 (TGFB1) in the A549 cell line. Single -cell RNA sequencing was performed at 778 five original time points (0d, 8h, 1d, 3d, and 7d). After quality control, 3,133 cells were retained for 779 downstream analysis. To simplify trajectory modeling and ensure sufficient cell coverage at each 780 time point, we merged the late-stage samples collected at 3 days and 7 days into a single time point, 781 reflecting a shared late EMT state. The resulting four time points were then indexed as 𝑡. to 𝑡[, 782 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint corresponding to 0d , 8h, 1d, and the merged 3d/7d samples . This discretization was used 783 consistently across all trajectory inference and evaluation analyses. 784 Mouse HSPC: This dataset profiles hematopoietic stem and progenitor cell (HSPC) differentiation 785 using a combination of single-cell RNA sequencing and clonal lineage tracing, originally introduced 786 by Weinreb et al. Cells were clonally barcoded using the LARRY (lineage and RNA recovery) system, 787 enabling joint measurement of transcriptional state and lineage relationships across time s. Single-788 cell transcriptomes were collected from heterogeneous progenitor populations at multiple time points 789 during differentiation, both in vitro and following transplantation into mice. The data capture 790 continuous hematopoietic trajectories from multipotent progenitors toward myeloid, erythroid, and 791 lymphoid lineages, providing a well -characterized reference for evaluating trajectory inference and 792 dynamical reconstruction methods. For compatibility with human -pretrained scFMs, one -to-one 793 human–mouse ortholog mapping was applied prior to embedding. After filtering, a total of 49,302 794 cells were included in the analysis. 795 Veres: This dataset profiles in vitro differentiation of human pluripotent stem cells toward pancreatic 796 endocrine lineages, originally reported by Veres et al. Using a high-resolution time-course single-cell 797 RNA-sequencing design, the study captures transcriptional dynamics underlying the emergence of 798 β-cells, α-like poly-hormonal cells, and other pancreatic and non-endocrine populations. Cells were 799 collected across multiple stages of directed differentiation, providing temporally ordered snapshots 800 that span early progenitor states through mature endocrine-like cell identities. For benchmarking, we 801 adopted the original stage and time-point annotations provided in the study, and 18099 cells passed 802 QC. 803 Embryoid bodies (EBdata): This dataset consists of a time -series single -cell RNA -sequencing 804 profiling of human embryonic stem cells differentiated as embryoid bodies over a 27 -day period, 805 capturing the emergence of multiple developmental lineages. Cells were collected at five sampling 806 time points spanning early pluripotent states through diverse differentiated populations, 18,204 cells 807 are retained after QC. 808 HSPC: This dataset characterizes the molecular and cellular dynamics of hematopoietic stem and 809 progenitor cells (HSPCs) during inflammatory response. Single-cell RNA sequencing was performed 810 at three time points (3, 24, and 72 hours) following treatment with the pro -inflammatory cytokine 811 IFNα. After minimal quality contro l, 9983 cells were retained for analysis. Original time -point 812 annotations and any lineage information were used as provided in the study. This dataset captures 813 heterogeneous and dynamic transcriptional responses of HSPCs to IFNα, providing a temporal 814 resolution suitable for benchmarking trajectory inference methods. 815 All datasets were subsequently processed to generate embeddings for six pipelines: five foundation 816 model embeddings (Geneformer, GeneCompass, scGPT, UCE, and scFoundation) and a highly 817 variable gene (HVG)-based baseline. 818 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint

819 1. Jovic, D. et al. Single‐cell RNA sequencing technologies and applications: A brief overview. Clin. Transl. 820 Med. 12, e694 (2022). 821 2. Theodoris, C. V. et al. Transfer learning enables predictions in network biology. Nature 618, 616–624 (2023). 822 3. Liu, Y., Huang, K. & Chen, W. Resolving cellular dynamics using single-cell temporal transcriptomics. Curr. 823 Opin. Biotechnol. 85, 103060 (2024). 824 4. MacLean, A. L., Hong, T. & Nie, Q. Exploring intermediate cell states through the lens of single cells. Curr. 825 Opin. Syst. Biol. 9, 32–41 (2018). 826 5. Moris, N., Pina, C. & Arias, A. M. Transition states and cell fate decisions in epigenetic landscapes. Nat. 827 Rev. Genet. 17, 693–703 (2016). 828 6. Ding, J., Sharon, N. & Bar-Joseph, Z. Temporal modelling using single-cell transcriptomics. Nat. Rev. Genet. 829 23, 355–368 (2022). 830 7. Luecken, M. D. & Theis, F. J. Current best practices in single‐cell RNA‐seq analysis: a tutorial. Mol. Syst. 831 Biol. 15, e8746 (2019). 832 8. Bunne, C., Schiebinger, G., Krause, A., Regev, A. & Cuturi, M. Optimal transport for single-cell and spatial 833 omics. Nat. Rev. Methods Primer 4, 58 (2024). 834 9. Zhao, R. et al. A systematic evaluation of highly variable gene selection methods for single -cell RNA-835 sequencing. Genome Biol. 26, 424 (2025). 836 10. Cui, H. et al. scGPT: toward building a foundation model for single-cell multi-omics using generative AI. 837 Nat. Methods 21, 1470–1480 (2024). 838 11. Yang, X. et al. GeneCompass: deciphering universal gene regulatory mechanisms with a knowledge -839 informed cross-species foundation model. Cell Res. 34, 830–845 (2024). 840 12. Bian, H. et al. scMulan: a multitask generative pre-trained language model for single-cell analysis. Preprint 841 at https://doi.org/10.1101/2024.01.25.577152 (2024). 842 13. Yang, F. et al. scBERT as a large-scale pretrained deep language model for cell type annotation of single-843 cell RNA-seq data. Nat. Mach. Intell. 4, 852–866 (2022). 844 14. Chen, H. et al. Quantized multi -task learning for context -specific representations of gene network 845 dynamics. 2024.08.16.608180 Preprint at https://doi.org/10.1101/2024.08.16.608180 (2024). 846 15. Hao, M. et al. Large-scale foundation model on single-cell transcriptomics. Nat. Methods 21, 1481–1491 847 (2024). 848 16. Wen, H. et al. CellPLM: Pre-training of Cell Language Model Beyond Single Cells. in (2023). 849 17. Rosen, Y. et al. Universal Cell Embeddings: A Foundation Model for Cell Biology. Preprint at 850 https://doi.org/10.1101/2023.11.28.568918 (2023). 851 18. Szałata, A. et al. Transformers in single-cell omics: a review and new perspectives. Nat. Methods 21, 852 1430–1443 (2024). 853 19. Boiarsky, R., Singh, N., Buendia, A., Getz, G. & Sontag, D. A Deep Dive into Single-Cell RNA Sequencing 854 Foundation Models. Preprint at https://doi.org/10.1101/2023.10.19.563100 (2023). 855 20. Wu, J. et al. Biology-driven insights into the power of single-cell foundation models. Genome Biol. 26, 1–856 39 (2025). 857 21. Kedzierska, K. Z., Crawford, L., Amini, A. P. & Lu, A. X. Zero-shot evaluation reveals limitations of single-858 cell foundation models. Genome Biol. 26, 101 (2025). 859 22. Csendes, G., Szalay, K. Z. & Szalai, B. Benchmarking a foundational cell model for post -perturbation 860 RNAseq prediction. bioRxiv 2024.09.30.615843 (2024) doi:10.1101/2024.09.30.615843. 861 23. Wei, Z. et al. Benchmarking algorithms for generalizable single-cell perturbation response prediction. Nat. 862

1–14 (2025) doi:10.1038/s41592-025-02980-0. 863 24. Bendidi, I. et al. Benchmarking Transcriptomics Foundation Models for Perturbation Analysis : one PCA 864 still rules them all. Preprint at https://doi.org/10.48550/arXiv.2410.13956 (2024). 865 25. Ahlmann-Eltze, C., Huber, W. & Anders, S. Deep-learning-based gene perturbation effect prediction does 866 not yet outperform simple linear baselines. Nat. Methods 22, 1657–1661 (2025). 867 26. Cao, J. et al. The single-cell transcriptional landscape of mammalian organogenesis. Nature 566, 496–868 502 (2019). 869 27. Bergen, V., Soldatov, R. A., Kharchenko, P. V. & Theis, F. J. RNA velocity -current challenges and future 870 perspectives. Mol. Syst. Biol. 17, e10282 (2021). 871 28. Gao, M., Qiao, C. & Huang, Y. UniTVelo: temporally unified RNA velocity reinforces single-cell trajectory 872 inference. Nat. Commun. 13, 6586 (2022). 873 29. Bergen, V., Lange, M., Peidli, S., Wolf, F. A. & Theis, F. J. Generalizing RNA velocity to transient cell states 874 through dynamical modeling. Nat. Biotechnol. 38, 1408–1414 (2020). 875 30. La Manno, G. et al. RNA velocity of single cells. Nature 560, 494–498 (2018). 876 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint 31. Schiebinger, G. et al. Optimal-Transport Analysis of Single-Cell Gene Expression Identifies Developmental 877 Trajectories in Reprogramming. Cell 176, 928-943.e22 (2019). 878 32. Lavenant, H., Zhang, S., Kim, Y. -H. & Schiebinger, G. Toward a mathematical theory of trajectory 879 inference. Ann. Appl. Probab. 34, 428–500 (2024). 880 33. Wang, L., Zhang, C. & Zhang, S. Batch Effects Remain a Fundamental Barrier to Universal Embeddings 881 in Single-Cell Foundation Models. 882 34. Klein, D. et al. Mapping cells through time and space with moscot. Nature 638, 1065–1075 (2025). 883 35. Zhang, Z. et al. Deciphering cell-fate trajectories using spatiotemporal single-cell transcriptomic data. Npj 884 Syst. Biol. Appl. 12, 2 (2025). 885 36. Zhang, Z., Sun, Y., Peng, Q., Li, T. & Zhou, P. Integrating Dynamical Systems Modeling with 886 Spatiotemporal scRNA-Seq Data Analysis. Entropy 27, 453 (2025). 887 37. Sha, Y., Qiu, Y., Zhou, P. & Nie, Q. Reconstructing growth and dynamic trajectories from single -cell 888 transcriptomics data. Nat. Mach. Intell. 6, 25–39 (2024). 889 38. Bunne, C. et al. Learning single-cell perturbation responses using neural optimal transport. Nat. Methods 890 20, 1759–1768 (2023). 891 39. Sun, Y., Zhang, Z., Wang, Z., Li, T. & Zhou, P. Variational Regularized Unbalanced Optimal Transport: 892 Single Network, Least Action. arXiv.org https://arxiv.org/abs/2505.11823v2 (2025). 893 40. Gangbo, W., Li, W., Osher, S. & Puthawala, M. Unnormalized optimal transport. J. Comput. Phys. 399, 894 108940 (2019). 895 41. Chizat, L., Peyré, G., Schmitzer, B. & Vialard, F. -X. Unbalanced optimal transport: Dynamic and 896 Kantorovich formulations. J. Funct. Anal. 274, 3090–3123 (2018). 897 42. Chizat, L., Peyré, G., Schmitzer, B. & Vialard, F.-X. An Interpolating Distance Between Optimal Transport 898 and Fisher–Rao Metrics. Found. Comput. Math. 18, 1–44 (2018). 899 43. Dong, X. et al. Data-driven selection of analysis decisions in single -cell RNA-seq trajectory inference. 900 bioRxiv 2023.12.18.572214 (2023) doi:10.1101/2023.12.18.572214. 901 44. Menon, R. et al. Single-cell analysis of progenitor cell dynamics and lineage specification in the human 902 fetal kidney. Dev. Camb. Engl. 145, dev164038 (2018). 903 45. Fischer, D. S. et al. Inferring population dynamics from single-cell RNA-sequencing time series data. Nat. 904 Biotechnol. 37, 461 (2019). 905 46. Dimitrov, D., Schrod, S., Rohbeck, M. & Stegle, O. Interpretation, extrapolation and perturbation of single 906 cells. Nat. Rev. Genet. 1–22 (2026) doi:10.1038/s41576-025-00920-4. 907 47. Liu, Y. et al. TemporalVAE: atlas-assisted temporal mapping of time -series single-cell transcriptomes 908 during embryogenesis. Nat. Cell Biol. 27, 1982–1992 (2025). 909 48. Liarou, M., Matthes, T. & Marchand-Maillet, S. TimeFlow: A Density‐Driven Pseudotime Method for Flow 910 Cytometry Data Analysis. https://doi.org/10.1002/cyto.a.24928 doi:10.1002/cyto.a.24928. 911 49. Ling, Y., Zhang, P., Zhang, Z. & Zhou, P. CellStream: Dynamical Optimal Transport Informed Embeddings 912 for Reconstructing Cellular Trajectories from Snapshots Data. arXiv.org https://arxiv.org/abs/2511.13786v1 913 (2025). 914 50. Gentil, I., Léonard, C. & Ripani, L. About the analogy between optimal transport and minimal entropy. Ann. 915 Fac. Sci. Toulouse Mathématiques 26, 569–600 (2017). 916 51. Léonard, C. A survey of the Schrödinger problem and some of its connections with optimal transport. 917 Discrete Contin. Dyn. Syst. 34, 1533–1574 (2013). 918 52. Zhang, Z., Li, T. & Zhou, P. Learning stochastic dynamics from snapshots through regularized unbalanced 919 optimal transport. in (2024). 920 53. Zhang, Z., Wang, Z., Sun, Y., Li, T. & Zhou, P. Modeling Cell Dynamics and Interactions with Unbalanced 921 Mean Field Schrödinger Bridge. Preprint at https://doi.org/10.48550/arXiv.2505.11197 (2025). 922 54. Zhang, Z., Li, T. & Zhou, P. Learning stochastic dynamics from snapshots through regularized unbalanced 923 optimal transport. Preprint at https://doi.org/10.48550/arXiv.2410.00844 (2025). 924 55. Conneau, A., Lample, G., Ranzato, M., Denoyer, L. & Jégou, H. Word Translation Without Parallel Data. 925 arXiv.org https://arxiv.org/abs/1710.04087v3 (2017). 926 56. Haghverdi, L., Büttner, M., Wolf, F. A., Buettner, F. & Theis, F. J. Diffusion pseudotime robustly 927 reconstructs lineage branching. Nat. Methods 13, 845–848 (2016). 928 57. Qiu, X. et al. Spatiotemporal modeling of molecular holograms. Cell 187, 7351-7373.e61 (2024). 929 58. Boufea, K., Seth, S. & Batada, N. N. scID Uses Discriminant Analysis to Identify Transcriptionally 930 Equivalent Cell Types across Single-Cell RNA-Seq Data with Batch Effect. iScience 23, (2020). 931 59. Korsunsky, I. et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat. Methods 932 16, 1289–1296 (2019). 933 60. Viñas Torné, R. et al. Systema: a framework for evaluating genetic perturbation response prediction 934 beyond systematic variation. Nat. Biotechnol. 1–10 (2025) doi:10.1038/s41587-025-02777-8. 935 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint 61. Liu, T., Li, K., Wang, Y., Li, H. & Zhao, H. Evaluating the Utilities of Foundation Models in Single-cell Data 936 Analysis. Preprint at https://doi.org/10.1101/2023.09.08.555192 (2023). 937 62. Lin, Z. et al. Evolutionary-scale prediction of atomic-level protein structure with a language model. Science 938 https://doi.org/10.1126/science.ade2574 (2023) doi:10.1126/science.ade2574. 939 63. Wolf, F. A., Angerer, P. & Theis, F. J. SCANPY: large -scale single-cell gene expression data analysis. 940 Genome Biol. 19, 15 (2018). 941 64. Benamou, J.-D. & Brenier, Y. A computational fluid mechanics solution to the Monge-Kantorovich mass 942 transfer problem. Numer. Math. 84, 375–393 (2000). 943 65. Kotlov, N. et al. Procrustes is a machine -learning approach that removes cross -platform batch effects 944 from clinical RNA sequencing data. Commun. Biol. 7, 392 (2024). 945 66. Dijksterhuis, G. B. & Gower, J. C. The interpretation of Generalized Procrustes Analysis and allied methods. 946 Food Qual. Prefer. 3, 67–87 (1991). 947 67. Cook, D. P. & Vanderhyden, B. C. Context specificity of the EMT transcriptional response. Nat. Commun. 948 11, 2142 (2020). 949 68. Weinreb, C., Rodriguez -Fraticelli, A., Camargo, F. D. & Klein, A. M. Lineage tracing on transcriptional 950 landscapes links state to fate during differentiation. Science 367, eaaw3381 (2020). 951 69. Veres, A. et al. Charting cellular identity during human in vitro β-cell differentiation. Nature 569, 368–373 952 (2019). 953 70. Moon, K. R. et al. Visualizing structure and transitions in high-dimensional biological data. Nat. Biotechnol. 954 37, 1482–1492 (2019). 955 71. Bouman, B. J. et al. Single-cell time series analysis reveals the dynamics of HSPC response to 956 inflammation. Life Sci. Alliance 7, (2024). 957 958 959 960

961 We thank Professor Weinan E for the helpful discussion. This work was supported by the National 962 Natural Science Foundation of China [Nos. 12288101 and T2321001 to P.Z., Nos. T2350003, 963 12131020, 42450084, 42450135, 12326614, 12426310, and T2542018 to L.C.], Zhejiang Province 964 Vanguard Goose-Leading Initiative [No.2025C01114 to L.C.], Shenzhen Medical Research Fund 965 [Nos. E250200620, E250200621 to L.C.], and Tianfu Jincheng Laboratory [No. TFJCPI20260001 to 966 L.C.]. We acknowledge the support from the High -performance Computing Platform of Peking 967 University and DP technology for computation. 968 969 970 Author contributions 971 P.Z., L.C. and X.Z. conceived the research. X.Z. and P.Z. designed the benchmark. X.Z. and Z.W. 972 performed benchmark and analysis. All authors discussed and interpreted the results. X.Z., Y.L. and 973 Q.T. drafted the manuscript with the input from all authors. P.Z. and L.C. supervised the research. 974 975 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint Competing interests 976 The authors declare no competing interests. 977 978 preprint (which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for thisthis version posted March 12, 2026. ; https://doi.org/10.64898/2026.03.10.710748doi: bioRxiv preprint