Adaptively biased molecular dynamics for free energy calculations

A boundary correction algorithm for metadynamics in multiple dimensions

Metadynamics is an efficient method for simulation of the free energy of many-particle systems. Over the last several years it has been applied to study a wide variety of systems, ranging from simple fluids to biological macromolecules. The method relies on uniform sampling along specified collective variables or order parameters. Such order parameters, however, are often bounded, and metadynamics algorithms as originally developed suffer from systematic errors at the corresponding boundaries. While several approaches have been proposed in the past to correct these errors for unidimensional systems, no method exists to fully correct these errors in multi-dimensional systems at points where multiple boundaries meet. Here we present a correction scheme that circumvents this limitation.

Unbinding free energy of acetylcholinesterase bound oxime drugs along the gorge pathway from metadynamics-umbrella sampling investigation

Because of the pivotal role that the nerve enzyme, acetylcholinesterase plays in terminating nerve impulses at cholinergic synapses. Its active site, located deep inside a 20 Å gorge, is a vulnerable target of the lethal organophosphorus compounds. Potent reactivators of the intoxicated enzyme are nucleophiles, such as bispyridinium oxime that binds to the peripheral anionic site and the active site of the enzyme through suitable cation-π interactions. Atomic scale molecular dynamics and free energy calculations in explicit water are used to study unbinding pathways of two oxime drugs (Ortho-7 and Obidoxime) from the gorge of the enzyme. The role of enzyme-drug cation-π interactions are explored with the metadynamics simulation. The metadynamics discovered potential of mean force (PMF) of the unbinding events is refined by the umbrella sampling (US) corrections. The bidimensional free energy landscape of the metadynamics runs are further subjected to finite temperature string analysis to obtain the transition tube connecting the minima and bottlenecks of the unbinding pathway. The PMF is also obtained from US simulations using the biasing potential constructed from the transition tube and are found to be consistent with the metadynamics-US corrected results. Although experimental structural data clearly shows analogous coordination of the two drugs inside the gorge in the bound state, the PMF of the drug trafficking along the gorge pathway point, within an equilibrium free energy context, to a multistep process that differs from one another. Routes, milestones and subtlety toward the unbinding pathway of the two oximes at finite temperature are identified.

Perspectives for hybrid ab initio/molecular mechanical simulations of solutions: from complex chemistry to proton-transfer reactions and interfaces

As a consequence of the ongoing development of enhanced computational resources, theoretical chemistry has become an increasingly valuable field for the investigation of a variety of chemical systems. Simulations employing a hybrid quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) technique have been shown to be a particularly promising approach, whenever ultrafast (i.e., picosecond) dynamical properties are to be studied, which are in many cases difficult to access via experimental techniques. Details of the quantum mechanical charge field (QMCF) ansatz, an advanced QM/MM protocol, are discussed and simulation results for various systems ranging from simple ionic hydrates to solvated organic molecules and coordination complexes in solution are presented. A particularly challenging application is the description of proton-transfer reactions in chemical simulations, which is a prerequisite to study acidified and basic systems. The methodical requirements for a combination of the QMCF methodology with a dissociative potential model for the description of the solvent are discussed. Furthermore, the possible extension of QM/MM approaches to solid/liquid interfaces is outlined.

Roadmaps through free energy landscapes calculated using the multi-dimensional vFEP approach

The variational free energy profile (vFEP) method is extended to two dimensions and tested with molecular simulation applications. The proposed 2D-vFEP approach effectively addresses the two major obstacles to constructing free energy profiles from simulation data using traditional methods: the need for overlap in the re-weighting procedure and the problem of data representation. This is especially evident as these problems are shown to be more severe in two dimensions. The vFEP method is demonstrated to be highly robust and able to provide stable, analytic free energy profiles with only a paucity of sampled data. The analytic profiles can be analyzed with conventional search methods to easily identify stationary points (e.g. minima and first-order saddle points) as well as the pathways that connect these points. These "roadmaps" through the free energy surface are useful not only as a post-processing tool to characterize mechanisms, but can also serve as a basis from which to direct more focused "on-the-fly" sampling or adaptive force biasing. Test cases demonstrate that 2D-vFEP outperforms other methods in terms of the amount and sparsity of the data needed to construct stable, converged analytic free energy profiles. In a classic test case, the two dimensional free energy profile of the backbone torsion angles of alanine dipeptide, 2D-vFEP needs less than 1% of the original data set to reach a sampling accuracy of 0.5 kcal/mol in free energy shifts between windows. A new software tool for performing one and two dimensional vFEP calculations is herein described and made publicly available.

