Reservoir pH replica exchange.

We present the reservoir pH replica exchange (R-pH-REM) method for constant pH simulations. The R-pH-REM method consists of a two-step procedure; the first step involves generation of one or more reservoirs of conformations. Each reservoir is obtained from a standard or enhanced molecular dynamics simulation with a constrained (fixed) protonation state. In the second step, fixed charge constraints are relaxed, as the structures from one or more reservoirs are periodically injected into a constant pH or a pH-replica exchange (pH-REM) simulation. The benefit of this two-step process is that the computationally intensive part of conformational search can be decoupled from constant pH simulations, and various techniques for enhanced conformational sampling can be applied without the need to integrate such techniques into the pH-REM framework. Simulations on blocked Lys, KK, and KAAE peptides were used to demonstrate an agreement between pH-REM and R-pH-REM simulations. While the reservoir simulations are not needed for these small test systems, the real need arises in cases when ionizable molecules can sample two or more conformations separated by a large energy barrier, such that adequate sampling is not achieved on a time scale of standard constant pH simulations. Such problems might be encountered in protein systems that exploit conformational transitions for function. A hypothetical case is studied, a small molecule with a large torsional barrier; while results of pH-REM simulations depend on the starting structure, R-pH-REM calculations on this model system are in excellent agreement with a theoretical model.

MDMS: Molecular Dynamics Meta-Simulator for evaluating exchange type sampling methods.

Replica exchange methods have become popular tools to explore conformational space for small proteins. For larger biological systems, even with enhanced sampling methods, exploring the free energy landscape remains computationally challenging. This problem has led to the development of many improved replica exchange methods. Unfortunately, testing these methods remains expensive. We propose a Molecular Dynamics Meta-Simulator (MDMS) based on transition state theory to simulate a replica exchange simulation, eliminating the need to run explicit dynamics between exchange attempts. MDMS simulations allow for rapid testing of new replica exchange based methods, greatly reducing the amount of time needed for new method development.

Evaluating the performance of the ff99SB force field based on NMR scalar coupling data.

Force-field validation is essential for the identification of weaknesses in current models and the development of more accurate models of biomolecules. NMR coupling and relaxation methods have been used to effectively diagnose the strengths and weaknesses of many existing force fields. Studies using the ff99SB force field have shown excellent agreement between experimental and calculated order parameters and residual dipolar calculations. However, recent studies have suggested that ff99SB demonstrates poor agreement with J-coupling constants for short polyalanines. We performed extensive replica-exchange molecular-dynamics simulations on Ala(3) and Ala(5) in TIP3P and TIP4P-Ew solvent models. Our results suggest that the performance of ff99SB is among the best of currently available models. In addition, scalar coupling constants derived from simulations in the TIP4P-Ew model show a slight improvement over those obtained using the TIP3P model. Despite the overall excellent agreement, the data suggest areas for possible improvement.

The unfolded state of the villin headpiece helical subdomain: computational studies of the role of locally stabilized structure.

The 36 residue villin headpiece helical subdomain (HP36) is one of the fastest cooperatively folding proteins, folding on the microsecond timescale. HP36's simple three helix topology, fast folding and small size have made it an attractive model system for computational and experimental studies of protein folding. Recent experimental studies have explored the denatured state of HP36 using fragment analysis coupled with relatively low-resolution spectroscopic techniques. These studies have shown that there is apparently only a small tendency to form locally stabilized secondary structure. Here, we complement the experimental studies by using replica exchange molecular dynamics with explicit solvent to investigate the structural features of these peptide models of unfolded HP36. To ensure convergence, two sets of simulations for each fragment were performed with different initial structures, and simulations were continued until these generated very similar final ensembles. These simulations reveal low populations of native-like structure and early folding events that cannot be resolved by experiment. For each fragment, calculated J-coupling constants and helical propensities are in good agreement with experimental trends. HP-1, corresponding to residues 41 to 53 and including the first alpha-helix, contains the highest helical population. HP-3, corresponding to residues 62 through 75 and including the third alpha-helix, contains a small population of helical turn residing at the N terminus while HP-2, corresponding to residues 52 through 61 and including the second alpha-helix, formed little to no structure in isolation. Overall, HP-1 was the only fragment to adopt a native-like conformation, but the low population suggests that formation of significant structure only occurs after formation of specific tertiary interactions.

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.

Secondary structure bias in generalized Born solvent models: comparison of conformational ensembles and free energy of solvent polarization from explicit and implicit solvation.

