Effects of long-range electrostatic forces on simulated protein folding kinetics

Molecular dynamics simulations are a useful tool for characterizing protein folding pathways. There are several methods of treating electrostatic forces in these simulations with varying degrees of physical fidelity and computational efficiency. In this article, we compare the reaction field (RF) algorithm, particle-mesh Ewald (PME), and tapered cutoffs with increasing cutoff radii to address the impact of the electrostatics method employed on the folding kinetics. We quantitatively compare different methods by a correlation of quantitative measures of protein folding kinetics. The results of these comparisons show that for protein folding kinetics, the RF algorithm can quantitatively reproduce the kinetics of the more costly PME algorithm. These results not only assist the selection of appropriate algorithms for future simulations, but also give insight on the role that long-range electrostatic forces have in protein folding.     

Dynamical motions of lipids and a finite size effect in simulations of bilayers

Molecular dynamics (MD) simulations of dipalmitoylphosphatidylcholine bilayers composed of 72 and 288 lipids are used to examine system size dependence on dynamical properties associated with the particle mesh Ewald (PME) treatment of electrostatic interactions. The lateral diffusion constant Dl is 2.92 x 10(-7) and 0.95 x 10(-7) cm2/s for 72 and 288 lipids, respectively. This dramatic finite size effect originates from the correlation length of lipid diffusion, which extends to next-nearest neighbors in the 288 lipid system. Consequently, diffusional events in smaller systems can propagate across the boundaries of the periodic box. The internal dynamics of lipids calculated from the PME simulations are independent of the system size. Specifically, reorientational correlation functions for the slowly relaxing phosphorus-glycerol hydrogen, phosphorus-nitrogen vectors, and more rapidly relaxing CH vectors in the aliphatic chains are equivalent for the 72 and 288 lipid simulations. A third MD simulation of a bilayer with 72 lipids using spherical force-shift electrostatic cutoffs resulted in interdigitated chains, thereby rendering this cutoff method inappropriate.     

The influence of simulation conditions in molecular dynamics investigations of model β-sheet peptides

The influence of simulation methods, cutoff based and particle mesh Ewald (PME) on the accuracy by which experimentally derived nuclear Overhauser effect (NOE) data are reproduced, has been investigated using 500-ns-long molecular dynamics simulations on a model -sheet peptide in explicit solvent. The structural and conformational features under the different conditions were evaluated in terms of flexibility, secondary structure content, hydrogen-bonding pattern and percent of native contacts as a function of time. It was found that the different simulation methods strongly influence the dynamics of the peptide, confirming previous observations based on ideal peptide models simulated for much shorter times. Moreover, the results of our simulations prove once more that it is necessary to reach extremely long time scales to obtain enough statistics to accurately reproduce experimental NOE restraints even in the case of the PME method, despite its tendency to the stabilization of conformations which are structurally closely related to the ones derived through experiment. Possible implications regarding the stabilization and folding mechanisms, together with their relationship to the experimental study of peptide models, are discussed.     

Beta-hairpin folding by a model amyloid peptide in solution and at an interface

The development of specific agents against amyloidoses requires an understanding of the conformational distribution of fibrillogenic peptides at a microscopic level. Here, I present molecular dynamics simulations of the model amyloid peptide LSFD with sequence LSFDNSGAITIG-NH2 in explicit water and at a water/vapor interface for a total time scale of approximately 1.8 micros. An extended structure was used as initial peptide configuration. At approximately 290 K, solvated LSFD was kinetically trapped in diverse misfolded beta-sheet/coil conformations. At 350 K, in contrast, the same type II' beta-hairpin in equilibrium with less ordered but also U-shaped conformations was observed for the core residues DNSGAITI in solution and at the interface in multiple independent simulations. The most stable structural unit of the beta-hairpin was the two residue turn (GA). The core residues exhibited a well-defined folded state in which the beta-hairpin was stabilized by a hydrogen bond between the side chain of Asn-385 and the main chain carbonyl group of Gly-387. My results suggest that beta-sheet conformations indicated from previous Fourier-transform infrared spectroscopy measurements immediately after preparation of the peptide solution may not arise from protofilaments as speculated by others but are a property of LSFD monomers. In addition, combined with previous results from X-ray scattering, my findings suggest that interfacial aggregation of LSFD implies a transition from U-shaped to extended peptide conformations. This work including the first simulations of reversible beta-hairpin folding at an interface is an essential step toward a microscopic understanding of interfacial peptide folding and self-assembly. Knowledge of the main conformation of the peptide core may facilitate the design of possible inhibitors of LSFD aggregation as a test ground for future computational therapeutic strategies against amyloid diseases.     