Unraveling the molecular structure of the catalytic domain of matrix metalloproteinase-2 in complex with a triple-helical peptide by means of molecular dynamics simulations

Herein, we present the results of a computational study that employed various simulation methodologies to build and validate a series of molecular models of a synthetic triple-helical peptide (fTHP-5) both in its native state and in a prereactive complex with the catalytic domain of the MMP-2 enzyme. First, the structure and dynamical properties of the fTHP-5 substrate are investigated by means of molecular dynamics (MD) simulations. Then, the propensity of each of the three peptide chains in fTHP-5 to be distorted around the scissile peptide bond is assessed by carrying out potential of mean force calculations. Subsequently, the distorted geometries of fTHP-5 are docked within the MMP-2 active site following a semirigid protocol, and the most stable docked structures are fully relaxed and characterized by extensive MD simulations in explicit solvent. Following a similar approach, we also investigate a hypothetical ternary complex formed between two MMP-2 catalytic units and a single fTHP-5 molecule. Overall, our models for the MMP-2/fTHP-5 complexes unveil the extent to which the triple helix is distorted to allow the accommodation of an individual peptide chain within the MMP active site.

Hidden Conformation Events in DNA Base Extrusions: A Generalized Ensemble Path Optimization and Equilibrium Simulation Study

DNA base extrusion is a crucial component of many biomolecular processes. Elucidating how bases are selectively extruded from the interiors of double-strand DNAs is pivotal to accurately understanding and efficiently sampling this general type of conformational transitions. In this work, the on-the-path random walk (OTPRW) method, which is the first generalized ensemble sampling scheme designed for finite-temperature-string path optimizations, was improved and applied to obtain the minimum free energy path (MFEP) and the free energy profile of a classical B-DNA major-groove base extrusion pathway. Along the MFEP, an intermediate state and the corresponding transition state were located and characterized. The MFEP result suggests that a base-plane-elongation event rather than the commonly focused base-flipping event is dominant in the transition state formation portion of the pathway; and the energetic penalty at the transition state is mainly introduced by the stretching of the Watson-Crick base pair. Moreover to facilitate the essential base-plane-elongation dynamics, the surrounding environment of the flipped base needs to be intimately involved. Further taking the advantage of the extended-dynamics nature of the OTPRW Hamiltonian, an equilibrium generalized ensemble simulation was performed along the optimized path; and based on the collected samples, several base-flipping (opening) angle collective variables were evaluated. In consistence with the MFEP result, the collective variable analysis result reveals that none of these commonly employed flipping (opening) angles alone can adequately represent the base extrusion pathway, especially in the pre-transition-state portion. As further revealed by the collective variable analysis, the base-pairing partner of the extrusion target undergoes a series of in-plane rotations to facilitate the base-plane-elongation dynamics. A base-plane rotation angle is identified to be a possible reaction coordinate to represent these in-plane rotations. Notably, these in-plane rotation motions may play a pivotal role in determining the base extrusion selectivity.

Driven Metadynamics: Reconstructing Equilibrium Free Energies from Driven Adaptive-Bias Simulations

We present a novel free-energy calculation method that constructively integrates two distinct classes of nonequilibrium sampling techniques, namely, driven (e.g., steered molecular dynamics) and adaptive-bias (e.g., metadynamics) methods. By employing nonequilibrium work relations, we design a biasing protocol with an explicitly time- and history-dependent bias that uses on-the-fly work measurements to gradually flatten the free-energy surface. The asymptotic convergence of the method is discussed, and several relations are derived for free-energy reconstruction and error estimation. Isomerization reaction of an atomistic polyproline peptide model is used to numerically illustrate the superior efficiency and faster convergence of the method compared with its adaptive-bias and driven components in isolation.

Solvent interactions stabilize the polyproline II conformation of glycosylated oligoprolines

