Sampling enhancement for the quantum mechanical potential based molecular dynamics simulations: a general algorithm and its extension for free energy calculation on rugged energy surface

An approach is developed in the replica exchange framework to enhance conformational sampling for the quantum mechanical (QM) potential based molecular dynamics simulations. Importantly, with our enhanced sampling treatment, a decent convergence for electronic structure self-consistent-field calculation is robustly guaranteed, which is made possible in our replica exchange design by avoiding direct structure exchanges between the QM-related replicas and the activated (scaled by low scaling parameters or treated with high "effective temperatures") molecular mechanical (MM) replicas. Although the present approach represents one of the early efforts in the enhanced sampling developments specifically for quantum mechanical potentials, the QM-based simulations treated with the present technique can possess the similar sampling efficiency to the MM based simulations treated with the Hamiltonian replica exchange method (HREM). In the present paper, by combining this sampling method with one of our recent developments (the dual-topology alchemical HREM approach), we also introduce a method for the sampling enhanced QM-based free energy calculations.     

Simulated scaling method for localized enhanced sampling and simultaneous "alchemical" free energy simulations: a general method for molecular mechanical, quantum mechanical, and quantum mechanical/molecular mechanical simulations

A potential scaling version of simulated tempering is presented to efficiently sample configuration space in a localized region. The present "simulated scaling" method is developed with a Wang-Landau type of updating scheme in order to quickly flatten the distributions in the scaling parameter lambdam space. This proposal is meaningful for a broad range of biophysical problems, in which localized sampling is required. Besides its superior capability and robustness in localized conformational sampling, this simulated scaling method can also naturally lead to efficient "alchemical" free energy predictions when dual-topology alchemical hybrid potential is applied; thereby simultaneously, both of the chemically and conformationally distinct portions of two end point chemical states can be efficiently sampled. As demonstrated in this work, the present method is also feasible for the quantum mechanical and quantum mechanical/molecular mechanical simulations.     

Computing Alchemical Free Energy Differences with Hamiltonian Replica Exchange Molecular Dynamics (H-REMD) Simulations

Alchemical free energy calculations play a very important role in the field of molecular modeling. Efforts have been made to improve the accuracy and precision of those calculations. One of the efforts is to employ a Hamiltonian replica exchange molecular dynamics (H-REMD) method to enhance conformational sampling. In this paper, we demonstrated that HREMD method not only improves convergence in alchemical free energy calculations but also can be used to compute free energy differences directly via the Free Energy Perturbation (FEP)algorithm. We show a direct mapping between the H-REMD and the usual FEP equations, which are then used directly to compute free energies. The H-REMD alchemical free energy calculation (Replica exchange Free Energy Perturbation, REFEP) was tested on predicting the pK(a) value of the buried Asp26 in thioredoxin. We compare the results of REFEP with TI and regular FEP simulations. REFEP calculations converged faster than those from TI and regular FEP simulations. The final predicted pK(a) value from the H-REMD simulation was also very accurate, only 0.4 pK(a) unit above the experimental value. Utilizing the REFEP algorithm significantly improves conformational sampling, and this in turn improves the convergence of alchemical free energy simulations.     

Free‐energy differences between states with different conformational ensembles

Multiple conformations separated by high‐energy barriers represent a challenging problem in free‐energy calculations due to the difficulties in achieving adequate sampling. We present an application of thermodynamic integration (TI) in conjunction with the local elevation umbrella sampling (LE/US) method to improve convergence in alchemical free‐energy calculations. TI‐LE/US was applied to the guanosine triphosphate (GTP) to 8‐Br‐GTP perturbation, molecules that present high‐energy barriers between the anti and syn states and that have inverted preferences for those states. The convergence and reliability of TI‐LE/US was assessed by comparing with previous results using the enhanced‐sampling one‐step perturbation (OSP) method. A linear interpolation of the end‐state biasing potentials was sufficient to dramatically improve sampling along the chosen reaction coordinate. Conformational free‐energy differences were also computed for the syn and anti states and compared to experimental and theoretical results. Additionally, a coupled OSP with LE/US was carried out, allowing the calculation of conformational and alchemical free energies of GTP and 8‐substituted GTP analogs. © 2013 Wiley Periodicals, Inc.     

