Reversible folding simulation by hybrid Hamiltonian replica exchange

Reversible foldings of a beta-hairpin peptide, chignolin, by recently invented hybrid Hamiltonian replica exchange molecular dynamics simulations based on Poisson-Boltzmann model in explicit water are demonstrated. Initiated from extended structures the peptide folded and unfolded a couple of times in seven out of eight replica trajectories during 100 nanoseconds simulation. The folded states have the lowest all-atom root mean squared deviation of 1.3 A with respect to the NMR structures. At T=300 K the occurrence of folded states was converged to 62% during 80 ns simulation which agrees well with experimental data. Especially, a detailed structural evolution map was constructed based on 800,000 structural snapshots and from where a unique folding doorway emerges. Compared with 130 ns standard replica exchange simulation using 24 replicas on the same system, the hybrid Hamiltonian replica exchange molecular dynamics simulation presents consistent results.     

Computational study on the structural diversity of amyloid Beta Peptide (abeta(10-35)) oligomers

We studied the oligomerization of Alzheimer amyloid beta peptide (Abeta) using a replica exchange molecular dynamics (REMD) simulation. The simulation was performed with Abeta(10-35) dimers, trimers, and tetramers. Extensive REMD simulations illustrated several possible oligomer conformations. As the size of the oligomer increased from a dimer to a tetramer, the number of possible configurations was reduced. We identified all the possible conformations for each oligomer and characterized their temperature dependence. It was found that the detailed structures of the oligomers, which may act as folding intermediates, are highly sensitive to the parameters of the simulation environment such as temperature and concentration. Structural diversities of Abeta oligomers suggest multiple pathways of the aggregation process.     

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.     

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)].     

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.     

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.     

Influence of water-protein hydrogen bonding on the stability of Trp-cage miniprotein. A comparison between the TIP3P and TIP4P-Ew water models

We report extensive replica exchange molecular dynamics (REMD) simulations on the folding/unfolding equilibrium of Trp-cage miniprotein using the Amber ff99SB all atom forcefield and TIP3P and TIP4P-Ew explicit water solvent models. REMD simulation-lengths in the 500 ns to the microsecond regime per replica are required to adequately sample the folding/unfolding equilibrium. We observe that this equilibrium is significantly affected by the choice of the water model. Compared with experimental data, simulations using the TIP3P solvent describe the stability of the Trp-cage quite realistically, providing a melting point which is just a few Kelvins above the experimental transition temperature of 317 K. The TIP4P-Ew model shifts the equilibrium towards the unfolded state and lowers the free energy of unfolding by about 3 kJ mol(-1) at 280 K, demonstrating the need to fine-tune the protein-forcefield depending on the chosen water model. We report evidence that the main difference between the two water models is mostly due to the different solvation of polar groups of the peptide. The unfolded state of the Trp-cage is stabilized by an increasing number of hydrogen bonds, destabilizing the α-helical part of the molecule and opening the R-D salt bridge. By reweighting the strength of solvent-peptide hydrogen bonds by adding a hydrogen bond square well potential, we can fully recover the effect of the different water models and estimate the shift in population as due to a difference in hydrogen bond-strength of about 0.4 kJ mol(-1) per hydrogen bond.     

Hybrid Hamiltonian replica exchange molecular dynamics simulation method employing the Poisson-Boltzmann model

A hybrid Hamiltonian replica exchange molecular dynamics simulation scheme based on explicit water model hybrided with Poisson-Boltzmann model is brought out. In this method the motions of atoms are governed by potential energy obtained from explicit water model. However, the exchanges between different replicas under different temperatures are controlled by the solvation energies of the solute calculated using the Poisson-Boltzmann model. In order to get the correct canonical ensembles, the van der Waals radii, which are used to define the dielectric boundary, have to be optimized. The conformational spaces of three distinct pentapeptides, Met-enkephalin, alanine 5, and glycine 5, are explored. We find that with the optimized radii the structural ensembles are nearly identical to those obtained by standard replica exchange simulations while the number of replica needed is reduced greatly.     

Simulating temperature jumps for protein folding