Absence of reptation in the high-temperature folding of the trpzip2 beta-hairpin peptide

We have carried out extensive all atom explicit solvent simulations of the high-temperature folding and unfolding of the trpzip2 beta-hairpin peptide and examined the resulting trajectories for evidence of folding via a reptation mechanism. Over 300 microcanonical simulations of 10 ns each were initiated from a Boltzmann ensemble of conformations at 425 K. Though we observed numerous folding and unfolding events, no evidence of reptation was found. The diffusional dynamics of the peptide are orders of magnitude faster than any observed reptation-like motion. Our data suggest that the dominant mechanisms for beta-hairpin folding under these conditions are hydrophobic collapse and turn formation, and that rearrangements occur via significant expansion of the polypeptide chain.     

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

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

Two-dimensional infrared spectroscopy as a probe of the solvent electrostatic field for a twelve residue peptide

The linear IR and two-dimensional (2D) IR spectra of the amide-I modes of the 12-residue beta-hairpin peptide tryptophan zipper-2 (SWTWENGKWTWK) and its two 13C isotopomers were simulated, with local mode frequencies evaluated by two solution-phase peptide amide-I frequency maps proposed recently: an electrostatic potential map and an electrostatic field map. Both maps predict a set of nondegenerate local amide-I mode transition energies for the hairpin. Spectral simulations using both maps predict the main spectral features of the linear IR and 2D IR experimental results of the (13)C-labeled and -unlabeled hairpin. The radial distribution functions obtained using trajectories from classical molecular dynamics simulations demonstrate different water distributions at different sites of the hairpin. Our results suggest that the observed difference of the (13)C-shifted band, including its peak position and frequency distributions for different isotopomers, in both linear IR and 2D IR spectra, is likely to be due to the difference in the local environment of the solvated peptide. Ab initio density functional theory calculations show a residue-independent (13)C shift of the amide-I mode, further supporting the result. The variations of these shifts are attributed to the residue level heterogeneity of the electrostatic environment of the peptide. Our results show that 2D IR of peptide with single (13)C isotopic labeling can be used to probe the electrostatic environment of the peptide local structure.     

Molecular dynamics study of peptides in implicit water: ab initio folding of beta-hairpin, beta-sheet, and beta beta alpha-motif

In this communication, we have demonstrated that molecular dynamics simulations using a GB implicit solvation model with the all-atom based force field (CHARMM19) can describe the spontaneous folding of small peptides in aqueous solution. The native structures of peptides with various structural motifs (beta-hairpin, beta-sheet, and betabetaalpha-moiety) were successfully predicted within reasonable time scales by MD simulations at moderately elevated temperatures. It is expected that the present simulations provide further insight into mechanism/pathways of the peptide folding.     

Effect of different treatments of long-range interactions and sampling conditions in molecular dynamic simulations of rhodopsin embedded in a dipalmitoyl phosphatidylcholine bilayer

The present study analyzes the effect of the simulation conditions on the results of molecular dynamics simulations of G-protein coupled receptors (GPCRs) performed with an explicit lipid bilayer. Accordingly, the present work reports the analysis of different simulations of bovine rhodopsin embedded in a dipalmitoyl phosphatidylcholine (DPPC) lipid bilayer using two different sampling conditions and two different approaches for the treatment of long-range electrostatic interactions. Specifically, sampling was carried out either by using the statistical ensembles NVT or NPT (constant number of atoms, a pressure of 1 atm in all directions and fixed temperature), and the electrostatic interactions were treated either by using a twin-cutoff, or the particle mesh Ewald summation method (PME). The results of the present study suggest that the use of the NPT ensemble in combination with the PME method provide more realistic simulations. The use of NPT during the equilibration avoids the need of an a priori estimation of the box dimensions, giving the correct area per lipid. However, once the system is equilibrated, the simulations are irrespective of the sampling conditions used. The use of an electrostatic cutoff induces artifacts on both lipid thickness and the ion distribution, but has no direct effect on the protein and water molecules.     

beta-Hairpins,alpha-helices,and the intermediates among the secondary structures in the energy landscape of a peptide from a distal beta-hairpin of SH3 domain

