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.     

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.     

Folding simulations with novel conformational search method

A novel scheme for fast conformational search has been developed by combining the replica exchange method (REM) with the generalized effective potential concept. The new method, referred to Q-REM [S. Jang et al. Phys. Rev. Lett. 91, 058305 (2003)], is expected to provide a useful alternative to the conventional REM for effective conformational sampling of complex systems. The authors have performed folding simulations of the Trp-cage miniprotein using Q-REM. All atom level simulations with generalized Born solvent access-area solvation model show that successful folding can be observed with much smaller number of replicas in Q-REM compared to the conventional REM. It can be concluded that the new method has potential to significantly improve sampling efficiency, allowing simulations of more challenging systems.     

Convergence of replica exchange molecular dynamics

Replica exchange molecular dynamics (REMD) method is one of the generalized-ensemble algorithms which performs random walk in energy space and helps a system to escape from local energy traps. In this work, we studied the accuracy and efficiency of REMD by examining its ability to reproduce the results of multiple extended conventional molecular dynamics (MD) simulations and to enhance conformational sampling. Two sets of REMD simulations with different initial configurations, one from the fully extended and the other from fully helical conformations, were conducted on a fast-folding 21-amino-acid peptide with a continuum solvent model. Remarkably, the two REMD simulation sets started to converge even within 1.0 ns, despite their dramatically different starting conformations. In contrast, the conventional MD within the same time and with identical starting conformations did not show obvious signs of convergence. Excellent convergence between the REMD sets for T>300 K was observed after 14.0 ns REMD simulations as measured by the average helicity and free-energy profiles. We also conducted a set of 45 MD simulations at nine different temperatures with each trajectory simulated to 100.0 and 200.0 ns. An excellent agreement between the REMD and the extended MD simulation results was observed for T>300 K, showing that REMD can accurately reproduce long-time MD results with high efficiency. The autocorrelation times of the calculated helicity demonstrate that REMD can significantly enhance the sampling efficiency by 14.3+/-6.4, 35.1+/-0.2, and 71.5+/-20.4 times at, respectively, approximately 360, approximately 300, and approximately 275 K in comparison to the regular MD. Convergence was less satisfactory at low temperatures (T<300 K) and a slow oscillatory behavior suggests that longer simulation time was needed to reach equilibrium. Other technical issues, including choice of exchange frequency, were also examined.     

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.     

On the Use of a Supramolecular Coarse‐Grained Model for the Solvent in Simulations of the Folding Equilibrium of an Octa‐β‐peptide in MeOH and H2O

Molecular dynamics (MD) simulation can give a detailed picture of conformational equilibria of biomolecules, but it is only reliable if the force field used in the simulation is accurate, and the sampling of the conformational space accessible to the biomolecule shows many (un)folding transitions to allow for precise averages of observable quantities. Here, the use of coarse‐grained (CG) solvent MeOH and H₂O models to speed up the sampling of the conformational equilibria of an octa‐β‐peptide is investigated. This peptide is thought to predominantly adopt a 3₁₄‐helical fold when solvated in MeOH, and a hairpin fold when solvated in H₂O on the basis of the NMR data. Various factors such as the chirality of a residue, a force‐field modification for the solute, coarse‐graining of the solvent model, and an extension of the nonbonded interaction cut‐off radius are shown to influence the simulated conformational equilibria and the agreement with the experimental NMR data for the octa‐β‐peptide.     

The beta-peptide hairpin in solution: conformational study of a beta-hexapeptide in methanol by NMR spectroscopy and MD simulation.

The structural and thermodynamic properties of a 6-residue beta-peptide that was designed to form a hairpin conformation have been studied by NMR spectroscopy and MD simulation in methanol solution. The predicted hairpin would be characterized by a 10-membered hydrogen-bonded turn involving residues 3 and 4, and two extended antiparallel strands. The interproton distances and backbone torsional dihedral angles derived from the NMR experiments at room temperature are in general terms compatible with the hairpin conformation. Two trajectories of system configurations from 100-ns molecular-dynamics simulations of the peptide in solution at 298 and 340 K have been analyzed. In both simulations reversible folding to the hairpin conformation is observed. Interestingly, there is a significant conformational overlap between the unfolded state of the peptide at each of the temperatures. As already observed in previous studies of peptide folding, the unfolded state is composed of a (relatively) small number of predominant conformers and in this case lacks any type of secondary-structure element. The trajectories provide an excellent ground for the interpretation of the NMR-derived data in terms of ensemble averages and distributions as opposed to single-conformation interpretations. From this perspective, a relative population of the hairpin conformation of 20% to 30% would suffice to explain the NMR-derived data. Surprisingly, however, the ensemble of structures from the simulation at 340 K reproduces more accurately the NMR-derived data than the ensemble from the simulation at 298 K, a question that needs further investigation.     

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     