We present a new computational methodology aimed to calculate the thermodynamics and kinetics of peptide folding. We focus in particular on temperature jump experiments of folding rates and show how a combination of replica exchange molecular dynamics (REMD) followed by multiplexed molecular dynamics starting from structures taken from the REMD runs can be used to extract properties in line with experiments. A model system, alanine20, is studied in this article as a proof of principle and used to describe the methodology.     

Velocity-scaling optimized replica exchange molecular dynamics of proteins in a hybrid explicit/implicit solvent

We propose a scheme for replica exchange molecular dynamics of proteins in explicit solvent that minimizes the number of required replicas using velocity rescaling. Our approach relies on a hybrid method where the protein evolves at each temperature in an explicit solvent, but replica exchange moves utilize an implicit solvent term. The two terms are coupled through the velocity rescaling. We test the efficiency of this approach for a common test case, the trp-cage protein.     

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.     

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.     

Solvophobic and steric effects of side groups on polymer folding: molecular modeling studies of amine-functionalized m-poly(phenyleneethynylene) foldamers in aqueous solution

The folding behavior of five different amine-functionalized m-poly(phenyleneethynylene) (m-PPE) oligomers containing 24 phenyl rings (12 residues, where a residue includes 2 phenyl rings) in water was examined by using a combination of molecular dynamics (MD) and replica exchange molecular dynamics (REMD) simulation techniques. The REMD method employed the highly parallelized GROMACS MD software and a modified OPLS-AA force field to simulate 44 replicas of each solvated system in parallel, with temperatures ranging from 300 to 577 K. Our results showed that the REMD method was more effective in predicting the helical conformation of the m-PPE in water, from an extended structure, than canonical MD methods in the same simulation time. Furthermore, we observed from canonical MD simulations of the explicitly solvated helical m-PPEs at 300 K that the radius of gyration, average helix inner diameter, and average helix pitch of the helical structure all pass through a minima when the side group is R = OC(2)H(5) as R is changed from R = H through OC(4)H(9).     

A bias-exchange approach to protein folding

By suitably extending a recent approach [Bussi, G.; et al. J. Am. Chem. Soc. 2006, 128, 13435] we introduce a powerful methodology that allows the parallel reconstruction of the free energy of a system in a virtually unlimited number of variables. Multiple metadynamics simulations of the same system at the same temperature are performed, biasing each replica with a time-dependent potential constructed in a different set of collective variables. Exchanges between the bias potentials in the different variables are periodically allowed according to a replica exchange scheme. Due to the efficaciously multidimensional nature of the bias the method allows exploring complex free energy landscapes with high efficiency. The usefulness of the method is demonstrated by performing an atomistic simulation in explicit solvent of the folding of a Triptophane cage miniprotein. It is shown that the folding free energy landscape can be fully characterized starting from an extended conformation with use of only 40 ns of simulation on 8 replicas.     

On the modeling of aggregates of an optically active regioregular polythiophene

The conformational properties of the optically active regioregular poly[(R)-3-(4-(4-ethyl-2-oxazolin-2-yl) phenyl) thiophene] (PEOPT) were explored by molecular dynamics on a single chain using several solvents of increasing polarity. Furthermore, their aggregate formation was studied over a wide range of temperatures using a replica exchange molecular dynamics simulation providing simulation data representative of the equilibrium behaviour of their aggregates. Results show a clear tendency of PEOPT to keep a syn-gauche conformation between continuous backbone thiophene rings favouring a bent chain structure in solvent. After studying their aggregation behaviour in acetonitrile, a strong tendency to pack stabilizing structures that reinforce the chirality of the polymer, in concordance with experimental data, was found. Two different aggregated structures were observed depending on oligomer length, a self-assembled helical aggregate based on stacked octamers and a bent double helix aggregate in large oligomers.     

Folding processes of the B domain of protein A to the native state observed in all-atom ab initio folding simulations

