Force-Guided Bridge Matching for Full-Atom Time-Coarsened Dynamics of Peptides

## Abstract

Abstract Molecular Dynamics (MD) is crucial in various fields such as materials science, chemistry, and pharmacology to name a few. Conventional MD software struggles with the balance between time cost and prediction accuracy, which restricts its wider application. Recently, data-driven approaches based on deep generative models have been devised for time-coarsened dynamics, which aim at learning dynamics of diverse molecular systems over a long timestep, enjoying both universality and efficiency. Nevertheless, most current methods are designed solely to learn from the data distribution regardless of the underlying Boltzmann distribution, and the physics priors such as energies and forces are constantly overlooked. In this work, we propose a conditional generative model called Force-guided Bridge Matching (FBM), which learns full-atom time-coarsened dynamics and targets the Boltzmann-constrained distribution. With the guidance of our delicately-designed intermediate force field, FBM leverages favourable physics priors into the generation process, giving rise to enhanced simulations. Experiments on two datasets consisting of peptides verify our superiority in terms of comprehensive metrics and demonstrate transferability to unseen systems.

## 1 Introduction

Figure: Figure 1: Illustration of how molecular conformations transfer from one state to another by MD (path in white) and time-coarsened dynamics (path in black).
Refer to caption: x1.png

Molecular Dynamics (MD), which simulates the physical movements of molecular systems at the atomic level via *in silico* methods, are widely applied in the fields of materials science, physics, chemistry, and pharmacology (Wolf et al., ; Durrant & McCammon, ; Salo-Ahen et al., ).
Accurate MD simulations enable the researcher to comprehend the equilibrium thermodynamics and kinetics of different molecular phases without the need for expensive wet-lab experiments.

Conventional MD software, like AMBER (Pearlman et al., ) and CHARMM (Vanommeslaeghe et al., ), mostly pre-defines an empirical force field of molecular systems and performs MD based on the numerical integration of Newtonian equations over the timestep $\Delta{t}$. To minimize discretization errors, $\Delta{t}$ should be chosen small enough, typically on the order of femtoseconds ($10^{-15}$s). Consequently, simulating the full duration of critical phase transitions that occur on the microsecond ($10^{-6}$s) or even millisecond ($10^{-3}$s) timescale becomes virtually impossible within a reasonable wall-clock time. To overcome the efficiency limitation in long-timescale simulation, a surge of approaches have been devised, including MD-like methods (Voter et al., ; Pang et al., ; Laio & Parrinello, ) and Monte Carlo-based methods (Sadigh et al., ; Sugita & Okamoto, ; Neyts & Bogaerts, ). However, all these methods share a common drawback: different molecular systems require customized simulations, despite the fact that many atomistic systems should, in principle, exhibit similar dynamic mechanisms.

Recently, data-driven approaches (Klein et al., ; Schreiner et al., ; Hsu et al., ) leveraging deep generative models have been used to enhance MD simulations. Unlike traditional numerical integration, these methods directly learn MD from diverse observed trajectories, offering superior transferability across different molecular systems. Additionally, they can accelerate simulations through the use of time-coarsened dynamics, where models learn to generate new states after a significantly larger timestep ($\tau\gg\Delta t$) starting from an initial state, as illustrated in Figure 1. Nevertheless, existing learning-based methods still face two significant issues. The first issue is that most current methods are designed solely to learn from the data distribution, which may be biased compared to the underlying Boltzmann distribution (Boltzmann, ) — a fundamental concept for describing the thermal equilibrium of molecular systems in relation to atomic positions and velocities. Although Timewarp (Klein et al., ) alleviates this issues by employing the Metropolis-Hastings algorithm (Metropolis et al., ) to resample from the equilibrium, the unbearably low acceptance rate still hinders its applicability.
The second issue is that physics priors (*e.g.*, energies and forces) are constantly overlooked during the learning process, which, yet, are critical for providing insights into the dynamics of molecular systems.

In this work, we delicately design a generative model to learn time-coarsened dynamics, and target the Boltzmann-constrained distribution $p(\vec{{\mathbf{X}}})$. This distribution incorporates a potential energy term $\exp(-k\varepsilon(\vec{{\mathbf{X}}}))$, into the original data distribution $q(\vec{{\mathbf{X}}})$, resulting in $p(\vec{{\mathbf{X}}})=\frac{1}{Z}q(\vec{{\mathbf{X}}})\exp(-k\varepsilon(\vec{
{\mathbf{X}}}))$. Here, $\vec{{\mathbf{X}}}$ stands for the molecular conformation, and $Z$ is the normalization factor. The definition of $p(\vec{{\mathbf{X}}})$ ensures that regions of high density in the generated distribution correspond to low potentials $\varepsilon(\vec{{\mathbf{X}}})$, thereby aligning more closely with thermodynamic principles.
To learn $p(\vec{{\mathbf{X}}})$, we utilize the bridge matching generative framework (Shi et al., ) in order to make the generation to be conditioned on the starting molecular conformations rather than Gaussian priors.
The most challenging part is that, the vector field to generate $p(\vec{{\mathbf{X}}})$ is strongly associated with the virtual “energies” and “forces” during the generation process, which, yet, are nontrivial to obtain. To address the issue, we derive an effective and rigorous way to interpolate a well-designed intermediate force field into the bridge matching framework based on reasonable assumptions. The proposed framework is termed as Force-guided Bridge Matching (FBM). It is also worth noting that FBM employs TorchMD-NET (Pelaez et al., ) as the backbone model, which enables full-atom modeling while preserving $\mathrm{SO}(3)$-equivariance. Our main contributions are summarized as follows:

- 1.
To our best knowledge, FBM is the first full-atom generative model for time-coarsened dynamics, by directly targeting the Boltzmann-constrained distribution without extra resampling steps.
- 2.
We interpolate a well-designed intermediate force field into the bridge matching framework, which effectively involves the guidance of physics priors for better MD simulation.
- 3.
We evaluate FBM on two datasets consisting of small peptides, where FBM exhibits transferability to unseen peptide systems and consistently showcases state-of-the-art results across various peptides.

## 2 Related Work

#### Boltzmann Generator

An important objective of MD research is to quickly sample from the Boltzmann distribution, thereby revealing the free energy landscape and collective variables of matter of molecular systems. Boltzmann generators (Noé et al., ; Köhler et al., ; ; Falkner et al., ; Klein et al., ; Klein & Noé, ) employ generative models to produce samples asymptotically from the Boltzmann distribution mainly by: (i) Apply reweighting techniques to i.i.d. generated samples. (ii) Inference in a Markov Chain Monte Carlo (MCMC) procedure. These approaches heavily rely on MC resampling techniques, which are the bottleneck of the sampling efficiency due to costly energy calculation and low acceptance rates. Most similar to our work, Wang et al. () propose ConfDiff that incorporate the energy and force guidance directly into score diffusion to target the Boltzmann-constrained distribution, yet it only works well on protein backbones and generate conformations in an unconditional way that fails to capture temporal dynamics. In contrast, FBM introduces the force guidance into the bridge matching framework, making it possible to sample straightforwardly from the Boltzmann-constrained distribution.

#### Time-Coarsened Dynamics

To overcome the instability of numerical integration of conventional MD simulations, many deep learning methods have adopted the fashion of time-coarsened dynamics, where models learn dynamics of diverse molecular systems over a long timestep. Fu et al. () proposes a multi-scale graph network to learn dynamics of polymers, which fails to operate in the all atom system. ITO (Schreiner et al., ) is devised to learn the transition probability over multiple time resolutions, yet its transferability across chemical space remains unknown. Recently, Timewarp (Klein et al., ) and TBG (Klein & Noé, ) utilize augmented normalizing flows and flow matching respectively for transferable time-coarsened dynamics of small peptides, while both need additional reweighting procedures to debias expected values to the Boltzmann distribution. Moreover, Score Dynamics (Hsu et al., ) applies score diffusion to capture dynamic patterns merely from the data distribution. On the contrary, FBM is designed to learn time-coarsened dynamics at the atomic level and target the Boltzmann-constrained distribution by directly inference without extra steps, which aligns well to thermodynamic principles and exhibits transferability to unseen peptide systems as well.

## 3 Method

In this section, we will present the overall workflow of our method. Specifically, in § 3.1, we will define our task formulation of time-coarsened dynamics, with providing necessary notations. Then in § 3.2, we introduce how to learn the time-coarsened dynamics from the data distribution, via a conditional generative model based on bridge matching, which is dubbed as FBM-BASE. Built upon this baseline, in § 3.3, we further propose a novel generative framework that targets the Boltzmann-constrained distribution by incorporating an intermediate force field into the FBM-BASE model, leading to our eventual model FBM.
All proofs of propositions are provided in § B.

### 3.1 Molecular Representation

We represent each molecule (namely peptide in our experiments) as a graph ${\mathcal{G}}=({\mathcal{V}},{\mathcal{E}})$ consisting of the node set ${\mathcal{V}}$ and the edge set ${\mathcal{E}}$. For a molecule with $N$ atoms (including hydrogens), ${\mathcal{V}}=\{v_{0},\cdots,v_{N-1}\}$ where $v_{i}$ ($0\leq{i}<N$) represents the $i$-th atom of the molecule. Each node $v_{i}$ is further attributed with Cartesian coordinates $\vec{{\bm{x}}}_{i}\in{\mathbb{R}}^{3}$ from the structural information and node features ${\bm{z}}_{i}\in{\mathbb{R}}^{H}$ from the embedding of atom types, where $H$ represents the hidden dimension. Particularly for peptides that are composed of 20 natural amino acids and exhibit similar features, we construct the atom type vocabulary based on the atom nomenclature of Protein Data Bank (Berman et al., ) to obtain more refined atomic representations. Formally, the features and Cartesian coordinates of all nodes are concatenated as node representations:

$$ $\vec{{\bm{X}}}=[\vec{{\bm{x}}}_{0},\cdots,\vec{{\bm{x}}}_{N-1}]^{\top}\in{ \mathbb{R}}^{N\times{3}},~{}{\bm{Z}}=[{\bm{z}}_{0},\cdots,{\bm{z}}_{N-1}]^{ \top}\in{\mathbb{R}}^{N\times{H}}.$ (1) $$

The edges are constructed with the cutoff radius $r_{\mathrm{cut}}$. For any node pair $(v_{i},v_{j})$, the connection is established iff $||\vec{{\bm{x}}}_{i}-\vec{{\bm{x}}}_{j}||<r_{\mathrm{cut}}$. Constructing cutoff graphs is a favorable choice for modelling molecular systems since forces and chemical bonds are highly related to interatomic distances.
With the aforementioned notations, we provide the task formulation below.

#### Task Formulation

Given each MD trajectory, we extract molecule pairs ${\mathcal{D}}=\{({\mathcal{G}}_{s},{\mathcal{G}}_{s+\tau})\mid s\in{\mathcal{S
}}\}$ to create a training dataset, where ${\mathcal{G}}_{s}$ denotes the starting state at time $s$ and ${\mathcal{G}}_{s+\tau}$ represents the future state after a temporal interval $\tau$.
Our goals are: (i) Train a baseline model that fits the conditional data distribution $\mu(\vec{{\bm{X}}}_{s+\tau}|\vec{{\bm{X}}}_{s})$ from ${\mathcal{D}}$, denoted as $q$. (ii) Based on the trained baseline model, we further train a new model that admits the Boltzmann constraint: $p(\vec{{\mathbf{X}}})=\frac{1}{Z}q(\vec{{\mathbf{X}}})\exp(-k\varepsilon(\vec{
{\mathbf{X}}}))$, where $k$, $\varepsilon$, and $Z$ denote the inverse temperature, the potential of the molecular system, and the partition function, respectively. During the training process, we require the queries of the MD energies and forces for each pair $({\mathcal{G}}_{s},{\mathcal{G}}_{s+\tau})$, namely, $(\varepsilon(\vec{{\bm{X}}}_{s}),\varepsilon(\vec{{\bm{X}}}_{s+\tau}))$ and $(\nabla\varepsilon(\vec{{\bm{X}}}_{s}),\nabla\varepsilon(\vec{{\bm{X}}}_{s+
\tau}))$. In the following context, we will train three neural networks $v_{\theta}$, $u_{\theta}$, and $w_{\theta}$ to approximate the vector fields that are used to generate the distributions $q$ and $p$.

