{'fraction_non_alphanumeric': 0.05385943279901356, 'fraction_numerical': 0.015215782983970408, 'mean_word_length': 4.042402387465804, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 22, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Trajectory Inference via Mean-field Langevin in Path Space

Trajectory inference aims at recovering the dynamics of a population from snapshots of its temporal marginals. To solve this task, a min-entropy estimator relative to the Wiener measure in path space was introduced by Lavenant et al. arXiv:2102.09204, and shown to consistently recover the dynamics of a large class of drift-diffusion processes from the solution of an infinite dimensional convex optimization problem. In this paper, we introduce a grid-free algorithm to compute this estimator. Our method consists in a family of point clouds (one per snapshot) coupled via Schr\"odinger bridges which evolve with noisy gradient descent. We study the mean-field limit of the dynamics and prove its global convergence to the desired estimator. Overall, this leads to an inference method with end-to-end theoretical guarantees that solves an interpretable model for trajectory inference. We also present how to adapt the method to deal with mass variations, a useful extension when dealing with single cell RNA-sequencing data where cells can branch and die.


Introduction
Trajectory inference aims at recovering the dynamics of a population of particles given samples from its temporal marginals at various time-points. This problem arises notably in the analysis of single-cell RNA-sequencing data [34,39,14], where one has access-via a destructive measurement process-to the cell state of samples from a large population of cells that evolve in time. In this application, one would like to recover the overall dynamics of the population as well as the trajectories of individual cells so as to improve our understanding of certain biological processes, such as embryonic development or tumor progression.
Trajectory inference can be cast as a regression problem where the unknown is the law of a continuous stochastic process. Let X be the ambient space containing the particles and Ω := C([0, 1], X) the path-space, i.e. the set of all possible continuous trajectories over the time interval [0, 1]. We consider the problem of recovering a probability distribution over paths P ∈ P(Ω), i.e. the law of a continuous stochastic process (X t ) t∈ [0,1] . For the inference to be well-behaved, one should look for a distribution P such that its time marginals are consistent with the observed snapshots and such that the overall dynamics it represents satisfies some notion of regularity.
The problem of trajectory inference has received significant attention over the past several years. However, while pioneering work has focused on creating new methods, the theoretical treatment has remained limited. One class of methods focuses on recovering a potential energy landscape that best fits the observed marginals. For example, Hashimoto et al. [17] encode the potential with a neural network, TrajectoryNet [39] encodes the potential as a neural ODE, and JKOnet [7] encodes the potential with a neural network architecture based on the JKO scheme [20] and input convex neural networks. However, these current approaches are all nonconvex, and therefore it can be difficult to establish rigorous guarantees.
Some guarantees have been established for the equilibrium case by Weinreb et al. [40], and recently consistency was established for the non-equilibrium case by Lavenant et al. [21], who, as a regularity prior, penalize the entropy of P relative to the Wiener measure on Ω (i.e. the law of the Brownian motion). Moreover, Bunne et al. have recently shown how to produce an improved reference process, beyond ordinary Brownian motion [8].
In this work we focus on the estimator introduced by Lavevant et al. [21], which reconstructs the process P by minimizing the entropy relative to the Wiener measure over Ω. By trading off data fitting with regularization, this approach generalizes the approach of Waddington-OT [34] to datasets with many time-points yet possibly few samples per time-point, which is required for consistency [21], yields better performance in practice [33], and is achievable through single-embryo profiling [27]. While the approach of Lavevant et al. leads to a consistent estimator, it also leads to a challenging infinite dimensional convex optimization problem. It is tackled in the original paper by discretizing the space and then applying convex optimization methods. The goal of the present paper is to show that the specific structure of this estimator makes it amenable to a newly introduced class of grid-free stochastic methods called Mean-Field Langevin (MFL) dynamics [19,25], a non-linear generalization of Langevin dynamics which enjoys quantitative global convergence guarantees [28,10]. Instantiated in our context, this method consists in a family of point clouds (one per snapshot) coupled via entropy regularized optimal transport, a.k.a. Schrödinger bridges [23], which evolve with noisy gradient descent. Intuitively, these dynamics can be interpreted as a (non-linear) Langevin diffusion over the path space Ω. Our approach leads to an inference method with end-to-end theoretical guarantees that solves an interpretable model for trajectory inference. We illustrated on Figure 1 the recovered estimator and on Figure 2 the MFL optimization dynamics.
Organization of the paper In § 2, we present the estimator of Lavenant et al. [21]. The heart of our theoretical contributions is in § 3 where we introduce the MFL dynamics and show its quantitative convergence towards the min-entropy estimator at an exponential rate. Numerical experiments are shown in § 4.
Notation and blanket assumptions For two probability measures µ, ν, their relative entropy is H(µ|ν) = log(dµ/dν)dµ if µ ν and +∞ otherwise. For n ∈ N, let [n] := {1, . . . , n}. Throughout, the ambient space X is a compact convex subset of R d and (B t ) t is a reflected Brownian motion on X with B 0 distributed uniformly over X. The path space is C([0, 1]; X) and laws on path space are noted by capital letters, P or R. We denote by P t 1 ,...,t T = (e t 1 , . . . , e t T ) # P the marginal of such laws under the (joint) evaluation maps e t : ω → ω(t). We use boldface letters for family of probability measures µ ∈ P(X) T .

Trajectory Inference as Stochastic Process Inference
Model of population dynamics The first step of any inference task is to determine a prior on the ground truth. For population dynamics with mass conservation, a natural prior is given by drift-diffusion processes. We thus model the population dynamics as where τ > 0 is the temperature/diffusivity, assumed known, and Ψ ∈ C 2 ([0, 1]×X) is unknown. As discussed in prior works [17,21], it is hopeless to recover the divergence-free component of the drift: it is thus natural to assume that the drift is given by the gradient of a function Ψ, called the Waddington potential or epigenetic landscape in the context of cell development. The noise level τ > 0 models the inherent randomness of the ground truth process. It turns out that this noise also has a favorable effect on the algorithm we develop here. This model can be extended to allow particles to branch and die [4], this extension is discussed in Section 4.2.
Model of measurements Let 0 ≤ t 1 < · · · < t T ≤ 1 be a family of measurement times. In contrast to the classical field of inference for stochastic processes [32] where each realization (X t ) t∈[0,1] is observed T times (with T large), in trajectory inference each realization is only observed once, i.e. we observe (X t i ,j ) i∈[T ],j∈[n i ] with independence for each couple (i, j). The measurements are summarized by the family of snapshotsμ t i = 1 Inference problem The general goal is to recover the law P = Law(X) ∈ P(Ω) of the stochastic process given the snapshots (μ t 1 , . . . ,μ t T ). This object P contains all the information about X: in particular it contains the marginals P t of the process and the transition probability between each family of marginals P t 1 ,...,t T . Note that, as discussed in the introduction, some works [39,17,7] focus on directly on recovering the Waddington potential Ψ within a parameterized class of functions. This is a different problem, with encouraging empirical results but weaker known theoretical guarantees.

Min-Entropy Estimator
Let W τ ∈ P(Ω) be the law of the Brownian motion with temperature τ . The min-entropy estimator 1 introduced in [21], is defined as the unique minimizer R * of the functional F : P(Ω) → R defined as where λ > 0 is the regularization strength, ∆t i := (t i+1 − t i−1 )/2 (with the convention t −1 = 0 and t T +1 = 1) and Fit σ is a divergence functional that quantifies how much R t i andμ t i differ, parameterized by a bandwidth σ > 0. For a particular choice of fitting functional (see below), the authors prove the following result: lim T →∞ R * = P weakly, almost surely.
In this paper, we consider the following data fitting term is the Gaussian density and C > 0 is a constant. Eq. (3) defines Fit σ and is always well-defined and finite, but Eq. (4) only makes sense when H(μ t ) < ∞. This is the negative log-likelihood under the noisy observation model X t i ,j = X t i ,j + σZ i,j whereX t i ,j is the observation and Z i,j iid ∼ N(0, I) the noise. Intuitively, when minimizing this soft-min objective over measures R t which are mixtures of Dirac masses, each of the observed particles fromμ t tends to attract the particles from R t that are the closest-even if they are far away in an absolute sense-but barely influences those that are the farthest.
This data-fitting functional Fit σ slightly differs from H(g σ * μ t |R t ), which is the one for which Thm. 2.1 has been established in [21], but it preserves its key properties, namely joint convexity in (R t ,μ t ) and linearity inμ t . Our choice for Fit σ is motivated by a more favorable behavior when R t is a discrete measure and by its natural statistical interpretation. In this paper, we focus on the optimization aspects and leave the statistical analysis of the estimator R * with this specific choice of Fit σ for future work. Let us also mention that Thm. 2.1 has been established for X a compact manifold without boundary, while under our assumptions X has boundaries.
We now turn to the presentation of the algorithm. The data model of Eq. (1) was introduced only to motivate the min-entropy estimator via Thm. 2.1, and plays no role in the rest of the paper.

Optimization Method: Mean Field Langevin Dynamics
In order to compute the min-entropy estimator R * ∈ P(Ω), one needs to minimize an entropy regularized functional over the space of paths F : P(Ω) → R which is of the form where we have posed Fit(R t 1 , . . . , R t T ) : . This is the sum of a convex functional Fit and the relative entropy with respect to the scaled Wiener measure. If, instead of P(Ω), the optimization space was P(X), a problem with this structure could be solved by Mean-Field Langevin (MFL) dynamics. More explicitly, MFL dynamics are designed to minimize problems of the form F (µ) = G(µ) + τ H(µ) where G : P(X) → R is "smooth" and H is minus the differential entropy. They enjoy global convergence guarantees when G is convex. Considering for some m ∈ N * the noisy gradient descent on the function (x 1 , . . . , x m ) → G( 1 m m i=1 δ x i ), MFL dynamics are obtained in the mean-field m → ∞ and vanishing step-size limit, see Section 3.2 for details.
Inspired by this parallel, we now design a drift-diffusion dynamics in path space that converges to the minimizer of F. The main idea is a reformulation of the problem as a "reduced" problem over P(X) T with a structure amenable to MFL dynamics.

Reduced Formulation
We first introduce the reduced functional F , and then state its connection to F in Thm. 3.1.
For µ, ν ∈ P(X), let Π(µ, ν) be the set of transport plans between µ and ν, that is, probability measures on X × X with respective marginals µ and ν. The entropy regularized optimal transport cost between µ and ν is defined, for some τ i > 0, as where p t (x, y) is the transition probability density of the Brownian motion on X over the time interval [0, t] and c τ i (x, y) := −τ i log(p τ i (x, y)). In numerical experiments we use the approximation suggested by Varadhan's formula when τ i is small: c τ i (x, y) ≈ 1 2 y − x 2 (up to an additive constant which is irrelevant when one is interested in minimizers only) [29]. This optimization problem is also known as the Schrödinger bridge problem: its solution gives the most likely evolution of a cloud of particles following a Brownian motion, conditioned on being distributed as µ at t = 0 and ν at t = 1 [23]. Some background on this problem is given in Appendix A, including the definition of the Schrödinger potentials (ϕ, ψ) used hereafter. Consider the function G : P(X) T → R defined for µ = (µ (1) , . . . , µ (T ) ) that represents a family of T reconstructed temporal marginals, by where where is minus the differential entropy of the family of measures µ.
The next result makes the link between minimizing F, the objective in path space (5), and F , the reduced objective (8). It can be interpreted as a Representer Theorem for the min-entropy estimator. 2 This theorem is straightforward to deduce from [5,Cor. 3.5], but the specific form of F that we use here is central for our subsequent algorithmic developments. 3 Theorem 3.1 (Representer Theorem). Let Fit : P(X) T → R be any function.
(ii) Conversely, if F admits a minimizer µ * ∈ P(X) T then a minimizer R * for F is built as where W τ (·|x 1 , . . . , x T ) is the law of W τ conditioned on passing through x 1 , . . . , x T at times t 1 , . . . , t T respectively and R t 1 ,...,t T is the composition of the transport plans γ i,i+1 which are optimal in the definition of The composition of the transport plans is obtained as where we have introduced the conditional probability (a.k.a. "disintegrations") characterized by The importance of the reduced problem comes from the following facts: • The optimization space has been "reduced" from P(Ω) to P(X) T . This reduction is enabled by the Markovian property of W τ . Moreover, Theorem 3.1 gives an explicit method to construct a minimizer for F from a minimizer for F and the associated optimal transport plans (γ i,i+1 ) i∈[T −1] . This fact was already exploited in Lavenant et al. [21].
• The objective function F is the sum of two terms: a "smooth" function G and a differential entropy term τ H. This is precisely the structure tackled by MFL dynamics.
Observe how the entropy in path space H(P |W τ ) is split into two parts: the Schrödinger bridges T τ , included in the "smooth" term G, and minus the differential entropy H.
Let us now describe some useful properties of G and F . Hereafter, the first-variation of G : Proposition 3.2. The function G is convex separately in each of its input (but not jointly), weakly continuous and its first-variation is given for µ ∈ P(X) T and i ∈ [T ] by , with the convention that the corresponding term vanishes when it involves ψ 1,0 or ϕ T,T +1 . The function F is jointly convex and admits a unique minimizer µ * , which has an absolutely continuous density (again denoted by µ * ) characterized by From now on, we thus focus on minimizing the reduced functional F in Eq. (8), keeping in mind that this is sufficient to minimize F, which is our main goal.

Mean-Field Langevin Dynamics
As previously mentioned, the MFL dynamics is natural when it comes to minimize functionals of the form F = G + (τ + )H. Here, we are increasing the entropy factor by > 0 (recalled as an index of F ) which will be useful to obtain convergence guarantees because G is not convex but F 0 = G + τ H is. Using the first-variation V [µ] of G given in Prop. D.1, the MFL dynamics is defined as the solution of the following non-linear SDE of McKean-Vlasov type, for s ≥ 0: with reflected boundary conditions in X. Beware that the pseudo-time s of the optimization dynamics in Eq. (10) should not be mistaken with the pseudo-time t of the process in Eq. (1), which is now represented by the discrete exponents i ∈ [T ]. The family of laws (µ s ) s≥0 of this stochastic process are characterized by the following system of PDEs (understood in the sense of distributions, with no-flux boundary conditions): which are coupled with the quantity ∇V (i) [µ s ]. The link between (10) and (11) follows from Ito's lemma, see e.g. [19]. This is a multi-species PDE where each of the species µ (i) attempts to minimize ∆t i λ Fit σ (·|μ t i ) + (τ + )H via a drift-diffusion dynamics, and is connected to µ (i−1) and µ (i+1) via Schrödinger bridges.

Quantitative Convergence
Let us now state the main convergence result proved in Appendix C. We note that extending this result to the non-compact case raises interesting challenges and is left as an open question. Theorem 3.3 (Convergence). Let µ 0 ∈ P(X) T be such that F (µ 0 ) < ∞. Then for ≥ 0, there exists a unique solution (µ s ) s≥0 to the MFL dynamics (11). For > 0 and moreover assuming that µ 0 has an absolute log-density bounded by A > 0, it holds where C = βe −α/ for some α, β > 0 independent of µ 0 and . Moreover, taking a smooth timedependent s that decays asymptotically asα/ log(s) for someα > α, it holds F 0 (µ s )−F 0 (µ * ) log(log(s))/ log(s) → 0 and µ s converges weakly to the min-entropy estimator µ * .
We can make the following comments: • Under these assumptions, and with fixed, the convergence of µ s to the minimizer of F in relative entropy and in Wasserstein distance, also holds, with the same rate [28,10].
• By Thm. 3.1, one can map (µ s ) s≥0 to a dynamics in P(Ω): it is this dynamics which we refer to as "Mean-Field Langevin in Path Space" in the title.
• This result gives a convergence rate for an optimization dynamics in a non-convex landscape, so it is not even obvious to have global convergence. The convergence rate depends exponentially on 1/ , a standard drawback of Langevin-like dynamics in absence of log-concavity. Summary of theoretical guarantees: Overall, the MFL dynamics of Eq. (11) with simulated annealing ( s → 0) converges to the unique minimizer µ * of the reduced problem in Eq. (8), from which we can build by Thm 3.1 the min-entropy estimator R * ∈ P(Ω) which is optimal for (5). In this sense, our method leads to an inference method with end-to-end theoretical guarantees for trajectory inference. However, this statement is about an idealized dynamics which in practice has to be discretized in time and with particles. In the next section ( §3.4), we will see that the resulting discretization error can be controlled, leveraging the long history of work on mean-field dynamics.

Discretization
In practice, an approximation of the MFL dynamics is obtained by running noisy gradient descent on the function G m : where m ∈ N is be number of particles used to discretized each of the time marginals µ (i) .
Since it can be shown that m∇ where η > 0 is a step-size. The convergence of such a numerical scheme towards the MFL dynamics is standard in the mean-field literature [37], and holds here thanks to the stability of ∇V which follows from [9].
In the related context of wide neural networks, quantitative bounds on the discretization error were studied in [26,Thm. 2]. It is shown that for s = kη, with high probability, |G(μ k ) − G(µ s )| =Õ(e Cs (m −1/2 + η 1/2 )) for some C > 0 independent of k, η and m. Inspecting their proof, it can be seen that it goes through as long as G(µ) and ∇V [µ] can be estimated at the "fast rate" from independent samples of µ, which is our case here [12,Thm. 4.5]. However, the exponential dependence in the pseudo-time s of these bounds make them not precise enough to obtain interesting long-time guarantees. We leave for future work a direct convergence analysis of the fully discretized MFL dynamics (note that [28,Cor. 2] gives a guarantee for the time-discretized MFL dynamics).
Each iteration of Eq. (13) requires to solve an entropic optimal transport problem in order to obtain ∇V . As discussed in Appendix E, this can be solved to precision˜ in time O(m 2 /˜ ) using Sinkhorn's algorithm.

Simulated data
In our first experiment, we compare the behaviour of the MFL dynamics to that of Global Waddington-OT (gWOT) [21] in a setting with few samples per time-point. Although both methods minimize a functional of the form (2), gWOT works on a fixed, discretized support: the union of all observed sample points. Thus, gWOT should perform poorly when the support set has missing regions or "gaps". To simulate this, we simulated a dataset with N i samples at time t i , with N 1 = N 10 = 64 and N i = N, 2 ≤ i ≤ 9 for N ∈ {2 0 , . . . , 2 6 }. The samples are drawn from the marginals of a bifurcating stochastic differential equation (see Appendix G for details).
In Figure 1(left) we illustrate two extreme cases: N = 1 (very few samples at intermediate timepoints) and N = 64 (uniform sampling over time) respectively, with λ = 0.05, λ gWOT = 0.0025. Visually, it seems that the output of MFL is robust to the few-sample regime, with relatively little qualitative difference between the reconstructed trajectories for N = 1, 64. On the other hand, the performance of gWOT degrades visibly once the set of support points is a poor reflection of the underlying process.
To examine this quantitatively, we applied both MFL dynamics and gWOT for various values of N and computed the root-mean-square (RMS) Energy Distance [36] over time to an approximate ground truth (see Appendix G). As shown in Figure 1(b), for large values of N we observe that both MFL and gWOT perform similarly, but with MFL outperforming gWOT for small N as expected.
Finally, Figure 2 illustrates the evolution of the reduced objective F

Dealing with Branching
In the context of single cell RNA-sequencing data it is crucial to be able to deal with the birth and death of cells, a.k.a. branching, which does not occur homogeneously in the domain X. Assume that we dispose of a function g ∈ C 1 ([0, 1] × X) such that g(t, x) is the prior knowledge about the instantaneous growth rate of the distribution of particles in x at time t in the model of Eq. (1). That is, in time dt the probability of a particle X t dividing is g(t, X t )dt [4]. We would like to incorporate this knowledge, and also allow for additional mass variations to account for the inaccuracy of our prior g.
To this end, we proceed heuristically by simply replacing T τ i (µ t i , µ t i+1 ) in Eq. (6) with the following "unbalanced" Schrödinger bridge problem whereμ are probability measures, ρ > 0 is a parameter, M + (X 2 ) is the set of nonnegative measures on X and γ 1 , γ 2 are the marginals of γ. This problem can be solved with a variant of Sinkhorn's algorithm [11,16]. The rationale behind this formula is as follows: (i) the modificationsμ t i of µ t i are intended to approximate the growth process ∂ tμt = g(t, ·)μ t over the time interval [t i , t i+1 ] and (ii)   we relax the marginal constraints, with a parameter ρ > 0, to account for the potential inaccuracy of our prior g. We illustrate on Figure 3 the practical advantage of taking into account branching in the algorithm.
Simulated data experiment Here we consider a modified version of the bistable process from [21, Figure 18] (see Appendix G for details). Figure 3 shows paths sampled from the ground truth, sample points, and the inferred trajectories. Since the lower potential well is closer to the initial condition, we observe in the ground truth that ∼ 72% of paths terminate in the lower branch. This effect is opposed by the high rate of branching in the upper branch, which ultimately results in more particles being sampled in the upper branch. By accounting for branching in the MFL dynamics we are able to reconstruct dynamics that isolate the effect of the potential, as evidenced by the comparable proportion (77%) of inferred paths ending near the lower branch. On the other hand, neglecting the presence of branching results in growth effects being confounded with the potential, and the proportion of paths ending near each branch are roughly equal.

Reprogramming dataset
We now consider the stem cell reprogramming dataset of [34], in which a population of differentiating cells were profiled over a time course using single-cell RNA sequencing. For the purpose of this example, we restrict our attention to days 2.5-6 inclusive making a total of 8 timepoints. As previously, we consider a regime where few samples are observed as snapshots of the time-series and we apply MFL dynamics to infer reconstructed marginals. We reason that if the underlying process is well described by a Schrödinger bridge (as was empirically validated in [34]), then MFL dynamics should be able to "improve" on the sample marginals. As a proxy for ground truth, we used the full dataset (consisting of 59154 cells over the 8 timepoints). From this, we subsampled timepoints consisting of 100 cells at days 2.5 and 6, and 10 cells at the remaining intermediate timepoints. After carrying out a series of preprocessing steps (see Appendix H for details), we applied MFL dynamics to the subsampled data to produce reconstructed marginals. Figure 4(left) shows an example of the sampled points, MFL reconstruction, and the full dataset.
A reconstruction "error" was then computed as the Energy Distance (20) of each reconstructed marginal to the corresponding snapshot in the full dataset. Figure 4(right) shows the error for each timepoint for each of MFL and gWOT, for the best performing parameters found in a parameter sweep. We observe that both MFL and gWOT improve upon the subsampled snapshots, but MFL by a larger margin.

Conclusion
We introduced a grid-free numerical method for trajectory inference that computes the minentropy estimator introduced in [21], with global and quantitative convergence guarantees in the mean-field limit. This method arises naturally when decomposing the optimization problem in a suitable way and, in practice, outperforms the fixed-grid method of [21].
Concerning limitations, our method shares those of the min-entropy estimator : since it does not incorporate fine structural prior on the structure of the Waddington potential (i.e. it is fully nonparametric), it may suffer from a limited statistical efficiency in high dimension. An interesting research direction is to quantify this statistical performance and to adapt our algorithm to learn the Waddington potential Ψ in a structured parameterized class of functions, jointly with the law on paths P .
On a more abstract level, our main insight is that min-entropy problems in Wiener space can be solved via a multi-species diffusion dynamics coupled with Schrödinger bridges. This point of view raises interesting theoretical questions, e.g. can we rigorously interpret our dynamics as a diffusion in path space?

A Background on the Schrödinger bridge problem
The following facts can be found for instance in the lecture notes [30,Sec. 4]. Consider µ, ν ∈ P(X) and the entropy-regularized optimal transport problem defining T τ (µ, ν) in Eq. (6) with a cost function c ∈ C ∞ (X × X). This problem admits a unique solution γ * and admits a dual formulation The dual problem admits a unique solution in L 1 (µ) × L 1 (ν) up to the transformation (ϕ + κ, ψ − κ) for κ ∈ R. Moreover, we have that T τ (µ, ν) = ϕdµ + ψdν and the primal-dual relation Also, the potentials satisfy the following equations, which are the first-order optimality conditions: These equations are a priori only satisfied µ (resp. ν) almost everywhere, but they can be used to extend ϕ and ψ as continuous (in fact infinitely differentiable) functions over X, which satisfy these equations everywhere and are unique in C ∞ (X), up to the additive invariance mentioned above. Throughout the paper, we refer to such functions (ϕ, ψ) as the Schrödinger potentials (they are also referred as EOT potentials in [30]). Let us conclude this section with two observations.
(i) One can bound the oscillation of ϕ as follows This is obtained by upper bounding c inside the integral which leads to ϕ(x) ≤ sup x,y c(x, y) − τ log e (ψ(y))/τ dν(y) . Subtracting the analogous lower bound, we observe that the log terms cancel and we get the bound on the oscillation.
(ii) One can differentiate (16) to see that we have for x ∈ X with Law(X, Y ) = γ * and γ * (dy|x) = g(x, y)ν(dy) is the conditional distribution of Y given X.

B Proof of Theorem 3.1
Let us first recall the statement of Thm. 3.1.
(ii) Conversely, if F admits a minimizer µ * ∈ P(X) T then a minimizer for F is built as where W τ (·|x 1 , . . . , x T ) is the law of W τ conditioned on passing through x 1 , . . . , x T at times t 1 , . . . , t T respectively and R t 1 ,...,t T is the composition of the transport plans γ i,i+1 which are optimal in the definition of T τ i (µ * (i) , µ * (i+1) ), for i = 1, . . . , T .
This theorem is a direct consequence of the following lemma, which is similar in spirit to [21,Prop. B.1], but the terms involved are different.
Lemma B.2. There exists C > 0 such that, for any R ∈ P(Ω) and t 1 , . . . , t T a collection of instants, it holds

The first inequality ( †) becomes an equality if and only if
where W τ (·|x 1 , . . . , x T ) is the law of W τ conditioned on passing through x 1 , . . . , x T at times t 1 , . . . , t T respectively. In addition, the second inequality ( ) becomes an equality if and only if R is Markovian.
Proof. The first inequality ( †) and the equality case follows from the additivity property of the relative entropy under conditioning, see [23,Eq. (A9)]. For the second inequality ( ), [5,Sec. 3.4] states that with equality if and only if R t 1 ,...,t T is Markovian. This is the formula of [21, Prop. B.1], but this expression is unsuitable for our purposes and we need to further reorganise the terms in E.
Without loss of generality, let us assume that R t i are absolutely continuous with density dRt i dx (x) = r i (x) (if this is not the case, both sides of the inequality are infinite) and let V X be the Lebesgue volume of X. On the one hand, since W τ t i is the normalized volume measure on X, it holds H(R t i |W τ t i ) = H(R t i ) + log(V X ). On the other hand, letting by definition of the transition probability density p of the Brownian motion on X. It follows that for any µ, ν ∈ P(X) with finite differential entropy and γ ∈ Π(µ, ν) it holds where we used the fact that γ ∈ Π(µ, ν) to simplify the two last terms (see [24,Lem. 1.6] for more details on the change of reference measure in regularized optimal transport). Putting everything together, and using the fact that which proves the formula.

C Proof of Theorem C.1
Let us recall the statement of Theorem C.1 and prove it, by an application of a convergence result proved independently in [28] and [10].
Proof. We verify the assumptions of [10], noticing that their proof can be adapted without difficulty to our context of families of T probability measures with T ≥ 1.
For the well-posedness of the dynamics, we split the objective function F as G + (τ + )H. Assumption 1, about the stability and regularity of the first-variation V , is guaranteed Prop. C.2 below (the stability and regularity of the component δFit/δµ is immediate) and leads to the well-posedness of the dynamics without requiring > 0 (in fact we could even take ∈ [−τ, +∞[ and have well-posedness).
For the global convergence guarantee we split the objective function F as (G + τ H) + H because we need the fact that F 0 = G + τ H is convex. Assumption 2, which requires convexity of F 0 and existence of a minimizer for F , is satisfied thanks to Prop. D.1. For the uniform Log-Sobolev Inequality (LSI) (Assumption 3), we first remark that [22,Thm. 7.3] states that the uniform distribution over X satisfies LSI with a constant ρ 0 that only depends on D the diameter of X. Now, the i-th component of the first-variation of F 0 is given by V (i) [µ] + log(µ i ). By the expression of Prop. D.1, the oscillation osc(V (i) [µ]) of V (i) [µ] (i.e. the difference between its maximum and minimum value over X) is bounded (independently of and µ 0 ). Indeed, the gradient formula for δFit σ /δµ (i) is nonnegative and bounded by e D 2 /(2σ 2 ) and by App A Eq. (17), the Schrödinger potential ϕ i,i+1 has an oscillation bounded by sup x,y c τ i (x, y) − inf x,y c τ i (x, y) which is D 2 /2 when c(x, y) = 1 2 y − x 2 or a less explicit constant that depends only on the domain X for the cost c τ i . Combining these estimates with Lemma C.3, we get that the oscillation of V (i) [µ] + log(µ (i) ) is bounded by some κ > 0 independent of and µ 0 (under the assumption that the log-density of µ 0 is absolutely bounded by A > 0).
Finally, the convergence result with simulated annealing is a direct application of [10,Thm. 4.1] which proves the convergence rate in log(log s)/ log s. Since in addition the sublevel sets of F 0 are compact and the minimizer µ * is unique, standard considerations imply that µ s converges weakly to µ * .
Let us report a stability result concerning the Schrödinger potentials, which is a consequence of a more general result in [9] and is used in the proof above.
where vol is the uniform distribution over X, and by Ito's formula Since, for S fixed, all the quantities in the exponential are uniformly bounded (in particular, for the term involving ∂ s V (i) [µ s ] this follows from Prop. C.2 which implies where v s is the Wasserstein derivative of (µ (i) s ) s and using facts from gradient flows theory [3]). This shows that the transition kernel of the process X (i) is upper and lower bounded by the heat kernel over X with Neumann boundary conditions. In particular, the density of µ Proposition D.1. The function G is convex separately in each of its input (but not jointly), weakly continuous and its first-variation is given for µ ∈ P(X) T and i ∈ [T ] by where (ϕ i,j , ψ i,j ) ∈ C ∞ (X) are the Schrödinger potentials for T τ i (µ (i) , µ (j) ), with the convention that the corresponding term vanishes when it involves ψ 1,0 or ϕ T,T +1 . The function F is jointly convex and admits a unique minimizer µ * , which has an absolutely continuous density (again denoted by µ * ) characterized by Proof. The properties and formulas for G follows from those for T τ which are well-known (see [15]) and for Fit which are immediate. In particular, T τ is convex separately in each of its variables as a supremum of linear forms but not jointly because of the product measure in Eq. (15). The joint convexity of F and the uniqueness of its minimizer is proved in [21] and the characterization of the minimizer can be formally deduced by writing the first order optimality condition V (i) [µ * ] + τ log(µ * (i) ) = 0. The rigorous argument, which is standard, can be found in a similar context e.g. in [25,Lem. 10.4].]

E Solving the Schrödinger Bridge Problems with Sinkhorn's algorithm
At each iteration, in order to compute V [μ[k]], one needs to compute the Schrödinger potentials (ϕ i,i+1 , ψ i,i+1 ) associated to the T − 1 Schrödinger bridges problems T τ (μ (i) ,μ (i+1) ) (see Prop. D.1). Among the various algorithms that can solve this problem [31], let us focus our discussion on the well-studied Sinkhorn's algorithm, which is alternate block maximization on the dual of Eq. (6).
Given two discrete probability measuresμ m = m i=1 p i δ x i andν m = m i=1 q i δ y i we define the cost matrix with entries c i,j = c(x i , y j ) (which we approximate with 1 2 x i − y j 2 using Varadhan's formula). The iterates u[ ], v[ ] ∈ R m , ≥ 1 of Sinkhorn's algorithm are defined as : This algorithm converges in value at a rate O(1/(τ k)), see [13,2]. In practice, Sinkhorn's iterations could be further sped up with non-linear acceleration methods for faster convergence [38]. Upon convergence, one can recover the Schrödinger potential ϕ, ψ ∈ C ∞ (X) via the formula Moreover the minimizer in Eq. (6), needed to recover P * by Thm. 3.1, is given by These alternatives are discussed in [1,Sec. 5.3] where (ϕ, ψ) are computed as a subroutine. There, the conclusion is that option (ii) is slightly more efficient, and this is the method we implemented.

F Simulated Annealing
A standard heuristic to accelerate diffusion-based algorithms with a temperature parameter τ is the so-called simulated annealing method [6], which consists in starting from a large value τ 0 and slowly decreasing it towards the desired value while the algorithm runs. In our context, a larger value for τ accelerates the convergence of both Sinkhorn's algorithm and the MFL dynamics.
In practice, we used simulated annealing with a geometric decay τ t = max{cr t , τ f } for some r ∈ ]0, 1[ and c > 0 in order to quickly reach the desired temperature τ f > 0. Empirically, we find also that allowing the step size η and the (squared) data-fitting bandwidth σ 2 to scale with the temperature leads to further acceleration of the MFL dynamics. We illustrate this for the example of Section 4.1 in Figure 5. In this experiment, particles of the MFL are started from N(0, 1.0). Consequently, some particles are distant from the data and require MFL to be run for an large number of iterations to converge. On the other hand, simulated annealing for the first 500 iterations with τ 0 = 5τ allows for the MFL dynamics to convergence much faster. As can be seen in Figure 5(b), the result of MFL with annealing is slightly less noisy. To prevent noise at high temperature from causing particles in the MFL dynamics to stray far away from the observed data, we add a confining potential to the objective (7) where we have written R = T −1 T i=1 R t i andμ = T −1 T i=1μ t i to be the mixtures (over time) respectively of the reconstructed and observed marginals. This has the effect of penalizing particles in the MFL dynamics which are far away from any observation. In the annealing examples, we used σ = 5.0 for the confining potential bandwidth.

G Simulation details
All numerical experiments were run using a CPU-based implementation of MFL dynamics. A copy of the code to reproduce the figures in this article is available at https://github.com/ zsteve/mfl.
The discretized MFL dynamics described in Section 3.4 was applied to each simulated dataset with m = 100 particles per timepoint and using the data fitting functional (3) with bandwidth σ = 0.5. Particles in the MFL dynamics were started from N(0, 0.1 2 ) and evolved following (13) with η = 0.1 for 2, 500 iterations. We repeated this for the following parameter values and random seeds: The approximate ground truth is taken to be a simulated dataset with 500 particles per timepoint. That is, we computed (T −1 T i=1 D 2 (µ t i , R t i )) 1/2 where µ t i and R t i are respectively the ground truth and reconstructed marginals at time t i , and the squared Energy Distance between two measures (α, β) is defined to be G.2 Additional details for Section 4.2 The ambient space X is taken as in [21], and dynamics follow (1)  . We prescribe a growth rate g(t, x) = 10(tanh(2x 0 ) + 1)/2 so that particles grow faster in the region x 0 > 0. Particles are started from X 0 ∼ N(0, 0.1 2 ), and 10 timepoints with 50 particles each were sampled on the interval t ∈ [0, 0.5]. We fit MFL dynamics both with and without the modification for branching described in Section 4.2 with λ = 0.025 and ρ = +∞ (since the growth rate is known exactly). Other hyperparameters were taken to be the same as in Section 4.1.

H Reprogramming dataset pre-processing details
For each set of subsampled snapshots, expression matrices were first centered and projected into 10 PCA dimensions. Similarly as in [21], a set of scaling factors were computed: • Pairwise scaling factors: σ 2 scale = Eμ t i ,μt i+1 Xt i+1 −Xt i The cost function for each transport term in (6) was divided by σ 2 scale so as to be of order one. Similarly, the pairwise squared distances in the data-fitting functional (3) were divided by η 2 scale . The timepoints were mapped to 0 = t 1 ≤ . . . ≤ t T = 1 and the value of τ was chosen such that the effective transport regularization level (τ i in (6)) was 0.1 for transport over 0.5-day interval. We applied MFL dynamics for λ ∈ {0.0125, 0.025, 0.05, 0.1, 0.2} and other parameters as in Section 4.1, and gWOT for λ ∈ {0.000625, 0.00125, 0.0025, 0.005, 0.01}. Of these parameter values, we found that λ = 0.025 performed best for MFL in terms of Energy Distance to the full dataset (projected onto the previously calculated principal components) and similarly λ = 0.01 for gWOT.