In nature, proline residues carry several post-translational modifications (PTMs), including 4R hydroxylation and glycosylation. A recent study synthesized contiguously hydroxylated and glycosylated nonaproline peptides and revealed that both PTMs lead to a significant increase in the thermal stability of PPII relative to the unmodified oligoproline. The increased stability of the hydroxylated peptide can be explained by increased stability of the trans isomer due to stereoelectronic effects. However, the effects of glycosylation cannot be completely explained by stereoelectronics since previous experimental results indicate that 4R-glycosylation does not produce observable changes in the trans preference compared to 4R-hydroxylation. We therefore used sophisticated molecular modeling techniques to determine the reason for the further increase in thermal stability upon glycosylation. Free energy estimates obtained from adaptively biased molecular dynamics calculations in implicit (explicit) solvent are -9 kcal mol(-1) (-20 kcal mol(-1)) for the hydroxylated compound and -9 kcal mol(-1) (-46 kcal mol(-1)) for the glycosylated compound, indicating that direct solvent-peptide interactions are vital for explaining the glycosylation effects on PPII stability. Our data reveals for the first time that interactions between the hydroxyl groups in the glycosylated compound and water act in a complementary fashion with stereoelectronic effects to stabilize the PPII conformation in these substituted oligoproline peptides.

Simulation of catalytic water activation in mitochondrial F1-ATPase using a hybrid quantum mechanics/molecular mechanics approach: an alternative role for β-Glu 188

The use of quantum mechanics/molecular mechanics simulations to study the free energy landscape of the water activation at the catalytic site of mitochondrial F(1)-ATPase affords us insight into the generation of the nucleophile OH(-) prior to ATP hydrolysis. As a result, the ATP molecule was found to be the final proton acceptor. In the simulated pathway, the transfer of a proton to the nucleotide was not direct but occurred via a second water molecule in a manner similar to the Grotthuss mechanism proposed for proton diffusion. Residue β-Glu 188, previously described as the putative catalytic base, was found to be involved in the stabilization of a transient hydronium ion during water activation. Simulations in the absence of the carboxylate moiety of β-Glu 188 support this role.

Recipes for free energy calculations in biomolecular systems

During the last decade, several methods for sampling phase space and calculating various free energies in biomolecular systems have been devised or refined for molecular dynamics (MD) simulations. Thus, state-of-the-art methodology and the ever increasing computer power allow calculations that were forbidden a decade ago. These calculations, however, are not trivial as they require knowledge of the methods, insight into the system under study, and, quite often, an artful combination of different methodologies in order to avoid the various traps inherent in an unknown free energy landscape. In this chapter, we illustrate some of these concepts with two relatively simple systems, a sugar ring and proline oligopeptides, whose free energy landscapes still offer considerable challenges. In order to explore the configurational space of these systems, and to surmount the various free energy barriers, we combine three complementary methods: a nonequilibrium umbrella sampling method (adaptively biased MD, or ABMD), replica-exchange molecular dynamics (REMD), and steered molecular dynamics (SMD). In particular, ABMD is used to compute the free energy surface of a set of collective variables; REMD is used to improve the performance of ABMD, to carry out sampling in space complementary to the collective variables, and to sample equilibrium configurations directly; and SMD is used to study different transition mechanisms.

A New Maximum Likelihood Approach for Free Energy Profile Construction from Molecular Simulations

A novel variational method for construction of free energy profiles from molecular simulation data is presented. The variational free energy profile (VFEP) method uses the maximum likelihood principle applied to the global free energy profile based on the entire set of simulation data (e.g from multiple biased simulations) that spans the free energy surface. The new method addresses common obstacles in two major problems usually observed in traditional methods for estimating free energy surfaces: the need for overlap in the re-weighting procedure and the problem of data representation. Test cases demonstrate that VFEP outperforms other methods in terms of the amount and sparsity of the data needed to construct the overall free energy profiles. For typical chemical reactions, only ~5 windows and ~20-35 independent data points per window are sufficient to obtain an overall qualitatively correct free energy profile with sampling errors an order of magnitude smaller than the free energy barrier. The proposed approach thus provides a feasible mechanism to quickly construct the global free energy profile and identify free energy barriers and basins in free energy simulations via a robust, variational procedure that determines an analytic representation of the free energy profile without the requirement of numerically unstable histograms or binning procedures. It can serve as a new framework for biased simulations and is suitable to be used together with other methods to tackle with the free energy estimation problem.

A sugar isomerization reaction established on various (βα)₈-barrel scaffolds is based on substrate-assisted catalysis

