Molecular enhanced sampling with autoencoders: On-the-fly collective variable discovery and accelerated free energy landscape exploration

Wei Chen; Andrew L Ferguson

Source-linked AI summary

Molecular enhanced sampling with autoencoders: On-the-fly collective variable discovery and accelerated free energy landscape exploration

Wei Chen, Andrew L Ferguson

arXiv:1801.00203v2Dec 30, 2017physics.bio-phphysics.comp-phstat.ML

TL;DR

Biomolecular simulations require effective collective variables to overcome high free-energy barriers, but nonlinear CV-discovery methods often lack the explicit coordinate mapping needed for direct biasing. The paper introduces MESA, an autoencoder-based approach that learns explicit, differentiable nonlinear CVs and biases sampling along them. In Trp-cage, MESA achieved more comprehensive configurational sampling and an approximately 20 kBT free-energy range, while the comparison sampled only the native basin over approximately 7 kBT.

Problem
Effective biased sampling requires good collective variables and explicit coordinate mappings, but existing nonlinear CV-discovery approaches typically do not provide the mapping needed for direct biasing.
Method
MESA uses autoencoders to discover nonlinear collective variables that are explicit and differentiable functions of atomic coordinates, enabling direct biasing along the learned coordinates.
Results
MESA achieved more comprehensive configurational sampling and an approximately 20 kBT free-energy range, compared with approximately 7 kBT and only the native basin for the comparison.
Takeaways & Limitations
The approach simultaneously discovers data-driven CVs and accelerates sampling directly in those coordinates.
Takeaways & Limitations
The force field prevented observation of unfolding of the Trp-cage N-terminal α-helix.

Abstract

from arXiv · show

Macromolecular and biomolecular folding landscapes typically contain high free energy barriers that impede efficient sampling of configurational space by standard molecular dynamics simulation. Biased sampling can artificially drive the simulation along pre-specified collective variables (CVs), but success depends critically on the availability of good CVs associated with the important collective dynamical motions. Nonlinear machine learning techniques can identify such CVs but typically do not furnish an explicit relationship with the atomic coordinates necessary to perform biased sampling. In this work, we employ auto-associative artificial neural networks ("autoencoders") to learn nonlinear CVs that are explicit and differentiable functions of the atomic coordinates. Our approach offers substantial speedups in exploration of configurational space, and is distinguished from exiting approaches by its capacity to simultaneously discover and directly accelerate along data-driven CVs. We demonstrate the approach in simulations of alanine dipeptide and Trp-cage, and have developed an open-source and freely-available implementation within OpenMM.

1 Introduction

Enhanced sampling of biomolecular folding is limited by slow motions, high free-energy barriers, and incomplete exploration, while effective CV biasing depends on difficult-to-identify collective variables. MESA uses autoencoders to discover nonlinear CVs that can be used directly for biasing, addressing the mapping problem faced by existing nonlinear approaches.

Sampling challenge: High free-energy barriers and slow conformational motions trap molecular dynamics simulations in incompletely sampled regions of phase space.These barriers arise from rough free-energy surfaces and limit access to thermally accessible configurations.
Sampling challenge: CV biasing can be more efficient than tempering when good CVs track important motions, but poor or unknown CVs can make it fail badly.Selecting suitable CVs is normally difficult except for simple systems.
Sampling challenge: Data-driven CV discovery faces a chicken-and-egg problem because good sampling is needed both to discover CVs and to use them for efficient sampling.Existing workflows therefore interleave accelerated sampling with repeated CV discovery.
Limits of prior approaches: Nonlinear manifold-learning methods can discover useful CVs but generally do not provide the atomic-coordinate mapping required to propagate biasing forces.Without an explicit differentiable mapping, direct biased molecular dynamics in the learned CVs is not possible.
MESA approach: MESA uses autoencoders to discover nonlinear CVs as explicit, differentiable functions of atomic coordinates and bias sampling directly along them.The CVs can be implemented in biasing methods such as umbrella sampling or metadynamics.
MESA approach: MESA is introduced as Molecular Enhanced Sampling with Autoencoders for simultaneous data-driven CV discovery and accelerated sampling.The approach is presented as distinct from existing nonlinear machine-learning methods because it supports direct acceleration in the discovered coordinates.

2.1 Nonlinear dimensionality reduction using autoencoders

Autoencoders learn low-dimensional nonlinear projections and approximate reconstructions of molecular data. Their bottleneck outputs define explicit, differentiable collective variables for molecular configurations.

