Curious structure in "canonical" alanine-based peptides

Configurational Entropy of Folded Proteins and Its Importance for Intrinsically Disordered Proteins

Many pairwise additive force fields are in active use for intrinsically disordered proteins (IDPs) and regions (IDRs), some of which modify energetic terms to improve the description of IDPs/IDRs but are largely in disagreement with solution experiments for the disordered states. This work considers a new direction-the connection to configurational entropy-and how it might change the nature of our understanding of protein force field development to equally well encompass globular proteins, IDRs/IDPs, and disorder-to-order transitions. We have evaluated representative pairwise and many-body protein and water force fields against experimental data on representative IDPs and IDRs, a peptide that undergoes a disorder-to-order transition, for seven globular proteins ranging in size from 130 to 266 amino acids. We find that force fields with the largest statistical fluctuations consistent with the radius of gyration and universal Lindemann values for folded states simultaneously better describe IDPs and IDRs and disorder-to-order transitions. Hence, the crux of what a force field should exhibit to well describe IDRs/IDPs is not just the balance between protein and water energetics but the balance between energetic effects and configurational entropy of folded states of globular proteins.

Molecular dynamics simulations identify the regions of compromised thermostability in SazCA

The present study examined the structure and dynamics of the most active and thermostable carbonic anhydrase, SazCA, probed using molecular dynamics simulations. The molecular system was described by widely used biological force-fields (AMBER, CHARMM22, CHARMM36, and OPLS-AA) in conjunction with TIP3P water model. The comparison of molecular dynamics simulation results suggested AMBER to be a suitable choice to describe the structure and dynamics of SazCA. In addition to this, we also addressed the effect of temperature on the stability of SazCA. We performed molecular dynamics simulations at 313, 333, 353, 373, and 393 K to study the relationship between thermostability and flexibility in SazCA. The amino acid residues VAL98, ASN99, GLY100, LYS101, GLU145, and HIS207 were identified as the most flexible residues from root-mean-square fluctuations. The salt bridge analysis showed that ion-pairs ASP113-LYS81, ASP115-LYS81, ASP115-LYS114, GLU144-LYS143, and GLU144-LYS206, were responsible for the compromised thermal stability of SazCA.

FF12MC: A revised AMBER forcefield and new protein simulation protocol

Specialized to simulate proteins in molecular dynamics (MD) simulations with explicit solvation, FF12MC is a combination of a new protein simulation protocol employing uniformly reduced atomic masses by tenfold and a revised AMBER forcefield FF99 with (i) shortened CH bonds, (ii) removal of torsions involving a nonperipheral sp(3) atom, and (iii) reduced 1-4 interaction scaling factors of torsions ϕ and ψ. This article reports that in multiple, distinct, independent, unrestricted, unbiased, isobaric-isothermal, and classical MD simulations FF12MC can (i) simulate the experimentally observed flipping between left- and right-handed configurations for C14-C38 of BPTI in solution, (ii) autonomously fold chignolin, CLN025, and Trp-cage with folding times that agree with the experimental values, (iii) simulate subsequent unfolding and refolding of these miniproteins, and (iv) achieve a robust Z score of 1.33 for refining protein models TMR01, TMR04, and TMR07. By comparison, the latest general-purpose AMBER forcefield FF14SB locks the C14-C38 bond to the right-handed configuration in solution under the same protein simulation conditions. Statistical survival analysis shows that FF12MC folds chignolin and CLN025 in isobaric-isothermal MD simulations 2-4 times faster than FF14SB under the same protein simulation conditions. These results suggest that FF12MC may be used for protein simulations to study kinetics and thermodynamics of miniprotein folding as well as protein structure and dynamics. Proteins 2016; 84:1490-1516. © 2016 The Authors Proteins: Structure, Function, and Bioinformatics Published by Wiley Periodicals, Inc.

Conformational Interconversions of Amino Acid Derivatives

Exhaustive conformational interconversions including transition structure analyses of N-acetyl-l-glycine-N-methylamide as well as its alanine, serine, and cysteine analogues have been investigated at the MP2/6-31G** level, yielding a total of 142 transition states. Improved estimates of relative energies were obtained by separately extrapolating the Hartree-Fock and MP2 energies to the basis set limit and adding the difference between CCSD(T) and MP2 results with the cc-pVDZ basis set to the extrapolated MP2 results. The performance of eight empirical force fields (AMBER94, AMBER14SB, MM2, MM3, MMFFs, CHARMM22_CMAP, OPLS_2005, and AMOEBAPRO13) in reproducing ab initio energies of transition states was tested. Our results indicate that commonly used class I force fields employing a fixed partial charge model for the electrostatic interaction provide mean errors in the ∼10 kJ/mol range for energies of conformational transition states for amino acid conformers. Modern reparametrized versions, such as CHARMM22_CMAP, and polarizable force fields, such as AMOEBAPRO13, have slightly lower mean errors, but maximal errors are still in the 35 kJ/mol range. There are differences between the force fields in their ability for reproducing conformational transitions classified according to backbone/side-chain or regions in the Ramachandran angles, but the data set is likely too small to draw any general conclusions. Errors in conformational interconversion barriers by ∼10 kJ/mol suggest that the commonly used force field may bias certain types of transitions by several orders of magnitude in rate and thus lead to incorrect dynamics in simulations. It is therefore suggested that information for conformational transition states should be included in parametrizations of new force fields.

Protein hydrogen exchange: testing current models

