Replica state exchange metadynamics for improving the convergence of free energy estimates

Metadynamics (MTD) is a powerful enhanced sampling method for systems with rugged energy landscapes. It constructs a bias potential in a predefined collective variable (CV) space to overcome barriers between metastable states. In bias‐exchange MTD (BE‐MTD), multiple replicas approximate the CV space by exchanging bias potentials (replica conditions) with the Metropolis–Hastings (MH) algorithm. We demonstrate that the replica‐exchange rates and the convergence of free energy estimates of BE‐MTD are improved by introducing the infinite swapping (IS) or the Suwa‐Todo (ST) algorithms. Conceptually, IS and ST perform transitions in a replica state space rather than exchanges in a replica condition space. To emphasize this, the proposed scheme is called the replica state exchange MTD (RSE‐MTD). Benchmarks were performed with alanine polypeptides in vacuum and water. For the systems tested in this work, there is no significant performance difference between IS and ST. © 2015 Wiley Periodicals, Inc.     

Acceleration of Markov chain Monte Carlo simulations through sequential updating

Strict detailed balance is not necessary for Markov chain Monte Carlo simulations to converge to the correct equilibrium distribution. In this work, we propose a new algorithm which only satisfies the weaker balance condition, and it is shown analytically to have better mobility over the phase space than the Metropolis algorithm satisfying strict detailed balance. The new algorithm employs sequential updating and yields better sampling statistics than the Metropolis algorithm with random updating. We illustrate the efficiency of the new algorithm on the two-dimensional Ising model. The algorithm is shown to identify the correct equilibrium distribution and to converge faster than the Metropolis algorithm with strict detailed balance. The main advantages of the new algorithm are its simplicity and the feasibility of parallel implementation through domain decomposition.     

Hamiltonian replica‐permutation method and its applications to an alanine dipeptide and amyloid‐β(29–42) peptides

We propose the Hamiltonian replica‐permutation method (RPM) (or multidimensional RPM) for molecular dynamics and Monte Carlo simulations, in which parameters in the Hamiltonian are permuted among more than two replicas with the Suwa‐Todo algorithm. We apply the Coulomb RPM, which is one of realization of the Hamiltonian RPM, to an alanine dipeptide and to two amyloid‐β(29–42) molecules. The Hamiltonian RPM realizes more efficient sampling than the Hamiltonian replica‐exchange method. We illustrate the protein misfolding funnel of amyloid‐β(29–42) and reveal its dimerization pathways. © 2013 Wiley Periodicals, Inc.     

Multiparticle moves in acceptance rate optimized monte carlo

Molecular Dynamics (MD) and Monte Carlo (MC) based simulation methods are widely used to investigate molecular and nanoscale structures and processes. While the investigation of systems in MD simulations is limited by very small time steps, MC methods are often stifled by low acceptance rates for moves that significantly perturb the system. In many Metropolis MC methods with hard potentials, the acceptance rate drops exponentially with the number of uncorrelated, simultaneously proposed moves. In this work, we discuss a multiparticle Acceptance Rate Optimized Monte Carlo approach (AROMoCa) to construct collective moves with near unit acceptance probability, while preserving detailed balance even for large step sizes. After an illustration of the protocol, we demonstrate that AROMoCa significantly accelerates MC simulations in four model systems in comparison to standard MC methods. AROMoCa can be applied to all MC simulations where a gradient of the potential is available and can help to significantly speed up molecular simulations. © 2015 Wiley Periodicals, Inc.     

Transition path sampling for discrete master equations with absorbing states

Transition path sampling (TPS) algorithms have been implemented with deterministic dynamics, with thermostatted dynamics, with Brownian dynamics, and with simple spin flip dynamics. Missing from the TPS repertoire is an implementation with kinetic Monte Carlo (kMC), i.e., with the underlying dynamics coming from a discrete master equation. We present a new hybrid kMC-TPS algorithm and prove that it satisfies detailed balance in the transition path ensemble. The new algorithm is illustrated for a simplified Markov State Model of trp-cage folding. The transition path ensemble from kMC-TPS is consistent with that obtained from brute force kMC simulations. The committor probabilities and local fluxes for the simple model are consistent with those obtained from exact methods for simple master equations. The new kMC-TPS method should be useful for analysis of rare transitions in complex master equations where the individual states cannot be enumerated and therefore where exact solutions cannot be obtained.     

