A self-learning algorithm for biased molecular dynamics

Introduction

• Barriers between low-energy states restrict short molecular-dynamics simulations to configurations near the input structure.The visited region is described as small and locally ergodic.Source: Introduction · p. 1
• Metadynamics generates bias from the previously visited trajectory, among approaches that bias selected collective variables to enhance sampling [14–22].The paper describes metadynamics as introduced by the authors [23,24] and applied to chemical problems [25].Source: Introduction · p. 1
• Choosing collective variables that span reactive motion is difficult, and bias volume in metadynamics grows exponentially with dimensionality.Barriers in transverse degrees of freedom can also leave sampling incomplete.Source: Introduction · p. 1
• Reconnaissance Metadynamics targets frequently occupied low-free-energy regions with a self-learning bias and is demonstrated on a Lennard-Jones cluster and a short protein.The authors implement the algorithm in PLUMED.Source: Introduction · p. 1

Background

• Variable transformations can reduce barriers for many collective coordinates, but the cited method does not work when hydrogen bonds are present [30, 31].Higher-temperature collective variables avoid that particular problem, but do not prevent revisiting configurations [11–13, 32].Source: Background · p. 1
• Principal components or path-based collective variables can reduce dimensionality, but path methods require insight to choose an informative initial path [33–35].The cited path method uses distance along and from the path as two collective variables.Source: Background · pp. 1–2
• Zero-temperature minimum-based methods use local basin modes as descriptors, while GAMUS fits visited configurations with a Gaussian mixture and adaptively biases kinetic traps [36, 37, 29].GAMUS speeds filling over conventional metadynamics with three or four collective variables, but filling time still grows exponentially with their number.Source: Background · p. 2
• The Fig. 2 schematic contrasts one PCA analyzer for a single basin with two for hopping between sub-basins, and illustrates spurious basins when hills are added or basin sizes fixed.The spurious-basin problem worsens as dimensionality increases.Source: Figure 2 · p. 2

Reconnaissance metadynamics algorithm

• A patchwork of basin-specific, locally valid one-dimensional coordinates lets the bias compensate depressions rapidly even when the total number of collective variables is large.Basins are identified dynamically by periodically analyzing the trajectory.Source: Reconnaissance metadynamics algorithm · p. 2
• The algorithm’s flow chart summarizes its sequence, with component details provided in the Methods and supplementary information.The figure is identified as a flow chart for reconnaissance metadynamics.Source: Figure 1 · p. 2, Reconnaissance metadynamics algorithm · p. 3
• Regular clustering of stored trajectory configurations uses combined dimensionality reduction and clustering because PCA alone cannot separate hopped basins and Gaussian-mixture EM can become unstable at high dimensionality [29, 38–40].The clustering returns basin centers and covariance matrices; low-weight or previously encountered similar basins can be discarded.Source: Reconnaissance metadynamics algorithm · p. 2
• Assuming spherical symmetry in the covariance-induced metric, the method adds small Gaussian hills along a radial coordinate to avoid incomplete compensation that can create spurious low-energy features.A single hill may fail to compensate a basin fully, with the resulting spurious features becoming more severe as dimensionality increases.Source: Reconnaissance metadynamics algorithm · pp. 2–3
• Each radial hill acts on a hyper-ellipsoidal crust shaped by the covariance matrix, so bias follows free-energy anisotropy and basin filling time no longer grows exponentially with coordinate count.The crust-like hills increase the integrated bias as radial distance increases.Source: Reconnaissance metadynamics algorithm · p. 3
• Basins expand when the system reaches their rim outside other basins’ influence, with a probabilistic criterion that makes expansion harder as simulation progresses.The expansion probability uses a user-defined parameter D and the time between checks.Source: Reconnaissance metadynamics algorithm · p. 3, Figure 3 · p. 3

Results