### 3.2 Baseline Bridge Matching

We will present the baseline BM model (*i.e.*, FBM-BASE) to fit the conditional data distribution $\mu(\vec{{\bm{X}}}_{s+\tau}|\vec{{\bm{X}}_{s}})$ from dataset ${\mathcal{D}}$. For the sake of convenience and consistency with the continuous diffusion time between $[0,1]$, we refer to the prior graph as ${\mathcal{G}}_{0}$ and the target graph as ${\mathcal{G}}_{1}$ in what follows, and denote their corresponding coordinates as $\vec{{\bm{X}}}_{0}$ and $\vec{{\bm{X}}}_{1}$ accordingly.

#### Bridge Matching Framework

We leverage the generative framework of bridge matching (Shi et al., ), which learns a mimicking diffusion process between two arbitrary distributions, allowing for more flexible choices of priors. Let $q_{0},q_{1}$ denote the prior and target data distributions, respectively, and $t\in[0,1]$ be the continuous diffusion time. We define the forward process:

$$ $d\vec{{\mathbf{X}}}_{t}=f_{t}(\vec{{\mathbf{X}}}_{t})dt+\sigma_{t}d{\mathbf{B} }_{t},~{}\vec{{\mathbf{X}}}_{0}\sim{q_{0}},$ (2) $$

where $\vec{{\mathbf{X}}}_{t}$ and ${\mathbf{B}}_{t}$ represent the random variable and the Brownian motion of the diffusion process at time $t$, respectively. With the process pinned down at an initial point $\vec{{\bm{X}}}_{0}$ and terminal point $\vec{{\bm{X}}}_{1}$, the conditional distribution $q_{t}(\cdot|\vec{{\bm{X}}}_{0},\vec{{\bm{X}}}_{1})$ will be realized as a diffusion bridge in the form of $d\vec{{\mathbf{X}}}_{t}=\{f_{t}(\vec{{\mathbf{X}}}_{t})+\sigma_{t}^{2}\nabla
\log{q_{t}(\vec{{\bm{X}}}_{1}|\vec{{\mathbf{X}}}_{t})}\}dt+\sigma_{t}d{\mathbf
{B}}_{t}$ with $\vec{{\mathbf{X}}}_{0}=\vec{{\bm{X}}}_{0}$, where Doob $h$-transform theory (Rogers & Williams, ) guarantees $\vec{{\mathbf{X}}}_{1}=\vec{{\bm{X}}}_{1}$. For simplification, by considering $f_{t}=0$ and $\sigma_{t}=\sigma$, the process will degenerate to the following Brownian bridge:

$$ $d\vec{{\mathbf{X}}}_{t}=\frac{\vec{{\bm{X}}}_{1}-\vec{{\mathbf{X}}}_{t}}{1-t} dt+\sigma{d{\mathbf{B}}_{t}},~{}\vec{{\mathbf{X}}}_{0}=\vec{{\bm{X}}}_{0},$ (3) $$

which yields the conditional distribution $q_{t}(\vec{{\mathbf{X}}}_{t}|\vec{{\bm{X}}}_{0},\vec{{\bm{X}}}_{1})$ at time $t\in[0,1]$, with its marginal defined as $q_{t}(\vec{{\mathbf{X}}}_{t})$. The core of bridge matching is to find a Markov diffusion governed by a vector field $v$:

$$ $d\vec{{\mathbf{X}}}_{t}=v(\vec{{\mathbf{X}}}_{t},t)dt+\sigma{d{\mathbf{B}}_{t}},$ (4) $$

which admits the same marginal $\vec{{\mathbf{X}}}_{t}\sim{q_{t}}$, such that $\vec{{\mathbf{X}}}_{1}\sim{q_{1}}$ holds. To achieve this, we can learn a parametric vector field $v_{\theta}$ via the following regression loss:

$$ ${\mathcal{L}}_{\mathrm{fwd}}=\mathbb{E}_{t\sim\mathrm{Uni}(0,1),({\mathcal{G}} _{0},{\mathcal{G}}_{1})\sim{\mathcal{D}},\vec{{\bm{X}}}_{t}\sim{q_{t}(\cdot| \vec{{\bm{X}}}_{0},\vec{{\bm{X}}}_{1})}}[||\frac{\vec{{\bm{X}}}_{1}-\vec{{\bm{ X}}}_{t}}{1-t}-v_{\theta}(\vec{{\bm{X}}}_{t},t)||^{2}],$ (5) $$

where $\mathrm{Uni}(0,1)$ denotes the uniform distribution of $[0,1]$, and $v_{\theta}$ is implemented by a neural network parameterized by $\theta$. The calculation of Eq. 5 requires sampling from $q_{t}(\vec{{\mathbf{X}}}_{t}|\vec{{\bm{X}}}_{0},\vec{{\bm{X}}}_{1})$, which usually needs additional SDE simulations. Fortunately, Eq. 3 enables a closed-form solution as follows:

$$ $q_{t}(\vec{{\mathbf{X}}}_{t}|\vec{{\bm{X}}}_{0},\vec{{\bm{X}}}_{1})={\mathcal{ N}}(t\vec{{\bm{X}}}_{1}+(1-t)\vec{{\bm{X}}}_{0},t(1-t)\sigma^{2}\bm{I}),$ (6) $$

where we can sample $\vec{{\bm{X}}}_{t}$ efficiently during each training step with any given $t$.

The optimal vector field $v_{\theta}$ to minimize the loss in Eq. 5 is actually equal to the expectation $v^{*}(\vec{{\mathbf{X}}}_{t},t)=\mathbb{E}_{q_{t}(\cdot,\cdot|\vec{{\mathbf{X}
}}_{t})}[\frac{\vec{{\bm{X}}}_{1}-\vec{{\mathbf{X}}}_{t}}{1-t}]$. Then the distribution $q_{t}$ can be estimated by performing SDE sampling in Eq. 4, using $v_{\theta}$ as the vector field $v$. In our experiments, $v_{\theta}$ is built upon TorchMD-NET (Pelaez et al., ) for full-atom modeling while satisfying $\mathrm{SO}(3)$-equivariance. Details of the architecture of neural networks used in our model are further elucidated in § C.

### 3.3 Force-guided Bridge Matching

In most cases, the training datasets are biased from the underlying Boltzmann distribution, leading to defective prediction even for superior generative models. In this section, we will introduce how to train a force-guided generative model FBM, which admits the Boltzmann-constrained distribution denoted as $p_{1}(\vec{{\mathbf{X}}}_{1})=\frac{1}{Z}q_{1}(\vec{{\mathbf{X}}}_{1})\exp(-k
\varepsilon(\vec{{\mathbf{X}}}_{1}))$. Here, the exponential term $\exp(-k\varepsilon(\vec{{\mathbf{X}}}_{1}))$ serves as the physics prior from thermodynamic principles. The overall framework of FBM as well as FBM-BASE is illustrated in Figure 2.

Figure: Figure 2: The overall framework of FBM-BASE and FBM. A. Firstly, FBM-BASE leverages the bridge matching framework to learn time-coarsened dynamics from the data distributions $q_{0}$ and $q_{1}$. B. With the guidance of the intermediate force field $\nabla\varepsilon_{t}$ at diffusion time $t$, the marginal distribution admits $p_{t}(\vec{{\mathbf{X}}}_{t})\propto{q_{t}(\vec{{\mathbf{X}}}_{t})\exp(-k \varepsilon_{t}(\vec{{\mathbf{X}}}_{t}))}$, thereby the target distribution of FBM is debiased to the Boltzmann-constrained distribution $p_{1}$.
Refer to caption: x2.png

#### Force-guided Bridge Matching Framework

In order to learn $p_{1}$ under the bridge matching framework, our key idea is to construct a new probabilistic path $p_{t}$ based on the existing probabilistic path $q_{t}$ from Eq. 3, such that the following condition is satisfied for $t\in[0,1]$:

$$ $p_{t}(\vec{{\mathbf{X}}}_{t})=\frac{1}{Z_{t}}q_{t}(\vec{{\mathbf{X}}}_{t})\exp (-k\varepsilon_{t}(\vec{{\mathbf{X}}}_{t})),~{}\varepsilon_{0}(\cdot)= \varepsilon_{1}(\cdot)=\varepsilon(\cdot).$ (7) $$

Here $Z_{t}$ is the partition function and $\varepsilon_{t}$ is an artificially-designed intermediate potential of the process, which should converge to the real MD potential $\varepsilon$ when $t\to 0^{+}$ and $t\to 1^{-}$ for consistency. Further, we assume $p_{t}(\vec{{\mathbf{X}}}_{t}|\vec{{\bm{X}}}_{0},\vec{{\bm{X}}}_{1})=q_{t}(\vec
{{\mathbf{X}}}_{t}|\vec{{\bm{X}}}_{0},\vec{{\bm{X}}}_{1})$ as in Eq. 6, thereby the stochastic process governed by $p_{t}$ shares the same form with Eq. 3, which can be modeled by the following Markov diffusion process associated with a vector field $v^{\prime}$:

$$ $d\vec{{\mathbf{X}}}_{t}=v^{\prime}(\vec{{\mathbf{X}}}_{t},t)dt+\sigma{d{ \mathbf{B}}_{t}}.$ (8) $$

Therefore, we are able to inference $p_{t}$ if we know how to learn $v^{\prime}$ from the dataset. Interestingly, we prove that $v^{\prime}(\vec{{\mathbf{X}}}_{t},t)$ can be expressed in terms of the vector field $v^{*}(\vec{{\mathbf{X}}}_{t},t)$ generating $q_{t}$ and the intermediate force field $\nabla\varepsilon_{t}(\vec{{\mathbf{X}}}_{t})$ allowing for the Boltzmann constraint, which will be formally demonstrated in Proposition 3.2. In prior to showing this proposition, we first derive the form of the intermediate force field below:

###### Proposition 3.1 .

Assume that the joint distributions $q(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1})$ and $p(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1})$ satisfy $p(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1})\propto{q(\vec{{\mathbf{X}}}_{
0},\vec{{\mathbf{X}}}_{1})\exp(-k(\varepsilon(\vec{{\mathbf{X}}}_{0})+
\varepsilon(\vec{{\mathbf{X}}}_{1})))}$, the intermediate force field $\nabla\varepsilon_{t}$ is given by:

$$ $\displaystyle\nabla\varepsilon_{t}(\vec{{\mathbf{X}}}_{t})$ $\displaystyle=\frac{\mathbb{E}_{q(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1 })}[q_{t}(\vec{{\mathbf{X}}}_{t}|\vec{{\bm{X}}}_{0},\vec{{\bm{X}}}_{1})\exp(-k (\varepsilon(\vec{{\bm{X}}}_{0})+\varepsilon(\vec{{\bm{X}}}_{1})))\zeta(\vec{{ \bm{X}}}_{0},\vec{{\bm{X}}}_{1},\vec{{\mathbf{X}}}_{t})]}{k\mathbb{E}_{q(\vec{ {\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1})}[q_{t}(\vec{{\mathbf{X}}}_{t}|\vec{{ \bm{X}}}_{0},\vec{{\bm{X}}}_{1})\exp(-k(\varepsilon(\vec{{\bm{X}}}_{0})+ \varepsilon(\vec{{\bm{X}}}_{1})))]},$ (9) $$

where we denote $\zeta(\vec{{\bm{X}}}_{0},\vec{{\bm{X}}}_{1},\vec{{\mathbf{X}}}_{t})=\nabla\log
{q_{t}}(\vec{{\mathbf{X}}}_{t})-\nabla\log{q_{t}}(\vec{{\mathbf{X}}}_{t}|\vec{
{\bm{X}}}_{0},\vec{{\bm{X}}}_{1})$ for brevity.

###### Proposition 3.2 .

Given the prerequisites in Proposition 3.1, we have $v^{\prime}(\vec{{\mathbf{X}}}_{t},t)=v^{*}(\vec{{\mathbf{X}}}_{t},t)-\frac{
\sigma^{2}}{2}k\nabla\varepsilon_{t}(\vec{{\mathbf{X}}}_{t})$ under some mild assumptions.

###### Proposition 3.1 .

###### Proposition 3.2 .

#### Estimation of Immediate Force Field

However, it is challenging to calculate the intermediate force field $\nabla\varepsilon_{t}(\vec{{\mathbf{X}}}_{t})$ owing to its nontrivial form in Eq. 9. Note that $q_{t}(\vec{{\mathbf{X}}}_{t}|\vec{{\bm{X}}}_{0},\vec{{\bm{X}}}_{1})$ can be directly calculated by Eq. 6, $\varepsilon(\vec{{\bm{X}}}_{0}),\varepsilon(\vec{{\bm{X}}}_{1})$ are both known as MD potentials, and the the expectation in the denominator can be estimated with samples during a training mini-batch instead of the entire data distribution (Lu et al., ). The most challenging part is the computation of the term $\nabla\log{q_{t}}(\vec{{\mathbf{X}}}_{t})$ in $\zeta(\vec{{\bm{X}}}_{0},\vec{{\bm{X}}}_{1},\vec{{\mathbf{X}}}_{t})$. To provide an unbiased estimation of this score, we investigate its relation with the vector fields.
Given ${\mathcal{G}}_{0},{\mathcal{G}}_{1}$, we first define the conditional score as $s_{t}(\vec{{\mathbf{X}}}_{t}|\vec{{\bm{X}}}_{0},\vec{{\bm{X}}}_{1})=\nabla\log
{q_{t}(\vec{{\mathbf{X}}}_{t}|\vec{{\bm{X}}}_{0},\vec{{\bm{X}}}_{1})}$. Based on Eq. 6, the closed-form of $s_{t}$ is given by:

$$ $s_{t}(\vec{{\mathbf{X}}}_{t}|\vec{{\bm{X}}}_{0},\vec{{\bm{X}}}_{1})=-\frac{ \vec{{\mathbf{X}}}_{t}-[t\vec{{\bm{X}}}_{1}+(1-t)\vec{{\bm{X}}}_{0}]}{t(1-t) \sigma^{2}}=\frac{1}{\sigma^{2}}[\frac{\vec{{\bm{X}}}_{1}-\vec{{\mathbf{X}}}_{ t}}{1-t}-\frac{\vec{{\mathbf{X}}}_{t}-\vec{{\bm{X}}}_{0}}{t}].$ (10) $$

Note that the first term of Eq. 10 has the same form as in the training objective of Eq. 5. Similarly, we train another network $u_{\theta}$ to imitate the second term:

$$ ${\mathcal{L}}_{\mathrm{rev}}=\mathbb{E}_{t\sim\mathrm{Uni}(0,1),({\mathcal{G}} _{0},{\mathcal{G}}_{1})\sim{\mathcal{D}},\vec{{\bm{X}}}_{t}\sim{q_{t}(\cdot| \vec{{\bm{X}}}_{0},\vec{{\bm{X}}}_{1})}}[||\frac{\vec{{\bm{X}}}_{t}-\vec{{\bm{ X}}}_{0}}{t}-u_{\theta}(\vec{{\bm{X}}}_{t},t)||^{2}],$ (11) $$

where the expectation of $u_{\theta}$ can be expressed as $u^{*}(\vec{{\mathbf{X}}}_{t},t)=\mathbb{E}_{q_{t}(\cdot,\cdot|\vec{{\mathbf{X}
}}_{t})}[\frac{\vec{{\mathbf{X}}}_{t}-\vec{{\bm{X}}}_{0}}{t}]$. Then we take the expectation over $\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1}$ conditioned on $\vec{{\mathbf{X}}}_{t}$ in Eq. 10, yielding:

$$ $s_{t}^{*}(\vec{{\mathbf{X}}}_{t})=\mathbb{E}_{q_{t}(\cdot,\cdot|\vec{{\mathbf{ X}}}_{t})}[s_{t}(\vec{{\mathbf{X}}}_{t}|\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X }}}_{1})]=\frac{v^{*}(\vec{{\mathbf{X}}}_{t},t)-u^{*}(\vec{{\mathbf{X}}}_{t},t )}{\sigma^{2}}.$ (12) $$

We present Proposition 3.3 to reveal that $s_{t}^{*}$ is identical to the marginal score $\nabla\log{q_{t}}$ of interest:

###### Proposition 3.3 .

We have $\nabla\log{q_{t}}(\vec{{\mathbf{X}}}_{t})=s_{t}^{*}(\vec{{\mathbf{X}}}_{t})$, where $\nabla\log{q_{t}}$ is the score of the Brownian bridge defined in Eq. 3 and $s_{t}^{*}(\vec{{\mathbf{X}}}_{t})$ is the expectation of the conditional score given by Eq. 12.

In practice, $\nabla\log{q_{t}}(\vec{{\mathbf{X}}}_{t})$ is estimated as $\nicefrac{{(v_{\theta}(\vec{{\mathbf{X}}}_{t},t)-u_{\theta}(\vec{{\mathbf{X}}}
_{t},t))}}{{\sigma^{2}}}$, by replacing the vector fields with the learned neural networks $v_{\theta},u_{\theta}$ in Eq. 12.
Now, since all quantities related to the intermediate force field $\nabla\varepsilon_{t}(\vec{{\mathbf{X}}}_{t})$ are commutable, we then train a neural network $w_{\theta}(\vec{{\mathbf{X}}}_{t},t)$ as its unbiased estimator:

$$ ${\mathcal{L}}_{\mathrm{iff}}=\mathbb{E}_{t,({\mathcal{G}}_{0},{\mathcal{G}}_{1 }),q_{t}(\cdot|\vec{{\bm{X}}}_{0},\vec{{\bm{X}}}_{1})}[||\frac{\exp(-k( \varepsilon(\vec{{\bm{X}}}_{0})+\varepsilon(\vec{{\bm{X}}}_{1})))\zeta(\vec{{ \bm{X}}}_{0},\vec{{\bm{X}}}_{1},\vec{{\bm{X}}}_{t})}{k\mathbb{E}_{q(\vec{{\bm{ X}}}_{0},\vec{{\bm{X}}}_{1})^{B}}[q_{t}(\vec{{\bm{X}}}_{t}|\vec{{\bm{X}}}_{0}, \vec{{\bm{X}}}_{1})\exp(-k(\varepsilon(\vec{{\bm{X}}}_{0})+\varepsilon(\vec{{ \bm{X}}}_{1})))]}-w_{\theta}(\vec{{\bm{X}}}_{t},t)||^{2}],$ (13) $$

where $B$ denotes the mini-batch size of each training step.

###### Proposition 3.3 .

### 3.4 Full Training Processes and Inference

In addition to the aforementioned regression losses, we introduce an auxiliary loss from Yim et al. (), which promotes predictions of pairwise atomic relations. The loss is defined as:

$$ ${\mathcal{L}}_{\mathrm{aux}}=(1-t)\cdot\frac{||\bm{1}_{{\bm{D}}_{0}<6\mathring {\mathrm{A}}}({\bm{D}}_{0}-\hat{{\bm{D}}}_{0})||^{2}}{\sum{\bm{1}_{{\bm{D}}_{0 }<6\mathring{\mathrm{A}}}}-N}+t\cdot\frac{||\bm{1}_{{\bm{D}}_{1}<6\mathring{ \mathrm{A}}}({\bm{D}}_{1}-\hat{{\bm{D}}}_{1})||^{2}}{\sum{\bm{1}_{{\bm{D}}_{1} <6\mathring{\mathrm{A}}}}-N},$ (14) $$

where $D_{0},D_{1}\in{\mathbb{R}}^{N\times{N}}$ denote pairwise distances between all atoms of ${\mathcal{G}}_{0}$ and ${\mathcal{G}}_{1}$, and $\hat{D}_{0},\hat{D}_{1}$ are defined in the same way based on the estimated starting point $\hat{\vec{{\bm{X}}}}_{0}=\vec{{\bm{X}}}_{t}-tu_{\theta}(\vec{{\bm{X}}}_{t},t)$ and terminal point $\hat{\vec{{\bm{X}}}}_{1}=\vec{{\bm{X}}}_{t}+(1-t)v_{\theta}(\vec{{\bm{X}}}_{t}
,t)$. The full loss of FBM-BASE is given by:

$$ ${\mathcal{L}}_{\mathrm{base}}={\mathcal{L}}_{\mathrm{fwd}}+{\mathcal{L}}_{ \mathrm{rev}}+\lambda_{\mathrm{aux}}\cdot{\mathcal{L}}_{\mathrm{aux}},$ (15) $$

where $\lambda_{\mathrm{aux}}$ is a hyper-parameter to balance the weight of different training objectives.

Empirically, large variances are noticed during training FBM with Eq. 13 when $t$ is close to 0 and 1. To address the issue, we find that the intermediate force field converges to the MD force field $\nabla\varepsilon$ at $t=0,1$, which is guaranteed by Prop. 3.4:

###### Proposition 3.4 .

Given $\varepsilon_{0}=\varepsilon_{1}=\varepsilon$ and the intermediate force field described in Eq. 9, the continuity condition $\lim_{t\to{0^{+}}}\nabla\varepsilon_{t}(\vec{{\mathbf{X}}}_{t})=\nabla
\varepsilon(\vec{{\mathbf{X}}}_{0}),~{}\lim_{t\to{1^{-}}}\nabla\varepsilon_{t}
(\vec{{\mathbf{X}}}_{t})=\nabla\varepsilon(\vec{{\mathbf{X}}}_{1})$ holds.

Therefore, we leverage two separate networks $w^{(1)}_{\theta},w^{(2)}_{\theta}$ to learn the boundary force fields:

$$ ${\mathcal{L}}_{\mathrm{bnd}}=\mathbb{E}_{t,({\mathcal{G}}_{0},{\mathcal{G}}_{1 }),q_{t}(\cdot|\vec{{\bm{X}}}_{0},\vec{{\bm{X}}}_{1})}[||\nabla\varepsilon( \vec{{\bm{X}}}_{0})-w^{(1)}_{\theta}(\vec{{\bm{X}}}_{t},t)||^{2}+||\nabla \varepsilon(\vec{{\bm{X}}}_{1})-w^{(2)}_{\theta}(\vec{{\bm{X}}}_{t},t)||^{2}].$ (16) $$

We construct the network $w_{\theta}$ in the interpolation form with another network $w^{(3)}_{\theta}$, similar to Máté & Fleuret (): $w_{\theta}(\vec{{\mathbf{X}}}_{t},t)=(1-t)w^{(1)}_{\theta}(\vec{{\mathbf{X}}}_
{t},t)+tw^{(2)}_{\theta}(\vec{{\mathbf{X}}}_{t},t)+t(1-t)w^{(3)}_{\theta}(\vec
{{\mathbf{X}}}_{t},t)$. The ultimate loss for training FBM is given by:

$$ ${\mathcal{L}}_{\mathrm{FBM}}={\mathcal{L}}_{\mathrm{iff}}+{\mathcal{L}}_{ \mathrm{bnd}}.$ (17) $$

###### Proposition 3.4 .

#### Full Training Processes

We first perform training to derive $v_{\theta}$ and $u_{\theta}$ under the base loss in Eq. 15, and then continue the training process to attain $w_{\theta}$ under the FBM loss in Eq. 17.
The pseudo codes for training FBM-BASE and FBM are in Alg. 1 and Alg. 2, respectively.

#### Force-guided Inference

For inference with FBM, we estimate the vector field $v^{\prime}$ based on Proposition 3.2, by replacing $v^{*}$ and $\nabla\varepsilon_{t}$ with the trained neural networks $v_{\theta}$ and $w_{\theta}$. Force-guided inference is then performed following the SDE process of Eq. 8, where the diffusion time $t$ is discretized equidistantly. Additional details and pseudo codes for inference are provided in § D.

## 4 Experiments

#### Dataset Generation