Energy landscape of a peptide, extracted from a distal beta-hairpin of src SH3 domain, in explicit water was obtained with the multicanonical molecular dynamics. A variety of beta-hairpins with various strand-strand hydrogen bonds were found in the energy landscape at 300 K. There was no energy barrier between random-coil and hairpins. Thus, the peptide conformation can easily change from the random-coil to the hairpins in the thermal fluctuations at 300 K. The landscape also included two clusters of alpha-helices, among which an energy barrier existed, and besides, these helix clusters were separated from the other conformations. Thus, the free-energy barrier exists among the helices and the other conformations. Intermediate clusters were found between the helix and the hairpin clusters. The current study showed that the isolated state of this peptide in water fluctuates among random-coil, beta-hairpin, and alpha-helix. In SH3 domain, which has a topology of mainly beta-protein, the whole-protein folding may proceed when the segment is folded in the beta-hairpin and the other parts of the protein are coupled with the beta-hairpin in an energetically or kinetically favorite way.     

Structural, dynamic, and electrostatic properties of fully hydrated DMPC bilayers from molecular dynamics simulations accelerated with graphical processing units (GPUs)

We present results of molecular dynamics simulations of fully hydrated DMPC bilayers performed on graphics processing units (GPUs) using current state-of-the-art non-polarizable force fields and a local GPU-enabled molecular dynamics code named FEN ZI. We treat the conditionally convergent electrostatic interaction energy exactly using the particle mesh Ewald method (PME) for solution of Poisson's Equation for the electrostatic potential under periodic boundary conditions. We discuss elements of our implementation of the PME algorithm on GPUs as well as pertinent performance issues. We proceed to show results of simulations of extended lipid bilayer systems using our program, FEN ZI. We performed simulations of DMPC bilayer systems consisting of 17,004, 68,484, and 273,936 atoms in explicit solvent. We present bilayer structural properties (atomic number densities, electron density profiles), deuterium order parameters (S(CD)), electrostatic properties (dipole potential, water dipole moments), and orientational properties of water. Predicted properties demonstrate excellent agreement with experiment and previous all-atom molecular dynamics simulations. We observe no statistically significant differences in calculated structural or electrostatic properties for different system sizes, suggesting the small bilayer simulations (less than 100 lipid molecules) provide equivalent representation of structural and electrostatic properties associated with significantly larger systems (over 1000 lipid molecules). We stress that the three system size representations will have differences in other properties such as surface capillary wave dynamics or surface tension related effects that are not probed in the current study. The latter properties are inherently dependent on system size. This contribution suggests the suitability of applying emerging GPU technologies to studies of an important class of biological environments, that of lipid bilayers and their associated integral membrane proteins. We envision that this technology will push the boundaries of fully atomic-resolution modeling of these biological systems, thus enabling unprecedented exploration of meso-scale phenomena (mechanisms, kinetics, energetics) with atomic detail at commodity hardware prices.     

Aggregation of polyalanine in a hydrophobic environment

The dimerization of polyalanine peptides in a hydrophobic environment was explored using replica exchange molecular dynamics simulations. A nonpolar solvent (cyclohexane) was used to mimic, among other hydrophobic environments, the hydrophobic interior of a membrane in which the peptides are fully embedded. Our simulations reveal that while the polyalanine monomer preferentially adopts a beta-hairpin conformation, dimeric phases exist in an equilibrium between random coil, alpha-helical, beta-sheet, and beta-hairpin states. A thermodynamic characterization of the dimeric phases reveals that electric dipole-dipole interactions and optimal side-chain packing stabilize alpha-helical conformations, while hydrogen bond interactions favor beta-sheet conformations. Possible pathways leading to the formation of alpha-helical and beta-sheet dimers are discussed.     

Simulation of infrared spectra for beta-hairpin peptides stabilized by an Aib-Gly turn sequence: correlation between conformational fluctuation and vibrational coupling

