PaperlayerView exact revision on arXiv

Source-linked AI summary

Deep learning the slow modes for rare events sampling

Luigi Bonati, GiovanniMaria Piccini, Michele Parrinello

arXiv:2107.03943v2physics.comp-phphysics.chem-ph

TL;DR

The paper addresses how to identify slow collective variables without prior long-term dynamics. It combines neural-network variational learning with OPES reweighting to extract slow modes from biased simulations and accelerate them. The method increases transition sampling and rapidly converges free-energy estimates in the demonstrated alanine-dipeptide case.

  • Problem

    Transfer-operator eigenfunctions provide attractive collective variables, but determining them ordinarily requires long-term dynamics that are generally unavailable.

  • Method

    Deep-TICA uses a neural-network ansatz to extract transfer-operator eigenfunctions from reweighted enhanced-sampling trajectories, then accelerates the leading mode with OPES.

  • Results

    200-fold increase in transitions per unit time was obtained for alanine dipeptide, with the free-energy difference converging within 0.1 kBT in 1 ns.

  • Takeaways & Limitations

    Deep-TICA can extract and accelerate slow modes from generalized-ensemble simulations or simulations using approximate collective variables.

Abstract

from arXiv · show

The development of enhanced sampling methods has greatly extended the scope of atomistic simulations, allowing long-time phenomena to be studied with accessible computational resources. Many such methods rely on the identification of an appropriate set of collective variables. These are meant to describe the system's modes that most slowly approach equilibrium. Once identified, the equilibration of these modes is accelerated by the enhanced sampling method of choice. An attractive way of determining the collective variables is to relate them to the eigenfunctions and eigenvalues of the transfer operator. Unfortunately, this requires knowing the long-term dynamics of the system beforehand, which is generally not available. However, we have recently shown that it is indeed possible to determine efficient collective variables starting from biased simulations. In this paper, we bring the power of machine learning and the efficiency of the recently developed on-the-fly probability enhanced sampling method to bear on this approach. The result is a powerful and robust algorithm that, given an initial enhanced sampling simulation performed with trial collective variables or generalized ensembles, extracts transfer operator eigenfunctions using a neural network ansatz and then accelerates them to promote sampling of rare events. To illustrate the generality of this approach we apply it to several systems, ranging from the conformational transition of a small molecule to the folding of a mini-protein and the study of materials crystallization.

Time lagged independent component analysis

The method replaces fixed TICA descriptors with neural-network-learned basis functions and applies variational analysis to reweighted enhanced-sampling trajectories. An OPES-based workflow then extracts and accelerates the slow modes as Deep-TICA collective variables.

  • TICA basis functions: TICA selects linear combinations of descriptors with maximal autocorrelation, while the variational formulation becomes a generalized eigenvalue problem.Mean-zero functions are constrained to be orthogonal to the trivial constant solution.
  • Neural-network ansatz: Neural networks learn nonlinear basis functions from descriptors in a lower-dimensional space, extending the variational approach to many descriptors.This increases the flexibility of the trial functions beyond a predetermined linear descriptor basis.
  • Neural-network ansatz: Time-lagged descriptor pairs are mapped to latent variables, whose mean-free covariance matrices yield eigenvalues and eigenfunctions optimized by maximizing the leading eigenvalues.The loss minimized by gradient descent corresponds to the VAMP-2 score.
  • Enhanced-sampling data: Reweighted enhanced-sampling trajectories can be analyzed in scaled time t′, where the VAC procedure is applied to recover equilibrium correlation functions.The scaled time is defined through the bias potential, and the resulting transfer-operator spectrum reflects the accelerated dynamics of the initial simulation.
  • Enhanced-sampling data: OPES estimates P(s) on the fly and biases the system toward a chosen target distribution, including well-tempered, uniform, or multithermal ensembles.At convergence, the free-energy surface is obtained from F(s) = −kBT log P(s).
  • Protocol: The recommended workflow explores transitions with trial-CV or generalized-ensemble OPES, trains Deep-TICA on the resulting trajectories, then biases the leading eigenfunction.The initial bias V*(s0) is retained in the Hamiltonian during subsequent sampling because the learned modes reflect convergence under that biased dynamics.

Results and discussion