Stochastic potential switching algorithm for Monte Carlo simulations of complex systems

This paper describes a new Monte Carlo method based on a novel stochastic potential switching algorithm. This algorithm enables the equilibrium properties of a system with potential V to be computed using a Monte Carlo simulation for a system with a possibly less complex stochastically altered potential V. By proper choices of the stochastic switching and transition probabilities, it is shown that detailed balance can be strictly maintained with respect to the original potential V. The validity of the method is illustrated with a simple one-dimensional example. The method is then generalized to multidimensional systems with any additive potential, providing a framework for the design of more efficient algorithms to simulate complex systems. A near-critical Lennard-Jones fluid with more than 20,000 particles is used to illustrate the method. The new algorithm produced a much smaller dynamic scaling exponent compared to the Metropolis method and improved sampling efficiency by over an order of magnitude.     

Convergence of folding free energy landscapes via application of enhanced sampling methods in a distributed computing environment

We have implemented the serial replica exchange method (SREM) and simulated tempering (ST) enhanced sampling algorithms in a global distributed computing environment. Here we examine the helix-coil transition of a 21 residue alpha-helical peptide in explicit solvent. For ST, we demonstrate the efficacy of a new method for determining initial weights allowing the system to perform a random walk in temperature space based on short trial simulations. These weights are updated throughout the production simulation by an adaptive weighting method. We give a detailed comparison of SREM, ST, as well as standard MD and find that SREM and ST give equivalent results in reasonable agreement with experimental data. In addition, we find that both enhanced sampling methods are much more efficient than standard MD simulations. The melting temperature of the Fs peptide with the AMBER99phi potential was calculated to be about 310 K, which is in reasonable agreement with the experimental value of 334 K. We also discuss other temperature dependent properties of the helix-coil transition. Although ST has certain advantages over SREM, both SREM and ST are shown to be powerful methods via distributed computing and will be applied extensively in future studies of complex bimolecular systems.     

Interior segment regrowth configurational-bias algorithm for the efficient sampling and fast relaxation of coarse-grained polyethylene and polyoxyethylene melts on a high coordination lattice

We demonstrate the application of a modified form of the configurational-bias algorithm for the simulation of chain molecules on the second-nearest-neighbor-diamond lattice. Using polyethylene and poly(ethylene-oxide) as model systems we show that the present configurational-bias algorithm can increase the speed of the equilibration by at least a factor of 2-3 or more as compared to the previous method of using a combination of single-bead and pivot moves along with the Metropolis sampling scheme [N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, J. Chem. Phys. 21, 1087 (1953)]. The increase in the speed of the equilibration is found to be dependent on the interactions (i.e., the polymer being simulated) and the molecular weight of the chains. In addition, other factors not considered, such as the density, would also have a significant effect. The algorithm is an extension of the conventional configurational-bias method adapted to the regrowth of interior segments of chain molecules. Appropriate biasing probabilities for the trial moves as outlined by Jain and de Pablo for the configurational-bias scheme of chain ends, suitably modified for the interior segments, are utilized [T. S. Jain and J. J. de Pablo, in Simulation Methods for Polymers, edited by M. Kotelyanskii and D. N. Theodorou (Marcel Dekker, New York, 2004), pp. 223-255]. The biasing scheme satisfies the condition of detailed balance and produces efficient sampling with the correct equilibrium probability distribution of states. The method of interior regrowth overcomes the limitations of the original configurational-bias scheme and allows for the simulation of polymers of higher molecular weight linear chains and ring polymers which lack chain ends.     

The temperature intervals with global exchange of replicas empirical accelerated sampling method: parameter sensitivity and extension to a complex molecular system

