Graph Fourier Neural ODEs: Modeling Spatial-temporal Multi-scales in Molecular Dynamics
Authors: Fang Sun
Zijie Huang
Haixin Wang
Huacong Tang
Xiao Luo
Wei Wang
Yizhou Sun
Venue: Accepted by TMLR ()

Abstract: Accurately predicting long-horizon molecular dynamics (MD) trajectories remains a significant challenge, as existing deep learning methods often struggle to retain fidelity over extended simulations. We hypothesize that one key factor limiting accuracy is the difficulty of capturing interactions that span distinct spatial and temporal scales—ranging from high-frequency local vibrations to low-frequency global conformational changes. To address these limitations, we propose **Graph Fourier Neural ODEs (GF-NODE)**, integrating a graph Fourier transform for spatial frequency decomposition with a Neural ODE framework for continuous-time evolution. Specifically, GF-NODE first decomposes molecular configurations into multiple spatial frequency modes using the graph Laplacian, then evolves the frequency components in time via a learnable Neural ODE module that captures both local and global dynamics, and finally reconstructs the updated molecular geometry through an inverse graph Fourier transform. By explicitly modeling high- and low-frequency phenomena in this unified pipeline, GF-NODE more effectively captures long-range correlations and local fluctuations alike. We provide theoretical insight through heat equation analysis on a simplified diffusion model, demonstrating how graph Laplacian eigenvalues can determine temporal dynamics scales, and crucially validate this correspondence through comprehensive empirical analysis on real molecular dynamics trajectories showing quantitative spatial-temporal correlations across diverse molecular systems. Experimental results on challenging MD benchmarks, including MD17 and alanine dipeptide, demonstrate that GF-NODE achieves state-of-the-art accuracy while preserving essential geometrical features over extended simulations. These findings highlight the promise of bridging spectral decomposition with continuous-time modeling to improve the robustness and predictive power of MD simulations. Our implementation is publicly available at https://github.com/FrancoTSolis/GF-NODE-code .

Published in Transactions on Machine Learning Research (06/2025)
Graph Fourier Neural ODEs: Modeling
Spatial-temporal Multi-scales in Molecular Dynamics
Fang Sun
fts@es. uela..llaaaoo edu
University of Colifornia, Los Angeles
Zijie Huang
zijichuang ucla edu
University of California, Los Angeles
Haixin Wang
who Catees..ucla.coo edu
University of California, Los Angeles
Huacong Tang
hetang&ucla. edu
University of California, Los Angeles
Xiao Luo
ziooluoBcs. www.:/essss......oo edu
University of Colifornia, Los Angeles
Wei Wang
weiwangBes. ucla edu
University of California, Los Angeles
Yizhou Sun
yzsum ucla. edu
University of California, Los Angeles
Reviewed on OpenReview: /ttps://openreveeee net/ forum? id= XX7cIdj6Fs
Abstract
Accurately predicting long-horizon molecular dynamics (MD) trajectories remains a significant challenge, existing deep learning methods often struggle to retain fidelity over extended simulations. We hypothesize that key factor limiting accuracy is the difficulty of capturing interactions that span distinct spatial and temporal scales rilosiggnging from high- frequency local vibrations to low-frequency global conformational changes. To address these limitations, we propose Graph Fourier Neural ODEs (GF-NODE), integrating graph Fourier transform for spatial frequency decomposition with a Neural ODE framework for continuous- time evolution. Specifically, GF-NODE first decomposes molecular configurations into multiple spatial frequency modes using the graph Laplacian, then evolves the frequency components in time via learnable Neural ODE module that captures both local and global dynamics, and finally reconstructs the updated molecular geometry through inverse graph Fourier transform. By explicitly modeling high- and low-frequency phenomena in this unified pipeline, GF-NODE more effectively captures long-range correlations and local fluctuations alike. We provide theoretical insight through heat equation analysis on a simplified diffusion model, demonstrating how graph Laplacian eigenvalues can determine temporal dynamics scales, and crucially validate this correspondence through comprehensive empírical analysis on real molecular dynamics trajectories showing quantitative peatial-tempora correlations across diverse molecular systems. Experimental results on challenging MD benchmarks, including MD17 and alanine dipeptide, demonstrate that GF-NODE achieves state-of- the-art accuracy while preserving essential geometrical fcatures over extended simulations. These findings highlight the promise of bridging spectral decomposition with continuous-Lime modeling to improve the robustness and predictive power of MD simulations. Our implementation is publicly available at httpee//gttub..oo//aaaoooooisse NODEE-codd

---

Published in Transactions on Machine Learning Research (06/2025)
1 Introduction
Molecular dynamics (MD) simulations indispensable tools for investigating the behavior of molecular systems at the atomic level, offering profound insights into physical (Bear & Blaisten-Barojas 1998) chemical (Wang et al., 2011), and biological (Salo-Ahen et al. 2020) processes. These simulations must capture interactions occurring across a wide range of spatial and temporal scales- from localized bond vibrations to long-range non-bonded interactions- interactions-poosing significant computational challenges. Accurately modeling these multiscale interactions is crucial for uncovering the mechanisms underlying complex molecular phenomena but remains prohibitively expensive for large systems and long trajectories (Vakis et al.) 2018
In recent years, Graph Neural
Time Scale
ODEs Zang & Wang 2020; Huang
et al. 2023a have gained traction for
modeling continuous-timm dynamics in multi-agent systems, including MD. By learning an Ordinary Differential Equation (ODE) function via Graph Neural Networks and solving it mu- merically, these methods allow flexi- ble sampling at arbitrary time points. Such a continuous-line formulation is well-suited for capturing multiple temporal scales inherent in molecu- lar simulations. However, significant challenges persist in accounting for the rich spatial multiscale effects that span localized bond vibrations to ex- lended nonbonded interactions. On another front, Fourier Neural Oper- ators (FNOs) (Li et al., 2020) have demonstrated success learning op- erators by decomposing signals into different frequency modes, thereby capluring various spatial scales effec- tively. Yet, they are not tailored for graph-structured molccular data or continuous- time temporal evolution, limiting their direct applicability to general MD simulations.
Conformation changes (10xNanometers, Nanoseconds)
Hydrogen bonding and van der Waals forces (Nanometers, Picoseconds)
Bond vibrations and rotations, angle bending (Angstroms, Femtoseconds)
Spatial Scale
Figure 1: Illustration of spatial and temporal multiscale interactions in molecular dynamics. The x-axis indicates spatial scales ranging from local bond vibrations to global conformational changes. The y-axis represents timescales ranging from femtoseconds (bond oscillations) to nanoscconds (conformational rearrangements) Local interactions, such as bond vibrations and angle bending, occur over short spatial scales (angstroms) and fast timescales (femtoseconds). Non-local interactions, including hydrogen bonding and van der Waals forces, span larger spatial scales (nanometers) and intermediate timescales (picoseconds). Conformational changes involve the longest spatial and temporal scales, potentially up to nanoseconds or longer.
Despite the promise of Graph Neural ODEs for handling multiple temporal scales via continuous time modeling, they alone are inadequate for fully capturing the complex spatial frequency components of molecular systems. Conversely, while approaches based on Fourier transforms can model multiple spatial scales, they do not naturally handle the intricacies molecular graphs or continuous time dynamics. To address these limitations, we introduce Graph Fourier Neural ODEs (GF-NODE). As demonstrated Figure our framework explicitly integrates a graph Fourier transform-ffe for dccomposing and encoding spatial multiscale interactions Neural ODE framework for continuous-time modeling of each spatial frequency. By leveraging an inverse graph Fourier transform the end of the pipeline, GF-NODE reconstructs the molecular state in physical space, thereby enabling unified approach to spatial and temporal multiscale simulation.
a
We conduct extensive experiments on benchmark molecular dynamics datasets, including MD17 and alanine dipeptide. Empirical results show that GF-NODE achieves state-of the--rr-tt--ttttt accuracy in predicting molecular trajectorics over long-horizon simulations, preserves essential geometric propertics such as bond lengths and angles, and demonstrates stable performance over temporal super-resolution tasks. These findings underscore the importance of explicitly decomposing molecular configurations into spatial frequency modes and evolving
2

---

