Modified replica exchange simulation methods for local structure refinement

Parallel tempering, also known as replica exchange molecular dynamics (REMD), has recently been successfully used to study the structure and thermodynamic properties of biomolecules such as peptides and small proteins. For large systems, however, applying REMD can be costly since the number of replicas needed increases as the square root of the number of degrees of freedom in the system. Often, enhanced sampling is only needed for a subset of atoms, such as a loop region of a large protein or a small ligand binding to a receptor. In such applications, it is often reasonable to assume a weak dependence of the structure of the larger region on the instantaneous conformation of the smaller region of interest. For these cases, we derived two variant replica exchange methods, partial replica exchange molecular dynamics (PREMD) and local replica exchange molecular dynamics (LREMD). The Hamiltonian for the system is separated, with replica exchange carried out only for terms involving the subsystem of interest while the remainder of the system is maintained at a single temperature. The number of replicas required for efficient exchange thus depends on the number of degrees of freedom in the fragment needing refinement rather than on the size of the full system. The method can be applied to much larger systems than was previously practical. This also provides a means to preserve the integrity of the structure outside the refinement region without introduction of restraints. LREMD takes this weak coupling approximation a step further, employing only a single representation of the large fragment that simultaneously interacts with all of the replicas of the subsystem of interest. This is obtained by combining replica exchange with the locally enhanced sampling approximation (LES), reducing the computational expense of replica exchange simulations to near that of a single standard molecular dynamics (MD) simulation. Use of LREMD also permits the use of LES without requiring the specification of a single temperature, a known difficulty for standard LES simulations. We tested these two methods on the loop region of an RNA hairpin model system and find significant advantages over standard MD and REMD simulations.     

Synergistic approach to improve "alchemical" free energy calculation in rugged energy surface

The authors present an integrated approach to "alchemical" free energy simulation, which permits efficient calculation of the free energy difference on rugged energy surface. The method is designed to obtain efficient canonical sampling for rapid free energy convergence. The proposal is motivated by the insight that both the exchange efficiency in the presently designed dual-topology alchemical Hamiltonian replica exchange method (HREM), and the confidence of the free energy determination using the overlap histogramming method, depend on the same criterion, viz., the overlaps of the energy difference histograms between all pairs of neighboring states. Hence, integrating these two techniques can produce a joint solution to the problems of the free energy convergence and conformational sampling in the free energy simulations, in which lambda parameter plays two roles to simultaneously facilitate the conformational sampling and improve the phase space overlap for the free energy determination. Specifically, in contrast with other alchemical HREM based free energy simulation methods, the dual-topology approach can ensure robust conformational sampling. Due to these features (a synergistic solution to the free energy convergence and canonical sampling, and the improvement of the sampling efficiency with the dual-topology treatment), the present approach, as demonstrated in the model studies of the authors, is highly efficient in obtaining accurate free energy differences, especially for the systems with rough energy landscapes.     

Coupling of replica exchange simulations to a non-Boltzmann structure reservoir

Computing converged ensemble properties remains challenging for large biomolecules. Replica exchange molecular dynamics (REMD) can significantly increase the efficiency of conformational sampling by using high temperatures to escape kinetic traps. Several groups, including ours, introduced the idea of coupling replica exchange to a pre-converged, Boltzmann-populated reservoir, usually at a temperature higher than that of the highest temperature replica. This procedure reduces computational cost because the long simulation times needed for extensive sampling are only carried out for a single temperature. However, a weakness of the approach is that the Boltzmann-weighted reservoir can still be difficult to generate. We now present the idea of employing a non-Boltzmann reservoir, whose structures can be generated through more efficient conformational sampling methods. We demonstrate that the approach is rigorous and derive a correct statistical mechanical exchange criterion between the reservoir and the replicas that drives Boltzmann-weighted probabilities for the replicas. We test this approach on the trpzip2 peptide and demonstrate that the resulting thermal stability profile is essentially indistinguishable from that obtained using very long (>100 ns) standard REMD simulations. The convergence of this reservoir-aided REMD is significantly faster than for regular REMD. Furthermore, we demonstrate that modification of the exchange criterion is essential; REMD simulations using a standard exchange function with the non-Boltzmann reservoir produced incorrect results.     

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.     

Walking freely in the energy and temperature space by the modified replica exchange molecular dynamics method

