Delay Flow Matching
Authors: Bolin Zhao
Xiaoyu Zhang
Yuting Dong
Xin Lu
Wei Lin
Qunxi Zhu
Venue: ICLR 2026 Poster (2026)
Keywords: Generative Models
Flow Matching
Delay Differential Equations
Trajectory Intersection
Heterogeneous Distribution Transfer
Primary area: generative models
TL;DR: We propose Delay Flow Matching, a generative framework based on delay differential equations that overcomes flow matching's limits-trajectory intersection, delay dynamics, heterogeneous transfer- while providing universal approximation guarantees.

Abstract: Flow matching (FM) based on Ordinary Differential Equations (ODEs) has achieved significant success in generative tasks. However, it faces several inherent limitations, including an inability to model trajectory intersections, capture delay dynamics, and handle transfer between heterogeneous distributions. These limitations often result in a significant mismatch between the modeled transfer process and real-world phenomena, particularly when key coupling or inherent structural information between distributions must be preserved. To address these issues, we propose Delay Flow Matching (DFM), a new FM framework based on Delay Differential Equations (DDEs). Theoretically, we show that DFM possesses universal approximation capability for continuous transfer maps. By incorporating delay terms into the vector field, DFM enables trajectory intersections and better captures delay dynamics. Moreover, by designing appropriate initial functions, DFM ensures accurate transfer between heterogeneous distributions. Consequently, our framework preserves essential coupling relationships and achieves more flexible distribution transfer strategies. We validate DFM's effectiveness across synthetic datasets, single-cell data, and image-generation tasks.

Published as a conference paper at ICLR 2026
DELAY FLOW MATCHING
Xin Bolin 12, Xiaoyu Zhange.? Yuting Dong"?, Xin Lu Wei Lin 124,5,4,51 Qunxi Zhu
I School of Mathematical Sciences, SCMS, SCAM, and CCSB, Fudan University, China Researcch Institute of Intelligent Complex Systems, Fudan University, China
College of Systems Engineering, National University of Defense Technology. China *State Key Laboratory of Medical Neurobiology and MOE Frontiers Center for Brain Science, Institutes of Brain Science, Fudan University, China SShnnghai Artificial Intelligence Laboratory. China.
[blzhao21, ytdong22] 8m. fudan. edu. on {wlin, qxzhu)@fudan edu cn kxyzhangphd, xin. lu. lou
ABSTRACT
Flow matching (FM) based on Ordinary Differential Equations (ODEs) has achieved significant success in generative tasks. However, faces several in- herent limitations, including an inability to model trajectory intersections, capture delay dynamics, and handle transfer between heterogeneous distributions. These limitations often result in significant mismatch between the modeled transfer process and real-world phenomena, particularly when key coupling or inherent structural information between distributions must be preserved. To address these issues, we propose Delay Flow Matching (DFM), new FM framework based on Delay Differential Equations (DDEs). Theoretically, universal approximation capability for continuous transfer maps. By incorporating delay terms into the vector field, DFM enables trajectory intersections and better captures delay dynamics. Moreover, by designing appropriate initial functions, DFM ensures accurate transfer between heterogeneous distributions. Consequently, our framework preserves essential coupling relationships and achieves more flexible distribution transfer strategies. We validate DFM's effectiveness across synthetic datasets, single-cell data, and image-generation tasks.
show that DFM possesses
INTRODUCTION
Generative modeling a key and rapidly advancing field in machine learning, focusing on learning transformations between different distributions. It underpins wide range of applications across many tasks in diverse fields, including image generation (Nichol et al., 2021). molecule design
(Sanchez-Lengeling & Aspuru-Guziz 2018), and single-cell trajectory inference (Saelens al. 2019). Traditional approaches, such as Variational Autoencoders (Kingma 2013 Rezende et al.
2014), Generative Adversarial Networks (Goodfellow et al.) (2014), and Normalizing Flows (Dinh
et al. 2014 2022 Papamakarios et al., 2021). have achieved notable success, yet they still face challenges like instability and computational inefficiency.
Recently, diffusion models (Sohl-Dickstein et al. 2015: Song & Ermon 2019 Ho et al. 2020
Song et al. 2020aa Yang et al. (202)) have garnered considerable attention for their ability to model complex data distributions via a diffusion process. These models incrementally add noise to data in forward process and then learn to reverse this process order regenerate data samples, which can be understood as learning stochastic dynamics that interpolate between the prior distribution and the target data distribution (Song et al. 2020b 202I). Additionally, Flow Matching (FM), simulation-free method for training continuous normalizing flows (Chen et al. 2018). is proposed to model the transformation between distributions as the flow maps Neural ODES (Chen et al.
2018 Grathwohl et al. 2018). FM achieves efficient training by directly regressing on explicitly
"Equal contribution. Correspnnning authors.

---

Published as a conference paper at ICLR 2026
constructed conditional vector field. Concurrently, the stochastic interpolant (Albergo & Vanden Eijnden )022) and rectified flow (Liu et al. 2023) are introduced, both of which employ flow maps for matching distributions, though from distinct conceptual frameworks.
As discussed above, most existing continuous generative models for distribution transfer rely on ODEs or SDEs. However, ODE-based models have limited representational capacity, restricting their ability to capture a broad range of distribution transfer strategies (Dupont et al.] 2019). In contrast, Delay Differential Equations (DDEs) have been widely adopted model various real-world systems, including neural dynamics (Campbell 2007), electro-optical systems (Chembo Kouomou al. 2005). population dynamics (Lotka] 1925). ). and many biological network motifs (Glass et al.,
2021). where delayed feedback mechanisms naturally give rise to DDE-based formulations. To address the limitations of ODEs in dynamic modeling, Neural DDEs (Zhu et al., (020) and their variants (Ji & Orosz. 2024 Zhu et al. 2022) are proposed, which explicitly incorporate the effects of historical states and thus exhibit significantly enhanced representational capacity. Despite these advancements, however, there is currently no approach that utilizes the probability flows of DDEs for distribution transport.
Contributions. We present Delay Flow Matching (DFM), generalized framework that enables more flexible and precise distribution transport strategies using DDEs. The key contributions of this study are summarized as follows:
Development of DFM: We introduce DFM, novel generative model framework based on Neural DDEs, which overcomes the inherent limitations of ODE-based models by incorporating delay terms the vector field and designing appropriate initial functions. DFM can model a broader range of transport strategies, including those with trajectory intersections, and achieves more precise transport between heterogeneous distributions. It also naturally adapts to the probability flow generated by delay dynamical systems.
2. Theoretical insights: We rigorously prove that DFM can universally approximate any continuous transport map between source and target distributions. In contrast, ODE-based models cannot represent certain simple transport maps, such as those involving trajectory intersections, and cannot achieve exact transport between heterogeneous distributions.
3. Integration with advanced techniques: DFM can be seamlessly integrated with cxisting methods, such as keypoint-guided optimal transport, enabling more effective alignment with known coupling information during transport.
4. Empirical validation: We validate DFM's 's effectiveness on both synthetic and real-world datasets. DFM accurately recovers underlying delay dynamics from snapshot data, sig- nificantly outperforms existing methods in single-cell trajectory inference, and surpasses ODE-based FM in image generation tasks.
Coepbee Insad Pandker
2 PRELIMINARIES
2.1 DDES AND PROBABILITY FLOWS
KPG-OT Cospling
ODE
DBE Dowes Justa Fradion
DDEs are widely used model systems in which the time evolution depends not only on the current state but also on past states, such as those governed by delayed feedback. For time- dependent DDEs with a single delay term, the general formulation is given by: dt)(t) dt
OT Couplina
ODE Coupling
DDE Cospling
Figure 1: Comparison of ODE and DDE coupling. (a) A Gaussian distribution is mapped to itself via KPG-OT, where the red lines and arrows denote the prescribed source-target keypoint pairs, enforc- ing x - -r. The resulting trajectory intersections cannot be realized by ODE flows, whereas the DDE à = - -2x(1 - 1) with constant initial func- tions achieves exact coupling. (b) The transport of ((-2,2) 2) to [ -1) + (1, 3) violates connectivity preservation of ODE flows, forcing mass near 0 into 1). DDEs achieve precise transport by assigning distinct initial functions.
u1(1, a: ar(1 - 7)), (O,
z(h) = (
0-,
where ult, =(t) (t) [0.7] Rdd x Rd Rd is smooth vector field which is abbreviated as utt, x, xx) in the following. vb(h) represents the continuous initial function. We further de- note zb h; wo) as the initial function that takes the value wo time 0, i.e. <b; mo) 26n
2

---

Published as a conference paper at ICLR 2026
For given initial distribution zu (zo) and the corresponding initial functions (t; æu), the above DDE (1) induces the associated probability flows p(a,t/1)) R x
[0. T
+ R satisfying
the delay Fokker-Planck equation (refer to Appendix [B) (Guillouzic et al. 1999
Frank 2005):
0p(z, - { (ul,(a((, 20)) p(a, ar, 7 [bb)]
where p(z, / ( is the probability density given the initial function, p(æ, (; Tr, / 7/45) is the joint probability density of being at at time and at a at time (-T. The specific relationships and distinctions between the probability flows of ODEs and DDEs are elaborated in Appendix A
2.2 OPTIMAL TRANSPORT
Optimal transport (OT) provides a mathematical framework for transforming one probability distri- bution to another in the most cost-efficient manner (Villani et al. 2009; Santambrogio 2015 Chen 2016; Peyré et al. 2019). . It seeks the optimal strategy for reallocating mass while minimizing the transportation cost. Formally, given two separable metric spaces X' and Y with probability measures on X' and // on V, the objective is to determine a transport plan T that minimizes the total transport cost, defined in the Kantorovich problem (Kantorovich) 1942):
[ c(x, y)d y), H1,,,)) v) { (T) P (m)= =v].
C(H,D)=
inf *Ell(u.)) Jxxx
(3)
where c(x, y) denotes the transportation cost of moving a unit mass from x to y. (T) and p (T)) denote the marginal distributions æ with respect to æ and y. In most cases, the transport cost between two points is defined as the squared Euclidean distance. The minimum total transport cost in this scenario corresponds to the squared 2- Wasserstein distance between the two distributions:
(1, inf =EIT(u.)) / | - yd y).
2.3 KEYPOINT-GUIDED OPTIMAL TRANSPORT
Traditional OT, which focuses solely on minimizing transport costs, often neglects other crucial constraints, resulting in suboptimal matching strategies. Thus, it is important to use small set of well-matched source-target keypoint pairs, K = [(Tk } to semi-supervise the transport strategy, especially when inherent structures and features of data need be preserved.
k==
To ensure correct transport with these keypoints, Reypoint-guidde Optimal Transport (KPG-OT) is proposed (Gu et al. 2022 2023), which guarantees that the relationships between each data point and the keypoints are preserved during transport. Formally, the objective can be expressed as:
inf / #Ell(H.N) Jxxy (x. y)d(w (x. ii(p v) = { P (w o म)) = μ, Pll (w o R) = =w] (5)
where म wxx y) (xx y) is the teypoint-masked transport plan, denotes the mask function.
For pair of keypoints (æx, yh E K. the mask function is defined as: w w yı ) = 1, w(x;, y) 0 for y j yk, ww. yı) Obviously, the mask function defined above exactly preserves the matching of the keypoints in K.
for x of 2h. and w(x y) if neither æ nor y matches any keypoint.
The cost function in Eq. 5 is defined g( y) ([("(2) Rd (y)). where d denotes the Jensen- Shannon divergence. (x) E 1) (resp. R² (y) E (0, 1)x) captures the relationship between æ (resp. y) and all source points (resp. target points) in K. with its k-th dimension representing the relation score between æ (resp. y) and the k-th source keypoint D& (resp. target keypoint yı:):
e o((0,)))
R (x)
R (v)
(6)
DK 0// DEE
K
EK e-y(yum)/
where is the temperature hyperparameter. Note that g quantifies the similarity distance between the two relationship vectors outlined above. Therefore, minimizing the cost function in Eq. (5) effectively promotes the preservation of the relationship between each point and the keypoints.
3

---

Published as a conference paper at ICLR 2026
2.4 ODE-BASED FLOW MATCHING
Flow Matching (FM) (Lipman et al., 2022) is a simulation-free training method for CNF, which trains the parameterized vector field v(I, x; 0) by regression on the target vector field u(I, x)) that generates the probability paths p(æ, 1) with marginals Po = 4o and P1 = 91- The training objective is
CFM = (.,)( ||
[=( x) - at)x
(7)
where ( ~ ((0,) 1) follows the uniform distribution.
However, the explicit forms of =tt, ax) and p(ar, are intractable to compute. To address this, Conditional FM (CFM) (Tong et al. 2023a; Pooladian et al.] (023) introduces a latent variable construct the target probability path as a mixture of conditional probability paths p(æ, t) By(=) [p(,, []=)] and modifies the training objective to:
(8) where ult,, #|]= denotes the conditional vector field generating the conditional probability path p(a, t|z). Then, it can be shown that the gradient w.r.t & of the CFM objective ( is the same as
(O) (r..(()) [|;((, -((,,=;;))
that of the FM objective 77). Typically, we set the latent condition == (za, 21), with g(z) being
the coupling between distributions yo (2o) and q (x)1 ). based on a coupling strategy such as OT or
KPG-OT. The conditional probability path is then modeled
Gaussian flow:
p(æ, t|z> = (x|tax) + (1 t) o²),
(9)
with the corresponding conditional vector field:
utt, x (2) = 31 TD.
(10)
LIMITATIONS OF ODE-BASED FLOW MATCHING
Current continuous-time generative models utilize the parameterized vector field of ODEs to facilitate distribution transport. However, the inherent constraints of ODEs limit their ability to model specific transport strategies, leading significant discrepancies between the modeled and true transport processes. In this section, we theoretically demonstrate the limitations of the ODF-bbaed FM frameworks. Note that the following propositions are stated within the classical FM framework Lipman et al. 2022 Liu et al. 2023). J. Recent ODE-based FM extensions enhance representational capacity and mitigate these limitations; we discuss these variants further in Section
Proposition 3.1 (Restriction on trajectory intersections). Suppose that the target transport strat- egy. guided by certain key points specific coupling constraints, inevitably induces trajectory intersections in the augmented phase space (t, x)) during the transport process. If the ODE-based FM is governed by vector field that Lipschitz continuous in and continuous in then the corresponding flow map cannot precisely preserve the transport strategy.
Remark 3.2. As demonstrated by Liu et al.) (2023), ) when trajectory intersections occur in the transport strategy, FM tends to learn a flow, with the targets of the trajectories being "rewired" at the intersection points. As result, the given transport strategy cannot be accurately preserved. A simple example is shown in Fig. (a).
Proposition 3.3 (Heterogeneity in distributions). Assume that the source (resp., target) distribution qu (resp., 91) is supported on M (resp., N) disjoint compact sets {U'IM, J (resp., {U, N Tz= 1). where
each Up E RRd (resp., U} Rdd) is path-connected set with non-zero measure and go (0)))) > 0
(resp., Q1 (an) > for all an E U. UP (resp., E U IfM < N, then for any transport map T: UM U + U U} satisfying T490- = q1, the following holds: if the flow map induced by a Neural ODE is equal to the transport map T, then the associated vector field cannot be Lipchitz continuous.
Remark 3.4. The conclusion emphasizes that accurate transport using an ODE-based FM requires regularity assumptions on distributions, which are often violated in real-world scenarios. For example, in single-cell dynamics, cells differentiate from one type into multiple types over time, leading to heterogeneous distributions before and after differeetiiaion. A simple example is shown in Fig. (b).
Moreover, when snapshot data are generated by delay dynamical system, based FM fails to recover the true vector field with the delay term, resulting in inaccurate distribution transfer, interpolation, and extrapolation predictions.
4

---

Published as a conference paper at ICLR 2026
4 DELAY FLOW MATCHING
To overcome the limitations of ODE-based FM, we introduce the Delay FM (DFM) framework, a new class of generative models based on Neural DDEs for distribution transport. Since the vector field incorporates delay terms, DFM allows trajectory intersections, enabling more accurate transport for tasks requiring keypoint-guided strategies. Additionally, by designing appropriate initial functions, DFM addresses singularities caused by distributional heterogeneity. A comparison between DDE- based and ODE-based FM under various scenarios shown in Table
Table I:
Comparison of ODE- and DDE- based FM. DFM can handle not only KPG-OT coupling and trajectory intersections but also distribution heterogeneity by employing diverse initial functions. Furthermore, DFM can accurately model and recover the dynamics of delay dynamical systems.
Methods
KPG-OT
Intersection
Heterogeneity
Delay Dynamics
FM
CFM
X
DFM(C) DFM(D)
/
Initial distribution
Target distribution Keypoint Pairs
Generated distribution
4.1 FORMULATION
Consider a target probability flow p(x, t), gen- erated by vector field with single delay term ulI, æx, x. and initial functions xb ~ (1b). The vector field and initial functions naturally define a joint probability flow at time t and t-7, denoted as p(x, t; ær, r), which satisfies
I pla, li 7) r)dx, p(xx,t) 1). (11)
(c)
(d)
Note that f(x,t) and p(æ, t; t r) can be modeled as mixture of conditional distribu- tions p(ær, f/ab) and p(x, w. 10 T(b)), respec-
tively, as follow:
Figure 2: (a. b) show the generation results of CFM (a) and DFM (C) (b) under the KPG-OT coupling Gaussian mixture data. (c. d) show the
p(a,t) t) = E (2) t/bb,
p(a, t; xr, 1- = ( 1n x., r //bb) generation results of CFM (c) and DFM (D) (d) for
(12) heterogeneous source and target distributions.
where p(æ, 1/bb) and p(æ, 4; (r, (1bb)) represent the probability flow generated by eul,, æe, xx)) and the initial function 1b and they satisfy the delay Fokker- Pekker-Planck equation (2). Based on the joint probability flow, we aim to learn a parameterized vector field vtt, a, æ,; 8). which can generate the target probability flow, by minimizing the following regression objective against the target vector field:
LDFM (0) = Ee,, (t) (/ x, ær; 0) ull, α, x111
(13)
In most cases, both p(x, t; xt,, - - ab) and u(t, x, x7) are computationally intractable. To address this, inspired by CFM, we introduce latent variable and further decompose the target probability path and joint probability path conditioned on the initial function into mixtures of simple probability paths and simple joint probability paths conditioned on both 2b and z, respectively, as follow:
p(z, t/bb) = E,,,) t|z, 4b), p(x, x,, t - 7 (bb) = (( t; xr, T|2, 1b),
(14)
where p(x, 2, xb) = p(a, x,, - T|2, dd We can define the marginal vector field by marginalizing over the conditional vector field u(t, x, 2z (z) follow:
2a æc, xx 2 p))pppp(p))))))) x m. t T|., bb)
u(t, x, x,,, 1g(z))
(15)
p(æ, ,,, - xb))
where u(t, x, xx)22) represents the conditional vector field that generates p(z, f/2, 2b) and p( 1r, L 7/2, 4b). Under this setup, we can derive the following result.
5