In the course of tryptophan biosynthesis, the isomerization of phosphoribosylanthranilate (PRA) is catalyzed by the (βα)₈-barrel enzyme TrpF. The reaction occurs via a general acid-base mechanism with an aspartate and a cysteine residue acting as acid and base, respectively. PRA isomerase activity could be established on two (βα)₈-barrel enzymes involved in histidine biosynthesis, namely HisA and HisF, and on a HisAF chimera, by introducing two aspartate-to-valine substitutions. We have analyzed the reaction mechanism underlying this engineered activity by measuring its pH dependence, solving the crystal structure of a HisF variant with bound product analogue, and applying molecular dynamics simulations and mixed quantum and molecular mechanics calculations. The results suggest that PRA is anchored by the C-terminal phosphate-binding sites of HisA, HisF and HisAF. As a consequence, a conserved aspartate residue, which is equivalent to Cys7 from TrpF, is properly positioned to act as catalytic base. However, no obvious catalytic acid corresponding to Asp126 from TrpF could be identified in the three proteins. Instead, this role appears to be carried out by the carboxylate group of the anthranilate moiety of PRA. Thus, the engineered PRA isomerization activity is based on a reaction mechanism including substrate-assisted catalysis and thus differs substantially from the naturally evolved reaction mechanism used by TrpF.

Reaction path ensemble of the B-Z-DNA transition: a comprehensive atomistic study

Since its discovery in 1979, left-handed Z-DNA has evolved from an in vitro curiosity to a challenging DNA structure with crucial roles in gene expression, regulation and recombination. A fundamental question that has puzzled researchers for decades is how the transition from B-DNA, the prevalent right-handed form of DNA, to Z-DNA is accomplished. Due to the complexity of the B–Z-DNA transition, experimental and computational studies have resulted in several different, apparently contradictory models. Here, we use molecular dynamics simulations coupled with state-of-the-art enhanced sampling techniques operating through non-conventional reaction coordinates, to investigate the B–Z-DNA transition at the atomic level. Our results show a complex free energy landscape, where several phenomena such as over-stretching, unpeeling, base pair extrusion and base pair flipping are observed resulting in interconversions between different DNA conformations such as B-DNA, Z-DNA and S-DNA. In particular, different minimum free energy paths allow for the coexistence of different mechanisms (such as zipper and stretch–collapse mechanisms) that previously had been proposed as independent, disconnected models. We find that the B–Z-DNA transition—in absence of other molecular partners—can encompass more than one mechanism of comparable free energy, and is therefore better described in terms of a reaction path ensemble.

Enhanced sampling of molecular dynamics simulation of peptides and proteins by double coupling to thermal bath

Low sampling efficiency in conformational space is the well-known problem for conventional molecular dynamics. It greatly increases the difficulty for molecules to find the transition path to native state, and costs amount of CPU time. To accelerate the sampling, in this paper, we re-couple the critical degrees of freedom in the molecule to environment temperature, like dihedrals in generalized coordinates or nonhydrogen atoms in Cartesian coordinate. After applying to ALA dipeptide model, we find that this modified molecular dynamics greatly enhances the sampling behavior in the conformational space and provides more information about the state-to-state transition, while conventional molecular dynamics fails to do so. Moreover, from the results of 16 independent 100 ns simulations by the new method, it shows that trpzip2 has one-half chances to reach the naive state in all the trajectories, which is greatly higher than conventional molecular dynamics. Such an improvement would provide a potential way for searching the conformational space or predicting the most stable states of peptides and proteins.

Statistical analysis of ion mobility spectrometry. II. Adaptively biased methods and shape correlations

Following a recent effort [J. Am. Soc. Mass Spectrom. 23, 386-396 (2012)], we continue to explore computational methodologies for generating molecular conformations to support collisional cross sections suggested by ion mobility measurements. Here, adaptively biased molecular dynamics (ABMD) simulations are used to sample the configuration space and to achieve flat-histogram sampling along the reaction coordinates of the first two moments of the gyration tensor. The method is tested and compared with replica-exchange simulations on triply-protonated bradykinin and on a larger 25-residue peptide. It is found to have a much higher efficiency for producing large sets of conformations in a broad range of diffusion cross-sections, whereas it does not compete with conventional replica-exchange molecular dynamics in locating the lowest-energy structure. Nevertheless, the broad sampling obtained from the ABMD method allows to quantitatively correlate the diffusion cross-section Ω with other geometric order parameters that have simpler interpretation. The strong correlations found between the diffusion cross-section and the radius of gyration, the surface area and the volume of the convex hull suggest an optimal template for accurately mimicking the variations of Ω in a broad range of conformations, using only geometrical information and doing so at a very moderate computational cost. The existence of such a correlation is confirmed on the much larger protein α-lactalbumin.

