A multi-grained symmetric differential equation model for learning protein-ligand binding dynamics

Domain: physics
Published: 2025-12-30
Editorial outcome: major_revision

Abstract
Molecular dynamics (MD) simulation is a key tool in drug discovery for predicting protein-ligand binding affinities, transport properties, and pocket dynamics. While advances in numerical and machine learning (ML) methods have improved MD efficiency, accurately modeling long-timescale dynamics remains challenging. We introduce NeuralMD, an ML surrogate that accelerates and enhances MD simulations of protein-ligand binding. NeuralMD employs a physics-informed, multi-grained, group-symmetric framework comprising (1) BindingNet, which enforces symmetry via vector frames and captures multi-level protein-ligand interactions, and (2) an augmented neural differential equation solver that learns trajectories under Newtonian mechanics. Across ten single-trajectory and three multi-trajectory tasks, NeuralMD achieves up to 15× lower reconstruction error and 70% higher validity than existing ML baselines. The predicted oscillations closely align with ground-truth dynamics, establishing NeuralMD as a foundation for next-generation protein-ligand simulation research.

Introduction
The molecular dynamics (MD) simulation for protein-ligand binding is one of the fundamental tasks in drug discovery 1–4. Such simulations of binding systems are a key component of the drug discovery pipeline to select, refine, and tailor the chemical structures of potential drugs to enhance their efficacy and specificity. To simulate the protein-ligand dynamics, numerical MD methods have been extensively developed 5,6. However, the numerical MD methods are computationally expensive due to the expensive force calculations on individual atoms in a large protein-ligand system.

To alleviate this issue, machine learning (ML) surrogates have been proposed to either augment or replace numerical MD methods to estimate the MD trajectories. However, all prior ML approaches for MD simulation are limited to single-system and not protein-ligand complex 7–10. A primary reason is the lack of large-scale datasets for protein-ligand binding (the first large-scale dataset with binding dynamics was released in May 2023 11). Further, prior ML-based MD approaches limit to fitting the energy surfaces so as to study the MD dynamics on a small timestep (e.g., 1e-15s) 12–19, while simulation on a larger timestep (e.g., 1e-9s) is needed for specific tasks, such as detecting the transient and cryptic states in binding dynamics 20. However, current longer-time ML MD simulations are challenging due to the catastrophic buildup of errors over longer rollouts 21.

Another critical aspect that needs to be considered in ML-based modeling is the group symmetry present in the protein-ligand geometry. Specifically, the geometric function over molecular systems should be equivariant to rotation and translation, i.e., SE(3)-equivariance. One principled approach to satisfy equivariance is to use vector frames, which have been previously explored for single molecules 22, but not yet for the binding complexes. The vector frame basis achieves SE(3)-equivariance by projecting vectors (e.g., positions and accelerations) to the vector frame basis, and such a projection can maintain the equivariant property with efficient calculations 17.

To this end, we propose NeuralMD, a multi-grained physics-informed approach designed to handle extended-timescale MD simulations in protein-ligand binding. Our multi-grained method explicitly decomposes the large molecular complexes into three granularities to obtain an efficient approach for modeling the large molecular system: the atoms in ligands, the backbone structures in proteins, and the residue-atom pairs in binding complexes. We achieve group symmetry in BindingNet through the incorporation of vector frames, and include three levels of vector frame bases for multi-grained modeling, from the atom and backbone level to the residue level for binding interactions.

Additionally, our ML approach, NeuralMD, leverages data-driven techniques to learn Newtonian mechanics. In MD, the movement of atoms is determined by Newton’s second law, F = m ⋅ a, where F is the force, m is the mass, and a is the acceleration of each atom. By integrating acceleration and velocity with respect to time, we can obtain the velocities and positions, respectively, as shown in Fig. 1. Thus, in NeuralMD, we formulate the trajectory simulation as a second-order ordinary differential equation (ODE) or second-order stochastic differential equation (SDE) problem. We augment derivative space by concurrently calculating the accelerations and velocities, allowing simultaneous integration of velocities and positions.

To evaluate the effectiveness and efficiency of NeuralMD, we design ten single-trajectory and three multi-trajectory binding simulation tasks. For quantitative assessment, we employed two reconstruction metrics and two validity metrics. The results demonstrate that NeuralMD consistently outperforms other ML methods, achieving up to 15 × reduction in reconstruction error and 70% increase in validity. Additionally, for qualitative analysis, we introduced a fluctuation measure to characterize ligand oscillation, supported by the Supplementary Video 1. These results highlight that while NeuralMD is not perfect, it generalizes more effectively than other ML baselines. Finally, in terms of efficiency, NeuralMD delivers up to over 1K × speedup compared to numerical methods.

Results
Preliminaries
Ligand data structure. In this work, we consider binding complexes involving small molecules as ligands. Small molecules can be treated as sets of atoms in the 3D Euclidean space, {f(l), x(l)}, where f(l) and x(l) represent the atomic numbers and 3D Euclidean coordinates for atoms in each ligand, respectively.

Protein data structure. Proteins are macromolecules, which are essentially chains of amino acids or residues. There are 20 natural amino acids, and each amino acid is a small molecule. Noticeably, amino acids are made up of three components: a basic amino group (-NH2), an acidic carboxyl group (-COOH), and an organic R group (or side chain) that is distinct to each amino acid. Additionally, the carbon that connects all three groups is called Cα. Due to the large number of atoms in proteins, this work proposes a multi-grained method for modeling protein-ligand complexes. In this regard, the backbone-level data structure for each protein is {f(p), x(p), x(p), x(p)}, for the residue type and the coordinates of N, Cα, and C in the residue, respectively. We may omit the superscript in the coordinates of backbone atoms, as these backbone structures are different from protein residues. In addition to the backbone level, as a coarser-grained view, we further consider residue-level information for modeling binding interactions, {f(p), x(p)}, where the coordinate of Cα is taken as the residue-level coordinate, i.e., x(p) ≜ x(p)Cα.

MD simulations. Generally, MD describes how each atom in a molecular system moves over time, following Newton’s second law of motion:

F = m ⋅ a = m ⋅ d^2x / dt^2, (1)

where F is the force, m is the mass, a is the acceleration, x is the position, and t is the time. Then, an MD simulation will take a second-order integration to get the trajectories for molecular systems like a small molecule, a protein, a polymer, or a protein-ligand complex. The numerical MD methods can be classified into classical MD and ab-initio MD, where the difference lies in how the force on each atom is calculated: classical MD uses force field approaches to predict the atomic forces 5, while ab-initio MD calculates the forces using quantum mechanical methods, such as density functional theory (DFT) 6. More recently, ML MD methods have opened another perspective by utilizing the group symmetric tools for geometric representation and energy prediction 12–19, as well as the automatic differential tools for trajectory learning 23–28. Please check Supplementary A and Supplementary Table S1 for a more detailed discussion.

Newtonian dynamics and langevin dynamics. Numerical MD methods can be additionally divided into two categories: using Newtonian dynamics or Langevin dynamics. Newtonian dynamics is suitable for idealized systems with negligible thermal effects or when deterministic trajectories are required, while Langevin dynamics is adopted where thermal effects play a significant role and when the system is being studied at a finite temperature. In the ML-based MD simulations, adopting Newtonian dynamics or Langevin dynamics can be treated as an option for introducing different inductive biases. In this work, we propose two versions: an ODE solver and a SDE solver concerning Newtonian dynamics and Langevin dynamics, respectively. Noticeably, our experiential dataset, MISATO 11, uses Newtonian dynamics with Langevin thermostats, and the information on solvent molecules is not provided.

Problem setting: MD simulation in protein-ligand binding. In this work, we are interested in learning the MD simulation in the protein-ligand binding system, and we consider the semi-flexible setting 29, i.e., proteins with rigid structures and ligands with flexible movements. Thus, the problem is formulated as follows: suppose we have a fixed protein structure {f(p), x(p), x(p), x(p)} and a ligand with its initial structure and velocity, {f(l), x(l), v(l)}. We want to predict the trajectories of ligands following the Newtonian dynamics, i.e., the movement of x(l) over time. We also want to clarify two critical points about this problem setting: (1) Our task is trajectory prediction, i.e., positions as labels, and no explicit energy and force are considered as labels. ML methods for energy prediction followed by a numerical ODE/SDE solver require a smaller timestep (e.g., 1e-15s), while trajectory prediction, which directly predicts the coordinates over time, is agnostic to the magnitude of the timestep. This is appealing for tasks with larger timestep (e.g., 1e-9s), as will be discussed below. (2) Each trajectory is composed of a series of geometries of molecules, and such geometries are called snapshots. We avoid using frame since we will introduce the vector frame in modeling the binding complex in the following sections.

Pipeline: neuralMD
In this section, we briefly introduce NeuralMD, our framework for learning MD simulations in protein-ligand binding. It has two main phases: (1) BindingNet: a multi-grained SE(3)-equivariant geometric model that represents the protein-ligand complex across three granularities. It utilizes vector frame bases at the atom level for ligands, the backbone level for proteins, and the residue level for protein-ligand complexes. (2) Dynamics Solver: either a second-order ODE solver or a second-order SDE, used to model trajectories as Newtonian or Langevin dynamics, respectively. Figure 2 depicts the key components and pipeline of NeuralMD. Please check the Methods Section for more details.

(1) Multi-grained vector frames. We consider three data perspectives: atom level for ligands, backbone level for proteins, and residue level for protein-ligand complexes. These perspectives allow us to construct vector frames corresponding to each level, serving as reference bases for SE(3)-geometric modeling.

(2) Multi-grained SE(3)-equivariant binding modeling: bindingNet. To model protein-ligand binding prediction, we introduce BindingNet, which leverages the three vector frames established in the first step. The core idea is to achieve rotation-equivariant modeling by projecting vector variables onto the reference frames.

(3) MD modeling of binding complexes: NeuralMD. Finally, we leverage the binding energy predictions from BindingNet to simulate protein-ligand binding dynamics. This involves using a second-order ODE for Newtonian dynamics and a SDE for Langevin dynamics. These three steps collectively form the foundation of NeuralMD for modeling binding dynamics.

Experiment Setting
Datasets. One of the main bottlenecks of studying ML for MD simulation in protein-ligand binding is insufficient data. Recently, the community has put more effort into gathering the datasets, and we consider MISATO in our work 11. It is built on 16,972 experimental protein-ligand complexes extracted from the protein data bank (PDB) 30. Such data is obtained using X-ray crystallography, Nuclear Magnetic Resonance (NMR), or Cryo-Electron Microscopy (Cryo-EM), where systematic errors are unavoidable. This motivates the MISATO project, which utilizes semi-empirical quantum mechanics for structural curation and refinement, including regularization of the ligand geometry. For each protein-ligand complex, the trajectory comprises 100 snapshots in 8 nanoseconds under fixed temperature and pressure. In Supplementary D (Supplementary Figs. S1 and S2), we list the basic statistics of MISATO, e.g., the number of atoms in small molecule ligands, and the number of residues in proteins.