Deep-TICA extracts slow collective variables from limited or poor exploratory simulations and uses them to accelerate rare-event sampling across molecular folding, conformational transitions, and crystallization. The resulting simulations show faster transitions and accurate free-energy estimates without requiring prior system-specific understanding.

  • Alanine dipeptide: Deep-TICA recovered efficient alanine dipeptide variables from limited exploratory sampling, remained robust to lag time and configuration count, and promoted both transition pathways.Sampling focused on minima and connecting transition regions.
  • Alanine dipeptide: 200-fold more transitions per unit time reduced alanine dipeptide inter-state intervals from 3 ns to about 25 ps.The free-energy difference converged within 0.1 kBT in 1 ns.
  • Alanine dipeptide: A poor ψ-driven exploratory simulation was followed by diffusive Deep-TICA dynamics and free-energy convergence in 1 ns.The authors present this as evidence that substantial speedups can result from identifying and accelerating slow modes.
  • Chignolin folding: For chignolin, Deep-TICA 1 described the slow folded–unfolded transition, while Deep-TICA 2 captured fine structure within the folded state.The analysis found no evidence of stable misfolded states along the dominant folding variables.
  • Chignolin folding: A Deep-TICA 1 bias increased chignolin folding-event rates 20-fold and produced free-energy profiles with about 0.5 kJ/mol average statistical error.At 340 K, the profiles agreed excellently with a 106 µs unbiased reference trajectory.
  • Silicon crystallization: For silicon crystallization, Deep-TICA increased solid–liquid transitions and enabled free-energy convergence after 20 ns, while its smoother transition representation aligned with crystalline-atom fraction.Silicon crystallization remains difficult because directional bonding permits defective and glassy structures.

Conclusions

The paper develops Deep-TICA as a general protocol that analyzes biased trajectories, extracts slow modes, and accelerates them for rare-event sampling.

  • Deep-TICA combines variational analysis of enhanced-sampling data, neural networks, and advanced sampling to construct a general protocol.
  • The method extracts slowly converging modes from biased trajectories and subsequently accelerates them.
  • Deep-TICA can analyze generalized-ensemble simulations and complement physically or data-driven approximate collective variables.
  • The alanine dipeptide test shows applicability even when the initial enhanced-sampling simulation uses a very poor collective variable.

Materials and methods

The study trains Deep-TICA collective variables from reweighted enhanced-sampling trajectories and uses them to bias subsequent simulations. It evaluates this workflow across alanine dipeptide, chignolin, and silicon crystallization simulations.

  • Time-lagged covariance matrices: Time-reweighted configuration pairs at lag time τ are used to construct symmetrized time-lagged covariance matrices, although τ is not a physical time.The matrices are symmetrized to enforce detailed balance, introducing a stated bias.
  • Deep-TICA CVs training: Deep-TICA CVs are trained in PyTorch by converting the covariance generalized eigenvalue problem into a standard problem through Cholesky decomposition.This formulation permits gradient backpropagation through the eigenvalue problem using automatic differentiation.
  • Deep-TICA CVs training: The neural network uses two feed-forward layers with hyperbolic-tangent activations and ADAM optimization at a learning rate of 1e-3.Training uses standardized inputs, training/validation splitting, and early stopping with a patience of 10 epochs; Deep-TICA outputs are scaled to [-1,1].
  • PLUMED-Pytorch interface: A modified PLUMED2–LibTorch interface loads Python-trained models to evaluate collective variables and derivatives and apply bias potentials during simulations.This connects neural-network CV training with enhanced-sampling production runs.
  • Alanine dipeptide and chignolin simulations: Alanine dipeptide and chignolin are simulated with OPES-based protocols, using multithermal ensembles or trial CVs before Deep-TICA-based biasing.The protocols use 300–600 K for alanine dipeptide and 270–700 K with eight shared-bias replicas for chignolin.
  • Silicon simulations: Silicon crystallization simulations use a 216-atom 3x3x3 supercell with the Stillinger–Weber potential, while structure-factor peaks initialize a Deep-LDA CV for OPES sampling.The protocol includes solid- and liquid-state standard-MD simulations followed by a 50 ns OPES run, excluding its first 25 ns from neural-network training.

Deep-TICA training

Deep-TICA training examines how lag time and data volume affect extraction of slow eigenfunctions. The first eigenfunction remains consistent across lag times, while useful approximations can emerge after only a few transitions.

  • Lag-time dependence: The first Deep-TICA eigenfunction remains consistent across the studied lag times, whereas the second loses signal as its eigenvalue becomes too small.The analysis therefore favors lag times for which all desired eigenvalues remain nonzero.
  • Configuration dependence: Already after a couple of transitions, Deep-TICA extracts a very good approximation of the eigenfunctions.Training configurations were sampled every 1 ps after 15 ns in the alanine dipeptide multithermal example.

Multithermal simulation - 2D free energy

The multithermal simulation produces a two-dimensional free-energy representation in Deep-TICA coordinates. Deep-TICA 1 captures the C7eq-to-C7ax transition, while Deep-TICA 2 resolves two substates within C7eq.

  • Free-energy landscape: Deep-TICA 1 describes the conformational transition between C7eq and C7ax.The central free-energy panel is plotted in the two Deep-TICA CVs, with one-dimensional projections shown along the top and left.
  • Free-energy landscape: Deep-TICA 2 highlights two substates within C7eq.This distinguishes internal structure within one of the conformational basins.