---

Published as a conference paper at ICLR 2026
Proposition 4.1. The marginal vector field u(t, x,,;,) given in Eq. (5), together with the lected initial function +, generates the probability path p(a,t|pb and the joint probability path p(a,t; t, r|pb) in Eq. (4).
Then, we can train the parameterized vector field by minimizing the following Delay Conditional FM (DCFM) objective:
(D) = "tt.:///.aaaaa................/. 2x;;0))))) 21,
(16)
because the following proposition holds.
Proposition 4.2. Given that pf t; 7|5)) for all æ, œ, Id and [0. 1}, up to a constant independent of 0, LDCMM(0) and (DFM(() are equal, which further implies that VoLocrM VoLoru (
4.2
SELECTION OF THE LATENT VARIABLE
4.2.1
TRANSPORT BETWEEN TWO MEASURES
As outlined in Section 2.4, ODE-based FM typically selects latent variables as two endpoints = (Fo, Ø1) jointly sampled from the source and target distributions. In contrast, DFM requires the construction conditional joint probability density path for (t)) and a 7) based on the latent variable. Hence, we define the latent variable as entire path ((; æo. æ1) connecting and æ1-
Formally, we define the probability distribution of the latent variable g/r((;; øo. æ1)] := o (Co, æi (r)) No, 21). where π (xo, 21) represents the OT KPG-OT coupling between the source and target distributions, and P(y; Fo, x1) denotes the path measure pinned at 2ig and æ1. In practice, we can simply construct the path measure Dirac Delta distribution at given path Y connecting 2o and Ø1:
(17) where 7 is constructed using specific interpolation method. For instance, it can be taken as the linear interpolation 7;- (1 t) + fæ or alternatively, as geodesic interpolation on the data manifold, based on appropriate manifold learning techniques.
P(y; To, x1) 8(> >"), % zo. Si 11.
After sampling a path (t: æo. 21) ~ q(z), the conditional joint probability density naturally degen- erates into the following form:
p(z, t; æ,, t +|>) = 8[æ (t)]][[z, 7(t - -T]]
(18)
which is exactly generated by the conditional vector field:
@y(1 x0
ult, x. By (nt;
(19)
Ōt
4.2.2 TRANSPORT BETWEEN MULTIPLE PROBABILITY MEASURES OVER TIME
For tasks with multiple target probability distributions over time {9,, 1 where to = 0, t. = T. define the latent variable sampled from
(t: { 1) trajectory passing through Øi, time ty.
J-1 n ( Rtmn):
a by (t: {xi, 1)) = P (7s f2i, D
II
(20)
where (xx, 1/=0 I π (x R...). ). with ot (Zt,, :t Rt,.)) representing the OT KPG-OT coupling between adjacent probability distributions at time t; and f.41, and P(y; (x1, } denotes the path measure pinned at (mx,, 1/ Similarly, we construct the path measure as:
P ( (Zi,, (/o)) = 6(7 - -y"), 7i, = Zty.
(21)
where 7 represents a path passing through Tiy at time ty. which can be constructed using the cubic spline (CSpline) interpolation method. With this path as the latent variable. the resulting conditional joint probability path and conditional vector field are identical to those in Eq. (18) and Eq. (19). respectively.

---

Published as a conference paper at ICLR 2026
4.3 SELECTION OF THE INITIAL FUNCTION
4.3.1 DFM WITH CONSTANT INITIAL FUNCTIONS
In general, the initial function can be simply chosen as a constant (DFM(C)), i.e. where 26 (t; zo) = Zy for & T-,,0 0) and 2 ~ 9o- Under this setting, we can theoretically prove that DFM can approximate any continuous transport strategy to arbitrary precision.
q 5( -bb")
Proposition 4.3 (Universal approximating capability of DDE's flow map). For any given > 0 and continuous transport map F _ Rd, which push-forward the source distribution qo to the targer distribution qı i.e. $# = q1, if there exists a neural network f(x; f) such that [1; 0) ( term (t,, x 2r; (), where the corresponding flow map G ; 0), under the constant initial function condition, satisfies [|G( 0) - F(z)]l < for all a E RR..
æ]ll < for all a E R then we can construct vector field with single delay
This indicates the exceptional representational capacity of DFM, enabling the modeling of wider range of transport processes than ODE-based FM.
4.3.2 DFM WITH DIVERSE INITIAL FUNCTIONS
In Section 3 we rigorously prove that ODE-based FM cannot effectively handle tasks with heteroge- neous source and target distributions. Here, we demonstrate that DFM addresses this by designing diverse initial functions (DFM(D)).
Specifically, we employ clustering methods, such as Gaussian Mixture Model, DBSCAN, to partition the source dataset Xo 40 (resp. target dataset X1 q1) into M (resp. N) mutually exclusive
subsets, denoted as X x (resp. x X (N)). We can assign a normalized mass each
subset (m) (resp. X x()) as Pa (m) = X 0 1/ /uu (resp.. (m)
| 1/ l/) If the endpoints
(aro 1 of sampled trajectory y are drawn from x-()) and x x-(). we assign an initial function
2b which has constant time derivative Cmn.
dab (t; zu) Cmn, tinn 20; = 2o.
(22)
dt
where t E [-T,0], OJ. In this case, (xb) is discrete distribution which satisfies:
M
(v mmn) Po
(m)
Σ
qd (vb mm) = (n)
(23)
p(æo. 0/2b = (. = 20 = vin) represents empirical data distribution available as data
points in x(m²) whose corresponding transport target is in X (()) By coupling the source and target
data through OT or KPG-OT, we can obtain the transport target for any initial point. This enables the
construction of qo (v) as described above, from which we can sample obtain the corresponding
initial function.
In summary, we design distinct initial functions for different subsets of the source and target data, guid- ing the vector field from different initial subsets to corresponding target subsets, thereby effectively handling distributional heterogeneity.
4.4 GENERATION PROCESS BASED ON NEURAL DDES
After training, we sample data from @o and initial functions from ( (v)), then generate the target data by solving the forward pass of the trained vector field with a single delay term using a piecewise ODE solver, in Neural DDEs (Zhu al. 2020).
RELATED WORKS
A number of recent methods have been proposed to address the limitations of ODE-based FM. To enable crossing trajectories, Constant Acceleration Flow (Park et al, 2024) models acceleration and
conditions it initial velocity; Hierarchical Rectified Flow (Zhang al. 2025) achieves trajectory
7

---

Published as a conference paper at ICLR 2026
crossings by hierarchically coupling multiple ODEs; and Augmented Bridge Matching (De Bortoli
ct al. 2023) incorporates information about the initial point directly into the vector field. Diversified
Flow Matching (Shrestha & Fu 2025 takes different perspective by training interpolation strategies
that construct non-intersecting target trajectories, thereby indirectly mitigating trajectory-crossin issues. To address the singularities caused by distributional heterogeneity, Switched Flow Matching
(Zhu & Lin 2024) and Variational Rectified Flow Matching (Guo & Schwing, (025) introduce additional latent variables that allow the model to represent multi-modal transport paths, alleviating the degeneracies that arise under heterogeneous target distributions.
The methods discussed above introduce various extensions classical FM, enabling the model to accommodate more transport strategies. While these approaches either allow trajectory intersections during the transport process or alleviate singularities arising from distributional heterogeneity, most of them cannot simultaneously resolve both issues. Moreover, they remain confined to the ODE-based generative framework, typically resolving limitations by introducing additional conditioning variables into the vector field or by redesigning the target transport process.
In contrast, our proposed DFM departs from the ODE paradigm and introduces a fundamentally new DDE-based generative framework. By modeling the vector field directly in the original phase space, without introducing auxiliary latent variables or specially engineered transport paths, DFM inherently accommodates both trajectory intersections and distributional heterogeneity within unified mechanism. Furthermore, DFM naturally captures transport dynamics governed by intrinsic delays, a capability that is entirely out of reach for those ODE-based extensions.
6 EXPERIMENTS
We demonstrate the advantages of DFM over ODE-based FM across various tasks, including recon- structing delay dynamics from snapshots, inferring differentiation trajectories from single-cell RNA sequencing (scRNA-seq) data, and image generation. The experimental details and additional results can be found in Appendix D A systematic ablation study is shown in Appendix
6.1 SIMPLE SYNTHETIC ILLUSTRATIONS
To more clearly illustrate the DFM workflow, we present two simple and concrete running examples, demonstrating the procedure of DFM(C) and DFM(D), respectively.
2-d Gaussian mixture dataset. As shown in Fig. 2 (a.b), aim transport two-component Gaussian mixture to its reflection-symmerric counterpart while respecting 6 keypoint pairs (i.e., the optimal mapping + =x). Using KPG-OT, we first obtain pointwise couplings and construct the linear interpolation paths = x 0 + (1 - t) x œ1, which unavoidably produce trajectory crossings. 1, and train the delay vector field. Notably, this toy problem admits an exact DDE solution: & = -2x(t - 1) -2zg. The learned DFM(C) exhibits the same behavior: trajectories select different directions at the crossing by leveraging historical states, yielding accurate keypoint-respecting transport.
We select the constant initial function with
1-d example with distributional heterogeneity. As shown in Fig. (c.d), we aim to realize the transport 1) - (43, 2) + (2, 3). Notably, we can use J. = xx() xtt (t 1) to achieve precise OT transport by assigning different initial functions: x(t) = xo - t, xo [-1, 0j and @ at() xo + f, 20 (0, 1]. During training, we employed the OT coupling. which maps [-1, to
[-3.-2] and (0, 1] to [2, 3), and construct the linear interpolation paths. In applying DFM(D), initial function ① is selected for 21 E [-3, 2], while @ is chosen for T1 E 3). with 7 = 1. The trained
DFM(D) achieves accurate transfer between heterogeneous distributions (Fig. 2 (d))
6.2 SYNTHETIC DATASETS OF DELAY DYNAMICAL SYSTEMS
Biological autoregulation motif. DDEs provide simplified framework for modeling biological network motifs (Glass et al. 2021). We consider recovering the delay dynamics from snapshot data
generated by the autoregulation model: i(t) [+="([--) with T 1, n. 2, 17 = 5, which produces damped oscillations (Fig. 3 (a.b)). Using 1, 000 initial values sampled from 20,2,, 1.2) and constant
initial functions, we generate trajectories and collect snapshots at /=0, 1,..., for training. During
8

---

Published as a conference paper at ICLR 2026
20
30
15
90 8
10
10
63
615
Figure 3: Comparison of predicted trajectories of CFM and DFM trained on snapshots of delay dy- namical systems. (a-b) Interpolation and extrapolation (insets) results on the biological autoregulation motif dataset using KP-CFM (a) and KP-DFM(C) (b), both with two keypoints. (c-d) Interpolation results on the spiral DDE using OT-CFM (c) and OT-DFM(CC) (d).
training, KPG-OT (KP-) matches minibatches at adjacent time steps using two known keyoints. After training, forward integration from (=0 shows that KP-CFM fails to recover the dynamics, while KP-DFM(C) with 7=1 accurately captures and extrapolates the oscillatory behavior (insets in Fig. (a,b)). Results with other delays and error metrics are provided in Appendix D.I
and
Spiral DDE. We next consider a 2-d DDE (Zhu et al. 2020 ((t) A (an)(((t) + xtt-r)] with 7=0.5 and A E 2 which produces crossing spiral trajectories (Fig. 3 (c.d)). Snapshots are taken every 0.05 in t E [o, 0.6] and coupled using minibatch-OT (OT-) between adjacent steps. After training. OC-DFM(() with 7=0.5 successfully reproduces the dynamics, while OT-CFM fails around
the crossing region. DFM also generalizes well across range of values (Appendix D.2 Fb.
6.3 TRAJECTORY INFERENCE SINGLE-CELL
We investigate the inference of differentiation trajectories from real scRNA-seq data, where het- erogeneity increases as cells transition from a single type to multiple types during development, complicating modeling with ODE-based approaches. To evaluate performance, we train on all time points except one intermediate time, then sample from the initial distribution and perform forward integration to predict distributions at each time point. Predictions are assessed via unsupervised leave-one-out validation (L) against the held-out intermediate distribution and supervised final-time validation (F) against the final distribution. Trajectory accuracy is measured using the 2- 2-Wasserstein distance (W2) and Maximum Mean Discrepancy with a Gaussian kernel (MMD(G)).
Table 2: Trajectory inference results on the scRNA-seq dataset in mouse hematopoiesis and the single-cell qPCR iPSC dataset with bifurcation. All results are averages of 10 runs.
Mouse hemulopoiesis
qPCR iPSC
Methods
W2(L)
MMD())
W2(F)
MMD(F)
(2())
MMD(L)
W2(F)
MMD(F)
TIGON
0.519
0.563
0,264
0.155
0.733
0.791
0.695
0.405
MIOFlow
0.514
0.629
0.220
0.056
0.770
1.039
0.345
0.155
OT-CFM
0.378
0.357
0.192
0.047
0.579
0.492
0.226
0.030
OT-DFMM()
0.379
0.384
0.136
0.021
0.553
0.447
0.234
0.041
OT-DFM(D)
0.372
0.341
0.095
0.010
0.532
0.399
0.213
0.027
Mouse hematopoiesis dataset. We evaluate DFM on scRNA-seq data from mouse hematopoiesis (Weinreb et al. 2020), focusing on cells differentiating into neutrophils (Neu) and monocytes (Mo) Day 2, 4, and 6. Models are trained on Day 2 and data; after training, Day 2 samples are integrated forward to predict distributions at Days 4 and 6. Geodesic interpolation is used to construct trajectories. Due to the heterogeneity between pre- and most-differentiation distributions, trajectories inferred by ODE-based models (TIGON (Sha et al. 2024), MIOFlow (Huguet et al. 2022 OT-CFM) partially deviate from the data manifold, falling between Neu and Mo fates (Fig. 4 (b-e)). In contrast,
OT-DFM(D), using distinct initial functions for Neu and Mo (Appendix D.3). preserves alignment
with the data manifold (Fig. (D) and achieves more accurate predictions (Table
qPCR iPSC dataset. We further apply DFM to model the bifurcation process of iPSCs in cardiomy- ocytes (Bargaje et al. 20. Data from Day 2, 3, 4, and 5 are selected, with bifurcation observed from Day 3, where progenitor cells differentiate into mesodermal (M) and endodermal (En) fates. Models are trained on all time points except Day 3. After training, samples from Day 2 are forward-intend -nwrrddintegratedd to infer trajectories and predict distributions Day 3 and 5. Due to distributional heterogeneity the

---

Published as a conference paper at ICLR 2026
bifurcation point, ODE-based methods and OT-DFM(C) produce trajectories that partially misalign with the data manifold (Fig. 4 (h-k)). In contrast, OT-DFM(D), by assigning separate initial functions for each fate (Appendix D.4), preserves trajectory alignment with the corresponding manifolds and yields more accurate predictions at both Day and 5 compared to other methods (Table
60
828
000
(0
2
Figure 4: Comparisons between predicted trajectories of MIOFlow (b. h), TIGON (c. i), OT-CFM (d, i). OF-DM(()) (c. k), and OT-DFM(D) (f. D) on scRNA-seq dataset in mouse hematopoiesis (a-f) and qPCR iPSC dataset (g-l). (a) and (g) illustrate the true data distribution.
6.4
IMAGE GENERATION
We design a Semi-paired Image-to-Image Translation task on the MNIST dataset, where the source domain consists of original images, and the target domain includes their negative images (the normalized pixel val- ues transformed via o + 1- To provide
MNIST dataset.
Table 3: Comparison of FID between KP-CFM
(T 0) and KP-DFM(C) with different time delay on the MNIST dataset.
(CFM)
0.125
0.250
0.500
1.000
HDD
45.020
28.497
11.747
12.653
12.031
partial supervision, 10% of the training data are paired with their negative counterparts as keypoints. During training, minibatches are independently sampled from the source and target distributions, and coupled via KPG-OT (KP-). The task aims transform source images into their negatives, where, under linear interpolation, all transport paths intersect at an image with uniform pixel values of 0.5. While KP-CFM struggles with this transformation, KP-DFM(C) effectively handles it. As shown in across varying delays т, KP-DFM(C) consistently achieves lower FID scores than KP-CFM.
Table 3
CIFAR-10 dataset. We further evaluate DFM on the CIFAR-10 10 dataset under a more general generation setup (Zhu & Lin 2024), where the source distribution is a two-component Gaus- sian mixture. DFM employs trainable initial functions with distinct constant time derivatives, enabling generation from different mixture com- ponents to specific image classes. We compare CFM and DFM(D) under two coupling strategies: independent coupling (I-) and OT coupling (OT-).
Table 4: Comparison of FID) between CFM and
DFM on the CIFAR- 10 dataset.
NFF
20
40
Adap.
I-CFM
108.291
94.629
91.404
90.254
88.306
OT-CFM
78.165
27.512
16:109
12.026
6.162
I-DFM(D)
54.064
18.248
11.429
0.008
4.980
OT-DFM(D)
51.222
18.598
11.891
9.287
5.191
As shown in Table 4 I-CFM struggles with mode heterogeneity, while I-DFM(D) generates higher- quality images. With OT coupling. OT-DFM(D) consistently outperforms OT-CFM, especially when the number of function evaluations (NFE) is small. Additional results are provided in Appendix D.6
7 CONCLUSION
We introduce DFM, a novel continuous- timimuouu--imeeeeee generative framework based on Neural DDEs. Through theoretical analysis, we highlight the limitations of ODE-based generative models, particularly their inability to capture certain transport strategies and preserve critical coupling information. In contrast, DFM offers universal approximation for arbitrary continuous transport strategies, addressing these shortcomings effectively. DFM also overcomes the challenge of transport between heterogeneous distributions by incorporating rask-specific initial functions. Furthermore, it is naturally suited for modeling delay dynamical systems, feature beyond the capability of ODEs. Extensive experiments on both synthetic and real-world datasets demonstrate that DFM achieves significantly more precise and versatile distribution transport strategies compared to FM.
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGMENTS
Q. Zhu supported by the National Science and Technology Major Project for Brain Science and Brain-like Intelligence Technology (No. 20252D0215700), and the National Natural Science Foundation of China (Nos. T2541024, 62406072, 12171350). W. Lin is supported by the NSFC (Nos. 19225103, 12531018), the STCSM (Nos. 2021SHZDZX0103 22IC1402500, 22JC1401402, 25JS28 10400), and the SMEC (Nos. 2023ZKZI044, 2023KFJ105-72), X. National Science and Technology Major Project for Brain Science and Brain-like Intelligence Tech- nology (No. 20252D0215700), the National Natural Science Foundation China (Nos. 72025405, 72421002, 92467302, 72474223, 72301285), the Hunan Science and Technology Plan Project (Nos. 2023J140685, 2024RC3133), and the Major Program of Xiangjiang Laboratory (No. 24XJJ- CYJ01001).
supported by the
REFERENCES
Michael Albergo and Eric Sunden-Eijnden Building normalizing flows with stochastic interpolants. Arxiv Preprini Arxiv:2209. 15571, 2022.
Uwe an der Heiden. Delays in physiological systems. Journal of Mathematical Biology, 8(4):
345-364, 1979.
Rhishikesh Bargaje, Kalliopi Trachana, Martin N Shelton, Christopher S McGinnis, Joseph X Zhou, Cora Chadick, Savannah Cook, Christopher Cavanaugh, Sui Huang, and Leroy Hood. Cell population structure prior to bifurcation predicts efficiency of directed differentiation in human induced pluripotent cells. Proceedings of the National Academy of Sciences, 114(9):277--2276 2017.
Charlotte Bunne, Geoffrey Schiebinger, Andreas Krause, Aviv Regev, and Marco Cuturi. Optimal transport for single-cell and spatial omics. Nature Reviews Methods Primers, 4(1):58, 2024.
Sue Ann Camphell. Time delays in neural systems. In Handbook of brain connectivity, pp. 65-90. Springer, 2007.
Yanne Chembo Kouomou, Pere Colet, Laurent Larger, and Nicolas Gastaud. Chaotic breathers in delayed electro-optical systems. Physical Review Letters, 95(20):203903, 2005.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances Neural Information Processing Systems, 31, 2018.
Yongxin Chen. Modeling and control of collective dynamics: From Schrödinger bridges optimal mass transport. PhD thesis, University of Minnesota, 2016.
Valentin De Bortoli, Guan-Horng Liu, Tianrong Chen, Evangelos A Theodorou, and Weilie Nie. Augmented bridge matching. arXiv preprint 1.Xiv:231L.0697,, 2023.
Laurent Dinh, David Krueger, and Yoshua Bengio. Nice: Non-linear independent components estimation. arXiv preprint arXiv:1410.8516 2014.
Laurent Dinh, Jascha Sohl-Dickstein, and Samy Bengio. Density estimation using real nvp. In International Conference on Learning Representations, 2022.
Emilien Dupont, Arnaud Doucet, and Yee Whye Teh. Augmented neural odes. Advances Neural Information Processing Systems, 32, 2019.
Patrick Esser, Sumith Kulal, Andreas Blattmann, Rahim Entezari, Jonas Müller, Harry Saini, Yam Levi, Dominik Lorenz, Axel Sauer, Frederic Boesel, et al. Scaling rectified flow transformers for high-resolution image synthesis. In International Conference on Machine Learning. 2024.
TD Frank. Delay fokker-planck equations, novikov's theorem, and boltzmann distributions small delay approximations. Physical Review Statistical, Nonlinear, and Soft Matter Physics, 72(1): 011112, 2005.
11

---

Published as a conference paper at ICLR 2026
David S Glass, Xiaofan Jin, and Ingmar H Kiedel-Kruse Nonlinear delay differential equations and their application to modeling biological network motifs. Nature Communications, 12(())1788, 2021.
Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. Generative adversarial nets. Advances in Neural Information Processing Systems, 27, 2014.
Will Grathwohl, Ricky TQ Chen, Jesse Bettencourt, Ilya Sutskever, and David Duvenaud. Ffjord: Free-form continuous dynamics for scalable reversible generative models. Arxiv Preprint Arxiv:1810.01367, 2018.
Xiang Gu, Yucheng Yang, Wei Zeng, Jian Sun, and Zongben Xu. Keypoint-guided optimal transport with applications in heterogeneous domain adaptation. Advances in Neural Information Processing Systems, 35:14972-11998 14985, 2022.
Xiang Gu, Liwei Yang. Jian Sun, and Zongben Xu. Optimal transport-guided conditional basore-based diffusion model. Advances in Neural Information Processing Systems, 36:36540-36552, 2023.
Steve Guillouzie, Ivan L'Heureuxx and André Longtin. Small delay approximation of stochastic delay differential equations. Physical Review E, 59(4):3970, 1999.
Pengsheng Guo and Alex Schwing. Variational rectified flow matching. Forty-second International Conference on Machine Learning, 2025.
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. Advances in Neural Information Processing Systems, 33:6840-6851, 2020.
Guillaume Huguet, Daniel Sumner Magruder, Alexander Tong, Oluwadamilola Fasina, Manik Kuchroo, Guy Wolf, and Smita Krishnaswamy. Manifold interpolating optimal-transport flows for trajectory inference. Advances in neural information processing systems, 35:29705-297188 2022.
Xunbi A. Ji and Gábor Orosz. Trainable delays in time delay neural networks for learning delayed dynamics. IEEE Transactions on Neural Networks and Learning Systems, pp. 1-11, 2024. doi: N0.11009TTNNLSSS222.33707000
Leonid V Kantorovich. On the translocation of masses. In Dokl. Akad. Nauk. USSR (NS), volume 37, pp. 199-201, 1942.
Kacper Kapušniak, Peter Potaptchik, Teodora Reu, Leo Zhang, Alexander Tong, Michael Bronstein, Avishek Joey Bose, and Francesco Di Giovanni. Metric flow matching for smooth interpolations on the data manifold. arXiv preprint arXiv:2405. 14780, 2024.
Diederik P Kingma. Auto-encoding variational bayes. arXiv preprint arXiv:1312.6114.
Yaron Lipman, Ricky TQ Chen, Heli Ben-Hamu, Maximilian Nickel, and Matt Le. Flow matching for generative modeling. arXiv preprint arXiv:2210.027477 2022.
Xingchao Liu, Chengyue Gong, et al. Flow straight and fast: learning to generate and transfer data with rectified flow. In International Conference on Learning Representations, 2023.
AJ Lotka. Elements physical biology. Williams and Wilkins, 1925.
Huanfei Ma, Bing Xu, Wei Lin, and Jianfeng Feng. Adaptive identification of time delays nonlinear dynamical models. Physical Review E. Stattisstiaal,, Nonlinear, and Soft Matter Physics, 82(6): 066210, 2010.
Huanfei Ma, Siyang Leng, Chenyang Tao, Xiong Ying, Jürgen Kurths, Ying- Cheng Lai, and Wei Lin. Detection of time delays and directional interactions based on time series from complex dynamical systems. Physical Review E, 96(1):01 ::21122,, 2017.
Alex Nichol, Prafulla Dhariwal. Aditya Ramesh, Pranav Shyam, Pamela Mishkin, Bob McGrew, Ilya Sutskever, and Mark Chen. Glide: Towards photorealistic image generation and editing with text-guided diffusion models. arXiv preprint arXiv:21 2021.
12

---

Published as a conference paper at ICLR 2026
George Papamakarios, Eric Nalisnick, Danilo Jimenez Rezende, Shakir Mohamed, and Balaji Lakshminarayanan. Normalizing flows for probabilistic modeling and inference. Journal of Machine Learning Research, 22(57):1-644 2021.
Dogyun Park, Sojin Lee, Sihyeon Kim, Taehoon Lee, Youngjoon Hong, and Hyunwoo Kim. Constant acceleration flow. Advances in Neural Information Processing Systems, 37:90030-90060, 2024.
Gabriel Peyré, Marco Cuturi, et al. Computational optimal transport: With applications data science. Foundations and Trends in Machine Learning, (1(5-6):355-06 2019.
-eam-Alexandre Pooladian, Heli Ben-Hamu, Carles Domingo-Enrich, Brandon Amos, Yaron Lipman, and Ricky TQ Chen. Multisample flow matching: Straightening flows with minibatch couplings. arXiv preprint arXiv:2304. 44772, 2023.
-laa-Alexandre Pooladian, Carles Domingo-Eninn Enrrgo,-Errch,,,, Ricky TQ Chen, and Brandon Amos. Neural optimal transport with lagrangian costs. In Uncertainty Artificial Intelligence, pp. 2989-3003. PMLR, 2024.
Danilo Jimenez Rezende, Shakir Mohamed, and Daan Wierstra. Stochastic backpropagation and approximate inference in deep generative models. In International Conference on Machine Learning, pp. 1278-1286. PMLR, 2014.
Fathalla A Rihan et al. Delay differential equations and applications to biology. 2021.
Wouter Saelens, Robrecht Cannoodt, Helena Todorov, and Yvan Saeys. A comparison of single-cell trajectory inference methods. Nature Biotechnology, 374555447-554, 2019.
Benjamin Sanchez-Lengeling and Alán Aspuru-Guzik Inverse molecular design using machine learning: Generative models for matter engineering. Science, (61(0)))) 305, 2018.
Filippo Santambrogio. Optimal transport for applied mathematicians. Birkäuser, NY, 55(58-33):94, 2015.
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.
Sagar Shrestha and Xiao Fu. Diversified flow matching with translation identifiability, In Forty-second International Conference Machine Learning, 2025.
Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan. and Surya Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. International Conference on Machine Learning, pp. 2266-2265. PMLR, 2015.
Jiaming Song, Chenlin Meng, and Stefano Ermon. Denoising diffusion implicit models. Arxiv Preprint Arxiv:2010.02502, 2020a.
Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. Advances in Neural Information Processing Systems, 32, 2019.
Yang Song. Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole. Score-based generative modeling through stochastic differential equations. arXiv preprint arXiv:2011.13456. 2020b. 2020b..
Yang Song, Conor Durkan, Iain Murray, and Stefano Ermon. Maximum likelihood training of score-based diffusion models. Advances in neural information processing systems, 34:1415-1428, 2021.
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. arXiv preprint arXiv:2302. 230.00200.. 2023a.
Alexander Tong, Nikolay Malkin, Kilian Fatras, Lazar Atanackovic, Yanlei Zhang, Guillaume Huguet, Guy Wolf, and Yoshua Bengio. Simulation- schrödinger bridges via score and flow matching. arXiv preprint arXiv:2307.0888 03672, 2023b.
13

---

Published as a conference paper at ICLR 2026
Cédric Villani al. Optimal transport: old and new, volume 338. Springer, 2009.
Caleb Weinreb, Alejo Fodriiuuz--raticelll) Fernando D Camargo, and Allon M Klein. Lineage tracing on transcriptional landscapes links state to fate during differentiation. Science, 367(6479): eaaw3381, 2020.
Ling Yang. Zhilong Zhang, Yang Song, Shenda Hong, Runsheng Xu, Yue Zhao, Wentao Zhang, Bin Cui, and Ming-Hsuan Yang. Diffusion models: A comprehensive survey of methods and applications. ACM Computing Surveys, 56(4):1 2023.
Yichi Zhang, Yici Yan, Alex Schwing, and Zhizhen Zhao. Towards hierarchical rectified flow. The Thirteenth International Conference Learning Representations, 2025.
Qunxi Zhu and Wei Lin. Switched flow matching: Eliminating singularities via switching odes. In Forty-first International Conference on Machine Learning. 2024.
Qunxi Zhu, Yao Guo, and Wei Lin. Neural delay differential equations. In International Conference on Learning Representutions, 2020.
Qunxi Zhu, Yifei Shen, Dongsheng Li, and Wei Lin. Neural piecewise-constant delay differential equations. In Proceedings of the AAAI Conference Artificial Intelligence, volume 36, pp. 9242-9250, 2022.
14

---

Published as a conference paper at ICLR 2026
The structure of the appendix is follows:
Appendix
provides the connection between the probability flows of ODEs and DDEs.
Appendix
B provides the formal proofs for the theoretical results presented in the main text.
Appendix provides pseudocode for DFM(C) and DFM(D).
Appendix includes the experimential setup details and additional experimental results.
Appendix provides the comparison of the training and inference time required by DFM and CFM.
Appendix F provides the sensitivity anaysis of time delay parameter on both synthetic and real-world datasets.
Appendix G compares DFM with other ODE-based generative methods designed to address trajectory intersections and distributional heterogeneity, further highlighting the advantages of the DDE-based framework.
A CONNECTION BETWEEN THE PROBABILITY FLOWS OF ODES AND DDES
Based on the vector field it (t), x(t - r) in Eq. we can define the conditional average drift (CAD) without delay terms follows (Guillouzic et 1999):
x/b)) / dæ, [u(t, æ )p(,, 7/ +x
(24)
where p(æ,, T T|x, t; 2b) denotes the conditional probability given æ(t) p(æ, 1; 2., -|bb) = p( T lac, t; (b) t/ab), ) so Eq. (2) equivalent to:
Note that
0p(æ, V [ma(t, (/ (1b)))
(25)
Öt
which is precisely the continuity equation satisfied by the probability flows of ODEs. Remark A. This implies that a non-delayed vector field can constructed match the probability flow of the delayed vector field in Eq. [ However, the integration in Eq. (24) climinates the coupling information between the states at times & and т. Consequently, while the ODE preserves the marginal probability flow, it fails to maintain the coupling relationship during distribution transfer, as illustrated by a simple example Fig. (a).
B PROOFS OF THEORETICAL RESULTS
B. FORKER-PLANCE EQUATION OF DDES
Theorem B.1 (Fokker-Plance equation of DDEs, Guillouzic et al. 1999 Frank 2005
Consider
dime-dependent DDE with a single delay term:
dx(()) = ult, (x)) xx1 - -]] 1 [0,77]
dt
(26)
x<h)= = (h)) h [
where ult, (x), x0t at 00, T] Rd Rdd R2d is smooth vector field which is abbreviated as ult, x. xx) in the following, (h) represents the continuous initial function. The associated
probability flows p(z, t/pb IRd 00, + E satisfy the delay Fokker-Plamt equation:
0p(x, t at
{ { /dær [u(1, x, x (x),, 1; 1; (0r,,--(y66) -|/1]]
(27)
where p(æ, t|xb) is the probability density time t given the initial function, while pl , t; xr, t- T ab) the joint probability density representing the likelihood of the system being at æ at time and at x.
at time L
15

---

Published as a conference paper at ICLR 2026
Proof. The proof follows similar approach as presented in Guillouzic et al. (1999). Without loss of generality, we assume that and the state variable x E a. bj in Eq. (26). Consider an arbitrary Ca function F(x) defined on the interval (a, b), i.. F C ([a, 6]). which satisfies the following conditions:
lim F(x) lim F(x) = 0,
(28)
z-b
d F(z) O.
lim F() lim
(29)
dx
z-6 dæ
Then, by applying the Taylor expansion, we obtain:
dF(x(t)) = F(x (t) + dz(t)) - F[*(t)] - ultt att - -) d dr
Fla(0)) dt.
(30)
The ensemble average (average over realizations) of (1F(x(t)] can be written as
<
Flo0ol)) = ( (wl,, (1) dr dr
d
(=((()))
(31)
We denote p(a, 4; xr, 1 T/y dxdd, as the probability that sr(4) E [ar, + dæj and x(1 - r) E
[T, dr.,) given the initial function vb. Then, Eq. (31) is equivalent to
[da [ dx d dx Flo(()) [ dr,u(t, x. or )o(r) t,, T//).
doF(()) J." dær oe 1,,,
(32)
By applying the integration by parts formula to the right-hand side, we obtain:
f. deF((s) S[" dr. 2
per (,,,, - -l/b)
(33)
=f -('dn((t))
dxF((x) 1 dr, {
lu(t, {
x,)( 1.)
tz -|/6]]
Ox
where the surface terms are neglected due to Eq. leads to
28 and Eq.
29
Since F(x) is arbitrary, Eq. 33
0 or p(x,
0 -
Ör { dx, de (,,
(34)
-rr
B.2 PROOFS OF PROPOSITIONS THE MAIN TEXT
Proposition 3.1 (Restriction on trajectory intersections). Suppose that the target transport strat- egy. guided by certain key points or specific coupling constraints, inevitably induces trajectory intersections in the augmented phase space (t,z) during the transport process. the ODE-based FM is governed by a vector field that Lipschitz continuous in x and continuous in then the corresponding flow map cannot precisely preserve the transport strategy.
Proof of Proposition 3./ Consider time-dependent Neural ODE corresponding to ODE-based FM with the following form:
dm(()) x(t; e [0,7],
dt
(35)
x00) 2D.
where (tt:: : [0, Rdd - Rd is a parameterized vector field. Since the vector field is Lipschitz continuous in æ and continuous in f, by the Picard-Lindelöf Theorem, the corresponding initial value problem (35) has unique solution on the interval [0, 7]. Suppose that there exist two distinct solutions x (t) and (t) corresponding to different initial values æ and wz. which intersect at (t", x)). By the uniqueness of the solution, these two trajectories must correspond to the same solution, leading to contradiction. Therefore, given transport strategy, which inevitably leads to trajectory intersections during the transport process, the flow map corresponding to the ODE- based FM cannot precisely ensure the transport strategy.
16

---

Published as a conference paper at ICLR 2026
Proposition 3.3 (Heterogeneity in distributions). Assume that the source (resp., target) distribution
qu (resp.. 91) is supported on M (resp., N) disjoint compact sets {U]) (resp.. { NN 1j-1). where
each U. E IRd (resp., U] E R2) is path-connected set with non-zero measure and qo go( > 0 (resp.. ạı (111) > for all Zo E U U (resp., 21 E - If M < N, then for any transport U3 satisfying T#90- = qu, the following holds: if the flow map induced by a Neural ODE is equal to the transport map T, then the associated vector field cannot be Lipchitz continuous,
map : U 0 + أ
Proof of Proposition 3.3 Assume there exists parameterizee vector field vtt, a; 0) whose flow map Ør coincides with the transport map T. i.e., ΦT(x)) T(x) for all x E U UT.
First, consider the case M < N. Then there must exist m E (1, and n1, n2 E (1,2, ..., , N such that q: (T(U)) n Un > 0 and (1 (T(U%) n Ung > 0. Otherwise, for each m E {1,2, .M there exists a unique
[1,2, ... ,N such that q [T(U)) n nU = 0 for
all n' * n. Since M N, there exists n² E (1,2,-.. N) such that ΣM (T() n nU-)] = 0.
which implies
M U
M = UT(U))NU
9i (U -
T
D
nu
=0,
(36)
a contradiction.
Hence, we can choose Lm,1, 2m,2 E U such that 1 E U T(xm,2) E Un- By the path-connectedness. of U there exists continuous path 7 [0 1] - U with 70) = m, and
r(1) = Tm.2. Then, for any >> O, there exist Tg, Tg E 7 such that [z, - žel < E, and T(x,) and
T(..) lie in distinct subsets of {U]}z= Since each U3 is a compact subset of R2d, it follows that for any J1 E {1,2, .............................. ,N). J1 # j2.
n. . E U > 0. dist (U2, U12) = inf {1o- vl x E U2,
(37)
Let
d = min dist (U), U//o) jn:b2 >0.
(38)
f1
Then we have
(x.) - (π. [(x.)) T(xs)l 0. T.
(39)
[. x.
6
On the other hand, if v x;0) 0) is Lipschitz continuous with constant L 2 0, then for any x, y E Rd
wut,, x; 0) - vtt, y; 0; < L x - y). By Grönwall's incquality, the corresponding flow satisfies
(40)
(x) xb(3)) x y +
uis, x(s): 0] vis, y(s); dlds
ya
|0 [s, as):: uis. (s): 0|
(41)
< x yl + (z) - (z)lds.
Thus,
(T(x)) (T( S LT yl
(42)
holds for any x, y E Rdd. Now, choose < pir Then from 39 there exist By Fg E U such that
d
d - -ar, n. > xr. -Ralp
(pr(x.)) - (r >
(43)
which contradicts the Lipschitz bound.
Therefore, any parameterized vector field realizing such a transport map cannot be Lipschitz continu- ous. This completes the proof.
17

---

Published as a conference paper at ICLR 2026
Proposition 4.1. The marginal vector field u(t, x, x,,;,) given in Eq. (5), together with the lected initial function s, generates the probability path ((x)t and the joint probability path p(a, 2; T| in Eq. (4).
Proof of Proposition 4.1 Given that utt, æ, ær [2) is the conditional vector field that generates p(x. z|z, 4) and p( 1; xr, T|2, 2)), it evident that the following delay Bokkeer-Plancee equation holds:
0p(æ, 12, --V. { dx,/uu(t, æ, ar )p(z tt xr, T|z, r|2,,11] Ot
(44)
To verify this proposition, we only need to check that the marginal vector field u11, x, x.)) given by Eq. [15), the marginal probability density path p(æ, t/ab) and the joint probability density
p(a, t; Zz, - -|1b) satisfy the following delay fokker-Plance equation:
0p(æ, l = - dzz,
f [u(t, x. ),(, f; t; r|0)]
(45)
Üt
Assuming that all the functions involved satisfy the regularity conditions necessary for the interchange of integration and differentiation, we have that
Op(x, t|pb
Öt
0
/
p(æ, Uz b))(()dzz
01
/
De ple, 1/2 (2, +b) g(z)dz
/
attp 7|2, {b)] g(2)d=
{ der I /u(1,,,,=1 2))
V { derr
((=)dz p(a,, 1; F/abb) }
T/6b)
{ dee, [e(tt, we )p(2, 1; / [0>1]
(46)
Proposition 4.2. Given that p(æ, nu,t;a t,, T|1b) for all x, Rd and t [0, 1]. up to a constant independent of 0, LDcrM(() and (0rM(0) are equal, which further implies that VoLucru = (oLoru(
To guarantee the existence of all integrals and the validity of changing the order integration (as justified by Fubini's Theorem), we assume that p(æ, t|z, 4)), p(z, t; x,, r|2, 4b) decay to zero sufficiently fast as ||æ||| - and ||=,|| -> co, and further assume that (t, æz, xz, ()), vtt, x xz; 8), Vov(t, xz: are bounded. Note that (0, 1) and sp (xb) are both independent of 0, so fixed t and va in the following analysis. Using the bilinearity of the Euclidean norm and the independence of ull, æ, 22/22 from O, we have:
Proof of Proposition 4.2
V o Fo (trl)) (1, n, 2:; 0) ull,, 212
Vo -111))) (1 ær 0)]l 2 (v x. x, 0), utt, + ||u(t, (47)
V (t (1 (t, 0)|2 (v(t, x, 0), u(t,,,,,)))
and
V (eBg((), u |0,,
(v (v(t,
(48)
(t, x, x
00 x. æz; 1| 2 (w(1,, a, x,: ult x x2))))
18

---

Published as a conference paper at ICLR 2026
Next,
= /f (1(1, 1,,,, - daade
(49)
fiʃ 10(t, =|z, #))((=((((((()(()))))) (x))xzdade,
Finally,
E Bp (-1)) (v x. u),
=// v(t, xz x,: htt.t/t.......//.tt..t..tttt---- (2)
pla,ti -lod)
N0,tia,,t-11
- -// v(t, 2., xri ult,,,,z,_1))p(,tt:
))) (z)dz
dæda
a0,,, 0) ult, x. xr x-
7 (xx))dd
= -oo (wtt, -.; 0), u x, xr,
(50)
From Eq. 449 and Eq. (50). it follows that Eq. 47
is equal to Eq. 448 for any latent variable CDMM
which can be further deduced that VoLDCFM
Proposition 4.3 (Universal approximating capability of DDE's flow map). For any given and continuous transport map F Rd + Rd which push-forward the source distribution go the target distribution 9ı, i.e. F# q1 there exists neural network f(x; f) such that [1f (): [F(x) æ|I| < for all E Fd. then term v(t, x. xz; 0), where the corresponding flow map aa; 0), under the constant initial function condition, satisfies ||G((;; (x; F(x))ll 6. for all x E RR..
construct vector field with single delay
The proof is straightforward. We consider the following DDE with constant initial function and a parameterized vector field with a single delay term:
Proof of Proposition 4.3
dxc(() 01 0);
[0, 1),
(51)
x(t) = Co. t 0),
where the time delay 1, and the vector field depends solely on the delay term, independent of the current state and time, which is special degenerate case of DFM. In this case, for any point sampled
from the initial distribution ao 4o- the corresponding vector field v|ar(1 - -); 0) = v ( 0) remains constant over the time interval /. [0, 1), which means that the flow map G Rdd - R
associated with the above DDE will map Do to 2b + v[z0; 0 1. For any given e > suppose
there exists a neural network f(z; 0) such that [ 0) [F(x) x]]| < e, for all E _d Let
v vxx;; () = f (xo; 0). Then it follows that G(æn; 0) F(x ||| = Do + < 0) 1 - F'(x))] s
[ıo; 8) f( I)ll [= 8) [F(x.) æull c. This implies that the flow map of the
Neural DDEs approximately peserves the target transport strategy, and naturally, can approximately
push-forward go qı.
C
ALGORITHMS
To more clearly illustrate the overall workflow of our proposed DFM, we provide pseudocode for different DFM methods: Algorithm forr DFM(C) and Algorithm for DFM(D).
D EXPERIMENTAL SETUP DETAILS AND ADDITIONAL RESULTS
In this section, provide detailed explanation of the experimental settings for different datasets described in the main text and present additional experimental results. The experimental details for the delay dynamical systems and single-cell datasets summarized in Table 5).
19

---

Published as a conference paper at ICLR 2026
Algorithm 1 Delay Flow Katching with Constant Initial Functions
Input: ( the source data X go and the target data X1 rig Time delay 7
Output: Model 0(t, 2,2;;, 0) for the vector field with a single delay term 2, ===((-17
Initialize Model: 2(t,2,2,;0) # Arbitrary neural network structure
Select Coupling Strategy: KPG-OT/OT/Indeeende
Select interpolation Method: Lánear/Geoéesic//S//
# Consider
the linear interpolation below
Initialize Initial Function: (t; = 2o for t [-
While not converge:
Sample 2p Co. C X1 # Generate batch data
aetimizer.zero.gran
#0,21 = Coupling Soupling.Stratery (#n, œ1) Construct coupling
torch.rand (batchsize)
Randomly sample le [0,1] Sample state at I based on the linear
(+() (1-)=11-4))
2 +1x
interpolation
=( (1- =(t-r)=(l-tt+)x<m 2o + (t-r) x for te [T, i Sample delay state
at t- based on the linear interpolation
att zn for t (0, # Sample delay state at t-r based the constant initial function
Loss Loss. bass.backward() Optimizer.step return Model
J.pow (2) mean ()
Table 5: The experimental setup details for delay dynamical systems and single-cell datasets.
Sctup
Autoregulation
Spiral DDE
Mouse hematopoiesis
iPSC's
Data
Dimension
2
Hidden layer
64
Structure
Activation
Tanh
Tanh
SELU
SELU
Time input
False
False
True
Tree
Delay tern inpur
Thee (DFM) / Palse (CIM)
Thue (DPM) / Palse (CPM)
True (DFM) Palse (CPM)
Thue (DFM) Palse (CFM)
Coupling
KPG-OT
OT
OT
OT
Initical functiio (DFM)
Lateni variable
CSpline
CSpline
Goodesic
CSpline
Training
Batch size
256
256
128
128
lierationn
2h
2k
102:
26
Optimiser
Adam
Adam
Adam
Adam
Learning
Results
Figures
Figs
Figs
Fig
Fig
D.1 BIOLOGICAL AUTOREGULATION MOTIF
Dataset generation. The autoregulation motif, which is one of the most common biological network motifs, can be modeled the following 1-d DDE:
11
a(t) + x T)
(52)
where 7
As shown in Fig. 3 the dynamics exhibit damped oscillations. We sample 1, 000 initial points from the uniform distribution i/.22 1.2) and perform forward integration based on constant initial functions 26 (b;; æx) æo, h E [-1, Oj. to obtain trajectories. After that, select snapshots at 0, as the training dataset.
2, 2
Details for training. For each training iteration, we randomly and independently sample 256 data points from each snapshot. Pairing between minibatch data points from adjacent time steps is performed using KPG-OT. For this task, utilize the true coupling of only two keypoint pairs between adjacent time steps as guidance, although incorporating more keypoint pairs would enable more accurate pairing. Once data points across time steps are fully matched, construct 256
20

---

Published as a conference paper at ICLR 2026
Algorithm 2 Delay Flow Katching with Diverse Initial Functions
Input: ( the source data X 9o and the target data X1 q1 Time delay 7
Output: Model 0(t, (1,,; 0) for the vector field with a single delay term 2, ===(t-1)
Initialize Model: 2(t,2,2z;0) # Arbitrary neural network structure
Select Coupling Strategy: KPG-OT/OT/Indepedee
Select Interpolate interpolation Method: Láhear/Geodesic//55g # Consider the linear interpolation below
Select Clustering Method: Gaussian Mixture Mode1/DBSCAN
Cluster the source and target dataset:
x : x(M) ClusteringgMethhod...
X 11)
.... X
(N) ClusteringgMeethodd. (Xi)
Initialize Diverse Initial Functions for Different Clusters:
Select different constants Cmnr for zu E eX"",, 21 x
())(())))) = Zo + Cmn x for LE [-T,0
While le not converge:
Sample 2o C Xo. æ1 C C # Generate batch data Detimizer.zzro.gga.
2; = Coupling
(æa, σ1) Construct coupling Identify the cluster membership of zu: Imo Identify the cluster membership of æ1: Im corch.rand (batchsize)
Randomly sample LE [0,1] Sample state at I based on the linear
a(()= (1-)=111))) #1 interpolation
(
2o + for te [T. Sample delay state
at t- based on the linear interpolation
Chay x (t for te (0, # Sample delay state at axt a(t-r)=
t t- based on the diverse initial functions
Loss =(1-) (01-20) pow (2) mean ()
Loss boss.backward()
Optimizer.step return Model
transport trajectories using cubic spline (CSpline) interpolation. These trajectories provide the states, delayed states, and vector fields at various time steps, enabling the computation of the training objective Eq. (16). which is then optimized via backpropagation. Further details on the experimental setup can be found in Table
Details for testing. After training, we evaluate the performance of KP-DFM(C) and KP-CFM on both interpolation and extrapolation tasks. Specifically, randomly select 100 initial points from the initial distribution 2(..2, 1.2) and generate predicted trajectories by performing forward integration on the learned vector field, obtaining the predicted trajectories and snapshots at various time steps. For interpolation, we compute the mean 2- Wassersstein distance W2 between the predicted and ground- 0.5, 1.5. calculate the enror between the predicted and ground-truth trajectories at these points, quantified by the Mean Squared Error (MSE) and Mean Relative Error (MRF). For extrapolation, we similarly compute the W2 distance between the predicted and ground-truth distributions at 8, 9, 10, along with the trajectory- wise error using MSE and MRE.
truth distributions
6.5. Additionally, since the true dynamics are known, we
Results. As shown in Table and g KP-DFM(C) significantly outperforms CFM in both interpolation and extrapolation tasks, achieving much lower prediction errors for both distributions and individual trajectories, demonstrating superior representational capability. Furthermore, illustrated in Table o CFM fails to reconstruct the underlying dynamics, whereas KP-DFM(C) successfully captures the damped oscillatory dynamics, when the time delay deviates from the ground truth.
21

---

Published as a conference paper at ICLR 2026
Table 6: Additional results biological autoregulation motif based CFM with keypoint-guided coupling (KP-CFM), as well as DFM with constant initial functions and keypoint-guis coupling (KP-DFM(C)). For DFM, the generation results illustrated with various time delays, specifically 0.6, 0.8, 1.0 (ground truth), 1.2, 1.4.
Method
KP-CFM
KP-DEMIC) (r = 0.6)
KP-DFMIC) (r = 0.8)
Results
Method
KP-DFMIC) (r 1.0)
KP-DFMIC) (r = 1.2)
KP-DFMIC) (r = 1.4)
Results
Table Comparison of the interpolation prediction error between KP-CFM and KP-DFM(C) with different time delay on the biological autoregulation dataset.
((CFM)
0.6
0.8
1.0
1.2
1.4
W2
0.1842
2.82e-2
1.55e-2
1.53e-2
2.27e-2
4.37e-2
MSE
9.21e-2
1.51e-3
4.42e-4
3.45e-4
5.36e-4
2.68e-3
MRE
0.117
2.22e-2
1.23e-2
9.47e-3
1.14e-2
2.57e-2
D.2 SPIRAL DDE
Dataset generation. We consider the following 2-d DDE:
((() = Atanh ( +
[ (+(
(53)
where 7 0.5 and A = [[-,,000, [-30, -1] E 2 As shown in Fig. the dynamics exhibit
trajectory crossings. We sample 1, 000 initial points from the uniform distribution L/(0., 1) x ((0.0.1) and perform forward integration based on constant initial functions (h; xo) = 2p. E
[-0.5, 0j. to obtain trajectories. After that, snapshots intervals of 0.05 within & E [0, 0.6]
selected as training dataset.
Details for training. In each training iteration, 256 data points randomly and independently sampled from each snapshot. Pairing between the minibatch data points at adjacent time steps is performed using OT. After matching data points across time steps, 256 transport trajectories are
Table 8: Comparison of the extrapolation prediction error between KP-CFM and KP-DFM with different time delay on the biological autoregulation dataset.
((CFM)
0.6
0.8
1.0
1.2
1.4
W2
2.13e-2
2.25e-2
8.67e-3
3,42e-3
4.55e-3
1.23e-2
MSE
5.51e-4
6.31e-4
1.28e
1.46e-5
4.38e-5
2. 2.13e-4
MRE
1.59e-2
1.68e-2
6.75e-3
2.26e-3
3.72e-3
9.37e-3
22

---

Published as a conference paper at ICLR 2026
constructed using CSpline interpolation. These trajectories provide the states, delayed states, and vector fields at various time steps, enabling the computation of the training objective in Eq. (16). which is subsequently optimized via backpropagation. Additional details on the experimental setup are provided in Table
Details for testing. After training, we assess the performance of OT-DFM(C) and OT-CFM on the interpolation task. Specifically, 100 initial points are randomly sampled from the initial distribution 4(0.8, 1) (4, 0.1). and predicted trajectories are generated via forward integration of the learned vector field, yielding trajectories and snapshots at various time steps. The mean W₂ distance between the predicted and ground-truti distributions is computed at intervals of 0.01 within [0, 0.6]. Moreover, we calculate the eπror between the predicted and ground-true -ruth trajectories at these points, quantified by MSE and MRE.
Results. As shown in Table TO, OT-DFM(C) significantly outperforms CFM in the interpolation prediction task, achieving substantially lower errors for prediction of both distributions and individual trajectories. Additionally, as illustrated in Table g CFM collapses around the trajectory crossing area, while OT-DFM(C) successfully captures the spiral dynamics, even when the time delay т deviates from the ground truth.
Table 9: Additional results on Spiral DDE based on CFM with OT coupling (OT-CFM), as well as DFM with constant initial functions and OT coupling (OT-DFM(C)). For DFM, the generation results are illustrated with various time delays, specifically T = 0.35, 0.40, 0.45, 0.50 (ground truth), 0.55.
Method
OT-CFM
OT-DFM(C) (r 0.35)
OT-DFM(C) (r = 0.40)
Podido
Parskin
Results
Gomnd
-2.5
Method
OT-DFM(C) (T D.45)
OT-DFM(C) (T D.50)
OT-DFM(C) (T 0.55)
Daliske Orendown
Palido Oweddww
Fadiiies Dranciaa
Results
Table 10: Comparison of the interpolation prediction error between OT-CFM and OT-DFM with different time delay on the spiral DDE dataset.
((CFM)
0.35
0.40
0.45
0.50
0.55
W
3.95
0.33
0.43
0.42
0.20
0.45
MSE
15.81
0.07
0.13
0.12
0.03
0.14
MRE
13.14
1.25
1.66
1.80
0.68
1.96
D.3 MOUSE HEMATOPOIESIS DATASET
Data preprocessing. To efficiently apply various trajectory inference methods to scRNA-sec data from mouse hematopoiesis, we first project the data into a low-dimensional space. Specifically, use the reduced two force-directed layouts (SPRING) space for this dataset with batch correction. following the same preprocessing pipeline TIGON (Sha et al. 2024). Besides, clustering algorithm to categorize all cells at Day into two distinct clusters, corresponding to the Neu fate and Mo fate, respectively.
the DBSCAN
23

---

Published as a conference paper at ICLR 2026
Selection of initial functions. During training, we match the sampled minibatch cells at Day 2 and 6 using OT. For each matched cell pairs (æo, Ø1). where Dy is sampled from data of Day 2 and æ1 is sampled from data of Day 6, the initial function is defined as:
deb" (t; æu) dt
[[1, 0), if æ is in Neu cluster, [0, 1], if æ1 is in Mo cluster.
(54)
Details of training. During training, we align Day 2 and Day with == 0 and = 1, respectively. For testing, the distribution of Day is predicted at t 0.5. In this task, the time delay is set In each training iteration, we randomly and independently sample 128 cells from each snapshot at Day 2 and Pairing between minibatch data points adjacent time steps is performed using OT. After matching data points across time steps, 128 transport trajectories are constructed using geodesic interpolation (refer to the next paragraph). These trajectories provide the states, delayed states, and vector fields at various time steps, enabling the computation of the training objective in Eq. (16), which is subsequently optimized via backpropagation. Additional details on the experimental setup provided in Table
1.
Construct geodesic interpolation yy. For each matched cell pairs, geodesic interpolation is employed to generate the interpolation trajectory (Kapušniak et al. 2024). Here, we first build a data-induced Riemannian metric on the data manifold M. Formally, Riemannian metric g on M is a smooth family of inner products the tangent spaces of M. Specifically, 9 provides each .M a positive defined symmetric bilinear form on T M,
Ix T&M X T-M-R
(55)
which induces a norm 11.1 T2M - R defined by ||o|la. V Im (v, Based on this norm, we can calculate the length of any path connecting two points on the manifold and construct the geodesic between them. We aim to equip the ambient space I' with an appropriate metric that constrains geodesics remain close to the data manifold. Specifically, given coordinates, the Riemannian metric g can be equivalently represented as a state dependenn d x d symmetric positive definite matrix (zx), such that Im (n( v) = u (z)) y for any u, E T.M. Intuitively, impose manifold constraints on the transport path, || |G())l is supposed to take smaller values around data points and larger values in regions far from all points in the dataset D. Therefore, we employ the LAND metric,
where Gæ; D) (diag (aong(((x))) (z: D) + eI)- -1- Here, the components of h(æ; D) are defined
Σ(nf 2x) Pxx (- 20 I (-20²
hg (x; D)
x²
(56)
where 1 15 k S d denotes the k-th component, the superscript i represents the i-th data point in the training dataset 'D and a is the kernel bandwidth. Subsequently, we utilize neural network o, to parameterize the geodesic interpolation path between 2g and æ1 as
Yn wo, 1 (1 t)20 + tæ + f(1 t) 3o, i1).
(57)
The trajectory Yy approximates the geodesic between any two points sampled at adjacent time steps on the Ricmannian manifold (M. g) implies learning parameters 7) to minimize the following loss function:
J (t,tt, Do, x1) G(m,:: (t,, Tu. x dt,
Ln(n) == -E((
(58)
where π represents the OT coupling. The optimal parameter obtained as:
1 arg min (n),
(59)
By substituting 17 into Y, in Eq. (57). the resulting trajectory 7n provides an approximation of the geodesic interpolation on the data manifold.
D.4 IPSCs DATASET
Data Preprocessing. To apply trajectory inference methods to the iPSC scRNA-seq data efficiently, first reduce the data to 4-dimensional space using Principal Component Analysis (PCA),
24

---

Published as a conference paper at ICLR 2026
following the preprocessing pipeline of TIGON (Sha et al. 2024). For cells after differentiation (from Day 3), we categorize them into two groups corresponding to the M and En fates using a Gaussian Mixture Model (GMM) with two components. For Day 2. we model the distribution with single-component GMM and denote its mean as Ma- As descrihed above, for Day 5, the distribution is modeled with two-component GMM, where the means of the M and En fates denoted J411 and 112, respectively,
Selection of initial functions. During training, we match the sampled minibatch cells at Day 2, and 5 using OT. For each matched triplet (æo. 11 Ø2), where 2o is sampled from the Day 2 dataset, 21 from the Day dataset, and ø2 from the Day 5 dataset, the initial function is defined as:
dep (t; 2o)) di
|1n 140, if 22 is in M component,
(60)
(112 po, if Ø2 is En component.
Details of training. During training, we align Day 2, 4 and 5 with & = 0. = 2 and t 3. respectively. For testing, the distribution of Day is predicted at t = 1. In this task, the time delay is set to 3. In each training iteration, we randomly and independently sample 128 cells from each snapshot at Day 2, 4, and 5. Pairing between minibatch data points at adjacent time steps is performed using OT. After matching data points across time steps, 128 transport trajectories 7 constructed using CSpline interpolation. These trajectories provide the states, delayed states, and vector fields at various time steps, enabling the computation of the training objective in Eq. (I6). which is subsequently optimized via backpropagation. Additional details on the experimental setup are provided in Table
D.5 MNIST DATASET
Dataset Generation. We propose Semi-paired Image-to-Image Translation task using the MNIST dataset. The source domain consists of the original MNIST images, while the target domain comprises their corresponding negative images. Specifically, the pixel values X of each image in the original MNIST dataset are inverted by replacing them with - X, which flips the brightness. The objective to map each image to its corresponding negative counterpart, as demonstrated in Table []]
Minibatch coupling. To provide partial supervision, 10% of the training data are paired with their negative counterparts keypoints. During training, minibatches are independently sampled from the source and target distributions, and KPG-OT coupling (KP-) is applied. Specifically, for each batch, we sample the same number of images from the source distribution and the target distribution. Out of these, around 10% image pairs are selected from the keypoint set, while the remaining images are independently sampled from the respective distributions. It important to note that the negative counterparts of the remaining images sampled from the initial distribution are not guaranteed to appear in the target batch. Subsequently, we perform KPG-OT OT pairing on the remaining images in the batch based on the selected keypoint pairs, which are then used for further training.
D.6 CIFAR-10 DATASET
Source distribution. The source distribution for zu is defined as follows: with a probability of 50%.
2p is sampled Iu corch.randn(33 32. /2 0.5. and with the remaining 50% probability, it is sampled as zo = cor 32, 32) + 0.5, following the setting in Zhu & Lin (0024). A sample drawn from the source distribution is shown in Table T2
Trainable initial functions. Designing appropriate initial values for high-dimensional images in image generation tasks is a particularly challenging problem. To address this, employ torch.nn. .ech.nn..mbbdding. to automatically learn the initial functions. Specifically, the source distribution consists of two modules, corresponding to two Gaussian distributions, while the target distribution consists of ten modules, corresponding to ten image categories. Therefore, for each pair of modules (m, where m E {0, 1} and n E {0, 1.. 9), we need to design the time derivative Cmn of the initial function (Eq. (22) ). To achieve this we begin by defining Emb=torch.nn wtw.rroccc.....medinennngggg Embeddinggnnn 3 x 32 32), which maps 20 discrete indices to corresponding tensors of the same shape as the image. Then, for each pair (m, rt), we define the initial function as:
dgb,n. (t; 2o)
Emb(10m + n), 3tnn (0; wo) = Zu.
(61)
dt
25

---

Published as a conference paper at ICLR 2026
Table 1:: Samples of generation process and results MNIST data based CFM with keypoint- guided coupling (KP-CFM), as well as DFM with constant initial functions and prrpoint-guider coupling (KP-DFM(C)). For DFM, the generation results are illustrated with various time delays, specifically 0. 0.125, 0.25, 0.5, and 1.0. In each image, the first and last columns represent samples obtained from the souree data and their negative counterparts, respectively, while the 10 intermediate columns depict the generation process.
Method
KP-CFM
KP-DFM(C) (T 0.125)
KP-DFM(C) (r 0.25)
0
000 0 0
000
0000
7
2
2
Generation
9
9
9
4
5
5
5
5 3
6
6
666
6
Method
KPACFM
KP-DEM 0.5)
KP-DFM(C) (r 1.0)
0
o
0
0
1
2
2
Generation
4
6
where E [-r,0], This approach allows the derivative of the initial functions across different modules to adaptively adjust during training, thereby enabling more flexible and efficient learning.
Experimental setup details. For the image generation tasks, both CFM and DFM are trained following the configurations outlined in Tong et al. (2023a]b) In particular, we utilize UNet architecture with the following structures and training parameters:
channels 128,
depth 2,
channels multiple = 2,
heads 4,
heads channels 64,
attention resolution 16,
dropout 0.1,
batch size per gpu 128, gpus = 1,
Adam optimizer with 31 0.9, 32 = 0.999, e = 10-8, and no weight decay,
learning rate = 2 x 10-4
gradient clipping with norm 1.0,
exponential moving average weights with decay 0.9999.
26

---

Published as a conference paper at ICLR 2026
Table 12: Samples from the source distribution and generation results based CFM and DFM(D) using independent coupling (I-) or OT coupling (OT-) for adaptive NFE on transporting the Gaussian mixture model with components the CIFAR-10 dataset.
IDDFM(D)
OT-DFM(D)
Sourc
Table 13: Samples from the source distribution and generation results based CFM and DFM(D) using independent coupling (I-) or OT coupling (OT-) for adaptive NFE on transporting the Gaussian mixture model with components the CIFAR- 10 dataset.
NFE-30
NPE-ll
LCRM
LDFM(D)
OT-DEMIT)
E TRAINING AND INFERENCE TIME OF DFM
In this section, we compare the training and inference times DFM and CFM. Through detailed analysis and illustrative examples, show that the training cost of DFM is only marginally higher
27

---

Published as a conference paper at ICLR 2026
than that of CFM. In terms of inference efficiency, when facing trajectory intersections or distributional heterogeneity, DFM achieves substantially faster inference because it produces transport paths that are effectively "straighter" than those of CFM.
Training time. Similar to CFM, we formulate a regression objective to efficiently train the vector field of DFM (see Eq. (T6)). The key difference is that DFM additionally requires specifying the initial functions and sampling delayed states. Note that the initial functions are selected and fixed before training and therefore introduce no additional computational cost. The only extra cost arises from delay sampling: in CFM, given sampled initial and target points, we construct an interpolation path and sample the state x at time t. In DFM, after constructing the same interpolation path, we must sample both the state x at time and the state 2, at time t - т. This additional sampling makes DFM slightly more expensive to train than CFM, though it still remains highly efficient in practice. On the 2-d Gaussian mixture dataset, we record the per-iteration training time of both models. Table [4] reports the mean training time per iteration together with the corresponding standard deviation (over 1. 1,000 iterations). DFM requires a slightly longer training time than CFM, but the difference is very small (DFM: 2.014+0. 109 s/ierr vs CFM: 1.986 sfiter).
Inference time. Compared to CFM, DFM possesses stronger representational capacity, enabling it to more accurately approximate the target transport map and the associated trajectories. As result, during inference, DFM typically follows much "straighter" transport path, leading to significantly improved inference efficiency. This advantage becomes particularly pronounced in scenarios involving trajectory intersections or distributional heterogeneity. We illustrate this difference using two examples. As shown in Fig. [5] on the 2-d Gaussian mixture dataset, the keypoint- -guided optimal transport induces trajectory intersections along the transport process. DFM can directly approximate the linear transport paths and therefore requires only a very small number of Euler discretization steps to achieve accurate generation. In contrast, CFM inevitably induces substantial trajectory bending near the crossing regions, necessitating much more Euler steps to obtain reasonable results. To quantify the inference efficiency, we employ an ODE solver with adaptive step size to compute the generation trajectories of both CFM and DFM, As reported in Table 14] DFM requires only 0.0091+#0000 0062 $ of average inference time per sample, substantially lower than the 0. 010053 required by CFM. The 1-d example in Fig. 6 Gillusttrtes the inference behavior of DFM and CFM under distributional heterogeneity. Near the singularity, CFM exhibits pronounced trajectory bending, which forces the model to take more discretization steps during inference. In contrast, DFM produces substantially straighter trajectories, achieving accurate generation with essentially a single step. Consequently, DFM attains higher inference efficiency along with improved accuracy in this setting.
Curvature of the transport paths. To assess the straightness of transport trajectories produced by DFM and CFM, we generate 1,000 transport paths the 2-d Gaussian mixture dataset using the trained DFM and CFM, and compute the absolute curvature at the midpoint of each trajectory. The results, reported in Table [4] show that DFM yields trajectories with consistently small curvature (0.010220.0133). indicating substantially straighter paths, whereas CFM produces trajectories with much larger curvature (4.3303 13.9777), reflecting pronounced bending. This observation further supports the higher inference efficiency of DFM compared with CFM.
we report generation performance of CFM and DFM(C) on the 2-d Gaussian mixture dataset under different numbers of Euler discretization steps. We evaluate generation quality using the W2 distance between the generated and target distributions. Moreover, because the keypoint-guided coupling strategy induces optimal transport map (X - -X), we additionally compute the mean squared error (MSE) between the generated target points and their ground- -ruth counterparts. The results show that DFM(C) consistently outperforms CFM and achieves high-quality generation with only a small number of steps. On the CIFAR- 10 image generation task, we further verify that DFM achieves accurate generation with small NFE (see Table 13). and that its generation quality at low NFE is significantly higher than that of CFM (see Table
Generation quality across different steps. In Table 14
Additional remarks. It is important to emphasize that our DFM framework is not designed solely for image generation, where faster inference and straighter transport paths are clear advantages. DFM is also motivated by a broader range of scientific problems- -ucch as inferring single-cell differentiation trajectories from snapshot data or recovering delay-driven dynamical systems where overly straight trajectories are not only unhelpful but can even be detrimental. Consider single-cell trajectory inference as an example. A fundamental requirement in this setting is that the inferred trajectories
28

---

Published as a conference paper at ICLR 2026
supposed to remain the data manifold, ensuring that the recovered differentiation paths reflect biologically plausible intermediate states. Excessively straight transport trajectories tend to leave the manifold, producing intermediate cell states that lack biological meaning and leading large inference errors. As illustrated Fig. 4 although MIOFlow (Fig. (b)) infers straighter trajectories than OT--FO (Fig. A (f), its paths deviate substantially from the gene-expression manifold, resulting in significantly lower trajectory inference accuracy (see Table 2
Table 14:
Comparison of CFM and DFM(C) the 2-d Gaussian mixture dataset in terms of training time, inference time, cunsport-trajector curvature, and generation quality across different discretization steps.
Method
CFM
DFM
Training time Inference time
1.986+0.11 s/iter
2.014=0.109 s/iter
0.0442.10.053 s/sample
0.0091 10.0062 s/sample 0.0102:100033
Curvature
1..903
1-Step W2
0. 0.15
0. 1569
1-Step MSE
0.0820
0.0211
2-Step W2
0.5107
0.1325
2-Step MSE
0.4030
0.0127
4- 4-Step W2
0.4531
0.0830
4-Step MSE
0.4123
0.0057
Adaptive-Stee W2
0.2884
0.0494
Adaptive-Step MSE
1.0729
0.0024
Initial distribution
Target distribution
Generuted distribution
Keypoint Puirs
XXX
(b)
X
Step
Steps
4 Steps
Adaptive Steps
Figure 5: Comparison of CFM (a) and DFM (b) on the Gaussian mixture dataset under different numbers of sampling steps. Since DFM allows trajectory crossings, its generated paths are straighter, enabling accurate generation with very few sampling steps and substantially lower inference time. In contrast, CFM exhibits pronounced bending near trajectory-crossing regions, which requires more sampling steps, leads to slower inference, and results in lower-quality generation.
29

---

Published as a conference paper at ICLR 2026
Initial distribution
Targel distribution
Generated distribution
Step
Steps
Steps
Adapiive Steps
Figure 6: Comparison CFM (a) and DFM (b) on the I-d synthetic dataset with distributional heterogeneity under different numbers of sampling steps. CFM exhibits singularity near O, resulting in inaccurate generation and pronounced trajectory bending in the local neighborhood around the singularity. This necessitates more sampling steps and slows down inference. In contrast, DFM, by employing diverse initial functions, produces straighter trajectories for the two subsets from the source to the target distribution, achieving accurate generation with very few sampling steps.
F ABLATION STUDY
In this section, we conduct systematic ablation study. Recall that, relative the FM, our DFM framework additionally incorporates the delay term and the initial functions. We conduct sensitivity analyses and discuss selection strategies for both components to demonstrate the robustness of method their choice. Moreover, our framework admits different path interpolation schemes, coupling strategies, and the clustering procedure used to determine initial functions, for which we also provide detailed examination.
F.1 SENSITIVITY ANALYSIS AND SELECTION STRATEGY FOR THE DELAY
To evaluate the effect of time delay 7 selection on DFM performance, we perform sensitivity analysis using both synthetic and real-world datasets.
First, we examine how different values of impact the reconstruction of delay dynamical systems. For the biological autoregulation motif, where the true time delay experiments with values ranging from 0.6 1.4. As shown in Table 7 g while the prediction errors for both the distribution and trajectories increase as the delay deviates from the true value, KP-DFM(C) consistently captures the damped oscillation pattern across wide range of delay values. In contrast, CFM fails to capture this behavior accurately. For the spiral DDE, with a true delay of 0.5, we experiment with delays from 0.35 0.55. Although the prediction errors become larger when the delay does not match the true value, OT-DFM(C) consistently outperforms CFM, achieving significantly better results, as shown in Table TO Additionally, illustrated in Table 9 OT-DFM(C) successfully captures the spiral dynamics across different delay values, whereas CFM exhibits divergence the intersections the trajectories.
1, we conduct extensive
Furthermore, we explore the effect of different time delay values on image generation performance. Using the MNIST dataset, we compare the generated images for 0.125, 0.25, 0.5, and 1.0. As shown in Table II] the generated images from DFM significantly outperform those generated by CFM, with FID scores also notably lower than CFM's, as summarized Table 3 Upon examining the impact of varying time delays, we find that performance is least favorable can be attributed to DFM's behavior approaching that of CFM as 7 tends towards 0. However, for values of 0.25, the generation quality improves substantially.
0.125, which
For the choice of the delay term, we consider two scenarios: (1) the underlying system has no inherent delay; (2) the underlying system is itself delay dynamical system.
the system has no explicit delay structure- for example, in image generation tasks the mere introduction of delay term is sufficient to mitigate issues such as trajectory intersections, thereby improving generation quality and yielding results that significantly outperform ODE-based FM models. The sensitivity analysis of the delay term the MNIST dataset further supports this claim
30

---

Published as a conference paper at ICLR 2026
In such cases, the key factor is simply the presence of delay term, whereas the exact numerical value of T plays a relatively minor role. If an optimal value of T is desired, standard hyperparameter selection techniques such grid search can be employed.
(see Table 3)
If the underlying data correspond to a genuinely delay system, then the choice of becomes critical for accurately capturing the system's 's dynamics,
illustrated in the biological autoregulation example (see Table [ 7 and 8) Here, the closer T is to its true value, the more accurate the trajectory reconstruction. For certain biophysical systems, prior domain knowledge often provides approximate
estimates of (an der Heiden 1979; Glass et al. 2021; Rihan et al. 202T). When such prior
information is unavailable, we can resort to delay identification methods from dynamical systems theory to infer the delay structure directly from data (Ma et al.) 22000] 2017).
F.2 SENSITIVITY ANALYSIS AND SELECTION STRATEGY FOR THE INITIAL FUNCTIONS
For generative frameworks based on DDEs, the initial function is an essential component. In the main text, we construct different forms of initial functions for different tasks, such as constant initial functions and diverse initial functions. For the diverse initial function, we further specify different constant time derivatives depending on the scenario (see Section 4.3.2 Eq. (22 In this section, demonstrate that DFM exhibits strong robustness with respect to both the choice of form and the selection of constant derivatives.
We first examine how to select the form of the initial functions. Note that diverse initial-function designs introduced to better address transport across heterogeneous distributions. Distinet initial functions guide samples toward different target subsets. Therefore, if the source and target distributions are known decompose into multiple modules (e.g.. when single cells differentiate into several cell types), it is preferable to employ the DFM with diverse initial functions to model the generative process. If no distributional heterogeneity is present, a simple constant initial function suffices. In fact, we have theoretically established that the DFM with constant initial functions already possesses universal approximation capability for continuous transport maps between arbitrary distributions. Moreover, experiments on synthetic delay-dynann dynamiccal data and the MNIST dataset further confirm that a constant initial function enables DFM to model transport maps with trajectory intersections and to capture the underlying delay-driven dynamical behavior.
Next, we examine how the selection of constant derivatives affects the performance of DFM when diverse initial functions are used. We conduct series of comparative experiments on the real iPSCs single-cell trajectory inference task. In the main text, the initial functions employed correspond to Eq. (60), which are specifically designed based the data structure. Here, we instead adopt randomized procedure to generate distinct initial functions for different cell types, as follows:
dub* (1; wo) dt
Vec1,, if 22 in M component, Rendom if 302 is in En component.
(62)
where Random independently from the uniform distribution. We randomly sampled 10 sets of initial functions, trained separate DFM for each, and conducted the single-cell trajectory inference task. The mean and standard deviation of the inference errors are reported in Table 15 (corresponding to DFM(DR)). The results show that, compared with the specially designed initial functions (Eq. (60 ), using randomized initial functions yields only a slight increase in the inference error, while still outperforming the DFM with constant initial functions and substantially surpassing the performance of CFM. This demonstrates the strong robustness of the DFM framework to the choice of initial functions. Moreover, the comparison between DFM(D) and DFM(C) further indicates that, in the presence of distributional heterogeneity, adopting diverse initial functions leads to improved generative performance, thereby reinforcing the selection strategy discussed above.
and Random veabbommee--rrhhrraadddd/ are 4-dimensional vectors drawn
F.3 ABLATION STUDY ON THE QUALITY OF CLUSTERING
For DFM with diverse initial functions, the design of these functions depends on the clustering results of the source and target samples. We examine two scenarios in which the clustering results are suboptimal and analyze how this affects model performance: (1) the number of clusters is smaller than the true number of underlying classes, and (2) the number of clusters exceeds the true number of classes.
31

---

Published as a conference paper at ICLR 2026
Table 15: Trajectory inference results the qPCR iPSC dataset with different initial functions. DFM(DR) denotes the DFM model equipped with diverse initial functions, where the corresponding constant time derivatives of the initial functions are randomly sampled.
Method
W2(L)
MMD(L)
W2(F)
MMD(F)
OT-CFM
0.579
0.492
0.226
0.030
OT-DFM(C)
0.553
0.447
0.234
0.041
OT-DFM(D)
0.532
0.399
0.213
0.027
OT-DFM(DR)
0.54510..00
1..411
0.218 10.008
0.028
When the number of clusters is smaller than the number of classes, multiple distinct subsets share the same initial function. In this case, the distributional heterogeneity among these subsets cannot be fully resolved. Nonetheless, the presence of the delay term mitigates the resulting singularities, and the overall generation quality remains substantially better than that of CFM. We illustrate this with 1-d example in which the task is transport U( ((1-,11) the two-component mixture - 22) + (2, 3). When the two modules are correctly identified and assigned distinct initial functions, the generated distribution (Fig. (c)) closely matches the target. However, the two modules are incorrectly grouped into single cluster, the resulting generation (Fig. 7 (b)) exhibits the singularity problem near in the source domain, causing a degradation in the generation quality. Nevertheless, the delay term still alleviates this issue, yielding results that remain superior to those of CFM.
When the number of clusters exceeds the true number of classes, single true class may be split into two subclusters. In this situation, the over-segmentation provides sufficient refinement to resolve the heterogeneity induced singularities. As shown in Fig. (d), further subdividing module of the target distribution does not negatively impact the final generation quality.
(a)
CFM
(b)
DFM (One2one)
(e)
DFM (One2two)
(d)
DFM (Onc2three)
12-0.111
W2-0.036
W2-0.013
W2-0.014
Figure 7: The generated results of CFM and DFM under different numbers of clusters, with the squared 2- Wasserstein distance between the generated and target distributions.
F.4 SELECTION OF OTHER COMPONENTS THE DFM FRAMEWORK
In this section, outline the rationale behind our choices of path-interpollaion methods and coupling strategies. Across different experiments, we adopt different interpolation schemes and coupling strategies, highlighting the flexibility of the DFM framework and its ability to accommodate variety of modeling choices. Both the selection of interpolation methods and the choice of coupling strategies are determined by the characteristics of each task and align with prior work related domains.
Specifically, for path interpolation, we use linear interpolation when no prior geometric constraints are present, yields the most direct transport path and improves generative efficiency (Liu et 2023: Esser et al., 2024). For tasks such as single-cell trajectory inference, where paths are required to remain on the underlying data manifold, we adopt geodesic interpolation, consistent with prior work (Huguet et al.) 2022 Kapušniak al. 2024). For continuous transport involving three or more distributions, paths connecting multiple intermediate points must be constructed, for which we employ cubic spline interpolation (Pooladian al. 2024
For the choice of coupling strategy, independent coupling or aptimal--ransppee coupling is typically used when prior constraints available (e.g., in image generation tasks). In contrast, for certain biological systems such as single-cell gene expression, many work has shown that their evolution tends follow energy-efficient transport patterns (Bunne al. 2024). In such cases, employing
32

---

Published as a conference paper at ICLR 2026
the optimal-transpor coupling strategy is the most appropriate choice and has been widely adopted in single-cell trajectory inference (Kapušniak et al. 2024 Sha et al.] 2024). Moreover, in many scenarios, partial keypoint correspondences are known prion, and the coupling strategy must respect these constraints. A line of research has proposed keypoint-guidee optimal transport (Gu et al. [022] 2023), which preserves keypoint correspondences while still minimizing the overall transport cost. Therefore, when partial prior coupling information is available, keypoint-guided optimal transport constitutes the most suitable strategy.
G COMPARISON OF DFM WITH OTHER ODE-BASED GENERATIVE METHODS
In Section we discuss a series of works that modify ODE-based generative frameworks to overcome FM's limitations in handling trajectory intersections and distributional heterogeneity. Here, we focus on two representative methods: Augmented Bridge Matching (AugBM) (De Bortoli et al. 2023) and Switched Flow Matching (SFM) (Zhu & Lin 2024). AugBM allows trajectory intersections during transport by b augmenting the velocity field with information from the initial point. SFM introduces an additional switching signal as conditioning variable and employs switching ODEs to enable transport between different modules of the source and target distribution, thereby alleviating singularities caused by distributional heterogencity.
As discussed Section 5 these ODE-based modifications typically address only one type of limitation and cannot simultancously resolve both. We illustrate this with two simple examples. In Fig, (a-c), we compare the generation results of CFM, DFM(C), and SFM trained under the keypoint-guide optimal transport strategy the 2-d Gaussian mixture dataset. While SFM can alleviate issues caused by distributional heterogeneity, its representation capacity is still weaker than DFM and cannot perfectly resolve trajectory intersections. Notably, although SFM allows trajectories cross between different mixture components, trajectories within the same component still cannot intersect. This is because points sampled from the same component share the same switching signal during generation, resulting in a single ODE that forbids intra- component trajectory crossing. Consequently, SFM cannot precisely realize the optimal transport mapping for this task (X + -X). In Fig. (d-f), we compare CFM, DFM(D), and AugBM on the 1-d example with distributional heterogeneity. Although AugBM allows trajectory intersections during transport, it still suffers from singularities when transporting between heterogeneous distributions, preventing it from achieving precise distribution transport.
We further highlight the distinctions between these ODE-based extensions and our proposed DFM framework. The ODE-based methods remain confined to the ODE generative paradigm, typically ad- dressing their limitations either by augmenting the vector field with additional conditioning variables or by designing highly specialized transport paths and procedures. In contrast, DFM departs from the ODE framework and leverages the vector field with a delay term to realize distributional transport. Notably, the inclusion of the delay term naturally overcomes the limitations inherent in ODE-based approaches within the original phase space, without requiring auxiliary conditioning variables or carefully engineered transport trajectories. Moreover, for transport tasks whose underlying dynamics involve intrinsic delays, such as the gene autoregulation system in Section ODE-based generative frameworks are incapable of recovering the true dynamics, whereas DFM, based on DDEs, can accurately model such delayed dynamics.
33

---

Published as a conference paper at ICLR 2026
Intiial distribution
Targel distribation Keypoint Pairs
Geneeaatdd distribution
(b)
(c)
CFM
DFM(C)
SFM
(d)
(c)
(
CFM
DFM())
AugBM
Figure 8: (a-c) Comparison of generation results from CFM (a), DFM(C) (b), and SFM (c) on the 2-d Gaussian mixture dataset. (d-f) Comparison of generation results from CFM(d), DFM(D) (e), and AugBM (f) on the dataset with distributional heterogeneity.
34


Reviews:
[{"review_id": "zLdgFx0c9l", "reviewer": "Reviewer_KnhQ", "summary": "The paper proposed to learn vector field that uses dynamics with memory (a dependence on a past state) to overcome two core limitations of standard flow matching: difficulty to handle trajectory crossings and difficulty with heterogeneous source/target supports. This “delay” preserves sample identity at intersections, enabling mappings like Gaussian point-wise negation and keypoint-consistent image translations that ODE flows struggle with. Both theoretical support and experimental validations are provided for the claims.", "questions": "1.\tCould the authors discuss realizations of initial functions and the path gamma for the two examples in Figure 1 and (intuitive) explanations of how the learnt vector field circumvent the described issues? \n\n2.\tCould the authors provide algorithms to summarize different DFM methods in the appendix ? \n\n3.\tWould it be more precise to say that the universal approximation theorem states the existence of delayed vector field that can approximate any mapping, rather than DFM can approximate any mapping? Because the flow matching objective for identifying those vector fields might still need to be modified (with different choice of coupling and gamma) for different cases ?", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "1.\tThe proposed approach is novel and and enables learning more flexible vector fields using flow matching than is possible with existing approaches.\n2.\tExperimental Observations are also interesting and clearly show the benefits of DFM.", "weaknesses": "1.\tThe writing could be made clearer. The paper is dense in terms of information content. A simple and concrete running example to illustrate each of the section’s content would be helpful for following the paper.\n\n2.\tThe claim that the described limitations of ODE based FM are inherent and unavaoidable is too strong. With the right choice of interpolant and coupling, those limitations could be addressed as in [1,2]\n\n[1] Shrestha, Sagar, and Xiao Fu. \"Diversified Flow Matching with Translation Identifiability.\" Forty-second International Conference on Machine Learning.\n\n[2] Albergo, Michael S., and Eric Vanden-Eijnden. \"Building normalizing flows with stochastic interpolants.\" arXiv preprint arXiv:2209.15571 (2022).\nHowever, ODE based FMs do require special design for different types of constraints and might not be as universal as the proposed DFM. The authors are recommended to include such discussions.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "R0MBpxwmkE", "reviewer": "Reviewer_vrRC", "summary": "The paper extends Flow Matching to Delay Differential Equations via Delay Flow Matching (DFM), where the drift $u(t,x,x_\\tau)$ depends on both the current state and a delayed past state. This lets the model capture trajectory intersections, delay dynamics, and heterogeneous source/target supports while remaining simulation-free to train. The authors instantiate variants with constant vs diverse initial functions and provide empirics on synthetic delay systems, single-cell trajectories, and image tasks, showing consistent gains over ODE-FM baselines.", "questions": "1) How does the required time-step (or number of function evaluations) for stable and accurate integration of DDEs compare to ODEs in your experiments? Do DDEs typically demand finer discretization?\n\n2) Many practical generative systems rely on guidance (e.g., classifier-free guidance). How does your framework integrate with these mechanisms? Does introducing delays change how guidance is applied?\n\n3) Your approach allows different initial function choices. How sensitive is performance to this choice in practice?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The paper is original in extending flow matching to delay differential equations and in using tailored initial-function designs, including constant and diverse variants, to cope with heterogeneous supports and keypoint-guided transport. The theoretical development is solid as the paper establishes a universal-approximation result for continuous transport maps together with a well-motivated conditional training objective. The empirical evaluation is broad, spanning synthetic delayed systems, single-cell trajectory inference, and image generation, with consistent improvements on standard metrics. The writing is generally clear and well structured.", "weaknesses": "1) The practical significance is less clear in settings where the sole goal is high-quality final samples and intermediate trajectories are irrelevant. In such cases, the added machinery of delays and historical states may not translate into noticeably better samples versus strong ODE/SDE baselines. \n\n2) Relatedly, some benefits attributed to the framework (e.g., handling heterogeneous classes) might also be achieved with common conditioning mechanisms such as classifier-free guidance, where the network ingests class or conditioning vectors directly. It is unclear whether the proposed approach would outperform such simpler alternatives.\n\n3) The paper is a little imprecise at times, using phrases like \"approximately push-forward\" without a precise definition, which makes it hard to interpret the guarantees.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "rbTGMbYGfK", "reviewer": "Reviewer_qSLq", "summary": "This paper proposes Delay Flow Matching (DFM), a framework that incorporates a delay term and appropriate initial functions to enable the model to capture delay dynamics, handle trajectory intersections, and transfer between heterogeneous distributions. The paper theoretically proves the universal approximation capability of DFM and demonstrates its effectiveness on a series of synthetic and real-world tasks through various experiments.", "questions": "1. For the experiments, various adaption with KP, OT and Conditional FM are used, the corresponding  cost function is not precised. \n2.  note that the performance of I-DFM(D) in Table 4 outperforms that of OT-DFM(D). However, only OT-DFM was utilized in the preceding sections of the paper. Could you explain the reason for this inconsistency? Additionally, how should we reasonably tune each component of the proposed algorithm in practice?\n    \n2. The diverse initial function appears to be built on the clustering of the dataset. If the clustering results are not particularly optimal (e.g., an excessive or insufficient number of clusters), to what extent would this affect the model’s performance?\n    \n3. All image generation tasks in the paper are designed for scenarios with trajectory intersections and heterogeneous data. If DFM were applied to a standard image generation benchmark task, what performance would you anticipate it to achieve? Would the \"Delay\" mechanism still provide advantages in such a scenario?", "limitations": "", "rating": 4, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "This work seems to be the first to combine DDE with FM. The delay term enables the model to capture delay dynamics and handle trajectory intersections, while the design of appropriate initial functions allows it to handle heterogeneous distributions. The paper provides a compelling theoretical analysis, demonstrating the limitations of FM and proving the expressive power of DFM. The effectiveness of DFM is validated through various experiments.", "weaknesses": "1. Lack of a Systematic Ablation Study: The DFM framework introduces multiple components, including the presence of a delay term, the choice of initial function (constant vs. diverse), the path interpolation method (linear vs. CSpline vs. geodesic), and the coupling strategy (OT-, KP-, I-). While the paper does conduct some comparisons across different experiments, these analyses are scattered and do not form a single, comprehensive ablation study.\n    \n2. Missing discussion and contextualization of related work: The paper fails to discuss the recent research progress in handling trajectory intersections of FM as well as FM for heterogeneous data. A prime example is Switched Flow Matching (SFM), whose core idea is highly similar to this work's: to handle heterogeneous distributions, SFM \"switches\" between different ODEs based on clustering, whereas this paper \"switches\" the initial function. However, the paper lacks both a discussion of this conceptual relationship and a necessary experimental comparison.\n    \n3. Lack of specific guidance for the selection of tau: Tau is a critical parameter of the proposed algorithm, exerting a significant impact on the algorithm’s performance. While the paper presents experiments with substantially different tau values (e.g., tau = 1 or tau = 3), it fails to offer specific guidance on how to select tau in practice (although provides a sensitivity analysis for tau)\n    \n4. Lack of details on efficiency and resource consumption: A major pain point of FM lies in its slow inference speed. While the paper designs experiments primarily focusing on the algorithm’s accuracy, it fails to provide detailed comparisons with other models regarding training and inference time, which are valuable for practical implementation considerations.\n\n5. For image translation experiments, the generated images are in wrong contrast and rather blurry, for which FID is not sufficient to justify the performance. More image quality index should be considered, especially for paired case. Furthermore, the image translation experiments are performed in low-dimensional (such as MINST and Cifar), how the method can be scalable to  high dimension (>=256) \nimage generation problems?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "U7U6SECYkw", "reviewer": "Reviewer_7LHR", "summary": "The paper proposes a new generative-model framework based on the idea of using delay‐differential equations (DDEs) instead of ordinary differential equations (ODEs) for the “flow‐matching” paradigm. Traditional flow matching (FM) uses an ODE to learn a vector field transporting samples from a simple base distribution (e.g., Gaussian) to a target data distribution. The authors identify three limitations of standard FM: inability to model trajectory intersections, and difficulty transferring between heterogeneous distributions. To address this, they introduce Delay Flow Matching (DFM): rather than modelling $\\dot{x}(t)=v(x(t),t)$, they incorporate delay terms $\\dot{x}(t)=v(x(t),x(t-\\tau),t)$ in a neural-parametrised vector field. They further provide a theoretical universal approximation guarantee.", "questions": "1. Could you provide the training time and inference time for DFM and other baseline models?\n\n2. Could you add more baselines (using the papers mentioned above) in the experiments? \n\n3. Could you add a section to discuss related work?", "limitations": "", "rating": 4, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. It clearly articulated the limitation of classic FM models. \n\n2. The authors introduce the delay differential equation (DDE) dynamics. \n\n3. They provide the universal approximation claim that DFM can approximate any continuous tansfer map, providing the theoretical pinning to the method. \n\n4. The experiments include a wide range of data, including synthetic, biological, and natural images.", "weaknesses": "1. There are several papers that addressed the problem of non-intersecting trajectories in rectified flow / flow matching [1, 2, 3, 4]. I think it would be good to review and compare with those related works.  \n\n2. It seems to me that the proposed framework might be computationally more expensive than the classical FM. What is the trade-off between efficiency and quality? \n\n3. I feel the FID value is too high in Table 3 for $\\tau = 0$ (CFM). \n\n4. The FID values for OT-CFM on CIFAR-10 seem to be higher than what earlier works reported. \n\nReferences: \n[1] Park et al. Constant Acceleration Flow, 2024\n[2] Zhang et al. Towards Hierarchical Rectified Flow, 2025\n[3] Chen et al. Gaussian Mixture Flow Matching Models, 2025\n[4] Guo et al. Variational Rectified Flow Matching, 2025", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "z61NvRoX1b", "reviewer": "Reviewer_LNpe", "summary": "The paper proposes Delay Flow Matching (DFM) which is a framework that replaces the role of ODEs within standard flow matching, with flows derived from delay differential equations. \nThis is motivated from the fact that ODEs can fail to represent certain distribution transport plans as they cannot represent flows with intersections (relevant to real world problems such as reconstructing trajectories of single-cell data)\nFrom this, the paper shows how the standard flow machine recipe can be generalised to work with flows from delay differential equations that overcome the previous issues. The paper then provide empirical validation of their framework on different datasets of image data and single-cell data, demonstrating compelling results.", "questions": "* Can you provide further details about the computational profile of DFM compared with flow matching\n\n* How restrictive is the design of the initial function. For example, this requires prior information about the structural of the problem. For cases such as text to image generation, how should we design the initial function - i.e thinking about where we want to generate based on a prompt instead of the fixed classes in the dataset.\n\n* Moreover, do you have experiments where you compare DFM with flow matching where the setup is closer to the standard setup - i.e. the prior is a isotropic Gaussian.", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* The core idea of replacing flow from ODEs with flows from delay differential equations is very interesting due to the properties such flows are afforded. In particular, this address key limitations of standard flow matching around trajectories not being able to intersect.\n\n* The derivation and setup of DFM is nicely presented and formulated, as well as the theoretical justification that standard flow matching cannot handle common transport plans\n\n* The paper provides good empirical evidence to support DFM across a range of modalities.", "weaknesses": "* The paper makes very limited evaluation of the profile of computational cost of training and inference (e.g. FLOPs and wall-clock time instead of just NFEs) for DFM in comparison with flow matching which makes it hard to see whether the increase in performance is worth it in cases such as image generation\n\n* The experimental setup in the image generation differ from the standard setup - i.e. the choice of the prior distribution is made to highlight regimes where DFMs are expected to perform better, which makes it hard to mentally compare against other flow matching approaches.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)