Baselines. Using ML for energy and force prediction, followed by trajectory prediction using numerical integration method, has been widely explored in the community, e.g., HDNNPs 7, DeePMD 9, TorchMD 31, and Allegro-LAMMPS 32. Here we extend this paradigm for binding dynamics and propose VerletMD, which utilizes BindingNet for energy prediction on each snapshot and velocity Verlet algorithm to get the trajectory. Additionally, we mainly focus on ML methods for trajectory prediction in this work, i.e., no energy or force is considered as labels. GNN-MD is to apply geometric graph neural networks (GNNs) to predict the trajectories in an auto-regressive manner 11,33. More concretely, GNN-MD takes the coordinates and other molecular information as inputs at time t and predict the coordinates at time t+1. Denoising LD (denoising diffusion for Langevin dynamics) 33–35 is a baseline method that models the trajectory prediction as denoising diffusion task 36, and the inference for binding trajectory essentially becomes the Langevin dynamics. CG-MD learns a dynamic GNN and a score GNN 33, which are essentially the hybrid of GNN-MD and Denoising LD. Here, to make the comparison more explicit, we compare these two methods (GNN-MD and Denoising LD) separately. We want to highlight that we keep the same backbone model, BindingNet, for force or position prediction for all the baselines and NeuralMD.

Experiments Settings. We consider two experimental settings. The first type of experiment is the single-trajectory prediction, where both the training and test data are snapshots from the same trajectory, and they are divided temporally. The second type of experiment is the multi-trajectory prediction, where each data point is the sequence of all the snapshots from one trajectory, and the training and test data correspond to different sets of trajectories. All the experiments are conducted with three random seeds (0, 42, 123), and details are provided in Supplementary E and Supplementary Table S2. The mean results are reported in the main article, and the standard deviations are reported in Supplementary F and Supplementary Tables S3-S6.

Evaluation Metrics. For quantitative evaluation, we use two reconstruction metrics and two validity metrics: (1) Reconstruction. For MD simulation, the evaluation is a critical factor for evaluating trajectory prediction. For both experiment settings, the trajectory reconstruction is the most straightforward evaluation metric. To evaluate this, we take both the mean absolute error (MAE) and root mean squared error (RMSE) between the predicted coordinates and ground-truth coordinates over all snapshots in the test set. (2) Validity. Molecules possess certain inherent properties, and we would like to measure the validity of the sampled positions. The first validity metric is the Matching metric. It is defined as the mean squared error between the atom pairwise distance from ground-truth conformations and our sampled conformations. This Matching metric reveals whether the model correctly reflects the true distance relationships in physical space. The second validity metric is Stability. It is defined as P i,j (|true distance - pred distance| ≤ Δ), where we take Δ = 0.5 Å. The underlying intuition is that, over long-time simulations, the predicted trajectory could enter an unrealistic state, such as bond breaking. Stability quantifies how often the predicted distances between atoms stay within a reasonable range of the true distances, indicating that the trajectory does not exhibit pathological behavior.

Generalization Among Multiple Trajectories
The first task is to test the generalization ability of NeuralMD among different trajectories. The MISATO dataset includes 13,765 protein-ligand complexes 11, and we first create two small datasets by randomly sampling 100 and 1k complexes, respectively. Then, we take 80%-10%-10% for training, validation, and testing on both datasets. We also consider the whole MISATO dataset, where the data split has already been provided. After removing the complexes with peptide ligands, we have 13,066, 1357, and 1357 complexes for training, validation, and testing, respectively.

The quantitative results are in Table 1. The first observation is that VerletMD has worse performance on all three datasets, and the performance gap with other methods is even larger compared to the single-trajectory prediction, as will be introduced below. When comparing VerletMD, NeuralMD reaches up to 15 × reduction in reconstruction error and 70% increase in validity. The second observation is that the other two baselines, GNN-MD and Denoising LD, show similar performance, while NeuralMD outperforms in all datasets and all four metrics. The last observation is that the two validity metrics (Matching and Stability) are more distinguishable than the two trajectory reconstruction metrics (MAE and RMSE), suggesting that ligand oscillation may contain more informative details. We will explore this further in a qualitative study below.

Generalization of one single trajectory
The second paradigm is the generalization of a single trajectory. The datasets are from the same protein-ligand trajectory, with the first 80 snapshots for training and the last 20 for testing.

The quantitative results are in Table 2. The first observation is that the baseline VertletMD has a clear performance gap compared to the other methods. There are two possible reasons: (1) Using ML models to predict the energy (or force) at each snapshot, and then using a numerical integration algorithm can fail in the long-time simulations; (2) ML for energy prediction methods require more data to train than the ML for coordinate prediction methods, thus they can perform worse in the low-data regime. The second observation is that, in general, both variants of NeuralMD demonstrate competitive performance across all ten tasks and four metrics, except for the reconstruction metrics on 3B9S and the validity metrics on 4K6W. The overall generalization performance reveals the potential of NeuralMD compared to existing ML baselines. Specifically, Stability (%) serves as a distinctive factor in method comparisons, with the two variants of NeuralMD outperforming on nine tasks by up to 70%.

Qualitative analysis on oscillation
In addition to the quantitative measures like reconstruction and validity, we would also like to discuss the qualitative performance of ML models for MD simulation. We begin by analyzing the Root Mean Square Fluctuation of the ligand (RMSF-Ligand) along single trajectories. RMSF-Ligand quantifies the mean square deviation of each ligand atom’s position from its average position over the trajectory, providing a measure of positional fluctuations. Formally, it is defined as:

RMSF-Ligand = sqrt(1/N * Σ_{i=1}^N ||r_{i,t} - <r_i>||^2), (2)

where r_{i,t} is the i-th atom position at time t (ground truth or sampled by ML methods) and <r_i> is the average position for the i-th atom across time. To capture localized oscillations more effectively, we employ a sliding window approach, calculating RMSF-Ligand over every five snapshots, such as snapshots [81, 85] and [96, 100].

RMSF-Ligand is a metric that measures the flexibility of individual atoms or residues during a simulation, reflecting how much they deviate from their average positions over time. It provides valuable insights into the extent of fluctuations in the system. Specifically, RMSF-Ligand describes the oscillations of individual trajectories. The objective here is to plot the RMSF-Ligand curves for the ground truth and four ML methods, enabling a comparison to determine which ML method aligns most closely with the ground truth curve. Notice that VerletMD is not considered as it cannot converge on the single-trajectory experiments (according to Table 2). The RMSF-Ligand curves are illustrated in Fig. 3. First, we observe that, overall, the two RMSF-Ligand curves produced by NeuralMD align better with the ground truth compared to those from ML baselines. Notably, NeuralMD demonstrates a significant improvement on datasets PDB 4ZX0, 3EOV, 4K6W, 4G3E, and 3B9S, and shows slightly better performance on 5WIJ, 1XP6, and 6B7F. Second, we note that almost all ML methods perform poorly on 1KTI and 4YUR. For these two trajectories, the ground truth exhibits sudden positional changes, representing out-of-distribution movements that remain an open challenge for the current ML community.

In addition to RMSF-Ligand, a more intuitive approach is to visualize the trajectories directly. We provide a 20-snapshot video comparing the trajectories generated by ML methods and the ground truth in the Supplementary Video 1. Additionally, we briefly present the sampled trajectories of four complexes using three ML methods alongside the ground truth in Fig. 3b. From these visualizations, it is evident that GNN-MD occasionally collapses, while Denoising LD maintains a comparatively structured trajectory. Notably, NeuralMD demonstrates the highest stability across all cases.

Efficiency Analysis
Efficiency is another important metric, and we consider measuring it with frames per second (FPS) on a single Nvidia-V100 GPU card. One main benefit of using NeuralMD for binding simulation is its efficiency. To show this, we list the computational time in Table 3. Since all ML methods share a similar backbone architecture (BindingNet and its variants), their efficiency differences arise from their computational approaches and hyperparameters. GNN-MD is the most efficient, as it requires no integration steps. NeuralMD incorporates an augmented integration step, with the optimal step size set to 0.5 of the snapshot interval, i.e., two integration steps are performed to predict the next snapshot, making it slightly slower than VerletMD. In contrast, Denoising LD is the slowest method, as it requires smaller integration timesteps to generate the next snapshot during inference.

We further approximate the wall time using the numerical method for MD simulation (PDB 5WIJ). Concretely, we can get an estimated speed of 1 nanosecond of dynamics every 0.28 h. This is running the simulation with GROMACS 37 on 1 GPU with 16 CPU cores and a moderately sized water box at the all-atom level (around 64,000 atoms) with the step size of 2 femtoseconds. This shows that NeuralMD is approximately 25K × faster than numerical methods under optimal conditions. However, as noted earlier, the model’s efficiency is highly sensitive to hyperparameter choices, so we consider an efficiency improvement factor of at least 1K to be a conservative estimate.

Case Study on 4G3E: biological meaning of neuralMD
In this section, we provide a detailed discussion of PDB entry 4G3E, explore the potential advantages of using NeuralMD for drug development, and examine its biological implications.

Nuclear factor kappa B (NF-κB) is a type of transcriptional factor that regulates transcription, cell survival, immune responses, and the development of the immune system 38–42. NF-κB inducing kinase (NIK), also known as MAP3K14, is a crucial signaling molecule that participates in the non-canonical NF-κB pathway, which regulates immune responses and lymphoid organ development 40,42,43. NIK is a serine/threonine kinase that activates the pathway by phosphorylating IKKα (IκB kinase α), which in turn leads to the processing of NF-κB2 (p100) into its active form, p52. This activation allows the p52-RelB complex to translocate to the nucleus and initiate gene transcription that regulates various immune functions, including B cell survival, osteoclastogenesis, and lymphoid tissue development 42,43.

NIK is a promising drug target due to its pivotal role in the non-canonical NF-κB signaling pathway 40,44–46. Aberrant NIK activation and expression have been implicated in a range of pathological conditions, including autoimmune diseases, solid cancers, hematologic malignancies, cardiovascular disease, obesity, and type 2 diabetes 46. Consequently, small-molecule inhibitors targeting NIK have demonstrated therapeutic potential by modulating immune responses and reducing inflammation, making NIK an attractive target for the treatment of various diseases 45,46.

Many drugs targeting NIK have been developed, with several structures of NIK bound to inhibitors resolved, providing valuable insights for drug design 40,47,48. However, many existing drugs still require optimization to improve efficacy and specificity. Recent efforts used metabolite identification and structure-based drug design to improve pharmacokinetics, leading to inhibitors with reduced clearance and enhanced kinase selectivity, significantly lowering predicted human doses 45. These advancements demonstrate how targeted optimization strategies can address key limitations. Our NeuralMD approach can complement such strategies by offering dynamic insights into inhibitor binding, guiding the rational design of improved protein-targeted therapies.

The example of the complex structure of the inhibitor CMP1, a 6-alkynylindoline, bound to its target protein, NIK (Mus musculus; PDB: 4G3E), highlights the final binding state of the NIK inhibitor at its target region of NIK 40. However, the dynamic process leading to this interaction remains poorly understood. Despite the development of numerous approaches for simulating these inhibitor binding processes (Fig. 3), their effectiveness is limited. For example, these methods, GNN-MD and Denoising LD, frequently produce unrealistic structural distortions of CMP1 during the simulation (Fig. 3).

In contrast, our NeuralMD model not only accelerates the simulation of the binding trajectory between CMP1 and NIK but also avoids such exaggerated deformations. NeuralMD simulations reveal that a specific structural region within the CMP1 molecule maintains stability or exhibits consistent, regular fluctuations during the binding process (Fig. 3b and Supplementary Video 1). These insights are invaluable for optimizing the NIK inhibitors, particularly by targeting molecular regions that enhance binding affinity and selectivity.

