Learning residue level protein dynamics with multiscale Gaussians

## Abstract

Abstract Many methods have been developed to predict static protein structures, however understanding the dynamics of protein structure is essential for elucidating biological function. While molecular dynamics (MD) simulations remain the in silico gold standard, its high computational cost limits scalability. We present DynaProt , a lightweight, SE(3)-invariant framework that predicts rich descriptors of protein dynamics directly from static structures. By casting the problem through the lens of multivariate Gaussians, DynaProt estimates dynamics at two complementary scales: (1) per-residue marginal anisotropy as 3 × 3 3\times 3 covariance matrices capturing local flexibility, and (2) joint scalar covariances encoding pairwise dynamic coupling across residues. From these dynamics outputs, DynaProt achieves high accuracy in predicting residue-level flexibility (RMSF) and, remarkably, enables reasonable reconstruction of the full covariance matrix for fast ensemble generation. Notably, it does so using orders of magnitude fewer parameters than prior methods. Our results highlight the potential of direct protein dynamics prediction as a scalable alternative to existing methods.

## 1 Introduction

Proteins rarely exist in static conformations. Due to interactions with ligands, other biomolecules, and external factors such as temperature and pH, protein structures continuously fluctuate. Many enzymes rely on loop motions or domain rearrangements to form catalytically active sites (Zinovjev et al., 2024), allostery often involves shifting backbone or side chain conformations propagating signals over long distances (Yu and Koshland Jr, 2001), and even membrane proteins such as GPCRs switch between inactive and active conformational states essential for signal transduction (Zhang et al., 2024). Clearly, protein structures are dynamic. Understanding these dynamics is central for mechanistic insight and, potentially, the design of functions (Guo et al., 2025; McCammon, 1984).

Capturing this ensemble-level behavior computationally has long been the domain of molecular dynamics (MD). MD simulates the time evolution of atoms under a force field, generating high-resolution conformational trajectories from which fluctuations, covariances, and time-dependent observables can be derived. MD remains the in silico gold standard for protein dynamics, offering fine-grained, physically grounded insights (Shaw et al., 2010; Hollingsworth and Dror, 2018; Childers and Daggett, 2017). However, it comes with an enormous computational cost: simulating 100 ns of dynamics for a single protein can take days or weeks on specialized hardware. This limits its scalability, especially for proteome-wide applications or tasks requiring real-time dynamics estimates.

Recent work has explored the use of deep learning to approximate and accelerate this process. Generative modeling-based methods like AlphaFlow (Jing et al., 2024) repurpose AlphaFold2 (Jumper et al., 2021) under a flow matching paradigm to sample protein conformations. Along this vein, Lewis et al. (2024) recently introduced BioEMU as a large-scale diffusion model pretrained on PDB (Burley et al., 2017) and AFDB structures, and fine-tuned on 200 ms of MD data, to efficiently generate protein conformations. Other methods, like MSA subsampling, make inference-time adjustments to the MSA input of AlphaFold2, yielding the structural ensembles (Del Alamo et al., 2022; Wayment-Steele et al., 2024; Stein and Mchaourab, 2022). Still, all of these approaches necessitate large scale PDB pretraining and suffer from inference-time computational overhead, requiring multiple stochastic forward passes to generate meaningful structural diversity. Moreover, while these ensembles can be used to approximate protein dynamics, generating them remains time intensive, and the full ensemble of diverse conformations may not always be necessary. In many practical settings, compact and interpretable representations of dynamics often suffice. This motivates the need for models that can explicitly predict such dynamics descriptors without relying on implicit dynamics learners like expensive ensemble generation methods.

Current explicit dynamics predictors, like FlexPert3D (Kouba et al., 2024), resort to predicting simple collective variables like per-residue RMSF, a scalar quantifying each residue’s positional fluctuation. RMSF is widely used due to its simplicity and interpretability, but it is fundamentally limited: it captures only the magnitude of local motion and discards directionality and residue-residue coupling. Similarly, Wayment-Steele et al. (2025) trained Dyna-1 to predict labels of $\mu$s–ms motion by cleverly exploiting missing chemical shift assignments as hidden observables in NMR ensembles, but these predictions also remain scalar and lack directionality. A different example of an explicit dynamics predictor is Normal Mode Analysis (NMA), a classic technique that approximates dynamics by identifying low frequency eigenmodes to describe the largest movements (Cui and Bahar, 2005; Skjaerven et al., 2009). NMA does not learn from data however, and instead operates solely on the input PDB structure. It can estimate the principal global directions of motion and offers insights into collective flexibility, but is sensitive to input structure quality and fails to adequately capture local anisotropy or conformational heterogeneity (Ma, 2005). This raises a natural question:

Figure: Figure 1: Dynamics methods information content vs. efficiency.
Refer to caption: x1.png