We evaluate FBM on two datasets consisting of small peptides: *Alanine-Dipeptide (AD)* that is commonly studied previously and *PepMD* which is created by us.
The AD dataset contains a simple peptide with only 22 atoms. The initial structure and reference MD trajectories of AD are all obtained from mdshare(^1^11https://github.com/markovmodel/mdshare) without post-processing.
As for PepMD, we first screen valid peptides between 3-10 residues from the sequence data provided by PDB (Berman et al., ). Next, we perform data cleaning according to the following criterion: each peptide must contain only the 20 natural amino acids, and the number of any type of residue should not exceed 50% of the sequence length. We then cluster the data with a sequence identity threshold of 60% by MMseq2 (Steinegger & Söding, ), and randomly sample one peptide from each cluster to obtain a non-redundant dataset. Considering the computing resource constraints, we select 136/14 peptides for constructing the training-validation/test set, respectively. The structures of all 150 peptides are predicted by open-source tools RDKit () and PDBfixer(^2^22https://github.com/openmm/pdbfixer), which are sent as initial states to generate MD trajectories using OpenMM (Eastman et al., ) afterwards. Finally, the peptide pairs for training are then sampled from trajectories in the way depicted in § 3.1. The MD simulation setups and the statistical details of our curated dataset are both illustrated in § E.1.

#### Baselines

We compare our FBM with the following generative models that learns time-coarsened dynamics: (i) Timewarp (Klein et al., ), the current state-of-the-art model targeting the Boltzmann distribution by MCMC resampling, which exhibits superior transferability to unseen peptide systems. (ii) ITO (Schreiner et al., ), a conditional diffusion model that learns multiple time-resolution dynamics. (iii) Score Dynamics (Hsu et al., ), a score-based diffusion model that learns discrete transitions of the dynamic variables. All models are trained on PepMD from scratch for fair comparison.

#### Metrics

Following Wang et al. (), we evaluate generated conformation ensembles against the full MD trajectories as to their validity, flexibility, and distributional similarity. We provide brief descriptions of the metrics in this part and further details are illustrated in § E.2:

- •
Validity. We regard a molecular conformation as valid when it is governed by certain physics constraints. Following Lu et al. (), we judge whether a conformation is valid by the criterion: no bond clashes between any residue pairs and no bond breaks between adjacent residues, based on coordinates of $\alpha$-carbons. This metric, named as VAL-CA, represents the fraction of valid conformations in the full generated conformation ensembles.
- •
Flexibility. The generated structures are further required to exhibit flexibility to capture dynamic characteristics. Following Janson et al. (), we report the root mean square error of contact maps between generated conformation ensembles and reference MD trajectories as a measure of flexibility, termed as CONTACT.
- •
Distributional similarity. We focus on the similarity between the sample distribution and the Boltzmann distribution. Instead of the costly computation of the Boltzmann density, we project the molecular conformations onto following low-dimensional feature spaces and calculate the Jensen-Shannon (JS) distance as a substitute (Lu et al., ): (i) pairwise distances between $\alpha$-carbons of residues (PWD); (ii) radius-of-gyration (RG) that measures the distribution of $\alpha$-carbons to the center-of-mass; (iii) the time-lagged independent components (Pérez-Hernández et al., ) (TIC) based on dihedrals and pairwise distances of $\alpha$-carbons, where only the slowest components, TIC 0 and TIC 1, are taken into consideration. For each metric, the mean JS distance along all feature dimensions are reported.

### 4.1 Metastable States Exploration for AD

We first investigate how well the generated conformations can travel across different metastable states of AD. Due to the simple structure of AD with only one peptide bond, some metrics are not applicable except for TIC and TIC-2D (*i.e.*, the joint distribution of TIC 0 and TIC 1). In particular, the backbone dihedrals of AD, psi and phi, are commonly considered as two challenging variables for state transitions. Therefore, we include the similarity measurement of the joint distribution of psi and phi, *i.e.* the Ramachandran plot (Ramachandran et al., ), denoted as RAM.

In Table 1 we show evaluation results on AD, where models sample in the time-coarsened manner from the same initial state for a chain length of $10^{3}$. According to Table 1, FBM outperforms existing baselines on both RAM and TIC metrics, and with the introduction of physics priors, it shows considerable improvements in distribution similarity across various feature spaces compared to FBM-BASE. Although Timewarp surpasses FBM in the TIC-2D metric, we will explain later that it comes at the cost of generating invalid conformations.

Further, Ramachandran plots of generated ensembles are illustrated in Figure 3, where three known metastable states are recognized based on MD trajectories and labeled in order. Apparently, ITO and Score Dynamics fail to capture the dynamics of AD with samples randomly allocated. Moreover, Timewarp cannot rapidly traverse through metastable states, resulting in a great portion of invalid samples. Despite both FBM-BASE and FBM show relatively “clean” plots with fewer unreasonable conformations, FBM pays more attention to high density regions including cluster 2 and cluster 3, confirming a strong guidance of the intermediate force field to align more closely with thermodynamic principles.

**Table 1: Results on alanine dipeptide. Values of each metric are averaged over three independent runs. The best result for each metric is shown in bold and the second best is \ulunderlined.**
| MODELS | JS DISTANCE ($\downarrow$) |  |  |
| --- | --- | --- | --- |
| RAM | TIC | TIC-2D |  |
| TIMEWARP | \ul0.722 | 0.546 | 0.719 |
| ITO | 0.740 | 0.696 | 0.833 |
| SD | 0.731 | 0.673 | 0.807 |
| FBM-BASE | 0.727 | \ul0.533 | 0.749 |
| FBM | 0.711 | 0.525 | \ul0.733 |

Figure: Figure 3: Ramachandran plots of alanine dipeptide with conformation ensembles generated by models. The initial state is indicated with the red cross. Contours represent the kernel densities estimated by the MD trajectory and the generated conformations are shown in scatter.
Refer to caption: x3.png

### 4.2 Transferability to Unseen Peptides of PepMD

We then explore the transferability of models to unseen peptides with various sequence lengths of PepMD. We use all metrics in § 4 for evaluation and the results on 14 test peptides of PepMD are demonstrated in Table 2, where all samples are generated for a chain length of $10^{3}$. We find Timewarp achieves a good performance on distributional similarity and flexibility, yet at the cost of only a small portion of valid samples. In contrast, our FBM showcases superiority across all metrics and achieves significant improvement on the validity of generated conformations in particular. The results indicate that by introducing the force guidance, generated ensembles of FBM are more akin to the underlying Boltzmann distribution and better comply with physical constraints. Additional experimental results can be found in § F.

**Table 2: Results on the test set of PepMD. Values of each metric are first averaged over 3 independent runs for each peptide and then shown in mean/std of all 14 test peptides. The best result for each metric is shown in bold and the second best is \ulunderlined.**
| MODELS | JS DISTANCE ($\downarrow$) | VAL-CA ($\uparrow$) | CONTACT ($\downarrow$) |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| PWD | RG | TIC | TIC-2D |  |  |  |
| TIMEWARP | \ul0.575/0.082 | 0.561/0.124 | \ul0.633/0.069 | \ul0.804/0.025 | 0.115/0.121 | \ul0.197/0.128 |
| ITO | 0.833/0.000 | 0.829/0.012 | 0.789/0.067 | 0.833/0.000 | 0.001/0.000 | 0.940/0.081 |
| SD | 0.823/0.030 | 0.818/0.041 | 0.773/0.032 | 0.832/0.001 | 0.006/0.016 | 0.824/0.095 |
| FBM-BASE | 0.576/0.066 | \ul0.560/0.153 | 0.639/0.061 | 0.807/0.020 | \ul0.367/0.173 | 0.208/0.142 |
| FBM | 0.573/0.064 | 0.542/0.140 | 0.631/0.077 | 0.801/0.032 | 0.616/0.188 | 0.188/0.127 |

In Figure 4 we provide TIC-2D results of the test peptide 1e28:C (*i.e.*, chain C of PDB 1e28) as a case, which is composed of 8 residues and reveals 3 major metastable states by MD simulations. For such complex peptide, ITO and Score Dynamics even hardly generate conformations close to any high density regions. Meanwhile, a large number of samples generated by Timewarp and FBM-BASE fall into low density regions and moreover, Timewarp scarcely accesses region 2, signaling a lack of validity and transferability. In contrast to FBM-BASE, samples generated by FBM are more concentrated in the highest density region 3 with relatively correct weights, and also consist of fewer invalid conformations that get trapped in low density regions.

Figure: Figure 4: TIC plots of peptide 1e28:C with conformation ensembles generated by models.
Refer to caption: x8.png

Furthermore, we explore the ability of FBM to recover different equilibrium structures, which is measured by the lowest sample Root-Mean-Square-Deviation (RMSD) to each cluster center (Wang et al., ). To be specific, we first select the reference structure with the lowest kernel density to each cluster center from MD trajectories, then based on the measurement of all-atom RMSD, we search for the most similar samples generated by FBM to each of the reference structures. The visualization results are demonstrated in Figure 5, where FBM shows a near-exact recovery of known, representative conformations. As for the generation efficiency according to visiting certain metastable states, FBM achieves 5-30x faster inference time than that of MD simulations.

Figure: Figure 5: Comparison between the reference equilibrium structures (blue) and the generated structures of FBM (yellow) with the lowest all-atom RMSD to each cluster center of peptide 1e28:C. Clusters are arranged in order based on the time of reaching the corresponding metastable state in MD simulations. Both MD simulations and FBM inference start with the same initial state (grey).
Refer to caption: x13.png

## 5 Conclusion and Future Work

In this work, we present a novel generative model FBM for time-coarsened dynamics in a full-atom fashion. We first leverage the bridge matching framework to construct the baseline model FBM-BASE for learning dynamics from the data distribution. Based on FBM-BASE, we further introduce physics priors and interpolate a well-designed intermediate force field, which is theoretically guaranteed to target the Boltzmann-constrained distribution that aligns with thermodynamic principles. Experiments on alanine dipeptide and our curated dataset PepMD showcase superiority of FBM on comprehensive metrics and demonstrate transferability to unseen peptide systems.

As the first attempt to incorporate the intermediate force field to bridge matching for full-atom time-coarsened dynamics, our method has considerable room for improvement. Firstly, our experiments have been conducted on small peptides with fewer than 10 residues. Further exploration on more complex molecular systems (*e.g.*, proteins) is warranted. Secondly, since the training labels for FBM depend on the marginal score calculations provided by FBM-BASE, we had to adopt a two-stage training process rather than an end-to-end one, which increases the training overhead. Lastly, Finally, the transitions between metastable states are still not fast enough, resulting in unreasonable conformations along the generated paths. Therefore, methods for rapid and jump-like state transition are of great importance.

## Appendix

## Appendix A Reproducibility

Our source code and the curated dataset PepMD are all available at [https://github.com/yaledeus/FBM](https://github.com/yaledeus/FBM).

## Appendix B Proofs of Propositions

Proof of Proposition 3.1. Given the joint distributions of random variables $\vec{{\mathbf{X}}}_{0}$ and $\vec{{\mathbf{X}}}_{1}$ following densities $q_{0},q_{1}$ and $p_{0},p_{1}$ satisfy:

$$ $p(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1})=\frac{1}{Z}q(\vec{{\mathbf{X} }}_{0},\vec{{\mathbf{X}}}_{1})\exp(-k(\varepsilon(\vec{{\mathbf{X}}}_{0})+ \varepsilon(\vec{{\mathbf{X}}}_{1}))),$ (18) $$

where $Z$ is the partition function to ensure $\iint{p(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1})}d\vec{{\mathbf{X}}}_{0}
d\vec{{\mathbf{X}}}_{1}=1$. It is easy to verify that the condition holds when $\vec{{\mathbf{X}}}_{0}$ and $\vec{{\mathbf{X}}}_{1}$ are independent variables.

Therefore, the marginal density $p_{t}$ is given by:

$$ $\displaystyle p_{t}(\vec{{\mathbf{X}}}_{t})$ $\displaystyle=\iint{p_{t}(\vec{{\mathbf{X}}}_{t}|\vec{{\mathbf{X}}}_{0},\vec{{ \mathbf{X}}}_{1})p(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1})d\vec{{ \mathbf{X}}}_{0}d\vec{{\mathbf{X}}}_{1}}$ (19) $\displaystyle=\iint{q_{t}(\vec{{\mathbf{X}}}_{t}|\vec{{\mathbf{X}}}_{0},\vec{{ \mathbf{X}}}_{1})q(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1})\frac{\exp(-k (\varepsilon(\vec{{\mathbf{X}}}_{0})+\varepsilon(\vec{{\mathbf{X}}}_{1})))}{Z} d\vec{{\mathbf{X}}}_{0}d\vec{{\mathbf{X}}}_{1}}$ (20) $\displaystyle=\iint{q_{t}(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1}|\vec{{ \mathbf{X}}}_{t})q_{t}(\vec{{\mathbf{X}}}_{t})\frac{\exp(-k(\varepsilon(\vec{{ \mathbf{X}}}_{0})+\varepsilon(\vec{{\mathbf{X}}}_{1})))}{Z}d\vec{{\mathbf{X}}} _{0}d\vec{{\mathbf{X}}}_{1}}$ (21) $\displaystyle=q_{t}(\vec{{\mathbf{X}}}_{t})\mathbb{E}_{q_{t}(\cdot,\cdot|\vec{ {\mathbf{X}}}_{t})}[\frac{\exp(-k(\varepsilon(\vec{{\mathbf{X}}}_{0})+ \varepsilon(\vec{{\mathbf{X}}}_{1})))}{Z}],$ (22) $$