Alchemical free energy calculations and multiple conformational substates

Thermodynamic integration (TI) was combined with (adaptive) umbrella sampling to improve the convergence of alchemical free energy simulations in which multiple conformational substates are present. The approach, which we refer to as non-Boltzmann TI (NBTI), was tested by computing the free energy differences between three five-atomic model systems, as well as the free energy difference of solvation between leucine and asparagine. In both cases regular TI failed to give converged results, whereas the NBTI results were free from hysteresis and had standard deviations well below +/-0.7 kcal/mole. We also present theoretical considerations that make it possible to compute free energy differences between simple molecules, such as the five-atomic model systems, by numerical integration of the partition functions at the respective end points.     

A hybrid recursion method to robustly ensure convergence efficiencies in the simulated scaling based free energy simulations

Recently, we developed an efficient free energy simulation technique, the simulated scaling (SS) method [H. Li et al., J. Chem. Phys. 126, 024106 (2007)], in the framework of generalized ensemble simulations. In the SS simulations, random walks in the scaling parameter space are realized so that both phase space overlap sampling and conformational space sampling can be simultaneously enhanced. To flatten the distribution in the scaling parameter space, in the original SS implementation, the Wang-Landau recursion was employed due to its well-known recursion capability. In the Wang-Landau recursion based SS free energy simulation scheme, at the early stage, recursion efficiencies are high and free energy regions are quickly located, although at this stage, the errors of estimated free energy values are large; at the later stage, the errors of estimated free energy values become smaller, however, recursions become increasingly slow and free energy refinements require very long simulation time. In order to robustly resolve this efficiency problem during free energy refinements, a hybrid recursion strategy is presented in this paper. Specifically, we let the Wang-Landau update method take care of the early stage recursion: the location of target free energy regions, and let the adaptive reweighting method take care of the late stage recursion: the refinements of free energy values. As comparably studied in the model systems, among three possible recursion procedures, the adaptive reweighting recursion approach is the least favorable one because of its low recursion efficiency during free energy region locations; and compared to the original Wang-Landau recursion approach, the proposed hybrid recursion technique can be more robust to guarantee free energy simulation efficiencies.     

Accelerated convergence of molecular free energy via superposition approximation-based reference states

The free energy of a molecular system can, at least in principle, be computed by thermodynamic perturbation from a reference system whose free energy is known. The convergence of such a calculation depends critically on the conformational overlap between the reference and the physical systems. One approach to defining a suitable reference system is to construct it from the one-dimensional marginal probability distribution functions (PDFs) of internal coordinates observed in a molecular simulation. However, the conformational overlap of this reference system tends to decline steeply with increasing dimensionality, due to the neglect of correlations among the coordinates. Here, we test a reference system that can account for pairwise correlations among the internal coordinates, as captured by their two-dimensional marginal PDFs derived from a molecular simulation. Incorporating pairwise correlations in the reference system is found to dramatically improve the convergence of the free energy estimates relative to the first-order reference system, due to increased conformational overlap with the physical distribution.     

Free Energy Calculations using a Swarm‐Enhanced Sampling Molecular Dynamics Approach

Free energy simulations are an established computational tool in modelling chemical change in the condensed phase. However, sampling of kinetically distinct substates remains a challenge to these approaches. As a route to addressing this, we link the methods of thermodynamic integration (TI) and swarm‐enhanced sampling molecular dynamics (sesMD), where simulation replicas interact cooperatively to aid transitions over energy barriers. We illustrate the approach by using alchemical alkane transformations in solution, comparing them with the multiple independent trajectory TI (IT‐TI) method. Free energy changes for transitions computed by using IT‐TI grew increasingly inaccurate as the intramolecular barrier was heightened. By contrast, swarm‐enhanced sampling TI (sesTI) calculations showed clear improvements in sampling efficiency, leading to more accurate computed free energy differences, even in the case of the highest barrier height. The sesTI approach, therefore, has potential in addressing chemical change in systems where conformations exist in slow exchange.     