Monte Carlo vs molecular dynamics for all-atom polypeptide folding simulations

An efficient Monte Carlo (MC) algorithm including concerted rotations is directly compared to molecular dynamics (MD) in all-atom statistical mechanics folding simulations of small polypeptides. The previously reported algorithm "concerted rotations with flexible bond angles" (CRA) has been shown to successfully locate the native state of small polypeptides. In this study, the folding of three small polypeptides (trpzip2/H1/Trp-cage) is investigated using MC and MD, for a combined sampling time of approximately 10(11) MC configurations and 8 micros, respectively. Both methods successfully locate the experimentally determined native states of the three systems, but they do so at different speed, with 2-2.5 times faster folding of the MC runs. The comparison reveals that thermodynamic and dynamic properties can reliably be obtained by both and that results from folding simulations do not depend on the algorithm used. Similar to previous comparisons of MC and MD, it is found that one MD integration step of 2 fs corresponds to one MC scan, revealing the good sampling of MC. The simplicity and efficiency of the MC method will enable its future use in folding studies involving larger systems and the combination with replica exchange algorithms.     

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.     

关于副本交换分子动力学模拟复杂化学反应的研究

本文研究用温度副本交换分子动力学(T-REMD)和哈密顿副本交换分子动力学(H-REMD)方法模拟复杂化学反应的问题。使用具有不同活化能和反应能的简单置换反应模型,我们检验了上述两种方法用来预测反应平衡产物的效率和应用范围。T-REMD方法对具有适度活化能(约< 20 kcal·mol^-1)或者反应能量(< 3 kcal·mol^-1)的放热反应是有效的。由于在相空间的不完整采样,对于同时具有高活化能和反应能量的反应其模拟效率有严重障碍,并且对于吸热反应问题更为显着。另一方面,H-REMD对一系列具有不同活化能的反应能的模型表现出色,与T-REMD相比,H-REMD可以使用更少的副本获得优异的结果。     

All-atom level direct folding simulation of a betabetaalpha miniprotein

We performed ab initio folding simulation for a betabetaalpha peptide BBA5 (PDB code 1T8J) with a modified param99 force field using the generalized Born solvation model (param99MOD5/GBSA). For efficient conformational sampling, we extended a previously developed novel Q-replica exchange molecular dynamics (Q-REMD) into a multiplexed Q-REMD. Starting from a fully extended conformation, we were able to locate the nativelike structure in the global free minimum region at 280 K. The current approach, which combines the more balanced force field with the efficient sampling scheme, demonstrates a clear advantage in direct folding simulation at all-atom level.     

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.     

Replica exchanging self-guided Langevin dynamics for efficient and accurate conformational sampling

This work presents a replica exchanging self-guided Langevin dynamics (RXSGLD) simulation method for efficient conformational searching and sampling. Unlike temperature-based replica exchanging simulations, which use high temperatures to accelerate conformational motion, this method uses self-guided Langevin dynamics (SGLD) to enhance conformational searching without the need to elevate temperatures. A RXSGLD simulation includes a series of SGLD simulations, with simulation conditions differing in the guiding effect and∕or temperature. These simulation conditions are called stages and the base stage is one with no guiding effect. Replicas of a simulation system are simulated at the stages and are exchanged according to the replica exchanging probability derived from the SGLD partition function. Because SGLD causes less perturbation on conformational distribution than high temperatures, exchanges between SGLD stages have much higher probabilities than those between different temperatures. Therefore, RXSGLD simulations have higher conformational searching ability than temperature based replica exchange simulations. Through three example systems, we demonstrate that RXSGLD can generate target canonical ensemble distribution at the base stage and achieve accelerated conformational searching. Especially for large systems, RXSGLD has remarkable advantages in terms of replica exchange efficiency, conformational searching ability, and system size extensiveness.     

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.     

Prediction of Protein Loop Conformations using the AGBNP Implicit Solvent Model and Torsion Angle Sampling