where in the first equality we use the assumption that probability paths $p_{t}$ and $q_{t}$ share the same conditional distribution given $\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1}$. Considering the assumption that $p_{t}$ admits the Boltzmann-constrained form as in Eq. 7, we can easily derive the formula for the intermediate potential:

$$ $\varepsilon_{t}(\vec{{\mathbf{X}}}_{t})=-\frac{1}{k}\log\mathbb{E}_{q_{t}(\vec {{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1}|\vec{{\bm{X}}}_{t})}[\exp(-k( \varepsilon(\vec{{\bm{X}}}_{0})+\varepsilon(\vec{{\bm{X}}}_{1})))]+\frac{1}{k} \log\frac{Z}{Z_{t}}.$ (23) $$

Take the gradient of Eq. 23 with regard to $\vec{{\mathbf{X}}}_{t}$, the intermediate force field is given by:

$$ $\nabla\varepsilon_{t}(\vec{{\mathbf{X}}}_{t})=-\frac{\iint\exp(-k(\varepsilon( \vec{{\mathbf{X}}}_{0})+\varepsilon(\vec{{\mathbf{X}}}_{1})))\nabla{q_{t}(\vec {{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1}|\vec{{\mathbf{X}}}_{t})}d\vec{{ \mathbf{X}}}_{0}d\vec{{\mathbf{X}}}_{1}}{k\mathbb{E}_{q_{t}(\cdot,\cdot|\vec{{ \mathbf{X}}}_{t})}[\exp(-k(\varepsilon(\vec{{\mathbf{X}}}_{0})+\varepsilon( \vec{{\mathbf{X}}}_{1})))]},$ (24) $$

where we assume that the integrals and gradients can be commuted. The numerator of Eq. 24 can be further expanded by:

$$ $\displaystyle\iint\exp(-k(\varepsilon(\vec{{\mathbf{X}}}_{0})+\varepsilon(\vec {{\mathbf{X}}}_{1})))\nabla{q_{t}(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1 }|\vec{{\mathbf{X}}}_{t})}d\vec{{\mathbf{X}}}_{0}d\vec{{\mathbf{X}}}_{1}$ (25) $\displaystyle=\iint\exp(-k(\varepsilon(\vec{{\mathbf{X}}}_{0})+\varepsilon( \vec{{\mathbf{X}}}_{1})))q_{t}(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1}| \vec{{\mathbf{X}}}_{t})\nabla\log{q_{t}(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X }}}_{1}|\vec{{\mathbf{X}}}_{t})d\vec{{\mathbf{X}}}_{0}d\vec{{\mathbf{X}}}_{1}}$ (26) $\displaystyle=\iint\exp(-k(\varepsilon(\vec{{\mathbf{X}}}_{0})+\varepsilon( \vec{{\mathbf{X}}}_{1})))q_{t}(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1}| \vec{{\mathbf{X}}}_{t})\nabla\log\frac{q_{t}(\vec{{\mathbf{X}}}_{t}|\vec{{ \mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1})q(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf {X}}}_{1})}{q_{t}(\vec{{\mathbf{X}}}_{t})}d\vec{{\mathbf{X}}}_{0}d\vec{{ \mathbf{X}}}_{1}$ (27) $\displaystyle=\mathbb{E}_{q_{t}(\cdot,\cdot|\vec{{\mathbf{X}}}_{t})}[\exp(-k( \varepsilon(\vec{{\mathbf{X}}}_{0})+\varepsilon(\vec{{\mathbf{X}}}_{1})))( \nabla\log{q_{t}(\vec{{\mathbf{X}}}_{t}|\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X }}}_{1})})-\nabla\log{q_{t}(\vec{{\mathbf{X}}}_{t})}].$ (28) $$

Substitute the numerator back into Eq. 24 and the conclusion holds.∎

Proof of Proposition 3.2. Similar to the definition of $v^{*}(\vec{{\mathbf{X}}}_{t},t)$ and $u^{*}(\vec{{\mathbf{X}}}_{t},t)$ under the probability path $q_{t}$, we define $v^{\prime}(\vec{{\mathbf{X}}}_{t},t)$ and $u^{\prime}(\vec{{\mathbf{X}}}_{t},t)$ under the probability path $p_{t}$ accordingly, which have the following form:

$$ $v^{\prime}(\vec{{\mathbf{X}}}_{t},t)=\mathbb{E}_{p_{t}(\cdot,\cdot|\vec{{ \mathbf{X}}}_{t})}[\frac{\vec{{\mathbf{X}}}_{1}-\vec{{\mathbf{X}}}_{t}}{1-t}], ~{}u^{\prime}(\vec{{\mathbf{X}}}_{t},t)=\mathbb{E}_{p_{t}(\cdot,\cdot|\vec{{ \mathbf{X}}}_{t})}[\frac{\vec{{\mathbf{X}}}_{t}-\vec{{\mathbf{X}}}_{0}}{t}].$ (29) $$

Further, the marginal scores of probability paths $q_{t},p_{t}$ at diffusion time $t$ (termed as $s_{t}^{*}$ and $s_{t}^{\prime}$ respectively) are connected based on Eq. 7 as follows:

$$ $s_{t}^{\prime}(\vec{{\mathbf{X}}}_{t})=s_{t}^{*}(\vec{{\mathbf{X}}}_{t})-k \nabla\varepsilon_{t}(\vec{{\mathbf{X}}}_{t}).$ (30) $$

Considering the linkage between scores and vector fields in Eq. 12, we obtain:

$$ $\displaystyle-k\nabla\varepsilon_{t}(\vec{{\mathbf{X}}}_{t})$ $\displaystyle=\frac{v^{\prime}(\vec{{\mathbf{X}}}_{t},t)-u^{\prime}(\vec{{ \mathbf{X}}}_{t},t)}{\sigma^{2}}-\frac{v^{*}(\vec{{\mathbf{X}}}_{t},t)-u^{*}( \vec{{\mathbf{X}}}_{t},t)}{\sigma^{2}}$ (31) $\displaystyle=\frac{1}{\sigma^{2}}[(v^{\prime}(\vec{{\mathbf{X}}}_{t},t)-v^{*} (\vec{{\mathbf{X}}}_{t},t))-(u^{\prime}(\vec{{\mathbf{X}}}_{t},t)-u^{*}(\vec{{ \mathbf{X}}}_{t},t))].$ (32) $$

Now we investigate the intrinsic properties of $v^{\prime}-v^{*}$ by expanding the term:

$$ $v^{\prime}(\vec{{\mathbf{X}}}_{t},t)-v^{*}(\vec{{\mathbf{X}}}_{t},t)=\frac{1}{ 1-t}\iint\vec{{\mathbf{X}}}_{1}\left(p_{t}(\vec{{\mathbf{X}}}_{0},\vec{{ \mathbf{X}}}_{1}|\vec{{\mathbf{X}}}_{t})-q_{t}(\vec{{\mathbf{X}}}_{0},\vec{{ \mathbf{X}}}_{1}|\vec{{\mathbf{X}}}_{t})\right)d\vec{{\mathbf{X}}}_{0}d\vec{{ \mathbf{X}}}_{1}.$ (33) $$

For convenience, we define:

$$ $\displaystyle f(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1},\vec{{\mathbf{X} }}_{t},t)$ $\displaystyle=p_{t}(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1}|\vec{{ \mathbf{X}}}_{t})-q_{t}(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1}|\vec{{ \mathbf{X}}}_{t})$ (34) $\displaystyle=\frac{p_{t}(\vec{{\mathbf{X}}}_{t}|\vec{{\mathbf{X}}}_{0},\vec{{ \mathbf{X}}}_{1})p(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1})}{p_{t}(\vec{ {\mathbf{X}}}_{t})}-\frac{q_{t}(\vec{{\mathbf{X}}}_{t}|\vec{{\mathbf{X}}}_{0}, \vec{{\mathbf{X}}}_{1})q(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1})}{q_{t} (\vec{{\mathbf{X}}}_{t})}$ (35) $\displaystyle=\frac{q_{t}(\vec{{\mathbf{X}}}_{t}|\vec{{\mathbf{X}}}_{0},\vec{{ \mathbf{X}}}_{1})q(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1})}{q_{t}(\vec{ {\mathbf{X}}}_{t})}[\frac{Z_{t}}{Z}\exp(-k(\varepsilon(\vec{{\mathbf{X}}}_{0}) +\varepsilon(\vec{{\mathbf{X}}}_{1})-\varepsilon_{t}(\vec{{\mathbf{X}}}_{t}))) -1].$ (36) $$

The two terms of Eq. 34 share the same form, where each term can be explained as the probability of states $(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1})$ to be the starting and terminal point respectively under probability path $p_{t}$ (*resp.* $q_{t}$), given the state $\vec{{\mathbf{X}}}_{t}$ at diffusion time $t$. Considering the symmetric role of the reverse-time process of the Brownian bridge (Anderson, ) and given the same state $\vec{{\mathbf{X}}}_{t}$ yet at diffusion time $1-t$, states $(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1})$ are now more likely to be the terminal and staring point, respectively. Formally, we assume $p_{t}(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1}|\vec{{\mathbf{X}}}_{t})=p_
{1-t}(\vec{{\mathbf{X}}}_{1},\vec{{\mathbf{X}}}_{0}|\vec{{\mathbf{X}}}_{t})$ and $q_{t}(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1}|\vec{{\mathbf{X}}}_{t})=q_
{1-t}(\vec{{\mathbf{X}}}_{1},\vec{{\mathbf{X}}}_{0}|\vec{{\mathbf{X}}}_{t})$, then the following equation holds based on Eq. 34:

$$ $f(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1},\vec{{\mathbf{X}}}_{t},t)=f( \vec{{\mathbf{X}}}_{1},\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{t},1-t).$ (37) $$

Denote $g(\vec{{\mathbf{X}}},\vec{{\mathbf{X}}}_{t},t)=\int{f(\vec{{\mathbf{X}}}_{0},
\vec{{\mathbf{X}}},\vec{{\mathbf{X}}}_{t},t)d\vec{{\mathbf{X}}}_{0}}$, Eq. 33 can be rewritten by first integrating over $\vec{{\mathbf{X}}}_{0}$:

$$ $v^{\prime}(\vec{{\mathbf{X}}}_{t},t)-v^{*}(\vec{{\mathbf{X}}}_{t},t)=\frac{1}{ 1-t}\int\vec{{\mathbf{X}}}g(\vec{{\mathbf{X}}},\vec{{\mathbf{X}}}_{t},t)d\vec{ {\mathbf{X}}}.$ (38) $$

Symmetrically, we can define $h(\vec{{\mathbf{X}}},\vec{{\mathbf{X}}}_{t},t)=\int{f(\vec{{\mathbf{X}}},\vec{
{\mathbf{X}}}_{1},\vec{{\mathbf{X}}}_{t},t)d\vec{{\mathbf{X}}}_{1}}$, while $u^{\prime}-u^{*}$ admits a similar form of Eq. 38, given by:

$$ $u^{\prime}(\vec{{\mathbf{X}}}_{t},t)-u^{*}(\vec{{\mathbf{X}}}_{t},t)=-\frac{1} {t}\int\vec{{\mathbf{X}}}h(\vec{{\mathbf{X}}},\vec{{\mathbf{X}}}_{t},t)d\vec{{ \mathbf{X}}}.$ (39) $$

Taking the assumption of Eq. 37 into consideration, an important observation is:

$$ $h(\vec{{\mathbf{X}}},\vec{{\mathbf{X}}}_{t},t)=g(\vec{{\mathbf{X}}},\vec{{ \mathbf{X}}}_{t},1-t).$ (40) $$

On the other hand, when $t$ approaches 1, the limit of function $g$ is given by:

$$ $\displaystyle\lim_{t\to 1^{-}}g(\vec{{\mathbf{X}}},\vec{{\mathbf{X}}}_{t},t)$ (41) $\displaystyle=\lim_{t\to 1^{-}}\frac{1}{q_{t}(\vec{{\mathbf{X}}}_{t})}\int \delta(\vec{{\mathbf{X}}}_{t}-\vec{{\mathbf{X}}})q(\vec{{\mathbf{X}}}_{0},\vec {{\mathbf{X}}})[\frac{Z_{t}}{Z}\exp(-k(\varepsilon(\vec{{\mathbf{X}}}_{0})+ \varepsilon(\vec{{\mathbf{X}}})-\varepsilon_{t}(\vec{{\mathbf{X}}}_{t})))-1]d \vec{{\mathbf{X}}}_{0}$ (42) $\displaystyle=\frac{1}{q_{1}(\vec{{\mathbf{X}}})}\int{q(\vec{{\mathbf{X}}}_{0} ,\vec{{\mathbf{X}}})}[\frac{Z_{1}}{Z}\exp(-k\varepsilon(\vec{{\mathbf{X}}}_{0} ))-1]d\vec{{\mathbf{X}}}_{0}$ (43) $\displaystyle=\frac{1}{q_{1}(\vec{{\mathbf{X}}})}Z_{1}\exp(k\varepsilon(\vec{{ \mathbf{X}}}))\int{\frac{1}{Z}}q(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}) \exp(-k(\varepsilon(\vec{{\mathbf{X}}}_{0})+\varepsilon(\vec{{\mathbf{X}}})))d \vec{{\mathbf{X}}}_{0}-1$ (44) $\displaystyle=\frac{1}{q_{1}(\vec{{\mathbf{X}}})}Z_{1}\exp(k\varepsilon(\vec{{ \mathbf{X}}}))\int{p(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}})}d\vec{{\mathbf {X}}}_{0}-1$ (45) $\displaystyle=\frac{1}{q_{1}(\vec{{\mathbf{X}}})}Z_{1}\exp(k\varepsilon(\vec{{ \mathbf{X}}}))p_{1}(\vec{{\mathbf{X}}})-1=0,$ (46) $$

where in the third and fourth equality we use the relationship between the marginal distribution and the joint distribution, and in the fifth equality we apply Eq. 7.

Here we suppose the function $g$ is separable with respect to the time variable $t$. Formally, there exist functions $\iota$ and $\Gamma$ such that the following identity holds:

$$ $g(\vec{{\mathbf{X}}},\vec{{\mathbf{X}}}_{t},t)\equiv\iota(t)\Gamma(\vec{{ \mathbf{X}}},\vec{{\mathbf{X}}}_{t}).$ (47) $$

Given $\lim_{t\to 1^{-}}g(\vec{{\mathbf{X}}},\vec{{\mathbf{X}}}_{t},t)=0$, we prescribe $\iota(t)=1-t$ for convenience, then we have $h(\vec{{\mathbf{X}}},\vec{{\mathbf{X}}}_{t},t)=t\Gamma(\vec{{\mathbf{X}}},\vec
{{\mathbf{X}}}_{t})$ and subsequently we can derive the closed-form of $v^{\prime}$ based on Eq. 38 and 39 by:

$$ $\displaystyle u^{\prime}(\vec{{\mathbf{X}}}_{t},t)-u^{*}(\vec{{\mathbf{X}}}_{t },t)$ $\displaystyle=-\frac{1}{t}\int\vec{{\mathbf{X}}}h(\vec{{\mathbf{X}}},\vec{{ \mathbf{X}}}_{t},t)d\vec{{\mathbf{X}}}$ (48) $\displaystyle=-\int\vec{{\mathbf{X}}}\Gamma(\vec{{\mathbf{X}}},\vec{{\mathbf{X }}}_{t})d\vec{{\mathbf{X}}}$ (49) $\displaystyle=-\frac{1}{1-t}\int\vec{{\mathbf{X}}}(1-t)\Gamma(\vec{{\mathbf{X} }},\vec{{\mathbf{X}}}_{t})d\vec{{\mathbf{X}}}$ (50) $\displaystyle=-(v^{\prime}(\vec{{\mathbf{X}}}_{t},t)-v^{*}(\vec{{\mathbf{X}}}_ {t},t)),$ (51) $$

which implies $v^{\prime}(\vec{{\mathbf{X}}}_{t},t)=v^{*}(\vec{{\mathbf{X}}}_{t},t)-\frac{
\sigma^{2}}{2}k\nabla\varepsilon_{t}(\vec{{\mathbf{X}}}_{t})$ according to Eq. 31.∎

Proof of Proposition 3.3. From the definition of $s_{t}^{*}$ in Eq. 12, the following equation holds:

$$ $\displaystyle s_{t}^{*}(\vec{{\mathbf{X}}}_{t})$ $\displaystyle=\mathbb{E}_{q_{t}(\cdot,\cdot|\vec{{\mathbf{X}}}_{t})}[\nabla \log{q_{t}(\vec{{\mathbf{X}}}_{t}|\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1 })}]$ (52) $\displaystyle=\iint\nabla\log{q_{t}(\vec{{\mathbf{X}}}_{t}|\vec{{\mathbf{X}}}_ {0},\vec{{\mathbf{X}}}_{1})q_{t}(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1} |\vec{{\mathbf{X}}}_{t})d\vec{{\mathbf{X}}}_{0}d\vec{{\mathbf{X}}}_{1}}$ (53) $\displaystyle=\frac{1}{q_{t}(\vec{{\mathbf{X}}}_{t})}\iint\nabla{q_{t}(\vec{{ \mathbf{X}}}_{t}|\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1})}q(\vec{{ \mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1})d\vec{{\mathbf{X}}}_{0}d\vec{{\mathbf{ X}}}_{1}$ (54) $\displaystyle=\frac{1}{q_{t}(\vec{{\mathbf{X}}}_{t})}\nabla\iint{q_{t}(\vec{{ \mathbf{X}}}_{t}|\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1})}q(\vec{{ \mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1})d\vec{{\mathbf{X}}}_{0}d\vec{{\mathbf{ X}}}_{1}$ (55) $\displaystyle=\frac{\nabla{q_{t}(\vec{{\mathbf{X}}}_{t})}}{q_{t}(\vec{{\mathbf {X}}}_{t})}=\nabla\log{q_{t}(\vec{{\mathbf{X}}}_{t})},$ (56) $$

where in the second equality we use the Bayesian rule of probability, $q(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1})$ denotes the joint distribution of random variables $\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1}$ following $q_{0},q_{1}$. Furthermore, the third equality is justified by assuming the integrands satisfy the regularity conditions of the Leibniz Rule.∎

Proof of Proposition 3.4. We first check whether the continuity condition holds when $t\to 0^{+}$. Note that under this condition, the density $q_{t}(\vec{{\mathbf{X}}}_{0},\vec{{\mathbf{X}}}_{1}|\vec{{\mathbf{X}}}_{t})$ involves into the Dirac mass $\delta(\vec{{\mathbf{X}}}_{t}-\vec{{\mathbf{X}}}_{0})$ at point $\vec{{\mathbf{X}}}_{t}$, subsequently we have:

$$ $\displaystyle\lim_{t\to 0^{+}}\nabla\varepsilon_{t}(\vec{{\mathbf{X}}}_{t})$ $\displaystyle=-\lim_{t\to 0^{+}}\frac{\iint\exp(-k(\varepsilon(\vec{{\mathbf{X }}}_{0})+\varepsilon(\vec{{\mathbf{X}}}_{1})))\delta(\vec{{\mathbf{X}}}_{t}- \vec{{\mathbf{X}}}_{0})d\vec{{\mathbf{X}}}_{0}d\vec{{\mathbf{X}}}_{1}}{k\iint \exp(-k(\varepsilon(\vec{{\mathbf{X}}}_{0})+\varepsilon(\vec{{\mathbf{X}}}_{1} )))\delta(\vec{{\mathbf{X}}}_{t}-\vec{{\mathbf{X}}}_{0})d\vec{{\mathbf{X}}}_{0 }d\vec{{\mathbf{X}}}_{1}}$ (57) $\displaystyle=-\lim_{t\to 0^{+}}\frac{\nabla\iint\exp(-k(\varepsilon(\vec{{ \mathbf{X}}}_{0})+\varepsilon(\vec{{\mathbf{X}}}_{1})))\delta(\vec{{\mathbf{X} }}_{t}-\vec{{\mathbf{X}}}_{0})d\vec{{\mathbf{X}}}_{0}d\vec{{\mathbf{X}}}_{1}}{ k\int\exp(-k(\varepsilon(\vec{{\mathbf{X}}}_{t})+\varepsilon(\vec{{\mathbf{X}} }_{1})))d\vec{{\mathbf{X}}}_{1}}$ (58) $\displaystyle=-\lim_{t\to 0^{+}}\frac{\nabla\int\exp(-k(\varepsilon(\vec{{ \mathbf{X}}}_{t})+\varepsilon(\vec{{\mathbf{X}}}_{1})))d\vec{{\mathbf{X}}}_{1} }{k\int\exp(-k(\varepsilon(\vec{{\mathbf{X}}}_{t})+\varepsilon(\vec{{\mathbf{X }}}_{1})))d\vec{{\mathbf{X}}}_{1}}$ (59) $\displaystyle=-\lim_{t\to 0^{+}}\frac{\nabla\exp(-k\varepsilon(\vec{{\mathbf{X }}}_{t}))}{k\exp(-k\varepsilon(\vec{{\mathbf{X}}}_{t}))}=\nabla\varepsilon( \vec{{\mathbf{X}}}_{0}).$ (60) $$

The case when $t$ approaches 1 is completely symmetrical and will not be elaborated further. Thus we have proven that the intermediate force field converges to the MD force field when $t$ approaches 0 and 1.

## Appendix C Model Architecture

In this work, we leverage the powerful TorchMD-NET (Pelaez et al., ) as the backbone model to process molecular graphs, which intrinsically satisfies $\mathrm{SO}(3)$-equivariance with the equivariant transformer (Thölke & Fabritiis, ) component. To adapt to our task setup, the inputs include not only the Cartesian coordinates $\vec{{\bm{X}}}$ and atom embeddings ${\bm{Z}}$, but also a one-dimensional continuous diffusion time $t$. TorchMD-NET will then output $\mathrm{SO}(3)$-equivariant vectors $\vec{{\bm{V}}}\in{\mathbb{R}}^{N\times{3}\times{H}}$ and node representations
${\bm{H}}\in{\mathbb{R}}^{N\times{H}}$. Formally, we have:

$$ $\vec{{\bm{V}}},{\bm{H}}=\texttt{TorchMD-NET}(\vec{{\bm{X}}},{\bm{Z}},t).$ (61) $$

To streamline the model, we only add lightweight output heads to a single TorchMD-NET module for the baseline model as well as FBM. For the baseline model, we use two separate two-layer Feed-Forward Networks (FFN) with no shared weights to transform the node representations into weights of vectors, which are then multiplied by
$\vec{{\bm{V}}}$ to obtain the final representation:

$$ $\displaystyle v_{\theta}(\vec{{\bm{X}}},t),u_{\theta}(\vec{{\bm{X}}},t) \coloneqq\vec{{\bm{V}}}\times\mathrm{FFN}({\bm{H}})\in{\mathbb{R}}^{N\times{3}},$ (62) $$