Architecture: Autoencoders project high-dimensional molecular data into a low-dimensional bottleneck and approximately reconstruct the original inputs.The architecture is a symmetric five-layer D-G-K-G-D network, with K specifying the intrinsic dimensionality.
Training: Training the network to reproduce its inputs discovers a bottleneck projection that captures salient features of the data.Parameters are optimized through reconstruction error, with early stopping used instead of weight decay for regularization.
Differentiable mapping: Smooth hyperbolic-tangent activations make the learned projection analytically differentiable with respect to its inputs.The input layer uses identity mappings, while the remaining layers use nonlinear activations with analytical first derivatives.
Dimensionality determination: The bottleneck dimensionality K is selected by training candidate architectures and locating a knee or plateau in the fraction of variance explained.The authors use cross-validation for G and the L-method to automate knee detection; K ≪ D with FVE above 70% supports a low-dimensional representation.
Molecular collective variables: The bottleneck outputs define collective variables, and the projection function explicitly maps high-dimensional inputs to those variables.For molecular applications, the inputs are Cartesian coordinates, with rotational and translational changes treated separately from internal structural reconfigurations.

2.2 Enhanced sampling in autoencoder CVs

The learned autoencoder variables can be used directly in biased molecular dynamics because their coordinate mappings provide analytical biasing forces. Umbrella sampling over these variables is combined with WHAM to estimate unbiased free-energy surfaces.

Biasing in learned variables: Autoencoders address a key limitation of nonlinear dimensionality reduction by furnishing an explicit mapping needed to bias directly in discovered collective variables.The mapping from high-dimensional data to CVs is used to propagate biasing forces to atomic coordinates.
Alternative biasing methods: The method can also use metadynamics or adaptive biasing force calculations after discovering the autoencoder CVs.The paper uses umbrella sampling, while identifying these other popular CV-biasing approaches as straightforward alternatives.
Umbrella sampling: Umbrella sampling augments the atomic Hamiltonian with restraining potentials defined on the learned CVs.Overlapping umbrella windows can collectively explore CV space and surmount free-energy barriers when the CVs track slow molecular motions and distinguish metastable states.
Force propagation: Biasing forces on the learned CVs are propagated to atomic forces through repeated application of the chain rule and analytical Jacobians.The final Jacobian acts as a filter, passing forces only to atoms whose coordinates are included in the autoencoder inputs.
Free-energy estimation: WHAM combines biased histograms from umbrella windows to estimate the unbiased probability distribution and corresponding free-energy surface.Free energies in other collective variables can also be estimated from the biased data and WHAM solution.

2.3 MESA: Interleaved on-the-fly CV discovery and enhanced sampling

MESA interleaves autoencoder-based CV discovery with boundary-targeted umbrella sampling, enabling explicit, differentiable CVs to guide successive exploration of configurational space.

MESA framework: MESA discovers explicit, differentiable collective variables from molecular coordinates and uses them directly for enhanced sampling.The framework combines nonlinear CV discovery with biasing in the learned coordinates.
Iterative protocol: The six-step protocol iteratively trains autoencoders, detects explored-region boundaries, launches umbrella runs, and assesses convergence.New configurations are incorporated into subsequent CV discovery and sampling rounds.
Computational cost: Autoencoder training is much cheaper than biased simulation: approximately 5 GPU-minutes versus 8 GPU-hours for the reported Trp-cage workload.The comparison used 32,000 harvested configurations and 15 × 2 ns biased simulations.
Boundary detection: Boundary detection identifies both the frontier and internal holes of the explored intrinsic manifold across shapes and dimensions.The grid-based procedure targets sparsely populated cells adjacent to dense regions.
Enhanced sampling: Enhanced sampling in identified CV-space cells advances the exploration frontier and fills internal holes through direct biasing.Harmonic umbrella potentials are centered in selected boundary cells, with force constants tuned to the system and local landscape.
Protocol benefits: MESA converges toward global data-driven CVs that parameterize accessible configurational space and provide interpretable slow collective motions.The explicit coordinate-to-CV mapping supports both direct biasing and physical interpretation.

2.4 Molecular dynamics simulations

The study simulates alanine dipeptide and Trp-cage under distinct molecular-dynamics conditions, then applies umbrella sampling through an OpenMM biasing-force plugin.

Alanine dipeptide: Alanine dipeptide umbrella calculations used 10–20 windows per round, harmonic restraints of κ=3000 kJ/mol·(unit of CV)^2, and 100 ps per window.Configurations were saved every 1 ps.
Performance: The simulations achieved approximately 6 µs/day·core for unbiased alanine dipeptide and 1.3 µs/day·core for biased calculations.These execution speeds were reported for Intel i7-5820K CPUs.

2.5 Implementation and open source availability

The authors released the MESA implementation and its OpenMM biasing-force components as free, open-source software for CPU and CUDA platforms.