To investigate the determinants of protein hydrogen exchange (HX), HX rates of most of the backbone amide hydrogens of Staphylococcal nuclease were measured by NMR methods. A modified analysis was used to improve accuracy for the faster hydrogens. HX rates of both near surface and well buried hydrogens are spread over more than 7 orders of magnitude. These results were compared with previous hypotheses for HX rate determination. Contrary to a common assumption, proximity to the surface of the native protein does not usually produce fast exchange. The slow HX rates for unprotected surface hydrogens are not well explained by local electrostatic field. The ability of buried hydrogens to exchange is not explained by a solvent penetration mechanism. The exchange rates of structurally protected hydrogens are not well predicted by algorithms that depend only on local interactions or only on transient unfolding reactions. These observations identify some of the present difficulties of HX rate prediction and suggest the need for returning to a detailed hydrogen by hydrogen analysis to examine the bases of structure-rate relationships, as described in the companion paper (Skinner et al., Protein Sci 2012;21:996-1005).

Residue-specific side-chain packing determines the backbone dynamics of transmembrane model helices

The transmembrane domains (TMDs) of membrane-fusogenic proteins contain an overabundance of β-branched residues. In a previous effort to systematically study the relation among valine content, fusogenicity, and helix dynamics, we developed model TMDs that we termed LV-peptides. The content and position of valine in LV-peptides determine their fusogenicity and backbone dynamics, as shown experimentally. Here, we analyze their conformational dynamics and the underlying molecular forces using molecular-dynamics simulations. Our study reveals that backbone dynamics is correlated with the efficiency of side-chain to side-chain van der Waals packing between consecutive turns of the helix. Leu side chains rapidly interconvert between two rotameric states, thus favoring contacts to its i±3 and i±4 neighbors. Stereochemical restraints acting on valine side chains in the α-helix force both β-substituents into an orientation where i,i±3 interactions are less favorable than i,i±4 interactions, thus inducing a local packing deficiency at VV3 motifs. We provide a quantitative molecular model to explain the relationship among chain connectivity, side-chain mobility, and backbone flexibility. We expect that this mechanism also defines the backbone flexibility of natural TMDs.

Salt dependence of an alpha-helical peptide folding energy landscapes

We used CD, UV resonance Raman spectroscopy, and molecular dynamics simulation to examine the impact of salts on the conformational equilibria and the Ramachandran Psi angle (un)folding Gibbs free energy landscape coordinate of a mainly polyalanine alpha-helical peptide, AP of sequence AAAAA(AAARA)(3)A. NaClO(4) stabilizes alpha-helical-like conformations more than does NaCl, which stabilizes more than Na(2)SO(4) at identical ionic strengths. This alpha-helix stabilization ordering is the reverse of the Hofmeister series of anions in their ability to disorder water hydrogen bonding. Much of the NaClO(4) alpha-helix stabilization results from ClO(4)(-) association with the AP terminal -NH(3)(+) groups and Arg side chains. ClO(4)(-) stabilizes 3(10)-helix conformations but destabilizes turn conformations. The decreased Cl(-) and SO(4)(2-) AP alpha-helix stabilization probably results from a decreased association with the Arg and terminal -NH(3)(+) groups. Cl(-) is expected to have a smaller binding affinity and thus stabilizes alpha-helical conformations intermediately between NaClO(4) and Na(2)SO(4). Electrostatic screening stabilizes pi-bulge conformations.

TIGER2: an improved algorithm for temperature intervals with global exchange of replicas

An empirical sampling method for molecular simulation based on "temperature intervals with global exchange of replicas" (TIGER2) has been developed to reduce the high demand for computational resources and the low computational efficiency of the conventional replica-exchange molecular dynamics (REMD) method. This new method overcomes the limitation of its previous version, called TIGER, which requires the assumption of constant heat capacity during quenching of replicas from elevated temperatures to the baseline temperature. The robustness of the TIGER2 method is examined by comparing it against a Metropolis Monte Carlo simulation for sampling the conformational distribution of a single butane molecule in vacuum, a REMD simulation for sampling the behavior of alanine dipeptide in explicit solvent, and REMD simulations for sampling the folding behavior of two peptides, (AAQAA)(3) and chignolin, in implicit solvent. The agreement between the results from these conventional sampling methods and the TIGER2 simulations indicates that the TIGER2 algorithm is able to closely approximate a Boltzmann-weighted ensemble of states for these systems but without the limiting assumptions that were required for the original TIGER algorithm. TIGER2 is an efficient replica-exchange sampling method that enables the number of replicas that are used for a replica-exchange simulation to be substantially reduced compared to the conventional REMD method.

Trapping a folding intermediate of the alpha-helix: stabilization of the pi-helix

We report the design, synthesis, and characterization of a short peptide trapped in a pi-helix configuration. This high-energy conformation was nucleated by a preorganized pi-turn, which was obtained by replacing an N-terminal intramolecular main chain i and i + 5 hydrogen bond with a carbon-carbon bond. Our studies highlight the nucleation parameter as a key factor contributing to the relative instability of the pi-helix and allow us to estimate fundamental helix-coil transition parameters for this conformation.

Comparison of protein force fields for molecular dynamics simulations

In the context of molecular dynamics simulations of proteins, the term "force field" refers to the combination of a mathematical formula and associated parameters that are used to describe the energy of the protein as a function of its atomic coordinates. In this review, we describe the functional forms and parameterization protocols of the widely used biomolecular force fields Amber, CHARMM, GROMOS, and OPLS-AA. We also summarize the ability of various readily available noncommercial molecular dynamics packages to perform simulations using these force fields, as well as to use modern methods for the generation of constant-temperature, constant-pressure ensembles and to treat long-range interactions. Finally, we finish with a discussion of the ability of these force fields to support the modeling of proteins in conjunction with nucleic acids, lipids, carbohydrates, and/or small molecules.