Vibrational spectra of a 12-residue beta-hairpin peptide, RYVEVBGKKILQ (HBG), stabilized by an Aib-Gly turn sequence (B = Aib) were investigated theoretically using a combination of molecular dynamics (MD) and density functional theory (DFT) calculations. Selected conformations of HBG were extracted from a classical MD trajectory and used for spectral simulations. DFT calculations, based on the Cartesian coordinate spectral property transfer protocol, were carried out for peptide structures in which all residues are replaced with Ala, except for the Aib and Gly residues, but the backbone (phi, psi, omega) structure of the original configuration is retained. The simulations provide a basis for interpretation of the HBG amide I infrared spectra in terms of structural variables such as detailed secondary structure and thermal conformational fluctuation as well as vibrational coupling as indicated by spectra of 13C isotope-labeled variants. The characteristic amide I band shape of such small beta-hairpin peptides appears to arise from the structure of the short antiparallel beta-sheet strands. The role of structural parameter fluctuation in vibrational coupling is evaluated by comparison of DFT-derived amide coupling constants for selected configurations and from transition dipole coupling calculations of coupling parameters between (13)C isotopically labeled residues for a MD-derived ensemble of configurations. Calculated results were compared with the experimentally obtained spectra for several (13)C isotope-labeled peptides of this sequence.     

Polypeptide folding using Monte Carlo sampling, concerted rotation, and continuum solvation

An efficient concerted rotation algorithm for use in Monte Carlo statistical mechanics simulations is applied to fold three polypeptides, U(1-17)T9D, alpha(1), and trpzip2, which exhibit native beta-hairpin and alpha-helix folds. The method includes flexible bond and dihedral angles, and a Gaussian bias is applied with driver bond and dihedral angles to optimize the sampling efficiency. Solvation in water is implemented with the generalized Born (GBSA) model. The computed lowest-energy manifolds for the folded structures of the two beta-hairpins agree closely with the corresponding NMR structures. In the case of the alpha(1) peptide, the folded alpha-helical state, which is observed as oligomers in concentrated solution and crystals, is not stable in isolation. The computed preference for random coil structures is in agreement with NMR experiments at low concentration. The fact that native states can be located on high dimensional energy surfaces starting from extended conformations shows that the present methodology samples all relevant parts of the conformational space. The OPLS-AA force field with the GBSA solvent model was also found to perform well in leading to clear energetic separation of the correctly folded structures from misfolded structures for the two peptides that form beta-turns.     

Comparison of Secondary Structure Formation Using 10 Different Force Fields in Microsecond Molecular Dynamics Simulations

We have compared molecular dynamics (MD) simulations of a β-hairpin forming peptide derived from the protein Nrf2 with 10 biomolecular force fields using trajectories of at least 1 μs. The total simulation time was 37.2 μs. Previous studies have shown that different force fields, water models, simulation methods, and parameters can affect simulation outcomes. The MD simulations were done in explicit solvent with a 16-mer Nrf2 β-hairpin forming peptide using Amber ff99SB-ILDN, Amber ff99SB*-ILDN, Amber ff99SB, Amber ff99SB*, Amber ff03, Amber ff03*, GROMOS96 43a1p, GROMOS96 53a6, CHARMM27, and OPLS-AA/L force fields. The effects of charge-groups, terminal capping, and phosphorylation on the peptide folding were also examined. Despite using identical starting structures and simulation parameters, we observed clear differences among the various force fields and even between replicates using the same force field. Our simulations show that the uncapped peptide folds into a native-like β-hairpin structure at 310 K when Amber ff99SB-ILDN, Amber ff99SB*-ILDN, Amber ff99SB, Amber ff99SB*, Amber ff03, Amber ff03*, GROMOS96 43a1p, or GROMOS96 53a6 were used. The CHARMM27 simulations were able to form native hairpins in some of the elevated temperature simulations, while the OPLS-AA/L simulations did not yield native hairpin structures at any temperatures tested. Simulations that used charge-groups or peptide capping groups were not largely different from their uncapped counterparts with single atom charge-groups. On the other hand, phosphorylation of the threonine residue located at the β-turn significantly affected the hairpin formation. To our knowledge, this is the first study comparing such a large set of force fields with respect to β-hairpin folding. Such a comprehensive comparison will offer useful guidance to others conducting similar types of simulations.     

The Thermodynamics of Folding of a β Hairpin Peptide Probed Through Replica Exchange Molecular Dynamics Simulations

Peptides that possess a well defined native state are ideal model systems to study the folding of proteins. They possess many of the complexities of larger proteins, yet their small size renders their study computationally tractable. Recent advances in sampling techniques, including replica exchange molecular dynamics, now permit a full characterization of the thermodynamics of folding of small peptides. These simulations not only yield insight into the folding of larger proteins, but equally importantly, they allow, through comparison with experiment, an objective test of the accuracy of force fields, water models and of different numerical schemes for dealing with electrostatic interactions. In this account, we present a molecular dynamics simulation of a small β-hairpin peptide using the replica exchange algorithm and illustrate how this enhanced sampling scheme enables a thorough characterization of the native and unfolded states, and sheds new light into its folding mechanism.     