where the dimensions of hidden and output layers of $\mathrm{FFN}$ are all $H$ and we use $\mathrm{SiLU}$ (Dong et al., ) for activation layers. Further, we construct the networks $\alpha_{\theta},\beta_{\theta},\gamma_{\theta}$ of FBM in the same way, while the only difference is that we add one $\mathrm{LayerNorm}$ (Ba et al., ) before the $\mathrm{FFN}$ layer due to the variance in scale of different targets.

## Appendix D Training and Inference Details

In this section, we provide additional details and pseudo codes for training and inference of the baseline bridge matching model FBM-BASE and the force-guided bridge matching model FBM.

### D.1 Normalization of Energies and Forces

In practice, we found that the unnormalized potential function is numerically unstable and its variance is positively correlated with the number of atoms $N$. For stable training, we need to perform certain pre-processing steps. Specifically, for potentials in kJ/mol, we divide by $3N$ to obtain the average potential energy per degree of freedom for the entire molecule, where $N$ varies with different peptides. For force fields in kJ/(mol$\cdot$nm), due to their relatively stable values across different molecular systems, we empirically multiply by the constant 0.002 for normalization.

### D.2 Guidance Strength η 𝜂 \eta italic_η

Similar to Wang et al. (), we introduce the guidance strength $\eta$ for better approximation of the Boltzmann distribution. Formally, for any positive constant $\eta>0$, we can define a new probability path based on $q_{t}$:

$$ $p_{t}(\vec{{\mathbf{X}}}_{t})=\frac{1}{Z_{t}}q_{t}(\vec{{\mathbf{X}}}_{t})\exp (-\frac{2\eta}{\sigma^{2}}k\varepsilon_{t}(\vec{{\mathbf{X}}}_{t})), \varepsilon_{0}=\varepsilon_{1}=\varepsilon.$ (63) $$

The only difference with Eq. 7 is the constant $\nicefrac{{2\eta}}{{\sigma^{2}}}$ in the exponential term, which can be interpreted as how well the probability path $p_{t}$ is guided by energies and forces. According to Prop. 3.2, it can be easily deduced that the vector field $v^{\prime}(\vec{{\mathbf{X}}}_{t},t)$ which generates $p_{t}$ has the following form:

$$ $v^{\prime}(\vec{{\mathbf{X}}}_{t},t)=v^{*}(\vec{{\mathbf{X}}}_{t},t)-\eta\cdot {k\varepsilon_{t}(\vec{{\mathbf{X}}}_{t})}.$ (64) $$

Thus practically, we regard $\eta$ as a hyperparameter during inference and enhance the similarity between $p_{t}$ and the Boltzmann distribution by selecting the proper guidance strength $\eta$.

### D.3 Refinement with Constrained Energy Minimization

We utilize the discrete form of the SDE process in Eq. 4 for inference with $T$ SDE steps, and the full conformation ensembles are generated in an autoregressive way, where the output from the previous step serves as the input for the next step. However, in the autoregressive fashion, errors at each inference step will accumulate, leading to out-of-distribution problem. Here we introduce an additional energy minimization procedure using OpenMM (Eastman et al., ) for refinement, which is performed for each generated conformation before sent to the next inference step. Note that we aim for the refinement to affect only the minor details (*e.g.*, X-H bonds) without altering the overall conformation; therefore, independent harmonic constraints are further applied on all heavy atoms with spring constant of 10 kcal/mol$\cdot\mathring{\mathrm{A}}^{2}$ and the tolerance of 2.39 kcal/mol$\cdot\mathring{\mathrm{A}}^{2}$ without maximal step limits (Wang et al., ).

### D.4 Algorithms for Training and Inference

We provide pseudo codes for training and inference with our models FBM-BASE and FBM in Algorithm 1,2,3 respectively.

Figure: Algorithm 1 Training with FBM-BASE

Figure: Algorithm 2 Training with FBM

Figure: Algorithm 3 Autoregressive inference with FBM/FBM-BASE

### D.5 Hyperparameters

The hyperparameters we choose are listed in Table 3.

**Table 3: Hyperparameter choice of FBM-BASE and FBM.**
| Hyperparameters | Values |
| --- | --- |
| Network |  |
| Hidden dimension $H$ of FBM-BASE | 128 |
| Hidden dimension $H$ of FBM | 176 |
| RBF dimension | 32 |
| Number of attention heads | 8 |
| Number of layers | 6 |
| Cutoff threshold $r_{\mathrm{cut}}$ | 5.0$\mathring{\mathrm{A}}$ |
| Training |  |
| Learning rate | 5e-4 |
| Optimizer | Adam |
| Warm up steps | 1,000 |
| Warm up scheduler | LamdaLR |
| Training scheduler | ReduceLRonPlateau(factor=0.8, patience=5, min_lr=1e-7) |
| Batch size of FBM-BASE | 16 |
| Batch size of FBM | 10 |
| SDE noise scale $\sigma$ | 0.2 |
| Inference |  |
| SDE steps $T$ | [25,30] |
| Guidance strength $\eta$ of FBM | [0.04,0.05,0.06,0.07,0.08] |

## Appendix E Experimental Details

### E.1 Dataset Details

As mentioned in § 4, all peptides of PepMD are simulated in explicit solvent model TIP3PFB using OpenMM (Eastman et al., ). The parameters we used for MD simulation are listed in Table 4 and the statistical information of PepMD is shown in Table 5.

Additionally, all 14 peptides of our test set are listed below with the format {pdb-id}:{chain-id}: 1hhg:C, 1k8d:P, 1k83:M, 1bz9:C, 1i7u:C, 1gxc:B, 1ar8:0, 2xa7:P, 1e28:C, 1gy3:F, 1n73:I, 1fpr:B, 1aze:B, 1qj6:I.

**Table 4: MD simulation setups using OpenMM.**
| Property | Value |
| --- | --- |
| Forcefield | AMBER-14 |
| Integrator | LangevinMiddleIntegrator |
| Integration time step | 1fs |
| Frame spacing | 1ps |
| Friction coefficient | 1.0$\text{ps}^{-1}$ |
| Temperature | 300K |
| Solvation | TIP3PFB |
| Electrostatics | NoCutoff |
| Constraints | all bonds between hydrogens |

**Table 5: Dataset statistics.**
| Dataset name | PepMD |
| --- | --- |
| Training set simulation time | 100ns |
| Test set simulation time | 100ns |
| MD integration time step $\Delta{t}$ | 1fs |
| coarsened predition time $\tau$ | $0.5\times{10^{6}}$fs |
| # Clusters | 2480 |
| # Training peptides | 136 |
| # Training pairs per peptide | $2\times{10^{3}}$ |
| # Validation pairs per peptide | $4\times{10^{2}}$ |
| # Test peptides | 14 |

### E.2 Details on Evaluation Metrics

In this part, we provide details for computing the evaluation metrics in § 4.

#### Flexibility

Following Janson et al. (), we compute the contact rates between residues as a measure of structural flexibility. For each residue pair $i,j$ ($1\leq{i}<j\leq{R}$) of a peptide with $R$ residues, the contact rate $r(i,j)$ of residue $i,j$ is defined as follows:

$$ $r(i,j)=\frac{1}{L}\sum_{l=1}^{L}\bm{1}_{d_{l}(i,j)<10\mathring{\mathrm{A}}},$ (65) $$

where $d_{l}(i,j)$ denotes the Euclidean distance between $\alpha$-carbons of residue $i,j$ of conformation $l$. Now we compute the root mean square error of contact maps between generated conformation ensembles and reference MD trajectories:

$$ $\text{C{ONTACT}}=\sqrt{\frac{2}{R(R-1)}\sum_{1\leq{i}<j\leq{R}}(r(i,j)-r_{ \mathrm{ref}}(i,j))^{2}}.$ (66) $$

#### Validity

We assess the structural validity by checking for bond breaks between adjacent residues and bond clashes between any residue pairs. The same as in Wang et al. (), bond clash occurs when the distance between $\alpha$-carbons of any residue pair is less than the threshold $\delta_{\mathrm{clash}}=3.0\mathring{\mathrm{A}}$, and bond break occurs when the distance between adjacent $\alpha$-carbons is greater than the threshold $\delta_{\mathrm{break}}=4.19\mathring{\mathrm{A}}$. Then the metric VAL-CA is assessed by the fraction of conformations without bond break and bond clash.

#### Distributional Similarity

Similar to Lu et al. (), we project peptide conformations onto the following three low-dimensional feature space: (i) Pairwise Distance (PWD) between $\alpha$-carbons excluding residue pairs within an offset of 3. (ii) Radius of gyration (RG) which computes the geometric mean of the distances from $\alpha$-carbons to the center-of-mass. (iii) Time-lagged Independent Components (TIC), where we featurize structures using backbone dihedrals $\psi,\phi,\omega$ and pairwise distances between $\alpha$-carbons (Klein et al., ), then TIC analysis is performed using Deeptime (Hoffmann et al., ). Only the slowest components, TIC 0 and TIC 1, are taken for further evaluation (Pérez-Hernández et al., ). (iv) the joint distribution of TIC 0 and TIC 1, termed as TIC-2D. (v) Specifically for the evaluation on AD, the joint distribution of backbone dihedrals $\psi$ and $\phi$, namely the Ramachandran plot (Ramachandran et al., ), is taken into consideration (RAM).

Afterwards we compute the Jensen-Shannon (JS) distance between generated samples and reference MD trajectories on the projection space. Features are discretized with 50 bins based on the reference ensembles, and a pseudo count 1e-6 is added for numerical stability. For each feature space, we report the mean distance along all dimensions.

## Appendix F Additional Experimental Results

### F.1 Ablation Study

In Table 6 we provide ablation results of SDE steps $T$ and the guidance strength $\eta$ on test peptides of PepMD. A clear pattern is that when $T$ is fixed, the greater the guidance strength $\eta$, the more likely it is to generate reasonable conformations, which demonstrates a strong correlation between the intermediate force field and real interatomic constraints of molecular systems.

**Table 6: Ablation results of SDE steps $T$ and the guidance strength $\eta$ on the test set of PepMD. Values of each metric are first averaged over 3 independent runs for each peptide and then shown in mean/std of all 14 test peptides.**
| Hyperparameters | JS DISTANCE ($\downarrow$) | VAL-CA ($\uparrow$) | CONTACT ($\downarrow$) |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| PWD | RG | TIC | TIC-2D |  |  |  |
| $T=25,\eta=0.04$ | 0.586/0.059 | 0.540/0.143 | 0.640/0.056 | 0.809/0.020 | 0.559/0.190 | 0.195/0.115 |
| $T=25,\eta=0.05$ | 0.578/0.069 | 0.550/0.151 | 0.639/0.075 | 0.805/0.027 | 0.606/0.192 | 0.179/0.132 |
| $T=25,\eta=0.06$ | 0.573/0.064 | 0.542/0.140 | 0.631/0.077 | 0.801/0.032 | 0.616/0.188 | 0.188/0.127 |
| $T=25,\eta=0.07$ | 0.585/0.084 | 0.548/0.177 | 0.644/0.089 | 0.800/0.044 | 0.647/0.182 | 0.201/0.135 |
| $T=25,\eta=0.08$ | 0.573/0.064 | 0.577/0.144 | 0.638/0.079 | 0.804/0.033 | 0.675/0.169 | 0.200/0.111 |
| $T=30,\eta=0.04$ | 0.585/0.074 | 0.587/0.151 | 0.644/0.077 | 0.804/0.029 | 0.553/0.208 | 0.208/0.131 |
| $T=30,\eta=0.05$ | 0.582/0.066 | 0.542/0.162 | 0.640/0.088 | 0.806/0.035 | 0.615/0.183 | 0.217/0.130 |
| $T=30,\eta=0.06$ | 0.596/0.085 | 0.591/0.139 | 0.637/0.089 | 0.803/0.041 | 0.604/0.205 | 0.231/0.116 |
| $T=30,\eta=0.07$ | 0.576/0.059 | 0.590/0.120 | 0.649/0.074 | 0.806/0.023 | 0.655/0.180 | 0.182/0.116 |
| $T=30,\eta=0.08$ | 0.590/0.077 | 0.575/0.151 | 0.614/0.101 | 0.789/0.065 | 0.661/0.176 | 0.208/0.116 |