The OPLS-AA all-atom force field and the Analytical Generalized Born plus Non-Polar (AGBNP) implicit solvent model, in conjunction with torsion angle conformational search protocols based on the Protein Local Optimization Program (PLOP), are shown to be effective in predicting the native conformations of 57 9-residue and 35 13-residue loops of a diverse series of proteins with low sequence identity. The novel nonpolar solvation free energy estimator implemented in AGBNP augmented by correction terms aimed at reducing the occurrence of ion pairing are important to achieve the best prediction accuracy. Extended versions of the previously developed PLOP-based conformational search schemes based on calculations in the crystal environment are reported that are suitable for application to loop homology modeling without the crystal environment. Our results suggest that in general the loop backbone conformation is not strongly influenced by crystal packing. The application of the temperature Replica Exchange Molecular Dynamics (T-REMD) sampling method for a few examples where PLOP sampling is insufficient are also reported. The results reported indicate that the OPLS-AA/AGBNP effective potential is suitable for high-resolution modeling of proteins in the final stages of homology modeling and/or protein crystallographic refinement.     

Molecular dynamics simulations of helix-forming, amine-functionalized m-poly(phenyleneethynylene)s

We present here the results of all-atom and united-atom molecular dynamics (MD) simulations that were used to examine the folding behavior of an amine-functionalized m-poly(phenyleneethynylene) (m-PPE) oligomer in aqueous environment. The parallelized GROMACS MD simulation code and OPLS force field were used for multiple MD simulations of m-PPE oligomers containing 24 phenyl rings in extended, coiled and helix conformations separately in water to determine the minimum energy conformation of the oligomer in aqueous solvent and what interactions are most important in determining this structure. Simulation results showed that the helix is the preferred minimum energy conformation of a single oligomer in water and that Lennard-Jones interactions are the dominant forces for the stabilization of the helix. In addition, these solvophobic interactions are strong enough to maintain the helix conformation at temperatures up to 523 K.     

Using enveloping distribution sampling to compute the free enthalpy difference between right- and left-handed helices of a β-peptide in solution

Recently, the method of enveloping distribution sampling (EDS) to efficiently obtain free enthalpy differences between different molecular systems from a single simulation has been generalized to compute free enthalpy differences between different conformations of a system [Z. X. Lin, H. Y. Liu, S. Riniker, and W. F. van Gunsteren, J. Chem. Theory Comput. 7, 3884 (2011)]. However, the efficiency of EDS in this case is hampered if the parts of the conformational space relevant to the two end states or conformations are far apart and the conformational diffusion from one state to the other is slow. This leads to slow convergence of the EDS parameter values and free enthalpy differences. In the present work, we apply the EDS methodology to a challenging case, i.e., to calculate the free enthalpy difference between a right-handed 2.7(10∕12)-helix and a left-handed 3(14)-helix of a hexa-β-peptide in solution from a single simulation. No transition between the two helices was detected in a standard EDS parameter update simulation, thus enhanced sampling techniques had to be applied, which included adiabatic decoupling (AD) of solute and solvent motions in combination with increasing the solute temperature, and lowering the shear viscosity of the solvent. AD was found to be unsuitable to enhance the sampling of the solute conformations in the EDS parameter update simulations. Lowering the solvent shear viscosity turned out to be useful during EDS parameter update simulations, i.e., it did speed up the conformational diffusion of the solute, more transitions between the two helices were observed. This came at the cost of more CPU time spent due to the shorter time step needed for simulations with the lower solvent shear viscosity. Using an improved EDS parameter update scheme, parameter convergence was five-fold enhanced. The resulting free enthalpy difference between the two helices calculated from EDS agrees well with the result obtained through direct counting from a long MD simulation, while the EDS technique significantly enhances the sampling of both helices over non-helical conformations.     

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.     

Folding of the 25 residue Abeta(12-36) peptide in TFE/water: temperature-dependent transition from a funneled free-energy landscape to a rugged one

The free-energy landscape of the Alzheimer beta-amyloid peptide Abeta(12-36) in a 40% (v/v) 2,2,2-trifluoroethanol (TFE)/water solution was determined by using multicanonical molecular dynamics simulations. Simulations using this enhanced conformational sampling technique were initiated from a random unfolded polypeptide conformation. Our simulations reliably folded the peptide to the experimental NMR structure, which consists of two linked helices. The shape of the free energy landscape for folding was found to be strongly dependent on temperature: Above 325 K, the overall shape was funnel-like, with the bottom of the funnel coinciding exactly with the NMR structure. Below 325 K, on the other hand, the landscape became increasingly rugged, with the emergence of new conformational clusters connected by low free-energy pathways. Finally, our simulations reveal that water and TFE solvate the polypeptide in different ways: The hydrogen bond formation between TFE and Abeta was enhanced with decreasing temperature, while that between water and Abeta was depressed.