Convergence and sampling efficiency in replica exchange simulations of peptide folding in explicit solvent

Replica exchange methods (REMs) are increasingly used to improve sampling in molecular dynamics (MD) simulations of biomolecular systems. However, despite having been shown to be very effective on model systems, the application of REM in complex systems such as for the simulation of protein and peptide folding in explicit solvent has not been objectively tested in detail. Here we present a comparison of conventional MD and temperature replica exchange MD (T-REMD) simulations of a beta-heptapeptide in explicit solvent. This system has previously been shown to undergo reversible folding on the time scales accessible to MD simulation and thus allows a direct one-to-one comparison of efficiency. The primary properties compared are the free energy of folding and the relative populations of different conformers as a function of temperature. It is found that to achieve a similar degree of precision T-REMD simulations starting from a random set of initial configurations were approximately an order of magnitude more computationally efficient than a single 800 ns conventional MD simulation for this system at the lowest temperature investigated (275 K). However, whereas it was found that T-REMD simulations are more than four times more efficient than multiple independent MD simulations at one temperature (300 K) the actual increase in conformation sampling was only twofold. The overall gain in efficiency using REMD resulted primarily from the ordering of different conformational states over temperature, as opposed to a large increase of conformational sampling. It is also shown that in this system exchanges are accepted primarily based on (random) fluctuations within the solvent and are not strongly correlated with the instantaneous peptide conformation raising questions in regard to the efficiency of T-REMD in larger systems.     

Overcoming entropic barrier with coupled sampling at dual resolutions

An enhanced sampling method is proposed for ab initio protein folding simulations. The new method couples a high-resolution model for accuracy and a low-resolution model for efficiency. It aims to overcome the entropic barrier found in the exponentially large protein conformational space when a high-resolution model, such as an all-atom molecular mechanics force field, is used. The proposed method is designed to satisfy the detailed balance condition so that the Boltzmann distribution can be generated in all sampling trajectories in both high and low resolutions. The method was tested on model analytical energy functions and ab initio folding simulations of a beta-hairpin peptide. It was found to be more efficient than replica-exchange method that is used as its building block. Analysis with the analytical energy functions shows that the number of energy calculations required to find global minima and to converge mean potential energies is much fewer with the new method. Ergodic measure shows that the new method explores the conformational space more rapidly. We also studied imperfect low-resolution energy models and found that the introduction of errors in low-resolution models does decrease its sampling efficiency. However, a reasonable increase in efficiency is still observed when the global minima of the low-resolution models are in the vicinity of the global minimum basin of the high-resolution model. Finally, our ab initio folding simulation of the tested peptide shows that the new method is able to fold the peptide in a very short simulation time. The structural distribution generated by the new method at the equilibrium portion of the trajectory resembles that in the equilibrium simulation starting from the crystal structure.     

Efficient evaluation of sampling quality of molecular dynamics simulations by clustering of dihedral torsion angles and Sammon mapping

A direct conformational clustering and mapping approach for peptide conformations based on backbone dihedral angles has been developed and applied to compare conformational sampling of Met-enkephalin using two molecular dynamics (MD) methods. Efficient clustering in dihedrals has been achieved by evaluating all combinations resulting from independent clustering of each dihedral angle distribution, thus resolving all conformational substates. In contrast, Cartesian clustering was unable to accurately distinguish between all substates. Projection of clusters on dihedral principal component (PCA) subspaces did not result in efficient separation of highly populated clusters. However, representation in a nonlinear metric by Sammon mapping was able to separate well the 48 highest populated clusters in just two dimensions. In addition, this approach also allowed us to visualize the transition frequencies between clusters efficiently. Significantly, higher transition frequencies between more distinct conformational substates were found for a recently developed biasing-potential replica exchange MD simulation method allowing faster sampling of possible substates compared to conventional MD simulations. Although the number of theoretically possible clusters grows exponentially with peptide length, in practice, the number of clusters is only limited by the sampling size (typically much smaller), and therefore the method is well suited also for large systems. The approach could be useful to rapidly and accurately evaluate conformational sampling during MD simulations, to compare different sampling strategies and eventually to detect kinetic bottlenecks in folding pathways.     

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.     

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.     