Replica Exchange Molecular Dynamics (REMD) method is a powerful sampling tool in molecular simulations. Recently, we made a modification to the standard REMD method. It places some inactive replicas at different temperatures as well as the active replicas. The method completely decouples the number of the active replicas and the number of the temperature levels. In this article, we make a further modification to our previous method. It uses the inactive replicas in a different way. The inactive replicas first sample in their own knowledge‐based energy databases and then participate in the replica exchange operations in the REMD simulation. In fact, this method is a hybrid between the standard REMD method and the simulated tempering method. Using different active replicas, one can freely control the calculation quantity and the convergence speed of the simulation. To illustrate the performance of the method, we apply it to some small models. The distribution functions of the replicas in the energy space and temperature space show that the modified REMD method in this work can let the replicas walk freely in both of the two spaces. With the same number of the active replicas, the free energy surface in the simulation converges faster than the standard REMD. © 2016 Wiley Periodicals, Inc.     

A divide-and-conquer strategy to improve diffusion sampling in generalized ensemble simulations

Generalized ensemble simulations generally suffer from the associated diffusion-sampling problem; the increased entropic barrier can greatly abolish sampling efficiency, in particular, with the increase of number of degrees of freedom in the target conformational space. Taking advantage of the recent simulated scaling method, we formulate a divide-and-conquer sampling strategy to solve this problem so as to robustly improve the sampling efficiency in generalized ensemble simulations. In the present method, the target conformational space sampling enhancement is decomposed to the sampling enhancements of several subconformational regions, and multiple independent SS simulations are performed to establish the individual sampling enhancement for each of the subconformational regions; in order to realize the global importance sampling, structure exchanges among these replicas are performed based on the Monte Carlo acceptance/rejection procedure. As demonstrated in our studies, the present divide-and-conquer sampling algorithm, named by us as "simulated scaling based variant Hamiltonian replica exchange method," has superior sampling capability so as to possibly play an essential role in dealing with the present bottleneck of generalized ensemble method developments: the system size limitations.     

Hamiltonian replica‐exchange simulations with adaptive biasing of peptide backbone and side chain dihedral angles

A Hamiltonian Replica‐Exchange Molecular Dynamics (REMD) simulation method has been developed that employs a two‐dimensional backbone and one‐dimensional side chain biasing potential specifically to promote conformational transitions in peptides. To exploit the replica framework optimally, the level of the biasing potential in each replica was appropriately adapted during the simulations. This resulted in both high exchange rates between neighboring replicas and improved occupancy/flow of all conformers in each replica. The performance of the approach was tested on several peptide and protein systems and compared with regular MD simulations and previous REMD studies. Improved sampling of relevant conformational states was observed for unrestrained protein and peptide folding simulations as well as for refinement of a loop structure with restricted mobility of loop flanking protein regions. © 2013 Wiley Periodicals, Inc.     

Replica exchange with nonequilibrium switches: enhancing equilibrium sampling by increasing replica overlap

We describe a replica exchange strategy where trial swap configurations are generated by nonequilibrium switching simulations. By devoting simulation time to the switching simulations, one can systematically increase an effective overlap between replicas, which leads to an increased exchange acceptance rate and less correlated equilibrium samples. In this paper, we derive our method for a general class of stochastic dynamics, and discuss various strategies for enhancing replica overlap through novel dynamical schemes and prudent choices of switching protocols. We then demonstrate our method on a model system of alanine dipeptide in implicit solvent, characterizing decreases in data correlations and gains in sampling efficiency.     

Smart resolution replica exchange: an efficient algorithm for exploring complex energy landscapes

A coarse-grained representation of a condensed phase system can significantly reduce the number of system degrees of freedom, making coarse-grained simulations very computationally efficient. Moreover, coarse graining can smoothen the free energy landscape of the system. Thus coarse-grained dynamics is usually faster than its fully atomistic counterpart. In this work, the smart resolution replica exchange method is introduced that incorporates the information from coarse-grained simulations into atomistic simulations in order to accelerate the sampling of rough, complex atomistic energy landscapes. Within this methodology, interactions between particles are defined by a potential energy that interpolates between a fully atomistic potential and a fully coarse-grained effective potential according to a parameter lambda. Instead of exchanging the configurations from neighboring resolutions directly, as has been done in the resolution replica exchange methods [E. Lyman et al., Phys. Rev. Lett. 96, 028105 (2006); M. Christen and W. F. v. Gunsteren, J. Chem. Phys. 124, 154106 (2006)], the configuration described at the coarser resolution is first relaxed before an exchange is attempted, similar to the smart walking method [R. Zhou and B. J. Berne, J. Chem. Phys. 107, 9185 (1997)]. This approach greatly increases the acceptance ratio of exchange and only two replicas, one at the atomistic level and one at the coarse-grained level, are usually required (although more can be implemented if desired). This new method can approximately obtain the correct canonical sampling if the exchange interval is sufficiently large to allow the system to explore the local energy landscape. The method is demonstrated for a two-dimensional model system, where the ideal population distribution can be recovered, and also for an alanine polypeptide (Ala(15)) model with explicit water, where its native structure, an alpha helix, is obtained from the extended structure within 1 ns.     