We introduce DynaProt, a lightweight, interpretable, and expressive framework for predicting protein dynamics through the lens of Gaussian distributions over structure (Section [2](https://arxiv.org/html/2509.01038v1#S2.F2)). Specifically, DynaProt predicts: (1) per-residue marginal anisotropy as $3\times 3$ covariance matrices capturing local dynamics while encompassing RMSF, and (2) joint scalar $N\times N$ covariances encoding pairwise dynamic coupling across residues. Remarkably, while DynaProt was not explicitly trained to directly model the full $3N\times 3N$ joint distribution, we find that its marginal and pairwise outputs can be composed into a reasonable approximation (Section [3.4](https://arxiv.org/html/2509.01038v1#S3.SS4)), enabling extremely fast ensemble generation in $\mathbb{R}^{3N}$. DynaProt is trained on only $\sim$1,000 MD-derived proteins, without large-scale pretraining on PDB structures, and improves upon Normal Mode Analysis (NMA) in both predictive accuracy and efficiency, while remaining dramatically smaller and faster than existing ensemble generation approaches.

To our knowledge, DynaProt is the first model to explicitly learn both marginal and pairwise Gaussian representations of protein dynamics, and the first to predict the full $3N\times 3N$ covariance structure—akin to NMA—in a data-driven, learnable fashion, rather than relying solely on analytical approximations or less informative per-residue fluctuations.

## 2 Gaussian representation of dynamics

Figure: Figure 2: Overview of protein dynamics under Gaussian view.
Refer to caption: figures/gaussianview.png

We propose a perspective for modeling protein dynamics through distributions over atomic coordinates, relying on tractable approximations such as Gaussians. Formally, we model a protein structure with $N$ residues as a random variable $\bm{X}\in\mathbb{R}^{3N}$, where each residue contributes the three-dimensional Cartesian coordinates of its Cα atom. While this coarse-grained representation omits side-chain flexibility, it enables scalable modeling of backbone dynamics, which is the scope of our work. We consider an ensemble to be $T$ independent samples after RMSD alignment. The distribution over conformational states is then represented as a multivariate normal distribution:

$$ $\displaystyle\bm{X}\sim\mathcal{N}(\bm{\mu},\bm{\Sigma_{\text{joint}}}),\quad\bm{\mu}\in\mathbb{R}^{3N},\;\bm{\Sigma_{\text{joint}}}\in\mathbb{R}^{3N\times 3N}$ (1) $$

Here, $\bm{\mu}$ corresponds to the average (or equilibrium) structure—typically the minimum energy conformation—and $\bm{\Sigma_{\text{joint}}}$ captures the full covariance across all Cα positions, encoding both local fluctuations and long-range correlated motions. This joint covariance matrix theoretically encodes all second-order information about the protein’s dynamics: from it, one can derive a wide range of collective variables including principal components (PCs) of motion, residue-residue distance variances, and global flexibility metrics.
The Gaussian formulation provides a principled way to decompose protein dynamics across different levels of granularity (Figure [2](https://arxiv.org/html/2509.01038v1#S2.F2) and Table [1](https://arxiv.org/html/2509.01038v1#S2.T1)), depending on the modeling objective. We highlight some of these below.

At the local level, the marginal distribution (Figure [2](https://arxiv.org/html/2509.01038v1#S2.F2).2) for a single residue $i$ is obtained by integrating out all other residue coordinates:
$p(\bm{x}_{i})=\int p(\bm{x}_{1},\dots,\bm{x}_{n})\,d\bm{x}_{\neg i},\text{where }d\bm{x}_{\neg i}:=\prod_{j\neq i}d\bm{x}_{j}$.
This results in a 3D Gaussian distribution over the Cα coordinates of residue $i$:

$$ $\displaystyle\bm{X}_{i}\sim\mathcal{N}(\bm{\mu}_{i},\bm{\Sigma}_{\text{marginal}}^{(i)}),\quad\bm{\mu}_{i}\in\mathbb{R}^{3},\;\bm{\Sigma}_{\text{marginal}}^{(i)}\in\mathbb{R}^{3\times 3}$ (2) $$

where $\bm{\Sigma}_{\text{marginal}}^{(i)}$ is the $3\times 3$ diagonal block of $\bm{\Sigma_{\text{joint}}}$. These marginals can be interpreted as Gaussian blobs encoding anisotropic local fluctuations—i.e., spatial variance of where each residue may reside.

Notably, this formulation allows for simple derivation of scalar flexibility metrics such as the root-mean-square fluctuation (RMSF) as $\text{RMSF}_{i}=\sqrt{\mathrm{Tr}(\bm{\Sigma}_{\text{marginal}}^{(i)}})$.
RMSF (Figure [2](https://arxiv.org/html/2509.01038v1#S2.F2).1) represents a simple notion of dynamics: a single scalar per residue quantifying positional fluctuation. However, it discards directional and covariance structure captured by the full marginal.

To capture dynamics beyond residue-local fluctuations, we also consider a covariance matrix $\bm{C}\in\mathbb{R}^{N\times N}$ of scalar pairwise coupling (Figure [2](https://arxiv.org/html/2509.01038v1#S2.F2).3). Each entry $\bm{C}_{ij}$ summarizes the dynamical coupling between residues $i$ and $j$, typically computed as a scalar projection of the corresponding $3\times 3$ block in the full joint covariance: $\bm{\Sigma}_{\text{joint}}[3i:3i+3,3j:3j+3]$. We choose MeanPooling as the scalar projection to compute each $\bm{C}_{ij}$. This compact representation enables efficient modeling of residue-residue coupling.

**Table 1: Taxonomy of protein dynamics representations under a Gaussian view.**
| Level | Description | Notation | Space | Captures |
| --- | --- | --- | --- | --- |
| 1 | Per-residue scalar (i.e. RMSF) | $\sqrt{\mathrm{Tr}(\bm{\Sigma}_{\text{marginal}}^{(i)}})$ | $\mathbb{R}^{N}$ | Magnitude of fluctuation per residue |
| 2 | Per-residue full (Gaussian blob) | $\bm{\Sigma}_{\text{marginal}}^{(i)}$ | $\mathbb{R}^{N\times 3\times 3}$ | Anisotropic local covariance per residue |
| 3 | Joint scalar (pairwise coupling) | $\bm{C}_{ij}$ | $\mathbb{R}^{N\times N}$ | Scalar covariance across all residues |
| 4 | Joint full covariance | $\bm{\Sigma}_{\text{joint}}$ | $\mathbb{R}^{3N\times 3N}$ | Full spatial covariance across all residues |

DynaProt focuses on levels 2 and 3 of this hierarchy—explicitly predicting both $3\times 3$ marginal Gaussians per-residue and a $N\times N$ matrix of residue-residue couplings. As noted before, from the $3\times 3$ marginals, we can easily derive RMSF (level 1). Interestingly, utilizing both the marginals and the pairwise coupling, we can retrieve a reasonable approximation of the full joint $3N\times 3N$ (level 4; Section [3.4](https://arxiv.org/html/2509.01038v1#S3.SS4)). This design strikes a balance between local interpretability and global coordination, while avoiding the intractibility of directly learning the full joint covariance.

## 3 Method

### 3.1 DynaProt Overview

DynaProt (Figure [3](https://arxiv.org/html/2509.01038v1#S3.F3)) consists of two models, each taking as input a protein structure but designed to explicitly capture different granularities of protein dynamics: (i) marginal Gaussian blobs per residue (Section [3.2](https://arxiv.org/html/2509.01038v1#S3.SS2)), and (ii) pairwise covariance across residues (Section [3.3](https://arxiv.org/html/2509.01038v1#S3.SS3)).
DynaProt is given the input structure as a set of local Cα residue frames. The frames are denoted $\{\bm{T}_{i}\}_{i=1}^{N}$, where each frame $\bm{T}_{i}\in\mathrm{SE}(3)$ is parameterized by a rotation matrix $\bm{R}_{i}\in\mathrm{SO}(3)$ and a translation vector $\bm{t}_{i}\in\mathbb{R}^{3}$. Simply put, $\bm{T}_{i}=(\bm{R}_{i},\bm{t}_{i})$ captures the local orientation and position of residue $i$. Additionally, an initial embedding layer is included to encode the amino acid sequence $\bm{s}\in\mathbb{R}^{N\times D}$.

Both models share a common architectural backbone composed of eight Invariant Point Attention (IPA) blocks from the structure module of AlphaFold2 (Jumper et al., 2021). These blocks are designed to encode geometric relationships between residues while maintaining invariance to SE(3) transformations (global transformations do not affect the learned residue-level representations). The IPA backbone processes the set of residue frames and the sequence representation, outputting a learned representation for each residue $\bm{h}\in\mathbb{R}^{N\times D}$. No pair representation is given as input to the model. The two models differ only in their readout layers, which we define in the proceeding sections.

### 3.2 Learning marginal Gaussians

After the input sequence representations and residue frames are processed through the IPA backbone, a simple MLP readout is used for marginal prediction. Given the hidden representation $h_{i}$ for each residue, the marginal readout outputs $\bm{\Sigma}_{\text{marginal}}^{(i)}\in\mathbb{R}^{3\times 3}$, modeling the local position (xyz) covariance of residue $i$. These outputs are trained to match empirical marginal distributions derived from the MD data. Note that the mean of each Gaussian is not learned. Instead, we take the input structure’s Cα coordinate $\bm{t}_{i}\in\mathbb{R}^{3}$ as the fixed mean $\bm{\mu}_{i}$ of the distribution: $\bm{\mu_{i}}:=\bm{t}_{i},\;\bm{X}_{i}\sim\mathcal{N}(\bm{\mu}_{i},\bm{\Sigma}_{\text{marginal}}^{(i)})$.
This assumption is motivated by the fact that the input structure usually corresponds to the experimentally determined (or AlphaFold-predicted) minimum energy conformation, and thus serves as a natural estimator of the ensemble mean. Consequently, the marginal prediction task reduces to learning the covariance matrices $\bm{\Sigma}_{\text{marginal}}^{(i)}$ alone.

Figure: Figure 3: DynaProt architecture.
Refer to caption: x2.png

##### Marginal dynamics module.

Recall that covariance matrices are required to be symmetric and positive definite (SPD). Predicting all 9 elements of a $3\times 3$ matrix would be overparameterized and does not guarantee SPD structure. Naively, one might consider symmetrizing an arbitrary matrix after predicting the 6 independent elements, but this only guarantees symmetry. Instead, we leverage the fact that any SPD matrix can be uniquely defined by its Cholesky factorization. Thus, we enforce SPD constraints directly by parameterizing the covariance via its Cholesky factor. Specifically, the model predicts the entries $\{a_{j}\}_{j=1}^{6}$ of a lower triangular matrix $\bm{L}_{i}\in\mathbb{R}^{3\times 3}$, enforces positivity along the diagonal with the Softplus activation function (Glorot et al., 2011), and recovers the covariance:

$$ $\displaystyle\bm{\Sigma}_{\text{marginal}}^{(i)}=\bm{L}_{i}\bm{L}_{i}^{\top},\quad\text{where }\bm{L}_{i}=\begin{bmatrix}\text{softplus}(a_{1})&0&0\\ a_{2}&\text{softplus}(a_{3})&0\\ a_{4}&a_{5}&\text{softplus}(a_{6})\end{bmatrix}$ (3) $$

This factorization ensures that the predicted covariance matrix is SPD by construction. Since SPD matrices lie on a Riemannian manifold with non-Euclidean geometry, using loss functions that respect this structure is critical for meaningful comparison. Standard Euclidean distances (e.g., MSE or Frobenius norm) ignore the curvature of this space and can lead to unstable or distorted gradients (see ablations in Appendix [A.3.2](https://arxiv.org/html/2509.01038v1#A1.SS3.SSS2)). We instead employ the log-Euclidean (or log-Frobenius) distance (Vemulapalli and Jacobs, 2015; Huang et al., 2015) that reflects the intrinsic geometry of the SPD manifold. The Bures-Wasserstein (Bhatia et al., 2019) distance can also be used, but we find the log-Frobenius distance to be more stable.

$$ $\displaystyle\mathcal{L}_{\text{LogFrob}}=\|\log(\bm{\Sigma}_{\text{pred}})-\log(\bm{\Sigma}_{\text{true}})\|_{F}^{2},\quad\text{where }\log(\bm{\Sigma})=\bm{Q}\log(\bm{\Lambda})\bm{Q}^{T}$ (4) $$

As Riemannian manifolds are “locally Euclidean", this loss applies the matrix logarithm mapping the SPD matrix to its tangent space where a Euclidean metric (canonical Frobenius norm) can be utilized.

### 3.3 Learning pairwise dynamics

Using the output representations $\bm{h}$ from the IPA backbone, the pairwise dynamics module produces a scalar-valued $N\times N$ covariance matrix $\bm{C}$, where each entry $\bm{C}_{ij}$ captures the dynamical coupling between residue pairs. These scalar couplings are derived from the full joint covariance matrix via averaging per block and trained to reproduce MD-derived pairwise fluctuations.

##### Pairwise Dynamics Module.

To predict the global pairwise covariance structure, we first construct pairwise features for all residue pairs. For each pair $(i,j)$, we concatenate their residue-level embeddings $[\bm{h}_{i}\,\|\,\bm{h}_{j}]\in\mathbb{R}^{2d}$ and project them into a lower-dimensional space: $f_{ij}^{(0)}=W_{\text{proj}}[\bm{h}_{i}\,\|\,\bm{h}_{j}]\in\mathbb{R}^{d^{\prime}}$.

We pass these features through a stack of AlphaFold-style pairwise attention blocks based on the Evoformer architecture (Jumper et al., 2021), which include triangle updates and residue-wise message passing. These operations are designed to model transitive and higher-order geometric dependencies across residue pairs, and have been shown to be highly effective in structure-aware tasks: $f_{ij}^{\text{attn}}=\text{PairwiseAttentionBlock}(f_{ij}^{(0)})\in\mathbb{R}^{d^{\prime}}$.
The output $f_{ij}^{\text{attn}}$ serves as a “learned basis"" over the space of residue-residue covariance structure. These basis features are then mapped to scalars through an MLP head, yielding a covariance for each pair of residues: $z_{ij}=\text{MLP}(f_{ij}^{\text{attn}}),\quad z_{ij}\in\mathbb{R},\quad\text{for }i\geq j$.
Following the same procedure as Section [3.2](https://arxiv.org/html/2509.01038v1#S3.SS2), we enforce SPD constraints on this covariance matrix by populating the lower-triangle entries of $\bm{L}\in\mathbb{R}^{N\times N}$ with the values of $z_{ij}$ and applying the Softplus activation when $i=j$. Finally, the pairwise covariance matrix is reconstructed via Cholesky composition $\bm{C}=\bm{L}\bm{L}^{T}$ and again equation [4](https://arxiv.org/html/2509.01038v1#S3.E4) is used for optimization.

### 3.4 Learning the full joint for ensemble sampling

##### Joint reconstruction heuristic.

Given a predicted scalar coupling matrix $\bm{C}\in\mathbb{R}^{N\times N}$ and a set of per-residue marginal covariances $\{\bm{\Sigma}_{\text{marginal}}^{(i)}\in\mathbb{R}^{3\times 3}\}_{i=1}^{N}$, we propose a heuristic to reconstruct an approximate full joint covariance matrix $\bm{\Sigma}_{\text{joint}}\in\mathbb{R}^{3N\times 3N}$.

Each marginal covariance $\bm{\Sigma}_{\text{marginal}}^{(i)}$ is SPD by construction, and thus admits a Cholesky factorization $\bm{\Sigma}_{\text{marginal}}^{(i)}=\bm{L}_{i}\bm{L}_{i}^{\top}$, where $\bm{L}_{i}\in\mathbb{R}^{3\times 3}$. We then define a block-diagonal matrix $\bm{L}_{\text{marginal}}\in\mathbb{R}^{3N\times 3N}$ as $\bm{L}_{\text{marginal}}=\bigoplus_{i=1}^{N}\bm{L}_{i}$. By construction, $\bm{L}_{\text{marginal}}$ is lower triangular with positive diagonal entries, since each $\bm{L}_{i}$ satisfies these properties.

Drawing from the univariate identity $\mathrm{Cov}(i,j)=\mathrm{Corr}(i,j)\cdot\sigma_{i}\sigma_{j}$, we define the multivariate cross-covariance block between residues $i$ and $j$ as $\bm{\Sigma}^{(i,j)}_{\text{joint}}=\bm{L}_{i}\widetilde{\bm{C}}_{ij}\bm{L}_{j}^{\top}$. Here, the Cholesky factor $\bm{L}_{i}$ serves as a matrix square root of the covariance $\bm{\Sigma}_{\text{marginal}}^{(i)}$, analogous to standard deviation in the univariate case. And, $\widetilde{\bm{C}}$ is a correlation matrix found by standardizing $\bm{C}$. Using the Kronecker product, we can denote this heuristic cleanly as follows,

$$ $\displaystyle\bm{\Sigma}_{\text{joint}}=\bm{L}_{\text{marginal}}\;(\widetilde{\bm{C}}\otimes\bm{I}_{3})\;\bm{L}_{\text{marginal}}^{\top}$ (5) $$

###### Proposition 3.1 (attach title to upper= ) .

x(SPD Closure). Given marginal covariances $\{\bm{\Sigma}_{\text{marginal}}^{(i)}\in\mathbb{R}^{3\times 3}\}_{i=1}^{N}$ and correlation matrix $\widetilde{\bm{C}}\in\mathbb{R}^{N\times N}$ to be symmetric and positive definite, then the reconstructed joint covariance
$\bm{\Sigma}_{\text{joint}}=\bm{L}_{\text{marginal}}(\widetilde{\bm{C}}\otimes\bm{I}_{3})\bm{L}_{\text{marginal}}^{\top}$
is also symmetric and positive definite.

We refer the reader to Appendix [A.1](https://arxiv.org/html/2509.01038v1#A1.SS1) for the proof.
This approximation combines local anisotropic uncertainty with global correlation structure. While not exact, we find it reconstructs the joint covariance to a reasonable degree and serves as a useful tool for downstream ensemble generation.

###### Proposition 3.1 (attach title to upper= ) .

##### Ensemble sampling.

Given the reconstructed joint covariance $\bm{\Sigma}_{\text{joint}}$ and our assumption that the mean $\bm{\mu}$ corresponds to the coordinates of the input structure (e.g., the PDB), we have now retrieved our Gaussian distribution over conformations $\mathcal{N}(\bm{\mu},\bm{\Sigma}_{\text{joint}})$. To sample from this distribution, we apply a multivariate generalization of the reparameterization trick used in univariate Gaussian sampling.

Note that this sampling relies directly on the Cholesky factor (similar to a matrix square root), mirroring the scalar case (Appendix [A.1](https://arxiv.org/html/2509.01038v1#A1.SS1)). Utilizing DynaProt predictions and this heuristic, ensemble sampling becomes extremely fast with minimal computational overhead.

## 4 Experiments

### 4.1 Setup

##### Preprocessing.

We construct ground-truth dynamics labels from the ATLAS molecular dynamics dataset, which comprises 1,390 proteins selected based on structural diversity using the ECOD domain classification (Vander Meersche et al., 2024). Following AlphaFlow (Jing et al., 2024) for preprocessing consistency, we concatenate each of the three replicate simulations of 100 ns per protein and extract the Cα coordinates. From each ensemble, we compute the empirical full joint covariance matrix over time and extract the relevant dynamics labels ($3\times 3$ marginals per residue and $N\times N$ residue coupling) as described in Section [2](https://arxiv.org/html/2509.01038v1#S2.F2).
We evaluate under two train/val/test split regimes. The primary matches AlphaFlow’s ($1265/39/82$), while comparisons to FlexPert3D use DynaProt trained on their topology-based split ($1112/139/139$). For naming, we refer to DynaProt-M for the model trained for marginals, DynaProt-J for the coupling predictions, and DynaProt for both.

##### Baselines.

For a faithful comparison, we mainly choose baseline methods that take a protein structure as input and predict dynamics descriptors either implicitly (AFMD+Templates, considered state-of-the-art) or explicitly (FlexPert3D, NMA). For NMA, we utilize the ProDy package (Zhang et al., 2021), specifically the Anisotropic Network Model instantiation. There are many methods worth noting that aim to predict ensembles or dynamics descriptors from sequence itself: standard AlphaFlow (AFMD), MSA-Subsampling, FlexPert-Seq, ESMDiff (Lu et al., 2024), ConfDiff (Wang et al., 2024), SeqDance (Hou and Shen, 2024), and SeaMoon (Lombard et al., 2024). With the advancements in structure prediction Abramson et al. (2024), we note that one could first predict a plausible structure from sequence and then predict dynamics using DynaProt. Thus, we compare against some of these sequence based methods as well (Appendix [A.3.3](https://arxiv.org/html/2509.01038v1#A1.SS3.SSS3)).

### 4.2 Predicting residue flexibility

**Table 2: RMSF Pearson correlation ($r$) against ATLAS MD-derived RMSF (FlexPert test split). Median and 75th percentile reported.**
| Method | RMSF $r$ ($\uparrow$) | # Params |
| --- | --- | --- |
| DynaProt-M | 0.865 / 0.930 | 955 K |
| FlexPert-3D | 0.830 / 0.899 | 1.2 B |
| NMA (ANM) | 0.697 / 0.784 | – |

Since DynaProt-M is trained to predict marginal Gaussians per residue, it inherently captures residue-level flexibility, as RMSF is defined as the square root of the trace of each marginal covariance (see Section [2](https://arxiv.org/html/2509.01038v1#S2.F2)). To evaluate DynaProt-M’s ability to recover this flexibility, we compare against what is, to our knowledge, the only method that *explicitly* predicts residue flexibility: FlexPert-3D. For fair comparison, we train and evaluate DynaProt-M under the same ATLAS train/val/test split defined in  Kouba et al. (2024). Despite solving the more challenging task of predicting full anisotropic marginals rather than scalar fluctuations alone, DynaProt-M achieves a substantially higher Pearson correlation with MD-derived RMSF (median $r=0.865$, 75th percentile $r=0.930$) than FlexPert-3D (Table [2](https://arxiv.org/html/2509.01038v1#S4.T2)), while using three orders of magnitude fewer parameters (955K vs. 1.2B) and without normal mode analysis as input. This allows DynaProt-M to generalize better while being more parameter efficient. See Appendix [A.4](https://arxiv.org/html/2509.01038v1#A1.SS4) for DynaProt-M additional RMSF plots.

### 4.3 Predicting residue full anisotropy

**Table 3: Comparison of methods on anisotropic blob prediction (ATLAS test split). Runtime for a length 271 protein (7lao_A). 25th %ile / Median reported ($\downarrow$ is better).**
| Method | RMWD Var | Sym. KL Var | # Params | Time |
| --- | --- | --- | --- | --- |
| DynaProt-M | 0.84 /  1.18 | 0.53 /  0.91 | 955 K | $\sim$0.02 s |
| AFMD+T | 0.87 /  1.10 | 0.37 /  0.60 | 95 M | $\sim$7000 s |
| NMA (ANM) | 1.14 /  1.45 | 3.03 /  4.56 | – | $\sim$5.37 s |

To assess the faithfulness of DynaProt-M’s predicted marginals, we compare against both physics-based and learned ensemble methods. In practice, only NMA (ANM) is a feasible baseline, as AlphaFlow is prohibitively slow: A single 271-length protein (7lao_A) requires $\sim$7000 s, compared to $\sim$0.02 s for DynaProt-M (Table [3](https://arxiv.org/html/2509.01038v1#S4.T3)). Note that DynaProt-M predicts this directly, but for AlphaFlow, we first sample 250 structures per protein and then calculate the empirical covariance to define the marginal Gaussians. With NMA, we retrieve the full joint via the normal modes and extract the marginal block diagonals.

To quantify the accuracy, we compute the variance contribution of the symmetric KL divergence (see Appendix [A.2](https://arxiv.org/html/2509.01038v1#A1.SS2)) and the root mean 2-Wasserstein distance (RMWD) as described in Jing et al. (2024), compared to the ground truth marginal Gaussians computed from the ATLAS test set (AFMD split). Despite being orders of magnitude faster and smaller (955k vs. 95M parameters), DynaProt-M achieves competitive accuracy. DynaProt-M attains a median RMWD of 1.18 and symmetric KL divergence of 0.91, both substantially better than NMA (1.45 and 4.56, respectively), and comparable to AFMD+Templates’s 1.10 and 0.60.

Figure: Figure 4: DynaProt marginal Gaussian and residue coupling analysis. A. Renderings of predicted marginal Gaussians compared to ATLAS MD constructed Gaussians (mean symmetric KL divergence and RMWD are reported). B. Joint distribution (within 75th percentile) of DynaProt performance vs. (AFMD+T, NMA). C. Band-wise Pearson correlation between predicted and ground-truth residue–residue coupling matrices as a function of sequence distance.
Refer to caption: x3.png

Moreover, rather than relying solely on summary statistics, we also visualize the distributions (75th percentile) of RMWD and mean symmetric KL variance contributions across test set proteins (Figure [4](https://arxiv.org/html/2509.01038v1#S4.F4)B). These plots compare DynaProt-M to both AFMD+T and NMA on a per-protein basis. Points below the diagonal (highlighted in pink) indicate that DynaProt-M outperformed the method in question on that particular protein. From this, we see that DynaProt-M achieves comparable performance to AFMD+T. Notably, within the 75th percentile it often outperforms AFMD+T (examples visualized in Figure [4](https://arxiv.org/html/2509.01038v1#S4.F4)A,B) on RMWD variance contribution. Moreover, DynaProt-M significantly outperforms NMA across both RMWD and symmetric KL. This further corroborates DynaProt-M’s ability to capture local anisotropic structure well despite being much smaller and faster than other methods.

Figure: Figure 5: DynaProt-M predicted residue Gaussians (ellipsoids) overlaid the apo form.
Refer to caption: x4.png

#### 4.3.1 DynaProt zero-shot cryptic pocket discovery of Adenylosuccinate Synthetase

Beyond accuracy, DynaProt-M’s marginals can also provide functional insight. Many proteins are considered to be undruggable as their apo form may not display a clear binding pocket. However, the druggable pocket may only become apparent after the drug is bound (holo form)–a so called “cryptic pocket." Identification of cryptic pockets is therefore an important task in drug discovery (Mou et al., 2025; Hollingsworth et al., 2019; Comitani and Gervasio, 2018).
As a case study, we sought to investigate DynaProt’s ability in cryptic pocket identification for the enzyme adenylosuccinate synthetase. As indicated in Meller et al., this enzyme is known to exhibit a cryptic pocket and both the apo and holo forms are available in the PDB (1ADE / 1CIB).

We applied DynaProt-M to zero shot predict the marginal Gaussians on the apo form. When we look at the predictions, we notice two features (Figure [5](https://arxiv.org/html/2509.01038v1#S4.F5)). First, the residues with the largest variance are exactly those that encompass the binding pocket.
Second, more than just the magnitude of the Gaussians, the directionality is consistent with a potential pocket opening motion (when compared with the holo form 1CIB). These early results suggest the potential of DynaProt’s utility in cryptic pocket discovery, but a systematic exploration is left for future work.

### 4.4 Predicting scalar coupling

To assess the performance of DynaProt-J in modeling residue–residue coupling, we evaluate its predicted $N\times N$ scalar covariance matrices against those derived from NMA, a classical method for capturing global structural dynamics. While DynaProt-J directly predicts these scalar coupling matrices, we construct comparable matrices for NMA by computing full $3N\times 3N$ anisotropic network models (using ProDy) from each test protein, and then projecting them into $N\times N$ scalar covariances as described in Section [2](https://arxiv.org/html/2509.01038v1#S2.F2). We obtain per-residue correlation matrices by normalizing the entries to be unit diagonal and constrained to the range $[-1,1]$.

We observe that in the ground truth $N\times N$ correlations, magnitude of entries diminish rapidly with distance from main diagonal, indicating weak long-range coupling. This is indicated by the dotted gray line in Figure [4](https://arxiv.org/html/2509.01038v1#S4.F4)C, which shows the mean absolute value of the entries from the principal diagonal up to the $k$th diagonal band. To focus on meaningful and prominent interactions, we define a *diagonal band* of width $k=50$ residues ($|i-j|\leq 50$), which captures local and medium-range interactions. This essentially measures residue-residue coupling as a function of sequence distance (how distal are $i$ and $j$ along the backbone). We compute the Pearson correlation between predicted and ground-truth residue–residue correlation matrices for the entries along each diagonal band $k$, by iteratively extracting the upper-triangular entries satisfying $|i-j|\leq k$ for $k=1$ to 50. This is repeated for each of the 82 test set proteins and the median Pearson correlation is reported (each point in Figure [4](https://arxiv.org/html/2509.01038v1#S4.F4)C). This band-wise analysis enables us to compare the accuracy of coupling signals at increasing residue distances, and we find that DynaProt-J (peak correlation of $r=0.71$) strongly outperforms NMA (peak correlation of $r=0.59$) particularly at short to mid-range coupling distances, where the coupling is the strongest.

### 4.5 Ensemble generation

As described in Section [3.4](https://arxiv.org/html/2509.01038v1#S3.SS4), given the output $3\times 3$ marginal covariances and $N\times N$ residue coupling from DynaProt-M and DynaProt-J respectively, we reconstruct a full joint covariance using the heuristic defined in Eq. [5](https://arxiv.org/html/2509.01038v1#S3.E5). This direct access to the joint distribution enables extremely fast sampling of diverse structures. For evaluation, we sample 250 structures with DynaProt, AFMD+Templates, and NMA to form ensembles for each of the 82 test set proteins in the AFMD split. Following the evaluation protocol of Jing et al. (2024), we assess these ensembles across a variety of metrics including flexibility accuracy, distributional similarity, and the ability to reproduce complex observables. For flexibility accuracy, we measure the pairwise RMSD to ground truth MD and RMSF correlation at the global and per-target level. For distributional coverage, we measure the 2-Wasserstein distance after projecting the ensembles onto the first two principle components derived from the MD trajectory (MD PCA $\mathcal{W}_{2}$) or the combined (MD+sampled) trajectory (Joint PCA $\mathcal{W}_{2}$).

Table [4](https://arxiv.org/html/2509.01038v1#S4.T4) summarizes the ensemble evaluation results across AFMD+Templates, DynaProt, and NMA. DynaProt achieves performance comparable to AFMD+Templates on key flexibility metrics such as pairwise RMSD and per-target RMSF correlation, while lagging slightly behind on distributional similarity and observable recovery. Some examples of where DynaProt outperforms AFMD+Templates on ensemble flexibility correlation are visualized in Figure [6](https://arxiv.org/html/2509.01038v1#S4.F6) and Appendix [A.4](https://arxiv.org/html/2509.01038v1#A1.SS4). For the visuals, all atom reconstruction is enabled by PULCHRA (Rotkiewicz and Skolnick, 2008).
Moreover, DynaProt consistently outperforms NMA across nearly all evaluations—except for transient contact prediction—particularly excelling in measures of local flexibility and pairwise distance preservation.

**Table 4: Comparison of Cα ensemble evaluation metrics on ATLAS MD Dataset between AFMD+Templates, DynaProt, and NMA. Underlined is second best value.**
| Metric | AFMD+Templates | DynaProt | NMA |
| --- | --- | --- | --- |
| Pairwise RMSD (=2.89) | 2.18 | 2.17 | 0.91 |
| RMSF (=1.48) | 1.17 | 1.10 | 0.52 |
| Global RMSF $r$ $(\uparrow)$ | 0.91 | 0.71 | 0.54 |
| Per-target RMSF $r$ $(\uparrow)$ | 0.92 | 0.86 | 0.76 |
| MD PCA $\mathcal{W}_{2}$ $(\downarrow)$ | 1.25 | 1.74 | 1.86 |
| Joint PCA $\mathcal{W}_{2}$ $(\downarrow)$ | 1.58 | 2.39 | 2.45 |
| Weak contacts $J$ $(\uparrow)$ | 0.62 | 0.51 | 0.43 |
| Transient contacts $J$ $(\uparrow)$ | 0.47 | 0.29 | 0.33 |
| # Parameters $(\downarrow)$ | 95 M | 2.86 M | – |
| Ensemble sampling time $(\downarrow)$ | $\sim 10,000$ s | $\sim 0.14$ s | $\sim 5.69$ s |

Notably, DynaProt requires only 2.86 million parameters (vs. 95 million for AFMD+Templates) and samples ensembles over 70,000$\times$ faster on average ($\sim$0.14 s vs. $\sim$10,000 s), all while being trained only to predict marginal and scalar covariances. This efficiency advantage is maintained when compared against sequence-based methods (i.e. ConfDiff, BioEmu, ESMDiff, see Appendix [A.3](https://arxiv.org/html/2509.01038v1#A1.SS3)). DynaProt even outperforms them on modeling flexibility and anisotropy.

As a final assessment of DynaProt’s generalization to longer timescale dynamics, we compare its zero-shot ensemble of BPTI (PDB: 5PTI) to the 1ms trajectory from Shaw et al. (2010). Even with these larger conformational changes, DynaProt performs reasonably well. It achieves RMSF correlation of 0.88 (c.f. 86 on ATLAS), local anisotropy with RMWD of 0.52 Å(c.f. 1.18 Å  on ATLAS), and strong recovery of transient contacts (Jaccard similarity 0.54, c.f. 0.29 on ATLAS). See appendix [A.3.1](https://arxiv.org/html/2509.01038v1#A1.SS3.SSS1) for more.

Figure: Figure 6: Comparison of DynaProt generated ensemble vs. AFMD+T to ATLAS MD simulation (PDB 7qsu_A) overlaid on reference. RMSF Pearson correlation $r$ and sample time reported.
Refer to caption: x5.png

## 5 Conclusion

Protein dynamics is critical for understanding biological function. Existing approaches to modeling dynamics often rely on complex generative models with large-scale PDB pretraining and expensive ensemble generation. In this work, we introduce DynaProt, a lightweight and data-driven alternative akin to Normal Mode Analysis (NMA), but designed to directly predict structured dynamics descriptors in the form of per-residue and pairwise Gaussian representations. This formulation enables extreme parameter efficiency while outperforming traditional baselines on key metrics, including flexibility estimation, marginal anisotropy, and residue–residue coupling. Remarkably, DynaProt’s outputs also support ultra-fast ensemble sampling with reasonable structural fidelity—offering a compelling alternative to conventional ensemble generation methods. While further scaling may be needed to match the full capabilities of state-of-the-art generative methods, our approach highlights a promising alternative grounded in explicitly learning structured representations of dynamics.

## Acknowledgements

We would like to thank Alexander Amini, Jeffrey Li, Yuseong (Nick) Oh, Ryan Mei, Manoj Niverthi, Renzo Soatto, and Kabir Doshi for their useful ideas and discussions.

This work was supported by the National Institute of General Medical Sciences of the National Institutes of Health under award number 1R35GM141861-01.

## References

- Abramson et al. [2024]
Josh Abramson, Jonas Adler, Jack Dunger, Richard Evans, Tim Green, Alexander Pritzel, Olaf Ronneberger, Lindsay Willmore, Andrew J Ballard, Joshua Bambrick, et al.
Accurate structure prediction of biomolecular interactions with alphafold 3.
*Nature*, 630(8016):493–500, 2024.
- Bhatia et al. [2019]
Rajendra Bhatia, Tanvi Jain, and Yongdo Lim.
On the Bures–Wasserstein distance between positive definite matrices.
*Expositiones Mathematicae*, 37(2):165–191, 2019.
- Burley et al. [2017]
Stephen K Burley, Helen M Berman, Gerard J Kleywegt, John L Markley, Haruki Nakamura, and Sameer Velankar.
Protein data bank (pdb): the single global macromolecular structure archive.
*Protein crystallography: methods and protocols*, pages 627–641, 2017.
- Childers and Daggett [2017]
Matthew Carter Childers and Valerie Daggett.
Insights from molecular dynamics simulations for computational protein design.
*Molecular systems design & engineering*, 2(1):9–33, 2017.
- Comitani and Gervasio [2018]
Federico Comitani and Francesco Luigi Gervasio.
Exploring cryptic pockets formation in targets of pharmaceutical interest with swish.
*Journal of chemical theory and computation*, 14(6):3321–3331, 2018.
- Cui and Bahar [2005]
Qiang Cui and Ivet Bahar.
*Normal mode analysis: theory and applications to biological and chemical systems*.
CRC press, 2005.
- Del Alamo et al. [2022]
Diego Del Alamo, Davide Sala, Hassane S Mchaourab, and Jens Meiler.
Sampling alternative conformational states of transporters and receptors with alphafold2.
*Elife*, 11:e75751, 2022.
- Glorot et al. [2011]
Xavier Glorot, Antoine Bordes, and Yoshua Bengio.
Deep sparse rectifier neural networks.
In *Proceedings of the fourteenth international conference on artificial intelligence and statistics*, pages 315–323. JMLR Workshop and Conference Proceedings, 2011.
- Guo et al. [2025]
Amy B Guo, Deniz Akpinaroglu, Christina A Stephens, Michael Grabe, Colin A Smith, Mark JS Kelly, and Tanja Kortemme.
Deep learning–guided design of dynamic proteins.
*Science*, 388(6749):eadr7094, 2025.
- Hollingsworth and Dror [2018]
Scott A Hollingsworth and Ron O Dror.
Molecular dynamics simulation for all.
*Neuron*, 99(6):1129–1143, 2018.
- Hollingsworth et al. [2019]
Scott A Hollingsworth, Brendan Kelly, Celine Valant, Jordan Arthur Michaelis, Olivia Mastromihalis, Geoff Thompson, AJ Venkatakrishnan, Samuel Hertig, Peter J Scammells, Patrick M Sexton, et al.
Cryptic pocket formation underlies allosteric modulator selectivity at muscarinic gpcrs.
*Nature Communications*, 10(1):3289, 2019.
- Hou and Shen [2024]
Chao Hou and Yufeng Shen.
Seqdance: A protein language model for representing protein dynamic properties.
*bioRxiv*, 2024.
- Huang et al. [2015]
Zhiwu Huang, Ruiping Wang, Shiguang Shan, Xianqiu Li, and Xilin Chen.
Log-euclidean metric learning on symmetric positive definite manifold with application to image set classification.
In *International conference on machine learning*, pages 720–729. PMLR, 2015.
- Jeffreys [1998]
Harold Jeffreys.
*The theory of probability*.
OuP Oxford, 1998.
- Jing et al. [2024]
Bowen Jing, Bonnie Berger, and Tommi Jaakkola.
AlphaFold meets flow matching for generating protein ensembles.
*arXiv preprint arXiv:2402.04845*, 2024.
- Jumper et al. [2021]
John Jumper, Richard Evans, Alexander Pritzel, Tim Green, Michael Figurnov, Olaf Ronneberger, Kathryn Tunyasuvunakool, Russ Bates, Augustin Žídek, Anna Potapenko, et al.
Highly accurate protein structure prediction with AlphaFold.
*Nature*, 596(7873):583–589, 2021.
- Kouba et al. [2024]
Petr Kouba, Joan Planas-Iglesias, Jiri Damborsky, Jiri Sedlar, Stanislav Mazurenko, and Josef Sivic.
Learning to engineer protein flexibility.
*arXiv preprint arXiv:2412.18275*, 2024.
- Kullback and Leibler [1951]
Solomon Kullback and Richard A Leibler.
On information and sufficiency.
*The annals of mathematical statistics*, 22(1):79–86, 1951.
- Lewis et al. [2024]
Sarah Lewis, Tim Hempel, José Jiménez-Luna, Michael Gastegger, Yu Xie, Andrew YK Foong, Victor García Satorras, Osama Abdin, Bastiaan S Veeling, Iryna Zaporozhets, et al.
Scalable emulation of protein equilibrium ensembles with generative deep learning.
*bioRxiv*, pages 2024–12, 2024.
- Lombard et al. [2024]
Valentin Lombard, Dan Timsit, Sergei Grudinin, and Elodie Laine.
Seamoon: Prediction of molecular motions based on language models.
*bioRxiv*, pages 2024–09, 2024.
- Lu et al. [2024]
Jiarui Lu, Xiaoyin Chen, Stephen Zhewen Lu, Chence Shi, Hongyu Guo, Yoshua Bengio, and Jian Tang.
Structure language models for protein conformation generation.
*arXiv preprint arXiv:2410.18403*, 2024.
- Ma [2005]
Jianpeng Ma.
Usefulness and limitations of normal mode analysis in modeling dynamics of biomolecular complexes.
*Structure*, 13(3):373–380, 2005.
- McCammon [1984]
JA McCammon.
Protein dynamics.
*Reports on Progress in Physics*, 47(1):1, 1984.
- [24]
Artur Meller, Michael Ward, Jonathan Borowsky, Meghana Kshirsagar, Jeffrey M. Lotthammer, Felipe Oviedo, Juan Lavista Ferres, and Gregory R. Bowman.
Predicting locations of cryptic pockets from single protein structures using the PocketMiner graph neural network.
14(1):1177.
ISSN 2041-1723.
[10.1038/s41467-023-36699-3](https:/doi.org/10.1038/s41467-023-36699-3).
URL [https://doi.org/10.1038/s41467-023-36699-3](https://doi.org/10.1038/s41467-023-36699-3).
- Mou et al. [2025]
Minyue Mou, Weicheng Yang, Guangyi Huang, Xiaoyan Yang, Xiao Zhang, Wasala Mudiyanselage Wishwajith Wickramabahu Kandegama, Charles R Ashby Jr, Gefei Hao, and Yangyang Gao.
The discovery of cryptic pockets increases the druggability of “undruggable” proteins.
*Medicinal Research Reviews*, 2025.
- Rotkiewicz and Skolnick [2008]
Piotr Rotkiewicz and Jeffrey Skolnick.
Fast procedure for reconstruction of full-atom protein models from reduced representations.
*Journal of computational chemistry*, 29(9):1460–1465, 2008.
- Shaw et al. [2010]
David E Shaw, Paul Maragakis, Kresten Lindorff-Larsen, Stefano Piana, Ron O Dror, Michael P Eastwood, Joseph A Bank, John M Jumper, John K Salmon, Yibing Shan, et al.
Atomic-level characterization of the structural dynamics of proteins.
*Science*, 330(6002):341–346, 2010.
- Skjaerven et al. [2009]
Lars Skjaerven, Siv M Hollup, and Nathalie Reuter.
Normal mode analysis for proteins.
*Journal of Molecular Structure: THEOCHEM*, 898(1-3):42–48, 2009.
- Stein and Mchaourab [2022]
Richard A Stein and Hassane S Mchaourab.
Speach_af: Sampling protein ensembles and conformational heterogeneity with alphafold2.
*PLOS Computational Biology*, 18(8):e1010483, 2022.
- Vander Meersche et al. [2024]
Yann Vander Meersche, Gabriel Cretin, Aria Gheeraert, Jean-Christophe Gelly, and Tatiana Galochkina.
ATLAS: protein flexibility description from atomistic molecular dynamics simulations.
*Nucleic Acids Research*, 52(D1):D384–D392, 2024.
- Vemulapalli and Jacobs [2015]
Raviteja Vemulapalli and David W Jacobs.
Riemannian metric learning for symmetric positive definite matrices.
*arXiv preprint arXiv:1501.02393*, 2015.
- Wang et al. [2024]
Yan Wang, Lihao Wang, Yuning Shen, Yiqun Wang, Huizhuo Yuan, Yue Wu, and Quanquan Gu.
Protein conformation generation via force-guided se (3) diffusion models.
*arXiv preprint arXiv:2403.14088*, 2024.
- Wayment-Steele et al. [2024]
Hannah K Wayment-Steele, Adedolapo Ojoawo, Renee Otten, Julia M Apitz, Warintra Pitsawong, Marc Hömberger, Sergey Ovchinnikov, Lucy Colwell, and Dorothee Kern.
Predicting multiple conformations via sequence clustering and alphafold2.
*Nature*, 625(7996):832–839, 2024.
- Wayment-Steele et al. [2025]
Hannah K Wayment-Steele, Gina El Nesr, Ramith Hettiarachchi, Hasindu Kariyawasam, Sergey Ovchinnikov, and Dorothee Kern.
Learning millisecond protein dynamics from what is missing in nmr spectra.
*bioRxiv*, pages 2025–03, 2025.
- Yu and Koshland Jr [2001]
Edward W Yu and Daniel E Koshland Jr.
Propagating conformational changes over long (and short) distances in proteins.
*Proceedings of the National Academy of Sciences*, 98(17):9517–9520, 2001.
- Zhang et al. [2024]
Mingyang Zhang, Ting Chen, Xun Lu, Xiaobing Lan, Ziqiang Chen, and Shaoyong Lu.
G protein-coupled receptors (gpcrs): advances in structures, mechanisms and drug discovery.
*Signal Transduction and Targeted Therapy*, 9(1):88, 2024.
- Zhang et al. [2021]
She Zhang, James M Krieger, Yan Zhang, Cihan Kaya, Burak Kaynak, Karolina Mikulska-Ruminska, Pemra Doruker, Hongchun Li, and Ivet Bahar.
Prody 2.0: increased scale and scope after 10 years of protein dynamics modelling with python.
*Bioinformatics*, 37(20):3657–3659, 2021.
- Zinovjev et al. [2024]
Kirill Zinovjev, Paul Guénon, Carlos A Ramos-Guzmán, J Javier Ruiz-Pernía, Damien Laage, and Iñaki Tuñón.
Activation and friction in enzymatic loop opening and closing dynamics.
*Nature Communications*, 15(1):2490, 2024.

## Appendix A Appendix

### A.1 Method details

##### SPD closure of joint reconstruction heuristic (restating Proposition 3.1).

Given marginal covariances $\{\bm{\Sigma}_{\text{marginal}}^{(i)}\in\mathbb{R}^{3\times 3}\}_{i=1}^{N}$ and a correlation matrix $\widetilde{\bm{C}}\in\mathbb{R}^{N\times N}$ that is symmetric and positive definite, then the reconstructed joint covariance
$\bm{\Sigma}_{\text{joint}}=\bm{L}_{\text{marginal}}\left(\widetilde{\bm{C}}\otimes\bm{I}_{3}\right)\bm{L}_{\text{marginal}}^{\top}$
is also symmetric and positive definite.

###### Proof A.1 .

Let $\bm{L}_{\text{marginal}}\in\mathbb{R}^{3N\times 3N}$ be the block-diagonal matrix defined as

$$ $\bm{L}_{\text{marginal}}=\bigoplus_{i=1}^{N}\bm{L}_{i},$ $$

where each $\bm{L}_{i}\in\mathbb{R}^{3\times 3}$ is the Cholesky factor (or any valid matrix square root) of the positive definite matrix $\bm{\Sigma}_{\text{marginal}}^{(i)}$. Since each $\bm{\Sigma}_{\text{marginal}}^{(i)}\succ 0$, it follows that each $\bm{L}_{i}$ is full rank, lower triangular, and has positive diagonal entries. Consequently, $\bm{L}_{\text{marginal}}$ is full rank and lower triangular with positive diagonal blocks.

Now consider the matrix $\widetilde{\bm{C}}\otimes\bm{I}_{3}\in\mathbb{R}^{3N\times 3N}$. Since $\widetilde{\bm{C}}\succ 0$ and $\bm{I}_{3}\succ 0$, the Kronecker product $\widetilde{\bm{C}}\otimes\bm{I}_{3}\succ 0$ as well (Kronecker product of two SPD matrices is also SPD). Finally, the product

$$ $\bm{\Sigma}_{\text{joint}}=\bm{L}_{\text{marginal}}\left(\widetilde{\bm{C}}\otimes\bm{I}_{3}\right)\bm{L}_{\text{marginal}}^{\top}$ $$

is a congruence transformation of the SPD matrix $\widetilde{\bm{C}}\otimes\bm{I}_{3}$ by the full-rank matrix $\bm{L}_{\text{marginal}}$. Since congruence preserves positive definiteness, we conclude:

$$ $\bm{\Sigma}_{\text{joint}}\succ 0$ $$

Moreover, $\bm{\Sigma}_{\text{joint}}$ is symmetric because it is of the form $\bm{A}\bm{B}\bm{A}^{\top}$.

###### Proof A.1 .

##### Multivariate Gaussian Sampling.

Let $\bm{\epsilon}\sim\mathcal{N}(\bm{0},\bm{I})$ be a standard multivariate normal in $\mathbb{R}^{d}$, and let $\bm{\mu}\in\mathbb{R}^{d}$, $\bm{\Sigma}\in\mathbb{R}^{d\times d}$ be a symmetric positive definite matrix. Suppose $\bm{L}\in\mathbb{R}^{d\times d}$ satisfies $\bm{\Sigma}=\bm{L}\bm{L}^{\top}$ (e.g., via Cholesky decomposition or matrix square root). Then the transformation $\bm{x}=\bm{\mu}+\bm{L}\bm{\epsilon}$ yields a random variable $\bm{x}\sim\mathcal{N}(\bm{\mu},\bm{\Sigma})$.

###### Proof A.2 .

Since Gaussian distributions are fully characterized by their first two cumulants (mean and covariance), it suffices to show that the transformed variable has the desired mean and covariance.

Mean of $\bm{x}$:

$$ $\mathbb{E}[\bm{x}]=\mathbb{E}[\bm{\mu}+\bm{L}\bm{\epsilon}]=\bm{\mu}+\bm{L}\cdot\mathbb{E}[\bm{\epsilon}]=\bm{\mu}$ $$

Covariance of $\bm{x}$:

$$ $\displaystyle\operatorname{Cov}[\bm{x}]$ $\displaystyle=\mathbb{E}\left[(\bm{x}-\bm{\mu})(\bm{x}-\bm{\mu})^{\top}\right]$ $\displaystyle=\mathbb{E}\left[(\bm{L}\bm{\epsilon})(\bm{L}\bm{\epsilon})^{\top}\right]$ $\displaystyle=\mathbb{E}\left[\bm{L}\bm{\epsilon}\bm{\epsilon}^{\top}\bm{L}^{\top}\right]$ $\displaystyle=\bm{L}\cdot\mathbb{E}[\bm{\epsilon}\bm{\epsilon}^{\top}]\cdot\bm{L}^{\top}$ $\displaystyle=\bm{L}\cdot\bm{I}_{d}\cdot\bm{L}^{\top}$ $\displaystyle=\bm{L}\bm{L}^{\top}=\bm{\Sigma}$ $$

###### Proof A.2 .

### A.2 Evaluation Metrics

##### RMWD Variance Contribution.

To evaluate the efficacy of the marginal Gaussian predictions, we adopt distributional similarity metrics used in Jing et al. [2024].
The first of these is the root mean 2-Wasserstein distance (RMWD), specifically its variance contribution term.
The 2-Wasserstein distance between two multivariate Gaussians has a closed-form expression as follows.

Let $\mathcal{N}_{0}=\mathcal{N}(\bm{\mu}_{0},\bm{\Sigma}_{0})$ and $\mathcal{N}_{1}=\mathcal{N}(\bm{\mu}_{1},\bm{\Sigma}_{1})$ be two $d$-dimensional Gaussian distributions. The squared 2-Wasserstein distance between them is given by:

$$ $\displaystyle\mathcal{W}_{2}^{2}(\mathcal{N}_{0},\mathcal{N}_{1})=\|\bm{\mu}_{0}-\bm{\mu}_{1}\|_{2}^{2}+\operatorname{Tr}\left(\bm{\Sigma}_{0}+\bm{\Sigma}_{1}-2\left(\bm{\Sigma}_{1}^{1/2}\bm{\Sigma}_{0}\bm{\Sigma}_{1}^{1/2}\right)^{1/2}\right)$ (6) $$

This expression consists of two additive components: a mean contribution and a covariance (variance) contribution.
This metric is also referred to as the Bures–Wasserstein distance [Bhatia et al., 2019].
Since our method predicts only the covariances, we isolate and evaluate only the second term.
We define the RMWD variance contribution across $N$ residues as follows:

$$ $\displaystyle\text{RMWD}_{\text{var}}(\mathcal{N}_{0},\mathcal{N}_{1})=\sqrt{\frac{1}{N}\sum_{i=1}^{N}\operatorname{Tr}\left(\bm{\Sigma}_{0,i}+\bm{\Sigma}_{1,i}-2\left(\bm{\Sigma}_{1,i}^{1/2}\bm{\Sigma}_{0,i}\bm{\Sigma}_{1,i}^{1/2}\right)^{1/2}\right)}$ (7) $$

##### Symmetric KL Divergence Variance Contribution.

Alongside the Wasserstein-based metric, we also evaluate the discrepancy between predicted and ground-truth marginal distributions using the symmetric Kullback–Leibler (KL) divergence, defined as the mean of the two directed KL divergences mentioned in [Kullback and Leibler, 1951, Jeffreys, 1998]:

$$ $\displaystyle\mathrm{KL}_{\text{sym}}(\mathcal{N}_{0}\,\|\,\mathcal{N}_{1})=\frac{1}{2}\left[\mathrm{KL}(\mathcal{N}_{0}\,\|\,\mathcal{N}_{1})+\mathrm{KL}(\mathcal{N}_{1}\,\|\,\mathcal{N}_{0})\right]$ $$

For two $d$-dimensional Gaussian distributions $\mathcal{N}_{0}=\mathcal{N}(\bm{\mu}_{0},\bm{\Sigma}_{0})$ and $\mathcal{N}_{1}=\mathcal{N}(\bm{\mu}_{1},\bm{\Sigma}_{1})$, the KL divergence from $\mathcal{N}_{0}$ to $\mathcal{N}_{1}$ is given by:

$$ $\displaystyle\mathrm{KL}(\mathcal{N}_{0}\,\|\,\mathcal{N}_{1})$ $\displaystyle=\frac{1}{2}\left[\operatorname{Tr}(\bm{\Sigma}_{1}^{-1}\bm{\Sigma}_{0})+(\bm{\mu}_{1}-\bm{\mu}_{0})^{\top}\bm{\Sigma}_{1}^{-1}(\bm{\mu}_{1}-\bm{\mu}_{0})-d+\log\frac{\det\bm{\Sigma}_{1}}{\det\bm{\Sigma}_{0}}\right]$ (8) $$

This expression consists of both a mean contribution—the Mahalanobis term—and a variance contribution, comprising the trace and log-determinant terms.
Since our method predicts only covariances (and uses the input structure coordinates as means), we isolate the variance terms by omitting ($\mu_{1}=\mu_{0}$) the mean term:

$$ $\displaystyle\mathrm{KL}_{\text{var}}(\mathcal{N}_{0}\,\|\,\mathcal{N}_{1})$ $\displaystyle=\frac{1}{2}\left(\operatorname{Tr}(\bm{\Sigma}_{1}^{-1}\bm{\Sigma}_{0})-d+\log\frac{\det\bm{\Sigma}_{1}}{\det\bm{\Sigma}_{0}}\right)$ (9) $$

To symmetrize the variance contribution of the divergence, we define the symmetric variance KL as:

$$ $\displaystyle\mathrm{KL}_{\text{symvar}}(\mathcal{N}_{0},\mathcal{N}_{1})$ $\displaystyle=\frac{1}{2}\left(\mathrm{KL}_{\text{var}}(\mathcal{N}_{0}\,\|\,\mathcal{N}_{1})+\mathrm{KL}_{\text{var}}(\mathcal{N}_{1}\,\|\,\mathcal{N}_{0})\right)$ (10) $\displaystyle=\frac{1}{4}\left(\operatorname{Tr}(\bm{\Sigma}_{1}^{-1}\bm{\Sigma}_{0})+\operatorname{Tr}(\bm{\Sigma}_{0}^{-1}\bm{\Sigma}_{1})-2d\right)$ (11) $$

### A.3 Additional Experiments

#### A.3.1 DynaProt zero-shot ensemble generation of BPTI

We note that DynaProt was trained on the ATLAS MD dataset comprising 100 ns per replicate trajectories. In contrast, D.E. Shaw Research performed simulations of BPTI (PDB: 5PTI) at millisecond-scale revealing structurally distinct conformational states [Shaw et al., 2010]. Thus, in an effort to understand DynaProt’s ability to generalize to long-timescale dynamics, we applied it to BPTI and compared to the DESRES trajectory.

Listed in Table [7(b)](https://arxiv.org/html/2509.01038v1#A1.F7.sf2), we compute the ensemble evaluation metrics from Jing et al. [2024] and observe that DynaProt performs remarkably well: e.g., RMSF correlation of 0.88 (c.f. 86 on ATLAS), local anisotropy with RMWD of 0.52 Å  (c.f. 1.18 Å  on ATLAS), and strong recovery of transient contacts (Jaccard similarity 0.54, c.f. 0.29 on ATLAS). These metrics emphasize that DynaProt is able to model larger conformational changes at high fidelity.

Figure: (a) Visualization of DynaProt zero shot BPTI (PDB: 5PTI) ensemble generation.
Refer to caption: x6.png

#### A.3.2 DynaProt ablations

**Table 5: DynaProt-M ablations of the log Frobenius loss loss and SE(3) invariance.**
| Metric | DynaProt | No LogFrob<br>Loss | No SE(3)<br>Invariance |
| --- | --- | --- | --- |
| RMWD Var (↓) | 1.18 | 2.70 | 1.92 |
| Sym KL Var (↓) | 0.91 | 9.26 | 4.46 |
| RMSF $r$ (↑) | 0.87 | 0.38 | 0.48 |

To test both the importance of DynaProt’s Riemannian aware loss (log Frobenius norm) and the SE(3) invariance from the IPA layers, we have performed the following ablations listed in Table [5](https://arxiv.org/html/2509.01038v1#A1.T5). Unsurprisingly, replacing the log Frobenius norm objective with standard Mean Squared Error loss significantly degrades performance as the optimization is over the space of positive definite covariance matrices, which lies on a well-studied Riemannian manifold. Replacing the IPA blocks with standard MLPs also degrades performance, suggesting that SE(3) invariance is crucial in this low-data, low-parameter regime.

#### A.3.3 Sequence input baselines

To further evaluate DynaProt’s ensemble generation, we compare against methods that take sequence as input following the AlphaFlow evaluations. DynaProt outperforms these methods on local RMSF correlation and marginal anisotropy prediction and is comparable with other distributional metrics. Moreover, the efficiency advantage is clear with DynaProt’s sub-second sample time.

Table: Table 6: Comparison of DynaProt generated with ensemble generation methods that take in sequence as input. ESMDiff, ESM3 entries reported from Lu et al. [2024].

### A.4 Supplementary Figures

Figure: Figure 8: DynaProt-M predicted RMSF correlations. Visualized test set examples of predicted RMSF per residue (derived from the predicted marginal Gaussians) compared to ground truth RMSF derived from MD trajectories. Pearson correlation coefficient ($r$) between predicted and ground truth RMSF is reported.
Refer to caption: x7.png

## Peer Reviews

confidence: 3
contribution: 3
initial_rating: 6
presentation: 4
questions: ### 1. Relationship between DynaProt-M, DynaProt-J, and DynaProt

In the Experiments “preprocessing” section, DynaProt-M, DynaProt-J, and the IPA backbone appear to be trained as three separate models. However, in Section 4.4, the final DynaProt results are described as being composed from the output of DynaProt-M and  DynaProt-J. Could the authors clarify which is the case? And if separate models are trained, why is it necessary (vs. a single model with two prediction heads)?
    
### 2. Cryptic pocket case study and sampling
In Section 4.2.1, the authors show that predicted marginal variances can indicate ligand-binding dynamics, such as cryptic pocket formation. Have the authors sampled from the full reconstructed Gaussian distribution for this protein, and verified whether both open and closed conformations can actually be generated?
    
Additionally, BioEmu (https://www.science.org/doi/10.1126/science.adv9817) includes several benchmarks with cryptic pocket motions. Evaluating DynaProt on those systems and reporting pocket coverage metrics would further strengthen the results.
    
### 3. Extending to atomic-level or protein-ligand complexes
    
DynaProt demonstrates impressive speed and potential scalability for residue-level modeling. Is it feasible to extend this approach to atomic-level resolution or to explicitly model protein-ligand complexes? Are there technical challenges such as dimensionality of the covariance matrix or parameterizing atom-level couplings?
soundness: 3
strengths: - The manuscript is well organized, clearly written, and easy to follow, with clean and informative figures that effectively support the main ideas.
- The model addresses a unique and practical perspective in conformational modeling, where not only accuracy but also the ability to efficiently generate large numbers of samples is essential. DynaProt improves upon classical normal mode analysis while being orders of magnitude faster than modern generative models.
- The method produces versatile outputs (e.g., RMSF, residue-level anisotropic Gaussian profiles, residue-residue coupling) and achieves competitive accuracy on some ensemble-level metrics compared to current state-of-the-art models.
summary: DynaProt is a Gaussian-based framework for ultra-fast protein conformation ensemble generation. Unlike traditional generative models that learn protein structure distributions via sampling, DynaProt learns a hierarchical Gaussian parameterization composed of residue-level anisotropic Gaussians and a protein-level isotropic residue–residue covariance matrix that captures coupling between residues.
This formulation enables DynaProt to analytically compute dynamics-related quantities such as RMSF without sampling; reconstruct the full 3N×3N covariance matrix using simple heuristics to support conformation sampling, generate large conformational ensembles extremely efficiently by sampling from reconstructed Gaussian distribution.
Compared with classical normal mode analysis, which also relies on Gaussian approximations, DynaProt learns variances and couplings directly from data rather than through simplified physical potentials, resulting in more accurate ensemble statistics.
Experimental results show that DynaProt achieves competitive performance on some ensemble-level metrics compared to diffusion generative models while being orders of magnitude faster. Although some accuracy is sacrificed due to the inherent Gaussian assumptions, the method offers a compelling balance between speed, accuracy, and physical interpretability, making it a strong and distinctive solution for conformation sampling and protein dynamics analysis.
weaknesses: - The core limitation lies in the Gaussian approximation and reliance on input structures, which restricts expressiveness compared to more flexible generative models - but this is aligned with the paper’s goal of trading some accuracy for major gains in speed and efficiency.
- Given this design philosophy, clearly demonstrating the cases where DynaProt is useful becomes important. Cryptic pocket modeling is a relevant application, but only one case study is provided.
- The current formulation models only coarse-grained residue-level dynamics. It would strengthen the paper if the authors could demonstrate potential extensions to atomic-resolution modeling or protein-ligand complexes.

---

confidence: 4
contribution: 2
initial_rating: 2
presentation: 3
questions: The same as weaknesses. In multi-conformation generation tasks, more attention is paid to the changes in system transitions, (rather than the local fluctuation). If the author is willing to explain this point, I am willing to increase the score
soundness: 3
strengths: 1. The paper is well organized and clearly stated.
2. The authors propose a protein dynamics framework based on a Gaussian hierarchical representation, which incorporates per-residue covariance tensors and the global coupling matrix into a Gaussian hierarchy.
3. Compared with a single scalar RMSF, the covariance matrix is richer: it captures both the shape and direction of fluctuations.
summary: The paper presents DYNAPROT, an innovative and lightweight SE(3)-invariant model that predicts protein dynamical features—marginal anisotropy and residue-level couplings—directly from static structures. It achieves accuracy comparable to leading generative models without large-scale PDB pretraining, while being more efficient. By constructing covariance matrices and introducing scalar covariances, the model can capture some key features, which are useful for fast ensemble generation.
weaknesses: 1. The core issue is model misspecification. A Gaussian hierarchy simplifies learning and sampling, and it can not reflect the thermodynamics of MD data. Molecular energy landscapes are multi-well with distinct basins. Gaussian models are locally valid only near a single basin (approximately harmonic vibrations) and cannot capture complex, anharmonic shapes or multimodality across multiple minima. This defect was masked in the atlas dataset. Proteins in Atlas are relatively large, so the existing MD trajectories predominantly capture local fluctuations rather than informative transitions. Could the authors validate the method on simpler systems—e.g., a 2D double-well potential and fast-folding proteins—to demonstrate effectiveness of model?
2. Although Eq. 5 and Proposition 3.1 guarantee positive semidefiniteness, the construction in Eq. 5 lacks clear physical motivation. The marginal covariance $\Sigma_{marginal}$ rightly emphasizes per-residue anisotropy, which is a strength. However, the coupling matrix $C$ enforces isotropic interactions—i.e., identical influence along all three spatial directions—between residues. This design oversimplifies residue–residue coupling and risks weakening key directional effects that are critical in complex biomolecular systems. Could the authors elaborate on the motivation and theoretical basis for this design choice?
3. When constructing $C$, the use of **mean pooling** as a scalar projection to compute each $C_{ij}$ appears casual and overly simplistic for modeling residue–residue coupling. A more completed ablation is needed to justify this design choice.
4. What is the $\mu$ of the BPTI results in Appendix A.3.1? Are the two main conformational changes caused by $\mu$ or by $\Sigma$?

---

confidence: 4
contribution: 2
initial_rating: 4
presentation: 3
questions: The model was trained on 100ns of data. Beyond the single example of BPTI, have the authors tested the model on additional proteins exhibiting known long-timescale (e.g. > 1us) dynamic characteristics? Is there a possibility that the model might overfit to short-timescale fluctuations, thereby underestimating genuine large-scale motions?
soundness: 2
strengths: Compared to models such as DYNAPROT (with 2.86 million parameters in total) and AlphaFlow (95 million parameters), the parameter count is reduced by several orders of magnitude. Its inference and sampling speeds far exceed those of other methods, rendering it exceptionally practical. Moreover, the model performs well and is trained solely on molecular dynamics data from approximately 1,000 proteins.
summary: This paper introduces DYNAPROT, a novel, lightweight, SE(3)-invariant deep learning framework for directly predicting protein dynamics from static protein structures. Unlike molecular dynamics (MD) simulations, which demand substantial computational resources, or generative models requiring extensive pre-training and costly sampling procedures, DYNAPROT is an efficient alternative.
weaknesses: The model was trained on the ATLAS dataset comprising 100 ns simulations. For numerous biologically significant processes—such as large-scale domain motions and allosteric transitions—100 ns is wholly inadequate, as these typically unfold over timescales ranging from microseconds to milliseconds or longer. Consequently, the model may have predominantly learned local, transient thermal fluctuations while exhibiting limited capacity to capture large-scale conformational changes. Although the paper's zero-shot testing results for BPTI appear favourable, this represents a single protein case study. One cannot conclude from this solitary example that the model generalises well to all long-timescale dynamics.

---

confidence: 4
contribution: 3
initial_rating: 8
presentation: 4
questions: 1. When sampling from the normal distribution of protein structures predicted by the model, how often do unphysical protein structures appear?  What prevents this approach from producing impossible torsion angles or bond lengths?  

2. How does the model handle proteins with pronounced conformational changes, such as large domain rotations, extensive motion cascades, or highly distinct alternate conformations like those observed in 4OLE?

3. How does the model behave when the provided protein is not in its lowest-energy state?  How would the model perform as part of a dynamical system, where a protein conformation is iteratively updated along the principal directions of each residue? Would the structure return toward the mean state, or would it break down entirely?
soundness: 4
strengths: This paper addresses protein dynamics prediction not implicitly through ensemble generation, but directly by modeling the protein’s residue C$\alpha$ covariance matrix.  This approach enables direct analysis of protein dynamics.  
Modeling the full covariance matrix using a heuristic that combines the cross-residue covariance matrix with the intra-residue position matrix provides an efficient and meaningful representation of dynamics.  
The experiments presented are diverse and interesting, demonstrating the versatility of the proposed approach.
summary: This paper proposes a pipeline for training a model to predict protein dynamics.  The authors model protein dynamics as sampling conformations from a normal distribution with a low-energy mean and a non-identity covariance matrix.  The overall covariance is represented as a combination of a local residue-wise 3×3 covariance matrix and a global residue-wise N×N covariance matrix, matching the dynamics learned from MD simulations.
weaknesses: 1. Throughout the paper, the examples presented focus on local dynamics, where the protein does not deviate significantly from the mean. Proteins that undergo large conformational changes may not be well represented by modeling their dynamics as a normal distribution with a single covariance matrix. For proteins that open and close or contain highly flexible regions, the full dynamic range might collapse into a single blurred representation under this parameterization. Nevertheless, this work represents an important step toward explicit modeling of protein dynamics.  

2. The authors assume that the provided protein structure is in a low-energy state, but there is no discussion of how the model behaves when the input structure is not in such a state.


Decision: Accept (Poster)