Molecular dynamics simulations of the catalytic pathway of a cysteine protease: a combined QM/MM study of human cathepsin K

Molecular dynamics simulations using a combined QM/MM potential have been performed to study the catalytic mechanism of human cathepsin K, a member of the papain family of cysteine proteases. We have determined the two-dimensional free energy surfaces of both acylation and deacylation steps to characterize the reaction mechanism. These free energy profiles show that the acylation step is rate limiting with a barrier height of 19.8 kcal/mol in human cathepsin K and of 29.3 kcal/mol in aqueous solution. The free energy of activation for the deacylation step is 16.7 kcal/mol in cathepsin K and 17.8 kcal/mol in aqueous solution. The reduction of free energy barrier is achieved by stabilization of the oxyanion in the transition state. Interestingly, although the "oxyanion hole" has been formed in the Michaelis complex, the amide units do not donate hydrogen bonds directly to the carbonyl oxygen of the substrate, but they stabilize the thiolate anion nucleophile. Hydrogen-bonding interactions are induced as the substrate amide group approaches the nucleophile, moving more than 2 A and placing the oxyanion in contact with Gln19 and the backbone amide of Cys25. The hydrolysis of peptide substrate shares a common mechanism both for the catalyzed reaction in human cathepsin K and for the uncatalyzed reaction in water. Overall, the nucleophilic attack by Cys25 thiolate and the proton-transfer reaction from His162 to the amide nitrogen are highly coupled, whereas a tetrahedral intermediate is formed along the nucleophilic reaction pathway.

Monte carlo folding of trans-membrane helical peptides in an implicit generalized Born membrane

An efficient Monte Carlo (MC) algorithm using concerted backbone rotations is combined with a recently developed implicit membrane model to simulate the folding of the hydrophobic transmembrane domain M2TM of the M2 protein from influenza A virus and Sarcolipin at atomic resolution. The implicit membrane environment is based on generalized Born theory and has been calibrated against experimental data. The MC sampling has previously been used to fold several small polypeptides and been shown to be equivalent to molecular dynamics (MD). In combination with a replica exchange algorithm, M2TM is found to form continuous membrane spanning helical conformations for low temperature replicas. Sarcolipin is only partially helical, in agreement with the experimental NMR structures in lipid bilayers and detergent micelles. Higher temperature replicas exhibit a rapidly decreasing helicity, in agreement with expected thermodynamic behavior. To exclude the possibility of an erroneous helical bias in the simulations, the model is tested by sampling a synthetic Alanine-rich polypeptide of known helicity. The results demonstrate there is no overstabilization of helical conformations, indicating that the implicit model captures the essential components of the native membrane environment for M2TM and Sarcolipin.

Direct UV Raman monitoring of 3(10)-helix and pi-bulge premelting during alpha-helix unfolding

We used UV resonance Raman (UVRR) spectroscopy exciting at approximately 200 nm within the peptide bond pi --> pi* transitions to selectively study the amide vibrations of peptide bonds during alpha-helix melting. The dependence of the amide frequencies on their Psi Ramachandran angles and hydrogen bonding enables us, for the first time, to experimentally determine the temperature dependence of the peptide bond Psi Ramachandran angle population distribution of a 21-residue mainly alanine peptide. These Psi distributions allow us to easily discriminate between alpha-helix, 3(10)-helix and pi-helix/bulge conformations, obtain their individual melting curves, and estimate the corresponding Zimm and Bragg parameters. A striking finding is that alpha-helix melting is more cooperative and shows a higher melting temperature than previously erroneously observed. These Psi distributions also enable the experimental determination of the Gibbs free energy landscape along the Psi reaction coordinate, which further allows us to estimate the free energy barriers along the AP melting pathway. These results will serve as a benchmark for the numerous untested theoretical studies of protein and peptide folding.

UV resonance Raman measurements of poly-L-lysine's conformational energy landscapes: dependence on perchlorate concentration and temperature

UV resonance Raman spectroscopy has been used to determine the conformational energy landscape of poly-L-lysine (PLL) in the presence of NaClO4 as a function of temperature. At 1 degree C, in the presence of 0.83 M NaClO4, PLL shows an approximately 86% alpha-helix-like content, which contains alpha-helix and pi-bulge/helix conformations. The high alpha-helix-like content of PLL occurs because of charge screening due to strong ion-pair formation between ClO4- and the lysine side chain -NH3+. As the temperature increases from 1 to 60 degrees C, the alpha-helix and pi-bulge/helix conformations melt into extended conformations (PPII and 2.51-helix). We calculate the Psi Ramachandran angle distribution of the PLL peptide bonds from the UV Raman spectra which allows us to calculate the PLL (un)folding energy landscapes along the Psi reaction coordinate. We observe a basin in the Psi angle conformational space associated with alpha-helix and pi-bulge/helix conformations and another basin for the extended PPII and 2.51-helical conformations.

UV Raman spatially resolved melting dynamics of isotopically labeled polyalanyl peptide: slow alpha-helix melting follows 3(10)-helices and pi-bulges premelting