Published in Transactions on Machine Learning Research (06/2025)
Graph Fourier Transform
Inverse Graph Fourier Transform
ODEs of Graph Fourier Coefficients
time
GF
9F
M wWA.nwmmnh
g
g
Figure 2: Overview of the Graph Fourier Neural ODE framework. The molecular graph G is first transformed into the spectral domain using a Graph Fourier Transform (GF), dccomposing the spatial structure into frequency components. Neural ODEs are then applied to evolve the Fourier coefficients over time. The evolved coefficients are finally transformed back into the physical domain using an inverse Graph Fourier Transform (GFF"),, reconstructing the molecular graph at future time tg-
them continuously in time. Our analysis suggests that this multiscale perspective is instrumental for capturing both rapid local fluctuations and slow global conformational changes.
Our contribution can be summarized follows.
(a) New perespective. We provide dynamics by jointly capturing spatial and temporal interactions within single framework.
perspective spatial-temporal multiscale modeling for molecular
(b) Novel architeclure. Building on this perspective, we introduce the Graph Fowrier Newral ODEs (GF- NODE) architecture, which combines graph Fourier transform for decomposing molecular interactions into distinct spatial frequencics with a Neural ODE formulation evolve these frequencics continuously in time.
(e) Theorelical foundations. We establish theoretical insight into the correspondence between spatial frequency modes and temporal dynamics scales through heat equation analysis simplified diffusion model, showing that graph Laplacian cigenvalues can intrinsically determine temporal evolution rates. Importantly, we provide comprehensive empirical validation through our Theoretical Framework for Joint Analysis, demonstrating quantitative spatial-temporal correlations on real molecular dynamics trajectories across diverse systems.
(d) Comprehensive evaluation. We conduct extensive experiments across multiple benchmarks, including the Revised MD17 dataset, larger molecular systems (20-326 heavy atoms), and extended temporal horizons up to 15,000 MD steps, consistently achieving state-of-the-art performance with detailed structural analysis confirming accurate preservation of molecular correlations.
2 Related Work
We review relevant works on multi-scale modeling in molecular dynamics, focusing on neural operator models and graph neural ODEs.
2.1 Classical Molecular Simulation Methods
Traditional molccular simulation methods, including force-field based MD (c.g., AMBER (Cornell al. 1995 Brooks et al. 2209),, GAMD (Li et al., (202bb)) and ab initio techniques such as Parrinello MD (Hutter 2012), have been foundational in exploring molecular behavior. However, these methods face significant.
3

---

Published in Transactions on Machine Learning Research (06/2025)
limitations: force-field simulations require extremely small time steps to accurately resolve high-frequency bond vibrations, which hampers long-term stability and computational efficiency; ab inilio MD, though offering first-principles accuracy, is computationally prohibitive for large systems and long trajectories; and while coarso-grained models (e.g., MARTINI (Souza et al. 2021)) enable more efficiont multiscale simulations, they often compromise on molecular detail and accuracy, particularly in reproducing local interactions and maintaining seamless force consistency at multiscale interfaces.
2.2 Neural Operator Models
Neural operator models (Kovachki et al. 2023) have emerged as powerful lools for learning mappings between infiite-dimenniional function spaces, demonstrating success in modeling complex dynamical systems. Among these, Fourier Neural Operators (FNOs) (Li et al. 2020 2022a Liu & Jafarzadeh 2023 Kovachki et al. 2021 Koshizuka et al.) are particularly notable for handling spatial multi-scale interactions in partial differential equation (PDE) data by learning representations in the Fourier domain. However, while FNOs efficiently capture spatial hierarchies, they do not inherently model temporal dynamics, making them suboptimal for time-evolving molecular systems.
In contrast, recent operator-based methods such as the Implicit Transfer Operator (TTO) (Schreiner al., 2024), Timewarp [Klein et al., 2024). and Equivariant Graph Neural Operator (EGNO) [u et al 2024) focus temporal multi-scale modeling in molecular dynamiis. ITO and Timewarp introduce coarse graining and adaptive time-stepping mechanisms to accelerate cong-oorizon simulations. EGNO employs neural operators with SE(3) equivariance to capture rotational and translational symmetries, yet it primarily addresses temporal evolution without explicitly handling spatial multi-scale effects. While these models successfully extend the applicability of neural operators to molecular simulations, none jointly addresses both spatial and temporal multi-scales in molecular dynamics.
2.3
Graph Neural ODE Models
Graph Neural ODEs combine Graph Neural Networks with Neural ODE frameworks (Chen et al. 2018
Kidger 2022
Goyal & Benner 2023 Holt et al. 2022 Luo et al. 2023 Luo al. to model continuous time
dynamics on graph-structured data Zang & Wang 2020 Kim et al. 2021 Huang et al. 2020 2021 2023a b).
These methods excel in capturing lemporal multi-scale behavior by a allowing flexible time integration, making them well-suited for systems with varying temporal resolutions. However, they primarily focus on modeling temporal dependencies, with little emphasis on explicitly handling spatial multi-scales in molecular dynamics.
A key limitation of existing Graph Neural ODE models is their reliance on local message passing, which inherently constrains their ability to capture long-range spatial dependencies within molecular systems. As result, they may fail to adequately represent the interplay between localized, high-frequency interactions (e.g., bond vibrations) and global, low-frequency effects (e.g., large-scale conformational changes). Unlike spectral-based approaches that can decompose spatial hierarchies, standard Graph Neural ODFs lack a mochanism to explicitly encode spatial multi-scale structures, limiting their effectiveness in modeling complex molecular dynamics.
3 The Proposed Approach
We propose GF-NODE, a framework specifically designed to address the limitations of existing methods in capturing both spatial and temporal multiscale dynamics in molecular systems. As discussed in Sections and 2 current approaches either handle spatial scales using Fourier-based methods focus on temporal scales using Graph Neural ODEs, but they do not jointly model these scales in unified framework. GF- NODE directly addresses this gap by integrating the Graph Fourier Transform (GFT) and Newral Ordinary Differential Equations (Neural ODFs) to simultaneously decompose spatial interactions and model their continuous-time evolution.
Specifically, as illustrated in Figure 3. GF-NODE first applies a Graph Fourier Transform to decompose the molecular graph into different frequency components, effectively separating localized, high-frequency
4

---

Published in Transactions on Machine Learning Research (06/2025)
ODE Solvers
time
tx
Spectral Domain
GF
GF
GF
GF
Encoder
Decoder
Decoder
Decoder
g
g
g
Physical Domain
Figure 3: An overview of the proposed GF-NODE architecture. The model first encodes the initial molecular graph G from the physical domain into the spectral domain via a graph Fourier transform (GF). Neural ODE solvers then propagate the dynamics the Fourier coefficients continuously across time points to, f,...., .k. The transformed coefficients are subsoquently decoded back to the physical domain using an inverse graph Fourier transform (GF-"), reconstructing the molecular graphs G g ,,.,,,, g al future time points. This design enables efficient modeling of spatial and temporal multiscale dynamics in molecular systems.
interactions from global, low-frequency patterns. This spectral representation allows the model to process molecular structures in frequency-sdaalive manner, capturing both fine-grained local interactions and large-scale conformational changes. The decomposed spectral coefficienss are then evolved continuously over time using Neural ODEs, ensuring flexible and adaptive modeling of multiscale temporal dynamics. Finally, the inverse Graph Fourier Transform reconstruets the molecular graph in the physical domain, preserving both local and global structural information over long-horizon molecular simulations.
The theoretical motivation for this design stems from heat equation analysis on simplified diffusion models, which suggests that spatial frequency modes (characterized by graph Laplacian eigenvalues) naturally correspond to different temporal evolution scales. Importantly, we provide comprehensive empirical validation of this correspondence through detailed spatial-temporal correlation analysis real molecular dynamics trajectories, confirming that the spectral decomposition elfectively separates distinct temporal dynamics scales in practice.
This design enables GF-NODE to overcome key challenges in molecular dynamics modeling: (1) explicitly encoding spatial multi-scale interactions via spectral decompasition, (2) leveraging continuous-Lime evolution capture complex temporal dependencies, and (3) integrating these spatial and temporal scales into single end-to-cnd framework that maintains SE(3) equivariance (formal proof provided in Appendix Dj. Crucially, approach is grounded in rigorous theoretical foundations: graph Laplacian eigenvalues intrinsically determine the characteristic time-scales of dynamics through heat equation analysis, establishing that low-frequency spatial modes naturally correspond to slow temporal evolution while high-frequency modes correspond to rapid dynamics (detailed analysis in Section 4.4 and Appendix). Below, we detail the core components and operations of GF-NODE.
3.1 Notation and problem setup
Let (V,8) E) represent a molccular graph, where N the number of nodes E (1,..., N) (atoms), and & is the of edges representing chemical bonds or interactions. Edges are identified by checking whether the distance between atoms falls below threshold. Each node has:
5

---

Published in Transactions on Machine Learning Research (06/2025)
Invariant (scalar) features h, R,, which include the velocity magnitude |vvll and normalized version of the atomic number Z;; Vector features z, E Rmxx,, which contain position x, and velocity v, (i.e. z, = (x,, v:) when m = 2). At time to, the molecular system's state is given as (h, z), and we aim to predict the future configuration g(tx) for any 1g lo.
3.2 Graph Neural Network encoder
We use Graph Neural Network (GNN) to encode scalar features h, E R and vector features z, E Rm for each node The initial scalar feature h for each atom is formed by concatenating ||v and Z where Z is the atomic number of atom i, and Zumx is a reference maximum. This concatenated vector is then mapped by linear embedding layer, producing the hidden dimension used by the GNN. The vector feature (0) contains x, and Va. Each GNN layer performs message passing, where node i's features are updated based its neighbors N():
z
(1) mij
= ge
h .
(1)
(1)
h
((+1)
Φh
h
(1) Σ m (1)
(2)
JEN(()
((+1)
x;
(1)
+
Σ
V m
(3)
[iii)
(EN(()
where K = =x - x and de, On. and & are neural networks. Similarly, we also update the velocity v. if present. Afler L layers, the GNN produces the encoded features h and (4)) capturing both local and global molecular information. These features are then passed to the Graph Fourier Transform (GFT) for spectral decomposition.
3.3 Graph Fourier Transform
Molecular dynamics exhibit behavior on multiple spatial scales: large-scale "global" deformations can be viewed as low-frequency modes, whereas fast "local" vibrations correspond to high frequency modes. In classical signal processing, Fourier analysis decomposes signals into sinusoids of different frequencies. Similarly, on graphs, we can decompose nodo-based signals into cigenmodes of a suitable operator (often the graph Laplacian). By retaining emphasizing certain frequency bands, one can explicitly model global vs. local patterns.
Graph Signal as Combination of Laplacian Eigenvectors. A scalar function V + R on the nodes can be seen as vector f E RN, The graph Laplacian L = D A (or a symmetrized variant) admits an eigen-decomposition::
L- UAU
(4)
where U = [uo, un,..., 1N is an orthonormal matrix of eigenvectors, i.e. us E RN,, and A = diag (Ao, , contains ascending eigenvalues Ao S A1 S ... AN-1-
In graph signal processing, the eigenvector uN is viewed as the "k-th frequency basis." Smaller eigenvalues ( (Aa,A1,.... correspond to low-frequency (more global, smooth) variations the graph, while larger eigenvalues correspond to high-frequency (local, rapidly changing) modes. Hence, any signal f e RN can be writlen as linear combination:
Σ
agus:
where on = u f.
(5)
This collection of coefficients {ax} is the Graph Fourier Transform (GFT) of f. We now denote og as Je-

---

Published in Transactions on Machine Learning Research (06/2025)
Truncation at M Bases. For large N, we often use only the first M eigenvectors {uo..... .u...... UM-1) to approximate f. This yields band-limited or iultiscale representation:
f
Og Uk.
(6)
where M S N. Sorting A in ascending order ensures the lowest-frequeney modes (i.c., global scales) appear first, and higuer-frequency (local scales) modes appear last. By choosing M appropriately, we focus on the most critical modes for modeling global vs. local behavior.
Applying the GFT Scalar vs. Vector Features. For scalor features H NxF, where F is the latent feature size, we substitute fo for Hg in the equation fx =uf (from Equation get:
M-1 H
k=0 = U. H,
(s.0:M)
where M) denotes the first M eigenvectors. For vector features Z E RVxxx3 where m is the feature size of the vector feature in the 3D space, we first remove the mean to ensure translational invariance (since the 0-th eigenvector un corresponds to the constant mode):
Z. Z-Z
where Z is the global mean over all nodes. We then apply the same truncated basis (..:.)) to each coordinate dimension, substituting f for Zg in the equation fi- = u (from Equation to get:
Ż = UT UT... (s.D:M) Z. E RRRM Mxmx3
(9)
Indices closer to = 0 indicate global motions, while larger i (up to M - 1) captures more local, high-froquency fluctuations. In practice, M can be hyperparameter that determines how many eigenmodes we keep, balancing elliciency (fewer modes to evolve) and accuracy (how many scales are caplured).
Theoretical Analysis. The Graph Fourier Transform (GFT) decomposes molecular dynamics into modes corresponding to global and local spatial patterns, as determined by the eigenvalues and cigenvectors of the graph Laplacian. As demonstrated in Proposition 3.1 low-frequency modes capture smooth, global deformations, while high-freeuency modes represent localized structural variations. This decomposition enables efficient modeling of multiscale spatial dynamics.
Proposition 3.1 (Global vs. Local Spatial Scales). Let ug be the k-th eigenvector of with eigenvalue dg- Suppose x encodes alomic coordinates or their latent features. Then:
1. If > is small, the corresponding mode un represents slowly varying (global) deformations across the molecule.
2. If \₂ is large, the corvesponding mode represents rapidly changing (local) structural variations.
Regarding the number of modes M need to use, Theorem 3.2 guarantees that use the first few modes to represent the entire system accurately:
Theorem 3.2 (Spectral Truncation Error). Truncating the spectral representation to the first M modes,
X(M) (:...M) U. ((.,,),,, yields an e-norm error:
N-1
x|- x(M) 12==
Σ K=M
u x
(10)
For a-bandlimited signals, choosing M such that AM-1 S guarantees exact recovery.
We provide detailed proofs and derivations in Appendix B where we analyze the spectral decomposition's efficacy and demonstrate its suitability for multiscale molecular modeling.
7

---

Published in Transactions on Machine Learning Research (06/2025)
3.4 Neural ODEs in the Spectral Domain
We propose novel approach that models the temporal evolution of the Fourier coefficients using Neural ODEs. By representing molecular interactions in the frequency domain, this approach enables the decomposition of dynamics across different spatial scales and provides a more compact representation for modeling temporal evolution. Neural ODEs then learn the dynamics of these Fourier coefficients over time, offering a continuous-time framework that captures multiscale spatial and temporal interactions simultaneously.
NeuralODE Preliminaries. Neural Ordinary Differential Equations (Neural ODEs) (Chen et al 2018
provide a framework for modeling continuous-time dynamics by learning the evolution of system as set of differential equations parameterized by a neural network. Specifically, for a system's 's state h(t), its temporal evolution is defined as:
dh(t) = fo th()).t)).
(11)
dt
where fo is a neural network parameterized by that learns the dynamics of h(t) over time. Given an initial state h(ta), the state at any future time t is computed by solving the ODE:
h(t) = h(6)) + S(d [o(((T)),)dTT.
(12)
This integral can be evaluated numerically using adaptive ODE solvers, such Runge-Kutta methods, allowing Neural ODFs to handle irregularly sampled data and continuously model temporal dynamics.
Joint Scalar-Vector Block-Diagonal Formulation. We update the spectral coefficients [H, 2] jointly.
We can write the combined spectral features as
NXF
+ R Nxmx3
(13)
and define Neural ODE dynamics over continuous time t:
H(t) 7(t)
((o(H()) (g)()
(14)
Mathematically, in the Fourier space, this equates to block-diagonal operator Mg. Let H. and Z. denote the coefficients for each frequency mode w. A single ODE step can be represented as:
H %
M (h)
O
H. Z...
(15)
M
Here, M(()) operates on the scalar channels and M(2) on the vector channels, allowing them to evolve in a coordinated but distinct manner. This is beneficial to preserving the 3D interactions while keeping track of the nuanced representations the latent space thus preserving the 3D interactions while capturing the muanced latent representations in block-diagonal design.
Implementation wisc, each ODE function fo go takes the current spectral modes (H Zw) and the time I, and produces their rate of change for the inlegrator. Specifically, we use multi-head self-atlention across these modes so that each frequency mode can attend to others. Formally, if (t) denotes the coefficients corresponding to frequency mode w at time f, the self-attention step can be written as
h (t) MHAttm (t)).
(16)
where MHAttn is the multi-head attention layer across different frequency modes. Next, the time & directly concatenated along the feature dimension via learned embedding >(t) Rd, yielding
h (1)
(17)
8

---

Published in Transactions on Machine Learning Research (06/2025)
Finally, each mode (t) is transformed by mode-wist lincar weight Wa. resulting in the derivative:
fo (t), = WW (t).
(18)
A similar procedure applies for go, handling any additional vector channels by suitably reshaping and performing attention across modes.
To predict time tk. we numerically integrate [.. Z] from to to is (e.g., via dopri5). Low-frequency components capture long-range, slower dynamics, while high-frequency components capture faster local fluctuations.
Theoretical connection between spatial and temporal modes. We also use simplified heat equation dynamics model to demonstrate the connection of spatial GFT modes to distinct temporal dynamics. Proposition 3.3 (Heat--Equatio Mode Dynamics). Let G = (V, E) be a graph with normalized Laplacian L admitting
L-UAU, A diag(Ao, .... -)); 0= Ao S .... AN--- Consider the graph-heat evolution scalar signal ((t N
f(((t) - L
(19)
dl
Write the kth Graph Fourier coefficient atl)) uh f(t). Then each mode evolves independently:
dax(t) du at())
(20)
dt
with closed-form solution on = (0) exp(- t). In particulor, modes with X small decay slowly (long time-scales), while modes with A large decay rapidly (short time-scales).
Although our GF-NODE dynamics are learned rather than the pure heat equation, Proposition 3.3 shows that under any diffusion-like operator the Laplacian cigenvalues X set intrinsic time-scales for each mode. The detailed proof and theoretical analysis can be found in Appendix C
3.5 From Spectral Domain to Physical Domain
Inverse GFT. After evolving the spectral coefficients (1x) and 2 (4x), we recover the node-level signals via the inverse GFT:
H (fx) = U (.. H (4x), Ze (fx) (G.0 (4x),
(21)
where U(,u:ar) (the matrix of the first M eigenvectors) is the same truncated basis used during the forward GFT. Finally, we restore the global translation by adding back the mean Z that was subtracted earlier:
Z (fx) Ze (tx) Z.
(22)
Graph Neural Network Decoder. We refine local interactions at each predicted time & with GNN decoder that again operates on [H()), Z(4)]_ This step can help capture short-range correlations that may not be fully resolved in the spectral update. The decoder GNN has a structure similar to the encoder:
m d' (hi(le), hū(lk), rij
(dec)
(dec)
(1x)
(23)
h
(dec) (dec)(tx) (tx) = 0% h((t)).
Σ JEN(()
m
(dec)
(24)
x (de))((+) = x1444) + [N(1]
(dec) (1x)
Σ
(mm (dec)
(25)
JEN())
and similarly for v,. Here, ri
x;(1x) (1x) (dec)(1x))
(dec)
(1(,,,, and &6 %ho 18 are learnable functions.
9

---

Published in Transactions on Machine Learning Research (06/2025)
4 Experiments
In this section, we provide a comprehensive evaluation of our proposed model on molecular dynamics datasets, comparing its performance against several basclines and conducting detailed ablation studies to assess the contributions of the different components. Our experimental evaluation designed to address the following key research questions:
RQ1: Prediction Accuracy. Does the proposed GF-NODE framework deliver improved molecular dynamics prediction accuracy compared to state-of-the-art methods?
RQ2: Continuous Time Modeling. How effectively does the continuous-tine time evolution component
capture multiscale temporal dynamics including long- Long-horizon forecasting and super-resolution compared to variants without continuous-time time propagation?
RQ3: Spatial Multiscale Modeling. How crucial the explicit spectral decomposition for capturing spatial multiscale interactions, and how do different Fourier mode interaction schemes and the number of retained modes affect overall performance?
RQ4: Theoretical Justification. Can rigorously establish the correspondence between spatial frequency modes (graph Laplacian cigenvectors) and temporal dynamics scales, providing thcoretical loundation for the GF-NODE architecture?
4.1 Dataset
We evaluate our model using both the Revised MD17 dataset (Christensen & von Lilienfeld 2020 and the original MD17 dataset (Chmiela et al. 2017), which contain molecular dynamics trajectories for small molecules including Aspirin, Benzene, Ethanol, Malonaldehyde, Naphthalene, Salicylic Acid, Toluene, and Uracil. The datasets provide atomic trajectories simulated under quantum mechanical forces, capturing realistic molecular motions. To demonstrate scalability, we further evaluate on five larger molecular systems with 20-326 heavy atoms: alanine dipeptide (Alaz) from MDShare, Ac-Alay-NHMe, AT-AT-CG-GG, Bucky- Catcher, and a doulle-walled carbon nanotube, from the MD22 dataset (Chmiela et al.) 2023 These larger systems exhibit complex multiscale dynamics spanning local bond vibrations to global conformational changes.
We also evaluated our model on the alaninc dipeptidc dataset (Shhriinrr et al. 2024 a standard benchmark for studying conformational dynamics in proteins. Our task is to predict the future positions of atoms given the initial state of the molecular system.
4.2 Experimental Setup
For each molecule, partition the trajectory data into training, validation, and test sets, using 500 samples for training, 2000 for validation, and 2000 for testing. The time scope AT for each piece of data is set to 3000 simulation steps, providing a challenging prediction horizon. To demonstrate long-term stability, we extend our evaluation to At = 10,000 steps and up to 15,000 MD steps for comprehensive temporal analysis.
A key aspect of our experimental sctup the use of irrogular timestep sampling, in contrast to the equi-timestep sampling used in some baseline models like EGNO, to better mimic the variable time intervals in weal-world physical systems. This setting tests the models' ability to handle irregular temporal data. Although each trajectory spans 3000 simulation steps, we randomly sample only 8 data points per instance. Because these samples are drawn from different points in the 3000-step window across data instances, this strategy enables the model to learn the dynamics over the entire time span without requiring training on every timestep, thereby significantly enhancing efficiency. Nevertheless, we also provide evaluations based on tquistim sampling in Appendix E for completeness.
4.3 Baseline Models
We compare our model against several state-of-the-art approaches:
NDCN ZZagg& & Wang 2020): Graph Neural ODE model that integrates graph neural networks into
the ODE framework to learn contimuous-time dynamics of networked systems.
LG-ODE (Huang al. 2020): A latent graph-based ODE model that integrates latent representations into continuous-tine evolution.
10

---

Published in Transactions on Machine Learning Research (06/2025)
EGNN (Satorras et al. 2021]: An Equivariant Graph Neural Network that models molecular systems
using 3D equivariant message passing but without explicit time propagation.
EGNO Xu et al. 2024] An Equivariant Graph Neural Operator that captures temporal dynamics using neural operators with regular timesteps.
ITO Schreiner al. 2024) An Implicit Time-stepping Operator that integrates differential equations into the learning process for temporal evolution.
These baselines represent range of approaches for modeling molecular dynamics, including methods that focus primarily spatial modeling (EGNN) temporal modeling (NDCN, LG-ODE, EGNO, ITO).
4.4 Results and Analysis (RQ1)
Units and Calculations: Note that the alanine dipeptide dataset operates nanometers (nm), whereas the MD17 dataset uses augstroms (A). When calculating bond lengths and bond angles, we consider all heavy aloms (excluding hydrogen) to focus on the core structure.
Our Prediction
Ground Truth
343.0 ps
867.0 ps
1324.0 ps
2735.0 ps
Figure 4: Representative snapshots of alanine dipeptide, comparing our model's predicted conformations (top row) against the ground-truth simulation (bottom row) at four different limestamps. The close agreement. illustrates the model's ability preserve key structural features over an extended trajectory.
To address RQ1, we evaluate the performance of our model and the basclines on both the Revised MD17 dataset and larger molecular systems, with Figure 4 providing visualization of our predictions. Our comprehensive evaluation includes results at multiple time horizons (At = 3000 and 10,000 steps) and extended evaluations up to 15,000 MD steps. The test Mean Squared Error (MSE) for our model and the baseline methods are summarized Tables and 3 for the original evaluation, under irregular timestep sampling. Comprehensive results on the Revised MD17 dataset and larger molecular systems
be found in Tables 11 14 in the appendix, demonstrating consistent improvements across all systems and time horizons. For completeness, we also provide results in the equi-timestep setting in Table 7 in Appendix E
Table MSE (x10-2 2 A2) on the MD17 dataset with irregular timestep sampling, Best results are in bold, and second-ber best are underlined.
Aspirin
Benzene
Ethanol
Malonaldehyde
Naphthalene
Salicylic
Toluene
Uracil
NDCN
29.75=0.02
70.13+0.98
10.05=0.02
12.28+0.07
2.30=0.00
3.43++0.05
12.33-0.00
2.39:+0.00
LG-ODE
51.65+0.01
68.29-0.21
12.32-0.05
43.95++0.07
2.38-0.02
2.85++0.08
18.11-0.09
2.38+0.077
EGNN
9.09-10.10
19.15.11.68
4.46:10.01
12.52.10.05
0.40=0.02
0.89
8.98-0.09
0.64 0.00
EGNO
10.60-0.01
52.53 40
4.52-10.06
12.80-0.006
0.46-0.01
1.07-10.00
9.31100100
0.67-10.01
ITO
12.74=0.10
07.8410.86
7.23:10.00
19.53.1.001
1.77-0.01
2.53:10.03
9.9641.04
1.7140.11
Ours
6.46++0.03
1.52+0.08
2.74+00.05
10.54+0.01
0.23+0.02
0.63+0.01
1.80+0005
0.41++0.1
11

---

Published in Transactions on Machine Learning Research (06/2025)
Table 2: Mean Absolute Errors (MAEs) and relative errors for bond lengths and bond angles on alanine dipeptide. Best results are bold.
Model
Bond Length MAE (nm)
Rel. Err. (%)
Bond Angle MAE ())
Rel. Err. (%)
EGNN 0.0209 + 0.0006
15.32 + 0.49
12.44 + 0.91
10.48 + 0.76
EGNO
0.0229 = 0.0018
16.75 + 1.23
10.54 + 0.11
8.89 = 0.11
Ours
0.0188 + 0.0022
13.74 74 + 1.66
10.47 + 1.03
8.80 = 0.89
From Table [ we observe that our model consistently outper- forms the bascline methods across all eight molecules under irregular timestep sampling. This demonstrates the effective- ness of our approach in jointly modeling spatial and temporal multiscale interactions. The performance gains are particu- larly pronounced for molecules such as Benzene and Aspirin, where our model significantly reduces the MSE compared to the baselines. Our comprehensive evaluation on the Revised MD17 dataset and five larger molecular systems (20-326 heavy atoms) shows that GF-NODE achieves the best accuracy on all cases across multiple time horizons (At = 3000 and 10,000 steps), with detailed results provided in the appendix (Tables
Table 3: Mean Squared Error (MSE) (x10-3 nm²) on the alanine dipeptide dataset. Best results are in bold.
Model
MSE (x10 mm?)
NDCN
12.27 + 0.19
ITO
26.95 + 0.19
EGNO
6.92 + 0.26
EGNN
5.67 = 0.08
Ours
4.48 0.07
14p. The model demonstrates superior long-horizon stability up to 15,000 MD steps (Figures 10 and [1]. capturing both short-range covalent vibrations and long-range collective modes such as radial breathing in carbon nanotubes.
Benzene Drifting. interestingly, we find that Benzene exhibits substantial drift during simulation. In our data, the maximum g-coordinate is 197.981A A, while the minimum z-coordinate is 1178.112 A, indicating large translations and rotations of the entire molccule. In contrast, other molecules exhibit minimal net translation, typically remaining within +3A from the origin. Methods that enforce strict invariance to translations and rotations (e.g., EGNN) may underperform on such drifting systems, since the global drift is part of the actual dynamics. Indeed, we discovered that replacing EGGNN-based layers with standard message passing layers (e.g., SAGEConv) can further boost performance on drifting molecules like Benzene. We provide details of the dataset isatis in Table 5 of Appendix A and experiment results using different GNN architectures in
Table g of Appendix
Experimental Results on Alanine Dipeptide. In addition, from Table 3) our model achieves the lowest MSE among all com- pared methods. This indicates a significant improvement on more complex molecular dynamics data and demonstrates the robustness of our approach.
70
Analysis of Molccular Structure Recovery. To further evalu- ale predictions al a structural level, we analyzed the bond lengths and bond angles for alanine dipeptide. Table 2 shows the Mean A Absolute Errors (MAEs) and relative errors for bond lengths and angles. Our method achieves the lowest errors, indicating superior recovery of internal molecular structures compared to baselines. Additionally, extensive radial distribution function (RDF) analysis for five representative systems (Figures 17 in Appendix F) con- firms that GF-NODE accurately reproduces both first-shell peaks and longer-range oscillations of the ab initio reference, demon- strating that predicted trajectories maintain realistic structural correlations well beyond pairwise coordinate accuracy.
60
10
Number Time-laged Companents
Figure 5: TICA scores for varying num-- bers of components. Lower scores indi- cale better alignment with slow modes.
TICA Analysis. Time-lagged Independent Component Analysis ( & Schuster 1994] is
used to extract slow collective motions from the trajectories. Figure 5 shows that, across various numbers
12

---

Published in Transactions on Machine Learning Research (06/2025)
of TICA components, model consistently achieves lower TICA scores (i.e., better alignment with the underlying slow modes) compared to EGNN and EGNO, indicating more effective capture of the underlying multiscale dynamics.
4.5
Continuous-Time Dynamics and Temporal Super-Resolution (RQ2)
A Ablaiion on continuous time components. To assess our model's ability to capture temporal multiscale
dynamics, we perform an ablation study on the continuousttime components. Specifically, we investigate:
1. No ODEs Evolution:
Removing the continuous-tinm evolution altogether.
2. No ODEs on Scalar Channels: Freeze the ODEs for among scalar features h during time propagation.
3. No ODEs Vector Channels: Freeze the ODEs among vector features x during time propagation.
Table 4 reports the test errors (MSE) for these variants across several molecules. The inferior performance of the ablated models confirms the importance of the contimuuus-inm dynamics and the effectiveness of the blok-diiagonal architecture that jointly propagates scalar and vector features.
Table 4: Ablation study on continuous-time evolution components (MSE x10-2 A values). Lower values indicate better performance. All results are inferior to our standard model (w/ ode).
Aspirin
Benzene
Ethanol
Malonaldehyde
Naphthalene
Salicylic
Toluene
Uracil
no ode 6.5640.03
18740008
3.09:1005
11.5810.02
0.4240.02
0.8620.01
3.0620.05
0.61=0.01
ode h 6.7740.05
1.994000
3.11+0003
10.70=10.0
0.4240.01
0.8820.02
3.48-0004
0.59=0.02
ode
6.95
1.7940.07
3.74220.05
10.61:0.03
0.4240.01
0.87
2.8920.01
0.59=0.01
w// ode 6.46=00.03
1.52=0.08
2.74+0.00
10.54+0.01
0.23+0.02
0.63+0.01
1.80+0.05
0.41#00011
Super-Resolution Task. To further validate temporal generalization, we perform a super-resolution experiment on the alanine dipeptide dataset by predicting trajectories at a 10x finer temporal resolution than the training samples, under the equi-timestep setting. Figure 0 compares the MSE for the original versus the super-resolved predictions. Our model maintains low error under super-resolution, demonstrating its ability interpolate continuous dynamics effectively.
Long-Term Prediction Stability. We evaluate the stability of long-term predictions across extended time horizons. We test models at 1000, 2000, 3000, 4000, and 5000 simulation steps on two representative molccules (Benzene and Malonaldehyde). Figure illustrates how errors evolve over these extended horizons. Our model maintains superior performance with slower error growth, indicating better capture of global low-frequency dynamics essential for accurate long-term predictions. In contrast, baselines lacking explicit multiscale modeling accumulate errors more rapidly. Extended evaluations up to 15,000 MD steps on both the Revised MD17 dataset (Figure 10) and larger molecular systems (Figure 11 further demonstrate the superior long-horizon stability of our approach, with computational efficiency analysis provided in Figure T2
Table in Appendix E provides additional ablations on the different types of temporal embeddings used in the ODE functions. A simple concatenation of the timestamp would work well enough.
min
1250
1207
103
2003
1ea
Figure 6: Comparison of MSE on the alanine dipep- tide dataset: Original vs. 10x Super-resolution.
Figure 7: Error growth over long-term forecasts (up to 5000 steps) for Benzene and Malonaldehyde. Mod- els with too high an error not presented here.
13

---

Published in Transactions on Machine Learning Research (06/2025)
4.6 Spatial Frequency Decomposition and Multiscale Analysis (RQ3)
To examine the role of spatial multiscale modeling, we perform ablations the spectral decomposition and mode interaction components. We consider two sets of modifications:
1. Ablation of the Spatial Decomposition:
No Fourier-based Decomposition: The model is run without any frequency-based transformation. The spectral decomposition is replaced by a multi-layer perception. Replacing GFT with FFT: The Graph Fourier Transform is substituted with a standard Fast Fourier Transform, disregarding the graph structure.
2. Fourier Mode Interaction Schemes: We compare different strategies for inter-mode communication:
Altention-bused Interaction: Each Fourier mode interacts with others via a multi-head self-attention mechanism (Used in our standard model).
Concatenation- based Interaction: Modes are concatenated before being processed.
No Interaction: Each mode is propagated independently.
Table 5 compiles the results for these spatial ablations. We observe that using the GFT for spectral decomposition rather than an MLP or FFT is most effective, underscoring the importance of capturing the inherent graph structure. Morcover, interaction schemes (whether via attention or concatenation) improve performance over treating modes independently.
Table 5: MSE (x10-2 2 A3) on the MD17 dataset with irregular timestep sampling for variants in spatial decomposition and Fourier modes interactions. Best results in bold, the standard model (GFT & Attn.).
Aspirin
Benzene
Ethanol
Malonaldehyde
Naphthaleme
Salicylic
Toluene
Uravil
FFT
6.53+0.03
1.94+0.08
3.36++0.05
10.83+0.01
0.43+0.02
0.88-0.01
3.36+0.05
0.60=0.01
No Fourier
6.6310.04
1.900.10
3.2810.06
10.68=0.07
0.4310002
0.88_0.01
3.7340.09
0.6020.05
No Interaction
6.51+0.08
1.7840.08
3,28:0.05
10.64=0.01
0.41:0002
0.87=0.02
3.08-0.05
0.59 0.00
Coucal.
6.55:10.03
1.7610007
3.13:10.05
10.58:10.01
0.43::0.02
0.88=0.01
2.80:10.05
0.58
GFT & Attn.
6.46++0.03
1.52+0.08
2.74+0.05
10.54+0.01
0.23+0.02
0.63++0.01
1.80+00005
0.41+0.01
4.70
NiF
Moan Sqqurved Emmor preic suar mm)
Impact of the Number of Fourier Modes. Finally, we investigate how the number of retained Fourier modes Eogureeploos the MSE against different
3.00
affects performance. numbers of modes used in the spectral decomposition. We observe that the performance improves as more modes are included up to a threshold, beyond which additional modes yield diminishing returns. This behavior is consistent with
4:00
4.4
12
20
the theoretical analysis presented earlier in Theorem 3.2 The number of modes used get the optimal results for each type of molecule can be found in Table of Appendix A
Figure 8: Effect of the number of Fourier modes on prediction performance. Perfor- mance platcans beyond n = 12.
4.7 Theoretical Justification and Opatial-Tempont Correspondence (RQ4)
To address concerns about the conceptual rigor of linking spatial GFT modes to distinct temporal dynamics, we provide both theoretical insight and comprehensive empirical validation. Our analysis demonstrates that. spatial frequency decomposition naturally corresponds to different temporal evolution scales in molecular systems.
Theoretical Insight through Simplified Model.
We establish theoretical foundation through heat equation analysis on a simplified diffusion model (detailed in Appendix CC This analysis shows that on a continuous domain, the graph Laplacian cigenvalucs 2₂ intrinsically determine the temporal evolution rates of different spatial frequency modes, with higher eigenvalues corresponding to faster decay rates. While this provides valuable theoretical insight, it represents simplified model of molecular dynamics.
14

---

Published in Transactions on Machine Learning Research (06/2025)
Spatial-emmprral Correation: 0.424
buckeballi-catcher
Sllemporl Correlation revisedd_elccylle
Correlation: 0.816
(Tempore Frequeno
alo (temppoed Prequercy
1.3
Aug (Spadal Frequency)
10ggsAAl (Spatial Frequency)
Spempollll Correlation Analysis revised nephttaleeee
Spetiel-emporal Correlation revised__ras
Correlation: 0.622
Correlation: 0,716
do (temporaal Frequarcyl
Wipul ctempore
(0g (Sp30a) Frequency)
M2gMMtrSppata Frequency)
Figure 9: Spatial-temporal correlation analysis across four molecular systems: (a) Buckybaal-catcher, r 0.424, a 0.552: (b) Salicylic acid, 0.816, a 2.642; (c) Naphthalene, 0.622, a 5.016; (d) Uracil, 0.716, 4.098. Each panel shows log-log plot of characteristic temporal frequency fx versus spatial frequency A (graph Laplacian eigenvalues) for individual cigenmodes. The consistently positive correlations validate the predicted correspondence between spatial and temporal scales across diverse molecular architectures.
Empirical Validation Real Molecular Systems. To validate this correspondence actual molecular trajectories, we developed comprehensive Theoretical Framework for Joint Spatial-Temmpral Analysis in Appendix Section C.Iphat quantifies spatial-temporal correlations across diverse molecular systems. Our empirical results on five representative systems (Naphthalene, Salicylic Acid, Uracil, Bucky-Catcher, and couble-walles nanotube), as demonstrated in Figure shows clear quantitative relationships between spatial frequency modes and temporal dynamics scales.
Key findings include: (1) Distinct spatial modes exhibit different temporal autocorrelation decay rates, with low-frequency modes showing slower decay dorg-term dynamics) and high-frequee feth-frequency modes showing faster decay (short-term fluctuations).
(2) Cross-correlation analysis reveals that spatial modes with similar frequencies exhibit. stronger temporal correlations. (3) The empirical correspondence between spatial eigenvalues and temporal scales aligns with theoretical predictions from the simplified heat equation model.
These results provide strong empirical support for our core hypothesis that spatial spectral docomposition effectively separates distinct temporal dynamics scales, justifying the GF-NODE architecture's design principle of processing different. frequency modes with tailored temporal evolution pathways.
15

---

Published in Transactions on Machine Learning Research (06/2025)
5 Conclusion
In this work, we presented Graph Fourier Neural ODEs (GF-NODE), a novel framework that unifics spatial spectral decomposition with continuous-lime evolution to effectively model the multiscale dynamics inherent. in molecular systems. By decomposing molecular graphs via the Graph Fourier Transform, our approach explicitly disentangles global conformational changes from local vibrational modes, and the subsequent Neural ODE-based propagation enables flexible, continuous-timg forocasting of these dynamics.
Our extensive evaluations span benchmark datasets including the Revised MD17 dataset, five larger molecular systems (20 326 heavy atoms), and extended temporal horizons up to 15,000 MD steps. GF-NODE achieves state-of-the-art performance across all evaluated systems and time horizons, demonstrating superior long- horizon trajectory prediction while accurately preserving essential molecular geometries. Comprehensive structural analyses including radial distribution functions confirm that our model maintains realistic molccular correlations well beyond coordinate-level accuracy. The ablation studies further highlight the critical role of both the spectral decomposition and the continuous- time dynamics in capturing complex temporalt interactions, with formal 5O(3) equivariance guarantees provided. Importantly, our theoretical insight through heat equation analysis simplified diffusion models, combined with comprchensive empirical validation on real molecular dynamics trajectories, establishes rigorous foundation for the spatial-temporal correspondence that underlies our approach. The empirical validation through our Theoretical Framework for Joint Spatial-Temporal Analysis demonstrates quantitative relationships between spatial frequency modes and temporal dynamics scales across diverse molecular systems, confirming that spectral decomposition effectively separates distinct temporal evolution patterns.
These findings suggest that incorporating thysics-informet spectral decomposition principles into neural architectures represents promising direction for advancing molecular dynamics modeling, with broader implications for understanding and predicting complex multiscale phenomena in chemical and biological systems.
6 Broader Impact Statement
Our work contributes to advancing molecular dynamics simulations, which have broad implications in scientific discovery, particularly in drug design, materials science, and biomolecular modeling. By improving the accuracy and efficiency of multiscale molecular predictions, our approach could accelerate the discovery of novel therapeutics and facilitate the design of functional materials with tailored properties.
While method relies data-driven modeling, it does not replace physics-based simulations but rather augments them, reducing computational costs while maintaining interpretability. As with any machine learning-driven approach in scientific domains, care must be taken ensure model reliabilily, particularly in high-stakes applications such as drug development. Further validation and collaboration with domain experts will be essential to maximize the positive socictal impact of our work while mitigating risks related to model uncertainty.
7 Acknowledgement
This work was partially supported by NSF 2211557, NSF 2119643, NSF 2303037, NSF 2312501; the NAIRR Pilot Program (NSF 2202693, NSF 2312501); the SRC JUMP 2.0 Center; Amazon Research Awards; Snapchat. Gifts; and C-CAS (Center for Computer Assisted Synthesis, NSF CHE-2202693).
16

---

Published in Transactions on Machine Learning Research (06/2025)
References
Michael Bear and Estela Blaisten- Barojas. Atomic and ion collisions in solids and at surfaces: Theory, simulation and applications. American Scientist, 86(2):188, 1998.
Bernard R Brooks, Charles Brooks III, Alexander D Mackerell Jr, Lennart Nilsson, Robert Petrella, Benoît Roux, Youngdo Won, Georgios Archontis, Christian Bartels, Stefan Boresch, al. Charmm: the biomolecular simulation program. Journal of computational chemistry, 30(10):1545-166) 1614, 2009.
Ricky TQ Chen, Yulia Rubanova, Jesse Bettencourt, and David K Duvenaud. Neural ordinary differential equations. Advances neural information processing systems, 31, 2018.
Stefan Chmiela, Alexandre Tkatchenko, Huziel E. Sauceda, Tgor Poltavsky, Kristof T Schütt, and Klaus-Robert Müller. Machine learning of accurate energy-conserving molecular force fields. Science advances, 3(5): c1603015, 2017.
Stefan Chmiela, Valentin Vassilev-Galinde Oliver T Unke, Adil Kabylda, Igor Poltavsky, Steven Vandenbrande, Michael Willatt, and Alexandre Tkatchenko. Machine learning force fields and coarse-grained variables in molecular dynamics: application to materials and biological systems. Journal of Chemical Theory and Computation, 19120)77000-7027, 2023.
Anders S Christensen and O Anatole von Lilienfeld. On the rolc of gradients for machine learning of molecular energies and forces. Machine Learning: Science and Technology, 1(4):045018, 2020.
Fan R. K. Chung. Spectral Graph Theory. American Mathematical Society, Providence, R.I., 1997.
Wendy D Cornell, Piotr Cieplak, Christopher Bayly, lan R. Gould, Kenneth M Merz, David M Ferguson, David C Spellmeyer, Thomas Fax, James W Caldwell, and Peter A Kollman. A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. Journal of the American Chemioal Society, 117((9):5119 5197, 1995.
Pawan Goyal and Peter Benner. Neural ordinary differential equations with irregular and noisy data. Royal Society Open Science, 10(7):221475, 2023.
Samuel I Holt, Zhaozhi Qian, and Mihacla van der Schaar. Noural laplace: Learning diverse classes of differential equations in the laplace domain. In International Conference on Machine Learning, pp- 8811 8832. PMLR, 2022.
Zijie Huang, Yizhou Sun, and Wei Wang. Learning continuous system dynamics from tregularly-samplle partial observations. Advances Neural Information Processing Systems, 33:16177-16187, 2020.
Zijie Huang, Yizhou Sun, and Wei Wang. Coupled graph ode for learning interacting system dynamics. In Proceedings of the 27th ACM SIGKDD conference on knowledge discovery & data mining, pp. 705 715, 2021.
Zijie Huang, Yizhou Sun, and Wei Wang. Generalizing graph ode for learning complex system dynamics across environments. In Proceedings of the 29th ACM SIGKDD Conference on Knowledge Discovery and Data Mining, pp. 798 809, 2023a.
Zijie Huang, Wanjia Zhao, Jingdong Gao, Ziniu Hu, Xiao Luo, Yadi Cao, Yuanzhou Chen, Yizhou Sun, and Wei Wang. Tango: Timo-reversal latent graphode for multi-agent dynamical systems. arXio preprint arXiv:2310.0644 0.47777, 2023b.
Jürg Hutter. Car parrinello molecular dynamics. Wiley Interdiscipliaary Reviews: Compudational Molecular Science, 244):604 612, 2012.
Patrick Kidger. On neural differential equations. arXiv preprint arXiv:2202.02.0 2022.
Suyong Kim, Weiqi Ji, Sili Deng, Yingbo Ma, and Christopher Rackauckas. Stiff neural ordinary differential equations. Chaos: An Interdisciplinary Journal of Nonlinear Science, 31(),, 2021.
17

---

Published in Transactions on Machine Learning Research (06/2025)
Leon Klein, Andrew Foong, Tor Fjelde, Bruno Mlodozeniec, Mare Brockschmidt, Sebastian Nowozin, Frank Noé, and Ryota Tomioka. Timewarp: Transferable acceleration of molecular dynamics by learning time-coarsened dynamics. Advances in Neural Information Processing Systems, 36, 2024.
Takeshi Koshizuka, Masahiro Fujisawa, Yusuke Tanaka, and Issei Sato. Understanding the expressivity and trainability of lourier neural operator: A mean-field perspective. In The Thirly-eighth Anveual Conference on Neural Information Processing Systems.
Nikola Kovachki, Samuel Lanthaler, and Siddhartha Mishra. On universal approximation and error bounds for fourier neural operators. Journal of Machine Learning Research, 22(290):1-76, 2021.
Nikola Kovachki, Zongyi Li, Burigede Liu, Kamyar Azizzadenesheli, Kaushik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Neural operator: Learning maps between function spaces with applications to pdes. Journal of Machine Learning Research, 24(89):1-99 17,, 2023.
Z Li, DZ Huang, B Liu, and Anandkumar. Fourier neural operator with learned deformations for pdes on general geometries. arxiv, 2022a.
Zijic Li, Kazem Mcidani, Prakarsh Yadav, and Amir Barati Farimani. Graph neural networks accelcrated molecular dynamics. The Journal of Chemical Physics, 156(14), 2022b.
Zongyi Li, Nikola Kovachki, Kamyar Azizzadenesheli, Burigede Liu, Kanshik Bhattacharya, Andrew Stuart, and Anima Anandkumar. Fourier neural operator for parametric partial differential equations. arXiv preprint Xiv::20010. 088899,, 2020.
Ning Liu and Siavash Jafarzadeh. Yue yu. domain agnostic fourier neural operators. ar preprint arXiv:2305.00477 2023.
Xiao Luo, Viyang Gu, Huiyu Jiang, Hang Zhou, Jinsheng Huang, Wei Ju, Zhiping Xiao, Ming Zhang, and Yizhou Sun. Pgode: Towards high-quality system dynamics modeling. In Forty-forst International Conference on Machine Learning.
Xiao Luo, Jingyang Yuan, Zijie Huang, Huiyu Jiang, Yifang Qin, Wei Ju, Ming Zhang, and Yizhou Sun. Hope: High-order graph ode for modeling interacting dynamics. In International Conference Machine Learning, pp. 23124-23139 PMLR, 2023.
Lutz Molgedey and Heinz Georg Schuster. Separation of a mixture of independent signals using time delayed corrolations. Physical review letters, (2(23)33634 1994.
Outi MH Calo-Ahen, Ida Alanko, Rajendra Bhadanc, Alexandre MJJ Bonvin, Rodrigo Vargas Honorato, Shakhawath Hossain, André H Juffer, Aleksei Kabedev, Maija Kakkoen- Anders Støttrup Larsen, et al. Molecular dynamics simulations in drug discovery and pharmaceutical development. Processes, 9(1):
71, 2020.
Victor Garcia Satorras, Emiel Hoogeboom, and Max Welling. E (n) equivariant graph neural networks. In International conference on machine learning, pp. 9323-9332 03-3233..222 PMLR, 2021.
Mathias Schreiner, Ole Winther, and Simon Olsson. Implicit transfer operator learning: Multiple time- resolution models for molecular dynamics. Advances in Neural Information Processing Systems, 36, 2024.
Paulo CT Souza, Riccardo Alessandri, Jonathan Barnoud, Sebastian Thallmair, Ignacio Faustino, Fabian Grünewald, Ilias Patmanidis, Haleh Abdizadeh, Bart MH Bruininks, Tsjerk A Wassenaar, et al. Martini 3: a general purpose force field for coarso-grained molecular dynamics. Nature methods, 18(4):382-3)) 3888228,..- 2021.
Antonis I Vakis, Vladislav A Yastrebov, Julien Scheibert, Lucia Nicola, Danicle Dini, Clotilde Minfray, Andreas Almqvist, Marco Paggi, Seunghwan Lee, Georges Limbert, al. Modeling and simulation in tribology across scales: An overview. Tribology International, 125:169- 2018.
18

---

Published in Transactions on Machine Learning Research (06/2025)
Quan-De Wang, Jing-Bo Wang, Juan-Qin Li, Ning- Xin Tan, and Xiang- Yuan Li. Reactive molecular dynamics simulation and chemical kinetic modeling of pyrolysis and combustion of n-dodecane. Combustion and Flame, 158(2):217 226, 2011.
Minkai Xu, Jiaqi Han, Aaron Lou, Jean Kossaifi, Arvind Ramanathan, Kamyar Azizzadenesheli,, Jure Leskovec, Stefano Ermon, and Anima Anandkumar. Equivariant graph neural operator for modeling 3d dynamics. orXiv preprint urXiv:2401.1100 11037, 2024.
Chengxi Zang and Fei Wang. Neural dynamics on complex networks. Proceedings of the 26th ACM SIGKDD International Conference Knowledge Discovery & Dala Mining, pp. 892-902, 2020.
19

---

Published in Transactions on Machine Learning Research (06/2025)
A
Experiment Setup
In this section, we describe key details about the datasets used (MD17 and alanine dipeptide), inchuding the simulation step sizes, and then outline the hyperparameter choices and training procedures.
A.1 Datasets
We use the MD17 dataset (Chmiela et al. 2017) for small-molecule dynamics, and the alanine dipeptide dataset (Schreiner et al. 2024) for conformational analysis in proteins. Statistios for MD17 are provided in Table ( From the table, we observe Benzene has very marked drifting that is much larger that the scale of the molecule itself. The MD17 dataset was generated using ab initio molecular dynamics (MD) simulations with a time slep of 0.5 femtoseconds (Is). The simulations were conducted in the NVT ensemble 500 K for total duration of 2000 picoseconds (ps). The final dataset was created by subsampling the full trajectory, preserving the Bolzell-Boltzmann distribution for the energies. For the alanine dipeptide dataset, we use the protein fragment's trajectory recorded at step size of 1.0 ps.
Table 6: Summary statistics for molecular structures, inchuding the number of atoms, position extrema (Xmin: Xmax; meea); and velocity extrema. (Vmin, V mmxx, V mean). The numbers of modes used get the optimal results are also listed.
Molecule
# Atoms
# #Modes used
Xmin
Xmax
Xmean
Vinn
V max
V mean
benzene
4
178.112
197.981
27.737
-0.004
0.003
-0.000
aspirin
13
6
-3.720
3.105
0.026
-0.111
0.012
0.000
ethanol
3
3
-1.398
1.417
-0.004
-0.011
0.010
-0.000
malonaldehyde
5
5
-2.397
2.370
0.000
-0.010
0.009
0.000
naphthalene
10
4
-2.997
2.593
0.000
-0.012
0.011
0.000
salicylic
8
-2.734
2.581
0.051
-0.013
0.012
-0.000
toluene
7
-1.90
2.630
0.015
-0.010
0.012
0.000
uracil
6
-2.338
2.558
0.012
-0.012
0.011
0.000
A.2 Training Setup
Hyperparameters. We train all models using the Adam optimizer at a learning rate of 1 x 10 and apply weight decay of 1 x 10-15 for regularization. Each experiment runs for 5000 epochs, processing batches of size 50 each training step. For molecular trajectory prediction, we set the sequence length (the number of timesteps) to 8, mcaning each training sample contains 8 frames from the overall simulation. In main experiments, we fix the hidden feature dimensionality to 64 and use maximum future horizon of 3000 simulation steps when constructing training samples. We also specify the dopri5 ODE solver with relative and absolute tolerances of 1 x 10-3 and 1 x 10-4, respectively, to integrate the continuous-tin model components. Although many of these choices (c.g., total layers, solver tolerances) can be altered, we find these particular settings maintain a good balance of accuracy and computational efficiency.
Training Procedure. During training, each mini-batch is formed by sampling short segments of length 8 from the molecule's dynamics trajectory. The model then predicts future positions of atoms after continuous-tim evolution, and the Mean Squared Error (MSE) between predicted coordinates and ground-truth coordinates is minimized. We checkpoint models whenever validation performance improves, and at the end of training, we report results using the best-performing checkpoint according to the validation set. In addition, regular evaluations on the test set help track the model's generalization to unseen trajectories. The overall process can be summarized as: (1) load training samples, (2) form mini-batches of molccular frames, (3) perform forward pass through the model to generate predictions, (4) compute the MSE loss, (5) update model parameters via backpropagation, repeating for each epoch until convergence.
20

---

Published in Transactions on Machine Learning Research (06/2025)
By following these procedures with our chosen hyperparameters, have observed stable convergence across both MD17 and alanine dipeptide datasets, as well as strong generalization to different. segments of the trajectory during test runs.
A.3 Loss Function
To train the model, we use the Mean Squared Error (MSE) loss, which measures the difference between the predicted atomic positions and the ground truth positions at each predicted time point. Given that the goal is to predict the molecular conformations at time points 4, 1,,..., the MSE is calculated as follows:
Let x E R be the ground truth 3D coordinates of atom at time iú. and x R be the predicted coordinates for the same atom at time ty. The MSE loss defined as:
K N
MMSE
ΣΣ ==
(26)
NK
where N is the number of atoms and K is the number of time points.
This loss encourages the model to minimize the Euclidean distance between the predicted and actual atomic positions across all time steps, ensuring accurate trajectory prediction for the molecular system.
The MSE loss is applied at each time point, thus aligning the predicted future states with the true molecular dynamics trajectory. The model is trained by minimizing LASE over all predicted time points.
21

---

Published in Transactions on Machine Learning Research (06/2025)
B Theoretical Guarantees on Spectral Decompositions
Below, present concise mathematical exposition on the theoretical underpinnings of the Graph Fourier Transform (GFT) decomposition used in our framework. We explain how the eigenaale--igeeveeoor structure of the graph Laplacian L induces a decomposition of graph signals into low-frequency (global) and high- froquency (local) modes, and we justify truncating to the first M modes. We follow standard nomenclature in spectral graph theory (Chung 1997).
B.1 Preliminaries and Definitions
Definition B.1 (Graph Laplacian). Let G = (V,E) be an undirected graph with N = M vertices. Let A E N be its adjacency matrix, and let D be the diagonal degree matrix, where
D(1,)) = Σ A(i,)).
(27)
1=1
The graph Laplacian is defined
== D A.
(28)
It is well known that L real symmetric and positive semidefinite. fact, the eigenvalues of the Laplacian matrix real and non-negative.
Definition B.2 (Graph Fourier Transform (GFT)). Given the eigen-decomposition
= UAUT
(29)
where
A diaggA AN-1). 0 0= do S S..... AN-1,
(30)
and U = [uo uuj ... uN- stores the corresponding orthonormal eigenvectors in columns. For graph signal
x (x1, "2,,.,,,,,,y IN) E RRV,
(31)
the Graph Fourier Transform (GFT) is given by
U
(32)
and the inverse GFT is
UR.
(33)
B.2 Truncation and Mode Selection
Lemma B.3 (Approximation Error for Spectral Truncation). Let E R be ony graph signal, and let X(M) be ils spectral apprecimation oblained by keeping the first M modes. Then.
N-1 Σ
Ixx
x
1 x
(34)
Moreover, x x is a-bundlimited the sense that.
u x = for all de > a.
(35)
then choosing M such that AM-1 S a yields an exact recovery x X(M)-
Proof. See the main text for details. We expand the signal in the Laplacian eigenbasis {ux}, and observe that discarding all modes with & > M removes the corresponding frequency components.
22

---

Published in Transactions on Machine Learning Research (06/2025)
B.3 Low-Frequency High-Frequency Modes
Because L is posilive semidefinite and the eigenvalues { increase with i, smaller eigenvalues correspond to slow, global variations, while larger eigenvalues capture more oscillatory, local phenomena. Proposition B.4 (Global vs. Local Spatial Scales). Let uk be the k-th eigenvector of L with eigenvalue Ag- Suppose x encodes atomic coordinates or their latent features. Then:
If dR small, the corresponding mode up represents slowly varying (global) deformations across the molecule.
If A large, the corresponding mode un represents rapidly changing (local) structural varialions.
Proof. From standard results in spectral graph theory (Chung 1997 The low-frequency (small X) modes vary smoothly across edges, whereas high-frequency (large X) modes exhibit large differences across edges.
B.4 Practical Mode Truncation Criteria
Definition B.5 (Mode Retention Threshold), For a desired tolerance e O, select M such that
N-1
1 x < 1
(36)
k= M
In practice, one may also pick M based on AM-1 S a, ignoring modes where d > a. Corollary B.6 (Error Control via Low-Pass Approximation). Under the some notation as above, if
Σ R=M
1 x S <|x11
(37)
then
- xto))l2 |x|12
(38)
Hence, discording high-frequency modes exceeding this threshold leads to bounded approzimation error.
23

---

Published in Transactions on Machine Learning Research (06/2025)
C Temporal Dynamics of Graph Fourier Modes
We illustrate on the canonical graph heat equation how Laplacian eigenvalues directly determine the time-scale of cach Fourier modc.
Proposition C.1 Eleat-Eqquaion Mode Dynamics). Let 9 (V. E) be a graph with normalized Laplacion L admitling
L=UAUT, A (tiag(No........................ AN-1), 0 = λo S AN-1- Consider the graph-heat evolution on scalar signal f(t) IR:;
d((()) - [())
(39)
dt
Write the kth Graph. Fourier coefficient
(1 us T.
S(1).
Then each mode evolces independently:
dax(t) ft(()
u f(t) dx uh f(t) da ax(t):
(40)
dt dt
with closed-form solution In particular:
ax ax) exp( (~x t).
exp(->x
Modes with D small decay slowly (long time-scales).
Modes with d large decay rapidly (short time-scales).
Proof. Starting from df f,, project onto the orthonormal eigenvector ug:
de
dak = U T df dt
L. (u, = = A = = dx (2F = ~1 akk
dt
The ODE dy - ag integrates immediately to ax au (0)e dst.
Remark C.2. Although our GF-NODE dynamics learned rather than the pure heat equation, Proposition 3.3 shows that under any diffusion-like operator the Laplacian cigenvalues A set intrinsic timo-scales for cach mode. In practice, our Neural-ODE learns a richer fe. but we still observe empirically that modes with larger A tend to exhibit faster temporal variation citiiiion-reciislly the behavior we exploit by evolving each og (and vector analogue) under separate ODE channels.
C.1 Empirical Validation of the Opatial-Temporal Scale Correspondence
Our central hypothesis that the spatial frequencies defined by the graph Laplacian eigenvectors correspond to the characteristic timescales of molecular motion. Specifically, we claim that low-frequency spatial modes (associated with small eigenvalues A&) capture slow, global dynamics, while high-frequency spatial modes (large A&) capture fast, local vibrations. To rigorously validate this claim, perform joint spatial-temporal frequency analysis on the ground-truth molecular dynamics trajectories.
The procedure is as follows:
Step Obtain the Eigendecomposition of the Graph Laplacian. For a given molecule, we construct. the graph and its Laplacian NxN, We then compute its eigendecomposition:
== UAU
where U
[ng u,...., ux-1] is the orthonormal matrix of cigenvectors, and A diag(>n; contains the corresponding real, non-negative eigenvalues sorted in ascending order (0 do I \ S ...).
These eigenvectors } form the spatial frequency basis.
24

---

Published in Transactions on Machine Learning Research (06/2025)
Step 2: Project the Ground-Truth Trajectory onto the Spatial Basis. Let a ground-truth trajectory be represented by sequence of atom positions (X(4)) where X(4) E NR is the matrix of coordinates for N atoms at time step 43. We first ensure translational invariance by mean-centering the coordinates at each step: Xx(tz) = x(t)) ( (tyy)
We then project these coordinates onto the spatial basis using the Graph Fourier Transform (GFT). For each spatial mode k; we obtain time series of its 3D spectral coefficient *x(t)) R². This is computed by projecting the centered coordinates onto the eigenvector un:
(t,) U Xc(t)) u X((tt)
This operation yields N distinct time series, {Et
(t))) Tuim, one for each spatial mode (0,............... (O. 1).
Step Analyze the Temporal Frequency of Each Spatial Mode. For each spatial mode k. we now have a signal (4) that describes how the amplitude of that spatial pattern evolves over time. To quantify its characteristic temporal frequency, we first compute the time series of its magnitude, sx(t,) (t))ll2
Next, we compute the power spectral density (PSD) of the signal s3 (t) using the Discrete Fourier Transform (DFT), commonly implemented via the Fast Fourier Transform (FFT). Let the PSD be Px(S), where s represents the temporal frequency.
To obtain single characteristic temporal frequency fx for each spatial mode k. we compute the power- power-weighte weighted average frequeney (i.e., the spectral centroid):
Jimax f df
J
Jiae Pa( df
where the integral is performed over the range relevant temporal frequencies up to the Nyquist limit fmar- This value fee represents the average timescale on which the spatial mode k is active.
Step 4: Visualize the Temporllll Correlation. Finally, we plot the characteristic temporal frequency fo against its corresponding spatial frequency (the Laplacian cigenvalue Aς). This creates set of points (Ax:s) for k = 1,....... . (we omit the k = 0 mode as it corresponds to the zero eigenvalue and has no dynamics). A clear positive correlation in this plot provides strong empirical evidence that spatially smoother modes (low As) indeed evolve more slowly (low fa), and spatially oscillatory modes (high Ax) evolve more quickly (high fi).
C.2 Theoretical Framework for Joint Spatial-Temporal Analysis
To rigorously ground our model's architecture, we introduce a formal framework for analyzing the joint spatial-temporal characteristics of molecular dynamics. We treat the molecular trajectory as time-varying signal defined on a graph and leverage tools from spectral graph theory and statistical signal processing to decompose and analyze its structure.
C.2.1 Molecular Configuration as a Signal in Hilbert Space
Let the molecular graph be G = EV, with | = N. The set of all possible scalar functions on the vertices, V:V V + R, forms an N-dimensional Hilbert space H. = (VV) equipped with the standard inner product
(f.g) = (()ry
The graph Laplacian L E NNN a self-adjoint, positive semi-definite operator on this space. By the spectral theorem, admits an orthonormal basis of eigenvectors [ N-11 for H. These eigenvectors are the natural analogues of the Fourier basis on Euclidean domains.
Definition C.3 (Spatial Frequency Basis). The eigenvectors {ux] of the graph Laplacian L, ordered by their corresponding eigenvalues 0 do S \ S .... S AN-1, are defined as the spatial frequency basis of the graph g. The eigenvalue A is interpreted the frequency of the mode ug.
25

---

Published in Transactions on Machine Learning Research (06/2025)
The frequency A quantifies the spatial variation of its mode. This is formalized by the concept of Total Variation. The Total Variation of graph signal f is given by the Laplacian quadratic form:
TV(f) = Lf Σ
w -
(F,)EE
For an eigenmode uk, its total variation is simply un Lue Ab. Thus, a small A implies the mode is spatially smooth, while a large d implics is highly oscillatory.
C.2.2 Spectral Representation of Time Varying Graph Signal
We formalize the molecular trajectory as a tiryiigrying, vector-valued graph signal, F: R + (E()))), which maps time & to the matrix of atomic coordinates ((t) E RNx3, To analyze this signal, we project it onto the spatial frequency basis using the Graph Fourier Transform (GFT).
The GFT of the trajectory F(t) is a set of time- time-varying spectral coefficients (t) where xx(() E R Each coefficient represents the projection of the molecular configuration onto specific spatial mode at time to
x (t) (t),, u == (u X. (t), u 2()) u (t))
The inverse GFT perfectly reconstructs the signal: =tt)= = E N-1 (t) From dynamical systems perspective, the evolution of the molecule, governed by complex coupled equations in the vertex domain, can be viewed as the superposition of the dynamics of these individual modes (x(() in the spectral domain.
C.2.3 Frequency Analysis of Temporal Mode Dynamics
Our core hypothesis is that the dynamics of these modes are not uniform; specifically, the characteristic temporal frequency of (t) should be correlated with its spatial frequency A To quantify this, we analyze the temporal signal associated with each spatial mode.
For each mode k we define a scalar time series representing its energetic contribution, sx(1) | (t)l.. The thcoretical foundation for analyzing the frequency content of such a signal is the theorem, which connects the power spectral density (PSD) of wide-sense stationary process to the Fourier transform of its autocorrelation function. Assuming local stationarity in the dynamics, we estimate the PSD, denoted (P)()) by computing the squared magnitude of the temporal Fourier transform (FFT) of the signal (!).
The PSD Px(f) reveals how the energy of the spatial mode k is distributed across different temporal frequencies f. To summarize this distribution into single charactoristic frequency, we compute the spectral centroid fe:
fx 10 f df Jo df
(41)
The centroid Fü represents the power-weighted average frequency, providing a robust measure the dominant timescale of the mode's dynamics.
C.2.4 The Joint Spatial- Temporal Spectrum
By performing this analysis for each spatial mode (I...., N -1} we can construct an empirical joint spatial-temporal spectrum of the molecular dynamics. This spectrum the set of points:
S = [(x.. No c R x R+
C.3 Empirical Validation Results
We validate our spatial-temporal correspondence hypothesis through comprehensive analysis of four distinct molocular systems: buckyball-aatcher (MD22), salicylic acid, naphthalcne, and uracil (rMD17). Following the four-step protocol outlined above, we analyze ground-truth molecular dynamics trajectories quantify the correlation wetween spatial frequencies x and characteristic temporal frequencies ft-
26

---

Published in Transactions on Machine Learning Research (06/2025)
Experimental Setup and Data Processing. For each molecular system, we apply the four-step protocol described in the preceding section. The molecular systems span diverse chemical environments: buckyball- catcher represents a large supramolecular complex with 120 heavy atoms, while salicylic acid, naphthalene, and uracil are smaller organic molecules with distinct aromatic and heterocyclic structures. This diversity enables assessment the generality of the spatial-temporal correspondence across different molecular archilectures and dynamical regimes.
Quantitative Correlation Analysis. We observe consistently positive correlations between spatial frequencies As and characteristic temporal frequencies f& across all molecular systems (Figure 9). The correlation coefficients range from moderate (r = 0.424) to very strong (r = 0.816), with all systems exhibiting positive slopes in the log-log representation. Specifically, buckyball-catcher exhibits r = 0.424 with slope a = 0.552, salicylic acid shows r = 0.816 with slope a = 2.642, naphthalene demonstrates r = 0.622 with slope a = 5.016, and uracil yields r = 0.716 with slope a = 4.098. These results provide quantitative validation of the predicted power-law relationship fa x Ag.
Physical Interpretation and Theoretical Validation. The consistently positive correlations provide empirical validation of our central hypothesis: spatial modes with larger Laplacian eigenvalues (high spatial frequency) exhibit faster characteristic temporal dynamics (high temporal frequency). This validates the theoretical prediction that the eigenspectrum of the molecular graph provides natural ordering of dynamical timescales, with smooth, collective modes evolving slowly and oscillatory, localized modes evolving rapidly.
The variation in correlation strength and slopes across molecules reflects differences in their structural complexity and dynamical behavior. The buckyball-catcher system exhibits moderate correlation (r 0.424) due to its complex multi-scale dynamics as large supramolecular assembly. In contrast, the smaller organic molecules demonstrate stronger correlations (r 0.622 0.816), consistent with more regular vibrational spectra and clearer timescale separation.
Connection to Heat Equation Dynamics. These empirical findings provide direct validation of our theoretical derivations based on heat equation mode dynamics. In the diffusion framework, solutions of du Lu take the form uy (t) = (0)e Ast, where the decay timescale тÔ- = 1/Ax is inversely proportional dt
the eigenvalue. This predicts positive correlation between spatial frequency d and temporal frequency fs x 1/Ts Av-
The observed positive slopes (a 0.552 5.016) confirm the power-law relationship fx Xg predicted by our theoretical analysis. The consistently positive slopes validate that the graph Laplacian eigenspectrun provides a natural frequency ordering for molecular dynamics, bridging spatial structure and temporal evolution through spectral graph theory.
Implications for Physics-Informes Neural Architectures. The demonstrated spatial-temporal corre- spondence validates graph Laplacian eigenmodes as a plysics-informed basis for neural network architectures in molecular dynamics. By aligning representational capacity with the natural frequency hierarchy, Fourier- based models can achieve more efficient and physically meaningful learning. The positive correlations further suggest that truncated representations retaining low-frequency spatial modes will preferentially capture slow, collective motions dominating cong-timescale behavior, providing theoretical justification for dimensionality reduction strategies in molecular machine learning.
27

---

Published in Transactions on Machine Learning Research (06/2025)
D Formal Proof of
for GF-NODE Pipeline
Below formal proof SO(3) (rotational) equivariance for our GF-NODE pipeline, closely following the style EGNO's Appendix proofs. We focus on the 3D rotational part of SE(3); translations can be handled by the separate mean-centering step (sce remarks below). Our proof is broken down into:
1. Defining the R-action,
2. Showing that each module (Fourier transforms, block-diagonal ODE, EGNN layers) is SO(3)- cquivariant, and
3. Composing these results to conclude overall equivariance.
D.1 Formal Statement of Equivaaanee
Let
f [n.
(42)
be a function describing the node features of 3D molecular system over some (possibly temporal) domain D. Concretely,
Jh D D Nx xx collects invariant (scalar) node fcatures,
fz: Rx collects equivariant (3D) features (positions, velocities, etc.).
Denote by R SO(3) 3D rotation matrix. The action of R on f is defined by
(R ()
(43)
which rotates only the Z-component in R and leaves the scalar h-component invariant. We claim that our overall GF-NODE operator To satisfies
To()) (R( f),
(44)
i.e., To (0(3)--quivariant, Formally:
Theorem D.1 (SO(3) Equivariance). Let To be the GF-NODE architecture composed of:
1. An EGNN encoder (mapping [fn. fz] + encoded features),
2. Mean-centering and Graph Fourier Transform (F),
3. A dick--iaaonal Neural ODE in the spectral domain,
4. Inverse GFT (F-1) plus adding back the mean, and
5. An EGNN decoder.
Then for any R E (O(3),, the pipeline satisfies
To T f) = TT(f..
(45)
In other words, rotating the input 3D features by R is equivalent to applying Te first and then rotating the result.
We prove this via the following steps:
1. Lemma 1: EGNN layers are (0(3)-
2. Lemma 2: GFT and its inverse are (3)))-equivart (dimension-wise linearity).
3. Lemma 3: The block-diagonal Neural ODE in spectral space preserves Reequivariance on the vector channels.
4. Conclusion:
Composing these yiclds the full pipeline's equivariance.
Below, we provide the details of each lemma and then the final proof of the Theorem.
28

---

Published in Transactions on Machine Learning Research (06/2025)
D.2 EGNN Equivariance
Lemma D.2 (EGNN layers are (Oui)eequiaariaant). Consider a generic EGNN luyer ④ which updales (h,, x,) + (h1 x(),
(46)
using message passing:
m = 2e (h,, h - X))
(47)
h Φb h,, Σ míj
(48)
X = x; +... (xi - Xz).
(49)
Then for any rotation E SO(3),
1 (Rxx, hi) (h', Rx.).
(50)
Hence & is -quiivranvva its 3D imputs.
A standard proof () shows that each update depends on x; Xj, which under a global rotation R(x; - -xy)
transforms consistently to yield R xx The same argument applies to 3D velocities (or any additional 3D vectors).
D.3 GFT Equivariance
We next show that the (inverse) Graph Fourier Transform is eqqiivaraa with respect to dimension- rotations of the 3D features.
Lemma D.3 (GFT and F-1 are equivaaraa) Let F be the dimension-wis wise GFT mapping a function
fz D + N (mmx3)
(51)
to its frequency coefficients F(ʃz) E ((modea)xmx3, Under R SO(3), define
R (Fiz) = Ffz but with each 3-D channel rotated by R.
(52)
Then
R F(Z)== 1z).
(53)
Similarly, F-1 is
in the sense that
RR F-
(54)
Proof Sketch. The GFT (and its inverse) act lincarly along each 3D axis. If R. rotates the 3D channels, we can commute R with the linear transform F. Precisely as in the EGNO proof, the multilinear expansions show that R = fz). The same argument applies to F-1 because it is also linear and dimension wise.
D.4 Block-Diagone Neural ODE Equivariance
In the GF-NODE pipeline, once we have GFT coefficients Ž, the Neural ODE acts as a diakonia operator:
H 7.
fo
H
(55)
go
Z
where Ž E ((modes)xmx3. Rotating on those 3D channels amounts to mixing the coordinate axes lincarly. Because the ODE is chosen to channelwise or "blockwise" linear or MLP-based, commutes with R. Hence:
Lemma D.4 (Block-Diagonal ODE is S0(()-equivariaat) For each frequency mode, the update on Ż is dimension-wste wise (like separate channel). A global rotation R that mixes Z 22 2 con be factored out of the ODE solution-So
R go (2) = go (R 2).
(56)
Integrating over preserves this properly.
29

---

Published in Transactions on Machine Learning Research (06/2025)
D.5 Proof of the Main Theorem (50(3) Equivariance)
Recall our overall operator To has the form:
1. EGNN Encode: (h, (h, 2) + (h(4),7()))
2. Mean-Center + GFT:: (44) + = + (()) = F(Zi
3. Block-Diagonal ODE: ( H Z() for any t.
4. Inverse GFT + Add Mean: 2(() + Zc(l)= = [ (2(1)) - 2(1).
5. EGNN Decode: (h, Z(t) + (h' (t), (t)).
To show R. Ta(s) = To f), we proceed step-by-step:
1. EGNN Encode: By Lemma A.1, if the input positions are replaced with Rx,, the output is Rx!")
2. Mean-Center Under global rotation, the centered coordinates also rotate, i.e, x; H Rx..
3. GFT: By Lemma A.2, dimension-wisc GFT on Xx yields the rotated spectral coefficients.
4. Block-Diagonal ODE: Lemma A.3 says the ODE in spectral space is equivariant w.r.t. 3D axis mixing, so R commutes with the ODE solution.
5. Inverse GFT: Again by Lemma A.2, inverse transforms are linear in each dimension, preserving R on the output.
6. Add Mean: The final global shift (if any) is consistent with R.
7. EGNN Decode: By Lemma A.1 again, if the input to the decoder is rotated, the output is the rotated version of the unrotated output.
Hence each sub-module respects the action of R. Composing them order yields the final statement
Ta(R f) = R Te(f).
(57)
This completes the proof of equivvaani..
D.6 Remarks Translations
In practice, SE(3) includes translations as well. Our pipeline removes the translational degree of freedom by mean-centering the positions (the DC mode). A global translation x; > x; simply shifts the mean x, so the centered coordinates x remain unchanged. This effectively "factors out" translation before the GFT steps. When we re-add the mean at the end, it ensures the final positions transform by x, H X; + μ. Thus the entire pipeline remains invariant to translations (i.e., translates its outputs accordingly). For brevity, SO(3); translation invariance follows from the mean-subtraction
the above proof focuses on rotations procedure plus the decoder's reliance relative positions.
30

---

Published in Transactions on Machine Learning Research (06/2025)
E
Additional Experiment Results
Below, we provide four tables corresponding to different experimental comparisons. The first table reports performance on the MD17 dataset with regular (equi-) timesteps. The remaining three tables focus on ablation experiments conducted under irregular sampling conditions: (i) ablations on the alanine dipeptide dataset, (ii) ablations comparing different GNN architectures, and (iii) ablations on different temporal embedding approaches.
Table 7: MD17 with Regular (Equi-) Timestep Sampling. MSE (x10-2 A2) on the MD17 dataset using regular timesteps. Best results are in bold, and second-best are underlined.
Aspirin
Benzene
Ethanol
Malonaldehyde
Naphthalene
Salicylic
Toluene
Uracil
NDCN LG-ODE
31.73:10.40
56.2110.30
10.7410.02
46.55:10.28
2.25:10.01
3.5810.11
13.92.1..2
2.38:10.00
19.3610.12
53.9211.32
7.0810.01
24.41110.03
1.7310.02
3.8210.04
11.1810.01
2.1110.02
EGNN
9.2410.07
57.8512.70
4.63:10.00
12.8110001
0.38110.1
0.8510.00
10.4110.04
0.5610.02
EGNO
9.41:10.09
55.1311322
4.6310.00
12.8110.01
0.40:10.01
0.9310.01
10.43:10.10
0.5910.01
ITO
20.56:10.03
57 57.2510.58
8.6011.27
28.4410.73
1.8210.17
2.4810.34
12.47:0.300
1.33:10.12
Ours
6.07-0.09
1.51+0.07
2.74+0.01
9.43+00.0
0.24+0.02
0.63++0.05
1.80++0.03
0.41+00.2
Explanation. Table shows results for MD17 under regular (evenly spaced) sampling scheme. Although the dataset inherently has fine-grained timesteps, constrain both training and evaluation to equidistant frames to compare methods fairly. Our approach demonstrates consistent improvements over baselines on nearly all molecules.
Table 8: Ablation Results on Alanine Dipeptide (Irregular Sampling). MSE (x10-3 nm²). Best results in bold.
standard
no ode
_o__ade
no_ode__ h
no mo__de___ no interaction
MSE
4.48++0.07
4.72+-0.05
4.60+0.05 time_mly
4.64+0.05
4.80+0.04
interaction_ _onnccat
posenc
FFT
no_fourrer
MSE
4.9810.06
4.62.10.05
4.60.10.5
4.57:10.05
4.5110.04
Explanation. Table organizes the ablation settings into two rows, each containing five columns. The first. row compares our "standard" model to variants that remove specific ODE blocks or modify ODE updates ("no_ "oto__dde""""" "no" "no__de______________ _dd__h, [___odd__")) and the second row compares different interaction modes, time embeddings, and Fourier settings. The "standard" configuration achieves the best overall MSE.
Table 9: Ablation on GNN Architectures (Irregular Sampling). MSE (x10-2 A2) on MD17 comparing different GNN layers (SAGEConv, GCNConv, EGNNConv). Best results in bold.
Aspirin
Beuzene
Ethanol
Malonaldehyde
Naphthalene
Salicylic
Toluene
Uracil
SAGEConY
6.46:00.03
1.79
2.74=0.05
10.54:0001
0.23+0.02
0.63:0.01
3.08-0.05
0.41#0.01
GCNCouv
6.002
10.08
3.09.10.06
10.05
0.42.10.01
0.8000.00
1.80.0.05
0.61.10.2
EGNNConV
8.85+0.02
40.86=0.98
4.41+00006
12.49+0.000
0.40=0.01
0.87=0.01
8.63-0.04
0.62=0.02
Explanation. Table 9 shows how our model performs with different GNN backbones on MD17 under irregular sampling. Overall, SAGECCovv yields robust performance for most molecules, whereas GCNConv provides better results specifically on Benzene and Toluene. EGNNConV performs well on some local metrics but struggles on large translations (i.e., Benzene).
Explanation. Table compares three different time-enbedding methods under irregular timestep sampling: positional encoding (posenc), small MLP (mlp), and direct concatenation of time tokens (concat).
31

---

Published in Transactions on Machine Learning Research (06/2025)
Table 10: Ablation on Time Embedding Approaches (Irregular Sampling). MSE (x10-2 A2) on
MD17 across different time encoding schemes.
Aspirin
Benzene
Ethanol
Malonaldehyde
Naphthalene
Salicylic
Toluene
Uracil
posenc
6.01=0.08
18110002
3.11=0.05
10.69=0.03
4.1910.04
0.87=0.01
3.5620.07
0.5910.02
mlp
6.61=0.06
1.6140.03
3.08=0.02
10.62=0.04
0.4140.01
0.87=0.02
3.25+0.06
0.56+0.01
concat
6.46:10.03
1.52=0008
2.7410.005
10.54=0.01
0.23
0.6310.01
1.8040.05
0.41-00001
Concalenation achieves the lowest MSE, suggesting that a straight.frwwrd inclusion of time in the feature vector can be beneficial, though the MLP variant also achieves competitive performance on several molecules.
Table Comparison of GF-NODE with baseline models on the revised MD17 dataset at At 3000. MSE (x10-² A values; best results in bold.
Model
Aspirin
Azobenzene
Ethanol
Malonaldehyde
Naphchalene
Paracetamol
Salicylic
Toluene
Uracil
NDCN
34.7820057
2.
21.67
39. 02.02
1. 1.282004
27.13:0.41
1.08
25. 23.99
0.88
LG ODE
33.4020.15
9.8820.31
23.1520.17
41 41.212001
1.42:0.00
26. 722222
1 1.3330005
24.75:2..7
0.9520.08
EGNN
31.45
11.03:0.41
22.95
80.0.
1.18
25.57=0.30
1.20+0.04
23.90+a.19
0.82+0.02
EGNO
32.01=0.8a
7.51+0.12
23.58
a7.90
1.37
26.02=0.58
0.88+a.02
24.82+0.66
0.78+0.04
38.50a1.z
10.87.0.53
25.3320.7
13.55
1.6020.09
1.6649.07
27.35
1.12:0.11
GF-NODE
30.27.0.0
7.0320.02
21.9220.05
37.9220.05
1.1020.01
24.46
0.81 20.01
28.18
0.62..001
Table 12: Comparison of GF-NODE with baseline models on the revised MD17 dataset at At 10000. MSE ( (x10-2 A2) values; best results in bold.
Model
Aspirin
Azobenzene
Ethanol
Malonaldehyde
Naphthalene
Paracetamol
Salicylic
Toluone
Uracil
NDON LG-ODE
42.67.0.91
113422072
29.4522.47
48 75
1.9022.05
83.332
1.95
34.1220.18
1.7220.00
16.12.0.87
9.88
31.0G20.a
20.80
2. 2.:300.07
29. 09.05
1.652020
30.48+0.19
1.22+0.04
EGNN
38.0920.16
13.67
36.1420.20
11.96
2.07=s0.1
28. 45
1.27
29.83
1.01
ECNO
40.9020.54
12.30
27. 27.88
42.33
2. 2.22:0.04
27. 12:0.47
1.58
12.1522022
0.95
ITO
49.77
15.03:0007
24.1
53.5020.73
2.5620.14
35.9820.05
2.2240.17
36.45sb0.51
1.832007
GF-NODE
33.18:0.05
7.29
22.31=0.04
38.74
1.27#0.01
27.20=0.04
0.93
27.92:00.4
0.72:0.01
Table 13: MSE (x10-2 A for Ala2 and larger molecules At 3000. Best results in bold, second best underlined.
Model
Ala2
Ac-Ala3-NHMe
AT-AT-CG-CG
Bucky-Catcher
DW Nanotube
NDCN
1877
22.34++.22
26.78=0.30
6.10=0.15
4.50=0.20
LG-ODE
90.15+0.90
30.12+1.00
33.50=1.10
8.25=0.40
5.80=2.30
EGNN
56. 70..0....84
18. 18.4511..1
20.75-0.45
7.10
5. 5.603
EGNO
69. 772.58
23. 1040.35
17.20-0.20
5.30+0.10
4.50
ITO
269.45 1.87
28. 10.95
32.00
8.60
3.80
GF-NODE
44.82=0.71
13.19=0.13
14.07+0.23
3.09+0.01
2.58+0.02
32

---

Published in Transactions on Machine Learning Research (06/2025)
Table 14: MSE (x10-2 A?) for Ala2 and larger molecules At 10000. Best results in bold, second best underlined.
Model
Ala2
Ac-Ala3-NHMe
AT-AT-CG-CG
Bucky-Catcher
DW Nanotube
NDCN
134.10=0.48
30.15+0.26
38.82=0.60
7.32=0.18
5.85=0024
LG-ODE
117.20
40.66.1..0
48.58 1.32
9.90
7.54
EGNN
88.63 10.36
20.11
30.09 0.54
8.52
7.28.10.2
EGNO
11:11
31.1940.42
24.94=0.24
6.36=0.12
5.85
ITO
297.21
39.02
46.40=1.50
10.32
4.94
GF-NODE
49.20=0.31
16.72=0.14
17.89+0029
4.37
3.22+0.05
Aspirin
Azobenzene
Ethanol
18:.5
18
60.0
55
37.5
5981
10
221
38.0
32.5
40
15
20.0
28
Malonaldehyde
Naphthalene
Paracetamal
10:004
42.5
10:232
40.0
37.5
30.0
32.5
MSS
30.0
27.3
25.0
3.0
22.5
Salicylic
Toluene
Uracil
1003
6541
00-000
40
609
(1400:
7.25
0081
5/400C
2.00
0400
UEE boto
3.50
20
3.25
125
1.00
1.00
0.75
9000 20
0.50
2000
4030
8000
10000
12000
14000
2000
4000
10000 12000
14000
2000
4000
1900
10000
1200U
14000
22
A8
61
Figure 10: Temporal error growth for GF-NODE and baseline models on nine molecules on the revised MD17
dataset. Each panel plots MSE (x10-2 Á?) versus integration horizon At 1000, 2000, .... 15000.
33

---

Published in Transactions on Machine Learning Research (06/2025)
Ala2
Ac AlaM--MMMe
AT-AT-CG-CO
300
15-000
15020
230
1655 10E0000-- 84
150E OCCE- 201
D902: 00:00-10000000000000 201
200
130
1t
Bucky-Catcher
DW Nanotube
12
15-00
1E-00C
102E (X20
8000 10000
2000
Figure 11: Long-horizon MSE trends for five larger molecules (Alaz, Ac-Alag-NHMe, AT-AT-CG-CG. Bucky- Catcher, DW Nanotube). Each panel shows MSE (x10-2 A?) for GF-NODE and five baselines over integration horizons Al- = 1000, 2000, ...> 15000.
Table 15: Number of heavy (non-H) atoms in each molecule.
Molecule
Dataset
Heavy Atoms
Revised MD17
Aspirin
13
Azobenzene
14
Ethanol
Malonaldehyde
Naphthalene
10
Paracetamol
11
Salicylic acid
10
Toluene
7
Uracil
Larger Molecular Systems
Alao
22
Ao-Alay-NHMe
20
Bucky-Catcher
120
AT-AT-CG-CG
76
DW Nanotube
326
34

---

Published in Transactions on Machine Learning Research (06/2025)
Epoch Time Number Modes
(NVIDIA L40 GPU, AMD Ryzen 7950X CPU)
100
Epoch (seconds)
70
100
150
200
250
300
Number of Modes
Figure 12: Epoch time (seconds) as function of the number of Fourier modes used, measured on an NVIDIA L40 GPU with an AMD Ryzen 7950X CPU. Error bars represent variability across three repeated timing runs at each mode count.
Ala2
Ac-Ala3-NHMe
AT-AT-CG-CG
N00: 10300- 80
1-30-
20
100
5
70.,de
FFT
10, Souner
Original
n0.0de
10, Soutrrrr
Original
n0.000
FFT
7O, Souree
Origiall
Varient
Bucky-Catcher
DW Nanotube
085E (K10 R1
821
1010
Griginal
Figure 13: Ablation study on hierarchical components for five larger molecules (Alaz, Ac-Alag- NHMe, AT-AT- CG-CG, DW Nanotube) at At = 10000. Variants shown are: no no_ode (red), FFT only (blue), no fourier (orange), and the full model (Original, green). Removing ODE or Fourier components degrades performance often exceeding baseline errors whereas the complete architecture attains the lowest MSE.
35

---

Published in Transactions on Machine Learning Research (06/2025)
F Radial Distribution Function Analysis
Section summary. This appendix contrasts the RDFs predicted by GF-NODE with reference ab initio data, both system averaged and element specific, to evaluate structural fidelity.
Radial distribution functions (RDF, ((r)) quantify how atomic density varies as a function of distance and therefore provide a stringent test of whether a Icarned model reproduces the local and intormediaao--aago structure of condensed phase systems. Below we compare the RDFss produced by GF-NODE with those obtained from reference ab inátio trajectories (blue).
Radal d6et tice Functionn reulaed
(a) Uracil
(b) Salicylic acid
(c) Naphthalene
Rodia sritution -attelalllcctttee
Roda Dicritudion Funciion
32
mm
8.
100
(d) Buckyball catcher
(c) Double-war wall NT
Figure 14: System averaged radial distribution functions g(r) for the five benchmark molecules/ molecules/comppexees complexes. Orange: GF GF-NODE; black: ab initio. The close match indicates that the model accurately reproduces both short- and medium-range order.
(a) Uracil: C-C
(b) Uracil: C-N
(c) Uracil: C-O
(d) Uracil: N-N
0. Betheehheeteeee...o www.ottooooooo uned
(c) Uracil: N-O
(f) Uracil: O-O
Figure 15: Element-spedtie RDFs (r) for uracil showing six unique heavy-atom pairs. GF-NODE reproduces both the peak positions and intensities of the reference curves.
36

---

Published in Transactions on Machine Learning Research (06/2025)
Hola nebdd
30
(a) Salicylic: C-C
(b) Salicylic: C-O
(e) Salicylic: O-O
Figure 16: Element-specific RDFs for salicylic acid.
C.C Radial Distriution Ruedioo repintalent
ww
Pali Podbllr
24-
(a) Naphthalene: C-C
(b) Bucky catcher: C-C
(c) DW NT. C-C
Figure 17: Carbon-carbon RDFs for thrce purely carbonaccous systems. GF-NODE captures the first
peak 1x4 Å) and the longer-range range oscillations characteristic of aromatic stacking and nanotube wall spacing.
37


Reviews:
[{"review_id": "C1OKPfSaZh", "reviewer": "Reviewer_MJAX", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n1. Novel Architecture: The proposed GF-NODE architecture presents an interesting and conceptually novel combination of Graph Fourier Transform, Neural ODEs, and GNN components for molecular dynamics modeling.\n2. Ablation Studies: The paper includes ablation studies investigating the contribution of different components (e.g., spatial decomposition method, ODE evolution, mode interactions, time embeddings), providing some insight into the model design.\n3. Detailed Analyses Attempted: The authors attempt several detailed analyses, such as long-term trajectory prediction, temporal super-resolution, TICA analysis for slow modes, and the effect of the number of Fourier modes.\n\nWeakness:\n\n1. Empirical Performance and Significance: One of the core weakness lies in the empirical results. The performance gains of GF-NODE over simpler baselines like EGNN (which lacks explicit temporal propagation and spectral decomposition) appear marginal on the tested datasets (Table 1, Table 3). This raises significant questions about the practical utility and added value of the proposed complex architecture (GFT + Neural ODE + Decoder) compared to existing, potentially simpler methods. The claim that the added components significantly improve long-term (there was no explanation whether the data exhibits specific long term behaviour and it was not properly cited what timescales are normally understood as long-term --> the times which were looked at are only double of the length of the times of the EGNO, which don't represent long-term behaviour) prediction is not convincingly supported by the quantitative results presented.\n2. Evaluation Scope and Rigor:  The evaluation is limited to relatively small molecules (MD17) and one peptide (alanine dipeptide). Demonstrating efficacy on more complex systems exhibiting clear multi-scale dynamics (e.g., larger proteins, solvation dynamics) is necessary to validate the approach's claims for practical MD simulations. The interpretation lacks comparison with ground truth  and clear metrics for \"preserving essential geometric properties.\" Definitions of \"long-horizon\" and \"stability\" in this context are vague.  The TICA analysis (Figure 5) lacks comparison to the TICA scores of the ground truth simulation, making it difficult to assess the absolute quality of slow-mode capture.\nThe justification for using MSE, especially under irregular sampling with potentially noisy MD data, is not provided. How baselines requiring regular steps were adapted for irregular sampling is unclear. As mentioned in previous work the data follows a stochastic propagation, why is MSE a proper metric to measure performance?\n3. Conceptual Clarity and Justification:\nThe link between spatial graph frequency modes (from GFT) and specific temporal dynamics (high/low frequency in time) is asserted rather than rigorously established or empirically verified, especially when GFT is taken in the latent space. Why should low graph Laplacian eigenvalues correspond to slow temporal evolution and high eigenvalues to fast temporal evolution? (See annotation re: Prop C.4). The proof provided for Prop C.4 merely restates standard spectral graph theory interpretations of smoothness, lacking deductive rigor regarding dynamic scales.\nThe motivation for evolving spectral modes independently via decoupled ODEs (Eq. 14-15) is not explained properly; the claim of \"coordinated but distinct manner\" needs better explanation, especially when interaction via attention is later introduced and shown to be beneficial (Sec 4.6), within the modes. The implementation details and justification for the attention mechanism across modes are sparse.\nWhy the Neural ODE component inherently captures temporal multi-scale behavior better than alternatives needs clearer articulation and evidence. How does evaluating the ODE at arbitrary points directly translate to capturing dynamics at different underlying physical timescales?\nThe motivation against prior work (Sec 1, Sec 2) often lacks specific citations or evidence for claims like inadequacy (e.g., \"fail to adequately represent the interplay...\"). Comparison to relevant SOTA diffusion models seems missing.\n4. Presentation and Methodological Details:\nAs noted in the overall assessment, the presentation can be confusing. Specific details are often unclear or imprecise: the exact nature of node features and embeddings (Sec 3.2, Fig 3 linkage), the graph construction (distance threshold vs. chemical bonds, Sec 3.1), the definition of \"interactions\" (Sec 3.1), the rationale for feature choices (Sec 3.1), normalization details (Sec 3.2), how FFT is applied to graphs (Sec 4.6), implementation of mode interactions (Sec 4.6), ODE solver details (dopri5 source/citation missing, Sec 3.4).\nThe term \"Neural Operator\" seems loosely applied; the learned component appears to be functions defining the ODEs, not necessarily operators mapping between function spaces in the canonical sense.\nThe analysis of the optimal number of modes (Figure 8, Table 6) is confusing: Why does variance seem to increase after the optimum? How can the optimal number of modes exceed the number of atoms for some molecules? Why weren't the purportedly optimal (higher) numbers used in the main experiments?", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper introduces Graph Fourier Neural ODEs (GF-NODE), a novel deep learning architecture designed to model the spatio-temporal dynamics of molecular systems, particularly addressing the challenge of capturing interactions across multiple scales. The core idea is to combine a Graph Fourier Transform (GFT), based on the graph Laplacian's eigenmodes, with a Neural Ordinary Differential Equation (Neural ODE) framework that models the continuous-time evolution of these spectral coefficients independently or interactively. The predicted future state in the spectral domain is then transformed back to the physical domain using an inverse GFT, potentially followed by refinement with a GNN decoder. The authors hypothesize that this explicit separation and modeling of spatial frequency modes coupled with continuous-time dynamics enable GF-NODE to better capture both high-frequency local fluctuations (e.g., bond vibrations) and low-frequency global effects (e.g., conformational changes) compared to existing methods. The approach is evaluated on the MD17 dataset and alanine dipeptide benchmark, comparing its predictive accuracy (MSE), preservation of geometric features (bond lengths/angles), and long-term stability against several baseline methods, including GNNs and other Neural ODE/Operator approaches.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Critical:\n1. Demonstrate Significant Empirical Advantage: Provide stronger evidence that the proposed GF-NODE architecture offers substantial benefits over simpler, relevant baselines (like EGNN, potentially others). This could involve evaluating on more complex datasets where multi-scale effects are prominent and where simpler models might fail, or showing significantly better performance/scaling on existing datasets. If gains are marginal, the paper needs to convincingly argue for other advantages (e.g., interpretability, specific physical regimes) supported by evidence.\n2. Improve Conceptual Rigor and Justification: Critically re-evaluate and clarify the core conceptual claims. Provide stronger theoretical arguments or empirical analysis linking spatial GFT modes to distinct temporal dynamics/scales. Justify the claim that the Neural ODE formulation inherently captures temporal multi-scales effectively in this context. Rigorously define and prove claims like Proposition C.4 or state them as assumptions/heuristics. Strengthen the motivation against prior work with specific evidence/citations.\nEnhance Clarity of Method and Implementation: Revise Section 3 significantly to provide precise details on graph construction, feature definition/normalization/embedding, the GFT process (especially for vector features), the exact structure of the Neural ODE (including interaction mechanisms like attention), and the decoder. Clarify notation (e.g., hidden dimensions). Provide necessary citations (e.g., for dopri5). Improve figure clarity (Fig 3).\n3. Refine Evaluation and Analysis: Re-evaluate or provide a much clearer interpretation of the long-term prediction results, explicitly addressing the large MSE values and comparing against ground truth behavior (including drift/rotation). Provide clear metrics and evidence for claims about preserving geometric features. Include ground truth TICA scores for comparison. Justify metric choices (MSE) and clarify baseline adaptation for irregular sampling. Clarify the analysis related to the number of modes (Fig 8, Table 6). Also the Figure 4 shows explicitly how Truth and Predicition don't align, the lower red atom exhibits rotation in the truth, but it is not captured in the in prediction.\n\nStrengthening:\n1. Expand Dataset Evaluation: Justify the choice of MD17/alanine dipeptide. Including revised MD17 and potentially larger/more complex systems (water, other biomolecules) would significantly strengthen the validation. Define \"simulation steps\" in physical time units.\n2. Clarify Sampling: Explain the random sampling strategy for training data (Sec 4.2) more clearly (e.g., are endpoints fixed?).\n3. Discuss SOTA: Briefly position the work relative to other potentially relevant SOTA approaches (e.g., diffusion models).", "broader_impact_concerns": "The authors include a Broader Impact Statement (Section 6) that appropriately identifies potential positive impacts (accelerating drug design, materials science) and acknowledges the need for model reliability and validation, particularly in high-stakes applications. This seems adequate.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "6JgIcCbsQu", "reviewer": "Reviewer_458H", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n- The GF-NODE pipeline of projecting node features into a Laplacian eigenbasis, evolving each mode with a Neural ODE, and reconstructing via inverse transform is an inuitive and  compelling approach to learning multiscale graph dynamics.\n- The reductions in MSE on MD17 and alainine dipeptide benchmarks are significant.\n- The paper is well-structured with clear algorithmic descriptions, illustrative figures of spectral modes and trajectories, and several ablation studies.\n\n**Weaknesses**\n-  It is very surprising that the EGNN is shown to be outperforming EGNO in the equal timestep sampling setting, which appears to be the same experimental setup as in the original EGNO paper). The EGNO paper reports improved results over EGNN [1]. Moreover, the ablation studies in Tables 4 and 5 show much better performance than existing baselines even for the naive settings. For example, replacing the spectral decomposition by a multi-layer perception still achieves almost the same performance, which indicates that the main architectural components of GF-NODE is nonessential.\n- The claim that GF-NODE explicitly handles multiscale spatial and temporal interactions as a novel contribution may be overstated, since many existing GNN/ODE models can implicitly capture such interactions via deep message-passing layers or learned temporal filters. \n- Experiments only focus on small molecules with less than roughly 20 atoms. Not only is it unclear whether the performance gains generalize to larger molecules or systems, but it is also unclear whether the MD17 and alanine dipeptide benchmarks require capturing _truly_ long-range spatial interactions, which is claimed to be one of the main benefits of GF-NODE.\n- There is no analysis of computational cost as the number of retained modes or graph size grows, which is critical for applying GF-NODE to larger molecular systems. One point of concern is that the time complexity of spectral decomposition is $O(n^3)$, so the method might not scale well with a larger number of nodes.\n\n**References**\n\n[1] Minkai Xu, Jiaqi Han, Aaron Lou, Kamyar Azizzadenesheli, Stefano Ermon, & Anima Anandkumar. (2024). Equivariant Graph Neural Operator for Modeling 3D Dynamics.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper introduces Graph Fourier Neural ODEs (GF-NODE), a framework for molecular dynamics modeling that combines graph Fourier transforms with Neural ODEs to effectively capture correlations across multiple spatial and temporal scales. GF-NODE first decomposes node features into Laplacian eigenmodes via a graph Fourier transform, evolves each mode through a block-diagonal Neural ODE augmented with multi-head attention, and finally reconstructs atomic geometries via the inverse transform. The graph spectral decomposition is based on the intuition that small Laplacian eigenvalues correspond to global spatial deformations while larger eigenvalues correspond to local deformations.\n\nGF-NODE is then benchmarked against state-of-the-art baselines on the MD17 and alanine dipeptide datasets. Empirical results demonstrate substantial gains in trajectory prediction,​ and ablation studies suggest that continuous-time ODE propagation and inter-mode communication help the model achieve these improvements.", "claims_and_evidence_yn": "No", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "**Changes critical for acceptance**:\n1. Ensure fair baseline comparisons using original papers’ implementations and hyperparameters. Re-implementations of EGNN, EGNO and ITO must reproduce their reported performance.\n2. Reconcile why all models in the ablations (including naive models that don't use Fourier-based decomposition, inter-mode communication, or any continuous-time evolution) outperform current SOTA models.\n3. Evaluate on large molecular systems which contain longer-range interactions. \n\n**Changes that would strengthen the paper**:\n1. Include a runtime and scalability study quantifying how training/inference time grows with node and mode count.\n2. Replace MD17 with the Revised MD17 dataset [1].\n\n**References**\n[1] Anders S. Christensen and O. Anatole von Lilienfeld (2020) \"On the role of gradients for machine learning of molecular energies and forces\" https://arxiv.org/abs/2007.09593", "broader_impact_concerns": "No concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "jn8WT1Jf8a", "reviewer": "Reviewer_BeZR", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths:**\n\n1. The integration of spectral methods with neural ODEs is a novel and potentially impactful approach that may appeal to a broad audience.  \n2. The manuscript presents compelling analyses and results concerning irregularly sampled time-series data.  \n3. The work includes an evaluation of long-term simulation stability, where the proposed method demonstrates improved robustness compared to MLFF-based approaches.\n\n**Weaknesses:**\n\n1. The manuscript lacks a thorough discussion of short-range versus long-range interaction modeling. Given that spectral methods are particularly well-suited for capturing long-range interactions, an explicit analysis along these lines would be beneficial. Moreover, evaluation on larger and more complex systems—such as the MD22 double-wall nanotube—could strengthen the empirical validation.  \n2. The manuscript does not include a formal proof of the model’s equivariance. While this may appear pedantic, providing a brief but rigorous proof would enhance the theoretical soundness and completeness of the work, and should be considered in the revised version.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The author propose a graph neural ODE for molecular dynamics simulation. The idea is simple, first maps the node features (both scalar and vectorial) to the spectral domain. Then fit a neural ODE on the spectral domain. Finally, the node features are mapped back via inverse Fourier transform. The author verified their methods on MD17 and Alanine Dipeptide with irregular time-series sampling.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "1. **Include Loss Function Equation**:  \n   Explicitly present the mathematical formulation of the loss function used during model training to clarify the optimization objective.\n\n2. **Add Proof of Equivariance**:  \n   Introduce a dedicated section that formally proves the equivariance properties of the proposed model, including the group actions involved and their preservation under the model’s operations.\n\n3. **Compare with MLFF Methods**:  \n   Add a section discussing the relationship between the proposed method and existing MLFF approaches. Highlight some theoretical distinctions.\n\n4. **Extend to MD22 Benchmark**:  \n   Demonstrate the applicability of the proposed model by evaluating its performance on the MD22 dataset.\n\n5. **Velocity Autocorrelation Function (VACF) Analysis**:  \n   Plot and compare the velocity autocorrelation functions across MLFF-based methods, the proposed approach, and DFT simulations to assess dynamical consistency.\n\n6. **Additional MD Metrics Beyond MAE**:  \n   Complement MAE with other molecular dynamics metrics such as the Radial Distribution Function (RDF)-based stability.\n\n7. **Discussion on Energy Conservation**:  \n   Address whether the proposed method conserves energy during MD simulations. If energy is not strictly conserved, include a discussion explaining why and its implications for stability and physical realism.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)