Application of MDGRAPE-3, a special purpose board for molecular dynamics simulations, to periodic biomolecular systems

We describe the application of a special purpose board for molecular dynamics simulations, named MDGRAPE-3, to the problem of simulating periodic bio-molecular systems. MDGRAPE-3 is the latest board in a series of hardware accelerators designed to calculate the nonbonding long-range interactions much more rapidly than normal processors. So far, MDGRAPEs were mainly applied to isolated systems, where very many nonbonded interactions were calculated without any distance cutoff. However, in order to regulate the density and pressure during simulations of membrane embedded protein systems, one has to evaluate interactions under periodic boundary conditions. For this purpose, we implemented the Particle-Mesh Ewald (PME) method, and its approximation with distance cutoffs and charge neutrality as proposed by Wolf et al., using MDGRAPE-3. When the two methods were applied to simulations of two periodic biomolecular systems, a single MDGRAPE-3 achieved 30-40 times faster computation times than a single conventional processor did in the both cases. Both methods are shown to have the same molecular structures and dynamics of the systems.     

Structure elucidation and 3D solution conformation of the antibiotic enduracidin determined by NMR spectroscopy and molecular dynamics

Enduracidin and ramoplanin belong to the large family of cyclodepsipeptide antibiotics, highly effective against Gram-positive bacteria. The primary and 3D solution structure of ramoplanin is already well known, and the primary structure of enduracidin has been determined by a combination of chemical and NMR spectroscopic methods. Both antibiotics share a similar peptide core of 17 amino acids and differ mainly in the length of the acyl chain and the presence of two D-mannose moieties in ramoplanin. Based on the high sequence homology with ramoplanin, the structure in solution of enduracidin is modeled as a cyclic peptide. The tertiary structure thus obtained was refined through molecular dynamics (MD) simulation, in which the interatomic NOE-derived distance restraints were imposed. MD simulations yielded a family of representative 3D structures (RMSD = 0.89), which highlighted a backbone geometry similar to that of ramoplanin in its beta-hairpin arrangement. In contrast, enduracidin displays a different arrangement of the side-chain and of the residues forming the hydrophobic core.     

Theoretical characterization of alpha-helix and beta-hairpin folding kinetics

By means of the conformational free energy surface and corresponding diffusion coefficients, as obtained by long time scale atomistic molecular dynamics simulations (mus time scale), we model the folding kinetics of alpha-helix and beta-hairpin peptides as a diffusive process over the free energy surface. The two model systems studied in this paper (the alpha-helical temporin L and the beta-hairpin prion protein H1 peptide) exhibit a funnel-like almost barrierless free energy profile, leading to nonexponential folding kinetics matching rather well the available experimental data. Moreover, using the free energy profile provided by Mu?oz et al. [Mu?oz et al. Nature 1997, 390: 196-199], this model was also applied to reproduce the two-state folding kinetics of the C-terminal beta-hairpin of protein GB1, yielding an exponential folding kinetics with a time constant (approximately 5 micros) in excellent agreement with the experimentally observed one (approximately 6 micros). Finally, the folding kinetics obtained by solving the diffusion equation, considering either a one-dimensional or a two-dimensional free energy surface, are also compared in order to understand the relevance of the possible kinetic coupling between conformational degrees of freedom in the folding process.     

Dependency of ligand free energy landscapes on charge parameters and solvent models

We performed replica-exchange molecular dynamics (REMD) simulations of six ligands to examine the dependency of their free energy landscapes on charge parameters and solvent models. Six different charge parameter sets for each ligand were first generated by RESP and AM1-BCC methods using three different conformations independently. RESP charges showed some conformational dependency. On the other hand, AM1-BCC charges did not show conformational dependency and well reproduced the overall trend of RESP charges. The free energy landscapes obtained from the REMD simulations of ligands in vacuum, Generalized-Born (GB), and TIP3P solutions were then analyzed. We found that even small charge differences can produce qualitatively different landscapes in vacuum condition, but the differences tend to be much smaller under GB and TIP3P conditions. The simulations in the GB model well reproduced the landscapes in the TIP3P model using only a fraction of the computational cost. The protein-bound ligand conformations were rarely the global minimum states, but similar conformations were found to exist in aqueous solution without proteins in regions close to the global minimum, local minimum or intermediate states.