BAR-based optimum adaptive steered MD for configurational sampling

The equilibrium and nonequilibrium adaptive alchemical free energy simulation methods optimum Bennett's acceptance ratio and optimum crooks' equation (OCE), based on the statistically optimal bidirectional reweighting estimator named Bennett's Acceptance Ratio or Crooks' equation, perform initial sampling in the staging alchemical transformation and then determine the importance rank of different states via the time-derivative of the variance. The method is proven to give speedups compared with the equal time rule. In the current work, we extend the time derivative of variance guided adaptive sampling method to the configurational space, falling in the term of steered MD (SMD). The SMD approach biasing physically meaningful collective variable (CV) such as one dihedral or one distance to pulling the system from one conformational state to another. By minimizing the variance of the free energy differences along the pathway in an optimized way, a new type of adaptive SMD (ASMD) is introduced. As exhibits in the alchemical case, this adaptive sampling method outperforms the traditional equal-time SMD in nonequilibrium stratification. Also, the method gives much more efficient calculation of potential of mean force than the selection criterion-based ASMD scheme, which is proven to be more efficient than traditional SMD. The OCE workflow is periodicity-of-CV dependent while ASMD is not. The performance is demonstrated in a dihedral flipping case and two distance pulling cases, accounting for periodic and nonperiodic CVs, respectively. © 2019 Wiley Periodicals, Inc.     

Free energy simulations: use of reverse cumulative averaging to determine the equilibrated region and the time required for convergence

A method is proposed for improving the accuracy and efficiency of free energy simulations. The essential idea is that the convergence of the relevant measure (e.g., the free energy derivative in thermodynamic integration) is monitored in the reverse direction starting from the last frame of the trajectory, instead of the usual approach, which begins with the first frame and goes in the forward direction. This simple change in the use of the simulation data makes it straightforward to eliminate the contamination of the averages by contributions from the equilibrating region. A statistical criterion is introduced for distinguishing the equilibrated (production) region from the equilibrating region. The proposed method, called reverse cumulative averaging, is illustrated by its application to the well-studied case of the alchemical free energy simulation of ethane to methanol.     

Improved convergence in dual-topology free energy calculations through use of harmonic restraints

When calculating free energy differences between two molecular systems by means of molecular dynamics simulation, accessory potential functions can help eliminate uninteresting configurational entropy contributions, improve convergence, and facilitate reversibility. In this work, we demonstrate that the use of a harmonic potential function to restrain key portions of a molecular system in a free energy perturbation dual-topology molecular dynamics approach dramatically improves convergence and precision of the calculation. Limitations of this technique are illustrated, and its use in conjunction with a fixed bond-length constraint is developed. © 1998 John Wiley & Sons, Inc. J Comput Chem 19: 1278–1283, 1998     

Combining the biased and unbiased sampling strategy into one convenient free energy calculation method

Constructing a free energy landscape for a large molecule is difficult. One has to use either a high temperature or a strong driving force to enhance the sampling on the free energy barriers. In this work, we propose a mixed method that combines these two kinds of acceleration strategies into one simulation. First, it applies an adaptive biasing potential to some replicas of the molecule. These replicas are particularly accelerated in a collective variable space. Second, it places some unbiased and exchangeable replicas at various temperature levels. These replicas generate unbiased sampling data in the canonical ensemble. To improve the sampling efficiency, biased replicas transfer their state variables to the unbiased replicas after equilibrium by Monte Carlo trial moves. In comparison to previous integrated methods, it is more convenient for users. It does not need an initial reference biasing potential to guide the sampling of the molecule. And it is also unnecessary to insert many replicas for the requirement of passing the free energy barriers. The free energy calculation is accomplished in a single stage. It samples the data as fast as a biased simulation and it processes the data as simple as an unbiased simulation. The method provides a minimalist approach to the construction of the free energy landscape. © 2019 Wiley Periodicals, Inc.     

Improving the efficiency and reliability of free energy perturbation calculations using overlap sampling methods