The recently developed "temperature intervals with global exchange of replicas" (TIGER2) algorithm is an efficient replica-exchange sampling algorithm that provides the freedom to specify the number of replicas and temperature levels independently of the size of the system and temperature range to be spanned, thus making it particularly well suited for sampling molecular systems that are considered to be too large to be sampled using conventional replica exchange methods. Although the TIGER2 method is empirical in nature, when appropriately applied it is able to provide sampling that satisfies the balance condition and closely approximates a Boltzmann-weighted ensemble of states. In this work, we evaluated the influence of factors such as temperature range, temperature spacing, replica number, and sampling cycle design on the accuracy of a TIGER2 simulation based on molecular dynamics simulations of alanine dipeptide in implicit solvent. The influence of these factors is further examined by calculating the properties of a complex system composed of the B1 immunoglobulin-binding domain of streptococcal protein G (protein G) in aqueous solution. The accuracy of a TIGER2 simulation is particularly sensitive to the maximum temperature level selected for the simulation. A method to determine the appropriate maximum temperature level to be used in a TIGER2 simulation is presented.     

Multiple Markov transition matrix method: Obtaining the stationary probability distribution from multiple simulations

We herein propose the multiple Markov transition matrix method (MMMM), an algorithm by which to estimate the stationary probability distribution from independent multiple molecular dynamics simulations with different Hamiltonians. Applications to the potential of mean force calculation in combination with the umbrella sampling method are presented. First, the performance of the MMMM is examined in the case of butane. Compared with the weighted histogram analysis method (WHAM), the MMMM has an advantage with respect to the reasonable evaluation of the stationary probability distribution even from nonequilibrium trajectories. This method is then applied to Met‐enkephalin nonequilibrium simulation. © 2008 Wiley Periodicals, Inc. J Comput Chem, 2009     

Cool walking: a new Markov chain Monte Carlo sampling method

Effective relaxation processes for difficult systems like proteins or spin glasses require special simulation techniques that permit barrier crossing to ensure ergodic sampling. Numerous adaptations of the venerable Metropolis Monte Carlo (MMC) algorithm have been proposed to improve its sampling efficiency, including various hybrid Monte Carlo (HMC) schemes, and methods designed specifically for overcoming quasi-ergodicity problems such as Jump Walking (J-Walking), Smart Walking (S-Walking), Smart Darting, and Parallel Tempering. We present an alternative to these approaches that we call Cool Walking, or C-Walking. In C-Walking two Markov chains are propagated in tandem, one at a high (ergodic) temperature and the other at a low temperature. Nonlocal trial moves for the low temperature walker are generated by first sampling from the high-temperature distribution, then performing a statistical quenching process on the sampled configuration to generate a C-Walking jump move. C-Walking needs only one high-temperature walker, satisfies detailed balance, and offers the important practical advantage that the high and low-temperature walkers can be run in tandem with minimal degradation of sampling due to the presence of correlations. To make the C-Walking approach more suitable to real problems we decrease the required number of cooling steps by attempting to jump at intermediate temperatures during cooling. We further reduce the number of cooling steps by utilizing "windows" of states when jumping, which improves acceptance ratios and lowers the average number of cooling steps. We present C-Walking results with comparisons to J-Walking, S-Walking, Smart Darting, and Parallel Tempering on a one-dimensional rugged potential energy surface in which the exact normalized probability distribution is known. C-Walking shows superior sampling as judged by two ergodic measures.     

Designing generalized statistical ensembles for numerical simulations of biopolymers

