BioMD: All-atom Generative Model for Biomolecular Dynamics Simulation
Authors: Bin Feng
Jiying Zhang
Xinni Zhang
Zijing Liu
Yu Li
Venue: ICLR 2026 Poster (2026)
Keywords: molecular dynamics
biomolecular trajectories generation
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)

Abstract: Molecular dynamics (MD) simulations are essential tools in computational chemistry and drug discovery, offering crucial insights into dynamic molecular behavior. However, their utility is significantly limited by substantial computational costs, which severely restrict accessible timescales for many biologically relevant processes. Despite the encouraging performance of existing machine learning (ML) methods, they struggle to generate extended biomolecular system trajectories, primarily due to the lack of MD datasets and the large computational demands of modeling long historical trajectories. Here, we introduce BioMD, the first all-atom generative model to simulate long-timescale protein-ligand dynamics using a hierarchical framework of forecasting and interpolation. We demonstrate the effectiveness and versatility of BioMD on the DD-13M (ligand unbinding) and MISATO datasets. For both datasets, BioMD generates highly realistic conformations, showing high physical plausibility and low reconstruction errors. Besides, BioMD successfully generates ligand unbinding paths for 97.1% of the protein-ligand systems within ten attempts, demonstrating its ability to explore critical unbinding pathways. Collectively, these results establish BioMD as a tool for simulating complex biomolecular processes, offering broad applicability for computational chemistry and drug discovery.

Published as a conference paper at ICLR 2026
BIOMD: ALL-ATOM GENERATIVE MODEL FOR BIOMOLECULAR DYNAMICS SIMULATION
Bin Feng's : Jiying Zhang! Xinni Zhang's Zijing Liu', / Yu Li
Interraaional Digital Economy Academy, "Siss Federal Technology Institute of Lausanne (EPFL)
ABSTRACT
Molecular dynamies (MD) simulations are essential tools in computational chem- istry and drug discovery, offering crucial insights into dynamic molecular behavior. However, their utility is significantly limited by substantial computational costs, which severely restrict accessible timescales for many biologically relevant pro- cesses, Despite the encouraging performance of existing machine learning (ML. methods, they struggle to generate extended biomolecular system trajectories, primarily due to the lack of MD datasets and the large computational demands of modeling long historical trajectories. Here, we introduce BioMD, the first all- atom generative model to simulate long-ting long-timescale ligand dynamics using a hierarchical framework of forecasting and interpolation. We demonstrate the effec- tiveness and versatility of BioMD on the DD-13M (ligand unbinding) and MISATO datasets. For both datasets, BioMD generates highly realistic conformations, show- ing high physical plausibility and low reconstruction errors. Besides, BioMD successfully generates ligand unbinding paths for 97.1% of the protein- systems within ten attempts, demonstrating its ability to explore critical unbinding pathways. Collectively, these results establish BioMD as a tool for simulating complex biomolecular processes, offering broad applicability for computational chemistry and drug discovery.
INTRODUCTION
Molecular dynamics (MD) simulations have emerged indispensable tool in computational chemistry and drug discovery, offering insights into the dynamic behavior of biomolecular systems. Through numerical integration of Newton's equations of motion, MD simulations directly produce atomic trajectories that reveal the time evolution of molecular structures (Hollingsworth & Dror 2018). These trajectories enable the exploration of conformational ensembles, optimization of small molecule structures, and identification of potential binding sites, significantly accelerating the design and development of novel therapeutics (Karplus & McCammon (002).
Despite their utility, traditional MD simulations face substantial computational limitations. Al- though algorithms like Particle Mesh Ewald (PME) have improved scaling to approximately O(N N) (Darden et al. 1993). these calculations remain the most computationally demand- Adeock
ing component (Dror al. Adccck & McCammon] 2006). Furthermore, accurately resolving high-frequency atomic vibrations necessitates extremely small time steps (on the order of femtosec- onds), severely limiting the accessible simulation timescales (Shaw et al. 2010 2009). Exploring biologically relevant processes, which often span microseconds to milliseconds, remains computa- tionally intensive, restricting the practical application of atomistic MD to obtain trajectories.
Recently, machine learning (ML) methods have emerged computational alternatives molecular dynamics (MD) simulations. Key advances include models for generating protein conformation en- sembles (Lewis et al. 2025) and neural network potentials trained on quantum mechanical data (Wang) et al. 2024a). For biomolecular systems, AlphaFold 3 (Abramson et al. (0224) has demonstrated promising accuracy in predicting protein-ligand interactions. Despite these achievements, generating full MD trajectories for complex protein ligand systems using ML remains major challenge. Exist- ing approaches tend fall into two categories: (i) methods that can generate protein conformation
"These authors contributed equally to this work.
1 Corresponning authors: liuzijing@ @idea.edu.cm liyu@idea.edu..c

---

Published as a conference paper at ICLR 2026
Average Ligand RMSD between frames
b
stop
step 10
stop 20
stop 40
step 80
stop 160
DD:13M
12.5
Averoge RMSD (4)
10.0
7.5
2.5
0.0
100
200
300
400
500
Frome Interval (steps)
Figure 1. Average Ligand RMSD between two frames. (a) I. Line plot showing that the average ligand RMSD between two frames the same trajectory increases with the frame interval. (b) Examples of ligand unbinding trjectoriis at time steps 0, 10, 20, 40, 80, and 160.
ensembles but cannot produce time-resole trajectories Jing et al. (2024b): Wang et al. (2024b), or (ii) methods that attempt trajectory modeling but struggle to capture protein protein-ligand ligand interactions. For example, NeuralMD (Liu et al., 2024) treats protein atoms as static and only models ligand dynamics, while MDGen (Jing et al.) 2024c) is specifically designed for peptides and proteins and does not handle small-molecule ligands. This limitation arises from both the complexity of -iggann-ligand energy landscapes and the scarcity of high-quality trajectory data for training generative models.
To address these limitations, we propose BioMD, hierarchical framework for generating all-atom biomolecular trajectories. Building upon the insight that short-timescale conformational changes exhibit little conformational change (Figure I.. BioMD decomposes long trajectory generation into two synergistic stages: forecasting of large-step conformations, followed by interpolation to refine intermediate steps. This strategy reduces sequence length by decoupling long-term evolution from local dynamics and helps manage the error accumulation problem for generating long trajectories. Crucially, BioMD unifies forecasting and interpolation within a conditional flow matching model, where we use the "noising-as-makking" methods following Diffusion Forcing (Chen et al. 2024) to our time-scale transformer. We apply independent noise to each frame, which enables flexible conditioning on partial trajectory segments, and we implement different tasks simply by using different masking schedules. Inspired by the success of AlphaFold 3, BioMD generates all-atom trajectories using velocity network that adapts its core transformer architecture, while employing an SE(3) -qqiivrruaa graph transformer to encode the initial conformation as conditional embeddings.
To evaluate the effectiveness of BioMD, we conducted experiments on two datasets: MISATO (Sieben) morgen et al. 2024) and DD-13M (Li et al. 2025). Our results show that BioMD generates highly realistic conformations with promising physical stability, evidenced by low energy and reconstruction errors across both benchmarks. On the MISATO dataset, which focuses on ligand dynamics within the binding pocket, our model accurately captures the system's 's conformational flexibility, outperforming existing methods. For the more challenging task of ligand unbinding on the DD-13M dataset, BioMD) successfully generates complete unbinding paths for up to 97.1 1% of the protein-ligand systems, demonstrating a robust ability to explore critical and tin-ttmmeccle biomolecular pathways. Col- lectively, these results establish BioMD as a powerful and efficient tool for simulating complex biomolecular processes, offering broad applicability for computational chemistry and drug discovery.
RELATED WORKS
Conformational Ensemble and Binding Pose Generation. One major line of research uses ML to generate biomolecule's conformational ensemble by modeling the equilibrium distribution of its dynamic structures. Early efforts like AlphaFold2 (Jumper al., 2021) produce set of diverse conformations primarily through multiple sequence alignment (MSA) subsampling and masking techniques (Stein & Mchaourab] 2022] del Alamo et al. 2022; Wayment- Steelle et al. 2024). More
advanced approaches now directly learn the conformational distribution from large-scale MD datasets using flow-based (Noé et al. 2019 Jing et al. 2024b) or diffusion-basee
(Wang et al. 2024b Jing et al. 2023) [u et al. 2024a Zheng et al. 2024 Lu et al. 2025a) generative models. Models such as BioEmu (Lewis et al. 2025) can effectively generate diverse and physically plausible conformations, providing powerful alternative to extensive MD sampling understand conformational space.
2

---