Protecting High Energy Barriers: A New Equation to Regulate Boost Energy in Accelerated Molecular Dynamics Simulations

Molecular dynamics (MD) is one of the most common tools in computational chemistry. Recently, our group has employed accelerated molecular dynamics (aMD) to improve the conformational sampling over conventional molecular dynamics techniques. In the original aMD implementation, sampling is greatly improved by raising energy wells below a predefined energy level. Recently, our group presented an alternative aMD implementation where simulations are accelerated by lowering energy barriers of the potential energy surface. When coupled with thermodynamic integration simulations, this implementation showed very promising results. However, when applied to large systems, such as proteins, the simulation tends to be biased to high energy regions of the potential landscape. The reason for this behavior lies in the boost equation used since the highest energy barriers are dramatically more affected than the lower ones. To address this issue, in this work, we present a new boost equation that prevents oversampling of unfavorable high energy conformational states. The new boost potential provides not only better recovery of statistics throughout the simulation but also enhanced sampling of statistically relevant regions in explicit solvent MD simulations.     

Comprehensive Characterization of the Self-Folding Cavitand Dynamics

The conformational equilibria and guest exchange process of a resorcin[4]arene derived self-folding cavitand receptor have been characterized in detail by molecular dynamics simulations (MD) and ¹H EXSY NMR experiments. A multi-timescale strategy for exploring the fluxional behaviour of this system has been constructed, exploiting conventional MD and accelerated MD (aMD) techniques. The use of aMD allows the reconstruction of the folding/unfolding process of the receptor by sampling high-energy barrier processes unattainable by conventional MD simulations. We obtained MD trajectories sampling events occurring at different timescales from ns to s: 1) rearrangement of the directional hydrogen bond seam stabilizing the receptor, 2) folding/unfolding of the structure transiting partially open intermediates, and 3) guest departure from different folding stages. Most remarkably, reweighing of the biased aMD simulations provided kinetic barriers that are in very good agreement with those determined experimentally by ¹H NMR. These results constitute the first comprehensive characterization of the complex dynamic features of cavitand receptors. Our approach emerges as a valuable rational design tool for synthetic host-guest systems     

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.     

Molecular dynamics simulations of proteins and peptides: from folding to drug design

Computer simulations of proteins, lipids and nucleic acids at equilibrium have become essentially routine. However, the fact remains that complete sampling of conformational space continues to be a bottle-neck in the field. The challenge for the future is to overcome such problems and use computational approaches to understand recognition and spontaneous self-organization in biomolecular systems (folding, aggregation and assembly of complexes), processes that cannot be directly observed experimentally. In this review, examples illustrating the extent to which simulations can be used to understand these phenomena in biomolecular systems will be presented along with examples of methodological developments to increase our physical understanding of the processes. The study cases will cover the problems of peptide-receptor recognition and the use of the information obtained for the design of new non-peptidic ligands; the study of the folding mechanism of small proteins and finally the study of the initial stages of peptide self-aggregation.     

Serial replica exchange

Parallel tempering (or the replica exchange method (REM)) is a powerful method for speeding up the sampling of conformational states of systems with rough energy landscapes, like proteins, where stable conformational states can be separated by large energy barriers. The usual implementation of the REM is performed on local computer clusters (or parallel processors) where the different replicas must be run synchronously. Here, we present serial replica exchange (SREM), a method that is equivalent to the standard REM in terms of efficiency yet runs asynchronously on a distributed network of computers. A second advantage is the method's greatly enhanced fault tolerance, which enables the study of biological systems on worldwide distributed computing environments, such as Folding@Home. For proof of concept, we apply the SREM to a single alanine dipeptide molecule in explicit water. We show that the SREM reproduces the thermodynamic and structural properties determined by the REM.     

Adaptive lambda square dynamics simulation: An efficient conformational sampling method for biomolecules