The effects of the use of three generalized Born (GB) implicit solvent models on the thermodynamics of a simple polyalanine peptide are studied via comparing several hundred nanoseconds of well-converged replica exchange molecular dynamics (REMD) simulations using explicit TIP3P solvent to REMD simulations with the GB solvent models. It is found that when compared to REMD simulations using TIP3P the GB REMD simulations contain significant differences in secondary structure populations, most notably an overabundance of alpha-helical secondary structure. This discrepancy is explored via comparison of the differences in the electrostatic component of the free energy of solvation (DeltaDeltaG(pol)) between TIP3P (via thermodynamic Integration calculations), the GB models, and an implicit solvent model based on the Poisson equation (PE). The electrostatic components of the solvation free energies are calculated using each solvent model for four representative conformations of Ala10. Since the PE model is found to have the best performance with respect to reproducing TIP3P DeltaDeltaG(pol) values, effective Born radii from the GB models are compared to effective Born radii calculated with PE (so-called perfect radii), and significant and numerous deviations in GB radii from perfect radii are found in all GB models. The effect of these deviations on the solvation free energy is discussed, and it is shown that even when perfect radii are used the agreement of GB with TIP3P DeltaDeltaG(pol) values does not improve. This suggests a limit to the optimization of the effective Born radius calculation and that future efforts to improve the accuracy of GB models must extend beyond such optimizations.

Comparison of multiple Amber force fields and development of improved protein backbone parameters.

The ff94 force field that is commonly associated with the Amber simulation package is one of the most widely used parameter sets for biomolecular simulation. After a decade of extensive use and testing, limitations in this force field, such as over-stabilization of alpha-helices, were reported by us and other researchers. This led to a number of attempts to improve these parameters, resulting in a variety of "Amber" force fields and significant difficulty in determining which should be used for a particular application. We show that several of these continue to suffer from inadequate balance between different secondary structure elements. In addition, the approach used in most of these studies neglected to account for the existence in Amber of two sets of backbone phi/psi dihedral terms. This led to parameter sets that provide unreasonable conformational preferences for glycine. We report here an effort to improve the phi/psi dihedral terms in the ff99 energy function. Dihedral term parameters are based on fitting the energies of multiple conformations of glycine and alanine tetrapeptides from high level ab initio quantum mechanical calculations. The new parameters for backbone dihedrals replace those in the existing ff99 force field. This parameter set, which we denote ff99SB, achieves a better balance of secondary structure elements as judged by improved distribution of backbone dihedrals for glycine and alanine with respect to PDB survey data. It also accomplishes improved agreement with published experimental data for conformational preferences of short alanine peptides and better accord with experimental NMR relaxation data of test protein systems.

Generating reservoir conformations for replica exchange through the use of the conformational space annealing method.

Temperature replica exchange molecular dynamics (T-REM) has been successfully used to improve the conformational search for model peptides and small proteins. However, for larger and more complicated systems, the use of T-REM is computationally intensive since the complexity of the free energy landscape and number of replicas required increase with system size. Achieving convergence of systems with slow transition kinetics is often difficult. Several methods have been proposed to overcome the size and convergence speed issues of standard T-REM. One of these is the Reservoir Replica Exchange Method (R-REM), in which the conformational search and temperature equilibration are separated by exchanging with a pre-existing reservoir of structures. This approach allows the integration of computationally efficient search algorithms with replica exchange. The Conformational Space Annealing (CSA) method has been shown to be able to determine the global energy minimum of proteins efficiently and has been used in structure prediction successfully. CSA uses a genetic algorithm to generate a diverse set of conformations to determine the minimum energy structure. We combine these methods by using conformations generated by the CSA method to build a reservoir. R-REM is then used to seed the top replica with the structures from the reservoir; fast convergence at every temperature is observed. The efficiency of this method is then demonstrated with model peptides and small proteins, and significant improvement of efficiency is observed while maintaining the overall shape of the free energy landscape.

Improving the description of salt bridge strength and geometry in a Generalized Born model.

The Generalized Born (GB) solvent model is widely used in molecular dynamics simulations because it can be less computationally expensive and it samples conformational changes more efficiently than explicit solvent simulations. Meanwhile, great efforts have been made in the past to improve its precision and accuracy. Previous studies have shown that reducing intrinsic GB radii of some hydrogen atoms would improve AMBER GB-HCT solvent model's accuracy on salt bridges. Here we present our finding that similar correction also shows dramatic improvement for the AMBER GB-OBC solvent model. Potential of mean force and cluster analysis for small peptide replica exchange molecular dynamics simulations suggested that new radii GB simulation with ff99SB/GB-OBC corrected salt bridge strength and achieved significantly higher geometry similarity with TIP3P simulation. Improved performance in 60 ns HIV-1 protease GB simulation further validated this approach for large systems.