Finite reservoir replica exchange to enhance canonical sampling in rugged energy surfaces

A "finite reservoir" replica exchange method is presented to further enhance sampling upon the regular replica exchange method (REM) in a rugged energy surface. The present method can facilitate important sampling more efficiently by exchanging structures with configurations randomly selected from a finite-sized reservoir; this finite reservoir is pregenerated and updated by a mechanism of replica exchange with neighboring "temperature" simulations. In practice, this proposal revises exchange schedule in REM simulations in order to make productive exchange for conformational "tunneling" more frequent.     

Lennard-Jones parameters for the combined QM/MM method using the B3LYP/6-31G*/AMBER potential

A combined DFT quantum mechanical and AMBER molecular mechanical potential (QM/MM) is presented for use in molecular modeling and molecular simulations of large biological systems. In our approach we evaluate Lennard-Jones parameters describing the interaction between the quantum mechanical (QM) part of a system, which is described at the B3LYP/6-31+G* level of theory, and the molecular mechanical (MM) part of the system, described by the AMBER force field. The Lennard-Jones parameters for this potential are obtained by calculating hydrogen bond energies and hydrogen bond geometries for a large set of bimolecular systems, in which one hydrogen bond monomer is described quantum mechanically and the other is treated molecular mechanically. We have investigated more than 100 different bimolecular systems, finding very good agreement between hydrogen bond energies and geometries obtained from the combined QM/MM calculations and results obtained at the QM level of theory, especially with respect to geometry. Therefore, based on the Lennard-Jones parameters obtained in our study, we anticipate that the B3LYP/6-31+G*/AMBER potential will be a precise tool to explore intermolecular interactions inside a protein environment.     

Adaptively biased molecular dynamics for free energy calculations

We present an adaptively biased molecular dynamics (ABMD) method for the computation of the free energy surface of a reaction coordinate using nonequilibrium dynamics. The ABMD method belongs to the general category of umbrella sampling methods with an evolving biasing potential and is inspired by the metadynamics method. The ABMD method has several useful features, including a small number of control parameters and an O(t) numerical cost with molecular dynamics time t. The ABMD method naturally allows for extensions based on multiple walkers and replica exchange, where different replicas can have different temperatures and/or collective variables. This is beneficial not only in terms of the speed and accuracy of a calculation, but also in terms of the amount of useful information that may be obtained from a given simulation. The workings of the ABMD method are illustrated via a study of the folding of the Ace-GGPGGG-Nme peptide in a gaseous and solvated environment.     

Two-dimensional replica exchange approach for peptide-peptide interactions

The replica exchange molecular dynamics (REMD) method has emerged as a standard approach for simulating proteins and peptides with rugged underlying free energy landscapes. We describe an extension to the original methodology--here termed umbrella-sampling REMD (UREMD)--that offers specific advantages in simulating peptide-peptide interactions. This method is based on the use of two dimensions in the replica cascade, one in temperature as in conventional REMD, and one in an umbrella sampling coordinate between the center of mass of the two peptides that aids explicit exploration of the complete association-dissociation reaction coordinate. To mitigate the increased number of replicas required, we pursue an approach in which the temperature and umbrella dimensions are linked at only fully associated and dissociated states. Coupled with the reweighting equations, the UREMD method aids accurate calculations of normalized free energy profiles and structural or energetic measures as a function of interpeptide separation distance. We test the approach on two families of peptides: a series of designed tetrapeptides that serve as minimal models for amyloid fibril formation, and a fragment of a classic leucine zipper peptide and its mutant. The results for these systems are compared to those from conventional REMD simulations, and demonstrate good convergence properties, low statistical errors, and, for the leucine zippers, an ability to sample near-native structures.     

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.     

Incorporation of a QM/MM buffer zone in the variational double self-consistent field method