A novel, efficient sampling method for biomolecules is proposed. The partial multicanonical molecular dynamics (McMD) was recently developed as a method that improved generalized ensemble (GE) methods to focus sampling only on a part of a system (GEPS); however, it was not tested well. We found that partial McMD did not work well for polylysine decapeptide and gave significantly worse sampling efficiency than a conventional GE. Herein, we elucidate the fundamental reason for this and propose a novel GEPS, adaptive lambda square dynamics (ALSD), which can resolve the problem faced when using partial McMD. We demonstrate that ALSD greatly increases the sampling efficiency over a conventional GE. We believe that ALSD is an effective method and is applicable to the conformational sampling of larger and more complicated biomolecule systems. © 2013 Wiley Periodicals, Inc.     

Force-momentum-based self-guided Langevin dynamics: a rapid sampling method that approaches the canonical ensemble

The self-guided Langevin dynamics (SGLD) is a method to accelerate conformational searching. This method is unique in the way that it selectively enhances and suppresses molecular motions based on their frequency to accelerate conformational searching without modifying energy surfaces or raising temperatures. It has been applied to studies of many long time scale events, such as protein folding. Recent progress in the understanding of the conformational distribution in SGLD simulations makes SGLD also an accurate method for quantitative studies. The SGLD partition function provides a way to convert the SGLD conformational distribution to the canonical ensemble distribution and to calculate ensemble average properties through reweighting. Based on the SGLD partition function, this work presents a force-momentum-based self-guided Langevin dynamics (SGLDfp) simulation method to directly sample the canonical ensemble. This method includes interaction forces in its guiding force to compensate the perturbation caused by the momentum-based guiding force so that it can approximately sample the canonical ensemble. Using several example systems, we demonstrate that SGLDfp simulations can approximately maintain the canonical ensemble distribution and significantly accelerate conformational searching. With optimal parameters, SGLDfp and SGLD simulations can cross energy barriers of more than 15 kT and 20 kT, respectively, at similar rates for LD simulations to cross energy barriers of 10 kT. The SGLDfp method is size extensive and works well for large systems. For studies where preserving accessible conformational space is critical, such as free energy calculations and protein folding studies, SGLDfp is an efficient approach to search and sample the conformational space.     

Calibration of force-field dependency in free energy landscapes of peptide conformations by quantum chemical calculations

The free energy landscapes of peptide conformations were calibrated by ab initio quantum chemical calculations, after the enhanced conformational diversity search using the multicanonical molecular dynamics simulations. Three different potentials of mean force for an isolated dipeptide were individually obtained by the multicanonical molecular dynamics simulations using the conventional force fields, AMBER parm94, AMBER parm96, and CHARMm22. Each potential of mean force was then calibrated based upon the umbrella sampling algorithm from the adiabatic energy map that was calculated separately by the ab initio molecular orbital method, and all of the calibrated potentials of mean force coincided well. The calibration method was also applied to the simulations of a peptide dimer in explicit water models, and it was shown that the calibrated free energy landscapes did not depend on the force field used in the classical simulations, as far as the conformational space was sampled well. The current calibration method fuses the classical free energy calculation with the quantum chemical calculation, and it should generally make simulations for biomolecular systems much more reliable when combining with enhanced conformational sampling.     

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.     

Simple continuous and discrete models for simulating replica exchange simulations of protein folding