MSCALE: A General Utility for Multiscale Modeling.

The combination of theoretical models of macromolecules that exist at different spatial and temporal scales has become increasingly important for addressing complex biochemical problems. This work describes the extension of concurrent multiscale approaches, introduces a general framework for carrying out calculations, and describes its implementation into the CHARMM macromolecular modeling package. This functionality, termed MSCALE, generalizes both the additive and subtractive multiscale scheme (e.g. QM/MM ONIOM-type), and extends its support to classical force fields, coarse grained modeling (e.g. ENM, GNM, etc.), and a mixture of them all. The MSCALE scheme is completely parallelized with each subsystem running as an independent, but connected calculation. One of the most attractive features of MSCALE is the relative ease of implementation using the standard MPI communication protocol. This allows external access to the framework and facilitates the combination of functionality previously isolated in separate programs. This new facility is fully integrated with free energy perturbation methods, Hessian based methods, and the use of periodicity and symmetry, which allows the calculation of accurate pressures. We demonstrate the utility of this new technique with four examples; (1) subtractive QM/MM and QM/QM calculations; (2) multi-force field alchemical free energy perturbation; (3) integration with the SANDER module of AMBER and the TINKER package to gain access to potentials not available in CHARMM; and (4) mixed resolution (i.e. coarse grain / all-atom) normal mode analysis. The potential of this new tool is clearly established and in conclusion an interesting mathematical problem is highlighted and future improvements are proposed.

Evaluation of Salt Bridge Structure and Energetics in Peptides Using Explicit, Implicit, and Hybrid Solvation Models.

Replica exchange or parallel tempering molecular dynamics (REMD) is widely used to enhance the exploration of free energy landscapes for complex molecular systems. However its application to large systems is hampered by the scaling of the number of required replicas with an increasing system size. We recently proposed an improved REMD method where the exchange probabilities were calculated using a hybrid explicit/implicit solvent model. We previously tested this hybrid solvent REMD approach on alanine polypeptides of 1, 3, and 10 residues and obtained very good agreement with fully solvated REMD simulations while significantly reducing the number of replicas required. In this study we continue evaluating the applicability of the hybrid solvent REMD method through comparing the free energy of formation of ion pairs using model peptides. In accord with other studies, pure GB simulations resulted in overstabilized salt bridges, whereas the hybrid models produced free energy profiles in close agreement with fully solvated simulations, including solvent separated minima. Furthermore, the structure of the salt bridge in explicit solvent is reproduced by the hybrid solvent REMD method, while the GB simulations favor a different geometry.

Improving Convergence of Replica-Exchange Simulations through Coupling to a High-Temperature Structure Reservoir.

Parallel tempering or replica-exchange molecular dynamics (REMD) significantly increases the efficiency of conformational sampling for complex molecular systems. However, obtaining converged data with REMD remains challenging, especially for large systems with complex topologies. We propose a new variant to REMD where the replicas are also permitted to exchange with an ensemble of structures that have been generated in advance using high-temperature MD simulations, similar in spirit to J-walking methods. We tested this approach on two non-trivial model systems, a ß-hairpin and a 3-stranded ß-sheet and compared the results to those obtained from very long (>100 ns) standard REMD simulations. The resulting ensembles were indistinguishable, including relative populations of different conformations on the unfolded state. The use of the reservoir is shown to significantly reduce the time required for convergence.

HIV-1 protease flaps spontaneously open and reclose in molecular dynamics simulations.

We report unrestrained, all-atom molecular dynamics simulations of HIV-1 protease that sample large conformational changes of the active site flaps. In particular, the unliganded protease undergoes multiple conversions between the "closed" and "semiopen" forms observed in crystal structures of inhibitor-bound and unliganded protease, respectively, including reversal of flap "handedness." Simulations in the presence of a cyclic urea inhibitor yield stable closed flaps. Furthermore, we observe several events in which the flaps of the unliganded protease open to a much greater degree than observed in crystal structures and subsequently return to the semiopen state. Our data strongly support the hypothesis that the unliganded protease predominantly populates the semiopen conformation, with closed and fully open structures being a minor component of the overall ensemble. The results also provide a model for the flap opening and closing that is considered to be essential to enzyme function.

