Joint Velocity-Growth Flow Matching for Single-Cell Dynamics Modeling
Authors: Dongyi Wang
Yuanwei Jiang
Zhenyi Zhang
Xiang Gu
Peijie Zhou
Jian Sun
Venue: NeurIPS 2025 poster ()
Keywords: flow matching
single-cell
optimal transport
dynamics modeling
Primary area: machine_learning_for_sciences
TL;DR: We propose a novel joint Velocity-Growth Flow Matching (VGFM) approach, driven by semi-relaxed optimal transport, for single-cell dynamics modeling.

Abstract: Learning the underlying dynamics of single cells from snapshot data has gained increasing attention in scientific and machine learning research. The destructive measurement technique and cell proliferation/death result in unpaired and unbalanced data between snapshots, making the learning of the underlying dynamics challenging. In this paper, we propose joint Velocity-Growth Flow Matching (VGFM), a novel paradigm that jointly learns state transition and mass growth of single-cell populations via flow matching. VGFM builds an ideal single-cell dynamics containing velocity of state and growth of mass, driven by a presented two-period dynamic understanding of the static semi-relaxed optimal transport, a mathematical tool that seeks the coupling between unpaired and unbalanced data. To enable practical usage, we approximate the ideal dynamics using neural networks, forming our joint velocity and growth matching framework. A distribution fitting loss is also employed in VGFM to further improve the fitting performance for snapshot data. Extensive experimental results on both synthetic and real datasets demonstrate that VGFM can capture the underlying biological dynamics accounting for mass and state variations over time, outperforming existing approaches for single-cell dynamics modeling.

Joint Velocity-Growth Flow Matching for Single-Cell Dynamics Modeling
Dongyi Wangh,, Yuanwei Jiang",, Zhenyi Zhang?*, Xiang Gu Peijie Zhou Jian Sun
'Scoool of Mathematics and Statistics, Xi an Jiaotong University, Xi an, China
LMMM and School of Mathematical Sciences, Peking University, Beijing, China.
Center for Machine Learning Research, Peking University, Beijing, China {dongyivang, xongyriaangg.j11578887777777 jyw.71117777777718att xjtu. edu. cn; {xianggu, jansun]@xxjuu edu. cn; zhenyizhangdatu, wtw.aangsstuuuuuuuuuuuuuu edu. cn; pj pjzhou0pku. edu.cr cn
Abstract
Learning the underlying dynamics of single cells from snapshot data has gained increasing attention in scientific and machine learning research. The destructive measurement technique and cell proliferation/dcath result in unpaired and unbal- anced data between snapshots, making the learning of the underlying dynamics challenging. In this paper, we propose joint (VGFM), a novel paradigm that jointly learns state transition and mass growth of single-cell populations via flow matching. VGFM builds an ideal single-cell dynamics containing velocity of state and growth of mass, driven by presented two-period dynamic understanding of the static semi-relaxed optimal transport, a mathematical tool that seeks the coupling between unpaired and unbalanced data. To enable practical usage, we approximate the ideal dynamics using neural networks, forming our joint velocity and growth matching framework. A dis- tribution fitting loss is also employed in VGFM to further improve the fitting performance for snapshot data. Extensive experimental results on both synthetic and real datasets demonstrate that VGFM can capture the underlying biological dynamics accounting for mass and state variations over time, outperforming ex- isting approaches for single-cell dynamics modeling. Our code is available at https//ggthub.coo/DDnggg//// Www.//DDnniiaggg6///GGFM. 66/VGGM.
Flow Matching
Introduction
Inferring the latent dynamics of complex systems from sparse and noisy data is a fundamental challenge in science and enginecring. In many domains, e.g., stock markets [1], climate systems [2], and biological processes [3-11], continuous trajectories are rarely fully observed by sensors. Instead, cross-sectional snapshot data collected at discrete time points are commonly provided. This challenge especially important in single-cell RNA sequencing [12-16], where destructive sampling yields unpaired population-level snapshots across time without having the tracked individual cell fates. Furthermore, due to the mass changes during cellular development or response processes, observed data often exhibits mass unbalancedness across time, violating mass conservation. Consequently, reconstructing the time-evolving, unnormalized density function from limited samples has become an important research problem and attracted increasing attention [14, 15, 17].
Deep farming-based models for dynamics inference have demonstrated great potential [3, 7, 18- 29]. These models typically employ ordinary or stochastic differential equations (ODEs or SDEs) parameterized by neural networks to approximate the velocity field governing density evolution.
'These authors contributed equally. 1 'Coresponding authors.
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

One class of methods, based on simulation [3, 18, 19, 23, 24, 26, 30-32], generates synthetic trajectories by feeding initial data through neural networks and solving ODEs or SDEs numerically, and compares simulated results with observations for loss computation. However, these simulation- based methods rely heavily on numerical solvers during training, which significantly increases computational cost. In high-dimensional settings, the enlarged search space further exacerbates instability, hindering scalability and convergence. Another line of research focuses on simulation-Fis free approaches [22, 29, 33-43J, where the velocity field is trained efficiently by constructing conditional probability paths without trajectory simulation, yielding better efficiency and stability in training compared with simulation approaches.
Nevertheless, the simulation-free approaches mentioned above rely solely on the velocity field v, igmoring the unbalancedness of observed data, which can result in incorrect reconstruction of underlying dynamies or unsatisfying generation performance [26, 37, 42]. Distributional imbalance is a common phenomenon in single-cell dynamics, where cellular proliferation and death occur, thus necessitating the introduction of a growth term g to allow for mass increasing or decreasing. To address this, [44, 45] proposed to jointly learn velocity and growth g by minimizing the Wasserstein- Fisher-Rao (WFR) metric [46, 47J. However, it mathematically enforces v = Vg. which lacks clear mechanistic justification in biological systems. Other WFR-inspired approaches [23, 26] design separate neural networks for e and and prioritize fitting the distributions. However, these methods require heavy simulation during training. limiting their scalability to high dimensions.
In this paper. we propose joint Velocity- Growth Flow Matching (VGFM), a novel approach for single-cell dynamics model- ing. VGFM aims to learn the joint state transition" and mass growth' single-cell evolution by flow matching, illustrated in Fig. I. VGFM is based on the static semi-relaxed optimal transport that allows mass variation in building coupling be- tween snapshots, for which we present dynamic understanding that enables mass growth and state transition accomplished in two time periods, respectively. Based on this dynamic understanding, we build a joint transition and growth dynamics, in- heriting the optimal properties of semi-relaxed optimal transport. We then approximate the velocity and growth in the built dynamics using neural networks with finite samples, yielding a joint velocity and growth flow matching framework. To further improve fitting performance for the single-cell snapshots, we employ distribution fitting loss based on the Wasserstein distance. Thanks to the joint velocity and growth flow matching, our approach enjoys more stable training and better scalability to high dimensions, compared with simulation- based approaches. Extensive experiments on both synthetic and real single-cell datasets are conducted to evaluate our proposed approach. Experimental results demonstrate accurate dynamics reconstruction and superior performance achieved by VGFM over recent baselines.
Mass|
Velacity manching
Growth matching
[|Po 8,001
1 () 0.00)
BA(N)
P200)
State Space
Figure 1: The goal of this paper is to learn a joint state transition (controlled by ve (x, t)) and mass growth (con- trolled by 9w (x, t)) dynamics for single-cell evolution.
2 Background and Related Works
This section revisits optimal transport, and reviews the Now-matching-basee methods for single-cell and umbalancedness-aware distribution learning methods, which are mostly related to this paper.
Optimal transport. Optimal transport, tool for transforming one distribution into another at minimal cost [48], has gained remarkable prominence in recent years across fields such as image processing [49-54] and bioinformatics [18, 21, 22, 26, 41]. Given two distributions Po and p1 on domain satisfying Jn Po(x) dr = fg Pi (x) dæ, the optimal transport [48] aims to find a joint distribution of Pu and Pl, named coupling, such that the transport cost is minimized, formulated as
2 State transition descrihes changes gene expression in a cell.
Mass growth means mass increasing or decreasing.

---

the following optimization problem:
[ c(xo, x1) (xo, x1), s.t. PP) = (" 2 0: # = po, P = pi). (1) TEU(Po)PL) Aqa
min
where c is the cost function. For map T, we define Txpo(x) = -a [x:T(x)=xXxx Po dx'. Let
Pi P'(xo x1) = x; for i 0, 1, then we have Pa (x,) Jn क (x0, xi) dx, for i 0,1, When c(xo, X1) = [xx X1 21 by Bremier's theorem [55], the optimal coupling can be expressed as T (Id, ('d =oo where T is called the Monge map.
Flow matching for learning single cell dynamics. Flow matching [33-35] is a simulation-free approach where the velocity field is trained efficiently by constructing conditional probability paths without trajectory simulation. Building on this foundation, further improvements have been achieved by incorporating optimal transport guidance [22, 36], considering unbalancedness [37], accounting for manifold structures [38], and approximating Schrödinger bridge via flow and score matching [41]] all of which has been applied to model single-cell dynamics [22. 37, 38, 41J. Different from these methods, our method explicitly models the simultaneous matching of both the velocity field and the growth function, driven from our developed ideal joint state transition and mass growth dynamics, which does not require the construction of conditional probability paths like the above approaches.
Unbalancedness-aware distribution learning methods. Data imbalance is prevalent across variety of domains, such as image synthesis and protein generation, motivating the development of flow matching models that account for unbalanced distributions [37, 42, 56]. However, these methods often do not model growth functions, which are essential in single-cell contexts. In the single-cell domain, several approaches have been proposed to model growth dynamics [5, 23, 44, 45, 57], yet they either fail to leverage the informative variations in cell abundance observable from snapshot data [44, 45], or rely heavily on computationally expensive simulations [5, 57], limiting their scalability and efficiency. Differently, our method explicitly learns the growth function from observed snapshot data, yielding better performance for single cells than these methods as shown in experiments.
3 Method
Given the snapshot population of single cells, this paper aims to build the dynamic trajectory of single cells. Considering the unbalancedness between snapshot populations due the undergo cell proliferation or death, we aim to build a model that can transform between unbalanced distributions by jointly learning a velocity field for controlling state transition and a growth function for controlling mass variation. Towards this goal, we propose joint Velocity- -Growhh Flow Matching (VGFM) based on semi-relaxed optimal transport to learn the dynamics of single cells. As illustrated in Fig. 2, we first present a decoupled understanding of state transition and mass growth for unbalanced dynamics based on semi-relaxed optimal transport (Figs. (a), (b)), upon which we build a dynamic process between unbalanced distributions (Fig. 2 (c) ). Finally, we propose the velocity and growth flow matching using the built dynamic process to learn the velocity and growth with samples (Fig. 2 (d)). Next, first describe the unbalanced dynamics of single-cell data, and then discuss the details of each component of our method.
3.1 Unbalanced Dynamics of Single Cell
Given two adjacent snapshots of populationsddisributtioss denoted as Po and P1. if the system is balanced, the dynamics (t E [0) between Pu and P1 can be generated by ODE dx, such that Xg ~ Po and X1 P1. Corresponding to this ODE, P. governed by continuity equation O,P, = - V (P.))). However, in biological systems, cells can proliferate and die. Such behaviors can be captured by a dimeddependent weight w, (x,) associated with the state transition, whose evolution is controlled by a growth function g, simulating cell proliferation or death processes:
de = (xx))
dx = (x)),, lo (x,) gt (x:).
dt
Under this formulation, the unbalanced distribution dynamics Pt is generated jointly by the position X₁ and the weight wy (xx). Specifically, wx(x)) models the variation of p: against Pt, i.e,, (x:))
3

---

(d) J oint velocity growth flow matching
Goal
(n) Semi-redaxed optimal transport
Iasptrol
Neural Netwarks
Po Be?
g
Prppotion
Uncssple to Ii:
(c) Joint dynamies
(b) Dynamic understanding
[|pp 1 P,001l
A.p -P- (9.00) AA [=
Therrm 1
Dium grovtul
(State trnellion
4,p. (w/))
(c(A.]
(5.00 8,000l
Figure 2: Illustration of our proposed VGFM, consisting of the velocity and growth flow matching deduced by the dynamic reformulation the semi-relaxed optimal transport.
px(x)//f(xx) for (x() + O. Corresponding to Eq. (2), Pt is governed by
O.p. = V (prv) + 9.Pr.
3.2 Building Unlabalnced Dynamics Based on Semi-Relaxed Optimal Transport
Though Eqs. (2) or (3) provide formulation of unbalanced dynamics between Po and Pt, the equations admit a wide range of admissible dynamics. We next seek the dynamics based on semi- relaxed optimal transport. Specifically, given the semi-relaxed optimal transport
min Toot (T) T20
A
Jgp= c(xo, x1) (xo, x1) + (L (lpo)
subject to
P =
between two unbalanced distributions po and Pl. we will first present dynamic understanding of Eq. (4) that decouples the state transition and mass growth into different periods. Based on this dynamic understanding, we will build reasonable unbalanced dynamics between po and P1. which naturally drives our velocity and growth matching to learn the dynamics with samples (see Sect. 3.3).
Dynamic understanding of semi-relaxed optimal transport. As mentioned above, the solution of the semi-relaxed optimal transport allows both state transition and mass variation/growth. From two-period dynamic process, in which
dynamic perspective, we understand the transport perform mass growth in the first period and state transition in the second period. Note that such an understanding is not practical in applications, but provides decoupled modeling of state transition and mass growth, which promotes a reasonable and friendly-to-matching unbalanced dynamics modeling, driving our flow matching algorithm for learning velocity and growth, as we will make clear later. More concretely, given À E (0, 1), we perform mass growth controlled by 9: for t E [0, A] through 0.P: gePt with initial and ending distributions of po and Pa respectively. For E (X, 1]. the state transition controlled by through å,p, - (pre,) is performed with initial and ending distributions of pa and pt, respectively. We then define the following two-period transport model:
1) So pr (x))) + H(v,, Pt,
min (Me..)) (Po.P))
T (ve, gx) -1 (1-
(5)
where H(vr, Pt) = Ja (x) (x())) ( (x Cx
= ((r, dep: =
gePt, t E [0, Eppt = -V (prv)), E (X, 1]}. and pu and p1 are the given distributions. The following proposition shows the relation between the semi-relaxed optimal transport model in Eq. (4) and the two-period transport model in Eq. (5).
| x1 and if we enforce PL and po share the same support for admissible solution π to problem (4), then we have mill Jsot (T) mine, The (ve, gr Vλ E (0, 1). Moreover, for any E (0, 1), given the optimal transport plan π to problem (4), let P & # n², then x" can be expressed as x- (Id, T 4P where T is the Monge map between pa and P1. Meanwhile, there exist a gt. such that p = po(x)elo ai (x)at, and a of given by _=(x)=x, satisfying (uī (g..)))) ) arg min Jipele (Ve, 9t).
Proposition 1. Assume c(xu,
v/ (T* - x))

---

The proof is provided in the Appendix A.2. Proposition 1 indicates that problems (4) and (5) share the same optimal objective function value. Meanwhile, the optimal objective function value of the two-period transport problem (5) is independent of X. Moreover, ri can be constructed using n²..
Building unbalanced dynamics. Although we have identified a scheme in which pa evolves into pi through two-period process governed separately by et and Ie respectively, this decoupling does not align with the behavior observed in biological systems, where transport and growth typically occur simultaneously. To address this, we next build an unbalanced dynamics based on the dynamic understanding discussed above that allows velocity and growth to jointly drive the state transition and mass growth while ensuring consistency with the target distribution. Specifically, given va (t E (A, 1)). gr. (t E [0, A]) of the two-period process, we define the joint velocity ün and growth 9i by
Te (x) = (1 A) ((2-))(+xx),,
(6)
(x)x = X gh (x))
where E [0, 1], and vh is the flow generated by U, i.e., Ve, = xo + Io Ux(x)ds.
Theorem Given the initial distribution po. denote the ending distribution of the two-period dynamics
a.p. = 9tPt. t E [0, A]; 0.p:= (pev), E (A, 1),
as Pt. and denote the ending distribution of the joint dynamics starting from po
a.p. - (p. grPr. t E 10,
Po Po:
(8)
as Pr. then it holds that P1
P1-
The proof is detailed in Appendix A.3. Theorem indicates that given the same initial distribution at 0, the two-period and the joint dynamics yield the same distribution at t 1. Since 6, Ør are defined on [0, 1). the joint dynamics is more practical in applications, e.g., biological systems. Meanwhile, if tt It is defined in Proposition I, é, ge will inherit the optimality properties of the semi-relaxed optimal transport, benefitting our flow matching algorithm, as shown in Sect. 3.3.
3.3 Velocity and Growth Flow Matching
Since the optimal growth function o in Proposition is required to evolve po to Po # satisfying the marginal constraint po exp (1o² g ds) = pa T possible choice for ge is not unique. For simplicity
log P " (xx log po(x)
and ease of implementation, we choose time-independent form: (x) =
which minimizes the - energy functional of g and can be explained by the Malthusian growth model [58], i.e., the exponential growth model (see Appendix A.4). By plugging the expression of of and consists with from Proposition I and ge into Eq. (6), we obtain
Uz (Wa, (xx)) T (xo) Xo,
Te (Vh, (xo)) log Po n (xx)) log (xx),, where xa Po- (9)
In practice, only finite samples of po and pi can be accessed. We then aim to learn neural networks us vx(xx f). Iw tx,t) approximate o (x), (x) by velocity and growth flow matching as
min Ex, E. [1110 (vr,, t) 0r 8,0
+ low (xo), f) ge (vr.) )|]
(10)
that can generalize to new samples. We first estimate π² and T using samples. Given samples (x) 1 of po and {xt Im of pr. We seek the optimal transport plan x 0 as 't" (Xo. x1) by solving the entropy-regularize
semi-relaxed transport problem using Sinkhorn-based algorithm [59, 60] as
0011 arg min Σ i.j
Ciytiy + (H)) + rKL(T1m)) [w), subject to n 1, Ims
(11)
+20
where Ciý |x x 1m is the all-one vector, (T)) denotes the negative entropy of and and are hyperparameters. With - T can be estimated using barycentric mapping I (x²))
+ E, =., I where Nr E, m Th-11 In implementation, for each x), we sample
5