The efficiency of temperature replica exchange (RE) simulations hinge on their ability to enhance conformational sampling at physiological temperatures by taking advantage of more rapid conformational interconversions at higher temperatures. While temperature RE is a parallel simulation technique that is relatively straightforward to implement, kinetics in the RE ensemble is complicated, and there is much to learn about how best to employ RE simulations in computational biophysics. Protein folding rates often slow down above a certain temperature due to entropic bottlenecks. This "anti-Arrhenius" behavior represents a challenge for RE. However, it is far from straightforward to systematically explore the impact of this on RE by brute force molecular simulations, since RE simulations of protein folding are very difficult to converge. To understand some of the basic mechanisms that determine the efficiency of RE, it is useful to study simplified low dimensionality systems that share some of the key characteristics of molecular systems. Results are presented concerning the efficiency of temperature RE on a continuous two-dimensional potential that contains an entropic bottleneck. Optimal efficiency was obtained when the temperatures of the replicas did not exceed the temperature at which the harmonic mean of the folding and unfolding rates is maximized. This confirms a result we previously obtained using a discrete network model of RE. Comparison of the efficiencies obtained using the continuous and discrete models makes it possible to identify non-Markovian effects, which slow down equilibration of the RE ensemble on the more complex continuous potential. In particular, the rate of temperature diffusion and also the efficiency of RE is limited by the time scale of conformational rearrangements within free energy basins.     

Enhanced conformational sampling method for proteins based on the TaBoo SeArch algorithm: Application to the folding of a mini‐protein,chignolin

The conformational samplings are indispensible for obtaining reliable canonical ensembles, which provide statistical averages of physical quantities such as free energies. However, the samplings of vast conformational space of biomacromolecules by conventional molecular dynamics (MD) simulations might be insufficient, due to their inadequate accessible time‐scales for investigating biological functions. Therefore, the development of methodologies for enhancing the conformational sampling of biomacromolecules still remains as a challenging issue in computational biology. To tackle this problem, we newly propose an efficient conformational search method, which is referred as TaBoo SeArch (TBSA) algorithm. In TBSA, an inverse energy histogram is used to select seeds for the conformational resampling so that states with high frequencies are inhibited, while states with low frequencies are efficiently sampled to explore the unvisited conformational space. As a demonstration, TBSA was applied to the folding of a mini‐protein, chignolin, and automatically sampled the native structure (C_α root mean square deviation < 1.0 Å) with nanosecond order computational costs started from a completely extended structure, although a long‐time 1‐µs normal MD simulation failed to sample the native structure. Furthermore, a multiscale free energy landscape method based on the conformational sampling of TBSA were quantitatively evaluated through free energy calculations with both implicit and explicit solvent models, which enable us to find several metastable states on the folding landscape. © 2015 Wiley Periodicals, Inc.     

New-generation amber united-atom force field

We have developed a new-generation Amber united-atom force field for simulations involving highly demanding conformational sampling such as protein folding and protein-protein binding. In the new united-atom force field, all hydrogens on aliphatic carbons in all amino acids are united with carbons except those on Calpha. Our choice of explicit representation of all protein backbone atoms aims at minimizing perturbation to protein backbone conformational distributions and to simplify development of backbone torsion terms. Tests with dipeptides and solvated proteins show that our goal is achieved quite successfully. The new united-atom force field uses the same new RESP charging scheme based on B3LYP/cc-pVTZ//HF/6-31g** quantum mechanical calculations in the PCM continuum solvent as that in the Duan et al. force field. van der Waals parameters are empirically refitted starting from published values with respect to experimental solvation free energies of amino acid side-chain analogues. The suitability of mixing new point charges and van der Waals parameters with existing Amber covalent terms is tested on alanine dipeptide and is found to be reasonable. Parameters for all new torsion terms are refitted based on the new point charges and the van der Waals parameters. Molecular dynamics simulations of three small globular proteins in the explicit TIP3P solvent are performed to test the overall stability and accuracy of the new united-atom force field. Good agreements between the united-atom force field and the Duan et al. all-atom force field for both backbone and side-chain conformations are observed. In addition, the per-step efficiency of the new united-atom force field is demonstrated for simulations in the implicit generalized Born solvent. A speedup around two is observed over the Duan et al. all-atom force field for the three tested small proteins. Finally, the efficiency gain of the new united-atom force field in conformational sampling is further demonstrated with a well-known toy protein folding system, an 18 residue polyalanine in distance-dependent dielectric. The new united-atom force field is at least a factor of 200 more efficient than the Duan et al. all-atom force field for ab initio folding of the tested peptide.     