OpenMM integration: Molecular simulations used the open-source OpenMM 7.0 package.The ANN Force plugin computes Jacobian elements from trained autoencoders and derives atomic biasing forces.
Software stack: Neural networks were constructed and trained with Keras running on Theano.The implementation used the Keras deep-learning Python API.
Released software: The enhanced-sampling framework automates autoencoder training, boundary detection, biased umbrella runs, and post-processing.The framework is available from the authors’ public repository.
Released software: ANN Force was implemented for both CPU and CUDA OpenMM platforms.The plugin is available under the authors’ public GitHub repository.

3 Results and Discussion

MESA rapidly expands sampling in alanine dipeptide and Trp-cage, discovers low-dimensional collective variables, and produces free-energy landscapes that agree with direct sampling or prior work in the reported comparisons.

Alanine dipeptide: The initial 800 ps alanine-dipeptide trajectory explored C5 and C7 but made no transition across the high barriers to αL.This illustrates the sampling difficulty that motivates MESA.
Alanine dipeptide: MESA reached the αL basin by iteration 6 and converged within 10 iterations, with CVs and the explored region stabilizing by iterations 9–10.The iteration-10 biasing run barely advanced the exploration frontier relative to iteration 9.
Alanine dipeptide: The MESA free-energy surface for alanine dipeptide reproduced the three metastable states but misidentified the αL well and introduced a spurious minimum.The artifacts were attributed to periodicity-related issues in the learned representation.
Alanine dipeptide: MESA identified four CVs for alanine dipeptide’s flat-torus embedding, avoiding projection problems from periodic φ and ψ variables.The projected four-dimensional free-energy surface agreed excellently with direct sampling in φ-ψ space.
Trp-cage: The two Trp-cage CVs correlated with global chirality (ρ=-0.92) and end-to-end distance (ρ=-0.91), identifying handedness and peptide extent.Both correlations had two-tailed p-values below 10^-15.
Trp-cage: The Trp-cage landscape resolved metastable states and a low-free-energy folding pathway consistent with several prior studies.The authors report good agreement between the free-energy surface, folding pathways, and earlier replica-exchange, transition-path, and related studies.
Trp-cage: The present Trp-cage calculations did not observe unfolding of the N-terminal α-helix, unlike one prior study’s alternative pathway.The authors relate the preserved helical content to favorable unfolded-state helical stability and note consistency with spectroscopy.
Trp-cage: For Trp-cage, MESA converged within 6 iterations and expanded the sampled free-energy range to approximately 20 kBT versus approximately 7 kBT for the unbiased run.Umbrella sampling barely advanced the frontier between iterations 5 and 6.

4 Conclusions

MESA combines iterative autoencoder-based nonlinear CV discovery with direct biased sampling in those coordinates. It recovers important free-energy features and folding pathways in alanine dipeptide and Trp-cage while substantially expanding sampled configurational space.

Method: MESA introduces data-driven accelerated order-parameter discovery using autoencoders to learn low-dimensional nonlinear collective variables.The framework is designed for accelerated sampling of macromolecular and biomolecular configurational space.
Method: Because the learned CVs are explicit and differentiable functions of atomic coordinates, biasing potentials can be converted into forces on atoms.This enables direct accelerated sampling in the learned nonlinear CVs rather than relying on proxy variables or basis projections.
Method: MESA iterates autoencoder CV discovery and biased sampling until the CVs and explored configurational region converge, requiring only an initial short trajectory to seed discovery.The framework is largely automated and requires minimal user intervention.
Validation: MESA recovered converged free-energy surfaces resolving important metastable states and folding pathways for alanine dipeptide and Trp-cage.For Trp-cage, the recovered free-energy surface and diffusion-collision folding mechanism agreed well with prior work.
Limitations: The approach remains limited by manifold periodicities, force-field-dependent sampling, network-training burden, training-data generation, regularization, and the lack of CV ranking.For Trp-cage, the force field prevented observation of unfolding of the N-terminal α-helix and the nucleation-condensation pathway was not resolved.

Appendix

The appendix evaluates the intrinsic-manifold dimensionality across MESA iterations for alanine dipeptide and Trp-cage, using autoencoder reconstruction variance and knee detection. Both systems identify a two-dimensional manifold, with Trp-cage maintaining K=2 across all six iterations.

Appendix figures: Figures A1 and A2 show the intrinsic-manifold dimensionality determination across iterations for the two systems.The appendix presents the iteration-wise dimensionality analysis for alanine dipeptide and Trp-cage.
Dimensionality analysis: Five independently initialized autoencoders are trained for each bottleneck size and iteration to obtain the FVE distributions.Boxplots summarize quartiles and whisker ranges, with outliers shown separately.
Alanine dipeptide: K=2 is determined for alanine dipeptide using the knee in the fraction-of-variance-explained analysis.The knee is identified with the L-method applied to piecewise linear fits.
Trp-cage: The intrinsic-manifold dimensionality is determined to be K=2 for Trp-cage in all six MESA iterations.The result is reported for MESA applied to Trp-cage in water.