We used UV resonance Raman (UVRR) to examine the spatial dependence of the T-jump secondary structure relaxation of an isotopically labeled 21-residue mainly Ala peptide, AdP. The AdP penultimate Ala residues were perdeuterated, leaving the central residues hydrogenated, to allow separate monitoring of melting of the middle versus the end peptide bonds. For 5 to 30 degrees C T-jumps, the central peptide bonds show a approximately 2-fold slower relaxation time (189 +/- 31 ns) than do the exterior peptide bonds (97 +/- 15 ns). In contrast, for a 20 to 40 degrees C T-jump, the central peptide bond relaxation appears to be faster (56 +/- 6 ns) than that of the penultimate peptide bonds (131 +/- 46 ns). We show that, if the data are modeled as a two-state transition, we find that only exterior peptide bonds show anti-Arrhenius folding behavior; the middle peptide bonds show both normal Arrhenius-like folding and unfolding. This anti-Arrhenius behavior results from the involvement of pi-bulges/helices and 3(10)-helix states in the melting. The unusual temperature dependence of the (un)folding rates of the interior and exterior peptide bonds is due to the different relative (un)folding rates of 3(10)-helices, alpha-helices, and pi-bulges/helices. Pure alpha-helix unfolding rates are approximately 12-fold slower (approximately 1 micros) than that of pi-bulges and 3(10)-helices. In addition, we also find that the alpha-helix is most stable at the AdP N-terminus where eight consecutive Ala occur, whereas the three hydrophilic Arg located in the middle and at the C-terminus destabilize the alpha-helix in these regions and induce defects such as pi-bulges and 3(10)-helices.

Helices and Sheets in vacuo

The structures and properties of unsolvated peptides large enough to possess secondary structure have been examined by experiments and simulations. Some of the factors that stabilize unsolvated helices and sheets have been identified. The charge, in particular, plays a critical role in stabilizing alpha-helices and destabilizing beta-sheets. Some helices are much more stable in vacuum than in aqueous solution. Factors like helix propensity, context, and the incorporation of specific stabilizing interactions have been examined. The helix propensities in vacuum differ from those found in solution. Studies of the hydration of unsolvated peptides can be performed one water molecule at a time. The first few water molecules only bind weakly to unsolvated peptides, and they bind much more strongly to some conformations than to others. The most favorable binding locations are not the protonation sites, but clefts or pockets where a water molecule can establish a network of hydrogen bonds. Non-covalent interactions between secondary structure elements leads to the formation of tertiary structure. Helical peptides assemble into complexes with a variety of intriguing structures. The intramolecular coupling of helices to make antiparallel coiled-coil geometries has also been investigated with model peptides.

How large is an alpha-helix? Studies of the radii of gyration of helical peptides by small-angle X-ray scattering and molecular dynamics

Using synchrotron radiation and the small-angle X-ray scattering technique we have measured the radii of gyration of a series of alanine-based alpha-helix-forming peptides of the composition Ace-(AAKAA)(n)-GY-NH(2), n=2-7, in aqueous solvent at 10(+/-1) degrees C. In contrast to other techniques typically used to study alpha-helices in isolation (such as nuclear magnetic resonance and circular dichroism), small-angle X-ray scattering reports on the global structure of a molecule and, as such, provides complementary information to these other, more sequence-local measuring techniques. The radii of gyration that we measure are, except for the 12-mer, lower than the radii of gyration of ideal alpha-helices or helices with frayed ends of the equivalent sequence-length. For example, the measured radius of gyration of the 37-mer is 14.2(+/-0.6)A, which is to be compared with the radius of gyration of an ideal 37-mer alpha-helix of 17.6A. Attempts are made to analyze the origin of this discrepancy in terms of the analytical Zimm-Bragg-Nagai (ZBN) theory, as well as distributed computing explicit solvent molecular dynamics simulations using two variants of the AMBER force-field. The ZBN theory, which treats helices as cylinders connected by random walk segments, predicts markedly larger radii of gyration than those measured. This is true even when the persistence length of the random walk parts is taken to be extremely short (about one residue). Similarly, the molecular dynamics simulations, at the level of sampling available to us, give inaccurate values of the radii of gyration of the molecules (by overestimating them by around 25% for longer peptides) and/or their helical content. We conclude that even at the short sequences examined here (< or =37 amino acid residues), these alpha-helical peptides behave as fluctuating semi-broken rods rather than straight cylinders with frayed ends.

Monte Carlo studies of folding, dynamics, and stability in alpha-helices

Folding simulations of polyalanine peptides were carried out using an off-lattice Monte Carlo simulation technique. The peptide was represented as a chain of residues, each of which contains two interaction sites: one corresponding to the C(alpha) atom and the other to the side chain. A statistical potential was used to describe the interaction between these sites. The preferred conformations of the peptide chain on the energy surface, starting from several initial conditions, were searched by perturbations on its generalized coordinates with the Metropolis criterion. We observed that, at low temperatures, the effective energy was low and the helix content high. The calculated helix propagation (s) and nucleation (sigma) parameters of the Zimm-Bragg model were in reasonable agreement with the empirical data. Exploration of the energy surface of the alanine-based peptides (AAQAA)(3) and AAAAA(AAARA)(3)A demonstrated that their behavior is similar to that of polyalanine, in regard to their effective energy, helix content, and the temperature-dependence of their helicity. In contrast, stable secondary structures were not observed for (Gly)(20) at similar temperatures, which is consistent with the nonfolder nature of this peptide. The fluctuations in the slowest dynamics mode, which describe the elastic behavior of the chain, showed that as the temperature decreases, the polyalanine peptides become stiffer and retain conformations with higher helix content. Clustering of conformations during the folding phase implied that polyalanine folds into a helix through fewer numbers of intermediate conformations as the temperature decreases.

A structural model of polyglutamine determined from a host-guest method combining experiments and landscape theory

