All-Atom GPCR-Ligand Simulation via Residual Isometric Latent Flow

Abstract: G-protein-coupled receptors (GPCRs), primary targets for over one-third of approved therapeutics, rely on intricate conformational transitions to transduce signals. While Molecular Dynamics (MD) is essential for elucidating this transduction process, particularly within ligand-bound complexes, conventional all-atom MD simulation is computationally prohibitive. In this paper, we introduce GPCRLMD, a deep generative framework for efficient all-atom GPCR-ligand simulation.GPCRLMD employs a Harmonic-Prior Variational Autoencoder (HP-VAE) to first map the complex into a regularized isometric latent space, preserving geometric topology via physics-informed constraints. Within this latent space, a Residual Latent Flow samples evolution trajectories, which are subsequently decoded back to atomic coordinates. By capturing temporal dynamics via relative displacements anchored to the initial structure, this residual mechanism effectively decouples static topology from dynamic fluctuations. Experimental results demonstrate that GPCRLMD achieves state-of-the-art performance in GPCR-ligand dynamics simulation, faithfully reproducing thermodynamic observables and critical ligand-receptor interactions.

• The paper introduces GPCRLMD, a generative ML framework for all-atom GPCR-ligand trajectory simulation, achieving state-of-the-art performance across multiple benchmarks.
• GPCRLMD employs a harmonic-prior VAE and residual latent flow matching, anchored to initial conformations, yielding efficient (2x speedup) and stable trajectory generation.
• Benchmarks show GPCRLMD outperforms BioMD and BioEmu in flexibility, electronic observables, and reconstruction fidelity, providing a framework for efficient and realistic simulation of GPCR-ligand interactions.

Motivation and problem setting

G protein-coupled receptors (GPCRs) are the target of more than one-third of approved drugs, and their signaling function depends on ligand-induced conformational transitions across the seven-transmembrane architecture. All-atom molecular dynamics (MD) remains the reference method for characterizing these transitions in ligand-bound complexes, but the sampling cost of step-wise numerical integration makes comprehensive ensemble generation impractical, even with enhanced-sampling techniques. The authors motivate the problem with an empirical observation: on PDB 6LW5, the root-mean-square fluctuation (RMSF) profile of the receptor differs substantially between apo and ligand-bound states (Pearson $r = 0.45$), so ligand-free conformational models are insufficient for the bound regime.

Existing generative approaches leave a gap. Protein-centric ensemble models (e.g., BioEmu) ignore non-protein components entirely; ligand-focused methods such as NeuralMD and UnbindingFlow do not resolve fine-grained protein fluctuations; BioMD attempts all-atom protein-ligand simulation but operates directly in Cartesian coordinates without GPCR-specific priors. The paper introduces GPCRLMD, a deep generative framework for trajectory forecasting—predicting the temporal evolution of a GPCR-ligand complex given an initial conformation—and claims state-of-the-art performance for all-atom GPCR-ligand dynamics.

Benchmark curation

The training corpus derives from GPCRmd, comprising 639 MD trajectories (3 replicas × 213 distinct complexes), each spanning 500 ns at 2500 frames. A strict split by receptor sequence identity (<50%) yields 531/57/51 trajectories for train/validation/test, corresponding to over 1.32 million distinct conformational states used in training. Training uses strided window sampling: sub-trajectories of 50 frames with stride 10 cover 100 ns of physical time per sample. Evaluation follows the protocol of Jing et al., covering flexibility correlation (pairwise RMSD, global and per-target RMSF Pearson correlations), distributional accuracy (root-mean 2-Wasserstein distance and its translation/variance components, PCA-based $W_2$ distances, PC-similarity fraction), and ensemble observables (Jaccard similarity of weak/transient contacts and exposed residues, Spearman correlation of exposed mutual-information matrices).

Harmonic-prior VAE

The first component replaces the standard Gaussian VAE prior $\mathcal{N}(\mathbf{0}, \mathbf{I})$ with an atom-centered harmonic prior $\mathcal{N}(\mathbf{x}_i, \sigma\mathbf{I})$, anchored to each atom's instantaneous coordinate $\mathbf{x}_i$. The KL term thus penalizes deviation of latent codes from atomic positions, which the authors justify on three grounds: it mirrors the harmonic approximation of potential wells used in normal-mode analysis; it aligns latent noise with thermal RMSF-like fluctuations; and it prevents posterior collapse while preserving explicit topology. Crucially, the latent space is isometric—the latent retains the same dimensionality as the input ($\mathbf{z}_t \in \mathbb{R}^{N \times 3}$)—so no compression bottleneck discards structural detail, contrary to most latent-space generative approaches.