Conformational properties of polymers, such as average dihedral angles or molecular alpha-helicity, display a rather weak dependence on the detailed arrangement of the elementary constituents (atoms). We propose a computer simulation method to explore the polymer phase space using a variant of the standard multicanonical method, in which the density of states associated to suitably chosen configurational variables is considered in place of the standard energy density of states. This configurational density of states is used in the Metropolis acceptance/rejection test when configurations are generated with the help of a hybrid Monte Carlo algorithm. The resulting configurational probability distribution is then modulated by exponential factors derived from the general principle of the maximal constrained entropy by requiring that certain average configurational quantities take preassigned (possibly temperature dependent) values. Thermal averages of other configurational quantities can be computed by using the probability distributions obtained in this way. Moments of the energy distribution require an extra canonical sampling of the system phase space at the desired temperature, in order to locally thermalize the configurational degrees of freedom. As an application of these ideas we present the study of the structural properties of two simple models: a bead-and-spring model of polyethylene with independent hindered torsions and an all-atom model of alanine and glycine oligomers with 12 amino acids in vacuum.     

Probability distributions of molecular observables computed from Markov models

Molecular dynamics (MD) simulations can be used to estimate transition rates between conformational substates of the simulated molecule. Such an estimation is associated with statistical uncertainty, which depends on the number of observed transitions. In turn, it induces uncertainties in any property computed from the simulation, such as free energy differences or the time scales involved in the system's kinetics. Assessing these uncertainties is essential for testing the reliability of a given observation and also to plan further simulations in such a way that the most serious uncertainties will be reduced with minimal effort. Here, a rigorous statistical method is proposed to approximate the complete statistical distribution of any observable of an MD simulation provided that one can identify conformational substates such that the transition process between them may be modeled with a memoryless jump process, i.e., Markov or Master equation dynamics. The method is based on sampling the statistical distribution of Markov transition matrices that is induced by the observed transition events. It allows physically meaningful constraints to be included, such as sampling only matrices that fulfill detailed balance, or matrices that produce a predefined equilibrium distribution of states. The method is illustrated on mus MD simulations of a hexapeptide for which the distributions and uncertainties of the free energy differences between conformations, the transition matrix elements, and the transition matrix eigenvalues are estimated. It is found that both constraints, detailed balance and predefined equilibrium distribution, can significantly reduce the uncertainty of some observables.     

Ingenuity in performing replica permutation: How to order the state labels for improving sampling efficiency

In the replica-permutation method, an advanced version of the replica-exchange method, all combinations of replicas and parameters are considered for parameter permutation, and a list of all the combinations is prepared. Here, we report that the temperature transition probability depends on how the list is created, especially in replica permutation with solute tempering (RPST). We found that the transition probabilities decrease at large replica indices when the combinations are sequentially assigned to the state labels as in the originally proposed list. To solve this problem, we propose to modify the list by randomly assigning the combinations to the state labels. We performed molecular dynamics simulations of amyloid-β(16–22) peptides using RPST with the “randomly assigned” list (RPST-RA) and RPST with the “sequentially assigned” list (RPST-SA). The results show the decreases in the transition probabilities in RPST-SA are eliminated, and the sampling efficiency is improved in RPST-RA.     

Steered transition path sampling

We introduce a path sampling method for obtaining statistical properties of an arbitrary stochastic dynamics. The method works by decomposing a trajectory in time, estimating the probability of satisfying a progress constraint, modifying the dynamics based on that probability, and then reweighting to calculate averages. Because the progress constraint can be formulated in terms of occurrences of events within time intervals, the method is particularly well suited for controlling the sampling of currents of dynamic events. We demonstrate the method for calculating transition probabilities in barrier crossing problems and survival probabilities in strongly diffusive systems with absorbing states, which are difficult to treat by shooting. We discuss the relation of the algorithm to other methods.     

Detailed balance condition and effective free energy in the primitive chain network model

We consider statistical mechanical properties of the primitive chain network (PCN) model for entangled polymers from its dynamic equations. We show that the dynamic equation for the segment number of the PCN model does not reduce to the standard Langevin equation which satisfies the detailed balance condition. We propose heuristic modifications for the PCN dynamic equation for the segment number, to make it reduce to the standard Langevin equation. We analyse some equilibrium statistical properties of the modified PCN model, by using the effective free energy obtained from the modified PCN dynamic equations. The PCN effective free energy can be interpreted as the sum of the ideal Gaussian chain free energy and the repulsive interaction energy between slip-links. By using the single chain approximation, we calculate several distribution functions of the PCN model. The obtained distribution functions are qualitatively different from ones for the simple slip-link model without any direct interactions between slip-links.     