Multithermal simulation - Convergence and 1D free energy profiles

Convergence analyses compare exploratory and Deep-TICA-based enhanced-sampling simulations using free-energy profiles and transition-rate measurements. The free-energy difference is defined from integrals over the C7eq and C7ax regions of the φ-based profile.

  • Convergence analysis: The convergence analyses track the free-energy difference between C7eq and C7ax over time and report one-dimensional free-energy profiles.Figure S4 compares multithermal exploratory sampling, Deep-TICA 1 biasing in the multithermal ensemble, and OPES biasing along Deep-TICA 1.
  • Transition-rate analysis: Average transition rates are obtained by dividing simulation time by the number of transitions detected from the running average of Deep-TICA 1.Transitions are counted using a 2 ps window and only quasi-static-bias portions are included.
  • Free-energy definition: The free-energy difference is computed between regions A = C7eq and B = C7ax, corresponding here to φ < 0 and φ > 0.The calculation uses the free-energy profile F(s) with s = φ.
  • Convergence analysis: Figure S5 compares ψ-biased OPES, Deep-TICA 1 with static V*(ψ), and OPES biasing only Deep-TICA 1.The same free-energy-difference and one-dimensional-profile diagnostics are used.

Deep-TICA CVs comparison

Deep-TICA CV isolines are compared in the Ramachandran plane using weights from multithermal and ψ-based OPES simulations. The analysis reports isolines for Deep-TICA 1 and Deep-TICA 2 under a uniform-sampling histogram.

  • Simulation comparison: The two columns compare Deep-TICA CV isolines weighted from the multithermal simulation and the ψ-based OPES simulation.Configurations come from an OPES simulation biasing φ-ψ toward a flat target distribution.
  • CV comparison: The two rows show isolines for Deep-TICA 1 and Deep-TICA 2 in the Ramachandran plane.The isolines are computed from a two-dimensional weighted histogram.

Replicas trajectories

Figure S7 compares Cα-RMSD trajectories across replicas in OPES Multi-T and OPES Multi-T* + Deep-TICA 1 simulations, distinguishing replicas with shifted times and dashed boundaries.

  • Figure S7 overlays Cα-RMSD evolution for OPES Multi-T and OPES Multi-T* + Deep-TICA 1 across system replicas.Each color represents a replica sharing the same bias potential.

Improving sampling of the transition region

The figures examine energy, Deep-TICA coordinates, and structural free energies across temperatures to assess transition-region sampling and compare it with reference behavior.

  • Figure S8 compares potential energy versus Deep-TICA 1 between OPES Multi-T and OPES Multi-T* + Deep-TICA 1.Points are colored by Cα-RMSD, and the energy range corresponds to temperatures from 280K to 500K.
  • Figure S9 compares 2D free energies in Deep-TICA coordinates at T=340K with a reference unbiased free-energy surface.It also evaluates projection onto Deep-TICA CVs trained using all heavy-atom distances rather than a limited descriptor subset.
  • Figures S5–S6-derived analyses use a single OPES Multi-T* + Deep-TICA 1 simulation to examine temperature-dependent free energies.The supplied captions specify temperature-dependent FES analyses in backbone hydrogen-bond/end-to-end coordinates and Deep-TICA coordinates.
  • The hydrogen-bond and end-to-end descriptors are defined using continuous switching functions and terminal-residue Cα distances, respectively.Hydrogen bonds use specified distance and switching-function parameters; end-to-end distance uses the two terminal tyrosine Cα atoms.
  • Figure S12 reports temperature-dependent free-energy profiles along Deep-TICA 1, including statistical uncertainty and a 340K reference profile.Deep-TICA 1 describes the folding–unfolding transition.

Characterization of the folded states

The folded ensemble contains multiple structurally distinct states that share common hydrogen bonds but differ in side-chain interactions and torsional configurations.

  • Three folded states are separated using weighted k-means clustering with reweighting and a Deep-TICA 1>0.65 cutoff.Backbone descriptors alone do not discriminate among the three folded states.
  • The folded states share common hydrogen bonds but differ in specific side-chain hydrogen bonds.The most populated state includes a hydrogen bond between the alcohol oxygens of the two threonines.
  • State 2 further divides into two states with different equilibrium values of the threonine side-chain torsional angles.The dihedral-angle analysis supports an ensemble rather than a single folded structure.
  • Relative folded-state populations are estimated by integrating P(s)=e^-βF(s) within basins defined along Deep-TICA 2.The analysis uses the condition Deep-TICA 1>0.65 when integrating the 2D free-energy surface.
  • Figure S11 compares free-energy profiles for Deep-LDA and Deep-TICA simulations, including free-energy differences over time.Shaded areas represent statistical uncertainties estimated with weighted block averages.
Loading 2107.03943v2…