The encoder is a DiT-based atom attention module that treats receptor residues and ligand atoms as one contiguous sequence, fusing allosterically coupled information into atom-level latents. The decoder is hierarchical: token-level temporal attention captures global dependencies across frames, then atom-level temporal attention reconstructs fine-grained coordinates conditioned on these tokens. Domain knowledge enters as conditioning: learnable embeddings of transmembrane helix indices (TM1–TM7) per residue and of ligand pharmacological class (agonist, antagonist, inverse agonist, unknown). Reconstruction uses rigidly aligned weighted MSE with ligand atoms upweighted ($\lambda_{\text{ligand}} = 10$), plus auxiliary losses covering ligand geometric centering, steric collision penalties, ligand bond lengths, side-chain bond/angle/torsion terms via cosine-parameterized angular losses, backbone torsions, and a smooth differentiable LDDT surrogate. An energy-guided extension performs test-time gradient minimization of a side-chain geometric potential after flow sampling, acting as a proximal operator toward chemically valid geometries.

Residual latent flow matching

The second component learns temporal evolution via rectified flow matching, but on the residual displacement field rather than absolute coordinates. Defining $$\mathbf{r}_{t,i} = \mathbf{z}_{t,i} - \mathbf{x}_{0,i}$$ relative to the initial frame, the flow transports isotropic Gaussian noise to the target displacement over $\tau \in [0,1]$, with straight-path interpolation and a velocity network regressing the clean residual. Because latents remain three-dimensional and anchored to valid geometry, generated samples stay within the VAE manifold's valid region, and model capacity concentrates on transition densities rather than reconstructing trivial covalent structure. Sampling integrates an ODE with only 10 Euler steps, producing all $T-1$ future frames in parallel; long trajectories are obtained autoregressively by chaining the final frame as the next initial condition. A Stochastic Token Propagator mixes initial-frame token context with time-scaled noise to bridge the information gap to later time steps.

Results

Against BioMD retrained on the same data, GPCRLMD reports consistent gains across 100–500 ns horizons: pairwise RMSD correlation 0.73 vs. 0.58 at 100 ns, global RMSF correlation 0.84 vs. 0.74, weak-contact Jaccard similarity 0.54 vs. 0.11, and exposed-residue Jaccard 0.54 vs. 0.31, while requiring 2.0 GPU seconds versus 10 per sample. Against BioEmu using official checkpoints (evaluated on receptor coordinates only), the gap is sharper: BioEmu achieves pairwise RMSD correlation of −0.10 versus GPCRLMD's 0.73, and requires roughly 1800 GPU seconds per sample—a reported 900× speedup. The authors interpret the BioEmu comparison as evidence that ligand-blind baselines fail to capture binding-induced conformational rearrangements, though this comparison involves a model not designed for the task and trained on different data, so it isolates the value of explicit ligand modeling rather than architectural superiority alone.

Qualitatively, the model reproduces the TM2–TM6 inter-helix distance dynamics that hallmark GPCR activation, maintains bounded ligand RMSD consistent with stable binding, recovers dual-cutoff contact maps, and preserves the conserved Asp$W_2$0 salt bridge when present in the initial frame (and correctly propagates its absence when not). TICA projection analysis shows the generated trajectories populate low-energy basins defined by held-out replicas, achieving lower Jensen-Shannon divergence than both BioMD and, on three of four metrics, a held-out 500 ns MD replica itself—an indication of broader conformational coverage than any single replica. Dihedral coupling distributions ($W_2$1-$W_2$2 backbone, $W_2$3-$W_2$4 side chain) and ligand torsion distributions closely track ground-truth statistics.

Ablations attribute performance to each component. Removing the residual anchoring is the most damaging change: weak-contact Jaccard collapses from 0.53 to 0.02 and the joint PCA $W_2$5 distance rises from 2.05 to 3.40. Replacing the latent flow with direct VAE sampling degrades flexibility correlations (0.84 → 0.76 global RMSF). The harmonic prior matters mainly for long-horizon generation: at 500 ns its removal raises the root-mean Wasserstein distance from 2.42 to 2.94, whereas at short horizons the difference is marginal or slightly favorable to the ablated variant—suggesting the dynamic prior primarily stabilizes chained autoregressive rollout rather than single-window fidelity.

Limitations and open questions

The paper concedes several constraints. Ligand chemical diversity in GPCRmd-derived training data is limited despite the million-plus conformations, raising questions about generalization to chemically novel ligands. Because all predictions are anchored to $W_2$6, performance depends on the physical plausibility of that starting frame; the authors suggest energy relaxation before inference but do not evaluate it here. The salt-bridge result also implies the model cannot form interactions absent from the initial structure, bounding its ability to discover new binding poses. Whether the harmonic prior's benefit persists beyond 500 ns horizons, whether the residual formulation transfers to systems with large induced-fit rearrangements far from $W_2$7, and whether the framework extends beyond GPCRs to general protein-ligand systems remain open questions explicitly deferred to future work.

Conclusion

GPCRLMD demonstrates that all-atom simulation of GPCR-ligand complexes benefits from two coupled design choices: an isometric latent space regularized by physics-informed harmonic priors, and flow matching on residuals anchored to the initial conformation. The approach challenges the prevailing reliance on dimensionality reduction in latent generative modeling, delivering second-scale generation of 100–500 ns ensembles with fidelity to thermodynamic observables, transmembrane helix mechanics, and receptor-ligand contact patterns. Its principal unresolved risks—limited ligand diversity and sensitivity to the initial frame—are clearly identified and define concrete next steps for validating the method's generality.