Modeling the structure of natively disordered peptides has proved difficult due to the lack of structural information on these peptides. In this work, we use a novel application of the host-guest method, combining folding theory with experiments, to model the structure of natively disordered polyglutamine peptides. Initially, a minimalist molecular model (C(alpha)C(beta)) of CI2 is developed with a structurally based potential and captures many of the folding properties of CI2 determined from experiments. Next, polyglutamine "guest" inserts of increasing length are introduced into the CI2 "host" model and the polyglutamine is modeled to match the resultant change in CI2 thermodynamic stability between simulations and experiments. The polyglutamine model that best mimics the experimental changes in CI2 thermodynamic stability has 1), a beta-strand dihedral preference and 2), an attractive energy between polyglutamine atoms 0.75-times the attractive energy between the CI2 host Go-contacts. When free-energy differences in the CI2 host-guest system are correctly modeled at varying lengths of polyglutamine guest inserts, the kinetic folding rates and structural perturbation of these CI2 insert mutants are also correctly captured in simulations without any additional parameter adjustment. In agreement with experiments, the residues showing structural perturbation are located in the immediate vicinity of the loop insert. The simulated polyglutamine loop insert predominantly adopts extended random coil conformations, a structural model consistent with low resolution experimental methods. The agreement between simulation and experimental CI2 folding rates, CI2 structural perturbation, and polyglutamine insert structure show that this host-guest method can select a physically realistic model for inserted polyglutamine. If other amyloid peptides can be inserted into stable protein hosts and the stabilities of these host-guest mutants determined, this novel host-guest method may prove useful to determine structural preferences of these intractable but biologically relevant protein fragments.

Structural and functional characterization of pi bulges and other short intrahelical deformations

We data-mined the Protein Data Bank for short intrahelical deformations, including pi bulges. These are defined as a contiguous stretch of intrahelical residues deviating from the standard alpha-helical i-->i-4 hydrogen bonding pattern, bilaterally flanked by at least one alpha-helical turn resulting in a helix kink of less than 40 degrees. We find that such motifs exist in 4.7% of a PDB subset filtered by quality metrics (resolution <2.5 A, R-factor <0.25, sequence identity <35%). These are typically characterized by at least one i-->i-5 main chain hydrogen bond, with energetically favorable main chain dihedral angles, followed by a variable number of main chain carbonyl groups that do not accept intrahelical main chain hydrogen bonds. Their stabilization commonly occurs via hydrogen bonding to water molecules or polar groups. Numerous deformations are implicated in basic yet vital functional roles, commonly as ligand binding site contributors.

Discriminating the Helical Forms of Peptides by NMR and Molecular Dynamics Simulation

The HNCO NMR pulse sequence was applied to three selectively labeled (15)N and (13)C isotopic homologues of the peptide Ac-WAAAH(AAARA)(3)A-NH(2) to probe directly for hydrogen bonds between residues 8 and 11 (characteristic of a 3(10)-helix), 8 and 12 (alpha-helix), and 8 and 13 (pi-helix). The experiments demonstrate conclusively, and in agreement with circular dichroism studies, that the center of the peptide is alpha-helical; there is no discernible 3(10)- or pi-helix at these specific positions. Molecular dynamics simulations of the preceding peptide and Ac-(AAAAK)(3)A-NH(2) in water using the potential energy parameter set CHARMM22/CMAP correctly yield an alpha-helix, in contrast to simulations with the set CHARMM22, which result in a pi-helix.

Extending the treatment of backbone energetics in protein force fields: limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations

Computational studies of proteins based on empirical force fields represent a powerful tool to obtain structure-function relationships at an atomic level, and are central in current efforts to solve the protein folding problem. The results from studies applying these tools are, however, dependent on the quality of the force fields used. In particular, accurate treatment of the peptide backbone is crucial to achieve representative conformational distributions in simulation studies. To improve the treatment of the peptide backbone, quantum mechanical (QM) and molecular mechanical (MM) calculations were undertaken on the alanine, glycine, and proline dipeptides, and the results from these calculations were combined with molecular dynamics (MD) simulations of proteins in crystal and aqueous environments. QM potential energy maps of the alanine and glycine dipeptides at the LMP2/cc-pVxZ//MP2/6-31G* levels, where x = D, T, and Q, were determined, and are compared to available QM studies on these molecules. The LMP2/cc-pVQZ//MP2/6-31G* energy surfaces for all three dipeptides were then used to improve the MM treatment of the dipeptides. These improvements included additional parameter optimization via Monte Carlo simulated annealing and extension of the potential energy function to contain peptide backbone phi, psi dihedral crossterms or a phi, psi grid-based energy correction term. Simultaneously, MD simulations of up to seven proteins in their crystalline environments were used to validate the force field enhancements. Comparison with QM and crystallographic data showed that an additional optimization of the phi, psi dihedral parameters along with the grid-based energy correction were required to yield significant improvements over the CHARMM22 force field. However, systematic deviations in the treatment of phi and psi in the helical and sheet regions were evident. Accordingly, empirical adjustments were made to the grid-based energy correction for alanine and glycine to account for these systematic differences. These adjustments lead to greater deviations from QM data for the two dipeptides but also yielded improved agreement with experimental crystallographic data. These improvements enhance the quality of the CHARMM force field in treating proteins. This extension of the potential energy function is anticipated to facilitate improved treatment of biological macromolecules via MM approaches in general.

π-Helix Preference in Unsolvated Peptides