Are long-range structural correlations behind the aggregration phenomena of polyglutamine diseases?

We have characterized the conformational ensembles of polyglutamine peptides of various lengths (ranging from to ), both with and without the presence of a C-terminal polyproline hexapeptide. For this, we used state-of-the-art molecular dynamics simulations combined with a novel statistical analysis to characterize the various properties of the backbone dihedral angles and secondary structural motifs of the glutamine residues. For (i.e., just above the pathological length for Huntington's disease), the equilibrium conformations of the monomer consist primarily of disordered, compact structures with non-negligible -helical and turn content. We also observed a relatively small population of extended structures suitable for forming aggregates including - and -strands, and - and -hairpins. Most importantly, for we find that there exists a long-range correlation (ranging for at least residues) among the backbone dihedral angles of the Q residues. For polyglutamine peptides below the pathological length, the population of the extended strands and hairpins is considerably smaller, and the correlations are short-range (at most residues apart). Adding a C-terminal hexaproline to suppresses both the population of these rare motifs and the long-range correlation of the dihedral angles. We argue that the long-range correlation of the polyglutamine homopeptide, along with the presence of these rare motifs, could be responsible for its aggregation phenomena.

Nine neurodegenerative diseases are caused by polyglutamine (polyQ) expansions greater than a given threshold in proteins with little or no homology except for the polyQ regions. The diseases all share a common feature: the formation of polyQ aggregates and eventual neuronal death. Using molecular dynamics simulations, we have explored the conformations of polyQ peptides. Results indicate that for peptides (i.e., just above the pathological length for Hungtington's disease), the equilibrium conformations were found to consist primarily of disordered, compact structures with a non-negligible -helical and turn content. We also observed a small population of extended structures suitable for forming aggregates. For peptides below the pathological length, the population of these structures was found to be considerably lower. For longer peptides, we found evidence for long-range correlations among the dihedral angles. This correlation turns out to be short-range for the smaller polyQ peptides, and is suppressed (along with the extended structural motifs) when a C-terminal polyproline tail is added to the peptides. We believe that the existence of these long-range correlations in above-threshold polyQ peptides, along with the presence of rare motifs, could be responsible for the experimentally observed aggregation phenomena associated with polyQ diseases.

Density of states-based molecular simulations

One of the central problems in statistical mechanics is that of finding the density of states of a system. Knowledge of the density of states of a system is equivalent to knowledge of its fundamental equation, from which all thermodynamic quantities can be obtained. Over the past several years molecular simulations have made considerable strides in their ability to determine the density of states of complex fluids and materials. In this review we discuss some of the more promising approaches proposed in the recent literature along with their advantages and limitations.

Dimerization of helical β-peptides in solution

Molecular simulations are used to examine the aggregation behavior of several β-peptides in explicit water. The particular peptides considered here adopt a helical, rodlike conformation in aqueous solution. Four distinct molecular sequences are considered. Earlier experimental studies have revealed the formation of ordered and disordered aggregates for such molecules, depending on sequence. The simulations reported here, which are conducted by resorting to metadynamics techniques, lead to free energy surfaces for dimerization of the peptides in water as a function of separation and relative orientation. Such surfaces are used to identify the molecular origins for the behaviors observed in the experiments.

Carbinolamine formation and dehydration in a DNA repair enzyme active site

In order to suggest detailed mechanistic hypotheses for the formation and dehydration of a key carbinolamine intermediate in the T4 pyrimidine dimer glycosylase (T4PDG) reaction, we have investigated these reactions using steered molecular dynamics with a coupled quantum mechanics–molecular mechanics potential (QM/MM). We carried out simulations of DNA abasic site carbinolamine formation with and without a water molecule restrained to remain within the active site quantum region. We recovered potentials of mean force (PMF) from thirty replicate reaction trajectories using Jarzynski averaging. We demonstrated feasible pathways involving water, as well as those independent of water participation. The water–independent enzyme–catalyzed reaction had a bias–corrected Jarzynski–average barrier height of approximately for the carbinolamine formation reaction and ) for the reverse reaction at this level of representation. When the proton transfer was facilitated with an intrinsic quantum water, the barrier height was approximately in the forward (formation) reaction and for the reverse. In addition, two modes of unsteered (free dynamics) carbinolamine dehydration were observed: in one, the quantum water participated as an intermediate proton transfer species, and in the other, the active site protonated glutamate hydrogen was directly transferred to the carbinolamine oxygen. Water–independent unforced proton transfer from the protonated active site glutamate carboxyl to the unprotonated N–terminal amine was also observed. In summary, complex proton transfer events, some involving water intermediates, were studied in QM/MM simulations of T4PDG bound to a DNA abasic site. Imine carbinolamine formation was characterized using steered QM/MM molecular dynamics. Dehydration of the carbinolamine intermediate to form the final imine product was observed in free, unsteered, QM/MM dynamics simulations, as was unforced acid-base transfer between the active site carboxylate and the N–terminal amine.

Rationale for the opposite stereochemistry of the major monoadducts and interstrand crosslinks formed by mitomycin C and its decarbamoylated analogue at CpG steps in DNA and the effect of cytosine modification on reactivity

Mitomycin C (MMC) is a potent antitumour agent that forms a covalent bond with the 2-amino group of selected guanines in the minor groove of double-stranded DNA following intracellular reduction of its quinone ring and opening of its aziridine moiety. At some 5'-CG-3' (CpG) steps the resulting monofunctional adduct can evolve towards a more deleterious bifunctional lesion, which is known as an interstrand crosslink (ICL). MMC reactivity is enhanced when the cytosine bases are methylated (5 MC) and decreased when they are replaced with 5-F-cytosine (5FC) whereas the stereochemical preference of alkylation changes upon decarbamoylation. We have studied three duplex oligonucleotides of general formula d(CGATAAXGCTAACG) in which X stands for C, 5MC or 5FC. Using a combination of molecular dynamics simulations in aqueous solution, quantum mechanics and continuum electrostatics, we have been able to (i) obtain a large series of snapshots that facilitate an understanding in atomic detail of the distinct stereochemistry of monoadduct and ICL formation by MMC and its decarbamoylated analogue, (ii) provide an explanation for the altered reactivity of MMC towards DNA molecules containing 5MC or 5FC, and (iii) show the distinct accommodation in the DNA minor groove of the different covalent modifications, particularly the most cytotoxic C1α and C1β ICLs.

The Role of Gln61 in HRas GTP hydrolysis: a quantum mechanics/molecular mechanics study

Activation of the water molecule involved in GTP hydrolysis within the HRas·RasGAP system is analyzed using a tailored approach based on hybrid quantum mechanics/molecular mechanics (QM/MM) simulation. A new path emerges: transfer of a proton from the attacking water molecule to a second water molecule, then a different proton is transferred from this second water molecule to the GTP. Gln(61) will stabilize the transient OH(-) and H(3)O(+) molecules thus generated. This newly proposed mechanism was generated by using, for the first time to our knowledge, the entire HRas-RasGAP protein complex in a QM/MM simulation context. It also offers a rational explanation for previous experimental results regarding the decrease of GTPase rate found in the HRas Q61A mutant and the increase exhibited by the HRas Q61E mutant.

Sucrose in aqueous solution revisited, Part 2: adaptively biased molecular dynamics simulations and computational analysis of NMR relaxation

We report 100 ns molecular dynamics simulations, at various temperatures, of sucrose in water (with concentrations of sucrose ranging from 0.02 to 4M), and in a 7:3 water-DMSO mixture. Convergence of the resulting conformational ensembles was checked using adaptive-biased simulations along the glycosidic Φ and ψ torsion angles. NMR relaxation parameters, including longitudinal (R₁) and transverse (R₂) relaxation rates, nuclear Overhauser enhancements (NOE), and generalized order parameter (S²) were computed from the resulting time-correlation functions. The amplitude and time scales of molecular motions change with temperature and concentration in ways that track closely with experimental results, and are consistent with a model in which sucrose conformational fluctuations are limited (with 80-90% of the conformations having ϕ-ψ values within 20° of an average conformation), but with some important differences in conformation between pure water and DMSO-water mixtures.

PPII propensity of multiple-guest amino acids in a proline-rich environment

There has been considerable debate about the intrinsic PPII propensity of amino acid residues in denatured polypeptides. Experimentally, this scale is based on the behavior of guest amino acid residues placed in the middle of proline-based hosts. We have used classical molecular dynamics simulations combined with replica-exchange methods to carry out a comprehensive analysis of the conformational equilibria of proline-based host oligopeptides with multiple guest amino acids including alanine, glutamine, valine, and asparagine. The tracked structural characteristics include the secondary structural motifs based on the Ramachandran angles and the cis/trans isomerization of the prolyl bonds. In agreement with our recent study of single amino acid guests, we did not observe an intrinsic PPII propensity in any of the guest amino acids in a multiple-guest setting. Instead, the experimental results can be explained in terms of (i) the steric restrictions imposed on the C-terminal guest amino acid that is immediately followed by a proline residue and (ii) an increase in the trans content of the prolyl bonds due to the presence of guest residues. In terms of the latter, we found that the more guests added to the system, the larger the increase in the trans content of the prolyl bonds, which results in an effective increase in the PPII content of the peptide.

Optimized Potential of Mean Force Calculations for Standard Binding Free Energies

The prediction of protein-ligand binding free energies is an important goal of computational biochemistry, yet accuracy, reproducibility, and cost remain a problem. Nevertheless, these are essential requirements for computational methods to become standard binding prediction tools in discovery pipelines. Here, we present the results of an extensive search for an optimal method based on an ensemble of umbrella sampling all-atom molecular simulations tested on the phosphorylated tetrapeptide, pYEEI, binding to the SH2 domain, resulting in an accurate and converged binding free energy of -9.0 ± 0.5 kcal/mol (compared to an experimental value of -8.0 ± 0.1 kcal/mol). We find that a minimum of 300 ns of sampling is required for every prediction, a target easily achievable using new generation accelerated MD codes. Convergence is obtained by using an ensemble of simulations per window, each starting from different initial conformations, and by optimizing window-width, orthogonal restraints, reaction coordinate harmonic potentials, and window-sample time. The use of uncorrelated initial conformations in neighboring windows is important for correctly sampling conformational transitions from the unbound to bound states that affect significantly the precision of the calculations. This methodology thus provides a general recipe for reproducible and practical computations of binding free energies for a class of semirigid protein-ligand systems, within the limit of the accuracy of the force field used.

A statistical analysis of the PPII propensity of amino acid guests in proline-rich peptides

There has been considerable debate about the intrinsic PPII propensity of amino-acid residues in denatured polypeptides. Experimentally, the propensity scale is based on the behavior of guest amino-acid residues placed in the middle of polyproline hosts. We have used classical molecular dynamics simulations, with state-of-the-art force fields to carry out a comprehensive analysis of the conformational equilibria of the proline-based host oligopeptides with single guests. The tracked structural characteristics include the PPII content, the cis/trans isomerization of the prolyl bonds, the puckering of the pyrrolidine rings of the proline residues, and the secondary structural motifs. We find no evidence for an intrinsic PPII propensity in any of the guest amino acids other than proline. Instead, the PPII content as derived from experiments may be explained in terms of: 1), a local correlation between the dihedral angles of the guest amino acid and the proline residue immediately preceding it; and 2), a nonlocal correlation between the cis/trans states of the peptide bonds. In terms of the latter, we find that the presence of a guest (other than proline, tyrosine, or tryptophan) increases the trans content of most of the prolyl bonds, which results in an effective increase of the peptide PPII content. With respect to the local dihedral correlations, we find that these are well described in terms of the so-called odds-ratio statistic. Expressed in terms of free energy language, the PPII content based on the odds-ratio of the relevant residues correlate well with the experimentally measured PPII content.