Reliable treatment of electrostatics in combined QM/MM simulation of macromolecules

A robust approach for dealing with electrostatic interactions for spherical boundary conditions has been implemented in the QM/MM framework. The development was based on the generalized solvent boundary potential (GSBP) method proposed by Im et al. [J. Chem. Phys. 114, 2924 (2001)], and the specific implementation was applied to the self-consistent-charge density-functional tight-binding approach as the quantum mechanics (QM) level, although extension to other QM methods is straightforward. Compared to the popular stochastic boundary-condition scheme, the new protocol offers a balanced treatment between quantum mechanics/molecular mechanics (QM/MM) and MM/MM interactions; it also includes the effect of the bulk solvent and macromolecule atoms outside of the microscopic region at the Poisson-Boltzmann level. The new method was illustrated with application to the enzyme human carbonic anhydrase II and compared to stochastic boundary-condition simulations using different electrostatic treatments. The GSBP-based QM/MM simulations were most consistent with available experimental data, while conventional stochastic boundary simulations yielded various artifacts depending on different electrostatic models. The results highlight the importance of carefully treating electrostatics in QM/MM simulations of biomolecules and suggest that the commonly used truncation schemes should be avoided in QM/MM simulations, especially in simulations that involve extensive conformational samplings. The development of the GSBP-based QM/MM protocol has opened up the exciting possibility of studying chemical events in very complex biomolecular systems in a multiscale framework.     

Enhanced Lipid Diffusion and Mixing in Accelerated Molecular Dynamics

Accelerated molecular dynamics (aMD) is an enhanced sampling technique that expedites conformational space sampling by reducing the barriers separating various low-energy states of a system. Here, we present the first application of the aMD method on lipid membranes. Altogether, ～1.5 μs simulations were performed on three systems: a pure POPC bilayer, a pure DMPC bilayer, and a mixed POPC:DMPC bilayer. Overall, the aMD simulations are found to produce significant speedup in trans-gauche isomerization and lipid lateral diffusion versus those in conventional MD (cMD) simulations. Further comparison of a 70-ns aMD run and a 300-ns cMD run of the mixed POPC:DMPC bilayer shows that the two simulations yield similar lipid mixing behaviors, with aMD generating a 2-3-fold speedup compared to cMD. Our results demonstrate that the aMD method is an efficient approach for the study of bilayer structural and dynamic properties. On the basis of simulations of the three bilayer systems, we also discuss the impact of aMD parameters on various lipid properties, which can be used as a guideline for future aMD simulations of membrane systems.     

Application of torsion angle molecular dynamics for efficient sampling of protein conformations

We investigate the application of torsion angle molecular dynamics (TAMD) to augment conformational sampling of peptides and proteins. Interesting conformational changes in proteins mainly involve torsional degrees of freedom. Carrying out molecular dynamics in torsion space does not only explicitly sample the most relevant degrees of freedom, but also allows larger integration time steps with elimination of the bond and angle degrees of freedom. However, the covalent geometry needs to be fixed during internal coordinate dynamics, which can introduce severe distortions to the underlying potential surface in the extensively parameterized modern Cartesian-based protein force fields. A "projection" approach (Katritch et al. J Comput Chem 2003, 24, 254-265) is extended to construct an accurate internal coordinate force field (ICFF) from a source Cartesian force field. Torsion crossterm corrections constructed from local molecular fragments, together with softened van der Waals and electrostatic interactions, are used to recover the potential surface and incorporate implicit bond and angle flexibility. MD simulations of dipeptide models demonstrate that full flexibility in both the backbone phi/psi and side chain chi1 angles are virtually restored. The efficacy of TAMD in enhancing conformational sampling is then further examined by folding simulations of small peptides and refinement experiments of protein NMR structures. The results show that an increase of several fold in conformational sampling efficiency can be reliably achieved. The current study also reveals some complicated intrinsic properties of internal coordinate dynamics, beyond energy conservation, that can limit the maximum size of the integration time step and thus the achievable gain in sampling efficiency.