Published as a conference paper at ICLR 2026
Beyond proteins, recent works have extended these generative frameworks to protein-ligand systems. DynamicBind (Lu et al. 2024b) and DynamicFlow (Zhou et al. (2225) leverage diffusion and flow matching models, respectively. recover key protein-ligand binding poses while accounting for receptor flexibility. However, these methods are fundamentally time-agnostic; they can sample what conformations possible (or the final binding states) but lack the temporal information to show the kinetic pathways between them.
Trajectory Learning for MD Simulation. To capture these kinetic pathways, complementary research direction aims to generate full, time-ordered trajectories. Approaches like EquiJump (Costal 2024 learn to sample future states based solely on the current conformation. To capture higher-order dependencies between the frames, MDGen (Jing et al. 2024c) models the joint dis- tribution of entire trajectories via masked frame modeling. ConíRover (Shen al. 2025 models these dependencies nuto-regrrssively by conditioning each frame on its entire history through causal transformer. While powerful, these methods are often specialized for protein-only dynamics. Conversely, methods that model protein-ligand interactions often introduce other simplifications. For instance, NeuralMD (Liu et al. 2024) treats the protein receptor as static. These simplifications limit their scope of accessible dynamics.
3 PRELIMINARIES
Notations. A complex C is composed a protein P and a ligand L. The trajectory of a complex contains T frames of coordinates, denoted as X, = (xo, X1, Xr) E R((++)xxxx3 where X₁ [x], J x3 represents the concatenation of protein coordinates x/" and ligand coordinates x at time-step t, and N the number of atoms in the complex. The complex trajectory prediction task is defined as generating subsequent conformations (coordinates) of complex trajectory given its initial conformation (i.e., the first frame).
Molecular dynamics. Molecular dynamics (MD) simulates the time evolution of particle system under classical mechanics. It leverages numerical schemes such as Verlet integration (Verlet) 1967) or Langevin dynamics to generate trajectories approximating the Boltzmann distribution. In the simplest deterministic case with no friction or noise, each particle evolves according to dæ; = Ps dt, dp: =
~Va, (z) dt, where ps and m, are the momentum and mass, and E(x) is the potential energy function. Metadynamics (Laio & Parrinello, 2002; Barducci et al. 2011 Li al.) 2025) extends MD by introducing a distory-dependent bias potential V(s, L), constructed over collective variables
x(x) as (s,,() = Evec 1 exp [((a(t)) a((()))) where Gaussians of height w and width
2n
are periodically added discourage revisiting explored states. This bias fills free-energy wells and enhances sampling of rare events and transition pathways beyond the reach of standard MD.
Flow matching based models. Flow matching (FM) (Lipman et al. 2023 an efficient and simulation-fii free method for training continuous normalizing flows (CNFs), class of generative mod- els based on ordinary differential equations (ODEs). Euclidean space, CNFs define transformation (.) IRN Nx3 + R.V.V via an ODE governed by tine-dependem velocity field e,:
dr の (x2) = e, (ør (x")), ex = x², TE [0, 1],
where x sampled from a simple distribution Pa. and o evolves it over time 7 E [0, 1] to match 1. Since u, is unknown, FM learns v, by regressing the conditional flow u(Ø, (x" (x") o (x" x²)), where or (x" x²) interpolates between x" Po and x ~ P1- In our setting, each conformation E N
the target distribution Pj at
represents frame in complex trajectory, and FM is used to generate future frames from an initial structure.
BIOMD METHOD
4.1
A UNIFIED GENERATIVE FRAMEWORK VIA FLOW MATCHING
Our model capitalizes on a fundamental insight into molecular dynamics: conformational changes typically subtle over short timescales but can involve significant global movements over longer
3

---

Published as a conference paper at ICLR 2026
b
Forecasting
schedule
Forecasting
Tg=1 T=U((011)
(x~U(0,1) (--U(011))
First stage
Interpolation
G-1 T..1~U(011) (v,,-(((111) Tap=
X
Xax
Interpolation
BioMD
Second
stage
Xiak
Ap-1
Flow vector field
Figure 2. Model framework. (a) The hierarchical framework, showing the two-stage process of coarse-gres coarse-grained forecasting followed by fine-grained interpolation. (b) The time scheduling mechanism for forecasting and interpolation tasks, where known frames are noise- free (r
generated frames are noised (r E [0. 1)..
This principle underpins our hierarchical prediction framework, which decom- poses the generation of long trajectories into two principal stages: coarse- consse-griiee forecasting and fine-grained interpolation (Figure 2).
timescales (Figure D
Notably, this entire framework is implemented within single model architecture that processes the sequence of the whole trajectory at once. We adopt a "noise as mask" strategy, where the distinction between the two stages is made simply by varying the input masking patterns (Figure 2b). In this unified framework, each frame in an input sequence is independently perturbed by noise according time variable т. Known or conditioning frames are kept clean (equivalent setting their corresponding т = 1, i..., "unmasked"), while frames to be generated are initialized from pure 0, i... "masked") and then iteratively denoised.
noise (equivalent
Nt, J. During training. we sample vector of independent time steps T (Tt, Ttg....., Tt. J. where each tr, (O, 1). The sequence is then noised to X = ixt,,... x7, ). where each frame is an interpolation between the real
Let trajectory sequence be denoted by X = ,
coordinates and Gaussian noise E; (O I): x7, = Xr, + (1-Tu, Ts, ),, The corresponding ground- truth velocity field for the sequence is U = {u ...., u. h. with u. (x,, - x1,)/(1- -1-11- -T.)).
Our velocity model uo takes the entire noisy sequence and conditioning information predict the velocities for all frames simultaneously. The training objective a Mean Squared Error loss over the entire sequence:
Lnow MSE(up (XT. T,, UT).
where Z contains static information including the first frame coordinate Xo, amino acid sequence and ligand atom types a. We explore two modeling approaches: BioMD-rel, which predicts coordinate changes relative to an anchor frame, and BioMD-abs, which predicts absolute atomic coordinates. For clarity, we focus on the absolute coordinate prediction task below.
4.2 HIERARCHICAL GENERATION WITH FORECASTING AND INTERPOLATION
The two stages of our hierarchical framework are realized simply by applying different masking schedules to our unified model during training and inference.
4.2. COARSE-GRAINED FORECASTING
The first stage generates a grarse-grained trajectory, constructed by sampling every k: steps (empirically chosen, ablation study in Appendix A.5) from the full trajectory, resulting in sequence Xc {xo, X3, .2...... }. This task is framed as a forecasting problem where, given the initial frame Xp, the model must generate all subsequent frames.
This is achieved by applying a specific masking schedule to our unified framework. During training, the time step for the initial frame is fixed the time steps for all other frames {Tk. T2k: ... are sampled independently from (O, 1). The model ug is trained predict the velocities for all frames in the sequence, conditioned the initial frame.
1 (making it known, "unmasked" condition), while
4

---

Published as a conference paper at ICLR 2026
First stage
Second stage
Flow
Frame
Single repe.
SE(3)-graph transf ooreer
(Rayem:: #head:15)
Atcm input repe.
rept.
Faracesting
Interpplation
stop
Atom pain
Global attention (tokens) (ffeedi11))
2.
local (atoms) (Phoar-)
local (atons) jenuss)))
X3
X3
Novey ccordinctes
Flow vector field
Kenzk
Forecasting: first frame is known, other frames begin with noise: Interpolation: first frome and last frame are known, other fromes begin with noise
Figure 3. Detailed architecture of BioMD. The model operates in modes, Forecasting and Interpolation, set up by the hierarchical framework (left). The core velocity network (right) processes noisy coordinates, conditioned on features from SE(3)-Graph Transformer. A local-global-local attention pathway generates the final flow vector field used for trajectory generation.
During inference, this setup supports multiple generation strategies:
All future frames (xk, X2k:... are generated concurrently. We set To 1, initialize all other frames from noise (i.e., their 7 values start at 0), and use an ODE solver like the Euler method to integrate all frames simultaneously to 1.
All-at-once:
Auto-regressive (AR): Frames are generated in sequential blocks of size j. To generate one such block. the model conditions on the previously generated history. This is controlled by the time variable r: the values for all frames in the history are set to be constant 1, making them clean, "unmasked" inputs. The values for j frames within the current target block are then jointly evolved from 0 to by the ODE solver. This process simultancously denoises all frames in the block, using the generated history as context. Once generated, this block is added to the history, and the process is repeated for the next block until the full trajectory is complete.
4.2.2 FINE-GRAINED INTERPOLATION
After obtaining the coarse-grained trajectory {xo, Xj:, X. J. the second stage replenishes the intermediate frames. This is interpolation task, where for each coarse interval, we generate the frames (x.tt1...., X(+1)) conditioned on the two "anchor" frames, X, and X(++))a-
This task uses the exact same velocity model ug and training framework, differing only in the data and masking schedule. The input sequence is now a fine-grained segment X (xa,, , )(+1)))) During training, the anchor frames are designated as known by fixing their time steps = and T(++))& = 1. The time steps for all intermediate frames are sampled independently from (0, 1). The model learns to generate the intermediate trajectory conditioned on the start and end conformations.
During inference, this task is always performed in an "all-at at-once" manner. The anchor frames X and ((++))& are provided as clean inputs (their 1), while all intermediate frames are initialized from noise (their = 0). The model then simultaneously generates all l: intermediate frames by integrating them 1. This process is described by:
+ Yi + 1o (T. Zueq T) AT,
where Yu represents the block of intermediate frames, and the velocity predictions extracted for only those frames. This hierarchical approach allows BioMD to efficiently generate long, physically plausible trajectories.
5

---

Published as a conference paper at ICLR 2026
Table Results on the MISATO test set. Comparison of methods on physical stability (first six metrics)
and conformational flexibility (last four metrics). Mean values the test samples are reported.
Geometry'
Angle Geometry"
Clashes
RMSF Comelation"
Value"
Method
MAE
MSE
MAE
MSE
Intra-Lig
Prui-Lgg
Ligund
Protein
Ligand (1.211)
Proiein (1.002)
Molecular Dynamics
.0377
.0023
.0555
.0033
Denoisingt.D
10
1027
10188
.0431
.0160
.0295
-0.0290
1002
GNNMD
.1133
1032
22115
1072
.3626
0028
-0.0103
.2165
NeuraMD.
.0883
.0776
0605
.086
.0114
.0578
.3405
.3220
NewnaMD--SDF
.083
.0076
.0604
0086
VerlerMD
19.73
1050
5847
5482
.183
3.111
3356
.3266
BioMD-rel
.0395
.0026
.0555
.0075
.0003
.0006
.8611
.5945
.369
.5177
BioMDDabs
.0495
.0.55
.0709
.0997
.0019
.0033
.789
.6854
.7023
.6242
Bond (bond length) and RMSF values are in ungsiroms (A). Angle geometry (bond angle) is in radians.
RMSF Correlation is reported using Pearson
coefficient.
RMSF values for reference trajectories are given in parentheses. Values closer those of the reference indicate better resalts.
4.3 VELOCITY MODEL ARCHITECTURE
BioMD is generative model that operates directly on all-atom Cartesian coordinates. In contrast to approaches that rely on internal coordinates such as coarse gaase-grained backbones and torsion angles, our method directly models all atoms, enabling to capture subtle structural variations that are critical for realistic biomolecular dynamics. The effectiveness of this all-atom modeling strategy has been demonstrated by state-of-the-ar biomolecular structure models like AlphaFold3 (Abramson et 2024). Notably, our unified model architecture is capable of performing both the forecasting and interpolation tasks (subsec. 4.2.1 and 4.2.2) within the same framework.
Our velocity model architecture is specifically tailored for generating trajectories from single initial structure (Figure 3). The model first employs an SE(3) Graph Transformer to encode the initial conformation, creating rich single and pair representations. Subsequently, our core generative module, the F FlowTra ject ooacotorrrrraa Trans former (Algorithm 6). operates on the entire trajectory sequence. To effectively capture complex biomolecular dynamics, each block of this transformer incorporates two primary attention mechanisms: Attent on rBias responsible for modeling intra- frame spatial interactions, while femporalAt:s dependencies by focusing on the same atom or token across different time steps. By stacking these two attention mechanisms, the model can simultaneously process spatial and temporal information, which is crucial for accurate trajectory prediction.
ion specifically addresses frame- temporal
4.4 AUXILIARY LOSSES
addition to the primary flow-matching objective, we incorporate several auxiliary losses to improve the physical plausibility of the generated structures. These losses are applied to the final predicted coordinates, which are obtained using the model's output velocity field.
Ligand Bond Loss: To preserve the ligand's local structure, we introduce a bond loss following 3
AlphaFold 3 (Abramson et al. 2024). For each bonded atom pair in the ligand, compute the mean squared error between the predicted inter-atomic distance and its ground-truth value, ensuring that the generated ligand structure maintains correct bond lengths.
To ensure physical plausibility and prevent steric clashes, we implement a collision loss that applies squared penalty to non-bonded atom pairs that unrcalistically close. This loss operates on both protein-ligane and intra-ligand interactions, and penalizes inter-atomic distance that falls below predefined threshold.
Collision Loss:
To penalize unrealistic rigid-body movements of ligands, we define geometric center loss. This loss calculates the mean squared error between the geometric center of the predicted ligand atoms and that the ground-truth ligand atoms, penalizing large and unrealistic movements of the entire molecule.
Ligand Geometric Center Loss:

---

Published as a conference paper at ICLR 2026
BioMD
MD
b
BioMD
MD
Ligand torsions
Ligand torsions
Figure 4. Conformation ensemble on the MISATO test set. A comparison of the distributions of conformations ligand torsion angles generated by BioMD MD simulation 6DGE (a) and 3FCF (b).
EXPERIMENTS
We evaluate BioMD three MD trajectory datasets: the MISATO Dataset (Siehenmorgen al. 2024), which comprises protein-ligand interaction trajectories focusing on ligand movement within the protein binding pocket; the DD-13M Dataset (Li et al. 2025), which contains trajectories ligand unbinding from protein binding pockets and ultimately reaching the protein surface; the ATLAS dataset (Vander Meersche et 2024). which contains 100 ns simulations for 1390 protein chains. Examples of predicted trajectories can be obtained from Zenodo.
To comprehensively evaluate our model's performance in generating all-atom biomolecular trajec- tories, we first evaluate the physical stability of the generated structures for both MISATO and DD-13M. enable direct comparison with methods that don't consider protein flexibility, physical stability metrics computed for all heavy atoms of ligands. For the MISATO dataset, given that this dataset provides conformational ensembles, we further evaluate our model's ability to predict conformational flexibility of both proteins and ligands. For the DD-13M ligand unbinding dataset, we also included several metrics to assess the accuracy of the predicted unbinding pathways. For experiments the Atlas dataset, (2024a). In this paper, we compare BioMD with several established ML methods, including methods developed for protein-ligand MD (e.g. NeuralMD (Liu al.) 2024)) and methods developed for protein conformation sampling (e.g. ESMFLow Jing et al. (2024a)). We also include a Static model a baseline, where the initial conformation of the system is held constant throughout the entire trajectory. For the ATLAS dataset, results of Str2Str (Lu et al. 2024a) are obtained using the released code from the official repositories with T3 set to 0.10 0.15 following their instructions, and evaluation results on other methods are provided by Lu et al.
evaluate all methods using metrics provided by Jing al.)
(2025b).
5.1 RESULTS MISATO
To evaluate BioMD's 's ability to generate realistic protein-ligand interaction trajectories, first conduct experiments on the MISATO dataset, which focuses on ligand dynamics within the protein binding pocket. MISATO comprises nearly 20,000 protein-ligand interaction trajectories, each con- taining 100 frames sampled from an 8 ns MD simulation. We evaluate all methods on the MISATO dataset using 1,031 targets, each with a protein sequence length < 800 and ligand size of < 100 heavy atoms. As shown in Table and Figure 7 both variants of our model, BioMD-rel and BioMD-abs, produce trajectories with promising physical stability. The bond and angle geometry errors closely approach the values of the static input structure, and the steric clash scores are orders of magnitude lower than all competing models. Further validation on relaxation consistency (Figure and pocket-ligand interactions (Figures 9 10) confirms the effectiveness of BioMD to generate physically plausible structures.
http://doi.org/10.55811 /2eoodo 16979768
7

---

Published as a conference paper at ICLR 2026
Input
BioMD
Metadynamics
Figure 5. Ligand unbinding path on 6EY8. (a) The input conformation. (b) The unbinding pathways generated by BioMD, with an unbinding path RMSD of 0.24 A(Metalynamics=0.22 Å). The novel pathway discovered by BioMD is highlighted in green. (c) The reference unbinding pathways obtained using metadynamics simulations.
Table 2. Results on the DD-13M test set. Comparison of methods on physical stability (first six metrics). ligand unbinding path reconstruction metric (Unbinding Path RMSD), and ligand unbinding success rates, Mean values on the test samples are reported.
Bond Geometry'
Angle Geometry"
Steric Clashes
Unhinding Path
Unhisding Success
Method
MAE
MSE
MAE
MSE
Intra
Prot-Lig
RMSD
@1
655
@10
Static
2778
.6544
Metadynamics
0246
.0112
.0452
.0030
2277
.217
BioMD-rel
2943
.004
.645
.0029
.0477
.0994
BioMD-abs
0369
.0026
.0545
.0061
.2911
.0003
.6802
.0176
.440
.0888
BioMD-rel (AR-5)
0580
.0100
0918
.0184
4021
.375
.7055
.7088
9295
9706
BioMD-abs (AR-5)
.0288
.0111
.0822
.0132
2943
0009
.5645
.676
.7199
.7911
Bond geometry (bond length) and unbinding path RMSD values are angstroms (A), and angle goometry (bond angle) is in radians. "The metadynamics trajectory serves as the lower-bound. The metrics are calculated among trajectories of multiple repeating simulations.
In terms of conformational flexibility, BioMD demonstrates superior ability to capture the system's ss dynamic behavior. We measure Pearson's correlation between the Root Mean Square Fluctuation (RMSF) of our generated trajectories and the reference MD trajectories. BioMD achieves the highest correlation score for ligand atoms, outperforming NeuralMD by 42.8%. Besides, BioMD achieves the correlation score of 0.685 for protein atoms, while other comparing methods fail to simulate protein conformation changes. Visual analysis in Figures 4 further corroborates these findings, showing that BioMD's predicted atomic fluctuations closely trace the ground truth profiles and that the generated conformational ensemble is qualitatively similar to that of traditional MD simulation. Collectively, these results indicate that BioMD can accurately simulate the flexibility of the entire protein-ligand complex.
5.2 RESULTS ON DD-13M
We further evaluate BioMD on the more challenging task of ligand unbinding using the DD-13M dataset, which comprises 26,612 dissociation trajectories across 565 complexes, each with an average of 480 frames. 36 complexes were held out as a test set for evaluation, while the remaining were used for training. As DD- is generated via metadynamics, the task focuses on reproducing sampling pathways, which do not necessarily represent true thermodynamic or kinetic behavior. A key advantage of our architecture is flexibility in supporting multiple generation strategies. concurrent denoising of all future frames, used on MISATO, results in minimal ligand movement because the model lacks historical guidance and averages over many potential paths. To overcome this, generate the trajectory nuto-regressively. which breaks the long-range prediction into steps and uses previously generated frames to help predict subsequent ones.
The results, summarized in Table 2 highlight the effectiveness of this regresseeee strategy. While maintaining high physical stability. the BioMD-abs (AR-5) model significantly improved path accuracy, reducing the Unbinding Path RMSD to 0.5645 Á and achieving high unhinding correct rate
8

---

Published as a conference paper at ICLR 2026
Table 3. Results on the ATLAS test set. Performance comparison of methods for modeling protein dynamics across different metrics. Bold values indicate the best performance for each metrie.
ESMFL Lw-MD
Metries
Disiiled
ConfDiff
BioEmu
MDGen
BioMI
Pairwise RMSD r 1
0.19
0.19
0.59
0.46
0.23
0.48
0.62
0.70
Predicting flexibility
Global RMSF
0.31
0.33
0.67
0.57
0.38
0.50
0.71
0.76
Por-target RMSF r
0.76
0.74
0.85
0.71
0.57
0,71
0.90
0.91
Root mean W2-dist,
3.60
4.23
2.76
4.32
4.05
2.09
2.43
2.18
+
Trans. contrib. 4
3.13
3.75
2:23
4,04
3.73
2.03
1.89
Distributional accuracy
+. Var. contrib, 4
1.74
1,90
1,40
1.77
1.43
1.20
1,10
MD PCA Wa-dist. 4
1,51
1.87
1.44
1,97
2.04
1.89
1.19
1.24
W2-dist. 4
3.79
3,98
% PC-sim >> 0.5
26
33
35
51
12
44
46
Weak contacts J
0.55
0.48
0.59
0.33
0.43
0.51
0.65
0.61
Ensemble ohservables
Transicut contacts J
0.34
0,30
0.36
0.41
0.46
Exposed residue 1
0.43
Exposed MI matrix
0.20
0.16
0.24
0.07
0.21
0.36
0.34
(Table()) Most importantly, the AR strategy enabled the successful generation of complete unbinding events. The BioMD-rel (AR-5) model achieved remarkable unbinding success rate, identifying valid path in 70.9% of cases with a single attempt @ 1), increasing to 97.1 % with ten attempts (@ 10). This demonstrates BioMD's reliability in exploring critical biomolecular pathways.
On the qualitative analysis for the 6EY8 system (Figure 5 our model not only reproduced the two distinct unbinding pathways found by metadynamies simulations with high fidelity but also discovered a novel third pathway, highlighting the exploratory power of our generative approach. Furthermore, BioMD achieves this with remarkable computational efficiency. While metadynamics required 2654 steps (approx. 1 hour on a single GPU) to find the first path, our model generated complete path in under 10 seconds using just 50 coarse- coarse-grained steps. Additionally, an out-of-distribution (OOD) case with low sequence similarity (Figure 12). BioMD successfully covers multiple escape clusters with an unbinding path RMSD of 0.77 Á (Metadynamics = 0.65 Å).
5.3 RESULT ON ATLAS
We finally used the Atlas dataset (Vander Meersche et al. 2024) to evaluate the model's ability to capture long-timescale protein dynamics with BioMD-rel. The Atlas dataset comprises 1,390
single-chain targets, each associated with three independent 100 ns Molecular Dynamics (MD) trajectories. We evaluated its performance on the Atlas test set, which comprises 82 targets. For each target, we sampled three 100 ns trajectories, where the initial frame for each trajectory was randomly selected from the corresponding MDD reference trajectories. From the resulting 300 of generated trajectories, we extracted 250 frames for evaluation. As shown in Table BioMD achieves of 13 metrics, demonstrating its superiority capturing both structural flexibility and distributional accuracy. Compared to MDGen, which shares the same input setting (sequence initial frame), BioMD exhibits substantial improvements across all metrics. Notably, we observe a 52% increase in the Global RMSF correlation coefficient (r), highlighting the effectiveness of our all-atom architecture in modeling residue-level fluctuations. Furthermore, BioMD surpasses EBA (Lu et al.) 2025a), the leading sequence-based method, on the majority of distributional metrics. This confirms that conditioning on the initial structure provides critical geometric guidance for producing high- tigh-fidelity and dynamically relevant ensembles.
state-of-the-ar performance
Beyond statistical metrics, BioMD successfully reproduces functional conformational transitions. A prime example is the domain motion of Adenylate Kinase (Figure 6). Starting from the closed state (IAKE, 0 ns), BioMD simulates the opening process, reaching a conformation at 23 ns that closely aligns with the experimental open state (4AKE) with an RMSD of 1.631 Á. This indicates that our model can capture major protein conformational transitions that are biologically significant.
To further analyze the exploration capability, we projected the generated conformations onto a 2D space using Time-lagged Independent Component Analysis (TICA). As shown in Figure 13, BioMD explores the conformational space more sufficiently than tiot--iimeeale MD simulations (1 ns and 10 ns), capturing diversity closer to long MD trajectories (100 ns). Remarkably, BioMD generates full 100 ns trajectory in approximately 56 seconds, offering orders-of of-mgnnitddee acceleration over traditional MD simulations.

---

Published as a conference paper at ICLR 2026
open
BioMD
IAKE (closed)
4AKE (apen)
BioMD 10
BioMD
closed
Input
RMSDCRAKE, HuMD 23Ta) 1633
Figure Case study the adenylate kinase. BioMD simulates the domain motion of the adenylate kinasc from close (0 ns) to open (23 ns). The test protein is not included in the training set.
5.4 ANALYSIS
The success of the auto-regressive (AR) strategy in modeling long-range dynamics simultaneously exposes a fundamental challenge in generative trajectory modeling: the error accumulation problem. As shown in Table I 2 and Figure 7 while the non-AR models produce local geometries with errors comparable to the metadynamics reference, the AR models exhibit notable increase in error. However, thanks our hierarchical framework, these errors remain manageable. The bond and angle MAEs for our AR models remain below 0.1 Á and 0.1 radians, respectively- threshold well within the range of thermal fluctuations for molecular systems. These geometrical errors can be readily corrected via a simple local refinement step with minor structural deviations (< 0.1 Å), similar to the relaxation procedure used in AlphaFold. In contrast, non-hierarchical approaches are trapped between two failure modes: large AR steps yield nearly static trajectories, while small AR steps cause significant етrоr accumulation that results in physically unrealistic structures.
Our results also reveal a distinct trade-off between predicting relative coordinate changes (BioMD-rel) and absolute coordinates (BioMD-abs). The absolute coordinate prediction method (BioMD-abs) demonstrates superior grasp of the global conformational landscape, evidenced by its higher protein RMSF correlation on MISATO and a more accurate centroid path RMSD on DD-13M, making it the preferred choice for tasks requiring the precise reproduction of specific dynamic pathways. In contrast, the relative coordinate prediction method (BioMD-rel) excels at encouraging more exploratory behavior while preserving local chemical fidelity. Its strength is highlighted by the significantly higher unbinding success rate on DD-13M, which makes it more effective for applications focused on sampling large-scale conformational changes and discovering novel dynamic events. This functional duality means BioMD can be flexibly adapted to the specific goals of a simulation, whether the priority is accuracy in reproducing known dynamics or exploration discover new ones.
CONCLUSION
In this work, we introduce BioMD, a novel all-atom generative model that overcomes the com- putational limitations of traditional molecular dynamics simulate long-timescale biomolecular events. Our hierarchical framework, which synergistically combines coarse-grained forecasting with fine-grained interpolation, effectively mitigates error accumulation and enables the generation of physically realistic trajectories. We demonstrated BioMD's capabilities on two challenging datasets, showing it can produce stable conformations that accurately capture protein-ligant flexibility on the MISATO dataset and successfully generate complete ligand unbinding pathways for up to 97.1 1% of systems on the DD-13M dataset. Notably, BioMD achieves this with remarkable computational effi- ciency, identifying unbinding paths in seconds compared to the hours required by traditional methods like metadynamics. By offering distinct modes optimized for either accurate pathway reproduction or broad exploratory sampling, BioMD provides powerful, flexible, and efficient tool poised to accelerate research computational chemistry and drug discovery. Nevertheless, BioMD's ability to generalize to substantially longer trajectories (uus ms) rare events beyond the training distribution remains limited, and we highlight this as an important direction for future work.
"This paper written with assistance from large language models (LLM) for proofreading and polishing.
10

---

Published as a conference paper at ICLR 2026
ACKNOWLEDGMENTS
This project was supported by Shenzhen Hetao Shenzhen-Hong Kong Science and Technology Innovation Cooperation Zone, under Grant No. HTTHZ0SW--KCYB--000000
REFERENCES
Josh Abramson, Jonas Adler, Jack Dunger, Richard Evans, Tim Green, Alexander Pritzel, Olaf Ronneberger, Lindsay Willmore, Andrew Ballard, Joshua Bambrick, et al. Accurate structure prediction of biomolecular interactions with alphafold 3. Nature, 630(8016)493-500,, 2024.
Stewart A Adcock and J Andrew McCammon. Molecular dynamics: survey of methods for simulating the activity of proteins. Chemical reviews, 106(5):1589-166 2006.
Alessandro Barducci, Massimiliano Bonomi, and Michele Parrinello. Metadynamics. Wiley Interdis- ciplinary Reviews: Computational Molecular Science, 1(5):826-843, 2011.
Boyuan Chen, Diego Martí Monsó, Yilun Du, Max Simchowitz, Russ Tedrake, and Vincent Sitzmann. Diffusion forcing: Next-token prediction meets full- full-sequence diffusion. Advances in Neural Information Processing Systems, 177:24811-24111 2024.
Allan dos Santos Costa, Ilan Mitnikov, Franco Pellegrini, Ameya Daigavane, Mario Geiger, Zhonglin Cao, Karsten Kreis, Tess Smidt, Emine Kucukbenli, and Joseph Jacobson. Equijump: Protein dy- namics simulation via so (3)-equivariant stochastic interpolants. arXiv preprint arXiv:2410.09660 2024.
Thomas A. Darden, Darrin M. York, and Lee G. Pedersen. Particle mesh ewald: An log(n) method for ewald sums in large systems. Journal of Chemical Physics, 98:10089-10092, 1993, URL
Diego del Alamo, Davide Sala, Hassane S Mchaourab, and Jens Meiler. Sampling alternative conformational states of transporters and receptors with alphafold2. elife, 11:775751, mar 2022. doi: 2t.7554/eLife.757555 URL http:s//oo..org11075544.........
Ron O Dror, Robert M Dirks, JP Grossman, Huafeng Xu, and David E Shaw. Biomolecular simulation: simmulaonnasiimulaaionn a computational microscope for molecular biology. Annual review of biophysics, 41:429 2012.
Scott Hollingsworth and Ron O Dror. Molecular dynamics simulation for all. Neuron, 99(6): 1129-1143, 2018.
Bowen Jing, Ezra Erives, Peter Pao-Huang, Gabriele Corso, Bonnie Berger, and Tommi Jaakkola. Eigenfold: Generative protein structure prediction with diffusion models. arXiv preprint arXiv:2304.02198, 2023.
Bowen Jing, Bonnie Berger, and T. Jaakkola. Alphafold meets flow matching for generating pro- tein ensembles. ArXiv, abs/2402.04845, 2024a. URL httpp::/appp.emmatttcchooooooooe ortparrooo/////lll/i
267522949
Bowen Jing, Bonnie Berger, and Tommi Jaakkola. Alphafold meets flow matching for generating protein ensembles. In Proceedings of the 41st International Conference on Machine Learning, pp. 22277-22303. 2024b.
Bowen Jing, Hannes Stärk, Tommi Jaakkola, and Bonnie Berger. Generative modeling of molecular dynamics trajectories. Advances in Neural Information Processing Systems, 37:40534 2024c.
John Jumper, Richard Evans, Alexander Pritzel, Tim Green, Michael Figurnov, Olaf Ronneberger, Kathryn Tunyasuvunakool, Russ Bates, Augustin Židek, Anna Potapenko, al. Highly accurate protein structure prediction with alphafold. Nature, 596(78733:583-558) 2021.
Martin Karplus and J Andrew McCammon. Molecular dynamics simulations of biomolecules. Nature structural biology, 9()):646-66 -6622 2002.
11

---

Published as a conference paper at ICLR 2026
Alessandro Laio and Michele Parrinello. Escaping free-energy minima. Proceedings of the national academy of sciences, 99(20):12562-2256) 2002.
Sarah Lewis, Tim Hempel, José Jiménez-Luna, Michael Gastegger, Yu Xie, Andrew YK Foong, Victor Gareía Satorras, Osama Abdin, Bastiaan S Veeling, Iryna Zaporozhets, et al. Scalable emulation of protein equilibrium ensembles with generative deep learning. Science, pp. cadv9817, 2025.
Maodong Li, Jiying Zhang, Bin Feng, Wenqi Zeng, Dechin Chen, Zhijun Pan, Yu Li, Zijing Liu, and Yi Isaac Yang. Enhanced sampling, public dataset and generative model for drug- protein dissociation dynamics, arXiv preprint arXiv:2504.18366 2025.
Yaron Lipman, Ricky T. Q. Chen, Heli Ben-Hamu, Maximilian Nickel, and Matthew Le. Flow matching for generative modeling. In ICLR, 2023.
Shengchao Liu, Weitao Du, Yanjing Li, Zhuoxinran Li, Vignesh Bhethanabotla, Nakul Rampal, Omar Yaghi, Christian Borgs, Anima Anandkumar, Hongyu Guo, et al. A multi multi-grained sym- metric differential equation model for learning protein-ligand binding dynamics. arXiv preprint 15X1: 2024.
Jiarui Lu, Bozitao Zhong, Zuobai Zhang, and Jian Tang. Str2str: A score-based framework for zero-shot protein conformation sampling. In The Twelfth International Conference on Learning Representations, 2024a.
Jiarui Lu. Xiaoyin Chen. Stephen Zhewen Lu. Aurelie Lozano, Vijil Chenthamarakshan, Payel Das, and Jian Tang. Aligning protein conformation ensemble generation with physical feedback. In Forty-second International Conference on Machine Learning. 2025a. URL. https:://opennrview.neei
forum"id==Ast955jaul
Jiarui Lu, Xiaoyin Chen, Stephen Zhewen Lu, Aurélie C. Lozano, Vijil Chenthamarakshan, Payel Das, and Jian Tang. Aligning protein conformation ensemble generation with physical feedback. ArXiv, abs/2505. 24203,, 2025b. URL
Wei Lu, Jixian Zhang, Weifeng Huang, Ziqiao Zhang, Xiangyu Jia, Zhenyu Wang, Leilei Shi, Chengtao Li, Peter G. Wolynes, and Shuangjia Zheng. Dynamicbind: predicting ligand-speciff protein-ligand complex structure with deep equivariant generative model. Nature Communica- lions, 15, 2024b. URL
Frank Noé. Simon Olsson, Jonas Köhler, and Hao Wu. Boltzmann generators: Sampling equilibrium states of many-body systems with deep learning. Science, 365(645)) 147, 2019.
David E Shaw, Ron O Dror, John K Salmon, JP Grossman, Kenneth M Mackenzie, Joseph A Bank, Cliff Young, Martin M Deneroff, Brannon Batson, Kevin Bowers, al. Millisecond-scale molecular dynamics simulations on anton. In Proceedings of the conference on high performance computing networking. storage and analysis, pp. 1-11, 2009.
David E Shaw, Paul Maragakis, Kresten Lindorff-Larsen Stefano Piana, Ron O Dror, Michael P Eastwood, Joseph A Bank, John M Jumper, John K Salmon, Yibing Shan, et al. Atomic-level characterization of the structural dynamics of proteins. Science, 33066002):411344. 2010.
Yuning Shen, Lihao Wang, Huizhuo Yuan, Yan Wang, Bangji Yang, and Quanquan Gu. Simul- taneous modeling of protein conformation and dynamics via autoregression. arXiv preprint arXiv:2505.174: 17478, 2025.
Till Siebenmorgen, Filipe Menezes, Sabrina Benassou, Erine Merdivan, Kieran Didi, André San- tos Dias Mourão, Radoslaw Kitel, Pietro Liò, Stefan Kesselheim, Marie Piraud, et al. Misato: machine learning dataset of mooeein-ligand complexes for structure-based drug discovery. Nature Computational Science, Science,,, pp. 1-12, 2024.
Richard Stein and Hassane S Mchaourab. Speach af: Sampling protein ensembles and conforma- tional heterogeneity with alphafold2. PLoS computational biology, 1810)2. 2022.
12

---

Published as a conference paper at ICLR 2026
Yann Vander Meersche, Gabriel Cretin, Aria Gheeraert, Jean- -han-Chrsttophe Gelly, and Tatiana Galochk- ina. Atlas: protein flexibility description from atomistic molecular dynamics simulations. Nucleic acids research, 52(D1):D384-D33 2024.
Loup Verlet. Computer "experiments" on classical fluids. i. thermodynamical properties of lennard- jones molecules. Physical Review, 159(1))98-13,, 1967.
Tong Wang, Xinheng He, Mingyu Li, Yatao Li, Ran Bi, Yusong Wang, Chaoran Cheng, Xiangzhen Shen, Jiawei Meng, He Zhang, al. Ab initio characterization of protein molecular dynamics with 2024a.
ai2bmd. Nature, 13580040):1019-0033
Yan Wang. Lihao Wang, Yuning Shen, Yiqun Wang, Huizhuo Yuan, Yue Wu, and Quanquan Gu. Pro- tein conformation generation via force-guided se (3) diffusion models. In Forty-first International Conference Machine Learning. 2024b.
Hannah K Wayment-Steele, Adedolapo Ojoawo, Renee Otten, Julia M Apitz, Warintra Pitsawong. Marc Hömberger, Sergey Ovchinnikov, Lucy Colwell, and Dorothee Kern. Predicting multiplc conformations via sequence clustering and alphafold2. Nature, 625(7996):832 829,,,9 2024.
Shuxin Zheng, Jiyan He, Chang Liu, Yu Shi, Ziheng Lu, Weitao Feng, Fusong Ju, Jiaxi Wang. Jianwei Zhu, Yaosen Min, et al. Predicting equilibrium distributions for molecular systems with deep learning. Nature Machine Intelligence, 6(5):558-567) 2024.
Xiangxin Zhou, Yi Xiao, Haowei Lin, Xinheng He, Jiaqi Guan, Yang Wang, Qiang Liu, Feng Zhou, Liang Wang, and Jianzhu Ma. Integrating protein dynamics into structure-based drug design via full-atom stochastic flows. ArXiv, abs/2503.03988) 2025. URL butps:/appi.smmmntiscclola.rrrrr 10rpusID:276813705
13

---

Published as a conference paper at ICLR 2026
A TECHNICAL APPENDICES AND SUPPLEMENTARY MATERIAL
A.1 DETAILED MODEL ARCHITECTURE
Hierarchical Generation Framework. As illustrated in Figure BioMD employs a hierarchical framework to perform both coarse-grained forecasting and fine-grained interpolation within a unified model. The specific task is controlled by applying noise selectively. For Forecasting, the initial frame Xy is provided without noise, while all subsequent frames are initialized from standard Gaussian distribution. For Interpolation. two anchor frames (e.g., and X21.) are kept clean, while the intermediate frames are initialized from noise. The model's 's objective to denoise the masked frames conditioned on the known ones.
Input Representation and Conditioning. The core of the model is the FlowModule (Algorithm 4). which processes three primary inputs. The main dynamic input is the set of Noisy Coordinates ((x) }), representing the current state of the trajectory. To provide structural con- text, the initial conformation (Frame 0 conf.) is processed by SB(3)-Graph Transformer, as detailed in the main inference loop (Algorithm D This produces static Single ((s)runk) }) and Pair ([2 }) representations. These representations, along with other atom features, are processed by the F FlowCondi t tioning module (Algorithm 5) generate the final conditioning signals.
Attention Pathway. The FlowModule uses a Docal-gioballlocat attention pathway to predict the velocity field. First, the noisy coordinates and conditioning features are passed to an At omAt tent i onEncoder, which models local atomic environments. The resulting repre- sentations are aggregated into tokens and fed into the central F FlowTra jectarycttrry Trans former (Algorithm 6). This module integrates spatial and temporal information using two key mechanisms: Attent i ConPairBiii captures inter-frame dynamics. The globally-aware token representations are then broadeast back to the atomic level, where At ttemAAtteniionoan
resolves intra-frame spatial relationships, while cemporalAtter i on
computes the final per-atom updates.
module is key component of our architecture, designed to capture dynamic dependencies across the time dimension of the trajectory. As shown in Algorithm 7 it operates the single representation, Ca. by treating the time axis the sequence length for the attention mechanism. This allows the model to integrate information from all frames of the trajectory for each residue or atom.
Ateeetioltttnnoon Module. The TemporalAttent
AttentionPairBias Module. The IntentionPaiiiie comes from AlphaFold. It projects the pair representation into bias term. This term is added di- rectly to the attention logits before the softmax operation, effectively steering the attention mechanism to focus on spatially relevant residue pairs.
module (Algorithm g) originally
Velocity Field Prediction and Trajectory Generation. The output of the FlowModule is the Flow vector field ((ū+)), which represents the predicted velocity for each atom. During training (Algorithm 2). the model is optimized via mean squared error loss between the predicted velocity and the true velocity. During inference (Algorithm 3). this vector field is used in an Euler integration step, x xi di ui, iteratively update the coordinates from noisy state fina trajectory.
A.2 AUXILIARY LOSSES
After we get the estimated vector field ug. we can get the predicted structure coordinates via
x + u T)2
and then we get the predicted protein and ligand structure XxP x = x]
Ligand geometric center loss. To stabilize the global placement of the ligand and prevent spurious rigid translations, we align the predicted and reference geometric centers of ligand atoms. Let
x [x.+). Jis and xf = ( Ji=1 denote ground-truth and predicted ligand coordinates at step I.
14

---

Published as a conference paper at ICLR 2026
The geometric center is
N₁ c(x( =
N₁
C(x))
N₂ im1
and the loss is the mean-squared discrepancy
Leenter ||C(x)) c(x)
This term softly anchors the ligand's global position while remaining agnostic to its internal geometry.
Collision loss. To penalize steric clashes, we define a collision loss between protein- ligand atoms and within ligand atoms. Let x and xf denote ligand and protein atom coordinates at step and x x? their predictions. We compute predicted distances
d' = P 12,
dit - 1|x)
and corresponding ground-truth minimal distances
dt min de
| p 1l.,
di off min |x/ x; 112
ligaad and ligand figandliggndd thresholds are set as PT. min 0.9 off..stt Spt) where Sar = 3.0 Å and Çu = 2.0 Å.
Sij min (0.9 Su).
The collision loss is then defined as
Coollision -E < 2 L dpp/2
11 ( (1-bb) (c/ dz.ly2
1
where 1(-) represents the indicator function and bij is the ligand bond mask to exclude bonded pairs.
Ligand bond loss. To preserve ligand bond lengths, we penalize deviations between predicted and ground-truth bonded atom distances. Let B denote the set of bonded atom pairs according to the ligand bond mask. For each bond (i, E B, we compute the predicted and ground-truth distances
K | 12. di |xx
f.j 1la.
The bond loss is then defined as the mean squared deviation:
1
Σ
( 2
B
(KJ)EB
Geometric constraint loss. We combine the above terms into a single geometric regularizer
Lyoom
Xenl Gcollision + Abond Lhond + Aetr Leenter,
where Acol, Aboud, Aetr > 0 balance steric clash avoidance, bond-length preservation, and global ligand anchoring, respectively.
A.3 EVALUATION METRICS
A.3.1 PHYSICAL STABILITY
This metric assesses whether the generated trajectories preserve physically stable conformations, which is essential to ensure chemical validity and avoid unrealistic molecular structures. We evaluate stability from two complementary perspectives:
1. Local Structure Stability. To assess whether the generated trajectories maintain chemically reasonable local geometries, we calculate the deviations of bond lengths and bond angles with respect the initial frame of the reference trajectories. Both the Mean Absolute Error (MAE) and Mean Squared Error (MSE) are reported. Lower values indicate that the generated conformations remain close to the idealized covalent structure and are thus more chemically stable.
15

---

Published as a conference paper at ICLR 2026
2. Steric Clashes. We further quantify the presence of steric conflicts, which occur when non-bonded atoms are unrealistically close to each other. Specifically, a clash is counted if the interatomic distance (excluding bonded pairs and angle-related atoms) is less than a threshold of 1.5 A. We compute clash scores for both intra-ligand and protein-ligand interactions, where the score corresponds to the average number of clashes per generated conformation. Lower clash scores indicate physically more plausible conformations.
A.3.2 CONFORMATIONAL FLEXIBILITY
In addition to stability, it is important that generated trajectories capture the dynamic flexibility of molecular systems. For the MISATO protein-ligand interaction dataset, we adopt the Root Mean Square Fluctuation (RMSF) to quantify the extent of atomic motion over time after trajectory alignment:
1 T Σ (m(()) - F.l1,,
RMSF;
where (t) is the position of atom at time t, and r, is its time-averaged position.
We evaluate flexibility from two perspectives: 1. Global Consistency. We compute the Pearson correlation coefficient between the RMSFs of generated and reference trajectories, where higher correlation indicates better agreement the fluctuation profiles. Magnitude Accuracy. We also report the average RMSF of the generated trajectories. Values closer to the reference average RMSF imply that the model produces realistic levels of conformational motion rather than being overly rigid or excessively flexible.
A.3.3 UNBINDING PATH DISTANCE
For the DD-13M ligand unbinding dataset, we evaluate whether generated unbinding trajectories follow realistic spatial pathways compared to reference simulations. We compute the Root Mean Square Deviation (RMSD) between generated and reference ligand centroid trajectories with the following procedure:
Trajectory Standardization. All ligand centroid trajectories are resampled uniform length (L = 100 frames) using linear interpolation, ensuring comparability between different sequences.
Best-Match Search. For each generated trajectory, we identify the reference trajectory that yields the minimum RMSD. This accounts for the possibility of multiple plausible unbinding pathways. Let Cgen = (Ci,..., C denote the sequence of ligand centroids in a generated trajectory, and let Tret (c(4) (k) K be the set of K reference trajectories. The distance metric calculated as:
f. c
RMSDrest min
Σ
|c
(k)
(
To ensure a fair comparison for the static baseline, we construct pseudo-trajectory by replicating the initial bound pose 7. times (i.e., C₁ for all t).
3. Final Score. The reported metric is the average of these best-match RMSDs across all generated trajectories. Lower RMSD values indicate that the model generates ligand motions more consistent with physically realistic unbinding paths.
A.3.4 UNBINDING SUCCESS
This metric evaluates whether the generated ligand trajectories successfully capture the unbinding event. Specifically, construct the convex hull of the protein heavy atoms in the initial bound state. If at least one predicted ligand centroid position lies outside this convex hull, the trajectory is considered successful unbinding case.
We report the Success @ k, which measures the probability that at least one out of l independently generated trajectories for the same protein-ligand complex achieves successful unbinding. A higher
16

---

Published as a conference paper at ICLR 2026
success rate indicates better capability of the model to reproduce realistic ligand unbinding processes. Formally, for each complex with k attempts, Success @k is defined as
N
1 Success @ k = N
Σ 7=
( =1]
max ISiC
where x/22 is the binary success indicator (1 if the j-th trajectory of complex n achieves unbinding, O otherwise), and N is the total number of complexes. We report Success@1 1, Success@5. and @uccess@ 10, which reflect performance under varying generation attempts, respectively.
A.3.5 UNBINDING CORRECT RATES
To provide more rigorous evaluation that directly addresses whether model identifies the correct exit pathway, we perform an analysis on the successful unbinding trajectories. We define the ground-truth exit position, e the ligand's centroid coordinates the final frame of the reference simulation. We introduce two metrics to quantify the spatial accuracy of the generated exit paths:
Endpoint Distance. This metric measures how closely the generated pathway approaches the known exit location. For a generated trajectory consisting of L frames, denoted as Cgen (ci.....Cl) CL), we calculate the minimum Euclidean distance between any centroid position along the predicted path and the ground- exit position:
D Dend
min
||c,
Cexit
GT
(E
We report the mean of this minimum distance across test samples. Lower values indicate that the generated trajectory passes closer to the true exit point.
Correct Rate. To provide an intuitive measure of success, we calculate the percentage trajec- tories that successfully navigate within specific proximity of the correct exit. Given distance threshold 7 (e.g., 0.5 Á 1.0 À), the correct rate is defined as:
Rate(T)
1
I(D
<T)
N
where N is the total number of samples, D is the endpoint distance for the i-th sample, and I(-) the indicator function.
A.4 EXPERIMENTAL DETAILS
A.4.1 SAMPLING SETTINGS
For the first caars--grained forecasting stage, set the sampling step k 10 for all three datasets, with different underlying simulation time intervals:
(i) For the MISATO dataset, the all-atom MD simulations were performed with 80 ps timestep. Therefore, our sampling of every & = 10 steps corresponds to a physical time interval of 0.8 ns (10 steps x 0.08 ns /step).
(ii) For the DD-13M dataset, the simulations utilized 0.1 ps timestep. Consequently, sampling every k=10 steps resulted in a time interval of ps between saved frames (10 steps 0.1 ps/step). For two case studies (6EY8 8 and 3PCU), we sampled 100 trajectories to compute the unbinding correct rate metric.
(iii) For the Atlas dataset, each trajectory has total simulation length of 100 ns and is provided as sequence of 10,000 frames. Consequently, sampling every k=10 steps resulted in a time interval of ns between saved frames (10 steps x 100 ps/step).
A.4.2 HYPERPARAMETERS
We use NVIDIA RTX A6000 GPUs for training. For each dataset, we use the Adam opti- mizer and train BioMD for 20,000 steps with fixed learning rate of 0.0001 and batch size of At AtomAttent i
32. The hyperparameters for
AtomAttent www.eeedee.....e..e......... and GraphTransformer is shown in Table 456
17

---

Published as a conference paper at ICLR 2026
Table 4. Key hyperparameters for the FlovTrajsccerrytraaatooreer model.
Hyperparameter
Symbol
Value
Description
Number of Transformer Blocks
Nylock
24
Number of sequential transformer layers.
Number of Attention Heads
Nbad
16
Number of heads in the multi-head attention.
Single Representation Dim.
Cs
384
Dimension of the per-residue representation.
Pair Representation Dim.
C2
128
Dimension of the residue-pat representation.
Table 5. Key hyperparameters for the At cmAttten con.etttodetinnneeeeeee. onEncode and AtomAtteet
model.
Hyperparameter
Symbol
Value
Description
Number of Transformer Blocks
Volock
3
Number of sequential transformer layers.
Number of Attention Heads
Nhad
4
Number of heads in the multi-head attention.
Single Atom Representation Dim.
Catom
128
Dimension of the per-atom representation.
Pair Atom Representation Dim.
Catompair
16
Dimension of the atom-pair representation.
Table
Key hyperparameters for the Graph rraphhrransfomme model.
Hyperparameter
Symbol
Value
Description
Number of Encoder Layers
L
6
Number of encoder layers.
Encoder Embedding Dimension
H
384
Encoder embedding dimension.
FFN Embedding Dimension
F
1536
Encoder embedding dimension for FFN.
Number of Attention Heads
A
16
Number of encoder attention heads.
Table 7. Ablation study on the MISATO dataset. Comparison of all methods on physical stability (first six metrics) and conformational flexibility (last four metries). Mean values the test samples reported.
Bond Geometry
Angle Geometry'
Sterie Clashes
RMSF Correlation
RMSF Value
Method
MAE
MSE
MAE
MSE
Intra Lig
Prot-Lig
Ligand
Protein
Ligand (1.211)
Protein (1.002)
Molecular Dynamics
0377
.0033
.0575
.0033
0
BioMD (k=1)
0356
.0021
.0568
.0053
0
.4526
.6291
.3717
.3825
BioMD (k=5)
0387
.0044
.0627
.066
0
.0001
,4982
.6233
A552
.4877
BioMD (k=10)
.0955
.0155
.0709
.097
.0019
.0033
.4789
.6854
.7023
.6242
BioMD (k=20)
.0421
.0330
.6988
.0877
.0002
.0008
4267
.5309
.5403
.4999
Bond geometry (bund length) and RMSF values are in angstroms (A). Angle geometry (bund angle) is in radians.
RMSF Correlation reported using the Pearson comelation coefficient.
RMSF values for reference traoectorics given in parentheses. Values closer to those of the reference indicate better results.
Table 8, Ablation study the DD-13M dataset. Comparison of methods on physical stability (first six metrics), ligand unbinding path reconstruction metric (Unbinding Path RMSD), and ligand unbinding success rates. Mean values on the test samples are reported.
Bond Geometry'
Angle Goometry'
Sterie Clashes
Unbinding Path
Unbinding Success
Method
MAB
MSE
MAE
MSB
Intra-Lig
Prot-Lig
RMSD
@1
@5
@10
Statie
.2778
6504
0
Metadynamics
JU446
.0012
.0452
.0030
.2777
4217
BioMD (AR-1. k=10) k01)
1880
.6122
.2011
.0761
1.727
4.079
1.131
.2855
.4288
5822
BioMD (AR-5,
.0111
.5676
.7419
BioMD (AR-10, k=10)
1067
0021
.1478
0450
7098
.5356
7465
.0088
.5719
6239
BioMD (w.o. AR, k=10)
.0308
.0018
.0606
.0077
.9444
.0004
.6845
.0029
.0147
.094
BioMD AA L
.1457
.0382
.1770
.U595
1.407
2.029
.7766
.5735
.7512
.8010
BioMD (AR 5, k=
.682
.0179
.1231
.0321
.3686
7,380
.7313
.3428
.5128
.6233
BioMD (AR-5, k=5)
0575
.0164
.1022
.0246
.4700
1.346
.5033
.5144
.6582
.7033
BioMD (AR-5. k=10)
0728
0111
.8022
.0122
.2943
.0009
.5645
.5676
.7419
.7941
BioMD (AR-5, k=20)
.0053
.3736
.0967
.7065
Bond goometry (bond longth) and unbinding path RMSD values are angstroms (A), and angle geometry (bond angle) is radians. The metadynamics trajectory serves the lower-bound, The metrics caleulated umong trajectories of multiple repeating simulations. "L=22" means that the generated irajectory length is doubled (i.e., 1000 frames).
represents forecasting-only" version of BioMD.
18

---

Published as a conference paper at ICLR 2026
A.5 ABLATION STUDY
We conducted ablation studies using BioMD-rel, generating 8 ns trajectories for MISATO and 500 frames for DD-13M. During forecasting, we employed moving window of the 50 most recent frames as historical context to mitigate GPU memory constraints. The quantitative results, presented in Tables and g demonstrate that the model's performance highly sensitive to both the auto- regressive block size (AR-z), which is the number of concurrently denoised frames during sampling, and the hierarchical step size (k), which the frame interval in the coarse-grained forecasting stage.
We first observe that the auto-regressive block size is a critical determinant of physical stability. A minimal block size (e.g., AR-1) necessitates a high number of sequential inference steps, leading to severe eπrоr accumulation. This accumulation manifests as catastrophic physical violations, with protein-ligand clashes increasing by over 4000 times compared to the baseline. Conversely, an excessively large block size (e.g., w.o. AR) appears to increase the task complexity, confusing the model and degrading predictive accuracy. Consequently, AR-5 emerges as the optimal trade-off, providing sufficient context for stable predictions without introducing unnecessary computational difficulty.
Furthermore, larger step size (k) proves essential for mitigating error propagation. Our experiments show that a non-hierarchical approach (k = 1) suffers from extreme instability due to the vast number of required inference steps, resulting in levels. While the model is generally less sensitive to & than to the AR setting, the results indicate that 10 yields superior efficacy for unbinding tasks. Therefore, the specific combination of AR-5 and k 10 was selected as the final configuration, it effectively balances unbinding metrics with the maintenance of physical realism in complex biomolecular dynamics.
clashes exceeding 7000 times the baseline
A.6 ADDITIONAL EXPERIMENTAL RESULTS
A.6.1 BOND AND ANGLE ERROR
To evaluate the geometric accuracy and temporal stability of our generative model, we analyzed the distributions of bond length and angle errors, as well as the error accumulation problem. All experiments are conducted with BioMD-rel. The results in Figure (soo that on the MISATO test set, 97.9% of bond length errors are below 0.2 A, and 95.3% of bond angle errors below 0.2 radians ( (11 11.5%). The performance on the DD-13M test set also demonstrates 99.7% of bond errors and 98.2% of angle errors falling below these thresholds.
Furthermore, we examined the trend of errors over time. On the MISATO dataset, due to the short generated trajectories (within 100 frames), the errors across different time spans are similar. In contrast, on the DD-13M dataset involving long trajectory generation, as the number of generated frames increases (up to 500 frames), significant increase in error is observed, indicating the occurrence of the eтror accumulation problem.
A.6.2 RELAXATION CONSISTENCY ANALYSIS
To assess the energetic quality of the generated conformations, we conducted a relaxation consistency analysis on them. For each trajectory generated by BioMD-rel on the MISATO test set, we took each frame and performed energy minimization using the Amber99SB force field. We then calculated the all-atom RMSD between the pre- and post-relaxation structures. A low RMSD indicates that our generated conformations are near a local energy minimum, confirming their physical plausibility. Our analysis reveals mean all-atom RMSD of 0.72 À (Figure & after relaxation across all generated frames (0.69 À for MD conformations, Figure 8), demonstrating the ability of BioMD generate energetically favorable structures and providing strong evidence for the physical plausibility of our generated trajectories beyond purely geometric metrics.
A.6.3 POCKET-LIGAND INTERACTION ANALYSIS
To evaluate the accuracy of the predicted binding interactions, we perform an analysis on the MISATO test set (1,031 complexes), focusing on the conformation of key binding residues and their interactions with the ligand. We identified residues located within 5 Á of the ligand and compared
19

---

Published as a conference paper at ICLR 2026
their conformational properties against ground truth MD simulations. All experiments conducted with BioMD-rel. As shown in Figure the chil dihedral angle distributions predicted by BioMD (red dashed line) closely align with those observed in the MD trajectories (blue solid line). Furthermore, we visualized the interaction patterns by calculating the contact probability between protein residues and ligand atoms using a 5 Å cutoff (Figure 10). The resulting contact probability map from BioMD exhibits high degree of consistency with the MD reference, confirming that the model accurately captures both the side-chain geometries and the dynamic contact profiles within the binding pocket.
A.6.4 TICA ANALYSIS
For each protein in the test set, we projected the generated conformations onto 2D space using Time-lagged Independent Component Analysis (TICA). We train the BioMD-rel model the Atlas dataset, which contains 100 ns MD trajectories for 1,390 protein chains. The first two independent components (tICI and tCC2), which capture the slowest and most significant motions of the system, were used to define a 2D conformational space. We first fit TICA on a 100 ns MD trajectory and then projected the conformations generated by BioMD onto this 2D space. Specifically, we used BioMD to sample three 100 ns trajectories and projected all generated conformations onto the TICA space. For comparison, we also projected conformations from other MD simulation replicas (100 ns, 10 ns, and ns) onto the same space. As shown in Figure 13 BioMD explores the sample space more sufficiently than ns and 10 MD simulations, capturing a broader range of conformational diversity. Notably, the average sampling time for one trajectory is only 56 seconds, making BioMD) much faster than MD simulations.
Table 9. Comparison of methods on endpoint distance and unbinding correct rates. Mean values on the test samples are reported.
Endpoiat Distance
Correct Rate (s X Á)
Method
RMSD (A)
<0.55 Á
<1
<2 Á
< Á
Statie
3,5852
0.1765
0.3235
0.4118
0.7353
Metadynamics
0.9641
BioMD-rel
2.3433
0.4265
0.5265
0.6971
0.8000
BioMD-abs
2.3897
0.3765
0.5324
0.6882
0,8059
BioMD-rel (AR-5)
1.8001
0.3471
0.5882
0.7529
0.8706
BioMD-abs (AR-5)
1.7253
0.4912
0.6824
0.7853
0.8588
The metadynamics trajectory trajecturies of multiple repeating simulations.
the lower-bound. The metries calculated among
Bon d and angle enrors MISATO tast
Bond Longth
Bond Angle (radians)
3.6
1.2
0.5
1.0
0.4
0.3
0-20
20-40
40-80
60-20
80-90
0-20
20-40
42-60
60-00
00-90
Trammes Segmennt
Tramees Segmennt
6
Bond and angle errors the DD-13M 13M test set
Bond Length troo (A)
Bond Angle Erroe (radlane)
Uie Enor
0.4
0.8
0.2
0.0
0:100
200
2:00 300
300-400
400 500
0:100
100-200
200 300
300 40-466
400 500
Prames
Prames
Figure 7. Bond and angle error distribution. Baxen plot showing error distribution the MISATO test set. b. Boxen plot showing error distribution on the DD-13M test set.
20

---

Published as a conference paper at ICLR 2026
MD
BioMD
30
Count
20
Count
20
10
10
-
0.5
0.6
0.7
0.8
0.9
1.0
1.1
0.5
0.6
0.7
0.8
0.9
1.0
1.1
RMSD (A)
RMSD (A)
Figure 8. All-atom RMSD after relaxation, All-atom RMSD when using the Amber99SB force field for
relaxation on MD conformations (mean=0.69 Á and BioMID-generated conformations (mean=0.72 Á
Dihedral distributions residues closest ligand (5XG5)
A)
(3.57
c1N1 (5.73 A)
0.04
JS0=2,002
0.2
|:0+0.005
|50+0 235
0.004
150+0.321
0.004
Denety
0.02
0.002
0.00
0.0
0.000
0.000
(3.57
chil (3.59 A)
ASP37
HI539
12,92
Denotty
0.2
Sundity
0.06
Dunally
0.30
1l0-0.152
0.04
0.05
0.02
1.00
-0
100
-100
-20
90
100
1.00
-00
92
100
-130
-00
20
100
ARG41 chil 15.21 A)
ASP125 chil (2.72 A)
0.03
J50m0.:52
0.075
J50+0.317
0.050
0.01
0.00
0.000
distributions residues ligand (6DYO)
PHE25 chil (3.57 A)
GLU32 chil (2.57
TYR44 chil [4.09 A)
LEUS3 -n133-cch111 14.03 At
|50+0
|S0+0.224
366
0.0075
J50=0.207
Dordty
0.04
Dereby
Derake
0.0050
0.02
0.00
0.00
0.0000
1.5
90
150
120
80
180
300
50
350
180
90
00
100
VAL96 chil 13.80 A)
PHE97 chil (4.23
LEU100
chil (3.76
PHE123
chil 13.08 A)
0.004
150+0 323
0.04
150+0.132
0.2
JSD=6.139
J50-0.329
Denstt
Denstty
Dendf
0.04
0.002
0.1
0.02
0.000
0.00
0.0
0.00
-383
180
|50+0.400
0.15
1S0+0.226
BC-0.375
0.004
0.10
0.002
0.05
0.000
0.00
0.0
Figure 9. Dihedral distributions for residues closest to the ligand. The dashed red line represents BioMD, while the solid blue line represents the ground truth MD simulations. The value in parentheses indicates the minimum distance between the considered residue atoms and all ligand atoms. Test cases (6DYO and 5XG5) from the MISATO test set.
21

---

Published as a conference paper at ICLR 2026
MD (8 ns)
BioMD
1.0
25
43
64
Residue
85
103
103
12
15
18
12
15
10
Ugand Atom Index
Upaad Abom Inoex
25
44
Residun
1RGL
55
55
60
69
0.0
18
21
15
18
21
Ligand Atom Index
Ligand Abom Index
23
73
74
688
75
Protein Residuue
80
0.6
6014
0.4
0.2
150
151
15L
0.0
Ugand Mom Index
Ugand Atom Index
61
1.0
0.8
64
64
6CWI
0.6
70
70
0.4
70
70
0.2
0.6
Ugand
Ugaad
Figure 10. Protein-ligand contact probability maps. The heatmaps compare the protein-ligane contact probability maps between MD reference (left) and BioMD (right), calculated using a 5 Á distance cutoff. Test cases (SGGK, IRGL, 6Q14, 6CWI) are from MISATO set.
22

---

Published as a conference paper at ICLR 2026
DereiShgLD
SNIND
NeuraMD_5DE
NeuraMD_G0E
BiuMD obs
BloMD-rel
60GE
BWXG
6M88
6Q36
RMSF (A)
25
576b
5NBV
3TD4
IGHB
RAKSF (A)
05
0.0
Figure 11. Ligand RMSF the MISATO dataset. Line plot showing Ligand RMSF for eight different protcin-ligand systems from the MISATO test set.
BloMD
Metadynamics
Cluster
Cluster
Cluster
Cluster B
Correct rete 2 0.92
Endpoint distance 1.18
Figure 12. OOD case study. The test case (3PCU) comes from the DD-13M test set, with an unbinding path RMSD of 0.77 A (Metadynamics=0.65 Å). Its largest sequence similarity with the training is 0.26.
MO 100me
MB 1064
MD Ins
RoMD
100ms
10ms
POR: 019
40
A00 GAAN
1000 CYPN
POR_
Figure 13. TICA plot the Atlas test set. Time-lagged Independent Component Analysis (TICA) plot showing 2D conformational space MD (100 ns, 10 ns, and ns) BioMD. we use dihedral angles of the backbone (g. v features for TICA. Test cases from the Atlas test set.
23

---

Published as a conference paper at ICLR 2026
Algorithm 1: Main Inference Loop
Input: } (xo,)), Neyele = 4, C₂ = 384, C2 = 128
{ inputss < ImputFeatureEmbeddder
s|nit LinearNoBBass,, (pputs)
init & (inearNoBBas(s (mputs) LinearNoBias(s, juputs).
by
{Zig), {s,}
< 0,0;
foreach e {1, [1,..... ] do
Zij < (i) + Linear\BBiasLLayeeNorrrii./
{z.z), 1, {s:} + GraphTransffrmmrr (Xo,) {s;}. (Zih. (s, (s]mputr)))
smmit + Limea\oBBis(LLaeeNNori]]))
(raj_list((
10 foreach [l..... do
[xt,] } - ((xxae) { (s) inputts J, (z,)))
11
traj lidiiaaddddd)
12
13
[Xhis] = traj__liaa list;
14
return traj_list
Algorithm TrainFlow
Input: (R J. {f*).
(s) imputs) (strunk), ). (2!runk)
# Indepentent noise levels:
(2/(0, (0, 1). ... 1)),
x } N8
Xi < CenteeRandmAugmmeetttooon/ (xx
[X]] [xx] (1 {x))
(ū7) + (lowMModulc((77) J: {Khis } [f {s; imputs J, (strunk), j. (z)junk) });
flow MSE((ū)) ), - });
return Liow
Algorithm 3: SampleFlow
Input: {Krish J. [f*). { imputs (strunk), {zly {ztrunkk
x] ~ IO,
foreach in (0,0.1, 0.2...., 0.9) do
{ūf Flow Flow ( T, {P*}, }, (simputs), J, (strunk), J, (2)zunk)
xi + dt ;f;;
return
{x]}}
Usage Note: The returned bias is added to the attention logits inside the main transformer block before the softmax activation: logits +- (Q KT)/ // + bias.
24

---

Published as a conference paper at ICLR 2026
Algorithm 4: Floule
Input: {xt (xuoiay I i, [f*) {s+ [mputs] j. {struak] { trunk]
Odata 16, Catom 128, Catompair 16, Croken 768
{si). [z,z} } + FlowConditionin((t, {f*), {s!mputa, }, (strunk), { data));
Sequence-local Atom Attention with history info and aggregation to coarse coarse-grained grained tokens fai), ( (p% (ta <
AlomAAtentioonnnoodrr(f (Xxce) {f*]. {x/ (ssss (2() Catom Catompair: Ctoken);
# Full self-attention token level.;
a; < BinearNoBBas(LayeNNooo a;));
{ax} } Flow TrajectoryTransfomme Transformer(faa)) (s:), (2i),, Rig 0, block 24, Nhead 16);
a; & LayerNorm(at.
Broadcast token activations atoms and Atom Attention.:
[ < AlommAttenioonDeoodr((44444 {4x Pk
skip skip skip });
10 return {ūr}
Algorithm 5: FlowConditioning
Input: i, {f*), {8) [s?"puts,, ), {strunk), (=]junk], Odatas Cz 128, C, 384
# Pair conditioning;
2ij < (aarrNBBBa(LyyyNrrmmaaa
foreach b E 2] do
L Fij += Transition(zy,
2);
Single conditioning:
8t concat([ strunk 8,
juputs });
8; é (anaarro foreach b E {1, 2} do
8: += Transition n = 2);
10 return (sth {
Algorithm 6: FlowTrajector hww.T/iijcccoryTTasffoomme Transformer
Input: {az}. {st]. { J. [B:g]. Nblock 24, Nnead 16
for n E [1. [1..... Nblock do
[b.] + AttentionPairrBias (fa,), {s;}, (.s,, (zig (B Nhead);
(b.) + Attertioottennnonn {a; + b:));
< by InndttionndTTannstiionBBeee
return {a.)
25

---

Published as a conference paper at ICLR 2026
Algorithm 7: TemporalAttention Module
Input: Single representation c, of shape (I' c.), where time, N residues. Output: Updated single representation output of shape (T, N, c.).
Permute dimensions to make time sequence axis for attention PermuteeC,,, dims (1,0,2))
d₂ 4
Project to Query, Key, Value for each residue independently;
Q Linearq (c,);
K f- Linear x (c)):
V Lineary (c));
# Calculate scaled dot-product attention scores across time;
dg dimension of K;
logits + (Q. //Q// Vdz:
weights < Softmax (logits, dim -));
10 Apply attention weights to values;
output' + weights V;
12
Permute back to the original dimension order;
output +- Permute Permute(output dims (1, 0)));
14 return outpul;
Algorithm
Module
Input: Pair representation of shape (N, N, c2), where N is residues. Output: Attention bias bias of shape (NNead, N).
# Project c2 to match the number of attention heads;
bias' <- Linearbias (c;);
# Permute dimensions to align with attention logits shape;
bias + Permute (aimst dims = (2, 1));
return bias;
26


Reviews:
[{"review_id": "7VoMgP1fV5", "reviewer": "Reviewer_TABm", "summary": "The paper introduces BioMD, an all-atom generative model for biomolecular dynamics simulation. It combines flow matching with a hierarchical forecasting–interpolation framework to generate long-timescale protein–ligand trajectories efficiently. The model is evaluated on MISATO and DD-13M, achieving high physical plausibility and impressive ligand-unbinding success rates.", "questions": "See above", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "* The work addresses a challenging and underexplored problem: generating full, time-resolved protein–ligand trajectories at all-atom resolution. The motivation is clearly articulated and relevant to molecular simulation and drug discovery.\n\n* The proposed hierarchical framework is conceptually elegant. Separating long-range forecasting from local interpolation aligns well with physical intuition and effectively mitigates error accumulation over long trajectories.\n\n* The unified conditional flow-matching formulation (“noising-as-masking”) is technically interesting and flexible, enabling both forecasting and interpolation within one architecture.\n\n* The experiments are comprehensive within the chosen scope.\n\n* The auxiliary geometric losses (bond, collision, and center penalties) are physically meaningful and help ensure chemically plausible outputs.", "weaknesses": "* The paper mainly reports geometric errors and RMSF correlations, but omits more direct physical checks (e.g., energy or temperature conservation, relaxation consistency). Without such validation, “physical plausibility” remains largely geometric.\n\n* The experimental baselines are narrow. Although the authors mention Str2Str (Lu et al., 2024) and BioEmu (Lewis et al., 2025) in the related work, these strong generative baselines are not included in the experiments. The justification that they are “time-agnostic” is somewhat artificial, since BioMD also measures static geometric stability and RMSF correlations—metrics directly comparable to those ensemble models. Moreover, both Str2Str and BioEmu are flow- or diffusion-based architectures closely related to BioMD’s conditional flow-matching design. Their absence weakens the empirical comparison and makes it unclear whether BioMD truly surpasses existing conformational generators.\n\n*  The architecture figure is dense but lacks quantitative configuration (model size, training schedule, runtime benchmarking). \n\n* It is unclear whether the hierarchical masking approach can generalize to unseen complexes or trajectories longer than those used in training. A discussion of scalability and failure cases would strengthen the work.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "HvqZcALRnM", "reviewer": "Reviewer_1cTN", "summary": "This paper introduces BioMD, a novel all-atom generative model designed to simulate long-timescale protein-ligand systems. The primary motivation is to overcome the significant computational cost of traditional Molecular Dynamics simulations, which severely limits the accessible timescales for many biologically relevant processes. The core contribution is a hierarchical framework that decomposes the complex task of long-trajectory generation into coarse-grained forecasting and fine-grained interpolation. The authors validate BioMD on two datasets. BioMD generates trajectories with high physical plausibility and successfully generates complete ligand unbinding paths.", "questions": "- Was any data cleaning performed on the MISATO dataset? For instance, were unstable simulation results (according to the RMSD result) filtered out before training?\n\n - Did you evaluate the difference between the generated protein conformations and the reference MD simulations, particularly regarding key conformational changes in the binding pocket? \n\n- Have you attempted to extend the model to even longer timescales (beyond those in DD-13M)? Given the acknowledged error accumulation, do you observe simulation collapse or instability over extended periods? Furthermore, the MISATO simulations are relatively short (8 ns), making it difficult to observe significant events like major conformational transitions. Could you provide representative case studies to demonstrate the model's ability (or limitations) in capturing such rare but significant events?\n\n- Regarding Figure 6, how consistent are the generated unbinding pathways with the reference metadynamics results?\n\n- I recommend including related works that do not directly recover the dynamics trajectories but recover the key binding pose and also use flow matching or diffusion models. There are some related works recommended. For example, DynamicBind [1] learns to transform unbound to bound states of protein-ligand complexes; DynamicFlow [2] uses a flow model to transform protein apo state and a noisy ligand to a holo state and a binding ligand to consider protein dynamics in structure-based drug design.\n\n\n**Note: I might raise my score according to the situation of rebuttal discussion and paper revision.**\n\nReferences:\n\n[1] DynamicBind: predicting ligand-specific protein-ligand complex structure with a deep equivariant generative model, Nature Machine Intelligence 2024\n\n[2] Integrating Protein Dynamics into Structure-Based Drug Design via Full-Atom Stochastic Flows. ICLR 2025.", "limitations": "", "rating": 4, "confidence": 2, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "1. The hierarchical  framework decomposes the long-sequence problem into \"capturing the backbone first, then filling the details,\"  effectively reduces the sequence length and the complexity of long-range dependencies in the generation process.\n2. Using a single flow matching model to achieve two different tasks (forecasting and interpolation) simply by changing the masking schedules.\n3. Supportive experimental results, especially the ligand unbinding experiment.\n4. The paper is clearly written and presented.", "weaknesses": "1. **Definition of \"Long-Timescale\"**: The paper claims to simulate \"long-timescale\" dynamics. However, the MISATO dataset's simulation time (8 ns) is relatively short, insufficient to observe major protein conformational transitions (e.g., open-to-closed states). While DD-13M is a dissociation task, it was generated by metadynamics, an enhanced sampling method, not a spontaneous process in equilibrium. Therefore, the model's ability to simulate equilibrium rare events at the microsecond (µs) to millisecond (ms) level has not been demonstrated\n2. **Insufficient Evaluation of Protein Conformational Changes**: The evaluation of protein flexibility relies mainly on RMSF, which is an aggregate metric. In drug discovery, we are often more concerned with specific conformational changes near the binding pocket, such as the flipping of key side-chains or the movement of \"gating\" residues. The paper lacks an in-depth case study on these critical biological events.\n3. The model's performance on DD-13M is highly dependent on the quality and coverage of the training data. The DD-13M dataset was generated using metadynamics, which means BioMD may have learned a biased distribution of \"how to dissociate quickly\" rather than the dynamics of the true Boltzmann distribution. Could the author provide some study on the OOD cases using BioMD?\n4. Lack of ablation studies, such as how the definition of \"large-step\" influences the final results and what is the effect of different sampling steps.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "otgedzudWK", "reviewer": "Reviewer_bQYZ", "summary": "The paper proposes BioMD, a model for all-atom molecular dynamics prediction focusing on protein–ligand complexes. BioMD is based on flow matching conditioned on the protein sequence, ligand atom types, and the initial complex structure. By combining coarse-grained prediction for long time intervals with fine-grained interpolation for the short intervals in between, BioMD can model long molecular dynamics sequences. In the model architecture, the decoupling of intra-frame spatial attention and inter-frame temporal dependencies enables the simultaneous prediction of velocities across multiple time frames at all-atom resolution. To enhance physical plausibility, auxiliary losses are introduced. Experimental results demonstrate superior physical stability, conformational flexibility, and ligand unbinding success rates with BioMD.", "questions": "1. In the formulation of the velocity model in Equation (2) (line 195), $\\mathbf{Z}$ is said to contain the “amino acid sequence” and “ligand atom types.” Since the model is all-atom, how and when is atom-wise information from the residues incorporated into the model? From Figure 3, it seems that the atom-wise information of residues is only used in the velocity module after the $SE(3)$ Transformer backbone.\n2. The paper focuses on protein–ligand complexes. However, the proposed techniques could, in principle, be applied to protein monomers and multimers as well. What is the motivation for focusing solely on protein–ligand complexes in this work?\n3. How is the time interval in the coarse-grained forecasting stage chosen? How does the scale of this interval influence the quality and efficiency of trajectory generation?", "limitations": "", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper addresses an important problem: long-timescale molecular modeling of protein–ligand binding complexes at all-atom resolution, which has significant applications in drug discovery.\n2. The efficient training and inference pipeline, along with the model architecture design, enables the prediction of multiple time frames simultaneously—crucial for modeling long MD trajectories. Experimental results show that the forecasting–interpolation workflow effectively reduces error accumulation.\n3. The evaluation of the unbinding success rate demonstrates the model’s capability to capture real biomolecular processes and the authenticity of the predicted MD trajectories.\n4. The paper is clearly written and well-presented.", "weaknesses": "1. Some implementation details are missing from both the main text and the appendix. For example, the architecture of the FlowTrajectoryTransformer is not fully described. In Algorithm 6, TemporalAttention and AttentionPairBias are used, but their definitions are missing, and only limited information about these modules is provided in the main text.\n2. The evaluation focuses primarily on low-dimensional metrics, including bond/angle geometry and RMSF. The unbinding process is described using the centroid coordinates of the ligand. These metrics can reflect trajectory quality to some extent, but higher-dimensional evaluations or demonstrations would be more informative. For instance, analyzing distributional similarity between generated and ground-truth trajectories projected onto two TICA components, or examining the transition processes between clusters of conformers in the projected space, would provide deeper insight.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "YhaPsjOpxF", "reviewer": "Reviewer_gM8W", "summary": "This paper develops a method BioMD, which generates time evolution of protein and ligand structures using a hierarchical architecture and flow matching. BioMD combined coarse grained time interval for structure ensemble generation conditioned on an initial frame, and fine grained interpolation conditioned on the starting and ending frames. It has demonstrated performance on MD datasets MISATO (unbiased protein-ligand simulations) and DD-13M (enhanced sampling MD simulations for ligand unbinding events), and shown advantage in the physicality, flexibility and unbinding success.", "questions": "- Section A.3.1 states: \"we calculate the deviations of bond lengths and bond angles with respect to the initial frame of the reference trajectories.\" However, in Section 5, the \"Static\" baseline is defined as \"the initial conformation of the system is held constant throughout the entire trajectory.\" If errors are computed relative to the initial frame, the Static model should have exactly zero error, yet Table 1 and Table 2 report non-zero errors for Static. The authors need to clarify: (1) what exactly is the reference for computing these errors, and (2) whether errors are computed for ligand only, protein only, or both.\n- The phrase \"sampling every k = 10 steps\" (Section 4.2.1) is imprecise without specifying the timestep used in the underlying simulations. The actual time interval (in ps or fs) should be stated explicitly.", "limitations": "", "rating": 4, "confidence": 5, "soundness": 3, "presentation": 4, "contribution": 2, "strengths": "- The model architecture is novel, proposing a hierarchical framework in combination with flow matching\n- The authors have tested multiple different ways of generation, including predicting relative conformational change v.s. absolute conformations, and predicting all frames v.s. auto-regressively, providing a thorough examination\n- The BioMD method has shown advantages in the MISATO and DD-13M benchmarks, in the physicality, flexibility and unbinding success.", "weaknesses": "- In the introduction, the paper states: \"The core bottleneck lies in the intensive calculation of non-bonded forces, particularly van der Waals and electrostatic interactions, which scale quadratically with the number of atoms.\" This is misleading, because with methods like Particle Mesh Ewald (PME) the scaling is approximately O(N log N), not O(N²). \n- The model is trained with a fixed time interval between frames (k=10 steps in coarse-graining). Can the model generalize to different time intervals, or is it restricted to this specific temporal resolution? If the latter, this should be explicitly stated as a limitation. \n- The DD-13M dataset is generated using enhanced sampling simulation, instead of unbiased MD. Therefore, BioMD is learning to reproduce biased trajectories, not true thermodynamic/kinetic unbinding behavior. The generated trajectories may not reflect realistic dissociation kinetics, and the model may be learning to mimic the \"push\" behavior by the biased potential in enhanced sampling. Currently, the paper does not clearly distinguish this subtlety.\n- The success criterion (Section A.3.4) only requires that \"at least one predicted ligand centroid position lies outside this convex hull.\" The success rate is thus measuring \"can find a way out\" rather than \"finds the correct way out.\" \n- The reported MAE and MSE metrics average over all frames and trajectories. This can mask catastrophic failures in individual frames. The authors should consider providing more comprehensive failure analysis, including per-frame error distributions, percentage of frames with physical violations etc.\n- The paper proposes the hierarchical architecture, but lacks the justification and ablation on the necessity of the hierarchy. Specifically, why the fine-grained generation is needed and how it will influence the performance should be demonstrated and discussed.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Poster)