Furthermore, similar results have been observed across numerous examples in our studies (Supplementary Video 1). This demonstrates the broad applicability of NeuralMD for drug design, emphasizing its potential to streamline the development of therapeutics for various targets.

Discussion
To sum up, we devise NeuralMD, an ML framework that incorporates a multi-grained group symmetric network architecture and second-order Newtonian mechanics, enabling accurate predictions of binding dynamics on a large timescale. Not only is such a timescale critical for understanding the dynamic nature of the ligand-protein complex, but our work marks the approach in developing a framework to predict coordinates for MD simulation in protein-ligand binding. We quantitatively and qualitatively verify that NeuralMD achieves superior performance on 13 binding dynamics tasks.

One potential limitation of our work is the dataset. Currently, we are using the MISATO dataset, a protein-ligand dynamics dataset with a large timestep. However, NeuralMD is agnostic to the timestep, and it can also be applied to binding dynamics datasets with a timestep of a femtosecond. This may require the effort of the whole community for the dataset construction.

The second limitation of our work is GPU memory efficiency. Currently, we run multi- and single-trajectory MD simulations in a semi-flexible setting. While the flexible setting, where proteins can move, is more practical, it requires up to 100 times more GPU memory in the multi-trajectory case, as all protein conformations across 100 snapshots per trajectory must be modeled. For future work, we aim to develop a memory-efficient, physics-informed foundation model for proteins, capable of rapidly adapting to MD simulation tasks, including MISATO.

Another limitation is the lack of verification through practical research. Although we have demonstrated the effectiveness of NeuralMD through computational simulations, we have not yet validated its biological relevance using experimental techniques, such as X-ray crystallization and cryogenic electron microscopy. While our model has shown promising results in terms of simulating protein-ligand interactions, its applicability in drug development requires further confirmation through practical research. To address this gap, we plan to synthesize the molecules based on the simulation results provided by NeuralMD. Such studies will help to enhance the reliability of NeuralMD as a tool for guiding the design of therapeutic agents.

Methods
Multi-grained vector frames
Proteins are macromolecules composed of up to thousands of residues (amino acids), where each residue is a small molecule. Thus, it is infeasible to model all the atoms in proteins due to the large volume of the system, and such an issue also holds for the protein-ligand complex. To address this issue, we propose BindingNet, a multi-grained SE(3)-equivariant model, to capture the interactions between a ligand and a protein. The vector frame basis ensures SE(3)-equivariance, and the multi-granularity is achieved by considering frames at three levels.

Vector frame basis for SE(3)-equivariance. Recall that the geometric representation of the whole molecular system needs to follow the physical properties of the equivariance w.r.t. rotation and translation. Such a group symmetric property is called SE(3)-equivariance. We also want to point out that the representation function should be reflection-equivariant for properties like energy, yet it is not for properties like chirality or ligand modeling in rigid protein structures. The vector frame basis inherently accommodates such reflection antisymmetry, and we leave a more detailed discussion in Supplementary C, along with the proof of group symmetry of the vector frame basis. In the following, we introduce three levels of vector frames for multi-grained modeling.

Atom-level vector frame for ligands. For small molecule ligands, we first extract atom pairs (i, j) within the distance cutoff c, and the vector frame basis is constructed using the Gram-Schmidt as:

F_{ij} = [ (x_i - x_j) / ||x_i - x_j||, (x_i × x_j) / ||x_i × x_j||, (x_i - x_j) × (x_i × x_j) / ||(x_i - x_j) × (x_i × x_j)|| ] (3)

where × is the cross product. Note that both x_i and x_j are for geometries at time t - henceforth, we omit the subscript t for brevity. Such an atom-level vector frame allows us to do SE(3)-equivariant message passing to get the atom-level representation.

Backbone-level vector frame for proteins. Proteins can be treated as chains of residues, where each residue possesses a backbone structure. The backbone structure comprises an amino group, a carboxyl group, and an alpha carbon, delegated as N - Cα - C. Such a structure serves as a natural way to build the vector frame. For each residue in the protein, the coordinates are x_N, x_{Cα}, and x_C, then the backbone-level vector frame for this residue is:

F_{protein} = [ (x_N - x_{Cα}) / ||x_N - x_{Cα}||, (x_{Cα} - x_C) / ||x_{Cα} - x_C||, (x_N - x_{Cα}) × (x_{Cα} - x_C) / ||(x_N - x_{Cα}) × (x_{Cα} - x_C)|| ] (4)

This is built for each residue, enabling the message passing for a residue-level representation.

Residue-level vector frame for protein-ligand complexes. It is essential to model the protein-ligand interaction to better capture the binding dynamics. We achieve this by introducing the residue-level vector frame. More concretely, proteins are sequences of residues, marked as

{(f_i^{(p)}, x_i^{(p)}), ..., (f_i^{(p)}, x_i^{(p)}), (f_{i+1}^{(p)}, x_{i+1}^{(p)}), ...}

Here, we use a cutoff threshold c to determine the interactions between ligands and proteins, and the interactive regions on proteins are called pockets. We construct the following vector frame on residues in the pockets sequentially:

F_{complex} = [ (x_i^{(p)} - x_{i+1}^{(p)}) / ||x_i^{(p)} - x_{i+1}^{(p)}||, (x_i^{(p)} × x_{i+1}^{(p)}) / ||x_i^{(p)} × x_{i+1}^{(p)}||, (x_i^{(p)} - x_{i+1}^{(p)}) × (x_i^{(p)} × x_{i+1}^{(p)}) / ||(x_i^{(p)} - x_{i+1}^{(p)}) × (x_i^{(p)} × x_{i+1}^{(p)})|| ] (5)

Through this vector frame, the message passing enables the exchange of information between atoms from ligands and residues from the pockets. The illustration of the above three levels of vector frames can be found in Fig. 2. Once we build up such three vector frames, we then design BindingNet, as will be introduced next. The key step involving vector frames is scalarization operation 49, which transforms the equivariant variables (e.g., coordinates) to invariant variables by projecting them to the three vector bases in the vector frame.

Multi-Grained SE(3)-equivariant binding force modeling: bindingNet
In this section, we introduce BindingNet, a multi-grained SE(3)-equivariant geometric model for protein-ligand binding. The input of BindingNet is the geometry of the rigid protein and the ligand at time t, while the output is the force on each atom in the ligand.

Atom-level ligand modeling. We first generate the atom embedding using one-hot encoding and then aggregate each atom’s embedding, z_i^{(l)}, by aggregating all its neighbors’ embedding within the cutoff distance c. Then, we obtain the atom’s equivariant representation by aggregating its neighborhood’s messages as (x_i^{(l)} - x_j^{(l)}) * z_i^{(l)}. A subsequent scalarization is carried out based on the atom-level vector frame as h_i^{(l)} = (h_i^{(l)} ⊕ h_j^{(l)}) * F_{ligand}, where ⊕ is the concatenation. Finally, it is passed through several equivariant message-passing layers (MPNN) defined as:

vec_i^{(l)} = vec_i^{(l)} + Agg_j vec_i^{(l)} * MLP(h_{ij}) + (x_i^{(l)} - x_j^{(p)}) * MLP(h_{ij}) (6)

where MLP(⋅) and Agg(⋅) are the multi-layer perceptron and mean aggregation functions, respectively. vec ∈ R^3 is a vector assigned to each atom and is initialized as 0. The outputs are atom representation and vector (h^{(l)} and vec^{(l)}), and they are passed to the complex module introduced below.

Backbone-level protein modeling. For the coarse-grained modeling of proteins, we consider three backbone atoms in each residue. We first obtain the atom embedding on three atom types, and then we obtain each atom’s representation z^{(p)} by aggregating its neighbor’s representation. Then, we obtain an equivariant atom representation by aggregating the edge information, (x_i^{(p)} - x_j^{(p)}) * z^{(p)}, within cutoff distance c. Following which is the scalarization on the residue frame h^{(p)} = (h^{(p)} ⊕ h^{(p)}) * F_{protein}. Recall that we also have the residue type, and with a type embedding z~^{(p)}, we can obtain the final residue-level representation using an MPNN layer as:

h^{(p)} = z~^{(p)} + (h_{N,Cα}^{(p)} + h_{Cα,C}^{(p)}) / 2 (7)

Residue-Level Complex Modeling. Once we obtain the atom-level representation (h^{(l)}, vec^{(l)}) from ligands, and backbone-level representation (h^{(p)}) from proteins, the next step is to learn the protein-ligand interaction. We first extract the residue-atom pair (i, j) with a cutoff c, based on which we obtain an equivariant interaction edge representation h_{ij} = (h_i^{(l)} + h_j^{(p)}) * (x_i^{(l)} - x_j^{(p)}). After scalarization, we can obtain invariant interaction edge representation h_{ij} = h_{ij} * F_{complex}. Finally, we adopt an equivariant MPNN layer to get the atom-level force as:

vec_i^{(pl)} = vec_i^{(l)} * MLP(h_{ij}) + (x_i^{(l)} - x_j^{(p)}) * MLP(h_{ij}) (8)

The ultimate force predictions for each atom include two parts: the internal force from the molecule vec_i^{(l)} and the external force from the protein-ligand interaction vec_{ij}^{(pl)}. To sum them up, we have:

F_i^{(l)} = vec_i^{(l)} + Agg_{j ∈ N(i)} vec_i^{(pl)} (9)

The three modules consist BindingNet. More completed descriptions can be found in Supplementary E and Supplementary Fig. S3.

BindingMD modeling: neuralMD
As previously clarified, MD follows Newtonian mechanics, and we solve it as either an ODE problem or an SDE problem. The BindingNet introduced above takes in the molecular system geometry at time t and outputs the forces. Then in this section, we describe how we use the neural differential equation solver to predict the coordinates at future snapshots. We mark this as “func” in Fig. 2. This means we can augment ODE or SDE derivative space by concurrently calculating the accelerations and velocities, allowing simultaneous integration of velocities and positions.

We want to highlight that one ML for the MD simulation research line is predicting the energy or force 9,31,32, which will be fed into the numerical integration algorithms for trajectory simulation. For accuracy, such an ML-based MD simulation must be at the femtosecond level (1e-15s). However, as shown in recent works 21,50,51, minor errors in the ML force field can lead to catastrophic failure for long-time simulations. For instance, there can be pathological behavior, such as extreme force predictions or bond breaking, within the low end of the distribution. Our experiments have yielded similar observations, as will be shown below.

In this paper, however, we overcome this issue by directly learning the extended-timescale MD trajectories (nanosecond level, 1e-9s).

Newtonian dynamics and Langevin dynamics for MD simulation. For MD simulation, if we assume the information of all particles in the molecular system (e.g., solvent molecules) is known and no thermal fluctuation is considered, then it can be modeled as Newtonian dynamics, which is essentially an ODE. For modeling, we consider using the BindingNet introduced above for force prediction at time τ, as:

F_τ^{(l)-ODE} = BindingNet(f^{(l)}, x_τ^{(l)}, f^{(p)}, x_N^{(p)}, x_{Cα}^{(p)}, x_C^{(p)}) (10)

On the other hand, Langevin dynamics introduces a stochastic component for large molecular systems with thermal fluctuations. Concretely, Langevin dynamics is an extension of the standard Newtonian dynamics with the addition of damping and random noise terms: F - γmv + sqrt(2mγk_B T)R(t), where γ is the damping constant or collision frequency, T is the temperature, k_B is the Boltzmann’s constant, and R(t) is a delta-correlated stationary Gaussian process with zero-mean. For our implementation, we learn this equation and arguments in a data-driven way with reparameterization as:

F_τ^{(l)-SDE} = BindingNet(f^{(l)}, x_τ^{(l)}, f^{(p)}, x_N^{(p)}, x_{Cα}^{(p)}, x_C^{(p)}) + BindingNet-Ligand(f^{(l)}, x_τ^{(l)}) * ϵ (11)

where ϵ ~ N(0, 1) is the standard Gaussian noise. Ultimately, adopting either F_τ^{(l)-ODE} or F_τ^{(l)-SDE} as the force prediction F_τ^{(l)} on each atom, the coordinates at time t can be obtained after integration as:

a_τ^{(l)} = F_τ^{(l)} / m, v^_t^{(l)} = v_0^{(l)} + ∫_0^t a_τ^{(l)} dτ, x^_t^{(l)} = x_0^{(l)} + ∫_0^t v^_τ^{(l)} dτ (12)

The training objective is the mean absolute error (MAE) or root mean squared error (RMSE) between the predicted coordinates and ground-truth coordinates along the whole trajectories:

L = E_t ||x^_t^{(l)} - x_t^{(l)}|| (13)

An illustration of NeuralMD pipeline is in Fig. 2.

Second-order differential equation. We also want to highlight that NeuralMD is solving a second-order differential equation, as in Newtonian and Langevin dynamics. The function, which returns the first-order derivative. Then, the outputs of the differential function will be integrated using algorithms like the Euler algorithm. To learn the MD trajectory following second-order ODE and SDE, we propose the following formulation of the second-order equation within one integration call:

[dx/dt; dv/dt] = [v; F/m] (14)

We mark this as “func” in Fig. 2. This means we can augment ODE or SDE derivative space by concurrently calculating the accelerations and velocities, allowing simultaneous integration of velocities and positions.