• The model demonstration used Langevin dynamics on a two-dimensional potential-energy surface.The setup is introduced as an illustration of algorithm operation.Source: Results · p. 3
• After bias was added, the particle quickly escaped its initial basin and rapidly explored the entire model energy surface as newly identified basins were compensated.The trajectory began trapped at low temperature before metadynamics was applied.Source: Results · p. 3
• The algorithm found only basins corresponding to true free-energy features and used small hills for sub-basins before larger hills compensated for the surrounding super-basin.The adaptive bias is shown handling small basins encompassed by larger depressions.Source: Results · p. 3

• LJ7 has four identified potential-energy minima and 19 saddle points, and at kBT = 0.1ϵ it usually oscillates around its symmetric central-atom minimum, with infrequent isomerizations.Source: Lennard-Jones 7 · p. 3
• Across seven coordination-number slices, basin positions and sizes remain consistent, suggesting the algorithm recovers LJ7’s interchange symmetry and explores configuration space correctly.Source: Figure 4 · p. 4, Lennard-Jones 7 · p. 4
• Because local CVs make the bias lack a straightforward connection to free energy, free-energy estimates require umbrella sampling; trajectories still provide basin landmarks for validating low-dimensional descriptions.Source: Lennard-Jones 7 · p. 4
• LJ7 basin centers cluster around minima on the well-tempered free-energy surface, with second and third coordination-number moments providing optimal in-plane separation.Source: Lennard-Jones 7 · p. 4, Figure 4 · p. 4
• For polyalanine-12, 24 periodic backbone dihedrals served as collective coordinates, with von Mises distributions and periodic differences and averages accounting for their periodicity.Source: Lennard-Jones 7 · p. 4
• After approximately 22 ns, polyalanine-12 folds into its global-minimum alpha helix, whereas a 1-µs conventional MD run from a random configuration did not fold.Minimizations every 20 ps show the helix is considerably lower in energy than the other encountered minima.Source: Lennard-Jones 7 · pp. 4–5, Figure 5 · p. 5

Conclusions

• Reconnaissance metadynamics patches locally valid one-dimensional coordinates into a bias, yielding a compact hierarchical landscape description that can inform umbrella-sampling bias potentials.The method identifies low-free-energy features dynamically and has no limit on the number of collective coordinates.Source: Conclusions · p. 5
• Across two demonstrated systems, trajectories extensively explore low-energy configuration space; the Lennard-Jones landscape also admits an effective two-dimensional description.Generic data such as torsional angles and coordination numbers can be used, including coordinates without large barriers.Source: Conclusions · p. 5

Mixture of probabilistic principle component analyzers.

• Clustering uses annealing and selects the cluster count by comparing fits with different counts using the largest Bayesian Information Criterion value.The annealing parameter σ^2 starts at the data covariance’s maximum eigenvalue and is lowered below 1% of its initial value.Source: Mixture of probabilistic principle… · pp. 5–6

Selecting novel basins.

• New basins are selected when their fit weight, discounted by similarity to existing basins, exceeds a user-defined tolerance.Similarity is measured with Matusita’s measure; expanded basins have their covariance scaled to account for their expanded volume.Source: Selecting novel basins · p. 6

1.1 Expectation-Maximization

• Expectation-maximization alternates point-probability estimation with parameter updates, improving likelihood toward a maximum-likelihood solution.The E-step calculates point probabilities under the current model; the M-step uses them to update its parameters.Source: §1.1 · p. 8

1.2 Gaussian mixture models

• Gaussian mixtures fit component centers, covariances, and weights with expectation maximization, but parameter counts grow quadratically with dimension and initialization can reach local maxima.The dimensionality burden limits their use to low-dimensional data; multiple starting fits are one way to address initialization sensitivity.Source: §1.2 · p. 9

1.3 Deterministic annealing

• Deterministic annealing lowers the Gaussian standard deviation as a temperature, reducing local-maximum trapping and starting from a single cluster at high temperature.As σ falls, clusters can split at successive phase transitions to reflect the data’s clustered structure.Source: §1.3 · p. 10

1.4 Probabilistic principle component analysis