Free-energy landscapes from adaptively biased methods: application to quantum systems

Several parallel adaptive biasing methods are applied to the calculation of free-energy pathways along reaction coordinates, choosing as a difficult example the double-funnel landscape of the 38-atom Lennard-Jones cluster. In the case of classical statistics, the Wang-Landau and adaptively biased molecular-dynamics (ABMD) methods are both found efficient if multiple walkers and replication and deletion schemes are used. An extension of the ABMD technique to quantum systems, implemented through the path-integral MD framework, is presented and tested on Ne38 against the quantum superposition method.

Balancing simulation accuracy and efficiency with the Amber united atom force field

We have analyzed the quality of a recently proposed Amber united-atom model and its overall efficiency in ab initio folding and thermodynamic sampling of two stable beta-hairpins. It is found that the mean backbone structures are quite consistent between the simulations in the united-atom and its corresponding all-atom models in Amber. More importantly, the simulated beta turns are also consistent between the two models. Finally, the chemical shifts on H alpha are highly consistent between simulations in the two models, although the simulated chemical shifts are lower than experiment, indicating less structured peptides, probably due to the omission of the hydrophobic term in the simulations. More interestingly, the stabilities of both beta-hairpins at room temperature are similar to those derived from the NMR measurement, whether the united-atom or the all-atom model is used. Detailed analysis shows high percentages of backbone torsion angles within the beta region and high percentages of native contacts. Given the reasonable quality of the united-atom model with respect to experimental data, we have further studied the simulation efficiency of the united-atom model over the all-atom model. Our data shows that the united-atom model is a factor of 6-8 faster than the all-atom model as measured with the ab initio first pass folding time for the two tested beta-hairpins. Detailed structural analysis shows that all ab initio folded trajectories enter the native basin, whether the united-atom model or the all-atom model is used. Finally, we have also studied the simulation efficiency of the united-atom model as measured in terms of how fast thermodynamic convergence can be achieved. It is apparent that the united-atom simulations reach convergence faster than the all-atom simulations with respect to both mean potential energies and mean native contacts. These findings show that the efficiency of the united-atom model is clearly beyond the per-step dynamics simulation of about 2 over the all-atom model. Thus, reasonable reduction of a protein model can be achieved with improved sampling efficiency while still preserving a high level of accuracy for applications in both ab initio folding and thermodynamic sampling. This study motivates us to develop more simplified protein models with sufficient consistency with the all-atom models for enhanced conformational sampling.