A challenge in free energy calculation for complex molecular systems by computer simulation is to obtain a reliable estimate within feasible computational time. In this study, we suggest an answer to this challenge by exploring a simple method, overlap sampling (OS), for producing reliable free-energy results in an efficient way. The formalism of the OS method is based on ensuring sampling of important overlapping phase space during perturbation calculations. This technique samples both forward and reverse free energy perturbation (FEP) to improve the free-energy calculation. It considers the asymmetry of the FEP calculation and features an ability to optimize both the precision and the accuracy of the measurement without affecting the simulation process itself. The OS method is tested at two optimization levels: no optimization (simple OS), and full optimization (equivalent to Bennett's method), and compared to conventional FEP techniques, including the widely used direct FEP averaging method, on three alchemical mutation systems: (a) an anion transformation in water solution, (b) mutation between methanol and ethane, and (c) alchemical change of an adenosine molecule. It is consistently shown that the reliability of free-energy estimates can be greatly improved using the OS techniques at both optimization levels, while the performance of Bennett's method is particularly striking. In addition, the efficiency of a calculation can be significantly improved because the method is able to (a) converge to the right answer quickly, and (b) work for large perturbations. The basic two-stage OS method can be extended to admit additional stages, if needed. We suggest that the OS method can be used as a general perturbation technique for computing free energy differences in molecular simulations.     

A convective replica‐exchange method for sampling new energy basins

Replica‐exchange is a powerful simulation method for sampling the basins of a rugged energy landscape. The replica‐exchange method's sampling is efficient because it allows replicas to perform round trips in temperature space, thereby visiting both low and high temperatures in the same simulation. However, replicas have a diffusive walk in temperature space, and the round trip rate decreases significantly with the system size. These drawbacks make convergence of the simulation even more difficult than it already is when bigger systems are tackled. Here, we present a simple modification of the exchange method. In this method, one of the replicas steadily raises or lowers its temperature. We tested the convective replica‐exchange method on three systems of varying complexity: the alanine dipeptide in implicit solvent, the GB1 β‐hairpin in explicit solvent and the Aβ_25–35 homotrimer in a coarse grained representation. For the highly frustrated Aβ_25–35 homotrimer, the proposed “convective” replica‐exchange method is twice as fast as the standard method. It discovered 24 out of 27 free‐energy basins in less than 500 ns. It also prevented the formation of groups of replicas that usually form on either side of an exchange bottleneck, leading to a more efficient sampling of new energy basins than in the standard method. © 2012 Wiley Periodicals, Inc.     

Peptide conformational equilibria computed via a single-stage shifting protocol

We study the conformational equilibria of two peptides using a novel statistical mechanics approach designed for calculating free energy differences between highly dissimilar conformational states. Our results elucidate the contrasting roles of entropy in implicitly solvated leucine dipeptide and decaglycine. The method extends earlier work by Voter and overcomes the notorious "overlap" problem in free energy computations by constructing a mathematically equivalent calculation with high conformational similarity. The approach requires only equilibrium simulations of the two states of interest, without the need for sampling transition states. We discuss possible extensions and optimizations of the approach.     

From A to B in free energy space

The authors present a new method for searching low free energy paths in complex molecular systems at finite temperature. They introduce two variables that are able to describe the position of a point in configurational space relative to a preassigned path. With the help of these two variables the authors combine features of approaches such as metadynamics or umbrella sampling with those of path based methods. This allows global searches in the space of paths to be performed and a new variational principle for the determination of low free energy paths to be established. Contrary to metadynamics or umbrella sampling the path can be described by an arbitrary large number of variables, still the energy profile along the path can be calculated. The authors exemplify the method numerically by studying the conformational changes of alanine dipeptide.     

Tinker‐OpenMM: Absolute and relative alchemical free energies using AMOEBA on GPUs

The capabilities of the polarizable force fields for alchemical free energy calculations have been limited by the high computational cost and complexity of the underlying potential energy functions. In this work, we present a GPU‐based general alchemical free energy simulation platform for polarizable potential AMOEBA. Tinker‐OpenMM, the OpenMM implementation of the AMOEBA simulation engine has been modified to enable both absolute and relative alchemical simulations on GPUs, which leads to a ∼200‐fold improvement in simulation speed over a single CPU core. We show that free energy values calculated using this platform agree with the results of Tinker simulations for the hydration of organic compounds and binding of host–guest systems within the statistical errors. In addition to absolute binding, we designed a relative alchemical approach for computing relative binding affinities of ligands to the same host, where a special path was applied to avoid numerical instability due to polarization between the different ligands that bind to the same site. This scheme is general and does not require ligands to have similar scaffolds. We show that relative hydration and binding free energy calculated using this approach match those computed from the absolute free energy approach. © 2017 Wiley Periodicals, Inc.     

A combined treatment of hydration and dynamical effects for the modeling of host–guest binding thermodynamics: the SAMPL5 blinded challenge

As part of the SAMPL5 blinded experiment, we computed the absolute binding free energies of 22 host–guest complexes employing a novel approach based on the BEDAM single-decoupling alchemical free energy protocol with parallel replica exchange conformational sampling and the AGBNP2 implicit solvation model specifically customized to treat the effect of water displacement as modeled by the Hydration Site Analysis method with explicit solvation. Initial predictions were affected by the lack of treatment of ionic charge screening, which is very significant for these highly charged hosts, and resulted in poor relative ranking of negatively versus positively charged guests. Binding free energies obtained with Debye–Hückel treatment of salt effects were in good agreement with experimental measurements. Water displacement effects contributed favorably and very significantly to the observed binding affinities; without it, the modeling predictions would have grossly underestimated binding. The work validates the implicit/explicit solvation approach employed here and it shows that comprehensive physical models can be effective at predicting binding affinities of molecular complexes requiring accurate treatment of conformational dynamics and hydration.     

Applying efficient implicit nongeometric constraints in alchemical free energy simulations

Several strategies have been developed for satisfying bond lengths, angle, and other geometric constraints in molecular dynamics simulations. Advanced variations of alchemical free energy perturbation simulations, however, also require nongeometric constraints. In our recently developed multisite λ-dynamics simulation method, the conventional λ parameters that are associated with the progress variables in alchemical transformations are treated as dynamic variables and are constrained such that: 0 ≤ λ(i) ≤ 1 and Σ(i = 1)(N) λ(i) = 1. Here, we present four functional forms of λ that implicitly satisfy these nongeometric constraints, whose values and forces are facile to compute and that yield stable simulations using a 2 fs integration timestep. Using model systems, we present the sampling characteristics of these functional forms and demonstrate the enhanced sampling profiles and improved convergence rates that are achieved by the functional form: λ(i) = e(c sinθ(i))/Σ(j = 1)(N) e(c sinθ(j)) that oscillates between λ(i) = 0 and λ(i) = 1 and has relatively steep transitions between these endpoints.     

Coulomb replica‐exchange method: Handling electrostatic attractive and repulsive forces for biomolecules

We propose a new type of the Hamiltonian replica‐exchange method (REM) for molecular dynamics (MD) and Monte Carlo simulations, which we refer to as the Coulomb REM (CREM). In this method, electrostatic charge parameters in the Coulomb interactions are exchanged among replicas while temperatures are exchanged in the usual REM. By varying the atom charges, the CREM overcomes free‐energy barriers and realizes more efficient sampling in the conformational space than the REM. Furthermore, this method requires only a smaller number of replicas because only the atom charges of solute molecules are used as exchanged parameters. We performed Coulomb replica‐exchange MD simulations of an alanine dipeptide in explicit water solvent and compared the results with those of the conventional canonical, replica exchange, and van der Waals REMs. Two force fields of AMBER parm99 and AMBER parm99SB were used. As a result, the CREM sampled all local‐minimum free‐energy states more frequently than the other methods for both force fields. Moreover, the Coulomb, van der Waals, and usual REMs were applied to a fragment of an amyloid‐β peptide (Aβ) in explicit water solvent to compare the sampling efficiency of these methods for a larger system. The CREM sampled structures of the Aβ fragment more efficiently than the other methods. We obtained β‐helix, α‐helix, 3₁₀‐helix, β‐hairpin, and β‐sheet structures as stable structures and deduced pathways of conformational transitions among these structures from a free‐energy landscape. © 2012 Wiley Periodicals, Inc.