Unbiased sampling of lattice Hamilton path ensembles

Hamilton paths, or Hamiltonian paths, are walks on a lattice which visit each site exactly once. They have been proposed as models of globular proteins and of compact polymers. A previously published algorithm [Mansfield, Macromolecules 27, 5924 (1994)] for sampling Hamilton paths on simple square and simple cubic lattices is tested for bias and for efficiency. Because the algorithm is a Metropolis Monte Carlo technique obviously satisfying detailed balance, we need only demonstrate ergodicity to ensure unbiased sampling. Two different tests for ergodicity (exact enumeration on small lattices, nonexhaustive enumeration on larger lattices) demonstrate ergodicity unequivocally for small lattices and provide strong support for ergodicity on larger lattices. Two other sampling algorithms [Ramakrishnan et al., J. Chem. Phys. 103, 7592 (1995); Lua et al., Polymer 45, 717 (2004)] are both known to produce biases on both 2x2x2 and 3x3x3 lattices, but it is shown here that the current algorithm gives unbiased sampling on these same lattices. Successive Hamilton paths are strongly correlated, so that many iterations are required between statistically independent samples. Rules for estimating the number of iterations needed to dissipate these correlations are given. However, the iteration time is so fast that the efficiency is still very good except on extremely large lattices. For example, even on lattices of total size 10x10x10 we are able to generate tens of thousands of uncorrelated Hamilton paths per hour of CPU time.     

On the acceptance ratio in Monte Carlo computer simulations

The functional form of acceptance probabilities in Monte Carlo algorithms bears a resemblance to the distance functions which are specifically defined to be bracketed by the unit interval. This observation led us to seek the average distance between any two points on the unit interval and this by analogy resulted in a suggestion of an upper and a lower bound of 1/2 and 1/3, respectively, for the acceptance ratio or the average acceptance probability in Monte Carlo computer simulations.     

Improved transition path sampling methods for simulation of rare events

The free energy surfaces of a wide variety of systems encountered in physics, chemistry, and biology are characterized by the existence of deep minima separated by numerous barriers. One of the central aims of recent research in computational chemistry and physics has been to determine how transitions occur between deep local minima on rugged free energy landscapes, and transition path sampling (TPS) Monte-Carlo methods have emerged as an effective means for numerical investigation of such transitions. Many of the shortcomings of TPS-like approaches generally stem from their high computational demands. Two new algorithms are presented in this work that improve the efficiency of TPS simulations. The first algorithm uses biased shooting moves to render the sampling of reactive trajectories more efficient. The second algorithm is shown to substantially improve the accuracy of the transition state ensemble by introducing a subset of local transition path simulations in the transition state. The system considered in this work consists of a two-dimensional rough energy surface that is representative of numerous systems encountered in applications. When taken together, these algorithms provide gains in efficiency of over two orders of magnitude when compared to traditional TPS simulations.     

Replica sub-permutation method for molecular dynamics and monte carlo simulations

We propose an improvement of the replica-exchange and replica-permutation methods, which we call the replica sub-permutation method (RSPM). Instead of considering all permutations, this method uses a new algorithm referred to as sub-permutation to perform parameter transition. The RSPM succeeds in reducing the number of combinations between replicas and parameters without the loss of sampling efficiency. For comparison, we applied the replica sub-permutation, replica-permutation, and replica-exchange methods to a β-hairpin mini protein, chignolin, in explicit water. We calculated the transition ratio and number of tunneling events in the parameter space, the number of folding–unfolding events, the autocorrelation function, and the autocorrelation time as measures of sampling efficiency. The results indicate that among the three methods, the proposed RSPM is the most efficient in both parameter and conformational spaces. © 2019 Wiley Periodicals, Inc.