• Ordinary PCA finds directions of maximal variation in the sampled region, but those directions are not guaranteed to remain valid in unexplored configuration space.Probabilistic PCA instead fits both principal components and the region over which they are valid.Source: §1.4 · p. 10

1.5 Mixtures of probabilistic principle component analyzers

• Mixtures of probabilistic principal component analyzers combine PPCA with Gaussian mixtures to reduce the parameter burden for high-dimensional clustering.The E-step uses the PPCA covariance rather than the data’s full covariance; component covariances are diagonalized to update their local models.Source: §1.5 · p. 11

1.6 A simulated annealing algorithm for PPCA clustering

• The clustering algorithm combines mixtures of PPCA analyzers with annealing, transitioning from spherical-covariance EM to PPCA fitting as σ decreases.The number of principal components at each EM step is the number of covariance eigenvalues above the current σ^2.Source: §1.6 · p. 12

2 Some useful properties of multivariate Gaussians

• Integrating out angular dependence gives the distribution’s dependence on distance from its mean in the covariance metric.The passages identify this distance as r.Source: §2 · p. 12
• The radial distribution is approximately Gaussian, centered at √(d − 1).This is stated for the distribution as a function of r.Source: §2 · p. 12
• Matusita’s measure for two multivariate Gaussians can be computed explicitly.The section presents the expression for this measure.Source: §2 · p. 12

3.1 Frequency of storage

• Collect coordinates at intervals long enough to reduce temporal correlation, using autocorrelation functions to estimate suitable timescales.Storing values at every timestep is not sensible because neighboring samples are highly correlated.Source: §3.1 · p. 13

3.2 Frequency of PPCA

• PPCA frequency trades fitting accuracy, which improves with more points, against the time required to collect statistics.The analysis time is typically insignificant compared with the cost of the dynamics.Source: §3.2 · p. 13

3.3 Weight tolerance

• Set the weight tolerance to determine what fraction of the trajectory a basin must contain to count as a significant kinetic trap.The tolerance applies to clusters not overlapping existing basins; overlap can be neglected when choosing it.Source: §3.3 · p. 13

3.4 Hill height and frequency of addition

• Frequent, low hills are preferred to avoid introducing too many discontinuities into the free energy over time.The authors often attempt hills lower than kBT rather than infrequently adding very large ones.Source: §3.4 · p. 13

3.5 Hill width

• Use the inverse covariance metric so locally adaptive hills are wide along high-fluctuation directions and narrow along low-fluctuation directions.The covariance also provides the average fluctuations in r(s), from which a sensible Gaussian width can be selected; hills should be a small multiple of its standard deviation.Source: §3.5 · pp. 13–14

3.6 Expansion parameter

• The basin-expansion probability compares the hill-addition interval with a diffusion-based time for basin size to grow, limiting growth to a diffusion-like speed.D must be set by trial and error: too small prevents expansion and leaves spurious basins, while too large lets basins keep expanding.Source: §3.6 · pp. 14–15

4 Supplementary information: Lennard-Jones 7

• Four well-tempered metadynamics runs were averaged, using Gaussian widths of 0.02 for both CVs and heights of 0.001 ϵ in two runs and 0.0025 ϵ in two.The notional temperature T + ∆T was 1.0 kBT; harmonic walls with force constants 1000 ϵCV−2 were placed at second-moment values 0.4 and 1.2 and third-moment value -0.3.Source: §4.1 · p. 15

5.1 Simulating ala12 with unbiassed MD

• Even during a 1 µs unbiased trajectory—20 times the reconnaissance-metadynamics timescale—the system stays in a small phase-space region and does not form an alpha helix.The figure includes 40 ps running averages of torsional angles, representative structures, and energies from annealing every 400 ps.Source: §5.1 · pp. 15–16

5.2 Quenched structures for ala12

• A supplied PDB file contains quenched structures and energies from a reconnaissance-metadynamics trajectory for benchmarking methods that explore potential-energy landscapes.Source: §5.2 · p. 17