Ion mobility measurements have been used to examine helix formations in the gas phase for a series of alanine/glycine-based peptides that incorporate a glutamic acid (E) and lysine (K) at various positions along the backbone. Incorporation of an EK pair lowers the percent helix for all positions (presumably because hydrogen bonding between the backbone and the E and K side chains stabilize the nonhelical globular conformations). The largest percent helix is found when the EK pair is in an i,i+5 arrangement, which suggests that the preferred helical conformation for these peptides is a pi-helix. This conclusion is supported by comparison of cross sections deduced from the ion-mobility measurements to average cross sections calculated for conformations obtained from molecular dynamics simulations. The glutamic acid and lysine may form an ion pair that is stabilized by interactions with the helix macro-dipole.

Reaction-path energetics and kinetics of the hydride transfer reaction catalyzed by dihydrofolate reductase

We have studied the hydride transfer reaction catalyzed by the enzyme dihydrofolate reductase (DHFR) and the coenzyme nicotinamide adenine dinucleotide phosphate (NADPH); the substrate is 5-protonated 7,8-dihydrofolate, and the product is tetrahydrofolate. The potential energy surface is modeled by a combined quantum mechanical-molecular mechanical (QM/MM) method employing Austin model 1 (AM1) and a simple valence bond potential for 69 QM atoms and employing the CHARMM22 and TIP3P molecular mechanics force fields for the other 21 399 atoms; the QM and MM regions are joined by two boundary atoms treated by the generalized hybrid orbital (GHO) method. All simulations are carried out using periodic boundary conditions at neutral pH and 298 K. In stage 1, a reaction coordinate is defined as the difference between the breaking and forming bond distances to the hydride ion, and a quasithermodynamic free energy profile is calculated along this reaction coordinate. This calculation includes quantization effects on bound vibrations but not on the reaction coordinate, and it is used to locate the variational transition state that defines a transition state ensemble. Then, the key interactions at the reactant, variational transition state, and product are analyzed in terms of both bond distances and electrostatic energies. The results of both analyses support the conclusion derived from previous mutational studies that the M20 loop of DHFR makes an important contribution to the electrostatic stabilization of the hydride transfer transition state. Third, transmission coefficients (including recrossing factors and multidimensional tunneling) are calculated and averaged over the transition state ensemble. These averaged transmission coefficients, combined with the quasithermodynamic free energy profile determined in stage 1, allow us to calculate rate constants, phenomenological free energies of activation, and primary and secondary kinetic isotope effects. A primary kinetic isotope effect (KIE) of 2.8 has been obtained, in good agreement with the experimentally determined value of 3.0 and with the value 3.2 calculated previously. The primary KIE is mainly a consequence of the quantization of bound vibrations. In contrast, the secondary KIE, with a value of 1.13, is almost entirely due to dynamical effects on the reaction coordinate, especially tunneling.

A point-charge force field for molecular mechanics simulations of proteins based on condensed-phase quantum mechanical calculations

Molecular mechanics models have been applied extensively to study the dynamics of proteins and nucleic acids. Here we report the development of a third-generation point-charge all-atom force field for proteins. Following the earlier approach of Cornell et al., the charge set was obtained by fitting to the electrostatic potentials of dipeptides calculated using B3LYP/cc-pVTZ//HF/6-31G** quantum mechanical methods. The main-chain torsion parameters were obtained by fitting to the energy profiles of Ace-Ala-Nme and Ace-Gly-Nme di-peptides calculated using MP2/cc-pVTZ//HF/6-31G** quantum mechanical methods. All other parameters were taken from the existing AMBER data base. The major departure from previous force fields is that all quantum mechanical calculations were done in the condensed phase with continuum solvent models and an effective dielectric constant of epsilon = 4. We anticipate that this force field parameter set will address certain critical short comings of previous force fields in condensed-phase simulations of proteins. Initial tests on peptides demonstrated a high-degree of similarity between the calculated and the statistically measured Ramanchandran maps for both Ace-Gly-Nme and Ace-Ala-Nme di-peptides. Some highlights of our results include (1) well-preserved balance between the extended and helical region distributions, and (2) favorable type-II poly-proline helical region in agreement with recent experiments. Backward compatibility between the new and Cornell et al. charge sets, as judged by overall agreement between dipole moments, allows a smooth transition to the new force field in the area of ligand-binding calculations. Test simulations on a large set of proteins are also discussed.

The role of alpha-, 3(10)-, and pi-helix in helix-->coil transitions

The conformational equilibrium between 3(10)- and alpha-helical structure has been studied via high-resolution NMR spectroscopy by Millhauser and coworkers using the MW peptide Ac-AMAAKAWAAKA AAARA-NH2. Their 750-MHz nuclear Overhauser effect spectroscopy (NOESY) spectra were interpreted to reflect appreciable populations of 3(10)-helix throughout the peptide, with the greatest contribution at the N and C termini. The presence of simultaneous alphaN(i,i + 2) and alphaN(i,i + 4) NOE cross-peaks was proposed to represent conformational averaging between 3(10)- and alpha-helical structures. In this study, we describe 25-nsec molecular dynamics simulations of the MW peptide at 298 K, using both an 8 A and a 10 A force-shifted nonbonded cutoff. The ensemble averages of both simulations are in reasonable agreement with the experimental helical content from circular dichroism (CD), the (3)J(HNalpha) coupling constants, and the 57 observed NOEs. Analysis of the structures from both simulations revealed very little formation of contiguous i --> i + 3 hydrogen bonds (3(10)-helix); however, there was a large population of bifurcated i --> i + 3 and i --> i + 4 alpha-helical hydrogen bonds. In addition, both simulations contained considerable populations of pi-helix (i --> i + 5 hydrogen bonds). Individual turns formed over residues 1-9, which we predict contribute to the intensities of the experimentally observed alphaN(i,i + 2) NOEs. Here we show how sampling of both folded and unfolded structures can provide a structural framework for deconvolution of the conformational contributions to experimental ensemble averages.