Implementation details on velocity modeling. As in Eq. (12) and (14), NeuralMD requires both initial coordinates and velocities to predict the trajectories of coordinates and velocities. However, the MISATO dataset we are using does not include velocity information. To handle this issue, we introduce two solutions: (1) We first estimate the initial velocity as the positional momentum, v_t^{(l)} = (x_t^{(l)} - x_{t-1}^{(l)}). Optionally, we can apply a velocity initial mapping function to get a surrogate initial velocity: v_t^{(l)} = BindingNet-Ligand(f^{(l)}, x_t^{(l)} - x_{t-1}^{(l)}). The decision to include a velocity mapping is treated as a binary hyperparameter. (2) During the integration step in (14), we can also take an extra module to refine the velocity. More specifically, we have v_t^{(l)} = v_t^{(l)} + α * BindingNet-Ligand(f^{(l)}, x_t^{(l)}), where α is the velocity refinement coefficient. This is not displayed from Fig. 2 for brevity, yet all critical details are provided to ensure reproducibility.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This paper introduces a machine learning approach for molecular dynamics (MD) simulations of protein-ligand complexes.\nSpecifically, the proposed method incorporates two main components: a multi-scale SE(3)-equivariant neural network and\ntwo different neural differential equation solvers inspired by Newton's second law and Langevin dynamics. The method is\nvalidated on the MISATO dataset, with a custom set of baselines and evaluation metrics demonstrating its advantages in\nboth speed and accuracy. However, the reviewers raised several concerns regarding the experimental setup, choice of\nbaselines, methodology, and experimental outcomes:\n1. Baselines\nThe study evaluates the proposed method against only three deep learning-based baselines, ensuring a fair comparison by\nusing the same backbone across different methods. However, could the authors include other MD methods, even those with\ndifferent types of backbones, to present more competitive results? While the paper highlights a 1000x speed improvement\nover standard numerical MD, it lacks direct comparisons of both accuracy and speed with traditional computational MD\napproaches.\n2. Model Comparisons and Ablation Studies:\nWhile the method employs multi-scale modeling of protein-ligand complexes, many other approaches exist for handling\nequivariant complex inputs. Comparisons with alternative network types, such as those developed for docking[1,2,3] or\nstructure-based drug design (SBDD)[4] scenarios, would strengthen the study. Additionally, ablation studies are absent,\nmaking it difficult to assess the contributions of different network components or the multi-scale design to the task.\n[1] EquiBind: Geometric Deep Learning for Drug Binding Structure Prediction\n[2] FABind: Fast and Accurate Protein-Ligand Binding\n[3] DiffDock-Pocket: Diffusion for Pocket-Level Docking with Side Chain Flexibility\n[4] MolCRAFT: Structure-Based Drug Design in Continuous Parameter Space\n3. Task Design and Protein Invariance Assumption\nThe evaluation metrics primarily focus on the accuracy of position and distance, without addressing binding affinity or\nenergy-related measures. This makes it challenging to understand how the simulation results directly contribute to drug\ndiscovery applications. The paper also assumes the protein structure to be invariant, which might oversimplify the dynamics\nand deviate from real-world MD scenarios involving protein-ligand complexes.\nFurthermore, the claim that \"all prior ML approaches for MD simulation are limited to single-system and not protein-ligand\ncomplex\" is not appropriate. In fact, AI2BMD [5] is a ML method to consider complex and consider both protein and ligand\nchanges in the MD process, which has been published in Nature. More discussion about the differences and benefits of the\nproposed method against AI2BMD should be included.\n[5] Ab initio characterization of protein molecular dynamics with AI2BMD\n4. Case Study on 4G3E\nWhile NIK is an interesting target, the case study lacks in-depth analysis. Additional results are necessary to demonstrate\nhow the proposed method enhances drug binding affinity and selectivity. Experimental validation (e.g., wet-lab experiments)\nor highly accurate computational methods (e.g., FEP, TI) should be included.\n5. About Simultaneous Calculation of Acceleration and Velocity\nThe paper mentions:\n\"We augment derivative space by concurrently calculating the accelerations and velocities, allowing simultaneous\nintegration of velocities and positions.\"\nHowever, it does not provide sufficient algorithmic details. Numerical integration of positions requires velocity values at\ncorresponding time steps, while velocity calculations require prior integration of accelerations. Please clarify how\n“simultaneous” calculations are achieved and describe the numerical integration method, including the choice of\ndiscretization step size.\n6. About Efficiency\nWhy does the NeuralMD algorithm reduce the need for discretization compared to the DenoisingLD method, enabling\nconformational prediction with fewer time steps? Please provide a detailed explanation.\n7. Pretraining\nCan pretraining on data annotated with energy/force labels further improve the performance?\n8. Suggestions for Improving Figures and Names\nThe information provided in Figure 1 is limited, as it only presents the formula for the Velocity Verlet algorithm. It is\nrecommended to enhance the figure by including additional information, such as comparisons with GNN-MD and\nDenoisingLD. This would help emphasize the innovative aspects of the proposed method. The design and color scheme of\nFigure 2a could be improved for better aesthetics and clarity.\nBesides, I suggest the authors consider using another name of BindingNet. Because it has been used by a protein-ligand\ncomplex data, which may easily cause confusions.", "reviewer_2": "The work by Liu et al. represents ongoing efforts to expedite MD and related computational methods for drug discovery\napplications using ML. It seems that they achieved significant improvement on efficiency of MD simulations. However, this\nwork lacks a few things for it to be useful in a drug discovery campaign.\n1. It is not clear that with the NeuralMD, one can perform interaction analysis throughout trajectories.\n2. As the authors point out, proteins being fixed during the simulation is a major drawback, as one usually resorts to MD\nsimulations for exactly the cases in which protein structures are flexible. This assumption severely limits the utility of the\nmethod in drug discovery.\n3. Training dataset is in the order of 10,000, which doesn't seem to be enough. I wonder if one can augment it by generating\nprotein-ligand complex structures by docking.\nI would like to see how these issues are addressed before recommending the manuscript for publication.", "reviewer_3": "In this study, the authors developed a multi-grained symmetric differential equation model for learning protein-ligand binding\ndynamics. The motivation of this work is that the conventional molecular dynamics method is too slow and too\ncomputationally expensive for predicting protein-ligand structures. The method, NeuralMD, outperforms other ML\napproaches, such as GNN-MD, DenoisingLD, etc, as for the accuracy, stability, and performances. The results look good\ndue to the comparisons. However, for this reviewer, who is not an expert of the ML field, there are many fundamental\nquestions about the manuscript in the current form.\n1) Purpose of this work\nWhat is the main purpose of this work? Do they want to predict the binding processes or just the binding poses? Since they\noften use the word, “trajectories”, like MD simulations, I was confused for their main research purposes.\n2) Initial structures\nAccording to the movies on their github, the motions of ligand are relatively small. Though they are flexible, their positions\nrelative to the protein do not change so much. What are the initial structures in the predictions?\nIf they separated a protein and a ligand, how they can predict the protein-ligand complexes?\nWhat are their target conformational spaces to be sampled by ML?\n3) Binder and non-binder distinctions\nCan they distinct binder ligands or non-binder ligands in their methods?\n4) RMSF in Figure 3\nWhat is the ground-truth RMSF? Is it MD simulation result or B-factor in X-ray structure?\nWhy the RMSF in ML is generally so small?\nWhat is the physical meaning of the RMSF?\nWhy the ground-truth RMSF in 1KTI is so large?\nThese questions are fundamental and useful to increase the readability of this manuscript to non-expert chemists. I want to\nhear the authors’ replies."}, "author_response": "Revisions made to manuscript NCOMMS-24-74647\nWe would like to extend our sincere appreciation to the editor and reviewers for their invaluable\ncomments and suggestions, which have significantly enhanced the quality of this manuscript. In\nresponse, we want to highlight four key points:\n1. In this community, many papers have been using E(3)-equivariant networks for modeling\nthe molecule geometry. However, they missed the key difference between E(3)-equivariance\nand SE(3)-equivariance, especially reflection anti-symmetry, which is very important in the\nbinding tasks. We highlighted how we handled this using the three granularities of vector\nframes in our NeuralMD pipeline, in Sections 4, B, and C (with detailed proof).\n2. We highlighted that there are two different research lines in AI for MD simulation, one is ML\nfor Potential Energy and Force Learning for MD Simulation, and the other one is\nML for Trajectory Learning for MD Simulation. We listed the details and comparison\nin Section A. Note that our current version of NeuralMD belongs to the second category.\n3. As discussed in Section 3, we currently adopt a semi-flexible setting due to computational\nlimitations. This setting has also been used in prior works [1,2]. While we acknowledge that\nthis is not a perfect solution, we view NeuralMD as an initial step in this direction. Our\ngroup is actively pursuing follow-up work to further advance this line of research.\n4. Meanwhile, NeuralMD still supports primary interaction analysis, as demonstrated in the\n4G3E case study in Section 2.8 and the 10 binding complex examples in the supplementary\nfiles(alsoonthewebsite). TheseobservationssuggestthatNeuralMDisapromisingdirection\nfor future exploration, especially from an AI and physics perspective.\nTo summarize, while we acknowledge the current limitations of NeuralMD, we believe it addresses\nimportant challenges in the community and offers a valuable step forward. We are grateful for\nthe constructive feedback, and we have meticulously revised the manuscript (marked in green)\nin accordance with the insightful comments. We would also be grateful if they could recognize\nthe contributions of our proposed multi-grained vector frame modeling approach. Please find our\ndetailed responses below.\nThank you for raising these detailed and insightful questions about our work. We have carefully\naddressed each point with additional explanations as outlined below.\nReviewer Comment1.1 — Baselines. The study evaluates the proposed method against only\nthree deep learning-based baselines, ensuring a fair comparison by using the same backbone across\ndifferentmethods. However,couldtheauthorsincludeotherMDmethods,eventhosewithdifferent\ntypes of backbones, to present more competitive results? While the paper highlights a 1000x speed\nimprovement over standard numerical MD, it lacks direct comparisons of both accuracy and speed\nwith traditional computational MD approaches.\nReply: Thank you for raising these questions. We would like to answer the key points below.\n• First would like to point out that this is a new research direction, on ML for trajectory learning\nfor MD simulation, as summarized in Section A.\n• In our work, the 1000× speedup is measured relative to the traditional numerical MD method,\nspecifically GROMACS. This reflects the speed comparison. Meanwhile, the MISATO dataset [3]\nprovides the binding dynamics data used in our experiments, which were generated using these\nnumerical MD simulations. In other words, while we achieve significant acceleration, the ground-\ntruth trajectories are obtained from accurate numerical methods.\n• By ‘... other MD methods, even those with different types of backbones...’ Do you mean other\nAI-MD methods? If so, then it is related to Comment1.2, and we would like to discuss the details\nbelow accordingly.\nReviewer Comment1.2 — Model Comparisons and Ablation Studies: While the method\nemploysmulti-scalemodelingofprotein-ligandcomplexes,manyotherapproachesexistforhandling\nequivariant complex inputs. Comparisons with alternative network types, such as those developed\nfor docking[1,2,3] or structure-based drug design (SBDD)[4] scenarios, would strengthen the study.\nAdditionally, ablation studies are absent, making it difficult to assess the contributions of different\nnetwork components or the multi-scale design to the task.\n[1] EquiBind: Geometric Deep Learning for Drug Binding Structure Prediction [2] FABind: Fast\nandAccurateProtein-LigandBinding[3]DiffDock-Pocket: DiffusionforPocket-LevelDockingwith\nSide Chain Flexibility [4] MolCRAFT: Structure-Based Drug Design in Continuous Parameter\nSpace\nReply:\nThank you for mentioning these works.\nFirstly, we would like to acknowledge that we are aware of this line of research, including [1, 2, 4].\nHowever, these works are not directly applicable to our setting due to a fundamental mathematical\nissue: They employ E(3)-equivariant networks, following EGNN [4], which are not suitable for bind-\ning tasks which require the model to be E(3)-equivariant and reflection-antisymmetry. We provide a\nantisymmetry.\ndetailed explanation of this limitation in the paragraph SE(3)-equivariance, E(3)-equivariance, and\nreflection-antisymmetry below.\n• EquiBind, E3-equivariant in Section 3, link.\n• FABind, E3-equivariant in Section 3.2, link. Besides, this is a parallel work.\n• MolCRAFT, E3-equivariant in Eq 25, link. Besides, this is also a parallel work.\nSecondly, from a modeling perspective, our approach—based on SE(3)-equivariance and reflection-\nantisymmetry—can be seen as an extension of the vector frame concept introduced in AlphaFold2 [5].\nIn AlphaFold2, the authors construct a local frame for each residue using three backbone atoms. We\nextend this idea by designing vector frames at three levels of granularity, as shown in Figure 2a of the\nmanuscript: a ligand frame, a protein frame, and a complex frame.\nContinuing the discussion on geometric modeling, we believe there are three reasons why adding the\nwork [3] as a baseline is not necessary: (1) The DiffDock-Pocket paper uses the Tensor Field Net-\nwork (TFN) [6], which was the first architecture to incorporate SE(3)-equivariance and reflection-\nantisymmetry. However, more expressive and effective approaches have since been developed. (2) In\nrecent top-tier ML conferences, vector frame–based modeling is the state-of-the-art for protein geo-\nmetric modeling, as demonstrated by works such as FrameDiff [7], FrameFlow [8], FoldFlow [9], and\nFoldFlow2 [10]. (3) DiffDock-Pocket is a parallel work that has not yet been peer-reviewed.\nSE(3)-equivariance, E(3)-equivariance, and SE(3)-equivariance and reflection-antisymmetry.\nSE(3)-equivariance encompasses translation and rotation in Euclidean space but does not account for\nreflection or chirality equivariance. Then, as shown in Figure 1, there are two cases:\n• SE(3)-equivariance and reflection-symmetry, or equivalently E(3)-equivariance. In the EGNN pa-\nper [4], the authors explicitly state: “It is often desired that predictions on these tasks are either\nequivariant or invariant w.r.t. E(3) transformations.” However, the paper does not provide a\nqualitative or quantitative assessment of how often this assumption holds in practice.\n• SE(3)-equivarianceandreflection-antisymmetry. ThisstandsincontrasttoE(3)-equivariance,and\nthere are two reasons that we should stick to this when applying the geometric modeling.\n1. In the supplementary information of the AlphaFold2 paper [5] – specifically, Supplementary\nFigure9(link)–theauthorsdemonstratethatnaturalmoleculesarenotreflection-symmetric,\nbut rather exhibit reflection-antisymmetry.\n3\n2. When modeling protein–ligand binding tasks, it is important to incorporate both SE(3)-\nequivarianceandreflection-antisymmetryfortheproteinandligand. UsingE(3)-equivariance\ninstead can lead to undesirable behavior. For example, if the protein is fixed and the small\nmolecule is reflected, an E(3)-equivariant model would produce reflected predictions—such\nas forces—which is physically incorrect and inconsistent with the nature of molecular inter-\nactions.\nBesides, in our supplementary information, we explicitly explain how we can guarantee the SE(3)-\nequivariance and reflection-antisymmetry in Section C. Please feel free to check if you are interested.\nReviewer Comment1.3 — 3. Task Design and Protein Invariance Assumption The evaluation\nmetricsprimarilyfocusontheaccuracyofpositionanddistance,withoutaddressingbindingaffinity\nor energy-related measures. This makes it challenging to understand how the simulation results\ndirectly contribute to drug discovery applications. The paper also assumes the protein structure\nto be invariant, which might oversimplify the dynamics and deviate from real-world MD scenarios\ninvolving protein-ligand complexes. Furthermore, the claim that ”all prior ML approaches for MD\nsimulation are limited to single-system and not protein-ligand complex” is not appropriate. In fact,\nAI2BMD [5] is a ML method to consider complex and consider both protein and ligand changes in\nthe MD process, which has been published in Nature. More discussion about the differences and\nbenefits of the proposed method against AI2BMD should be included.\n[5] Ab initio characterization of protein molecular dynamics with AI2BMD\nReply: Thank you for raising the questions on the problem setting, energy, and related works.\nProblem formulation. First of all, we would like to mention that the problem formulation is:\nL = E (cid:2) ||xˆ (l) −x (l) || (cid:3) . (1)\nt t t\nBTW. This is Equation 13 in the manuscript. We would like to point out a few key differences or\nattributes of our problem formulation:\n• This is different from the other research line on energy and force prediction for MD simulation.\n• There are two key differences between energy prediction for MD and trajectory prediction for\nMD: (1) Energy prediction takes each conformation and energy as IID, while trajectory learning\noptimizes the conformations along the whole trajectory, enforcing the temporal relation. (2) The\ntrajectory learning is agnostic to the magnitude of the timesteps, and energy prediction can be\nsensitive to longer-timestep MD simulations.\n• This research line includes CGDMS [11], DiffMD [12], DFF [13], CG-MD [14].\n• We listed the detailed comparison in Section A.2: the comparison between ML for Potential\nEnergy and Force Learning for MD Simulation and ML for Trajectory Learning for MD\nSimulation..\nAbout the problem formulation of the semi-flexible setting. Assuming the protein stays fixed\nbecause:\n• This is due to the large volume of proteins.\n• Modeling both the ligand and protein and applying the second-order integration for simulating\nthe dynamics goes above the GPU memory limitation.\n4\n• This assumption is also common in the top-tier community, e.g., DiffDock [1]. BTW, we also\nfound that the DiffDock-PP [2] shared by you using the inflexible setting, and they claimed that\n“This assumption is often realistic (Vakser, 2014) and even leads to improved results for most\ninteracting proteins (Desta et al., 2020).”\nThen on the analysis. We acknowledge that NeuralMD has certain limitations, as we discussed in\nSection 3 of the manuscript, including GPU limitations for the flexible setting. However, NeuralMD still\npossesses the ability for certain primary interaction analysis.\n• We provide a case study on 4G3E, as described in the manuscript Section 2.8, as well as in\nComment1.4.\n• We provide the demonstration for 10 binding complex examples in the supplementary files (also\non the website).\nThese observations reveal that NeuralMD is a promising direction for future exploration, especially from\nthe AI and physics perspective.\nAbout the related works. We would like to kindly point out that the first version of our work was\nreleased in September 2023, and AI2BMD is indeed a parallel work.\nBesides, we also checked if AI2BMD open-sourced the codes, and the authors stated on Jan 14th 2025\nthat “The current version support single chain protein simulations, while more functions will be added\nand released later”. Please check this link for more details.\nMeanwhile,wewouldalsoliketopointout,aswillbediscussedinComment1.7,theproblemformulation\nof AI2BMD is ML for Potential Energy and Force Learning for MD Simulation, while ours is in\nthe ML for Trejactory Learning for MD Simulation research line. There are two key differences\nbetween energy and trajectory prediction for MD (we also discussed them in the comments above): (1)\nEnergy prediction takes each conformation and energy as IID, while trajectory learning optimizes the\nconformations along the whole trajectory, enforcing the temporal relation. (2) The trajectory learning\nis agnostic to the magnitude of the timesteps, and energy prediction can be sensitive to longer-timestep\nMD simulations.\nPlease feel free to check more discussions in Section A.2. BTW. We have added this discussion to the\nrevised manuscript in Section A.2.\nReviewer Comment1.4 — Case Study on 4G3E. While NIK is an interesting target, the\ncase study lacks in-depth analysis. Additional results are necessary to demonstrate how the pro-\nposedmethodenhancesdrugbindingaffinityandselectivity. Experimentalvalidation(e.g., wet-lab\nexperiments) or highly accurate computational methods (e.g., FEP, TI) should be included.\nReply: Thank you for your valuable feedback.\nWe agree that NIK is an important kinase involved in the regulation of immunity and inflammation,\nand it has also been identified as a potential therapeutic target in glioma, melanoma, and prostate\ncancer [15,16]. Your suggestion to perform wet-lab validation or employ highly accurate computational\nmethods such as Free Energy Perturbation (FEP) and Thermodynamic Integration (TI) is indeed a\nstandard and rigorous approach for confirming improvements in the drug binding affinity and selectivity\nas simulated by our NeuralMD model. However, we plan to continue this work in a future project. Our\n5\ndecision is based on the following two reasons:\n• The application of FEP/TI does not provide stronger evidence for demonstrating the\nacceleration of molecular dynamics (MD) simulations achieved by our model compared\nto what is shown in Figure 3, and the current architecture of our model is not yet well-\nsuited for reliably executing FEP/TI simulations. We agree that FEP and TI are powerful\ntoolsforcomparing the relativebinding affinitiesof differentligands to thesame protein target. In\nour case, this would require identifying additional molecules that bind to NIK for comparison with\nCMP1, a6-alkynylindoline. However, the primaryfocusof thiswork istointroduceourNeuralMD\nmodel, which achieves over a 1000-fold improvement in MD simulation speed across various\nprotein–ligand systems (Figure 3). We believe this demonstration provides stronger evidence of\nthe model’s capability in accelerating MD simulations than a comparison of the relative binding\nfree energies of different ligands for a single target. Moreover, our current model architecture is\nnot well-suited to support the execution of FEP/TI simulations. However, we recognize the value\nof such comparative studies and plan to include them in future work. We are currently applying\nthe model to novel molecules not present in the MISATO dataset, and extending the model to\nintegrate FEP/TI represents a promising direction to further assess both MD simulation speedup\nand accuracy against known protein targets such as NIK.\n• The experimental validation results related to NIK have already been discussed in our\npaper, and additional wet-lab experiments require further resources. The binding between\nCMP1 and NIK has already been experimentally validated in the original publication that reports\nthe crystal structure 4G3E [17]. The structure itself is the result of wet-lab crystallographic work,\nand has been discussed in Section 2.8. We are currently in contact with wet-lab collaborators to\nscreen and evaluate additional molecules that bind to NIK, in parallel with MD simulation using\nour NeuralMD model.\nThe main goal of our paper is to introduce our NeuralMD model, which can be applied across various\nprotein–ligand systems and offers over a 1000-fold speedup in MD simulations. Therefore, we did\nnot devote extensive space to NIK-specific analysis. However, we sincerely appreciate your thoughtful\nsuggestion and fully intend to incorporate both highly accurate computational validations and wet-lab\nexperiments in our future research.\nReviewer Comment1.5 — 5. About Simultaneous Calculation of Acceleration and Velocity\nThe paper mentions: ”We augment derivative space by concurrently calculating the accelerations\nand velocities, allowing simultaneous integration of velocities and positions.” However, it does not\nprovide sufficient algorithmic details. Numerical integration of positions requires velocity values\nat corresponding time steps, while velocity calculations require prior integration of accelerations.\nPlease clarify how “simultaneous” calculations are achieved and describe the numerical integration\nmethod, including the choice of discretization step size.\nReply: Thank you for raising this question.\nThis formulation is indeed rooted in physics, as presented in Equation 14 in the manuscript. We are\nhappy to provide additional details here.\nF (l) (cid:90) t (cid:90) t\na(l) = τ , vˆ (l) = v (l) + a(l)dτ, xˆ (l) = x (l) + vˆ(l)dτ. (2)\nτ m t 0 τ t 0 τ\n0 0\nCase 1: non-simultaneous integration. In this case, we need to call the ODE/SDE function twice.\n(cid:2) (cid:3)\nThe first ODE/SDE function call corresponds to solving the acceleration integration with dv/dt =\n(cid:2) (cid:3)\nF/m . This means that if we are integrating for T timesteps, then we have\nv = ODE/SDE Integration(model(x )/m),\n1 0\nv = ODE/SDE Integration(model(x )/m),\n2 1\n(3)\n...\nv = ODE/SDE Integration(model(x )/m).\nT T−1\n(cid:2) (cid:3)\nThesecondODE/SDEcallfunctioncallcorrespondingtosolvingthevelocityintegrationwith dx/dt =\n(cid:2) (cid:3)\nv . Similarly, this means that if we are integrating for T timesteps, then we have:\nx = ODE/SDE Integration(v ),\n1 0\nx = ODE/SDE Integration(v ),\n2 1\n(4)\n...\nx = ODE/SDE Integration(v ).\nT T−1\nThus, we will call the ODE/SDE integration calls twice, as in Equations (3) and (4).\nCase 2: simultaneous integration. The case means that we can call the acceleration integration and\nvelocity integration simultaneously, as shown in Equation 14 of the manuscript, we utilize\n(cid:20) (cid:21) (cid:20) (cid:21)\ndx/dt v\n= . (5)\ndv/dt F/m\nThen we call the integration for T timesteps as:\n(cid:20) (cid:21) (cid:20) (cid:21)\nv model(x )/m\n1 = ODE/SDE Integration( 0 ),\nx v\n1 0\n(cid:20) (cid:21) (cid:20) (cid:21)\nv model(x )/m\n2 = ODE/SDE Integration( 1 ),\nx 2 v 1 (6)\n...\n(cid:20) (cid:21) (cid:20) (cid:21)\nv model(x )/m\nT = ODE/SDE Integration( T−1 ).\nx v\nT T−1\nThus, we will call the ODE/SDE integration calls just once, in comparison to twice in the first case\n(Equations (3) and (4)), resulting in a twofold increase in efficiency.\nWe added this in Section E.3 in the revised Supplementary Information.\nneed for discretization compared to the DenoisingLD method, enabling conformational prediction\nwith fewer time steps? Please provide a detailed explanation.\nReply: Thank you for raising the question on efficiency. We briefly described this in the previous\nmanuscript, and we revised it in Section 2.7 in the latest manuscript. Here is a polished revision:\nSinceallMLmethodsshareasimilarbackbonearchitecture(BindingNetanditsvariants),theirefficiency\ndifferences arise from their computational approaches and hyperparameters. GNN-MD is the most\nefficient, as it requires no integration steps. NeuralMD incorporates an augmented integration step,\nwith the optimal step size set to 0.5 of the snapshot interval, i.e., two integration steps are performed\nto predict the next snapshot, making it slightly slower than VerletMD. In contrast, DenoisingLD is the\nslowest method, as it uses a smaller optimal timestep of one-tenth the snapshot interval, requiring 10\nintegration steps to generate the next snapshot during inference.\nReviewer Comment1.7 — 7. Pretraining. Can pretraining on data annotated with energy/-\nforce labels further improve the performance?\nReply: Thank you for raising this question on pretraining. Firstly, a brief answer to your question is\nthat: yes, pretraining on energy/force prediction should be helpful. We have an ongoing project that is\ndoing this.\nMeanwhile, we would like to mention the key differences in energy/force prediction and trajectory\nprediction.\n• In our current problem formulation, we assume that only positional data is available. This follows\ntheML for Trejactory Learning for MD SimulationresearchlineofseveralDenoisingLDpapers\n(CGDMS [11], DiffMD [12], DFF [13], CG-MD [14]), as discussed in Section A.2.\n• Therefore, incorporating additional supervised labels such as energy or force would define a\ndifferent problem setting, as discussed in ML for Potential Energy and Force Learning for MD\nSimulation in Section A.2. That said, if the objective is to address the general problem of dy-\nnamics simulation, incorporating such additional data—like energy or force labels—would indeed\nbe beneficial.\nReviewer Comment1.8 — 8. Suggestions for Improving Figures and Names. The information\nprovided in Figure 1 is limited, as it only presents the formula for the Velocity Verlet algorithm. The\ninformation provided in Figure 1 is limited, as it only presents the formula for the Velocity Verlet algorithm. It is\nrecommended to enhance the figure by including additional information, such as comparisons\nwithGNN-MDandDenoisingLD.Thiswouldhelpemphasizetheinnovativeaspectsoftheproposed\nmethod. The design and color scheme of Figure 2a could be improved for better aesthetics and\nclarity.\nReply: Thank you for raising this comment. We have revised the figure, and the revision is in Figure 2\nbelow. We also uploaded the revised version in the manuscript. We would be happy to further refine\nthe figures based on any additional suggestions you may have.\ndynamic leads to an equilibrium state with lower energy. (a) Pipeline for second-order methods,\nVelocity Verlet and NeuralMD (ours). (b) Pipeline for GNN-MD. (c) Pipeline for Denoising-LD.\nReviewer Comment1.9 — Besides, I suggest the authors consider using another name of\nBindingNet. Because it has been used by a protein-ligand complex data, which may easily cause\nconfusions.(Remarks on code availability)\nReply: Thank you for your suggestion. We will rename BindingNet to FrameNet in the final version.\nWe would like to express our gratitude to Reviewer 2 for raising important questions regarding\ndataset clarification and problem formulation. These are indeed limitations of the current Neu-\nralMD work, as discussed in Section 3 of the manuscript. We consider them valuable directions for\nfuture exploration and provide further details in our responses below.\nReviewer Comment2.1 — 1. It is not clear that with the NeuralMD, one can perform\ninteraction analysis throughout trajectories.\nReply: Thank you for raising this comment on the interaction analysis.\nFirst, we would like to mention that NeuralMD has certain limitations, as we discussed in Section 3 of\nthe manuscript, including the dataset and GPU limitations for the flexible setting.\nHowever, NeuralMD still possesses the ability for certain primary interaction analysis:\n• We provide a case study on 4G3E, as described in the manuscript Section 2.8, as well as in\nComment1.4.\n• We provide the demonstration for 10 binding complex examples in the supplementary files (also\non the website).\nReviewer Comment2.2 — 2. As the authors point out, proteins being fixed during the\nsimulation is a major drawback, as one usually resorts to MD simulations for exactly the cases in\nwhich protein structures are flexible. This assumption severely limits the utility of the method in\ndrug discovery.\nReply: Thank you for raising this important question. In our current setup, we assume the protein\nremains fixed for the following reasons:\n• This is due to the large volume of proteins.\n• Modeling both the ligand and protein and applying the second-order integration for simulating\nthe dynamics goes above the GPU memory limitation.\n• This assumption is also common in the top-tier community, e.g., DiffDock [1]. BTW. We also\nfound that the DiffDock-PP [2] shared by Reviewer 1 using the inflexible setting, and they claimed\nthat“Thisassumptionisoftenrealistic(Vakser, 2014)andevenleadstoimprovedresultsformost\ninteracting proteins (Desta et al., 2020).”\nWe would like to highlight that memory efficiency is currently the primary limitation, and addressing\nit is a key focus of our future work. In particular, our group is actively exploring geometric modeling\ntechniques that better capture molecular rigidity, which can also be naturally extended to model protein\nresidues more effectively.\nto be enough. I wonder if one can augment it by generating protein-ligand complex structures by\ndocking.\nReply: Thank you for pointing out this issue. We acknowledge that the dataset insufficiency is one\ninherent limitation. To the best of our knowledge, MISATO [3] is by far the largest one with semi-\nquantum methods. Of course, we could try the docking algorithm, yet it will be another research\ndirection on dataset generation, along the MISATO. We mentioned this as one limitation for future\nexploration in Section 3 of the manuscript.\nWe would like to express our gratitude to Reviewer 3 for raising questions on the problem formula-\ntion and RMSF evaluation. We have carefully addressed each of these points in detail, as explained\nbelow.\nReviewer Comment3.1 — What is the main purpose of this work? Do they want to predict\nthe binding processes or just the binding poses? Since they often use the word, “trajectories”,\nlike MD simulations, I was confused for their main research purposes.\nReply: Thank you for raising this question. Our goal is to predict the trajectories of the protein-ligand\nbinding processes. We have mentioned this in Section 2.1 Problem Setting and Section 4 Methods\nof the manuscript.\nReviewer Comment3.2 — 2) Initial structures. According to the movies on their github,\nthe motions of ligand are relatively small. Though they are flexible, their positions relative to the\nproteindonotchangesomuch. Whataretheinitialstructuresinthepredictions? Iftheyseparated\na protein and a ligand, how they can predict the protein-ligand complexes? What are their target\nconformational spaces to be sampled by ML?\nReply: The initial positions are provided in the MISATO dataset.\nIn the MISATO dataset [3], they simulate the trajectories for 10ns, and drop the first 2ns. We take the\nposition at the first timestep as the input, and the goal is to predict the remaining positions along the\nMD trajectories.\nReviewer Comment3.3 — 3) Binder and non-binder distinctions. Can they distinct binder\nligands or non-binder ligands in their methods?\nReply: Thank you for raising this question. The problem definition of our NeuralMD paper is p(ligand\npositions at t+1, ...— initial positions at t).\nThe binder and non-binder distinction is another task: F(ligand, protein) = True/False.\nReviewer Comment3.4 — 4) RMSF in Figure 3 What is the ground-truth RMSF? Is it MD\nsimulation result or B-factor in X-ray structure? Why the RMSF in ML is generally so small?\nWhat is the physical meaning of the RMSF? Why the ground-truth RMSF in 1KTI is so large?\nReply: The ground-truth RMSF is provided by the MISATO dataset [3], which is the MD simulation\nresult.\nRMSF meaning. RMSF is short for Root Mean Square Fluctuation. It quantifies the mean square\ndeviation of each atom’s position from its average position over the trajectory, providing a measure of\n12\npositional fluctuations. Or formally, it is defined as:\n(cid:118)\n(cid:117) N\n(cid:117) (cid:116)E (cid:2) 1 (cid:88) ∥r −⟨r ⟩∥2 (cid:3) , (7)\nt i,t i\nN\ni=1\nwhere r is the i-th atom position at time t (ground truth or sampled by ML methods) and ⟨r ⟩ is the\ni,t i\naverage position for the i-th atom across time. We have this detailed in Section 2.6 in the manuscript.\nWhy RMSF in ML methods is small. ThatisthegeneralizationabilityofMLmodels. InFigure3, the\nexperiment setting corresponds to the Generalization of One Single Trajectory. In other words, for\neach trajectory, we train NeuralMD on the first 80 snapshots (positions) and generalize to the remaining\nsnapshots (the remaining positions along the MD trajectory). The RMSF values are small and tend to\nbe smaller because the first 80 snapshots fluctuations are smaller.\nRMSF in 1KTI. We have created demonstrations for the 10 binding complex examples in the supple-\nmentary files (also on the website). Notably, there is a significant jump for 1KT1 (ground-truth), which\nis due to the simulation results from the MISATO dataset [3].\n[1] G. Corso, H. St¨ark, B. Jing, R. Barzilay, and T. Jaakkola, “Diffdock: Diffusion steps, twists,\nand turns for molecular docking,” arXiv preprint arXiv:2210.01776, 2022.\n[2] M. A. Ketata, C. Laue, R. Mammadov, H. St¨ark, M. Wu, G. Corso, C. Marquet, R. Barzilay,\nand T. S. Jaakkola, “Diffdock-pp: Rigid protein-protein docking with diffusion models,” arXiv\npreprint arXiv:2304.03889, 2023.\n[3] T. Siebenmorgen, F. Menezes, S. Benassou, E. Merdivan, K. Didi, A. S. D. Mour˜ao, R. Kitel,\nP. Li`o, S. Kesselheim, M. Piraud et al., “Misato: machine learning dataset of protein–ligand\ncomplexes for structure-based drug discovery,” Nature Computational Science, pp. 1–12, 2024.\n[4] V. G. Satorras, E. Hoogeboom, and M. Welling, “E (n) equivariant graph neural networks,”\nin International conference on machine learning. PMLR, 2021, pp. 9323–9332.\n[5] J.Jumper,R.Evans,A.Pritzel,T.Green,M.Figurnov,O.Ronneberger,K.Tunyasuvunakool,\nR. Bates, A. Zˇ´ıdek, A. Potapenko et al., “Highly accurate protein structure prediction with\nalphafold,” nature, vol. 596, no. 7873, pp. 583–589, 2021.\n[6] N. Thomas, T. Smidt, S. Kearnes, L. Yang, L. Li, K. Kohlhoff, and P. Riley, “Tensor field\nnetworks: Rotation-and translation-equivariant neural networks for 3d point clouds,” arXiv\npreprint arXiv:1802.08219, 2018.\n[7] J. Yim, B. L. Trippe, V. De Bortoli, E. Mathieu, A. Doucet, R. Barzilay, and T. Jaakkola,\n“Se (3) diffusion model with application to protein backbone generation,” arXiv preprint\narXiv:2302.02277, 2023.\n[8] J. Yim, A. Campbell, A. Y. Foong, M. Gastegger, J. Jim´enez-Luna, S. Lewis, V. G. Satorras,\nB. S. Veeling, R. Barzilay, T. Jaakkola et al., “Fast protein backbone generation with se (3)\nflow matching,” arXiv preprint arXiv:2310.05297, 2023.\n[9] A. J. Bose, T. Akhound-Sadegh, G. Huguet, K. Fatras, J. Rector-Brooks, C.-H. Liu, A. C.\nNica, M. Korablyov, M. Bronstein, and A. Tong, “Se (3)-stochastic flow matching for protein\nbackbone generation,” arXiv preprint arXiv:2310.02391, 2023.\n[10] G. Huguet, J. Vuckovic, K. Fatras, E. Thibodeau-Laufer, P. Lemos, R. Islam, C.-H. Liu,\nJ. Rector-Brooks, T. Akhound-Sadegh, M. Bronstein et al., “Sequence-augmented se (3)-flow\nmatchingforconditionalproteinbackbonegeneration,”arXivpreprintarXiv:2405.20313,2024.\n[11] J. G. Greener and D. T. Jones, “Differentiable molecular simulation can learn all the pa-\nrameters in a coarse-grained force field for proteins,” PloS one, vol. 16, no. 9, p. e0256990,\n2021.\n[12] F.WuandS.Z.Li,“Diffmd: Ageometricdiffusionmodelformoleculardynamicssimulations,”\ninProceedings of the AAAI Conference on Artificial Intelligence, vol.37, no.4, 2023, pp.5321–\n5329.\n14\n[13] M. Arts, V. Garcia Satorras, C.-W. Huang, D. Zugner, M. Federici, C. Clementi, F. Noe,\nR. Pinsler, and R. van den Berg, “Two for one: Diffusion models and force fields for coarse-\ngrained molecular dynamics,” Journal of Chemical Theory and Computation, 2023.\n[14] X. Fu, T. Xie, N. J. Rebello, B. D. Olsen, and T. Jaakkola, “Simulate time-integrated\ncoarse-grained molecular dynamics with geometric machine learning,” arXiv preprint\narXiv:2204.10348, 2022.\n[15] K. M. Pflug and R. Sitcheran, “Targeting nf-κb-inducing kinase (nik) in immunity, inflam-\nmation, and cancer,” International journal of molecular sciences, vol. 21, no. 22, p. 8470,\n2020.\n[16] H. Teng, L. Xue, Y. Wang, X. Ding, and J. Li, “Nuclear factor κb-inducing kinase is a diag-\nnostic marker of gastric cancer,” Medicine, vol. 99, no. 5, p. e18864, 2020.\n[17] G. de Leon-Boenig, K. K. Bowman, J. A. Feng, T. Crawford, C. Everett, Y. Franke, A. Oh,\nM. Stanley, S. T. Staben, M. A. Starovasnik et al., “The crystal structure of the catalytic\ndomain of the nf-κb inducing kinase reveals a narrow but flexible active site,” Structure,\nvol. 20, no. 10, pp. 1704–1714, 2012.\n15"}, {"round": 2, "reviewer_comments": {"reviewer_1": "I appreciate the author's revisions to the paper and figures, as well as the renaming of the method. My questions regarding\nthe distinction between energy/force and MD tasks, the explanation of pretraining, and the subsequent experimental\nvalidation plan have been resolved.\nI also thank the author for the explanations regarding SE(3) equivariance, E(3) equivariance, and reflection-antisymmetry in\nthe network. However, even if these methods don't explicitly emphasize the concepts of SE(3) equivariance and reflection-\nantisymmetry, they also focus on multi-scale modeling of protein-ligand interactions, which is a central aspect of this work.\nWithout a comparison, it's difficult to assess the practical significance of the proposed network’s architectural innovations.\nTherefore, I encourage the authors to conduct such comparisons and provide corresponding discussions or empirical\nanalyses to better highlight the novelty and effectiveness of their design.", "reviewer_2": "As the authors admit, there are limitations to their proposed method. However, with their detailed response, I am convinced\nthat NeuralMD has utility in the current drug discovery practice, though there is quite a bit of room to improve. I recommend\nthe publication of the manuscript.", "reviewer_3": "In the discussion with reviewers, the limitation of this machine learning model becomes clearer. The rigid treatment of\nprotein is a significant approximation compared to the flexible molecular dynamics simulations. In this sense, their statement\nabout the comparison between the ML surrogate and MD is indeed an overstatement. For instance, in the abstract, they\nwrote:\nWe demonstrate the efficiency and effectiveness of NeuralMD, achieving over 1K× speedup compared to standard\nnumerical MD simulations.\nBut, considering the rigid treatment in the ML surrogate and the flexible treatment in the MD simulation, I cannot accept the\nstatement. The author needs to emphasize the limitation of the current model in the abstract and maintext and discuss the\nfuture perspectives about the flexible treatment of protein structures.\nAlso, the physical meaning of RMSF shown in Figure 3 is different from the RMSF in MD simulations, because of the\ndifferent treatment of protein flexibility. If such quantities are shown without reasonable explanations, the data will confuse\nreaders of this manuscript.\nAlso, the treatment of solvent is different between the ML surrogate model and the MD simulation. In the MD simulation,\nsolvent molecules are explicitly included, while the ML surrogate does not.\nI suggest that the author compared these two different methods with the same conditions and approximations. Otherwise,\njust the numerical comparison of the computational time is not meaningful."}, "author_response": "Revisions made to manuscript NCOMMS-24-74647\nWe would like to extend our sincere appreciation to the editor and reviewers for their invaluable\ncomments and suggestions, which have significantly enhanced the quality of this manuscript.\nIn response, we want to highlight four key points:\n1. We conducted ablation studies by changing the three granularities to two granularities.\n2. We modified the statements on rigidity and acceleration of using NeuralMD.\n3. We clarified the main usage of NeuralMD, which is complementary to the existing classic MD\nmethods.\nTo summarize, while we acknowledge the current limitations of NeuralMD, we believe it addresses\nimportant challenges in the community and offers a valuable step forward. We are grateful for\nthe constructive feedback, and we have meticulously revised the manuscript (marked in green)\nin accordance with the insightful comments. We would also be grateful if they could recognize\nthe contributions of our proposed multi-grained vector frame modeling approach. Please find our\ndetailed responses below.\nReviewer 1\nThank you for acknowledging our revisions to the statements, figures, and methods – particularly\nour clarifications regarding the distinction between energy/force and MD tasks, the explanation of\npretraining, and the subsequent experimental validation plan. Your remaining question on SE(3)-\nequivariance has been addressed below.\nReviewer Comment1.1 — I also thank the author for the explanations regarding SE(3)\nequivariance,E(3)equivariance,andreflection-antisymmetryinthenetwork. However,evenifthese\nmethodsdon’texplicitlyemphasizetheconceptsofSE(3)equivarianceandreflection-antisymmetry,\nthey also focus on multi-scale modeling of protein-ligand interactions, which is a central aspect of\nthis work. Without a comparison, it’s difficult to assess the practical significance of the proposed\nnetwork’s architectural innovations.\nTherefore, I encourage the authors to conduct such comparisons and provide corresponding discus-\nsions or empirical analyses to better highlight the novelty and effectiveness of their design.\nReply:\nThank you for raising this question on multi-scale modeling. In comparison, we list the main differences\nofthebaselinesbelow. Onethingwewouldliketohighlightisthatmostofthesebaselines(e.g.,FABind,\nMolCRAFT) are parallel work and their underlying physics is wrong. They applied a simple paradigm\nas the interaction modeling, which we added as an ablation study.\n• EquiBind is just applying E(3)-equivariant on ligands and proteins for the binding prediction.\nThere is no multi-scaling modeling. link.\n1\ni i\ni+1 i+1\nCa C Ca\nC C i-1 N N Ca C i-1 C C Ca\nC C C C\nC C Ca C Ca C C\nN\ndnagiL nietorP\n(a) Three granularities of vector frame basis in NeuralMD, Ours\n(a.1) Ligand Frame (a.2) Protein Frame (a.3) Complex Frame\ni i\ni+1 i+1\nCa C Ca\nC C i-1 N N Ca C i-1 C C Ca\nC C C C\nC C Ca C Ca C C\nN\ndnagiL nietorP\n(b) Two granularities of vector frame basis in NeuralMD, Ablation\nCa C\nN Ca C\nN\nCa C\nN\n(b.1) Ligand Frame (b.2) Protein Frame (b.3) Complex Frame = Ligand Frame + Protein Frame\nFigure 1: Fig. (a) corresponds to the three granularities of vector frame basis in NeuralMD (Ours). Fig.\n(b)liststhethreegranularitiesofvectorframebasisintheablationstudies, i.e., NeuralMD(Ablation). The\nmain difference lies in the third granularity (if we take (b.3) as one): Fig. (a.3) utilizes the backbone atoms\ninpocketstomodelthebindingcomplex,whileFig. (b.3)considersallthebackboneatomsandligandatoms\nfor frame construction.\n• FABind is applying E(3)-equivariant on two hierarchies, link. Besides, this is a parallel work.\n– It includes a pocket prediction module to determine the binding site. Note that this is not\nrequired in our case, since we include the ligand information as the initial input.\n– Then goes to the FABind layer, which is similar to the EquiBind layer, using E(3)-equivariant\nlayers on ligands and proteins for their interaction prediction.\n• MolCRAFT,directlymodelstheligandatomsandproteinatoms,i.e.,nohierarchicalormulti-scale\nmodeling, link. Besides, this is also a parallel work.\nBoth EquiBind and FABind perform message passing between ligand and receptor, using either joint-\nthen-separateorseparate-then-jointmulti-scalemodeling. Incontrast,ourNeuralMDintroducesathree-\nlevel hierarchical design. As shown in Figure 1, Fig. (a) illustrates the three vector frame granularities\nin NeuralMD (Ours), while Fig. (b) shows the ablation version. The key difference is in the third\ngranularity (if we take (b.3) as one): Fig. (a.3) uses pocket backbone atoms, whereas Fig. (b.3)\nincludes all backbone and ligand atoms for frame construction.\nWe conducted an ablation study of the new NeuralMD architecture using 10 single-traj experiments,\nand the results are listed in Table 1. We mark our proposed hierarchical structure (Fig. a in Figure 1)\nas NeuralMD (ODE), Ours and NeuralMD (SDE), Ours. We mark the structures in ablation\nstudies (Fig. b in Figure 1) as NeuralMD (SDE), Ablation and NeuralMD (SDE), Ablation. For\nthe general performance, the hierarchical structures obtain (ours) obtain better performance than the\nnon-hierarchical architectures (ablation studies).\nwith optimal training loss are reported. Four evaluation metrics are considered: MAE (˚A, ↓), MSE (↓),\nMatching(↓), and Stability (%, ↑).\nNeuralMD(ODE),Ours NeuralMD(ODE),Ablation NeuralMD(SDE),Ours NeuralMD(SDE),Ablation\nMAE(↓) 2.252 2.259 2.260 2.600\nMSE(↓) 1.514 1.506 1.514 1.766\n5WIJ\nMatching(↓) 0.464 0.509 0.615 0.735\nStability(↑) 82.046 77.125 67.464 57.138\nMAE(↓) 1.878 1.879 2.158 2.156\nMSE(↓) 1.263 1.255 1.455 1.455\n4ZX0\nMatching(↓) 0.428 0.475 0.696 0.886\nStability(↑) 81.401 77.489 59.109 52.255\nMAE(↓) 3.858 3.624 3.395 3.435\nMSE(↓) 2.651 2.471 2.309 2.350\n3EOV\nMatching(↓) 1.062 1.160 0.962 1.095\nStability(↑) 47.328 42.652 50.108 44.788\nMAE(↓) 3.656 7.933 3.765 4.049\nMSE(↓) 2.400 5.204 2.501 2.650\n4K6W\nMatching(↓) 0.928 10.751 1.076 1.352\nStability(↑) 49.438 41.319 49.700 42.515\nMAE(↓) 6.675 6.706 6.646 6.664\nMSE(↓) 4.176 4.181 4.141 4.163\n1KTI\nMatching(↓) 0.337 0.303 0.167 0.214\nStability(↑) 86.430 87.932 98.508 95.103\nMAE(↓) 1.924 1.929 2.061 2.119\nMSE(↓) 1.280 1.285 1.356 1.403\n1XP6\nMatching(↓) 0.537 0.588 0.615 0.768\nStability(↑) 75.533 70.409 69.423 62.079\nMAE(↓) 6.957 6.975 7.038 6.977\nMSE(↓) 4.597 4.607 4.679 4.593\n4YUR\nMatching(↓) 0.584 0.700 0.749 0.820\nStability(↑) 69.775 64.365 60.344 55.515\nMAE(↓) 2.191 2.672 2.345 2.365\nMSE(↓) 1.453 1.771 1.536 1.564\n4G3E\nMatching(↓) 0.505 0.698 0.521 0.594\nStability(↑) 71.436 65.011 68.729 62.840\nMAE(↓) 3.921 3.952 3.842 3.974\nMSE(↓) 2.652 2.660 2.601 2.691\n6B7F\nMatching(↓) 0.459 0.543 0.741 0.622\nStability(↑) 75.692 71.733 57.917 63.308\nMAE(↓) 3.039 2.990 3.132 3.390\nMSE(↓) 1.999 1.995 2.078 2.263\n3B9S\nMatching(↓) 0.659 1.055 0.444 0.453\nStability(↑) 76.065 56.667 77.801 75.093\nReviewer 2\nWe sincerely thank Reviewer 2 for the thoughtful evaluation and encouraging comments. We\nappreciate your recognition of NeuralMD’s potential utility in current drug discovery practices, as\nwell as your acknowledgment of our detailed response. We also agree that there remains room for\nimprovement, and we see this as an exciting direction for future research. Thank you again for\nrecommending our manuscript for publication.\nReviewer 3\nWe would like to express our gratitude to Reviewer 3 for raising concerns on the clarification of the\nlimitation, the RMSF annotation, as well as the detailed discussion on comparison with the classic\nMD methods. We have revised the statements and clarified them in the revised manuscript. Please\ncheck the details below.\nThe author needs to emphasize the limitation of the current model in the abstract and maintext\nand discuss the future perspectives about the flexible treatment of protein structures.\nReply: Thank you for highlighting the potential impact of the semi-flexible/rigidity assumption. We\nhave removed the unclear statement and highlighted the rigidity assumption in the revised manuscript\nand marked in green.\nReviewer Comment3.2 — Also, the physical meaning of RMSF shown in Figure 3 is different\nfrom the RMSF in MD simulations, because of the different treatment of protein flexibility. If\nsuch quantities are shown without reasonable explanations, the data will confuse readers of this\nmanuscript.\nReply: Thankyouforraisingthisquestion. Inourproject,wearecalculatingtheRMSFonlyconsidering\nthefluctuationsoftheligands,andwerenameittoRMSF-Ligandintherevisedmanuscriptandmarked\nin green.\nReviewer Comment3.3 — Also,thetreatmentofsolventisdifferentbetweentheMLsurrogate\nmodel and the MD simulation. In the MD simulation, solvent molecules are explicitly included,\nwhile the ML surrogate does not.\nI suggest that the author compared these two different methods with the same conditions and\napproximations. Otherwise, just the numerical comparison of the computational time is not mean-\ningful.\nReply:\nThank you for bringing this up. It’s an important point, and we appreciate the opportunity to clarify.\nAs mentioned, the training dataset for NeuralMD is based on the MISATO dataset [1]. The molecular\ndynamics (MD) simulations in MISATO were carried out using Amber20, where the protein–ligand\ncomplexes were solvated in the TIP3P water model and neutralized with Na+ and Cl- ions [1]. As\na result, NeuralMD was trained to model protein–ligand dynamics under aqueous conditions. When\napplied to other complexes in the test set, its predictions implicitly assume similar solvent environments,\ni.e., in an in-distribution setting (in contrast to out-of-distribution). This is important when interpreting\nthe results.\nIn other words, NeuralMD and Amber20/GROMACS are different in the underlying principles.\n• NeuralMD is a machine learning (ML) model trained on MD trajectories of protein–ligand com-\nplexes in explicit solvent, which are typically generated from traditional MD methods like Am-\nber/GROMACS. Its simulations/predictions reflect patterns learned from such data, with solvent\neffects implicitly captured from the training dataset rather than explicitly simulated.\n• Additionally, the NeuralMD SDE (Eq. 11 in the manuscript) includes a stochastic term that\nimplicitly captures the effects of the solvent environment, even though explicit solvent information\nis not available in the MISATO dataset. If such solvent information were provided, it could be\ndirectly incorporated into the modeling. However, this remains a limitation of the current dataset.\n4\n• In contrast, Amber20 and GROMACS are traditional physics-based MD engines that use defined\nforce fields and explicitly model solvent molecules under user-specified simulation conditions.\nGiven these differences, we believe it may not be appropriate to directly apply traditional MD parameter\nsettings to NeuralMD. Nevertheless, we fully agree with your suggestion: NeuralMD would benefit from\nadditional training data generated under well-controlled MD protocols to further improve its simula-\ntion/prediction accuracy.\nRegarding the comparison of different methods under the same conditions and approximations, we\nensuredthatthesettingswerealigned. WhilethetrainingdataofNeuralMDcomefromAmber20-based\nMD simulations (with TIP3P water, as well as Na+ and Cl- ions neutralization), we used GROMACS\nfor benchmark simulations under comparable conditions for the protein-ligand complex, applying the\nAMBER99SB-ILDN force field, solvating with TIP3P water, and neutralizing with ions.\nReferences\n[1] T. Siebenmorgen, F. Menezes, S. Benassou, E. Merdivan, K. Didi, A. S. D. Mour˜ao, R. Kitel,\nP. Li`o, S. Kesselheim, M. Piraud et al., “Misato: machine learning dataset of protein–ligand\ncomplexes for structure-based drug discovery,” Nature Computational Science, pp. 1–12, 2024.\n5"}, {"round": 3, "reviewer_comments": {"reviewer_3": "The authors have improved the manuscript significantly."}, "author_response": ""}]

Decision letter