HIV-1 protease flaps spontaneously close to the correct structure in simulations following manual placement of an inhibitor into the open state.

We report unrestrained, all-atom molecular dynamics simulations of HIV-1 protease (HIV-PR) with a continuum solvent model that reproducibly sample closing of the active site flaps following manual placement of a cyclic urea inhibitor into the substrate binding site of the open protease. The open form was obtained from the unbound, semi-open HIV-PR crystal structure, which we recently reported (Hornak, V.; et al. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 915-920.) to have spontaneously opened during unrestrained dynamics. In those simulations, the transiently open flaps always returned to the semi-open form that is observed in all crystal structures of the free protease. Here, we show that manual docking of the inhibitor reproducibly induces spontaneous conversion to the closed form as seen in all inhibitor-bound HIV-PR crystal structures. These simulations reproduced not only the greater degree of flap closure, but also the striking difference in flap "handedness" between bound and free enzyme. In most of the simulations, the final structures were highly accurate. Root-mean-square deviations (RMSD) from the crystal structure of the complex were approximately 1.5 A (averaged over the last 100 ps) for the inhibitor and each flap despite initial RMSD of 2-5 A for the inhibitors and 6-11 A for the flaps. Key hydrogen bonds were formed between the flap tips and between flaps and inhibitor that match those seen in the crystal structure. The results demonstrate that all-atom simulations have the ability to significantly improve poorly docked ligand conformations and reproduce large-scale receptor conformational changes that occur upon binding.

Improved Efficiency of Replica Exchange Simulations through Use of a Hybrid Explicit/Implicit Solvation Model.

The use of parallel tempering or replica exchange molecular dynamics (REMD) simulations has facilitated the exploration of free energy landscapes for complex molecular systems, but application to large systems is hampered by the scaling of the number of required replicas with increasing system size. Use of continuum solvent models reduces system size and replica requirements, but these have been shown to provide poor results in many cases, including overstabilization of ion pairs and secondary structure bias. Hybrid explicit/continuum solvent models can overcome some of these problems through an explicit representation of water molecules in the first solvation shells, but these methods typically require restraints on the solvent molecules and show artifacts in water properties due to the solvation interface. We propose an REMD variant in which the simulations are performed with a fully explicit solvent, but the calculation of exchange probability is carried out using a hybrid model, with the solvation shells calculated on the fly during the fully solvated simulation. The resulting reduction in the perceived system size in the REMD exchange calculation provides a dramatic decrease in the computational cost of REMD, while maintaining a very good agreement with results obtained from the standard explicit solvent REMD. We applied several standard and hybrid REMD methods with different solvent models to alanine polymers of 1, 3, and 10 residues, obtaining ensembles that were essentially independent of the initial conformation, even with explicit solvation. Use of only a continuum model without a shell of explicit water provided poor results for Ala3 and Ala10, with a significant bias in favor of the α-helix. Likewise, using only the solvation shells and no continuum model resulted in ensembles that differed significantly from the standard explicit solvent data. Ensembles obtained from hybrid REMD are in very close agreement with explicit solvent data, predominantly populating polyproline II conformations. Inclusion of a second shell of explicit solvent was found to be unnecessary for these peptides.

Using PC clusters to evaluate the transferability of molecular mechanics force fields for proteins.

The transferability of molecular mechanics parameters derived for small model systems to larger biopolymers such as proteins can be difficult to assess. Even for small peptides, molecular dynamics simulations are typically too short to sample structures significantly different than initial conformations, making comparison to experimental data questionable. We employed a PC cluster to generate large numbers of native and non-native conformations for peptides with experimentally measured structural data, one predominantly helical and the other forming a beta-hairpin. These atomic-detail sets do not suffer from slow convergence, and can be used to rapidly evaluate important force field properties. In this case a suspected bias toward alpha-helical conformations in the ff94 and ff99 force fields distributed with the AMBER package was verified. The sets provide critical feedback not only on force field transferability, but may also predict modifications for improvement. Such predictions were used to modify the ff99 parameter set, and the resulting force field was used to test stability and folding of model peptides. Structural behavior during molecular dynamics with the modified force field is found to be very similar to expectations, suggesting that these basis sets of conformations may themselves have significant transferability among force fields. We continue to improve and expand this data set and plan to make it publicly accessible. The calculations involved in this process are trivially parallel and can be performed using inexpensive personal computers with commodity components.