Importance of substrate and cofactor polarization in the active site of dihydrofolate reductase

By using a combined quantum-mechanical and molecular-mechanical potential in molecular dynamics simulations, we have investigated the effects of the enzyme electric field of dihydrofolate reductase on the electronic polarization of its 5-protonated dihydrofolate substrate at various stages of the catalyzed hydride transfer reaction. Energy decomposition of the total electrostatic interaction energy between the ligands and the enzyme shows that the polarization effect is 4% of the total electrostatic interaction energy, and, significantly, it accounts for 9kcal/mol of transition state stabilization relative to the reactant state. Therefore it is essential to take account of substrate polarization for quantitative interpretation of enzymatic function and for calculation of binding free energies of inhibitors to a protein. Atomic polarizations are calculated as the differences in the average atomic charges on the atoms in gas phase and in molecular simulations of the enzyme; this analysis shows that the glutamate tail and the pterin ring are the highly polarized regions of the substrate. Electron density difference plots of the reactant and product complexes at instantaneous configurations in the enzyme active center confirm the inferences made on the basis of partial atomic charges.

Breaking non-native hydrophobic clusters is the rate-limiting step in the folding of an alanine-based peptide

The formation mechanism of an alanine-based peptide has been studied by all-atom molecular dynamics simulations with a recently developed all-atom point-charge force field and the Generalize Born continuum solvent model at an effective salt concentration of 0.2M. Thirty-two simulations were conducted. Each simulation was performed for 100 ns. A surprisingly complex folding process was observed. The development of the helical content can be divided into three phases with time constants of 0.06-0.08, 1.4-2.3, and 12-13 ns, respectively. Helices initiate extreme rapidly in the first phase similar to that estimated from explicit solvent simulations. Hydrophobic collapse also takes place in this phase. A folding intermediate state develops in the second phase and is unfolded to allow the peptide to reach the transition state in the third phase. The folding intermediate states are characterized by the two-turn short helices and the transition states are helix-turn-helix motifs-both of which are stabilized by hydrophobic clusters. The equilibrium helical content, calculated by both the main-chain Phi-Psi torsion angles and the main-chain hydrogen bonds, is 64-66%, which is in remarkable agreement with experiments. After corrected for the solvent viscosity effect, an extrapolated folding time of 16-20 ns is obtained that is in qualitative agreement with experiments. Contrary to the prevailing opinion, neither initiation nor growth of the helix is the rate-limiting step. Instead, the rate-limiting step for this peptide is breaking the non-native hydrophobic clusters in order to reach the transition state. The implication to the folding mechanisms of proteins is also discussed.

Recent advances in helix-coil theory

Peptide helices in solution form a complex mixture of all helix, all coil or, most frequently, central helices with frayed coil ends. In order to interpret experiments on helical peptides and make theoretical predictions on helices, it is therefore essential to use a helix-coil theory that takes account of this equilibrium. The original Zimm-Bragg and Lifson-Roig helix-coil theories have been greatly extended in the last 10 years to include additional interactions. These include preferences for the N-cap, N1, N2, N3 and C-cap positions, capping motifs, helix dipoles, side chain interactions and 3(10)-helix formation. These have been applied to determine energies for these preferences from experimental data and to predict the helix contents of peptides. This review discusses these newly recognised structural features of helices and how they have been included in helix-coil models.

Helix-coil transitions re-visited

The thermally-induced helix-coil transition in polyamino acids is a good model for determining the helix-forming propensities of amino acids but not for the two-state folding/unfolding transition in globular proteins. The equilibrium and kinetic treatments of the helix-coil transition are summarized here together with a description of applications to various types of homopolymers and copolymers. Attention is then focused on the helix-coil transition in poly-L-alanine as an example of a non-polar polyamino acid. To render such a non-polar polymer water soluble, it is necessary to introduce polar amino acids such as lysines, but care must be taken as to the location of such polar residues. If they are attached as end groups, as in a triblock copolymer, they do not perturb the helix-forming tendency of the central poly-L-alanine block significantly, but if they are introduced within the sequence of alanine residues, then the hydration properties of the lysines dominate the behavior of the resulting copolymer, thereby leading to erroneous values of the parameters characterizing the helix-forming tendency of the alanines. Neutral but polar residues, such as glutamines, also exhibit hydration-dominating properties but less so than charged lysines. Some details of the calculations for an alanine/glutamine copolymer are presented here. It is concluded that random copolymers based on a neutral water-soluble host provide reliable information about the helix-forming tendencies of amino acid residues that are introduced as guests among such neutral host residues.

Measuring dynamic flexibility of the coil state of a helix-forming peptide

To investigate the dynamic flexibility of the coil state of a helix-forming peptide the end-to-end contact rate was determined. Nanosecond optical excitation of tryptophan at one end of a 22 residue, alanine peptide populates a long-lived triplet state which is quenched upon close contact with a cyclic disulfide attached to the opposite end. Analysis of the decay of the triplet population using a two-state model for helix formation yields the diffusion-limited end-to-end contact rate of the coil state of the peptide as well as the helix-->coil and coil-->helix rates. The helix-coil rates are very similar to those previously measured in laser temperature-jump experiments. The end-to-end contact rate of 1.1 x 10(7) s(-1) in the coil state is tenfold faster than the rate for a disordered peptide with threonine substituted for alanine and, somewhat surprisingly, is about twice the rate for a disordered glycine-containing peptide. These differences are discussed in terms of the theory of Szabo, Schulten and Schulten. The rates should provide important new benchmarks for testing the accuracy of atomistic molecular dynamics simulations.