Reaching the native states of small proteins, a necessary step towards a comprehensive understanding of the folding mechanisms, has remained a tremendous challenge to ab initio protein folding simulations despite the extensive effort. In this work, the folding process of the B domain of protein A (BdpA) has been simulated by both conventional and replica exchange molecular dynamics using AMBER FF03 all-atom force field. Started from an extended chain, a total of 40 conventional (each to 1.0 micros) and two sets of replica exchange (each to 200.0 ns per replica) molecular dynamics simulations were performed with different generalized-Born solvation models and temperature control schemes. The improvements in both the force field and solvent model allowed successful simulations of the folding process to the native state as demonstrated by the 0.80 A C(alpha) root mean square deviation (RMSD) of the best folded structure. The most populated conformation was the native folded structure with a high population. This was a significant improvement over the 2.8 A C(alpha) RMSD of the best nativelike structures from previous ab initio folding studies on BdpA. To the best of our knowledge, our results demonstrate, for the first time, that ab initio simulations can reach the native state of BdpA. Consistent with experimental observations, including Phi-value analyses, formation of helix II/III hairpin was a crucial step that provides a template upon which helix I could form and the folding process could complete. Early formation of helix III was observed which is consistent with the experimental results of higher residual helical content of isolated helix III among the three helices. The calculated temperature-dependent profile and the melting temperature were in close agreement with the experimental results. The simulations further revealed that phenylalanine 31 may play critical to achieve the correct packing of the three helices which is consistent with the experimental observation. In addition to the mechanistic studies, an ab initio structure prediction was also conducted based on both the physical energy and a statistical potential. Based on the lowest physical energy, the predicted structure was 2.0 A C(alpha) RMSD away from the experimentally determined structure.     

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.     

Molecular dynamics simulations of supramolecular polymer rheology

Using equilibrium and nonequilibrium molecular dynamics simulations, we studied the equilibrium and rheological properties of dilute and semidilute solutions of head-to-tail associating polymers. In our simulation model, a spontaneous complementary reversible association between the donor and the acceptor groups at the ends of oligomers was achieved by introducing a combination of truncated pseudo-Coulombic attractive potential and Lennard Jones repulsive potential between donor, acceptor, and neighboring groups. We have calculated the equilibrium properties of supramolecular polymers, such as the ring/chain equilibrium, average molecular weight, and molecular weight distribution of self-assembled chains and rings, which all agree well with previous analytical and computer modeling results. We have investigated shear thinning of solutions of 8- and 20-bead associating oligomers with different association energies at different temperatures and oligomer volume fractions. All reduced viscosity data for a given oligomer length can be collapsed into one master curve, exhibiting two power-law regions of shear-thinning behavior with an exponent of -0.55 at intermediate ranges of the reduced shear rate β and -0.8 (or -0.9) at larger shear rates. The equilibrium viscosity of supramolecular solutions with different oligomer lengths and associating energies is found to obey a power-law scaling dependence on oligomer volume fraction with an exponent of 1.5, in agreement with the experimental observations for several dilute or semidilute solutions of supramolecular polymers. This implies that dilute and semidilute supramolecular polymer solutions exhibit high polydispersity but may not be sufficiently entangled to follow the reptation mechanism of relaxation.     

Temperature weighted histogram analysis method, replica exchange, and transition paths

We analyzed the data from a replica exchange molecular dynamics simulation using the weighted histogram analysis method to combine data from all of the temperature replicas (T-WHAM) to obtain the room-temperature potential of mean force of the G-peptide (the C-terminal beta-hairpin of the B1 domain of protein G) in regions of conformational space not sampled at room temperature. We were able to determine the potential of mean force in the transition region between a minor alpha-helical population and the major beta-hairpin population and identify a possible transition path between them along which the peptide retains a significant amount of secondary structure. This observation provides new insights into a possible mechanism of formation of beta-sheet secondary structures in proteins. We developed a novel Bayesian statistical uncertainty estimation method for any quantity derived from WHAM and used it to validate the calculated potential of mean force. The feasibility of estimating regions of the potential of mean force with unfavorable free energy at room temperature by T-WHAM analysis of replica exchange simulations was further tested on a system that can be solved analytically and presented some of the same challenges found in more complex chemical systems.