Conformations and free energy landscapes of polyproline peptides

The structure of the proline amino acid allows folded polyproline peptides to exist as both left- (PPII) and right-handed (PPI) helices. We have characterized the free energy landscapes of hexamer, nanomer, and tridecamer polyproline peptides in gas phase and implicit water as well as explicit hexane and 1-propanol for the nanomer. To enhance the sampling provided by regular molecular dynamics, we used the recently developed adaptively biased molecular dynamics method, which describes Landau free energy maps in terms of relevant collective variables. These maps, as a function of the collective variables of handedness, radius of gyration, and three others based on the peptide torsion angle omega, were used to determine the relative stability of the different structures, along with an estimate of the transition pathways connecting the different minima. Results show the existence of several metastable isomers and therefore provide a complementary view to experimental conclusions based on photo-induced electron transfer experiments with regard to the existence of stable heterogeneous subpopulations in PPII polyproline.

Enhanced conformational sampling of nucleic acids by a new Hamiltonian replica exchange molecular dynamics approach

Although molecular dynamics (MD) simulations have been applied frequently to study flexible molecules, the sampling of conformational states separated by barriers is limited due to currently possible simulation time scales. Replica-exchange (Rex)MD simulations that allow for exchanges between simulations performed at different temperatures (T-RexMD) can achieve improved conformational sampling. However, in the case of T-RexMD the computational demand grows rapidly with system size. A Hamiltonian RexMD method that specifically enhances coupled dihedral angle transitions has been developed. The method employs added biasing potentials as replica parameters that destabilize available dihedral substates and was applied to study coupled dihedral transitions in nucleic acid molecules. The biasing potentials can be either fixed at the beginning of the simulation or optimized during an equilibration phase. The method was extensively tested and compared to conventional MD simulations and T-RexMD simulations on an adenine dinucleotide system and on a DNA abasic site. The biasing potential RexMD method showed improved sampling of conformational substates compared to conventional MD simulations similar to T-RexMD simulations but at a fraction of the computational demand. It is well suited to study systematically the fine structure and dynamics of large nucleic acids under realistic conditions including explicit solvent and ions and can be easily extended to other types of molecules.

Intrinsic Free Energy of the Conformational Transition of the KcsA Signature Peptide from Conducting to Nonconducting State

We explore a conformational transition of the TATTVGYG signature peptide of the KcsA ion selectivity filter and its GYG to AYA mutant from the conducting α-strand state into the nonconducting pII-like state using a novel technique for multidimensional optimization of transition path ensembles and free energy calculations. We find that the wild type peptide, unlike the mutant, intrinsically favors the conducting state due to G77 backbone propensities and additional hydrophobic interaction between the V76 and Y78 side chains in water. The molecular mechanical free energy profiles in explicit water are in very good agreement with the corresponding adiabatic energies from the Generalized Born Molecular Volume (GBMV) implicit solvent model. However comparisons of the energies to higher level B3LYP/6-31G(d) Density Functional Theory calculations with Polarizable Continuum Model (PCM) suggest that the nonconducting state might be more favorable than predicted by molecular mechanics simulations. By extrapolating the single peptide results to the tetrameric channel, we propose a novel hypothesis for the ion selectivity mechanism.