Occurrence, conformational features and amino acid propensities for the pi-helix

The most abundant helix type in proteins is the alpha-helix, accounting for about 31% of amino acid secondary structure states, while the 3(10)-helix accounts for about 4%. The pi-helix appears to be extremely rare and is considered to be unstable. Existing secondary structure definition methods find very few within the Protein Data Bank. Using an improved pi-helix definition algorithm to search a non-redundant subset of high-resolution and well-refined protein structures, we found that almost every tenth protein contained a pi-helix. This enabled us to show for the first time that the pi-helix has structural parameters that are different from the hypothesized model values. It also has distinctive amino acid preferences and it is conserved within functionally related proteins. Features that may contribute to the stability of the pi-helical structure have also been identified. In addition to hydrogen bonds, several other factors contribute to the stability of pi-helices. The pi-helix may have some functional advantages over other helical structures. Thus, we describe cases where the side chains of functionally important residues at every fourth position within a pi-helix could be aligned and brought close together in a way that would not be allowed by any other helix type.

Folding study of an Aib-rich peptide in DMSO by molecular dynamics simulations

To evaluate the ability of molecular dynamics (MD) simulations using atomic force-fields to correctly predict stable folded conformations of a peptide in solution, we show results from MD simulations of the reversible folding of an octapeptide rich in alpha-aminoisobutyric acid (2-amino-2-methyl-propanoic acid, Aib) solvated in di-methyl-sulfoxide (DMSO). This solvent generally prevents the formation of secondary structure, whereas Aib-rich peptides show a high propensity to form secondary structural elements, in particular 3(10)- and alpha-helical structures. Aib is, moreover, achiral, so that Aib-rich peptides can form left- or right-handed helices depending on the overall composition of the peptide, the temperature, and the solvation conditions. This makes the system an interesting case to study the ensembles of peptide conformations as a function of temperature by MD simulation. Simulations involving the folding and unfolding of the peptide were performed starting from two initial structures, a right-handed alpha-helical structure and an extended structure, at three temperatures, 298 K, 340 K, and 380 K, and the results are compared with experimental nuclear magnetic resonance (NMR) data measured at 298 K and 340 K. The simulations generally reproduce the available experimental nuclear Overhauser effect (NOE) data, even when a wide range of conformations is sampled at each temperature. The importance of adequate statistical sampling in order to reliably interpret the experimental data is discussed.

Hydrogen bonding in helical polypeptides from molecular dynamics simulations and amide hydrogen exchange analysis: alamethicin and melittin in methanol

Molecular dynamics simulations of ion channel peptides alamethicin and melittin, solvated in methanol at 27 degrees C, were run with either regular alpha-helical starting structures (alamethicin, 1 ns; melittin 500 ps either with or without chloride counterions), or with the x-ray crystal coordinates of alamethicin as a starting structure (1 ns). The hydrogen bond patterns and stabilities were characterized by analysis of the dynamics trajectories with specified hydrogen bond angle and distance criteria, and were compared with hydrogen bond patterns and stabilities previously determined from high-resolution NMR structural analysis and amide hydrogen exchange measurements in methanol. The two alamethicin simulations rapidly converged to a persistent hydrogen bond pattern with a high level of 3(10) hydrogen bonding involving the amide NH's of residues 3, 4, 9, 15, and 18. The 3(10) hydrogen bonds stabilizing amide NH's of residues C-terminal to P2 and P14 were previously proposed to explain their high amide exchange stabilities. The absence, or low levels of 3(10) hydrogen bonds at the N-terminus or for A15 NH, respectively, in the melittin simulations, is also consistent with interpretations from amide exchange analysis. Perturbation of helical hydrogen bonding in the residues before P14 (Aib10-P14, alamethicin; T11-P14, melittin) was characterized in both peptides by variable hydrogen bond patterns that included pi and gamma hydrogen bonds. The general agreement in hydrogen bond patterns determined in the simulations and from spectroscopic analysis indicates that with suitable conditions (including solvent composition and counterions where required), local hydrogen-bonded secondary structure in helical peptides may be predicted from dynamics simulations from alpha-helical starting structures. Each peptide, particularly alamethicin, underwent some large amplitude structural fluctuations in which several hydrogen bonds were cooperatively broken. The recovery of the persistent hydrogen bonding patterns after these fluctuations demonstrates the stability of intramolecular hydrogen-bonded secondary structure in methanol (consistent with spectroscopic observations), and is promising for simulations on extended timescales to characterize the nature of the backbone fluctuations that underlie amide exchange from isolated helical polypeptides.

Insight into signal transduction: structural alterations in transmembrane helices probed by multi-1 ns molecular dynamics simulations

The hypothesis of structural alteration in transmembrane helices for signal transduction process is viewed by molecular dynamics simulation techniques. For the c-erbB-2 transmembrane domain involved in oncogenicity, the occurrence of conformational changes has been previously described as transition from the alpha to pi helix. This dynamical feature is thoroughly analyzed for the wild phenotype and oncogenic sequences from a series of 18 simulations carried out on one nanosecond time scale. We show that these structural events do not depend upon the conditions of simulations like force field or starting helix coordinates. We demonstrate that the oncogenic mutations Val659 Glu, Gln and Asp do not prevent the transition. Furthermore, we show that beta branched residues, in conjunction with Gly residues in the c-erbB-2 sequence, act as destabilizers for the alpha helix structure, pi deformations are tightly related to other local structural motifs found in soluble and membrane proteins. These structural alterations are discussed in term of structure-activity relationships for the c-erbB-2 activating mechanism mediated by transmembrane domain dimerization.