### F.2 PepMD Additional Results

More examples of the PepMD test set are shown in Figure 6. Compared to the baselines, FBM presents more accurate depictions of the molecular free energy landscape in most cases. Due to potential issues such as insufficient generated samples or the limited scale of the training set, there is still room for improvement in the prediction accuracy and transferability of FBM.

Figure: Figure 6: Experiments on PepMD test peptides 1n73:I, 1gxc:B and 1qj6:I (top, middle and bottom). Samples were generated in the time-coarsened manner for a chain length of $10^{3}$. Free energies (*i.e.*, the relative log probability) along the first two TIC components of FBM, FBM-BASE and Timewarp are displayed in the left, middle, and right columns, respectively. The blue solid line represents MD-generated samples, while the yellow dashed line represents model-generated samples.
Refer to caption: x14.png

## Appendix G Computing Infrastructure

Our models, FBM-BASE and FBM, were trained on 4 NVIDIA GeForce RTX 3090 GPUs within a week. The inference procedure with baselines and our model were all performed on one single NVIDIA GeForce RTX 3090 GPU.

## Peer Reviews

confidence: 5
contribution: 1
initial_rating: 3
presentation: 3
questions: 1. Using bridge matching may not be ideal, as all samples along the same trajectory are expected to follow the same Boltzmann distribution. Ideally, we would achieve **a constant transformation**; however, in practice, we are performing point-to-point correspondence (from one delta distribution to another). This approach, therefore, appears theoretically unsound. What are your thoughts on this?
2. Experiments on peptides alone are too limited and should at least be performed on the smallest “protein” (backbone), Chignolin.
soundness: 2
strengths: 1. The paper is clearly written and easy to follow, also this task is indeed important.
2. The experiments are effectively support the conclusions.
3. The code and detailed hyperparameters are provided, making it easy for others to reproduce the results.
summary: The paper proposes a conditional generative model called Force-guided Bridge Matching (FBM) to improve Molecular Dynamics (MD) simulations. Traditional MD simulations, while accurate, are time-consuming due to the need for small time steps. FBM leverages a physics-informed approach to time-coarsened dynamics by incorporating an intermediate force field to guide the generation process. This allows the model to better approximate the Boltzmann distribution, essential for realistic molecular conformations, and improves simulation efficiency without sacrificing accuracy.
weaknesses: I believe the main weaknesses with this paper is its lack of novelty, as it appears to be a straightforward combination of [1], [2] and [3] in a slightly altered setting. It seems the author did not conduct a comprehensive literature review. My reasons are as follows:

1. Training with trajectory pairs was already introduced in [1]; the only difference here is that [1] focused on protein backbones and employed SE(3) FM. In addition, next frame prediction is also a subtask in [2].
2. The concept of force-guided training was discussed in [3], with the main difference being that [3] uses SE(3) diffusion.

In essence, this paper combines aspects of both [1] ([2]) and [3] without properly referencing [1], [2]. While it presents a solid engineering application, it lacks originality.

Additionally, I have another concerns, which I’ve outlined in the Questions section.

>[1] Li, S., Wang, Y., Li, M., Shao, B., Zheng, N., Jian, Z., & Tang, J. F $^ 3$ low: Frame-to-Frame Coarse-grained Molecular Dynamics with SE (3) Guided Flow Matching. In ICLR 2024 Workshop on Generative and Experimental Perspectives for Biomolecular Design.

>[2] Jing, B., Stark, H., Jaakkola, T., & Berger, B. Generative Modeling of Molecular Dynamics Trajectories. In ICML'24 Workshop ML for Life and Material Science: From Theory to Industry Applications.

>[3] Wang, L., Shen, Y., Wang, Y., Yuan, H., Wu, Y., & Gu, Q. Protein Conformation Generation via Force-Guided SE (3) Diffusion Models. In Forty-first International Conference on Machine Learning.

---

confidence: 4
contribution: 4
initial_rating: 8
presentation: 4
questions: The TIC plots in Figure 4 are not entirely clear in demonstrating the advantages of FBM. It might be more interpretable for readers to adjust the color scheme, explore alternative visualizations, or include additional quantifiable metrics.
soundness: 3
strengths: - The integration of force-guided priors within the bridge matching framework for MD simulations is innovative. This approach goes beyond standard generative modeling by incorporating Boltzmann-distributed physical priors, improving the model's ability to simulate realistic molecular transitions over extended time scales. 
- The authors conduct a comprehensive evaluation of FBM on two peptide datasets, demonstrating that it consistently outperforms baseline models across several metrics, including distribution similarity and validity of generated molecular conformations. 
- The paper is generally well-structured, with a clear presentation of the methodology and objectives. 
- The model’s ability to generate high-quality conformations with transferability across unseen systems is particularly valuable for applications requiring generalized and scalable MD simulations.
summary: The paper presents a novel model, FBM, which integrates force guidance into a bridge matching framework for time-coarsened dynamics. Traditional MD methods often struggle with achieving a balance between computational efficiency and accuracy, especially over long time scales. FBM addresses this challenge by directly targeting a Boltzmann-constrained distribution, effectively combining physical priors with generative modeling to capture full-atom time-coarsened dynamics.
weaknesses: - The experiments to evaluate the model’s performance are conducted mainly on datasets consisting of small peptides. It would provide a more comprehensive understanding of the model's performance to extend the evaluation to include larger peptides or more complex molecular structures.
- While the FBM model demonstrates promising results, a more extensive comparison with traditional MD simulations across additional cases would strengthen its credibility. It would be helpful to provide more detailed time comparisons between FBM and MD simulations.

---

confidence: 3
contribution: 3
initial_rating: 8
presentation: 3
questions: Overall all those figures and metrics are complementary and I believe for a few samples to have all of those would help a lot. Indeed we have the different JS values for the full test set but having them also for the samples presented, as well as having directly the distribution comparison to the ground truth as in Figure 6 would be super helpful. I  believe it can be hard to understand what a JS value means except that it is bigger for the output of one model compared to the other. 
It could be also nice to have more of those figures for your best-picked trajectories but also for those that didn’t work as well: again just to have a feel about what could be driving those numbers.
Moreover the study of hyperparameters effect (table 6) is a great add to show the robustness of the results. I believe it could be used to show the effect of force guidance on the trajectories generated and their metrics:  mainly following for a single or a few peptides how their trajectory change when the strength of guidance increases. This is touched upon in section f1 where it is pointed out that validity and contact metrics indeed nicely follow the guidance strength, but I feel like this is less true/obvious for the other metrics.
Finally, it would also be nice to see how easily this is applicable in more real case scenarios, with some scaling laws with peptide size, etc…

Overall I would say that all of my remarks are minors, but I still encourage  the authors to follow them as much as possible as I believe it will ease the overall understanding of the model ability. I would also like to state that even though I was able to follow all the proofs and didn’t see anything that I would deem to be problematic, I don’t believe I  have the level of expertise truly required to judge the math present in this paper with certainty.
soundness: 3
strengths: The paper offers a rather intuitive and straightforward way to add some physics-derived constraints to the trajectory produced at prediction time. Moreover, their model is state of the art when taking into account a rather diverse set of evaluation metrics: from comparing model output collective variables at the distribution level, in terms of validity through checking bond clashes and breaking, or finally in term structure flexibility compared to ground truth. The simplicity of adding or not the correction term in a tunable way through a guidance strength hyperparameter allows for a direct ablation study which shows the usefulness of the correction (guidance) term.
The care brought into splitting the dataset (clustering and splitting by  cluster), indeed allows them to claim generalizability, even though I  didn’t really understood why AD seems to be treated a bit separately  from the peptide in the test set.
Last but not least, it also proposes a robustness/explanatory study by adding the results of their hyperparameter tunning and some effects of those hyperparameters.
summary: In this paper the authors propose a model to produce relevant MD  trajectories from an initial structure, using bridge diffusion with a possibility to ensure that the recovered conformation statistic obeys a  Boltzmann distribution rather than blindly following the learned data distribution. They accomplish their goal by constraining the form of the conformation distribution to be the learned data distribution with each conformation being weighted by a learnable Boltzmann factor. So the constraints appeared from the beginning and are present in the trained model from the get-go: no need at prediction time to go through expensive re-weighting procedure, as it is part of the model. Hence their model relies on first learning the reverse and forward process that would generate the initial and final experimental conformations (conditional to those), then use this newly learned score function and use it within a dedicated network to calculate the associated Boltzmann weight (i.e. the conformation energy). At prediction time one can either just perform the simple diffusion, or decide to add the correction term to the simple diffusion-associated force field, which turns out to “simply” be the gradient of the energy term in the learned Boltzmann re-weighting coefficient: hence the name.
weaknesses: The paper is pretty strong and its only rather minor weakness has to do with easiness to have a direct interpretation/overall view of the metrics for the samples presented. Right now the metrics show their claims (their model is better) which is great, but they, I  believe, are hard to link to the model's performance and usefulness. 
For example, I think Figure 6 is nice because it shows a direct comparison with ground truth in a more obvious way than other figures. Hence the few following questions and propositions of improvement described in the next section will be about that.

Meta-review: This paper received mixed scores, with two reviewers leaning to a clear accept (score 8), and one reviewer initially learning towards recommending to reject this paper (score 3). Among other things, the paper was thought to have a clear presentation, and claims that were validated with their empirical results, which included evaluations on alanine dipeptide and small dipeptides. As areas for improvement, points were raised on the need for validation on more than just small peptides and the need for timing comparisons. Furthermore, the reviewer who recommended rejection indicated that they viewed the work as a straightforward combination of existing concepts in prior work, and that the work did not show an advantage through this combination. 

During the rebuttal, the authors have added additional experiments for the small protein Chignolin, where they compare their model against MD only. Additionally, the authors included comparisons of effective sample size per second of wall-clock time as an additional metric for efficiency, which was introduced in prior work.

Based on this rebuttal and following discussions the two reviewers in favor of acceptance maintained their score, while the reviewer with the most concerns raised their score from a 3 to a 5, based on the remaining concern of limited novelty.

After the author-reviewer discussion ended, I read through the reviews again, as well as the discussions with the authors and the newly computed results and the revised paper. Upon going through the revised paper to look through the added results, I found that the results on alanine dipeptide raise questions about their correctness. In particular, the distribution of samples from the models in figure 3 do not seem to resemble the distribution represented by the contours based on kernel density estimation of MD trajectories, and the baseline methods in particular do not perform well at all. This discrepancy between model distributions and the MD-based distribution is striking since one of the baseline methods, Timewarp, has shown in their original paper (https://arxiv.org/pdf/2302.01170) to be able to draw samples of alanine dipeptide according to a distribution that closely mimics the distribution of samples from MD trajectories. See figure 3 in the Timewarp paper. This raises doubts for me on whether or not the results in this current submission were calculated correctly. 
After consulting with the reviewers, the reviewers agreed that the results raise concerns about correctness. Unfortunately, the range of angles used to compute Ramachandran plots of alanine dipeptide in this submission were not displayed, nor are samples from the reference MD simulation included in the Ramachandran plots. This makes it hard to determine if the discrepancy between the results reported in this paper compared to the Timewarp paper are influenced by different data or plotting settings. Moreover, we received input from one additional reviewer with significant expertise in this field, who was initially unable to review the paper but looked through the results and other reviews. In addition to sharing concerns about the alanine dipeptide results, they noted that the TICA plots in figure 4 show worryingly bad results, and that the energies in Figure 6 are also not good. Unfortunately, during the reviewer-AC discussion phase it was no longer possible to ask for clarification from the authors, and the presented results did not have enough details to be able to determine where the discrepancies with previously reported results come from. Given the overall consensus among reviewers and AC that the presented results raise concerns, I recommend rejecting this paper. 

I would like to thank the authors for submitting the paper and for providing additional results during the rebuttal. I hope the above explanation can be helpful in rewriting the manuscript for a future submission.

Decision: Reject