---

from (1,2,--,,m) with probability Ni
1 m and approximate T (x²) = xi- This approximation is
accurate as c - Therefore, Vi, (xi) R xb [(x] - -xx) - Xx for E [0, 1]. and Eq. (10) becomes
LVGFM((0,0)
0 bj
Ee
[o((x.)) t) - (x] - -x)o|| + [x t) log( 0 1m)i)
(12)
where F]; is the i-th element. The last term Eq. (12) is because po (xb) 1, for 1, n. In
experiments, LVGPM can be implemented using mini-batch samples. Specifically, we sample (2, j) from -0-1 and from (0, 1), then calculate the loss in square brackets. To improve robustness, we add Gaussian noise to x; in experiments.
3.4 Training Process
Our ultimate goal is generate data close to the real data distribution. The current learned ve- locity field of flow matching is the expectation of the conditional velocity field. i.e.,ua(z,t) = Ezl(x.4) f|z). In practice, since we use limited samples learn the velocity field v, each nu- merical integration step introduces an approximation error, which can be accumulated during the numerical ODE solving process, resulting in a deviation from the true trajectory [61]. To address this bias, we employ distribution fitting loss besides the flow matching loss in Eq. (12), to improve performance further. Next, we introduce the distribution fitting loss and our training algorithm.
To achieve better performance, we incorporate the Wasserstein distance between generated samples and observed samples as part of the loss function. Specifically, let X; = [x]] i1
N. for 1 = 0, 1, 1.... T N₂
1 denote the observed samples at different time points, and denote p(X) Ei 5xj- We define on, RRd x T + T7 as the trajectory mapping function parameterized by the neural network uo. which takes given starting point as the initial condition and outputs particle coordinates time indices T according to ODE dynamics, where T denotes a set of time steps. Given an initial set No. the model &ve predicts particle positions at future time points in T. Specifically, the predicted samples are computed by applying the neural ODE to 2g over the time indices {1, T 1}, i.e.
Xr - den (No, (1, T - 1))) Similarly, we define &g., as the particle weight mapping function parameterized by the neural network Iw. which takes the initial weight of particle i as input and outputs the corresponding weight values at T under ODE dynamics. For simplicity of notation, use (x) for x E X denote the weight generated by 6 corresponding to the particles in set & which is generated by Φva The distribution fitting loss is defined the Wasserstein distance as
[oT(O,w)= Σ mi ( ExER ((x) Pio (8i. t=l
1
(
(13)
where Più (&)) Exex, (x)0xx Wi is the W-Wasserstein distance with Euclidean norm.
Total loss and algorithm. Note that the matching loss can be easily generalized multi-time snapshots by employing Eq. (12) among two consecutive snapshots pt. and pt+1- Combined with distribution fitting loss in Eq. (13), we obtain the final loss
[(0,w) = (VCFM(0,w) CoT(O, ())
(14)
During training, we employ parameter scheduling scheme in which we initially use LVGFM as warm-up stage, before enabling both loss terms for joint training. This design provides a notable advantage: after the warm-up stage, and 9w are well-initialized through conditional matching, which facilitates convergence when switching to Wi-based training. The detailed algorithm presented in Algorithm 1. For further details on the warm-up procedure and parameter scheduling, please refer to Appendix B.I.
4 Experiments
We evaluate VGFM on synthetic and real datasets, and compare it with state-of-the-art approaches, including methods that account for distributional unbalancedness [23, 26, 57] as well those that assume mass conservation [22, 38, 41]. Our code will be released online.

---

Algorithm 1:
Training algorithm of joint velocity-grovee growthh flow matching
Input:
Observed data Xo....; Xr-1; warm-up epochs M1; training epochs M2 Output: Trained velocity field vy and growth function I
Compute transport plans + for ,,...,
for i = to M2 do
for = 0 to T' 2 do
Sample batch (i,j) T ml0
Sample ~ 4(0, 1) + to, and x, ~ ~N(xx) + (1 t)x², T2T);
C L vo(xx,,) t) - (xl -x)i|| + (xw(x); - log(( 1-11 'mmli)
if > Mi then
Moot1 t Qvo (dio + 1), utoo (Riot1) Oo (cio to + 1):
C L L Wi
Nopt (Nigt (( Excr, Pix (clott)))
Update 0,w using
Synthetic datasets. Inspired by [23, 26], we adopt the Simulation Gene dataset that applies a three-gene regulatory network to produce a quiescent region and an area exhibiting both transition and growth that can be observed. We also use Dyngen [62] to simulate a scRNA-seq dataset from a dynamic cellular process, which exhibits pronounced branching unbalancedness, with significantly different cell abundances across divergent lineages. Following the experiment setup of [31], we use PHATE [63] to reduce its dimensions to 5. Additionally, inspired by [26] that employs high- dimensional Gaussian mixture model [64] to evaluate the scalability of models, we adopt a more challenging setting: 1000-dimensional Gaussian mixtures.
Real-world dataset. We conduct experiments on three real-world datasets, Embryoid Body (EB) [65], CITE-seq (CITE) [66] and Pancreas [67]. preprocessed following the procedures in [18, 22, 37]. EB and CITE dataset are evaluated under both 5- and 50-dimensional PCA projections, while Pancreas dataset evaluated under 2000 highly variable gene space assess VGFM's scalability on real- world data. Specifically, the EB (5D), CITE (5D), and CITE (50D) configurations are assessed using a hold-out strategy, same as [38], in which intermediate time point excluded during training. The model is then used to predict the distribution at the hold-out time, and we compute the Wi distance between the predicted and true distributions at that time point. The EB (50D) setting is further used for comparison against [26] and ablation study, evaluated using both Wi and the Relative Mass Error (RME) (see below). In the Panereas dataset, in addition to Wi and RME, we also present the mean-variance trends between generated and real gene expressions and analyze the interpretability of the learned growth function g.
Evaluation metrics. We follow [26] to use W1 to measure the distance between the predicted cell distribution and the true cell distribution. Additionally, to assess the accuracy of the growth function Ia. we normalize the total mass at time O to I, and denote the relative mass at time / w.s./ time O N. The predicted relative mass time based on the evolution induced by g denoted as m. = =x (x)/Noo where ExEx, vx()) is the predicted total mass at time t. We define the relative mass error, denoted RME, i.e., RME mull
as my i.e., my No
Implementation details. We employ a 3-layer (5-layer for dimensions greater than 50) MLP with 256 hidden units and LeakyReLU activation to parameterize both the velocity field Uy and the growth function 9w. Optimization is performed using the Adam optimizer with a learning rate of 10 at warm-up stage and 10 4 after applying distribution fitting loss. A warm-up stage of 500 iterations for synthetic datasets and 5000 iterations for real-world dataset is applied (only matching loss in Eq. (12)), ), after which the distribution fitting loss Lor defined in Eq. (13) is applied for an additional 30 training epochs. We will discuss the strategy to set c and in the Appendix B.1.
4.1 Results and Analysis
Ability to reconstruct cellular dynamics. As shown Tab. 1, our model achieves lowest mean Wi and RME on synthetic datasets, outperforming both unbalancedness-aweat methods [23, 26, 57] and
7

---

(a) Predicted dynamics
(b) Predicted growth rate
(c) True growth rate
22: Rico
2.25
Man
Oberarred Pau
sranl
3.75
Vill
3:35
30.
(d) Simulation Gene
(c) Dyngen
(1) EB (50D)
Figure (a) Predicted dynamics with trajectories of VGFM on Simulation Gene. (b) Predicted and (c) true growth rates Simulation Gene. Note that for Simulation Gene, the true growth rates could be accessed. (d), (e) and (f) respectively compare predicted relative mass by UDSB, TIGON, DeepRUOT, and VFGM, and observed relative mass from data, on three datasets. Note that in (f,, on 50D dataset, the mass predicted by TIGON deviates significantly from the observed trend, because of its difficulty in handling high-dimenn nigh-dimensional data as discussed in their paper [23].
Table I: Mean W1 and RME over all time on synthetic datasets. *OT-CFM and OT-MFM do not model the growth function, thus, RME is not computed. "N/C" means "not converge".
Simulation Gene (2D)
Dyngen (5D)
Gaussian (1000D)
Method
Wi
RME
Wi
RME
Wi
RME
OT-CFM* [22]
0.302
3.926
10.126
OT-MFM" [38]
0.311
3.976
11.008
UDSB [57]
0.665
0.192
1.914
0.658
N/C
N/C
TIGON [23]
0.099
0.065
1.029
0.542
N/C
N/C
DeepRUOT [26]
0.068
0.016
0.474
0.199
N/C
N/C
VGFM
0.046
0.006
0.420
0.053
3.010
0.037
umbalancednes--inooiie methods [22, 38], which demonstrate the ability of VGFM to model cellular dynamics for data and mass fitting. We also visualize the learned trajectory and growth function on the Simulation Gene dataset in Figs. 3 (a-c). Figs. 3 (a-c) show that VGFM not only reconstructs accurate trajectories but also produces faithful predictions of spatially varying growth rates.
Scalability and mass-matching performance compared with existing unbalanced models. While existing SOTA approaches of UDSB [57], TIGON [23] and DeepRUOT [26) incorporate unbalanced modeling, these models either introduce biologically implausible transitions [57] (as discussed in (26]) or face difficulties scaling to high-dimensions [23, 26, 57]. As result, these methods do not perform well or converge the 1000-dimensional Gaussian dataset, in Tab. 1. In contrast, VGFM integrates flow matching based on unbalanced transport, allowing it to scale more effectively to high-dimensional datasets. Our method more accurately recovers the ground-truth dynamics defined in [26] and achieves superior mass matssin accuracy, as shown Tab. and Figs. (d-f).
Hold-one-out results on real-world datasets. Table compares hold-one-out results of different methods on EB (5D) and CITE (SD and 50D) datasets. can be observed that VGFM outperforms most existing approaches in terms of Wi distance. We attribute this to the integration of complemen- tary strengths from both flow matching and simulation-basee (Neural ODE) frameworks. On the one hand, flow matching provides robust initialization for learning the velocity field and growth rate,

---

which helps maintain training stability after incorporating the distribution fitting loss Cor. leading to better performance than purely simulation- based models, as also shown in Fig. 4. On the other hand, since the flow matching objective serves only an upper bound on the true Wasserstein distance [68], the introduction of Lor allows direct optimization of the distribution fitting loss, further enhancing model performance. Moreover, by explicitly modeling cell growth through a time-varying weight function, our framework generalizes beyond the conventional setting where Wasserstein distance is computed between unweighted discrete measures, thereby offering greater flexibility in evaluation and enabling more possible transport solutions.
Pancreas dataset under 2000-dimensional gene space. To further explore the scalability of VGFM, we applied our model and compared methods explicitly modeling (x(x) [23, 26] the Pancreas dataset with 2000 highly variable genes. We first observed that VGFM is the only method showing a steadily decreasing training loss, both for LVGFM and LOT, as shown Fig.A-11, 1. Next, following [69]. we calculated the means and variances of the real and generated gene at day 15.5 and plotted the corresponding mean-variance trend and histograms. The results show that the generated samples closely follow the real data, as shown in Fig. A-15. Finally, we analyzed the learned 9u (w, f) and found that our model successfully reconstructs the unbalanced pattern and identifies key genes without being given any prior knowledge of the cell types. The results are reported in Fig. A-16, A-17 and A-18. For more details, please refer to Appendix B.8.
Ablation study on loss terms. To assess the contribution of each loss component in our framework, we perform an ablation study on the EB dataset with unwhitened 50-dimensional PCA features. Ta- ble reports the performance across four time points. VGFM in general achieves the lowest W1 and RME values (excluding time 1) compared with VGFM (w/o LOT) and VGFM (w/o LVGFM), demonstrat- ing the effectiveness of both our joint velocity- growth matching loss and distribution fitting loss. Notably, by removing LVGFM. both W1 and RME increase by significant margin, implying the im- portance of our velocity and growth matching to our framework. Consistent with our motivation in Sect. 3.4, adding in distribution fitting loss Lor re- sults in lower Wi and RME, improved fitting ability to observed distributions.
Table 2: Mean hold-one results on EB and CITE datasets over hold-out times.
EB
CITE
Method
5D
5D
50D
OT-CFM [22]
0.790
0.882
38.756
SF [41]
0.793
0.920
38.524
UDSB [57]
1.206
2.023
44.168
OT-MFM [38]
0.713
0.724
36.394
DeepRUOT [26]
0.774
0.845
38.681
VGFM
0.676
0.745
37.386
(a) Predicted dynamics (w/o Lot) (b) Predicted dynamics
(c) Predicted growth rate
Figure 4: Visualization of predicted dynamics by (a) VFGM (w/ LOT) and (b) VGFM, and growth
rate by (c) VFGM, on EB (5D) dataset, where the hold-out time is the first intermediate timepoint.
Comparison on computational efficiency. We compare the computational cost of our method and the SOTA method DeepRUOT [26] in Tab. 3. For fair comparison, all methods in Tab. 3 adopt the same networks 5-layer MLP with 256 hidden units, and are trained on the same GPU until convergence. As shown in Table 3, our full model requires slightly more training time than the variant using only LVGFM, due to the training with additional distribution fitting loss. However, it still achieves faster convergence compared with simulation-based methods of DeepRUOT [26] and the ablation variant without LVGEM. This demonstrates the critical role of the loss LveFM for facilitating training convergence.

---

To make Lor fully differentiable component in our training pipeline, we also employed the Sinkhorn divergence [70] to replace Wi as the distribution fitting loss Cor. forming a variant of VGFM denoted as VGFM((* in Tab. 3. It can be observed that using the Sinkhorn divergence leads to a further improvement in the accuracy of dynamics reconstruction, while reducing the computational cost for distribution fitting. For more details, please refer to Appendix B.6.
Table 3: Ablation study on EB (50D) dataset. For comparison, we also report the results of the SOTA approach DeepRUOT [26] in this table.
Method
t
f2
t4
Training time
W
RME
W4
RME
W1
RME
Wi
RME
DeepRUOT [26]
8.169
0.416
9.041
0.415
9.348
0.119
9.808
0.296
90 (mins)
VGFM (wfo LoT)
8.915
0.020
10.590
0.098
10.915
0.067
11.635
0.088
(mins)
VGFM (w/ LVGFM)
8.644
0.650
10.167
0.710
11.052
0.823
11.530
0.862
62 (mins)
VGFM
7.951
0.089
8.747
0.042
9.244
0.019
9.620
0.044
13 (mins)
VGFM (*)
7.902
0.018
8.767
0.013
9.063
0.083
9.507
0.096
(mins)
5
Conclusion, Limitation, and Future Work
The paper proposes the joint Velocity- Flow Matching (VGFM) method that jointly learns state transition and mass growth of single-cell populations via flow matching. VGFM designs an ideal dynamics containing velocity field and growth function, driven by our presented two-period dynamic understanding of the static semi-relaxed optimal transport model. Approximating the ideal dynamics using neural networks yields the velocity-growe joint flow matching framework.
Although our approach achieves better scalability and training efficiency compared to simulation- based methods, it is not entirely simulation-fiis free due to the incorporation of the distribution fitting loss Lor- Moreover, since the learned growth rate relies on the observed cell counts each time point, may inadvertently capture effects that are not solely attributable to biological growth, Future work could address them by integrating biological priors into the learning process of the growth function g. and by exploring fully simulation-fis free alternatives that can achieve comparable performance.
Acknowledgment
This work was supported by the National Key R&D Program 2021 YFA 1003002, NSFC (12125104 623B2084, 12501709, 12426313), China National Postdoctoral Program for Innovative Talents (BX20240276), China Fundamental Research Funds for the Central Universities (xzy022025047) and China Postdoctoral Science Foundation (2025M773058).
References
[1] V Gontis, J Ruseckas, and A Kononovičius. A long-range memory stochastic model of the return in financial markets. Physica A: Statistical Mechanics and its Applications, 180(1):000-166, 2010.
[2] Christian LE Franzke, Terence J O'nen Judith Berner, Paul D Williams, and Valerio Lucarini. Stochastic climate theory and modeling. Wiley Interdisciplinary Reviews: Climate Change, 6(1):63-78. 78, 2015.
[3] Grace Hui Ting Yeo, Sachit D Saksena, and David K Gifford. Generative modeling of single-cell time series with prescient enables prediction of cell trajectories with interventions. Nature Communications, 12(1):3222, 2021.
[4] Geoffrey Schiebinger, Jian Shu, Marcin Tabaka, Brian Cleary, Vidya Subramanian, Aryeh Solomon, Joshua Gould, Siyan Liu, Stacie Lin, Peter Berube, et al. Optimal-transport analysis of single-cell gene expression identifies developmental trajectories reprogramming. Cell, 176(4):928-4) 943, 2019.
10

---

[5] Stephen Zhang, Anton Afanassiev, Laura Greenstreet, Tetsuya Matsumoto, and Geoffrey Schiebinger. Optimal transport analysis reveals trajectories in steady-state systems. PLoS Computational Biology, 17(12)x1100446, 2021.
[6] Conrad H Waddington. Canalization of development and the inheritance of acquired characters. Nature, 150(381 1:563111)::36633 1942.
[7] Charlotte Bunne, Stefan Stark, Gabriele Gut, Jacobo Sarabia Del Castillo. Mitch Levesque, Kjong- Van Lehmann, Lucas Pelkmans, Andreas Krause, and Gunnar Rätsch. Learning single- cell perturbation responses using neural optimal transport. Nature Methods, 20111)1159-11768 2023.
[8] James Hay, Lee Kennedy-Shaffer Sanjat Kanjilal, Niall J Lennon, Stacey Gabriel, Marc Lipsitch, and Michael J Mina. Estimating epidemiologic dynamics from cross-sectional viral load distributions. Science, 373(6552))eabh06333 2021.
[9] Jim Oeppen and James W Vaupel. Broken limits to life expectancy. Science, 296(5570):1029- 1031, 2002.
[10] Zhenyi Zhang, Yuhao Sun, Qiangwei Peng, Tiejun Li, and Peijie Zhou. Integrating dynamical systems modeling with spatiotemporal scrna-seq data analysis. Entropy, 27(5), 2025.
[11] Alessandro Palma, Sergei Rybakov, Leon Hetzel, Stephan Günnemann, and Fabian J Theis. Enforeing latent euclidean geometry single-cell VAEs for manifold interpolation. In ICMI., 2025.
[12] Evan Macosko, Anindita Basu, Rahul Satija, James Nemesh, Karthik Shekhar, Melissa Goldman, Itay Tirosh, Allison R Bialas, Nolan Kamitaki, Emily M Martersteck, et al. Highly parallel genome-wide expression profiling of individual cells using nanoliter droplets. Cell, 161(5):1202-122 1214, 2015.
[13] Allon M Klein, Linas Mazutis, Ilke Akartuna, Naren Tallapragada, Adrian Veres, Victor Li, Leonid Peshkin, David A Weitz, and Mare W Kirschner. Droplet barcoding for single-cell transcriptomics applied to embryonic stem cells. Cell, 161(5):1177-100) 2015.
[14] Dragomirka Jovic, Xue Liang, Hua Zeng, Lin Lin, Fengping Xu, and Yonglun Luo. Single-cell ma sequencing technologies and applications: brief overview. Clinical and Translational Medicine, 12(3) 2022.
[15] Intoine-Emmanue Saliba, Alexander J Westermann, Stanislaw A Gorski, and Jörg Vogel. Single-cell rna-seq: advances and future challenges. Nucleic Acids Research, 42(14):8845- 8860, 2014.
[16] Ashraful Haque, Jessica Engel, Sarah A Teichmann, and Tapio Lömnberg. A practical guide to single-cell ma-sequencing for biomedical research and clinical applications. Genome Medicine, 9:1-12, 2017.
[17] David Lähnemann, Johannes Köster, Ewa Szczurek, Davis McCarthy, Stephanie C Hicks, Mark D Robinson, Catalina A Vallejos, Kieran R Campbell, Niko Beerenwinkel, Ahmed Mahfouz, et al. Eleven grand challenges in single-cell data science. Genome Biology, 21:1-35, 2020.
[18] Alexander Tong, Jessie Huang, Guy Wolf, David Van Dijk, and Smita Krishnaswamy. Tra- jectorynet: A dynamic optimal transport network for modeling cellular dynamics. In ICML, 2020.
[19] Charlotte Bunne, Laetitia Papaxanthos, Andreas Krause, and Marco Cuturi. Proximal optimal transport modeling of population dynamics. In AISTATS, 2022.
[20] Zhenyi Zhang, Zihan Wang, Yuhao Sun, Tiejun Li, and Peijie Zhou. Modeling cell dynamics and interactions with unbalanced mean field schrödinger bridge. arXiv preprint arXiv:2505. / 2025.
11

---

[21] Yuhao Sun, Zhenyi Zhang. Zihan Wang, Tiejun Li, and Peijie Zhou. Variational regularized unbalanced optimal transport: Single network, least action. arXiv preprint arXiv:2505./1823, 2025.
[22] Alexander Tong, Kilian Fatras, Nikolay Malkin, Guillaume Huguet, Yanlei Zhang, Jarrid Rector- Brooks, Guy Wolf, and Yoshua Bengio. Improving and generalizing flow-based generative models with minibatch optimal transport. Transactions on Machine Learning Research, pages 1-34, 2024.
[23] Yutong Sha, Yuchi Qiu, Peijie Zhou, and Qing Nie. Reconstructing growth and dynamic trajectories from single-cell transcriptomics data. Nature Machine Intelligence, 6(1):25-39, 2024.
[24] Qi Jiang and Lin Wan, A physics-informed neural sde network for learning cellular dynamics from time-series serna-seq data. Bioinformatics, 40:11120-i112) 2024.
[25] Antonio Terpin, Nicolas Lanzetti, Martín Gadea, and Florian Dorfler. Learning diffusion at lightspeed. In NeurlPS, 2024.
[26] Zhenyi Zhang, Tiejun Li, and Peijie Zhou. Learning stochastic dynamics from snapshots through regularized unbalanced optimal transport. In ICLR, 2025.
[27] Xi Zhang, Yuan Pu, Yuki Kawamura, Andrew Loza, Yoshua Bengio, Dennis Shung, and Alexander Tong. Trajectory flow matching with applications to clinical time series modelling. In NeurlPS, 2024.
[28] Anming Gu, Edward Chien, and Kristjan Greenewald. Partially observed trajectory inference using optimal transport and dynamics prior. In ICLR, 2025.
[29] Martin Rohbeck, Charlotte Bunne, Edward De Brouwer, Jan-Christian Huetter, Anne Biton, Kelvin Y. Chen, Aviv Regev, and Romain Lopez. Modeling complex system dynamics with flow matching across time and conditions. In ICLR, 2025.
[30] Ricky T. Q. Chen, Yulia Rubanova, Jesse Bettencourt, and David Duvenaud. Neural ordinary differential equations. In NeurIPS, 2018
[31] Guillaume Huguet, Daniel Sumner Magruder, Alexander Tong, Oluwadamilola Fasina, Manik Kuchroo, Guy Wolf, and Smita Krishnaswany, Manifold interpolating optimal-transport flows for trajectory inference. NewrIPS, 2022.
[32] Takeshi Koshizuka and Issei Sato. Neural lagrangian schrödinger bridge: Diffusion modeling for population dynamies. In ICLR, 2023.
[33] Yaron Lipman, Ricky TQ Chen, Heli Ben-Hamu, Maximilian Nickel, and Matthew Le. Flow matching for generative modeling. In ICLR, 2023.
[34] Xingchao Liu, Chengyue Gong, and qiang liu. Flow straight and fast: Learning to generate and transfer data with rectified flow. In ICLR, 2023.
[35] Michael S Albergo and Eric Vanden-Eijnden Building normalizing flows with stochastic interpolants. In ICIR, 2023.
[36] Aram-Alexandre Pooladian, Heli Ben-Hamu, Carles Domingo-Fanich, Brandon Amos, Yaron Lipman, and Ricky TQ Chen. Multisample flow matching: Straightening flows with minibatch couplings. In ICLR, 2023.
[37] Luca Eyring, Dominik Klein, Théo Uscidda, Giovanni Palla, Niki Kilbertus, Zeynep Akata, and Fabian Theis. Unbalancedness in neural monge maps improves unpaired domain translation. In ICLR, 2024.
[38] Kacper Kapusniak, Peter Potaptchik, Teodora Reu, Leo Zhang, Alexander Tong, Michael Bron- stein, Joey Bose, and Francesco Di Giovanni. Metric flow matching for smooth interpolations on the data manifold. In NeurIPS, 2024.
12

---

[39] Lazar Atanackovic, Xi Zhang, Brandon Amos, Mathieu Blanchette, Leo J Lee. Yoshua Bengio, Alexander Tong, and Kirill Neklyudov. Meta flow matching: Integrating vector fields on the wasserstein manifold. In ICMI, 2024 Workshop on Geometry- Representation Learning and Generative Modeling, 2024.
[40] Dominik Klein, Théo Uscidaa, Fabian J Theis, and marco cuturi. GENOT: Entropic (gromov) wasserstein flow matching with applications to single-cell genomics. In NeurlPS, 2024.
[41] Alexander Y Tong, Nikolay Malkin, Kilian Fatras, Lazar Atanackovic, Yanlei Zhang, Guillaume Huguet, Guy Wolf, and Yoshua Bengio. Simulation-free schrödinger bridges via score and flow matching. In AISTATS, 2024.
[42] Gabriele Corso, Vignesh Ram Somnath, Noah Getz, Regina Barzilay, Tommi Jaakkola, and Andreas Krause. Composing unbalanced flows for flexible docking and relaxation. ICLR, 2025.
[43] Alessandro Palma, Till Richter, Hanyi Zhang, Manuel Lubetzki, Alexander Tong, Andrea Dittadi, and Fabian Theis. Multi-modal and multi-attriihn generation of single cells with cfgen. In ICLR, 2025.
[44] Kirill Neklyudov, Rob Brekelmans, Daniel Severo, and Alireza Makhzani. Action matching: Learning stochastic dynamics from samples. In ICML, 2023.
[45] Kirill Neklyudov, Rob Brekelmans, Alexander Tong, Lazar Atanackovic, qiang liu, and Alireza Makhzani. A computational framework for solving wasserstein lagrangian flows. In ICML, 2024.
[46] Lenaic Chizat, Gabriel Peyré, Bernhard Schmitzer, and François-Xavics Vialard. Unbalanced optimal transport: Dynamic and kantorovich formulations. Journal of Functional Analysis, (74(11):3090-3) 2018.
[47] Matthias Liero, Alexander Mielke, and Giuseppe Savaré. Optimal entropy-transport problems and new hellingger-knttooovic distance between positive measures. Inventiones Mathematicae, 21 111(3):99-1117, 2018.
[48] Leonid Kantorovich. On the translocation of masses. In Dokl. Akad. Nauk. USSR (NS), 1942.
[49] Xiaole Tang, Xin Hu, Xiang Gu, and Jian Sun. Residuul-conditiinned optimal transport: towards structure-reeeerving unpaired and paired image restoration. In ICML, 2024.
[50] Xiaole Tang, Xiang Gu, Xiaoyi He, Xin Hu, and Jian Sun. Degradation-aware residual- conditioned optimal transport for unified image restoration. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2025.
[51] Xiang Gu, Liwei Yang. Jian Sun, and Zongben Xu. Optimal transport-guided conditional score-based diffusion model. In NeurIPS, 2023.
[52] Nikita Gushchin, Alexander Kolesov, Alexander Korotin, Dmitry P. Vetrov, and Evgeny Burnaev. Entropic neural optimal transport via diffusion processes. In NewrIPS, 2023.
[53] Alexander Korotin, Daniil Selikhanovych, and Evgeny Burnaev. Neural optimal transport. In The Eleventh International Conference on Learning Representations, 2023.
[54] Taoran Zheng, Xing Li, Yan Yang, Xiang Gu, Zongben Xu, and Jian Sun. Towards prospective medical image reconstruction via inowledge-informed dynamic optimal transport. arXiv preprint arXiv:2505. 2025.
[55] Yann Brenier. Décomposition polaire et réarrangement monotone des champs de vecteurs. CR Acad. Sci. Paris Sér. I Math., 305:805-808, 1987.
[56] Jaemoo Choi and Jaewoong Choi. Scalable simulation-fire entropic unbalanced optimal transport. arXiv preprint arXiv:4410002656, 2024.
[57] Matteo Pariset, Ya-Ping Hsieh, Charlotte Bunne, Andreas Krause, and Valentin De Bortoli. Unbalanced diffusion schrödinger bridge. arXiv preprint (tti::3306000999 2023.
13

---

[58]
Thomas Robert Malthus, An Essay on the Principle of Population. 1798.
[59] Marco Cuturi. Sinkhorn distances: Lightspeed computation of optimal transport. In NeurIPS,
2013.
[60] Gabriel Peyré, Marco Cuturi, et al. Computational optimal transport: With applicatiios to data science. Foundations and Trends@ in Machine Learning, 11(5-6):3556677, 2019.
[61] Kevin Frans, Danijar Hafner, Sergey Levine, and Pieter Abbeel. One step diffusion via shorteut models. In ICLR, 2025.
[62] Robrecht Cannoodt, Wouter Saelens, Louise Deconinck, and Yvan Saeys. Spearheading future omics analyses using dyngen, multi-mode simulator of single cells. Nature Communicattonss, 12(1):3942, 2021.
[63] Kevin R Moon, David Van Dijk, Zheng Wang, Scott Gigante, Daniel B Burkhardt, William S Chen, Kristina Yim, Antonia van den Elzen, Matthew J Hirn, Ronald Coifman, al. Visu- alizing structure and transitions in high- nigh-dimensional biological data. Nature Biotechnology. 37(12):1482-1492 2019.
[64] Lars Ruthotto, Stanley J Osher, Wuchen Li, Levon Nurbekyan, and Samy Wu Fung. A machine learning framework for solving high-dimensional mean field game and mean field control problems. Proceedings of the National Academy of Sciences, 17(17):1133-919, 2020.
[65] Kevin R Moon, Jay S Stanley III, Daniel Burkhardt, David van Dijk, Guy Wolf, and Smita Krishnaswamy, Manifold learning-basd based methods for analyzing single-cell ma-sequeneing data. Current Opinion in Systems Biology, 7:36-48 46, 2018.
[66] Christopher Lance, Malte D) Luecken, Daniel B Burkhardt, Robrecht Cannoodt, Pia Rauten- strauch, Anna Laddach, Aidyn Ubingazhibov, Zhi-Jie Cao, Kaiwen Deng, Sumeer Khan, et al. Multimodal single cell data integration challenge: results and lessons learned. BioRxiv, pages 2022-04, 2022.
[67] Aimée Bastidas- Ponce,, Sophie Tritschler, Leander Dony, Katharina Scheibner, Marta Tarquis- Medina, Ciro Salinno, Silvia Schirge, Ingo Burtscher, Anika Böttcher, Fabian J Theis, et al. Comprehensive single cell mrna profiling reveals a detailed roadmap for pancreatic endocrino- genesis. Development, 146(12):dev 1et.vvi178848 173849, 2019.
[68] Joe Benton, George Deligiannidis, and Arnau Doucet, Error bounds for flow matching methods. Transactions on Machine Learning Research, 2024.
[69] Alessandro Palma, Till Richter, Hanyi Zhang, Manuel Lubetzki, Alexander Tong, Andrea Dittadi, and Fabian Theis. Multi-modal and multi-attribute generation of single cells with efgen. In ICLR, 2025.
[70] Jean Feydy, Thibault Séjourné, François- François-Xavier Xavier Vialard, Shun-ichi Amari, Alain Trouve, and Gabriel Peyré. Interpolating between optimal transport and mmd using sinkhorn divergences. In AISTATS, 2019.
[71] Jean-David Benamou and Yann Brenier. A computational fluid mechanics solution to the monge-kantorovich mass transfer problem. Numerische Mathematik, 84(3):375-393, 2000.
[72] Alessio Figalli and Federico Glaudo. An invitation to optimal transport, Wasserstein distances, and gradient flows. EMS Press, 2021.
[73] Caleb Weinreb, Alejo Rodriguez- eage tracing on transcriptional landscapes links state to fate during differentiation. Science, 367764479)caaww338,, 2020.
Rodriguez-Fraticelln Fraticelli, Fernando D Camargo, and Allon M Klein. Lin-
14

---

Outline of Appendix
The appendix is organized into three main parts.
Appendix A presents the theoretical foundations. We begin by reviewing the Brenier-Benamou theorem [71]. followed by detailed proofs of Proposition I and Theorem 1 that are stated in the main text, as well as the justification for our choice of the reparameterized growth function Ò in Eq. (9).
Appendix B focuses experiments. We first describe general implementation setup, including the strategy for selecting the hyperparameters in Eq. (11). Next, we provide detailed settings, visualizations and analysis for each experiment.
Appendix C discusses the broader impacts of our work.
A Proofs
A.1 Background: Brenier- Bename Formulation
We first recall the dynamic formulation of balanced optimal transport, which is known as Brenier- Benamou formula [71]. In this paper, we assume po and P1 are absolutely continuous w.E.t Lebesgue measure, which omitted for convenience of description.
Theorem A-2 (Brenier-Be Benamo- formula). Given two probability measures Po,P1 E P2().. it holds that
W pı) = inf { (x)(())) dxdt
0.p: dxdet
- (pec)), Pi=u Po: Pt=1 = P1
Peite
(A-15)
and the optimal of (x) can be expressed by Monge map between po and pr. ie,
2x + (x) x)) = Tx) x, x
(A-16)
A.2 Proof of Proposition 1
For the convenience of reading, here we restate Proposition in the main text as Proposition A-2. Proposition A-2. Assume x(x,,x1) ||xo - x1I1² and if we enforce Pa== and po share the same support for admissible solution æ to problem (4), then we have min, [so min Jipo (Ve, gt), AX (0,1). Moreover, for any (0, 1). given the optimal transport plan n" to problem (4), let p po n", then x" can be expressed as n-" = (Id, T ) where T is the Monge map between p} and p1. Meanwhile, there exist S. such that p3 = (x(xee of (x)dt, and et given by
t 0i (x+1=>(r(<x) (TM((x) x))
(x)x) (x)
T
(A-17)
satisfying (vt,gi) E arg min Tigh (Ve,
Proof. Recall that when c(Xu, = ||x0 11 problems (4) and (5) are respectively
min Tsot (Tr A / ||xo - x1112 dT(xx, x1) + (L # (ppp) subject #20
P1 = P1,
Jgz
and
(pp (1-))) [/ (xxx)]loxxx
min (EAx).PPPP
where H(V,, gt, Pt) = Jn Po (x) (xx)de 9e (x(x)tt11++yy) Cx = {(v,, )= dep: =
DePr. t [0, A); O.pt V (p(v)), (A,1))
(1) We first prove min,, (TT) > min T (vs gt).
15

---

Let nT" be the optimal solution to problem (4). For any X E (0, 1), define
= # T,, and let
wx(x)= = Px(x)
Define the growth as
Po(x)
d
log wx(xx = go(x), E [].
dt
This implies the continuity equation Oept = SePr for t E [0, AJ, and hence we have Po(x) exp (Io (x()))))) = P ...
By Brenier's theorem [60], we have xt = (Id, #P> where (x) is the Monge map from px Define
gx(x) log Po - log (x)
0x 1- 1 (T x)) (x)
-->
(A-18)
1-
Plugging into the dynamic objective and using the definition of KL-divergence:
KT(PP)lpo) = 1/8 po(x)
(xxxx) po(xx Pu(x) -1 1) + +1) dx.
(log (x(x)
we obtain
(vr: g1) J. po(x)) afó se(x))ee oaxi) -1 1) +1) dx + (1 (2-2) / (T' -x112x(xxddh
-KI(PAllRo) + I ||xo x1ll s1l (Xo. x1)
-Toot(TT))
Thus, (T=) = Tipe (Vt, 9t) 2 mine... Tupt (Ve, 9t).
(2) To show the reverse inequality ming Jsoe S mine. Tepe (vt, 9t), we assume the contrary that min, (T) > min (v 9c). Let (vt, 98) be optimal solution problem (5), and define
(x)x==)oxx) (x) = po(x) exp (f. (f
Since the second- stage evolution Pt, for E (A, 1], is governed solely by the velocity field and does not involve mass creation or destruction, both PA and p1 have the same total mass. Thus, the Monge map 7 from px to P1 under quadratic cost is well-defined. Then, by the Benamou-Brenier formulation [71, 72], we have
> (T (T* (
x)) 7x((x)
of (x -
1-
and the corresponding coupling # = (Id, T 4DA satisfies
Toot ( -KL(px + / ||xo 1x-xx1l1 xxxx (xo, x1)
f po(x)
edo' of (x))) (x)a(x)) +1 dx + (1 (1-2))p
12
= JTe JTieeei (e g6)
min The (v. < min Toot(T)
ve.ge
which leads the contradiction.
Combining (1) and (2), we have Jnnt (x") = Thd (uf, gf))
Given T", we construct 9i, Vt as in Eq. (A-18), According to the proof of (1), we have J (VZ,gt) (ve,
Toot (T²) min.,.g gz J (v, (c). gc).
16

---

A.3 Proof and Empirical Evidence of Theorem
For the convenience of reading, here we restate Theorem in the main text as Theorem A-3. Theorem A-3. Gvien the initial distribution Po, denote the ending distribution of the two-period dynamics
Oep: 9.Pr, t E [O; dep: = - (pen)), E (2,
(A-19)
as Pt, and denote the ending distribution of the joint dynamics starting from po
O.p: - (p(o.) gcPt. t [0,0,]],
Pa Po.
(A-20)
as Pr. then it holds that P1
P1-
Proof. Given Xo ~ Po, consider the two systems below.
System I (original two-period transport dynamics);
dx:
= vx(xx) (2))
dt
(A-21)
d
log (x(xxx) (xx) (t),
dt
where we the time-dependent weight function and In is the indicator function, for any function
f. If E Q,
f.f TT(() = o, Otherwise.
System II (joint dynamics defined via reparameterization (6)):
dx: dt
(x()));
(A-22)
log (x,)) 9o (x+):
dt
with w, being the time-dependent weight under joint dynamics. We first recall the definition (x(x) (1 X) Ø(1 P(1=))++a(X) (x), = X - Ixe (v,,(())) To prove Theorem 1, it suffices to show that given the same initialization Xo, we (xo). the final state of system I (X1, wi (x1)) and II (x), wi (x)) are identical, i.e., x1 x and ūū (X1) w (xi).
x X0 +
(x) dt
xo +
0i1 0(1-1je-1 111-1111111 (1- A) dt
xx + / Us (x) ds (let (1 A)t + A)
X1-
Meanwhile,
log wi (x1) = log wi (xx) = log wx (Xo)
1' (x(x)) ds
log wo (xo)
/-² A gxe(xx) dt (lett= 1
= loo wx(o) +
x
log (o(()) + / x exe( (x= xo + ,'' (x)dss
f (x) dt
log wo(x))
= log úi (x1).
Hence, the final state and mass under both systems coincide. Applying the above analysis to all Xg po completes the proof.
17

---

Empirical evidence of Theorem 1. We provide the following example an evidence to Theorem
Set Po (2, 0.5), ex(x) = 2t, gx(x) log(x + 1) + ta, using (q.( we have (,(x) (1 -)) (x(x) ( (x)) Agxe ((1 Ajt + A)² + X²). By setting
0.4 and using numerical solvers solve the corresponding continuity equation, we validate our correctness of our theorem. As shown in Fig. A-5, the two dynamics ended in the same distribution.
3.-
35-
00-
P1
25-
34-
04-
0.4
635
511
0.1
30
50
60
(a) Dynamics geverned by Vo.rr
(b) Dynamics governed by 9:.
(c) D..Ps
0=
Figure A-5: Empirical evidence of Theorem 1.(a) Original two-period transport dynamics, (b) Joint dynamics defined via reparametcrization 6. (c) The two dynamics ended the same distribution.
A.4 Motivation for time-independent growth function J. according to Eq. (9)
Minimization of 1.2 energy potential. We aim prove the following statement: Given function (x) whose logarithmic derivative satisfies de log wx(x) gx(x) and the boundary condition
log wi (x) - log wo (x) = S (x) dt,
then among all functions ge (x) satisfying this constraint, the one minimizing the energy functional
E(g) = Sa 11 (x) dfdx
is the constant function fo 9r (x)dx = log w1 (x) log (x). To show this, we apply the method of Lagrange multipliers. Introduce a multiplier X(x) and consider the Lagrangian
C(ge) = - S.S (gx(x)² + A(x) 9x))) di X(x) (log wi (x) log wo (x) (x)x)xx
Taking the first variation of with respect to gx(x) yields the optimality condition
8C (x) 2xx(x) + A(x) = + gx(x)
A(x)
5g:
This shows that the optimal gr (x) is constant with respect to time. Plugging this into the constraint gives 9r (x) = log wh (x) - log wp (x). Hence, the constant function ge (x() is the unique minimizer of the energy functional under the given constraint.
Explanation from the Malthusian growth model [58], i.e. the exponential growth model. The Malthusian growth model is dp:
dt = gpt, where pt the population size and is a constant growth åp: rate. Our model on the growth ät JrPt is consistent with the above Malthusian growth model. In these models, g is treated as a constant based on the assumption that the resources are abundant and the environment is stable, causing growth rates to remain relatively stable over time. In the context of scRNA-seq experiments, these conditions are typically satisfied because cells are often cultured or sampled under controlled laboratory conditions, where nutrient supply, temperature, and other environmental factors are maintained at constant levels. Thus, choosing this form is biologically reasonable, especially when there is no prior knowledge about the growth rate.
18

---

100
10
312
166
20
10
180
102
079
(a) Simulation Gene data
(b) Dygen data
(c) Gaussian data
Figure A-6: Fixing 6, we plot the transport cost as a function of r. (a) In Simulation Gene, the curve flattens around 5; (b) In Dyngen, the transport cost also stabilizes near 5.
B Experimental Details
The model architecture has been described in Sect. 4. We next detail the strategy for selecting hyperparameters e and r in Eq. (11), as well as provide a comprehensive description of each experiment. All experiments performed single-core CPU without GPU acceleration, and all visualizations are based on projections onto the first two dimensions the high-dimensional data.
B.1 Determine the Entropy Regularization Parameter and Relaxation Coefficient
Our algorithm involves solving a static semi-relaxed optimal transport problem (11), which includes two critical hyperparameters and These parameters play pivotal role in the behavior and stability of the model. one hand, if is too small, it may lead to numerical instability; it too large, cause incorrect cross-branch matching in unbalanced data scenarios. To select e, first set to a moderately large value, such as т = 50, and then perform a grid search over increasing values of starting from 0 until numerical stability is achieved. The choice of is even more crucial. When r is too small, the KI. divergence term between P0 π and Po receives a weak penalty. As a result, transport mass becomes overly concentrated around points close to Pt, leading to highly uneven marginals 1m in the discrete setting. This causes instability in the learning of the growth function and may result erroneous modeling. Conversely, as > +o0, Eq. (11) reduces to balanced problem, which contradicts the unbalanced nature of our formulation. To determine suitable value for T, we fix the previously selected small and gradually increase T while observing the variation in Σ TiyCij Cij. This curve typically exhibits increasing then-flaee flattening trend. Similar the "elbow metho" used in clustering to select the optimal number of clusters, we identify the region where the curve becomes stable and choose accordingly, illustrated in Fig. A-6.
On the simulation gene dataset, we fix c 0.001 and evaluate the model with different values of {5, 10, 15, 20). all of which lie within the identified stable region. The corresponding model performance metrics are summarized in Tab. A-4. We observe that the models trained with 7 this range are stable, validating the effectiveness of our hyperparameter selection strategy.
Table A-4: Fix e 0.001, change different 7 on simulation gene data, we present the Wasserstein- distance of each timepoint.
t
t2
t4
Models
W1
RME
Wi
RME
Wi
RME
Wi
RME
VFGM (r 5)
0.046
0.007
0.062
0.001
0.053
0.003
0.063
0.003
VFGM (T 10)
0.041
0.005
0.053
0.005
0.038
0.007
0.040
0.008
VFGM (T -15)
0.045
0.003
0.056
0.003
0.046
0.005
0.052
0.004
VFGM (T =00)
0.043
0.007
0.057
0.004
0.045
0.007
0.058
0.011
19

---

B.2 Simulation of Gene Data
Following the setup of [26], the dynamics of the simulated gene regulatory network are governed by the following system of differential equations:
dx
04 x + B
81X1 7/5..
dt
1 X X2 + 13X3 +
dX2
X2 + B
82X₂ + 72
(A-23)
dt
1+ + 2 + 73X3 +
dX3
03X3
1+ 23x3 63X3 + naat:
dt
where (t) denotes the expression level of gene at time t. Genes X1 and X₂ mutually inhibit each other while self-activating, forming toggle-switch regulatory motif. An external signal B independently activates both X3 and X2, whereas Xg inhibits the expression of both X, and X2.
The parameters and "/ determine the strengths of self-activation and inhibition, respectively. The ó, terms denote degradation rates, and TıC represents additive Gaussian noise with intensity y/:- The cell division rate is positively correlated with X2 expression and is calculated as
X2 9 1+ X3 %.
At each cell division event, daughter cells inherit gene expression values (X1; X2) parent, subject to small perturbations yaN 1) per gene. Post-division, cells evolve independently according to the same stochastic dynamics.
X3) from the
Initial Conditions and Simulation Setup. Initial expression states are sampled from two normal distributions: ([2, 0.2, 0], 0.01) and A((O, 0, 2), 0.01). every step, negative values are clipped to zero. Gene expression data is recorded at time points / {0, 8, 16, 24,32).
By setting 10 and 0.003 according
selection scheme described in Appendix B.1 and illustrated Fig. A-6 (a). Our model achieve the best performance in terms of (weighted) Wi and relative mass error shown in Tab.
Table A-5: Simulation parameters on gene regulatory network [26].
Parameter
Value
Description
0.5
Strength of self-activation for X
0.5
Inhibition of X₂ by Xi
012
Strength of self-activation for X₂
7/2
Inhibition of X, by X₂
03
Strength of self-activation for X3 Half-saturation constant for inhibition Degradation rate for X1 Degradation rate for X₂ Degradation rate for X
73
81
0.4
62
0.4
63
0.4
1/1
0.05
Noise intensity for X₁ Noise intensity for X₂
172
0.05
0.05
Noise intensity for X3
1jd
0.014
Noise for perturbation during cell division
B
External activation signal Simulation time step Observation time points
dt
Time Points
{0, 8, 16, 24, 32}
Growth rate correlation analysis.
To validate the accuracy of our growth rate modeling and learning, we perform correlation analysis between the predicted growth rates and the ground truth. We utilize this dataset, which is generated
20

---

from the prescribed dynamical systems (Eq. (A-23). Since the ground truth growth rates can be computed analytically, we conduct correlation analysis on out-of-distribution (OOD) time points excluded from the training data assess the generalization capability of VGFM. The results are summarized in Table A-6.
Table A-6: Correlation analysis four OOD time points for Simulation Gene Data
OOD Time Point Pearson Correlation Coefficient
$0.5
0,980
$1.5
0,992
$2.5
0.995
f3.5
0,996
The growth rates predicted by our VGFM exhibit strong correlation with the ground truth values.
Uncertainty quantification. Although our current model deterministic, we incorporate dropout into the velocity and growth networks to quantify predictive uncertainty. During inference, compute the average variance for each dimension considering different branches (0 denotes the quiescent region i.e., the lower left corner of Fig. while denotes the region where exhibits mass variation and state transition) and time points. We conduct the experiment with dropout rate of 0.1, performing 5 stochastic forward passes through the velocity and growth networks for all data points. The results are presented in Table A-7.
velocity and growth across all data points,
Table A-7: Uncertainty estimation of predicted velocity (2-dimensional) and growth on Simulation Gene Data. var,, and varg, represent the average variance of the first and second velocity dimensions, respectively, while var, denotes the variance of growth. All variance values scaled by 10
Branch Time
varx, (x10~
varx, (x10~>)
varg (x10-3)
0.1080
0.1000
0.0550
0
0.0770
0,0600
0.0660
0
0.0660
0.0670
0.1100
0
0,0850
0,0820
0.1910
0
0.1220
0.1370
0,3500
1
2.9320
6,5010
0.2510
1
3.0230
4,2870
0.3060
1
2
2,0410
1.9780
0.5590
1
0,9530
0.9760
0.7660
0.2800
0.3100
1.1680
We analyze the uncertainty patterns from both spatial and temporal perspectives:
For spatial analysis, cells in branch 0 exhibit minimal movement with little variation in state and counts. Correspondingly, the variances of ei, and in this branch are consistently small, indicating low predictive uncertainty. In contrast, branch demonstrates substantially larger variances, reflecting higher uncertainty in predictions.
For temporal analysis, The Simulation Gene Data features significant state transitions and quantity changes during initial phases, with diminishing magnitudes over time. The variance patterns of vh, U2, and y accurately capture this temporal evolution, validating our uncertainty estimation approach. Specifically, for branch 1, the of u decrease monotonically over time, aligning with the reduced dynamics in later stages.
Effect of sampling bias. In cellular dynamics reconstruction, conventionally assumed that observed datasets accurately capture the temporal evolution of the true cell distribution. Consequently, the observed cell counts should align with the biological processes of cell proliferation and death. Evaluating the robustness of VGFM when trained biased observational data holds significant practical relevance.
To generate non-uniform or subsampled data, we resample the original dataset with perturbed branch ratios. Specifically, for each time point I, we randomly sample perturbation factor a, ~ NN(0, Given original branch ratio of s, st), the resampled ratio becomes (st -sal)))) :1 (1- s,+ (..).
21

---

where quantifies the bias level. We train our model on the resampled data and evaluate its performance on the original dataset to assess VGFM's robustness.
Table A-8: Performance evaluation (W/RMM metrics) on Simulation Gene Data with varying sampling bias levels, averaged over three random seeds.
Bias Factor o
t2
f3
4
0 (unbiased)
0.041/0.007
0.053/0.005
0.038/0.007
0.040/0.008
0.1
0.074/0.035
0.094/0.051
0.075/0.063
0.077/0.032
0.2
0.085/0.059
0.113/0.073
0.134/0.062
0.123/0. 125
B.3 Dyngen Data
Dyngen is multi-modal simulation engine for studying dynamic cellular processes at single-cell resolution [62]. We follow the setup of [31], which simulates gene expression time-series data that mimics cell proliferation processes, including branching and temporal progression. In the analyzed instance of the dyngen dataset, the samples span five discrete time points, labeled from 0 to 4, capturing the temporal evolution of the cell population.
We use PHATE [63] to reduce its dimensions to 5, the same as [31]. Starting from time point O, cells can be divided into two branches based on the sign of the second PHATE coordinate (12). Quantitative analysis reveals pronounced branch imbalance: for instance, at time point 4, 88 cells belong to one branch while 213 belong to the other, which poses a great challenge to the models that do not consider unbalancedness, resulting in cross-branch inference. Throughout this study, assume that the temporal change cell counts within each branch is entirely governed by a growth function, without contributions from migration observational noise.
By setting 7 = 5 and e = 0.03 according to our selection scheme described in Appendix B. and illustrated in Fig, A-6 (b). Our model successfully avoids cross-branch reconstructions in trajectory modeling, as illustrated in Fig. A-7. This leads to generated samples that remain well-aligned with the underlying manifold. In terms of growth rate estimation, our model also accurately captures the rapid expansion observed the branch below between time point 3 and time point 4, where the number of observed cells increases dramatically from 25 to 213, indicating the highest growth rate in the system.
338
(a) Predicted dynamics
(b) Predicted growth rate
Figure A-7: Predicted dynamics and growth rates by VGFM on Dyngen data.
B.4 Gausssian 1000D Data
DeepRUOT [26] employs a high-dimensional Gaussian mixture distribution (100D) [64] to evaluate the scalability of models. In this work, we adopt an even more challenging setting by testing scalability at 1000 dimensions. Following the setup of [26]. for the initial distribution, we generated 400 samples from the Gaussian located lower the (æ1, #2) plane, and 100 samples from the Gaussian positioned higher. For the final distribution, we generated 1,000 samples from the upper Gaussian, and 200
22

---

samples each from the two lower Gaussians, and assume cells in the upper region exhibit proliferation without transport [26]. We observe that simulation-bais -based methods [26] encounter training instability in the 1000-dimensional setting and flow matching-hased methods fail to learn reliable velocity field due to the unbalancedness of data, as shown in Tab. (Some works incorporating unbalancedness [37, 42] may as well reconstruct correct velocity field but fail to construct growth rate). In contrast, By setting 7 = 5 and e = 0.03 according to our selection scheme described in Appendix B.1 and illustrated in Fig, A-6 (c), Our method remains effective even in high-dimensional regimes for both velocity and growth field learning.
6N
(a) Predicted dynamics
(b) Predicted growth rate
Figure A-8: Predicted dynamics and growth rates by VGFM on Gaussian 1000D data.
B.5
Additional Experiment: Mouse Hematopoiesis Data
We also validate algorithm on the mouse hematopoicsis data previously analyzed in [23, 26, 73]. This dataset leverages lineage tracing to track differentiation trajectories. After applying batch correction integrate data across multiple experiments, the cells were embedded into two- dimensional force-directed layout (SPRING plot). The resulting visualization reveals pronounced bifurcation, where early progenitor cells diverge into distinct differentiation lineages.
-05
Predicted dynamics
(b) Predicted growth rate
Figure A-9: Predicted dynamics and growth rates by VGFM mouse hematopoiesis data.
Our model effectively captures the underlying branching structure in the data, and the predicted growth rates align well with those reported in [26], as well as with established biological priors and are consistent with known biological lineage patterns. Following the evaluation protocol of [26]. we also compute the Wasserstein-1 distance as a quantitative metric. Under this evaluation, by setting - 20 and 0.005 according to our selection scheme described in Appendix B. we train the model with warm-up stage of 200 iterations. After that, the distribution fitting loss Lot introduced for an additional 100 training epochs. Our method achieves superior performance, demonstrating improved trajectory inference. In addition, evaluate the RME (Relative Mass Error)
23

---

metric introduced in this paper. The results indicate that our model significantly outperforms the other baselines in terms of mass-matching reconstruction accuracy, as shown in Tab. A-9 and Fig. A-10.
Table A-9: Wi and RME over all time on mouse hematopoiesis data. Part of the results were adopted from [26].
t
t/
Models
W1
RME
W1
RME
SF [41]
0.167
0.190
uAM[44]
0.745
0.777
UDSB [57]
0.388
0.159
0.128
0.249
TIGON [23]
0.314
0.124
0.342
0.177
DeepRUOT [26]
0.145
0.140
0.132
0.202
VGFM
0.115
0.043
0.094
0.019
5
UDSB
Tigon
Relative Mass
4
DeepRuot
VGFM
Oberserved Mass
3
2
1
2
Figure A-10: Comparison of predicted relative mass (UDSB, TIGON, DeepRUOT, VFGM) with observed values from mouse hematopoiesis data.
B.6 EB Data
Extensive experiments the Embryoid Body (EB) [65] dataset with varying PCA dimensions conducted Sect. 4. We preprocess this dataset first by normalizing EB (5D) and leaving EB (50D) unnormalized. The proposed VFGM demonstrates strong performance in both cellular dynamics reconstruction and growth prediction, as illustrated in Fig. A-11. Specifically, we set the parameters 0.01, == 5 and normalize the cost matrix to ensure numerical stability when computing the semi-relaxed optimal transport plans (Eq. (11)).
Big batch strategy for large-scale datasets. When the number of observed data points is large, computing the transport plan between two consecutive time points using Eq. (11) becomes compu- tationally inefficient. To address this, adopt a Big-Batch Strategy. Specifically, we partition the data at each time point into n subsets, referred to as Big Batches. For each of the n Big Batches, we precompute and store the transport plans between adjacent time points. During training, randomly select one Big Batch and sample smaller mini-batches from it based on the precomputed transport plan for model optimization. In real-world datasets, set 5, while in synthetic datasets above we do not apply this strategy, i.e., n = 1, as the sample size is relatively smaller. The mini-batch size set to 256 for all experiments, except for Dyngen, where its limited number of samples.
smaller batch size of 60 due to
Comparison with sinkhorn divergence as distribution fitting loss LoT. We compute the distri- bution fitting loss Wi Sect. 3.4 using EMD distance by pot library, following MIOFlow [31] and DeepRUOT [26]. Particularly, compute the OT using the pot library, using the function

---

pot. emd O) This function solves for и using a network flow algorithm, which is not compatible with PyTorch's i automatic differentiation. However, the gradients of the cost function c can still be backpropagated.
We replace the EMD distance with Sinkhorn divergence [70], which is fully differentiable and enjoys many favorable properties, using CUDA/C+ based geomloss library. Tab. A-10 verified the effectiveness of Sinkhorn divergence for the distribution fitting loss.
Sinkhorn divergence between distributions Pu and P1 is defined as S. (Po. W. Pi) (Po. Pı) 2 (P1, P1), where (Po, P1) is the cetropy-regulaizzer optimal transport distance (Sinkhorn distance) and >> 0 is the entropy regularization parameter.
Table A-10: Comparison EMD and sinkhorn divergence (with different entropic regularization parameters) on EB 50D dataset with metrics (W>//RME at each timepoint).
LoT
t
t
t3
t4
Training time(min)
EMD
7.951/0.033
8.747/0.042
9.244/0.019
9.620/0.044 044
13
( 0.001)
7.902/0.018
8.767/0.013
9.063/0.083
9.507/0.096
9
(. 0.002)
7.904/0.032
8.791/0.033
9.111/0.104
9.523/0. 120
9
(. 0.005)
7.917/0.048
8.768/0.062
9.086/0.137
9.518/0.151
9
(. = 0.01)
7.904/0.035
8.801/0.038
9.102/0.113
9.511/0.127
9
( 0.05)
7.908/0.036
8.787/0.045
9.080/0.120
9.512/0.129
9
(E 0.1)
7.916/0.030
8.773/0.022
9.111/0.091
9.591/0.108
9
It can be found that the sinkhorn divergence achieves further improvement on the accuracy of dynamics reconstruction, and costs less time for computing the distribution fitting loss.
031
057
(a) Predicted dynamics
Predicted growth rate
Figure A-11: Predicted dynamics and predicted growth rates by VGFM on 5D EB data with timepoint 1 held out.
We also examine whether the distribution fitting loss might lead to overfitting on the given population snapshots, In the hold-one-out experiments (Tab. 2), we relabel the timestamps of the four timepoints 1, 2, 3, and 4. Following the protocol in [38], we evaluate the model by holding intermediate timepoint at time (i.e., 1, 2, and 3) and report the average Wasserstein-1 distance. While hold-one-out experiments have already validated the superiority of VFGM over other methods, we further analyze the loss curves of both the distribution fitting loss and the W1 distance on hold-out distribution (e.g., the distribution The results show that the Wi distance on the hold-out distribution decreases during the optimization of LOT, indicating the absence of overfitting and demonstrating the generalization ability of the distribution fitting loss in modeling unseen distributions.
one
time 1) that is not seen during training (Fig. A-12).
25

---

10
0.66
0.5
07
a6-
20
Figure A-12: Loss curves of the training loss after incorporating Lor (blue) and the Wi distance on hold-out distribution at timepoint 1 (red) 5D EB data.
B.7 CITE Data
CITE-seq (Cellular Indexing of Transcriptomes and Epitopes by Sequencing) [66] is advanced technique that enables the simultancous profiling of transcriptomes and surface protein expression the single-cell level through the use of antibody-derived tags. In this study, we utilize only the gene expression matrix from the CITE-seq dataset, and preprocess the data by normalizing CITE (5D) and leaving CITE (50D) unnormalized. The experimental settings for the CITE dataset are consistent with those used in the EB data experiments, employing PCA with 5 and 50 dimensions, and hyperparameters set as e = 0.01, т = We relabel the timestamps of the four timepoints as 2, and 3. Following the protocol in [38], we evaluate the model by holding out one intermediate timepoint at time (i.e., I and 2) and report the average Wasserstein- I distance.
To assess the model's capability to capture state trajectories (Fig. A-13 (a)) and predict mass growth (Fig. A-13 (b)), the distribution at time point is held out. Additionally, the distribution time point 2 is held out to evaluate the performance of the Cor (Fig. A- 14). The strong performance observed in the visualizations substantiates the potential the proposed VGFM framework.
(a) Predicted dynamics
(b) Predicted growth rate
Figure A-13: Predicted dynamics and predicted growth rates by VGFM on 5D CITE Data with timepoint held out.
26

---

Predicted dynamies (w/o LOT)
Predicted dynamics
(a)
(b) Wi hold-out timepoint: 36.462
Wi hold-out timepoint: 38.130
Figure A-14: Visualization of predicted dynamics by (a) VFGM (wlo Lor) and (b) VGFM on CITE
(50D) dataset, where the hold-out time is the second intermediate timepoint.
B.8 Pancreas Data Analysis
To further explore the scalability of VGFM, we applied our model and compared methods explicitly modeling 9r (z) [23, 26] to the Pancreas dataset [67] comprising measurements of the developing mouse pancreas. We select cells at day 14.5 and 15.5 and relabel them as 0, and select 2000 highly variable genes. As shown in [37], we assume the following cell type transitions are exclusively correct (denoted by -+), i.e., there is descending cell type (or set of cell types) other than the given one. We partition all considered cell type transitions into three branches.
Endocrine branch (ED) transitions.
Fev+ Alpha (FA) Alpha (A)
Fev+ Beta (FB) - Beta (B)
Fev+ Delta (FD) - Delta (D)
Fev+ Epsilon (FE) -> Epsilon (E)
A A
B + B
D - D
E E
Ngn3 EP transitions.
Ngn3 high early (NE) ED
Ngn3 high late (NL) -> ED
Non-endocrine branch (NEB) transitions.
Ductal (DU) + DU
Tip (T) - Acinar (AC)
AC AC
We observed that VGFM the only method showing steadily decreasing training loss, both for LyGFM and Lor- We report W1 and RME as shown in Tab. A-11.
This indicates that both loss terms still effective in high-dimensional real-world data.
Analysis on mean and variance. We calculated the means and variances of the real and generated gene at day 15.5 and plotted the corresponding mean-variance trend (Fig. A-15 (a), (b)) and histograms (Fig. A-15 (c), (d)). The results show that the generated samples closely follow the vean-variance trend of the real data, especially by incorporating Lot.
27

---

Table A-11: W1 and RME at day 15.5 between real data and generated data from VGFM and its variant.
Method
W1
RME
VGFM (w/o Lor)
27.002
0.025
VGFM
24.416
0.017
Tirepsint
Timepoinc
Real
Real Generated
1.75
Generated
1.50
1.25
VogVar
TogVar
0.75
0.50
0.25
0.00
0.00
1.00
1.50
2.00
0.50
1.25
logeMean
logiMean
(a) Mcan-Variance Variance trend (w/o Lot)
(b) Mcan- Veran-Varaance trend
Time Mean
Veriance
3.0
Cenerated
Cenerated
2.5
2.0
2.0
1.5
1.5
1.0
0.5
0.5
(c) Mean- Variance histogram (w/o Lor)
Tima Mean
Time 2 Mariania
3.0
Real
Generated
2.0
2.0
2.5
1.0
1.0
0.0
0.0
3
5
(d) Mean- Veran-Variance histogram
Figure A-15: mean- -eaannvaiancee trend ((a), (b)) and histograms ((c), (d)) of real and generated gene data.
28

---

Interpretable learned growth function. We highlight that our main contribution lies in modeling and training the growth function (x(x),, with VGFM in the flow matching framework, leverage snapshot mass change for learning (2. f). During training. the growth loss rapidly converges to negligible value, much faster than the velocity loss. To further investigate Iu (x; t), We first calculate growth rate of each cell at the initial time point (day 14.5) and visualize them in Fig. A-16. At this time, the proliferation observed during this developmental stage mainly originates from Acinar and Ductal cells. Remarkably, our model successfully recapitulates this pattern without being given any prior knowledge of the cell types, demonstrating its ability to infer biologically meaningful dynamics directly from the data.
Normalized growth Acinar and Ductal cells are known to have high growth rate at this time.
day on
brans
NEB
Mgn EP
cell Eype
hormalized arowth
Acthar
Feert Delia
0.6
Mgn late
Ducta
Fore Alpha
Cpsiion
0.2
Fov
Fave Epci on
Mon3 High early
Cell Index
Figure A-16: Normalized growth rate of each cell time point 0 (day 14.5).
We also demonstrate that the learned 9w (z, not only successfully recapitulates the relative growth rates of different cell types without any prior information about cell identities, but also that the branch-wise mass obtained from numerical integration of the ODE closely matches the total branch mass at time point! 1.
Mass (relative to total mass at day 14.5) of 3 different branches (EB, Ngn3 EP, NEB).
1.0
0.96
Mass at day 14.5
0.91
Predicted mass at day 15.5
0.8
day 15.5
0.6
Mass
0.4
0.29
0.2
0.15
0.13
0.10
0.08
0.0
EB
Ngn3 EP Branch
NEB
Figure A-17: Relative mass of different branches
0g(2.()
we then follow [23], for each time compute the element-wise absolute value gene space for every cell, and average across cells to quantify each gene's contribution to growth
in the
Ox,
29

---

dynamics. From our analysis of (2, t), We have identified Pnliprp1, Clps, and Ctrb1, as shown in Fig. A-18. These genes are well-known markers of Acinar cells and are highly expressed in the exocrine pancreas. Notably, shown above, cells from the non-endocrine branch, particularly Acinar and Ductal cells, are much more abundant later stages due to their high proliferation rates. This alignment between the learned growth-driving genes and known biological processes indicates that VGFM captures growth dynamics in a biologically interpretable manner, successfully linking the latent growth function to meaningful cell-type-speeific proliferation patterns.
091
Codl2
Key genes at time
Key genes at time
Figure A-18: Key genes identified by In (sr, 1) time (a) and (b).
C
Broader Impacts
Our method provides a scalable and efficient framework for modeling cellular dynamics, enabling trajectory reconstruction and growth rate inference in nigh-dimensional single-cell datasets. In biomedical and clinical settings, such capabilities can facilitate a deeper understanding of devel- opmental processes, disease progression, and response to treatment a single-cell resolution. For example, modeling differentiation trajectories and proliferation patterns of stem immune cells could inform therapeutic strategies in cancer, regenerative medicine, and immunotherapy.
However, the proposed method inherently data-driven and relies on statistical patterns learned from observational data. As such, it may produce biologically implausible trajectories or growth behaviors that conflict with known biological priors, especially when the training data is noisy, biased, or incomplete. This could potentially lead misleading interpretations or incorrect clinical hypotheses if not carefully validated by domain experts. Therefore, we emphasize that any downstream medical conclusions drawn from the model's output should be interpreted with caution and in conjunction with biological prior knowledge and experimental validation.
30

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification:
The main claims made in the abstract and introduction accurately reflect the paper's contributions and scope.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
is fine to include aspirational goals as motivation as long as is clear that these goals not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations the work performed by the authors?
Answer: [Yes]
Justification: Please Sect.
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model sell-specification, asymptotic approximations only holding locally). The authors should reflect how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution is low images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't It acknowledged judgment and recognize that individual actions in favor of transparency play impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
the paper. The authors should use their best
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
31

---

Justification: We provide the full set of assumptions and complete (and correct) proof for each theoretical result.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that it affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: The paper fully discloses all the information needed to reproduce the main experimental results of the paper.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is novel architecture, describing the architecture fully might suffice, or if the contribution a specific model and empirical evaluation, it may be necessary either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of large language model), releasing of a model checkpoint, or other means that appropriate to the research performed.
While NeurlPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome deseribe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
32

---

Does the paper provide open access to the data and code, with sufficient instruc- tions faithfully reproduce the main experimental results, as described in supplemental material?
Question:
Answer: [Yes]
Justification:
Our code is available at
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips.. cc/
http://gguuieeeoooetttiissssssss soooiiiiissioiicoiiiooo for more details.
While we encourage the release of code and data, we understand that this might not be possible, "No" an acceptable answer. Papers cannot be rejected simply for including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https:
/ //wpp...c//////// httpccccppppccppceeeeee wwwlcc/ggiss///eee /CodeSummm cesionPoll onPolicy) for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing much information possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: The paper specifies all the training and test details necessary to understand the results.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to a level of detail that is necessary appreciate the results and make sense of them.
The full details can provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Question: Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Answer: [No]
Justification: The current version of the paper does not report enror bars or statistical significance metries.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall with given experimental conditions).
33

---

The method for calculating the error bars should be explained (closed form formula, call a library function, bootstrap, etc.)
The assumptions made should be given (e.g.. Normally distributed errors).
It should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report 1-sigma error bars, but one should state it. The authors should preferably report a 2-sigma enror bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not to show in tables or figures symmetric error bars that would yield results that are range (e.g. negative enror rates).
If error bars reported in tables plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification: The paper provide sufficient information on the computer resources needed to reproduce the experiments.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs as well as estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make it into the paper).
9. Code of ethics
Question: Does the research conducted the paper conform, in every respect, with the NeurlPS Code of Ethics https: https::/ httpa:////e//sp///eeee cc http:////////iiiiccctteeeeeeeeee
Answer:[Yoe
Justification: The research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurlPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: Please see Appendix C.
Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
34

---

Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), faireess considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Does the paper describe safeguards that have been put in place for responsible release of data models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification:
The paper poses such risks.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Are the creators or original owners of assets (e.g.. code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Question:
Answer: [Yes]
Justification: The creators original owners of assets used in the paper are properly credited and the license and terms of use are explicitly mentioned and properly respected.
Guidelines:
The answer NA means that the paper does use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset is used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
35

---

If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. www.coo//a/oo/oo datasatssss has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information not available online, the authors are encouraged to reach out to the asset's s creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification: The paper does not release new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/moode as part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsoureing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well details about compensation (if any)?
Answer: [NA]
Justification:
The paper does not involve crowdsourcing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country institution) were obtained?
Answer: [NA]
Justification: The paper does not involve crowdsoureing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
36

---

We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Question:
Answer: [NA]
Justification: The core method development in this research does involve LLMs as any important, original, or non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to our LLM policy htops://nnerips... http://p//../oooreeeee////////// for what should or should not be described.
37


Reviews:
[{"review_id": "q5qthRrORy", "reviewer": "Reviewer_kHRn", "summary": "The paper presents a new approach, called Velocity-Growth Flow Matching (VGFM), for learning single-cell population dynamics from snapshot data. Unlike existing methods that either rely on simulation-heavy training or ignore the mass imbalance due to cell proliferation and death, VGFM jointly models cell-state transitions (velocity) and mass variations (growth). The method leverages a theoretical foundation based on a \"two-period\" interpretation of semi-relaxed optimal transport, separating the processes of mass growth and cell-state transitions to derive an ideal dynamics formulation. Neural velocity and mass-variation functions are then trained to approximate these dynamics via flow matching, augmented with a Wasserstein-based distribution fitting loss. Extensive experiments on synthetic datasets (e.g., Gaussian mixtures and Dyngen) and real-world single-cell datasets demonstrate that VGFM reconstructs trajectories and growth dynamics, achieving superior performance compared to multiple strong baselines.", "questions": "- Could the authors provide a discussion regarding the influence of batch size on their method? As both the growth matching loss and the OT fitting loss rely explicitly on mini-batch optimal transport computations, this aspect is crucial. Statistically, we should expect that increasing the batch size would improve performance, as larger batches yield more accurate OT approximations and stable gradient estimates. Providing empirical evidence or a sensitivity analysis on how batch size affects model performance and stability would significantly strengthen the evaluation.\n\n- Could the authors further elaborate on their motivation for choosing a time-independent growth function? Currently, the manuscript briefly states (lines 177–178): \"For simplicity and ease of implementation, we choose a time-independent form...\". However, it would be beneficial to explicitly discuss why this simplification—assuming that mass variations depend solely on spatial variables—is biologically or computationally justified. Clarifying whether this choice could potentially limit model expressiveness or generalization would improve the manuscript’s transparency and reader understanding.\n\n- Could the authors clarify whether they rely on the Sinkhorn algorithm to compute the OT fitting loss (involving the Wasserstein-1 distance)? If so, do you choose a small entropic regularization parameter to closely approximate the true Wasserstein distance? Additionally, have you explored or considered using the Sinkhorn divergence (as introduced by  [Feydy et al., 2019]) instead of the standard entropic Wasserstein distance, given its favorable metric properties?", "limitations": "The authors clearly acknowledge important limitations of their work (section 5, page 9).", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "**Strengths.**\n\nThis paper has several substantial strengths. First, it addresses an important problem—learning single-cell dynamics from unpaired snapshot data—while explicitly handling the biologically relevant challenge of mass imbalance. The proposed theoretical framing via semi-relaxed optimal transport is both elegant and well-motivated. Experimentally, the authors conduct thorough evaluations using synthetic and real-world datasets, demonstrating notable performance improvements over existing methods, especially in high-dimensional settings (up to 1,000 dimensions), where competing approaches significantly struggle. Additionally, the empirical validation provides compelling visualizations of the learned trajectories and growth rates.\n\n**Weaknesses.**\n\n- As acknowledged by the authors (lines 37–41), incorporating the OT-based distribution fitting loss renders the training procedure simulation-based. Specifically, this requires both simulating the underlying ODE dynamics and backpropagating through these simulations, making the method computationally cumbersome. This design choice diverges from recent advances in generative modeling, which explicitly aim for simulation-free training strategies to improve scalability and efficiency.\n\n- Moreover, the authors do not explicitly describe how they compute the distribution fitting loss. Presumably, they rely on the Sinkhorn algorithm with Euclidean cost computed over mini-batches. If this is indeed the case, it should be clearly stated. This approach presents scalability concerns, as mini-batch OT approximations suffer substantially from the curse of dimensionality. A practical suggestion is to replace the Wasserstein-1 distance with a Sinkhorn divergence [Feydy et al., 2019], which possesses more favorable statistical properties. Nonetheless, even this improvement may not fully resolve scalability issues, particularly when dimensionality exceeds $d > 1,00$.\n\n- The authors do not discuss the critical issue of batch size selection, despite the fact that two of the three loss components (growth function matching and OT distribution fitting) depend explicitly on mini-batch OT computations. Statistically, performance should improve as batch size increases. The current choice of a batch size of 256 for most experiments seems relatively small. Exploring batch sizes of 1,024–2,048 is recommended, especially given that Sinkhorn-type algorithms (including their unbalanced variants) remain computationally feasible at these scales. Furthermore, the authors propose partitioning computations into subsets (Appendix B.6) for larger datasets, which could negatively impact the quality and stability of the learned velocity and growth functions.\n\n- Additionally, concerns arise regarding the stability and scalability of the proposed growth matching loss. While the theoretical justification for its introduction is appealing, the practical implementation—which involves regressing against the right marginal of empirical unbalanced OT couplings computed via mini-batches—is potentially unstable, particularly for smaller batch sizes. Specifically, estimating the growth function effectively involves approximating potentials from the dual unbalanced OT problem [Séjourné et al., 2024], a procedure known to suffer considerably from the curse of dimensionality. Moreover, this instability is exacerbated by the potential inconsistency in how mini-batches handle outliers; some points might have their mass reduced or discarded in one mini-batch, yet retained or increased in another.\n\n- Finally, the notation throughout the paper, particularly in Sections 3.3 and 3.4, is often cumbersome and difficult to follow. Simplifying notation in these sections would significantly improve readability and accessibility.", "was_revised": true, "final_justification": "I appreciate the thoughtful and detailed responses provided by the authors. However, I remain not entirely convinced by the current methodological approach. I believe the paper would significantly benefit from additional revisions and resubmission in the future.", "rebuttal": "We sincerely appreciate your valuable and insightful suggestions. \n### **Q1: Choice of distribution fitting loss diverges from simulation-free training.**\n\n**Necessity of distribution fitting loss:** Training the flow matching model only at the level of velocity and growth regression without an explicit global distribution-level fitting constraint may lead to inaccuracy distribution learning. Specifically, since the theoretically optimal velocity $v\\_\\theta(x,t)=\\mathbb{E}\\_{z|(x,t)}v_\\theta(x,t|z)$ in flow matching, there may exist a possible velocity average due to the expectation, and the expected velocity may deviate from the real velocity. Introducing the distribution fitting loss will enforce the global distributional match constraint over the evolved population of cells.\n\n\n**Computational Efficiency Comparison:** As the ablation study (Table 3 in Sect. 4.1, Page 9) shows, the computation cost needed for VGFM (13 mins) is **much less** than the totally simulation-based method DeepRUOT (90 mins), although VGFM sacrifices some computation efficiency compared with the totally simulation-free VGFM (w/o $\\mathcal{L}\\_{\\rm OT}$). At the same time, VGFM achieves the **best results** on the dynamics reconstruction among the totally simulation-based and totally simulation-free methods.\n\n### **Q2: Influence of batch size.**\n\nWe clarify that **our VGFM approach did NOT use mini-batch OT** either in the computation of the OT plan or the distribution fitting loss. Instead, we propose **an offline sampling strategy for OT plan computation** and its improved version tailored for large-scale datasets.\n\n**Offline sampling strategy**: **The OT plan between the adjacent snapshots for the whole population is precomputed and cached before training**, as **Algorithm 1 (Page 7)** shows. Then, in training of the flow matching, we sample a batch of paired timestep samples using this precomputed OT plan, alleviating the sampling bias. This strategy needs the computation of OT for the whole dataset. We further propose the **Big Batch strategy** in Page 22, Line 668. Specifically, we  **uniformly partition** the data at each time point into $n$ subsets as Big Batches. Then we precompute and cache the $n$ OT plans for these Big Batches and randomly select one Big Batch for constructing batches in the training process. **During training the velocity field and growth function, we don't have to compute any ot plan**, one just need to sample batch $(x_0,x_1)$ from the precomputed ot plan(s) $\\pi$ (to train $v$) and sums up the corresponding rows of $x_0$ (to train $g$) **(See Eq.(12), Page5)**. Since the pre-computed plan(s) is on full sample or Big Batches and is fixed during training, we have better training stability compared with other mini-batch ot approaches.\n\nIn the following Tables r1-1,r1-2,r1-3,r1-4, we report the sensitivity analysis on $n$ (Table r1-1: number of big batches in precomputing the OT over dataset), batchsize (Table r1-2: batchsize in training loss $\\mathcal{L}\\_{\\mathrm{VGFM}}$ ), sample size (Table r1-3: sampling size in distribution fitting loss $\\mathcal{L}\\_{\\mathrm{OT}}$ ), $n$ (Table r1-4: number of batches in computing distribution fitting loss $\\mathcal{L}\\_{\\mathrm{OT}}$ ), on EB 50D dataset under the same experiment setting of ablation study (Table 3, Page 9). \n\nTable r1-1: Sensitivity analysis of $n$ on EB dataset with metrics (computation time and $W_1$/RME at each time).\n\n|$n$|Time for precomputing OT plan(s)|t1|t2|t3|t4|Average|\n|-|-|-|-|-|-|-|\n|1|10.67|7.951/0.039|8.747/0.042|9.244/0.019|9.620/0.044|8.890/0.036|\n|2|2.96|7.985/0.035|8.823/0.057|9.284/0.035|9.649/0.071|8.935/0.049|\n|5|2.54|7.997/0.038|8.969/0.031|9.436/0.009|9.799/0.042|9.050/0.030|\n|8|2.13|7.970/0.026|8.855/0.070|9.237/0.046|9.536/0.081| 8.899/0.055|\n\nTable r1-2: Sensitivity analysis of batch size in $\\mathcal{L}\\_{\\mathrm{VGFM}}$  on EB dataset with metrics ($W_1$/RME at each time t).\n\n|Batchsize in matching loss|t1|t2|t3|t4|Average|\n|-|-|-|-|-|-|\n|64|7.949/0.034|8.863/0.035|9.250/0.012|9.712/0.047|8.943/0.032|\n|128|7.995/0.023|8.942/0.060|9.227/0.036|9.685/0.069|8.962/0.047|\n|256|7.951/0.039|8.747/0.042|9.244/0.019|9.620/0.044|8.890/0.036|\n|512|7.997/0.037|8.743/0.034|9.150/0.013|9.685/0.045|8.893/0.032|\n|1024|7.989/0.035|8.819/0.059|9.296/0.023|9.738/0.058|8.960/0.043|\n|2048|7.904/0.029|8.800/0.069|9.277/0.041|9.580/0.076|8.890/0.053|\n\nTable r1-3: Comparison of full sampling and mini-batch sampling in $\\mathcal{L}\\_{\\mathrm{OT}}$ on EB dataset with metrics ($W_1$/RME at each time t).\n\n|Sample size|t1|t2|t3|t4|Average|\n|--|--|--|--|--|--|\n|64|8.251/0.032|9.221/0.046|9.745/0.028|10.366/0.057|9.395/0.040|\n|128|8.079/0.030|9.095/0.053|9.625/0.029|9.930/0.057|9.182/0.042|\n|256|8.117/0.032|9.125/0.052|9.758/0.029|9.892/0.057|9.223/0.042|\n|512|8.021/0.034|8.886/0.046|9.184/0.022|9.632/0.048|8.930/0.037|\n|1024|7.978/0.034|8.772/0.044|9.214/0.017|9.519/0.042|8.870/0.034|\n|2048|7.939/0.037|8.807/0.039|9.193/0.016|9.499/0.041|8.859/0.033|\n|full samplling|7.951/0.039|8.747/0.042|9.244/0.019|9.620/0.044|8.890/0.036|\n\n\nTable r1-4: Results of using Big Batch strategy in $\\mathcal{L}\\_{\\mathrm{OT}}$ to speed up computation on EB dataset with metrics ($W_1$/RME at each time t and the training time).\n\n|$n$|t1|t2|t3|t4|Average|Training time(min)|\n|-|-|-|-|-|-|-|\n|1|7.951/0.039|8.747/0.042|9.244/0.019|9.620/0.044|8.890/0.036|13|\n|2|7.938/0.035|8.806/0.040|9.134/0.014|9.545/0.043|8.855/0.033|11|\n|5|7.948/0.034|8.764/0.041|9.153/0.016|9.435/0.045|8.825/0.034|10|\n|8|7.968/0.033|8.774/0.043|9.172/0.017|9.446/0.046|8.840/0.034|10|\n\nThese results demonstrate the robustness of performance to the number of big batches in the offline sampling strategy (Table r1-1) and distribution fitting loss $\\mathcal{L}\\_{\\mathrm{OT}}$ (Table r1-4), and the effect of batch size in $\\mathcal{L}\\_{\\mathrm{VGFM}}$ and sampling size in computing distribution-fitting loss $\\mathcal{L}\\_{\\mathrm{OT}}$ (Tables r1-2/3) on the final performance.\n\n### **Q3: Clarification on if relying on the Sinkhorn algorithm and exploration of Sinkhorn divergence.**\n\nWe compute the distribution fitting loss $W_1$ in Section 3.4 using EMD distance by `pot` library, following MIOFlow [Guillaume et al., NeurIPS22] and DeepRUOT [Zhang et al., ICLR25]. Particularly, we compute the OT using the `pot` libary, using the function `pot.emd()` and **only backpropagate the gradient with respect to the cost $c$** in $W_1$ since $\\pi$ is solved using a network flow algorithm, which is not compatible with PyTorch’s automatic differentiation. We will make a clear statement in Appendix.\n\n\nAs suggested, we replace the EMD distance with Sinkhorn divergence using CUDA/C++ based `geomloss` library. Table r1-5 verified the effectiveness of Sinkhorn divergence for the distribution fitting loss. \n\nTable r1-5: comparison of EMD and sinkhorn divergence (with different entropic regularization parameters) on EB dataset with metrics ($W_1$/RME at each time step).\n|Distribution fitting loss|t1|t2|t3|t4|Average|Training time(min)|\n|-|-|-|-|-|-|-|\n|EMD|7.951/0.039|8.747/0.042|9.244/0.019|9.620/0.044|8.890/0.036|13|\n|$S_{\\varepsilon}(\\varepsilon=0.001)$|7.902/0.018|8.767/0.013|9.063/0.083|9.507/0.096|8.809/0.052|9|\n|$S_{\\varepsilon}(\\varepsilon=0.002)$|7.904/0.032|8.791/0.033|9.111/0.104|9.523/0.120|8.832/0.072|9|\n|$S_{\\varepsilon}(\\varepsilon=0.005)$|7.917/0.048|8.768/0.062|9.086/0.137|9.518/0.151|8.817/0.099|9|\n|$S_{\\varepsilon}(\\varepsilon=0.01)$|7.904/0.035|8.801/0.038|9.102/0.113|9.511/0.127|8.829/0.078|9|\n|$S_{\\varepsilon}(\\varepsilon=0.05)$|7.908/0.036|8.787/0.045|9.080/0.120|9.512/0.129|8.821/0.082|9|\n|$S_{\\varepsilon}(\\varepsilon=0.1)$|7.916/0.030|8.773/0.022|9.111/0.091|9.591/0.108|8.847/0.062|9|\n\nIt can be found that the sinkhorn divergence achieves further improvement on the accuracy of dynamics reconstruction, and costs less time for computing the distribution fitting loss. We will include this result of sinkhorn divergence in the appendix and cite it in the discussion section of the manuscript.\n\n### **Q4: Motivation for time-independent growth function.**\n\nThe selection of time-independent growth function is because of\n+ **First, it minimizes $L^2$ energy potential.** As we have discussed in Appendix A.4, time-independent growth function is not only easy to implement but also minimizes the $L^2$ energy potential, which satisfies the principle of least action.\n+ **Second, it can be explained by the Malthusian growth model**, i.e., the exponential growth model. The Malthusian growth model is $\\frac{\\mathrm{d}p_t}{\\mathrm{d}t}=gp_t$, where $p_t$ is the population size and $g$ is a constant growth rate. Our model on the growth $\\frac{\\partial p_t}{\\partial t}=g_tp_t$ (refer to Eq.5, Page 4, Line 145-146) is consistent with the above Malthusian growth model. **In these models, g is treated as a constant** based on the assumption that the resources are abundant and the environment is stable, causing growth rates to remain relatively stable over time. In the context of scRNA-seq experiments, these conditions are **typically satisfied** because cells are often cultured or sampled under controlled laboratory conditions, where nutrient supply, temperature, and other environmental factors are maintained at constant levels. Thus, choosing this form is biologically reasonable, especially when there is no prior knowledge, such as the maximum carrying capacity of the environment.\n\n### **Q5: On notations.**\n\nThanks for your suggestion. Upon revisiting our paper, we realized that the discrete entropic semi-relaxed optimal transport formulation in Eq. (11) can be removed, as it is a standard result. Accordingly, we can directly use $\\pi^*$ in place of $\\pi^{0\\to1}$. In addition, we can describe the barycentric mapping and $N_i$ in plain language (line 188) for better clarity. We will also carefully review the manuscript to refine the notations and formulations throughout."}, {"review_id": "37Ad713orf", "reviewer": "Reviewer_fzn6", "summary": "The paper proposes VGFM (Joint Velocity-Growth Flow Matching), a scalable, simulation-free approach to modeling single-cell dynamics from unpaired, unbalanced snapshot data. By framing the dynamics as a two-period semi-relaxed optimal transport problem, the authors decouple and then recombine velocity and growth fields into a joint dynamic system. The method uses neural networks to parameterize these fields and optimizes a flow-matching loss with an additional distribution-fitting term. Experiments demonstrate VGFM’s performance on synthetic Gaussian data, a 5D EB dataset, and CITE-seq data.", "questions": "- The scalability of VGFM is very interesting. Would training directly on gene-level counts (say 1k-2k highly variable genes) be feasible? Do the authors have an intuition how the sparsity of scRNA-seq would be handled by VGFM? \n\n- Could you elaborate on the motivation behind the distribution-fitting loss? I’m curious whether it was derived from a principled rationale or chosen empirically.\n\n- I would be interested to hear the authors’ perspective on the biological significance of explicitly modeling growth effects. Methodologically, this choice is well motivated and aligns well with the nature of the data. \n\n- Could the authors elaborate on how their model could be made conditional on covariates (e.g., cell type, batch)? Incorporating conditioning is important in single-cell modeling and would enable meaningful applications such as data augmentation and batch correction.\n\n&nbsp; \n_____\nI appreciate the authors’ effort in addressing an important research question and look forward to their clarifications and responses during the rebuttal phase. I am happy to raise my score if the authors adequately address my questions and concerns.", "limitations": "While the authors acknowledge some limitations of their approach, the discussion remains focused on technical aspects. It would strengthen the manuscript to also address broader challenges in single-cell modeling (noise structure, sparsity, overdispersion, data acquisition). A discussion of where VGFM is likely to succeed and where it may face limitations in the context of real single-cell data would provide valuable perspective (for example, the scalability claim is only demonstrated on synthetic Gaussian data rather than realistic gene-level scRNA-seq datasets).", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "### Strengths \n- [Originality] The formulation of the dynamics via semi-relaxed OT with separate velocity and growth fields is novel and well-grounded. The manuscript clearly explains the derivation of VGFM and its training objectives.\n- [Significance] The proposed method allows training on very high-dimensional synthetic data (up to 1000D) efficiently, outperforming prior work in terms of scalability.\n\n### Weaknesses \n- [Significance] My biggest concern is the misalignment of motivation and experiments. Although the method is explicitly motivated as being tailored for single-cell dynamics, the experiments largely fail to demonstrate its merits in this domain. The synthetic Gaussian and low-dimensional EB experiments are insufficient to validate biological relevance or single-cell-specific challenges. For this, it is essential to analyse the modeled samples (i.e. gene) space. I would suggest to investigate VGFM in the following contexts: \n     - A comparison of the gene-wise empirical mean-variance trend with that of real data and an analysis of the modeled sparsity in gene space, cf. 5.1 in [1]\n     - The inclusion of the MEF dataset [2] as well as the analysis of the decoded space as done in 5.2 of [3]\n     - An analysis on the Pancreas dataset [4], checking that with VGFM one can identify all terminal states as well as investigating the dynamics of linage drivers. For the linage drivers, in particular, it would be interesting if the growth modelling enables more accurate modeling of the real data, cf. [3].\n- [Quality] I encourage the authors are encouraged to discuss key characteristics of scRNA-seq data, such as sparsity, overdispersion, and its negative binomial noise structure. Since VGFM can readily be combined with alternatives to PCA that better reflect the data’s nature, I encourage the authors to consider and elaborate on this connection (incorporating an NB noise model may also be relevant for the additional experiments suggested above).\n- [Clarity] Key dataset details (e.g., number of cells, timepoints, gene dimensions; see L233–234) are missing. Including this information would take little space but greatly improve the clarity of the experimental setup for the reader.\n- [Minor] The additional distribution loss term (L68–69) appears a bit ad-hoc and lacks justification (beyond empirical performance) compared to the VGFM part.\n\n&nbsp; \n\n_____\n\n- [1] Palma et al. \"Multi-Modal and Multi-Attribute Generation of Single Cells with CFGen\", ICLR 2025\n- [2] Schiebinger et al. \"Optimal-transport analysis of single-cell gene expression identifies developmental trajectories in reprogramming\", Cell, 176(4):928–943, 2019.\n- [3] Palma et al. \"Enforcing Latent Euclidean Geometry in Single-Cell VAEs for Manifold Interpolation\", ICML 2025.\n- [4] Bastidas-Ponce et al. \"Comprehensive single cell mRNA profiling reveals a detailed roadmap for pancreatic endocrinogenesis\" Development, 146(12):dev173849, 2019.", "was_revised": true, "final_justification": "The authors addressed most of my concerns and put great effort into the rebuttal. I believe that the additional experiments really benefit the paper and demonstrate the potential of the VGFM method.", "rebuttal": "We sincerely appreciate your valuable and insightful suggestions. We have conducted additional experiments and will **include results and discussions in our manuscript and appendix**.\n\n### **Q1: The scalability of VGFM on real-world high-dimensional gene-space datasets.**\n\nTo further explore the scalability of VGFM, we applied our model and compared methods explicitly modeling $g\\_t(x)$ (Sha *et al*., 2024; Zhang *et al*., 2025) to the Pancreas dataset with 1000 highly variable genes.\n\n+ **Stable training loss curve.** We observed that VGFM is the only method showing a steadily decreasing training loss, both for $\\mathcal{L}\\_{\\mathrm{VGFM}}$ and $\\mathcal{L}\\_{\\mathrm{OT}}$. We compute $W\\_1$ at day 15.5 between real data and generated data from VGFM and its variant VGFM (w/o) $\\mathcal{L}\\_{\\mathrm{OT}}$, as shown in the table below:\n\nTable r2-1: $W_1$ at day 15.5 between real data and generated data from VGFM and its variant.\n\n|Method|$W_1$|\n|-|-|\n|VGFM (w/o $\\mathcal{L}_{\\mathrm{OT}})$|24.583|\n|VGFM|21.026|\n\nThis indicates that both loss terms are still effective in high-dimensional real-world data.\n\n+ **Analysis on mean and variance.** We calculated the means and variances of the real and generated data at day 15.5 and plotted the corresponding trend. The results show that the generated samples closely follow the mean–variance trend of the real data. To present our results, we provide a table showing binned statistics, summarizing the distribution of real and generated data across different mean and variance intervals.\n\n\nTable r2-2: Distribution of real and generated data across different mean and variance intervals.\n|Mean intervals|Generated data count|Real data count|Variance intervals|Generated data count|Real data count|\n|-|-|-|-|-|-|\n|[0.03,0.59)|512|514|[0.02,0.52)|829|765|\n|[0.59,1.15)|278|269|[0.52,1.01)|119|164|\n|[1.15,1.71)|87|88|[1.01,1.50)|25|32|\n|[1.71,2.28)|48|48|[1.50,1.99)|15|17|\n|[2.28,2.84)|46|50|[1.99,2.48)|5|9|\n|[2.84,3.40)|21|21|[2.48,2.97)|5|6|\n|[3.40,3.96)|2|4|[2.97,3.46)|1|2|\n|[3.96,4.53)|1|0|[3.46,3.96)|0|3|\n|[4.53,5.09)|4|5|[3.96,4.45)|0|1|\n|[5.09,5.65)|1|1|[4.45,4.94)|1|1|\n\nAs shown in the table, the distribution of means in the generated data exhibits a high degree of consistency with that of the real data, while the variances are generally slightly smaller than those of the real data but still remain reasonably accurate.\n\n+ **Analysis of the sparsity.**  Our current paper is applied to PCA-projected single-cell data, with the main contributions lying in the theoretical formulation and algorithmic design for the **joint dynamic modeling of $g_t$ and $v_t$**, from which we derive our flow-matching-based approach. In its present form, VGFM does not explicitly impose sparsity constraints on the data, generating data not as sparse as the real data. Nevertheless, recent studies have introduced encoder–decoder frameworks that operate in a latent space while explicitly modeling single-cell gene expression properties such as sparsity and the negative binomial distribution (Palma *et al.*, 2025). By integrating VGFM into such a framework, for example, using advanced encoders to map raw count data into a latent space and appropriately designing decoders to reconstruct gene counts, one could learn the growth function from OT plans in the latent space with more accurate metrics beyond the Euclidean cost, while simultaneously preserving biologically relevant statistical properties. We believe this **direct extension to raw data with explicit sparsity and negative binomial modeling** is a promising direction, deserving thorough investigation in future work. We will cite these suggested papers and discuss them in the manuscript.\n\n+ **Computing fate probabilities.** We follow Eyring *et al*. (2024) to compute the averaged fate probabilities of each ground truth transition. Notably, in almost each transition, VGFM ranked among the top three performers, despite operating directly in the 1,000-dimensional gene expression space, whereas all other methods were evaluated in a 10–50-dimensional PCA space, as done in Eyring *et al*. (2024). This might indicate that VGFM can achieve biological interpretability even in high-dimensional settings.\n\nTable r2-3:Averaged fate probabilities of each ground truth transition. \n|Model|dim|FA→A|A→A|FB→B|B→B|FD→D|D→D|FE→E|E→E|NL→ED|DU→DU|T→AC|AC→AC|\n|-|-|-|-|-|-|-|-|-|-|-|-|-|-|\n|TrajectoryNet|10 PCs|0.07|0.46|0.07|0.39|0.11|0.75(3)|0.10|0.52|0.01|0.32|0.79(3)|0.99(2)|\n|scVelo|50 PCs|0.79(1)|0.80(1)|0.30(1)|0.61|0.03|0.52|0.04|0.56(3)|1.00(1)|0.32|0.04|0.90|\n|WOT|50 PCs|0.37|0.62(2)|0.18(2)|0.44|0.19(3)|0.74|0.55(1)|0.49|0.50|0.82(2)|0.48|0.84|\n|OT-ICNN|50 PCs|0.51(2)|0.54|0.11|0.70(3)|0.39(2)|0.94(2)|0.48|0.62(2)|0.08|1.00(1)|1.00(1)|1.00(1)|\n|UOT-ICNN|50 PCs|0.25|0.60|0.29(2)|0.75(1)|0.78(1)|0.99(1)|0.36|0.55|0.73(2)|0.59|1.00(1)|1.00(1)|\n|VGFM|1000 HVGs|0.42(3)|0.80(1)|0.13|0.63(3)|0.14|0.70|0.64(1)|0.80(1)|0.73(2)|0.66(3)|0.39|0.91|\n\n+ **Interpretable learned growth function.** We highlight that our main contribution lies in modeling and training the growth function $g_t(x)$, with VGFM in the flow matching framework, to leverage snapshot mass change for learning $g_w(x,t)$. During training, the growth loss rapidly converges to a negligible value, much faster than the velocity loss. To further investigate $g_w(x,t)$, we follow Sha *et al.* (2024), for each time $t$, compute the element-wise absolute value $\\left|\\frac{\\partial g(x,t)}{\\partial x_i}\\right|$ in the gene space for every cell, and average across cells to quantify each gene’s contribution to growth dynamics. From our analysis of $g_w(x,t)$, the three genes with the largest absolute gradients are **Pnliprp1, Clps, and Ctrb1**. These genes are well-known markers of Acinar cells and are highly expressed in the exocrine pancreas. Notably, at this time point, cells from the non-endocrine branch, particularly Acinar cells, are much more abundant at later stages due to high proliferation rates. This alignment between the learned genes and known biological processes indicates that VGFM captures growth dynamics in a biologically interpretable manner, successfully linking the learned growth function to meaningful cell-type-specific proliferation patterns.\n\n### **Q2: The motivation behind the distribution-fitting loss.**\n\nOur ultimate goal is to generate data close to the real data distribution. Although flow matching is a powerful generative modeling framework, its learned velocity field is the expectation of the conditional velocity field, i.e.,$v\\_\\theta(x,t)=\\mathbb{E}\\_{z|(x,t)}v_\\theta(x,t|z)$， where $z=(x_0,x_1)$ is the entropic semi-relaxed ot plan in our setting, which is a **soft coupling**. For the same $x_0$, $x_1$ may come from different branches of the trajectory. This bias could accumulate during generation, possibly leading to deviation from the target distribution. To address this bias, we introduce a distribution-fitting loss that directly measures how well the generated samples match the empirical data distribution at each time point. As done in the ablation study (Page 9, Table 3) on loss terms in our paper, this complements the flow-matching objective and ensures that even if the velocity field is biased, the model still adjusts its dynamics to match the correct data distribution.\n\n### **Q3: The biological significance of explicitly modeling growth effects.**\nWe appreciate the reviewer’s positive feedback on our methodological choice of explicitly modeling growth effects. From a biological perspective, cell development is inherently a coupled process of **state transitions** and **mass changes**. Ignoring the growth component may lead to erroneous state transition inference, for example, reconstructing cross-branch trajectories that are biologically implausible, as we demonstrated in the simulation gene dataset experiment in our paper.  Moreover, modeling growth explicitly enables downstream biological interpretation. By computing element-wise absolute value $\\left|\\frac{\\partial g(x,t)}{\\partial x_i}\\right|$ in the gene space, one can **identify growth-related genes and proliferative cell populations**, as done in Sha *et al.*,2024. However, their approach relies on autoencoder-based dimensionality reduction. \nIn our response to Q1, we have analyzed our learned $g_w(x,t)$, which is trained directly on gene space, on the Pancreas dataset. We successfully identified key genes such as Pnliprp1, Clps, and Ctrb1, which are known to be associated with proliferating Acinar cells. This demonstrates that VGFM also yields biologically meaningful insights into cell population dynamics, even in high-dimensional settings where previous methods face limitations.\n\n### **Q4: How VGFM could be made conditional on covariates.**\nA natural approach is to encode cell type as a one-hot vector $c$ and integrate it into both $v_\\theta(x,t,c)$ and $g_w(x,t,c)$. The main challenge lies in determining the coupling between cells of the same type across two time points. This becomes considerably easier if prior knowledge of the transition fate for each cell type is available. While VGFM can identify proliferating cell types (as demonstrated in the Pancreas dataset), incorporating such prior knowledge could improve model accuracy. Moreover, VGFM could be extended to model perturbation responses conditioned on different drugs or gene perturbations. These embeddings could be generalized to unseen perturbations via the mode-of-action framework (Bunne *et al.*, 2022), enabling VGFM to predict which perturbations might facilitate or inhibit cell growth. \n\n### **Q5: Key dataset details.**\n\nThanks for your suggestions, we will add this table to our appendix.\n\nTable r2-4: Key dataset details.\n\n|Dataset|Timepoints|Cell counts|Dimensions|\n|-|-|-|-|\nSimulation Gene|5|3031|2|\nDyngen|5|728|5|\nGaussian|2|1900|1000\nMouse Hematopoiesis|3|10998|2\nEB|5|16819|5,50|\nCITE|4|31240|5,50|\nPancreas|2|18941|1000|"}, {"review_id": "GxlVfW0cp9", "reviewer": "Reviewer_DyDD", "summary": "This paper proposes a flow and growth model for modeling interpolation of unbalanced snapshots of an evolving distribution. The evolution is derived from a two-stage process of growth and then flow, with the assumption that $p_0$ and the marginal of the coupling for the flow part have the same support. The resulting growth factor and flow vector field is then mixed to occur simultaneously, and proven to result in $p_1$.\n\nImplementation of this framework is done via learning of the vector field and growth weight parameterized as neural networks. Training is done with sample access to the snapshot distributions and solving entropy-regularized semi-relaxed OT between the discrete distributions. An additional $W_1$ loss is also added to the result of the joint growth and flow model to improve alignment to the snapshots. Empirical performance is shown to outperform SOTA methods on synthetic and real datasets in achieving alignment to the snapshots.", "questions": "1. See point 1 in Weaknesses.\n2. Is there a problem with non-differentiability of $W_1$ in the pipeline?\n3. For Eq. 5, is there a reason you couldn't incorporate the energy from A.4 in order to get a unique specification for $g_t$? It would nice if it\nmirrored the uniqueness of the velocity field from the Benamou-Brenier formulation.\n4. Will an implementation be released?", "limitations": "Yes, though Weakness point 1 should be addressed in some way.", "rating": 5, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 4, "originality": 3, "strengths_and_weaknesses": "Strengths:\n1. Good theoretical justification of the approach is given.\n2. The method is relatively simple, but nontrivially combines existing tools for entropically-regularized transport and neural ODEs.\n3. Empirical performance and runtime seems to clearly beat existing methods.\n\nWeaknesses:\n1. It still seems unclear to me if the suggested model accurately reflects the potential behavior of genetic populations of cells. In particular, is it not possible for novel populations of cells to grow or appear outside the support of the current (or initial distribution) $p_0$? Or for whole populations to just disappear, and to have the remaining population be transported to the next snapshot? It would nice to comment on this or address this.\n2. As noted by the authors, the method is not fully simulation-free due to the additional $W_1$ loss against snapshots.", "was_revised": false, "final_justification": "", "rebuttal": "We sincerely thank the reviewer for the positive comments and insightful questions.\n\n### **Q1: It still seems unclear to me if the suggested model accurately reflects the potential behavior of genetic populations of cells. In particular, is it not possible for novel populations of cells to grow or appear outside the support of the current (or initial distribution) ? Or for whole populations to just disappear, and to have the remaining population be transported to the next snapshot?**\n\nCell differentiation naturally causes transitions in the state space, which can lead to changes in both the support set and the population size at the next time point. Our proposed VGFM aims to model such dynamics.\n\nIt is a particularly challenging scenario that part of the population in a future snapshot does not originate from the first snapshot, but gradually emerges at some intermediate time from outside the support of the earlier distribution. From our current biological understanding, such a “spontaneous appearance” of an entirely novel population without any lineage connection is highly unlikely, as new cell populations generally arise through differentiation from existing ones. That said, if such a situation were to occur, we acknowledge that VGFM in its current form could not handle it, as it assumes that the population at the next snapshot evolves entirely from the previous one. \n\nFor the second scenario mentioned, where a certain population disappears while the remaining population transitions to the next snapshot, VGFM can indeed identify such disappearing populations via the semi-relaxed OT plan, as these cells would have vanishing future weights. However, since such cells are not sampled during training, their growth rates and velocity fields could not be learned. Handling this case may require additional constraints, such as enforcing the velocity field to be zero and the growth rate to be very small in those regions, or incorporating biological priors to determine the likely fate and disappearance time of such populations.\n\n\n\n\n### **Q2: Is there a problem with non-differentiability in the pipeline?**\n\nWe appreciate the reviewer’s careful observation. We now clarify how the $W_1$ loss is computed. In practice, we use the `ot.emd()` function from the `pot` library to compute the transport plan $\\pi$ in $W_1$. This function solves for $\\pi$ using a network flow algorithm, which is not compatible with PyTorch’s automatic differentiation. However, the gradients of the cost function $c$ can still be backpropagated. This strategy is also adopted in the implementations of Zhang *et al.* (ICLR 2025) and Huget *et al.* (NeurIPS 2022). Therefore, the gradient of $W_1$ w.r.t. network parameters is partial rather than full. Nevertheless, this does not affect the effectiveness of our method nor the feasibility of computation, as demonstrated by our results (see sec. 4). Furthermore, in the rebuttal to reviewer kHRn (Please refer to Q3 of response to reviewer kHRn), we also explored replacing $W_1$ with the Sinkhorn divergence (Feydy *et al.*, AISTATS 2019), which is fully differentiable and enjoys many favorable properties. The Sinkhorn divergence and $W_1$ yield similar results. The corresponding results and discussion are provided in our response to Reviewer kHRn.\n\n### **Q3: For Eq. 5, is there a reason you couldn't incorporate the energy from A.4 in order to get a unique specification for $g_t$?**\n\nThe static semi-relaxed optimal transport in Eq. 4 is a **strictly convex** optimization problem: the first term is linear, the second term $\\mathrm{KL}(\\cdot \\| q)$ is convex, and the constraints are linear. Therefore, the optimal transport plan $\\pi^\\ast$ is **unique**, which in turn makes $p_\\lambda$ unique.  Given the additional constraint in A.4, we can then uniquely determine $g_t$. We aim develop an dynamic understanding in Eq. (5) of the static semi-relaxed OT in Eq. (4). Directly adding the $L^2$ energy term into the objective function of Eq. (5) would break its equivalence with Eq. (4). Instead, we treat it as a **constraint** to select $g_t$.\n\n\n### **Q4: Will an implementation be released?**\n\nYes. We are currently organizing and cleaning up the code, and it will be released publicly upon acceptance of the paper.\n\n### **Q5: Introducing distribution fitting loss needs simulations.**\n\nIndeed, our method is not fully simulation-free due to the additional distribution fitting loss, we stress the importance of distribution fitting loss:\n\n**Necessity of distribution fitting loss:** Training the flow matching model only at the level of velocity and growth regression without an explicit global distribution-level fitting constraint may lead to inaccuracy distribution learning. Specifically, since the theoretically optimal velocity $v\\_\\theta(x,t)=\\mathbb{E}\\_{z|(x,t)}v_\\theta(x,t|z)$ in flow matching, there may exist a possible velocity average due to the expectation, and the expected velocity may deviate from the real velocity. Introducing the distribution fitting loss will enforce the global distributional match constraint over the evolved population of cells.\n\n\n**Computational Efficiency Comparison:** As the ablation study (Table 3 in Sect. 4.1, Page 9) shows, the computation cost needed for VGFM (13 mins) is **much less** than the totally simulation-based method DeepRUOT (90 mins), although VGFM sacrifices some computation efficiency compared with the totally simulation-free VGFM (w/o $\\mathcal{L}\\_{\\rm OT}$). At the same time, VGFM achieves the **best results** on the dynamics reconstruction among the totally simulation-based and totally simulation-free methods."}, {"review_id": "uaUpDBVeAl", "reviewer": "Reviewer_G6tj", "summary": "This paper proposes VGFM (Velocity-Growth Flow Matching), a novel method for modeling single-cell dynamics from snapshot data that accounts for both state transitions and mass variation (e.g., due to cell proliferation or death). The authors introduce a dynamic reinterpretation of semi-relaxed optimal transport (SROT) as a two-phase process: first, mass growth occurs (modeled via a growth function), and then cells transition in state space (modeled via a velocity field). From this decomposition, the authors derive closed-form supervision targets for learning velocity and growth via neural networks using a flow matching approach. An additional distribution fitting loss is incorporated to improve the match between predicted and observed snapshots. The method outperforms existing state-of-the-art (SOTA) methods on both synthetic and real datasets, with superior accuracy, interpretability, and scalability.", "questions": "- Could the authors provide correlation analyses between predicted growth rates and known proliferation markers or synthetic ground truth?\n-  Have the authors considered estimating uncertainty in predicted velocity/growth (e.g., via dropout or ensembles)? \n-  Can the method be extended to scenarios where timepoint alignment is unknown or uncertain, possibly by integrating OT-based time alignment?\n- Effect of Sampling Bias: Since the method relies on observed cell counts for growth learning, have the authors tested robustness under non-uniform or subsampled cell distributions across time?", "limitations": "Yes\nAdditional suggestions:\n- Consider including diagnostics or visualizations of growth error spatially or per-cell type.\n- Investigate how growth modeling could be regularized or debiased using biological priors.", "rating": 6, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 4, "clarity": 4, "significance": 4, "originality": 4, "strengths_and_weaknesses": "Strengths:\n- The paper is technically sound and builds on a solid mathematical foundation (SROT and continuity equations).\n- The formulation of joint mass and state dynamics is well-posed, with provable equivalence to SROT objectives (Proposition 1, Theorem 1).\n- Empirical results are strong across diverse settings, including high-dimensional synthetic data and real scRNA-seq datasets.\n- Ablation studies and comparisons to unbalanced transport baselines (e.g., DeepRUOT, TIGON) are well done and validate all components.\n\nWeeknesses:\n- The estimation of the Monge map via barycentric projection introduces approximation error, especially when entropy regularization is non-negligible.\n- The growth function $g(x, t)$ is derived from empirical densities and may absorb technical artifacts (sampling bias, dropout), but this is not fully analyzed.\n- The method assumes aligned timepoint snapshots, limiting its direct application to asynchronous or pseudotime datasets.", "was_revised": true, "final_justification": "I appreciate the authors' effort to clarify their papers and address my concerns. As my concerns have been addressed, I would like to increase the rate by one score (from 5 to 6).", "rebuttal": "We thank the reviewer for the positive comments and insightful questions.\n\n### **Q1: Could the authors provide correlation analyses between predicted growth rates and known proliferation markers or synthetic ground truth?**\n\nAs suggested, we implement correlation analysis between the predicted growth rate and the ground truth to verify the accuracy of our growth rate modeling and learning as follows. We choose the Simulation Gene Data (Appendix B.2, Page 19) which is generated from the given dynamic systems (eq. (A-23), Page 19). The ground truth of growth rate can also be calculated analytically. Therefore, we can conduct **correlation analysis on the out-of-distribution (OOD) time points**, which are not included in the training data, to verify the generalization ability of VGFM. Results are reported in the following Table r4-1.\n\nTable r4-1: correlation analysis on four OOD time points on Simulation Gene Data\n\n\n| OOD time  | Pearson correlation coefficient |\n| :-------: | :-----------------------------: |\n| $t_{0.5}$ |              0.980              |\n| $t_{1.5}$ |              0.992              |\n| $t_{2.5}$ |              0.995              |\n| $t_{3.5}$ |              0.996              |\n\nThe growth rate predicted by our VGFM shows a high correlation with the ground truth, and we will add this analysis in the Appendix.\n\n### **Q2: Have the authors considered estimating uncertainty in predicted velocity/growth (e.g., via dropout or ensembles)?**\n\nWhile our current model is deterministic, dropout can be applied to our velocity/growth networks to obtain uncertainty.  During inference, we calculate the averaged variance of each dimension of velocity and growth over all the points across different branches (i.e., cell types) and times. We experiment on Simulation Gene Data (Appendix B.2), set dropout rate as 0.1 and execute 5 runs on the velocity/growth network during inference for all the points. Results are shown in the following Table r4-2. \n\nTable r4-2: uncertainty estimation of predicted velocity (2-dimensional) and growth on Simulation Gene Data. var_x1/var_x2  and var_g are the average variance of $1$-st/$2$-nd dimension of velocity and growth respectively.\n\n| branch | time | var_x1 (1e-3) | var_x2 (1e-3) | var_g (1e-3) |\n| ------ | ---- | ------------- | ------------- | ------------ |\n| 0      | 0    | 0.1080        | 0.1000        | 0.0550       |\n| 0      | 1    | 0.0770        | 0.0600        | 0.0660       |\n| 0      | 2    | 0.0660        | 0.0670        | 0.1100       |\n| 0      | 3    | 0.0850        | 0.0820        | 0.1910       |\n| 0      | 4    | 0.1220        | 0.1370        | 0.3500       |\n| 1      | 0    | 2.9320        | 6.5010        | 0.2510       |\n| 1      | 1    | 3.0230        | 4.2870        | 0.3060       |\n| 1      | 2    | 2.0410        | 1.9780        | 0.5590       |\n| 1      | 3    | 0.9530        | 0.9760        | 0.7660       |\n| 1      | 4    | 0.2800        | 0.3100        | 1.1680       |\n\n\nWe analyze the results at both spatial and temporal levels.\n\n**Spatial level:** In the Simulation Gene Data, cells of branch 0 relatively stand still with little variation on state and counts. The result shows that the variance of $v_1$, $v_2$ and $g$ in this branch is small, indicating low uncertainty, while in branch 1 the variances is larger, indicating relatively higher uncertainty.\n\n**Temporal level:** In the Simulation Gene Data, cells undergo significant state transitions and quantity changes in the initial phase, and the magnitude of the change decreases thereafter. Accordingly, the variance of $v_1,v_2$ and $g$ also reflects this trend, confirming our uncertainty estimation.\n\nWe will include these results as well as visulize the heatmap of uncertainty in the state space at each time points  in appendix.\n\n\n\n\n### **Q3: Can the method be extended to scenarios where timepoint alignment is unknown or uncertain, possibly by integrating OT-based time alignment?**\n\nOur current method assumes known time labels. In the absence of such labels, one possible extension is to first infer pseudotime using standard approaches (e.g., Monocle, Slingshot), then apply OT-based alignment between discretized pseudotime segments to estimate cell transitions and population growth. Alternatively, traditional pseudotime methods typically rely on low-dimensional embeddings (e.g., PCA, UMAP) followed by trajectory reconstruction. Inspired by this, we are also interested in future work that jointly learns cell evolution trajectories, velocity fields, and growth dynamics directly in the latent space, integrating pseudotime inference and dynamic modeling into a unified OT-based framework. \n\n### **Q4: Effect of Sampling Bias: Since the method relies on observed cell counts for growth learning, have the authors tested robustness under non-uniform or subsampled cell distributions across time?**\n\nIn the cell dynamics reconstruction, we usually assume that the observed dataset accurately reflects the evolution of true cell distribution across time. So the observed cell counts is consistent with the fact of cell proliferation and death. It's also practically meaningful to evaluate the robustness of VGFM trained on the biased observed dataset. \n\nFor the experiment settings, we use Simulation Gene Data (Appendix B.2) which has two spacial branches. To generate the non-uniform or subsampled data, we resample the original data with perturbed branch ratio. Specifically, we randomly select perturbation factor where $a_t\\sim\\mathcal{N}(0,\\sigma^2)$ for time $t$. Then if the original  branch ratio is $s_t:1-s_t$, the resampled branch ratios become $s_t-|a_t|:1-s_t+|a_t|$ . $\\sigma$ indicates the bias level. We train the model on the resampled data and evaluate on the original data, to testify the robustness of VGFM.\n\nTable r4-3: results of training on sampling bias on Simulation Gene Data with merics ($W_1$/RME) averaged over three random seeds.\n\n| bias factor $\\sigma$ |     t1      |     t2      |     t3      |     t4      |\n| :------------------: | :---------: | :---------: | :---------: | :---------: |\n|     0 (unbiased)     | 0.041/0.007 | 0.053/0.005 | 0.038/0.007 | 0.040/0.008 |\n| 0.1  | 0.074/0.035 | 0.094/0.051 | 0.075/0.063 | 0.077/0.032 |\n| 0.2  | 0.085/0.059 | 0.113/0.073 | 0.134/0.062 | 0.123/0.125 |\n\nWe also experiment on the real-world dataset EB (Appendix B.6). Since this dataset doesn't have obvious spatial branches, we take a different dropout method to construct the biased dataset and likewise **control the bias level with a factor $\\eta$**. Specifically, we sample several $a_t\\sim \\mathcal{U}([0,\\eta])$ with different bias factors, where $\\mathcal{U}([0,\\eta])$ is the uniform distribution on $[0,\\eta]$. Then, the points at time $t$ with a proportion of $a_t$ are randomly removed from the dataset. In such way, we can get a biased version of EB dataset. \n\nTable r4-4: results of training on the biased EB datasets with merics ($W_1$/RME) averaged over three random seeds.\n| bias factor $\\eta$ |     t1      |     t2      |     t3      |     t4      |\n| :----------------: | :---------: | :---------: | :---------: | :---------: |\n|    0 (unbiased)    | 7.951/0.039 | 8.747/0.042 | 9.244/0.019 | 9.620/0.044 |\n|        0.1         | 7.981/0.129 | 8.825/0.118 | 9.219/0.108 | 9.632/0.085 |\n|        0.2         | 8.019/0.072 | 8.908/0.165 | 9.330/0.111 | 9.723/0.181 |\n\nTables r4-3 and r4-4 indicates that when the bias level is relatively low ($\\sigma\\le 0.1$ or $\\eta\\le 0.1$), our VGFM seems robust to the sampling bias. But if the bias level is high, the perfromance degrades. For such cases, we believe that with more biology prior beyond the cell counts integrated into growth modeling, our method could be more robust to the bias obervation. We will include these experimental results and the discussions in the appendix.\n\n### **Q5: The estimation of the Monge map via barycentric projection introduces approximation error, especially when entropy regularization is non-negligible.**\n\nAs stated in Appendix B.1 (Line 574-576, Page 18), we try to choose the entropy regularization coefficient $\\varepsilon$ as small as possible to ensure it can approximate the exact optimal transport and monge map, while ensuring numerical stability and computational efficiency. Extensive experiments (Table 1, Table 2, Page 8) demonstrate the advantages of our VGFM despite the approximation error of exact optimal transport and Monge map.  We acknowledge that theoretical research of this approximation remains an open and intereting direction, and we thank the reviewer for highlighting it.\n\n### **Q6: On limitations**\n\nFor the first point, we will consider adding visualizations of growth error, both spatially and per cell type, in the appendix for the simulated gene dataset. However, for real datasets, growth rates are generally not directly observable, which limits such an analysis.\n\nFor the second point, certain existing databases, such as UniProtKB, record genes that are potentially related to growth. We can leverage this information by comparing the learned $g$ with these known growth-associated genes, and incorporate such comparison as a regularization term during training to regularize or debias the growth modeling process."}]

Author rebuttal:

Meta-review: This paper presents velocity-growth flow matching (VGFM) a method for modeling single-cell population dynamics that handles mass imbalance using semi-relaxed optimal transport. This is a crucial challenge in the single-cell modeling problem. The consensus among reviewers is to accept this paper which I agree with.

While one reviewer (kHRn) expressed concerns during the initial review, that reviewer did not engage in the discussion despite a large number of reminders and will therefore be discounted. Otherwise, the rebuttal discussion went smoothly and improved the opinion of two reviewers. I recommend acceptance at this time.
Decision: Accept (poster)