The explicit polarization (X-Pol) potential is an electronic-structure-based polarization force field, designed for molecular dynamics simulations and modeling of biopolymers. In this approach, molecular polarization and charge transfer effects are explicitly treated by a combined quantum mechanical and molecular mechanical (QM/MM) scheme, and the wave function of the entire system is variationally optimized by a double self-consistent field (DSCF) method. In the present article, we introduce a QM buffer zone for a smooth transition from a QM region to an MM region. Instead of using the Mulliken charge approximation for all QM/MM interactions, the Coulombic interactions between the adjacent fragments are determined directly by electronic structure theory. The present method is designed to accelerate the speed of convergence of the total energy and charge density of the system.     

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.     

Quantum mechanics/molecular mechanics minimum free-energy path for accurate reaction energetics in solution and enzymes: sequential sampling and optimization on the potential of mean force surface

To accurately determine the reaction path and its energetics for enzymatic and solution-phase reactions, we present a sequential sampling and optimization approach that greatly enhances the efficiency of the ab initio quantum mechanics/molecular mechanics minimum free-energy path (QM/MM-MFEP) method. In the QM/MM-MFEP method, the thermodynamics of a complex reaction system is described by the potential of mean force (PMF) surface of the quantum mechanical (QM) subsystem with a small number of degrees of freedom, somewhat like describing a reaction process in the gas phase. The main computational cost of the QM/MM-MFEP method comes from the statistical sampling of conformations of the molecular mechanical (MM) subsystem required for the calculation of the QM PMF and its gradient. In our new sequential sampling and optimization approach, we aim to reduce the amount of MM sampling while still retaining the accuracy of the results by first carrying out MM phase-space sampling and then optimizing the QM subsystem in the fixed-size ensemble of MM conformations. The resulting QM optimized structures are then used to obtain more accurate sampling of the MM subsystem. This process of sequential MM sampling and QM optimization is iterated until convergence. The use of a fixed-size, finite MM conformational ensemble enables the precise evaluation of the QM potential of mean force and its gradient within the ensemble, thus circumventing the challenges associated with statistical averaging and significantly speeding up the convergence of the optimization process. To further improve the accuracy of the QM/MM-MFEP method, the reaction path potential method developed by Lu and Yang [Z. Lu and W. Yang, J. Chem. Phys. 121, 89 (2004)] is employed to describe the QM/MM electrostatic interactions in an approximate yet accurate way with a computational cost that is comparable to classical MM simulations. The new method was successfully applied to two example reaction processes, the classical SN2 reaction of Cl-+CH3Cl in solution and the second proton transfer step of the reaction catalyzed by the enzyme 4-oxalocrotonate tautomerase. The activation free energies calculated with this new sequential sampling and optimization approach to the QM/MM-MFEP method agree well with results from other simulation approaches such as the umbrella sampling technique with direct QM/MM dynamics sampling, demonstrating the accuracy of the iterative QM/MM-MFEP method.     

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.     

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.     

Hamiltonian replica exchange molecular dynamics using soft-core interactions

To overcome the problem of insufficient conformational sampling within biomolecular simulations, we have developed a novel Hamiltonian replica exchange molecular dynamics (H-REMD) scheme that uses soft-core interactions between those parts of the system that contribute most to high energy barriers. The advantage of this approach over other H-REMD schemes is the possibility to use a relatively small number of replicas with locally larger differences between the individual Hamiltonians. Because soft-core potentials are almost the same as regular ones at longer distances, most of the interactions between atoms of perturbed parts will only be slightly changed. Rather, the strong repulsion between atoms that are close in space, which in many cases results in high energy barriers, is weakened within higher replicas of our proposed scheme. In addition to the soft-core interactions, we proposed to include multiple replicas using the same Hamiltonian/level of softness. We have tested the new protocol on the GTP and 8-Br-GTP molecules, which are known to have high energy barriers between the anti and syn conformation of the base with respect to the sugar moiety. During two 25 ns MD simulations of both systems the transition from the more stable to the less stable (but still experimentally observed) conformation is not seen at all. Also temperature REMD over 50 replicas for 1 ns did not show any transition at room temperature. On the other hand, more than 20 of such transitions are observed in H-REMD using six replicas (at three different Hamiltonians) during 6.8 ns per replica for GTP and 12 replicas (at six different Hamiltonians) during 8.7 ns per replica for 8-Br-GTP. The large increase in sampling efficiency was obtained from an optimized H-REMD scheme involving soft-core potentials, with multiple simulations using the same level of softness. The optimization of the scheme was performed by fast mimicking [J. Hritz and C. Oostenbrink, J. Chem. Phys. 127, 204104 (2007)].