Exploring peptide energy landscapes: a test of force fields and implicit solvent models

Investigation of the Effect of Peptide p5 Targeting CDK5-p25 Hyperactivity on Munc18-1 (P67) Regulating Neuronal Exocytosis Using Molecular Simulations

Munc18-1 is an SM (sec1/munc-like) family protein involved in vesicle fusion and neuronal exocytosis. Munc18-1 is known to regulate the exocytosis process by binding with closed- and open-state conformations of Syntaxin1, a protein belonging to the SNARE family established to be central to the exocytosis process. Our previous work studied peptide p5 as a promising drug candidate for CDK5-p25 complex, an Alzheimer's disease (AD) pathological target. Experimental in vivo and in vitro studies suggest that Munc18-1 promotes p5 to selectively inhibit the CDK5-p25 complex without affecting the endogenous CDK5 activity, a characteristic of remarkable therapeutic implications. In this paper, we identify several binding modes of p5 with Munc18-1 that could potentially affect the Munc18-1 binding with SNARE proteins and lead to off-target effects on neuronal communication using molecular dynamics simulations. Recent studies indicate that disruption of Munc18-1 function not only disrupts neurotransmitter release but also results in neurodegeneration, exhibiting clinical resemblance to other neurodegenerative conditions such as AD, causing diagnostic and treatment challenges. We characterize such interactions between p5 and Munc18-1, define the corresponding pharmacophores, and provide guidance for the in vitro validation of our findings to improve therapeutic efficacy and safety of p5.

Computational Study of the Allosteric Effects of p5 on CDK5-p25 Hyperactivity as Alternative Inhibitory Mechanisms in Neurodegeneration

The cyclin-dependent kinase (CDK5) forms a stable complex with its activator p25, leading to the hyperphosphorylation of tau proteins and to the formation of plaques and tangles that are considered to be one of the typical causes of Alzheimer's disease (AD). Hence, the pathological CDK5-p25 complex is a promising therapeutic target for AD. Small peptides, obtained from the truncation of CDK5 physiological activator p35, have shown promise in inhibiting the pathological complex effectively while also crossing the blood-brain barrier. One such small 24-residue peptide, p5, has shown selective inhibition toward the pathological complex in vivo. Our previous research focused on the characterization of a computationally predicted CDK5-p5 binding mode and of its pharmacophore, which was consistent with competitive inhibition. In continuation of our previous work, herein, we investigate four additional binding modes to explore other possible mechanisms of interaction between CDK5 and p5. The quantitative description of the pharmacophore is consistent with both competitive and allosteric p5-induced inhibition mechanisms of CDK5-p25 pathology. The gained insights can direct further in vivo/in vitro tests and help design small peptides, linear or cyclic, or peptidomimetic compounds as adjuvants of orthosteric inhibitors or as part of a cocktail of drugs with enhanced effectiveness and lower side effects.

Peptide and Protein Structure Prediction with a Simplified Continuum Solvent Model

A continuum solvent model based on screened Coulomb potentials has been simplified and parametrized to sample native-like structures in replica-exchange simulations of each of six different peptides and miniproteins. Low-energy, native, and non-native structures were used to iteratively refine 11 parameter values. The centroid of the largest cluster of structures sampled in simulations initiated from an extended conformation represents the predicted structure. The main-chain rms deviation of this prediction from the experimental structure was 0.47 Å for the 12-residue Trp-zip2, 0.86 Å for the 14-residue MBH12, 2.53 Å for the 17-residue U(1-17)T9D, 2.03 Å for the 20-residue BS1, 1.08 Å for the 20-residue Trp-cage, and 3.64 Å for the 35-residue villin headpiece subdomain HP35. The centroid of the sixth largest cluster sampled for HP35 deviated by 0.91 Å. The CHARMM22/CMAP force field was used, with an additional ψ torsion term for residues other than glycine and proline. Six parameters govern the dielectric response of the continuum solvent, and four values of surface tension approximate nonpolar effects. An atom's self-energy and interaction energies are screened independently, each depending on whether the atom is part of a charged group, a neutral hydrogen-bonding main-chain group, or any other neutral group. The parameters inferred result in strong main-chain hydrogen bonds, consistent with the view that protein folding is dominated by the formation of these bonds. (1,2) Conformations of MBH12 and BS1 were excluded from the energy-function refinement, suggesting the parameters, referred to as SCP18, are transferable. An efficient estimate of solvent-accessible surface area is also described.

Computational study of the inhibitory mechanism of the kinase CDK5 hyperactivity by peptide p5 and derivation of a pharmacophore

The hyperactivity of the cyclic dependent kinase 5 (CDK5) induced by the activator protein p25 has been linked to a number of pathologies of the brain. The CDK5-p25 complex has thus emerged as a major therapeutic target for Alzheimer's disease (AD) and other neurodegenerative conditions. Experiments have shown that the peptide p5 reduces the CDK5-p25 activity without affecting the endogenous CDK5-p35 activity, whereas the peptide TFP5, obtained from p5, elicits similar inhibition, crosses the blood-brain barrier, and exhibits behavioral rescue of AD mice models with no toxic side effects. The molecular basis of the kinase inhibition is not currently known, and is here investigated by computer simulations. It is shown that p5 binds the kinase at the same CDK5/p25 and CDK5/p35 interfaces, and is thus a non-selective competitor of both activators, in agreement with available experimental data in vitro. Binding of p5 is enthalpically driven with an affinity estimated in the low µM range. A quantitative description of the binding site and pharmacophore is presented, and options are discussed to increase the binding affinity and selectivity in the design of drug-like compounds against AD.

Detection and characterization of nonspecific, sparsely populated binding modes in the early stages of complexation

A method is proposed to study protein-ligand binding in a system governed by specific and nonspecific interactions. Strong associations lead to narrow distributions in the proteins configuration space; weak and ultraweak associations lead instead to broader distributions, a manifestation of nonspecific, sparsely populated binding modes with multiple interfaces. The method is based on the notion that a discrete set of preferential first-encounter modes are metastable states from which stable (prerelaxation) complexes at equilibrium evolve. The method can be used to explore alternative pathways of complexation with statistical significance and can be integrated into a general algorithm to study protein interaction networks. The method is applied to a peptide-protein complex. The peptide adopts several low-population conformers and binds in a variety of modes with a broad range of affinities. The system is thus well suited to analyze general features of binding, including conformational selection, multiplicity of binding modes, and nonspecific interactions, and to illustrate how the method can be applied to study these problems systematically. The equilibrium distributions can be used to generate biasing functions for simulations of multiprotein systems from which bulk thermodynamic quantities can be calculated.

Variational Optimization of an All-Atom Implicit Solvent Force Field to Match Explicit Solvent Simulation Data

The development of accurate implicit solvation models with low computational cost is essential for addressing many large-scale biophysical problems. Here, we present an efficient solvation term based on a Gaussian solvent-exclusion model (EEF1) for simulations of proteins in aqueous environment, with the primary aim of having a good overlap with explicit solvent simulations, particularly for unfolded and disordered states - as would be needed for multiscale applications. In order to achieve this, we have used a recently proposed coarse-graining procedure based on minimization of an entropy-related objective function to train the model to reproduce the equilibrium distribution obtained from explicit water simulations. Via this methodology, we have optimized both a charge screening parameter and a backbone torsion term against explicit solvent simulations of an α-helical and a β-stranded peptide. The performance of the resulting effective energy function, termed EEF1-SB, is tested with respect to the properties of folded proteins, the folding of small peptides or fast-folding proteins, and NMR data for intrinsically disordered proteins. The results show that EEF1-SB provides a reasonable description of a wide range of systems, but its key advantage over other methods tested is that it captures very well the structure and dimension of disordered or weakly structured peptides. EEF1-SB is thus a computationally inexpensive (~ 10 times faster than Generalized-Born methods) and transferable approximation for treating solvent effects.

Specific and non-specific protein association in solution: computation of solvent effects and prediction of first-encounter modes for efficient configurational bias Monte Carlo simulations

Weak and ultraweak protein-protein association play a role in molecular recognition and can drive spontaneous self-assembly and aggregation. Such interactions are difficult to detect experimentally, and are a challenge to the force field and sampling technique. A method is proposed to identify low-population protein-protein binding modes in aqueous solution. The method is designed to identify preferential first-encounter complexes from which the final complex(es) at equilibrium evolve. A continuum model is used to represent the effects of the solvent, which accounts for short- and long-range effects of water exclusion and for liquid-structure forces at protein/liquid interfaces. These effects control the behavior of proteins in close proximity and are optimized on the basis of binding enthalpy data and simulations. An algorithm is described to construct a biasing function for self-adaptive configurational-bias Monte Carlo of a set of interacting proteins. The function allows mixing large and local changes in the spatial distribution of proteins, thereby enhancing sampling of relevant microstates. The method is applied to three binary systems. Generalization to multiprotein complexes is discussed.

Using D-Amino Acids to Delineate the Mechanism of Protein Folding: Application to Trp-cage

Using the miniprotein Trp-cage as a model, we show that D-amino acids can be used to facilitate the delineation of protein folding mechanism. Specifically, we study the folding-unfolding kinetics of three Trp-cage mutants where the native glycine residue near the C-terminus of the α-helix is replaced by a D-amino acid. A previous study showed that these mutations increase the Trp-cage stability, due to a terminal capping effect. Our results show that the stabilizing effect of D-asparagine and D-glutamine originates almost exclusively from a decrease in the unfolding rate, while the D-alanine mutation results in a similar decrease in the unfolding rate, but it also increases the folding rate. Together, these results support a folding mechanism wherein the α-helix formation in the transition state is nucleated at the N-terminus, whereas those long-range native interactions stabilizing this helix are developed at the downhill side of the folding free energy barrier.

On the foldability of tryptophan-containing tetra- and pentapeptides: an exhaustive molecular dynamics study

Short peptides serve as minimal model systems to decipher the determinants of foldability due to their simplicity arising from their smaller size, their ability to echo protein-like structural characteristics, and their direct implication in force field validation. Here, we describe an effort to identify small peptides that can still form stable structures in aqueous solutions. We followed the in silico folding of a selected set of 8640 tryptophan-containing tetra- and pentapeptides through 15 210 molecular dynamics simulations amounting to a total of 272.46 μs using explicit representation of the solute and full treatment of the electrostatics. The evaluation and sorting of peptides is achieved through scoring functions, which include terms based on interatomic vector distances, atomic fluctuations, and rmsd matrices between successive frames of a trajectory. Highly scored peptides are studied further via successive simulation rounds of increasing simulation length and using different empirical force fields. Our method suggested only a handful of peptides with strong foldability prognosis. The discrepancies between the predictions of the various force fields for such short sequences are also extensively discussed. We conclude that the vast majority of such short peptides do not adopt significantly stable structures in water solutions, at least based on our computational predictions. The present work can be utilized in the rational design and engineering of bioactive peptides with desired molecular properties.

Water-exclusion and liquid-structure forces in implicit solvation

A continuum model of solvation is proposed to describe (i) long-range electrostatic effects of water exclusion resulting from incomplete and anisotropic hydration in crowded environments and (ii) short-range effects of liquid-structure forces on the hydrogen-bond interactions at solute/water interfaces. The model is an extension of the phenomenological screened coulomb potential-based implicit model of solvation. The developments reported here allow a more realistic representation of highly crowded and spatially heterogeneous environments, such as those in the interior of a living cell. Only the solvent is treated as a continuum medium. It is shown that the electrostatic effects of long-range water-exclusion can strongly affect protein-protein binding energies and are then related to the thermodynamics of complex formation. Hydrogen-bond interactions modulated by the liquid structure at interfaces are calibrated based on systematic calculations of potentials of mean force in explicit water. The electrostatic component of the model is parametrized for monovalent, divalent and trivalent ions. The conceptual and practical aspects of the model are discussed based on simulations of protein complexation and peptide folding. The current implementation is ~1.5 times slower than the gas-phase force field and exhibits good parallel performance.

Computational design of thermostabilizing D-amino acid substitutions

Judicious incorporation of D-amino acids in engineered proteins confers many advantages such as preventing degradation by endogenous proteases and promoting novel structures and functions not accessible to homochiral polypeptides. Glycine to D-alanine substitutions at the carboxy termini can stabilize α-helices by reducing conformational entropy. Beyond alanine, we propose additional side chain effects on the degree of stabilization conferred by D-amino acid substitutions. A detailed, molecular understanding of backbone and side chain interactions is important for developing rational, broadly applicable strategies in using D-amino acids to increase protein thermostability. Insight from structural bioinformatics combined with computational protein design can successfully guide the selection of stabilizing D-amino acid mutations. Substituting a key glycine in the Trp-cage miniprotein with D-Gln dramatically stabilizes the fold without altering the protein backbone. Stabilities of individual substitutions can be understood in terms of the balance of intramolecular forces both at the α-helix C-terminus and throughout the protein.

Molecular dynamics modeling of tubulin C-terminal tail interactions with the microtubule surface

Tubulin, an α/β heterodimer, has had most of its 3D structure analyzed; however, the carboxy (C)-termini remain elusive. Importantly, the C-termini play critical roles in regulating microtubule structure and function. They are sites of most of the post-translational modifications of tubulin and interaction sites with molecular motors and microtubule-associated proteins. Simulated annealing was used in our molecular dynamics modeling to predict the interactions of the C-terminal tails with the tubulin dimer. We examined differences in their flexibility, interactions with the body of tubulin, and the existence of structural motifs. We found that the α-tubulin tail interacts with the H11 helix of β-tubulin, and the β-tubulin tail interacts with the H11 helix of α-tubulin. Tail domains and H10/B9 loops interact with each other and compete for interactions with positively-charged residues of the H11 helix on the neighboring monomer. In a simulation in which α-tubulin's H10/B9 loop switches on sub-nanosecond intervals between interactions with the C-terminal tail of α-tubulin and the H11 helix of β-tubulin, the intermediate domain of α-tubulin showed more fluctuations compared to those in the other simulations, indicating that tail domains may cause shifts in the position of this domain. This suggests that C-termini may affect the conformation of the tubulin dimer which may explain their essential function in microtubule formation and effects on ligand binding to microtubules. Our modeling also provides evidence for a disordered-helical/helical double-state system of the T3/H3 region of the microtubule, which could be linked to depolymerization following GTP hydrolysis.

Comparative study of the energetics of ion permeation in Kv1.2 and KcsA potassium channels

Biological ion channels rely on a multi-ion transport mechanism for fast yet selective permeation of ions. The crystal structure of the KcsA potassium channel provided the first microscopic picture of this process. A similar mechanism is assumed to operate in all potassium channels, but the validity of this assumption has not been well investigated. Here, we examine the energetics of ion permeation in Shaker Kv1.2 and KcsA channels, which exemplify the six-transmembrane voltage-gated and two-transmembrane inward-rectifier channels. We study the feasibility of binding a third ion to the filter and the concerted motion of ions in the channel by constructing the potential of mean force for K(+) ions in various configurations. For both channels, we find that a pair of K(+) ions can move almost freely within the filter, but a relatively large free-energy barrier hinders the K(+) ion from stepping outside the filter. We discuss the effect of the CMAP dihedral energy correction that was recently incorporated into the CHARMM force field on ion permeation dynamics.

Kinetics and mechanism of the unfolding native-to-loop transition of Trp-cage in explicit solvent via optimized forward flux sampling simulations

The native-to-loop (N-L) unfolding transition of Trp-cage protein was studied via optimized forward flux sampling (FFS) methods with trajectories evolved using molecular dynamics. The rate constant calculated from our simulations is in good agreement with the experimental value for the native-to-unfolded transition of this protein; furthermore, the trajectories sampled a phase region consistent with that reported in previous studies for the N-L transition using transition path sampling and transition interface sampling. A new variant of FFS is proposed and implemented that allows a better control of a constant flux of partial paths. A reaction coordinate model was obtained, at no extra cost, from the transition path ensemble generated by FFS, through iterative use of the FFS-least-square estimation method [E. E. Borrero and F. A. Escobedo, J. Chem. Phys. 127, 164101 (2007)] and an adaptive staging optimization algorithm [E. E. Borrero and F. A. Escobedo, J. Chem. Phys. 129, 024115 (2008)]. Finally, we further elucidate the unfolding mechanism by correlating the unfolding progress with changes in the root mean square deviation from the α carbons of the native state, the root mean square deviation from an ideal α-helix, and other structural properties of the protein.

Conformational determinants of phosphotyrosine peptides complexed with the Src SH2 domain

The inhibition of specific SH2 domain mediated protein-protein interactions as an effective chemotherapeutic approach in the treatment of diseases remains a challenge. That different conformations of peptide-ligands are preferred by different SH2 domains is an underappreciated observation from the structural analysis of phosphotyrosine peptide binding to SH2 domains that may aid in future drug design. To explore the nature of ligand binding, we use simulated annealing (SA) to sample the conformational space of phosphotyrosine-containing peptides complexed with the Src SH2 domain. While in good agreement with the crystallographic and NMR studies of high-affinity phosphopeptide-SH2 domain complexes, the results suggest that the structural basis for phopsphopeptide- Src SH2 interactions is more complex than the “two-pronged plug two-hole socket” model. A systematic study of peptides of type pYEEX, where pY is phosphotyrosine and X is a hydrophobic residue, indicates that these peptides can assume two conformations, one extended and one helical, representing the balance between the interaction of residue X with the hydrophobic hole on the surface of the Src SH2 domain, and its contribution to the inherent tendency of the two glutamic acids to form an α-helix. In contrast, a β-turn conformation, almost identical to that observed in the crystal structure of pYVNV bound to the Grb2 SH2 domain, predominates for pYXNX peptides, even in the presence of isoleucine at the third position. While peptide binding affinities, as measured by fluorescence polarization, correlate with the relative proportion of extended peptide conformation, these results suggest a model where all three residues C-terminal to the phosphotyrosine determine the conformation of the bound phosphopeptide. The information obtained in this work can be used in the design of specific SH2 domain inhibitors.

Importance of the peptide backbone description in modeling the selectivity filter in potassium channels

A dihedral energy correction (CMAP) term has been recently included in the CHARMM force field to obtain a more accurate description of the peptide backbone. Its importance in improving dynamical properties of proteins and preserving their stability in long molecular-dynamics simulations has been established for several globular proteins. Here we investigate its role in maintaining the structure and function of two potassium channels, Shaker K(v)1.2 and KcsA, by performing molecular-dynamics simulations with and without the CMAP correction in otherwise identical systems. We show that without CMAP, it is not possible to maintain the experimentally observed orientations of the carbonyl groups in the selectivity filter in Shaker, and the channel loses its selectivity property. In the case of KcsA, the channel retains some selectivity even without CMAP because the carbonyl orientations are relatively better preserved compared to Shaker.

CHARMM: the biomolecular simulation program

CHARMM (Chemistry at HARvard Molecular Mechanics) is a highly versatile and widely used molecular simulation program. It has been developed over the last three decades with a primary focus on molecules of biological interest, including proteins, peptides, lipids, nucleic acids, carbohydrates, and small molecule ligands, as they occur in solution, crystals, and membrane environments. For the study of such systems, the program provides a large suite of computational tools that include numerous conformational and path sampling methods, free energy estimators, molecular minimization, dynamics, and analysis techniques, and model-building capabilities. The CHARMM program is applicable to problems involving a much broader class of many-particle systems. Calculations with CHARMM can be performed using a number of different energy functions and models, from mixed quantum mechanical-molecular mechanical force fields, to all-atom classical potential energy functions with explicit solvent and various boundary conditions, to implicit solvent and membrane models. The program has been ported to numerous platforms in both serial and parallel architectures. This article provides an overview of the program as it exists today with an emphasis on developments since the publication of the original CHARMM article in 1983.

In silico molecular engineering for a targeted replacement in a tumor-homing peptide

A new amino acid has been designed as a replacement for arginine (Arg, R) to protect the tumor-homing pentapeptide CREKA (Cys-Arg-Glu-Lys-Ala) from proteases. This amino acid, denoted (Pro)hArg, is characterized by a proline skeleton bearing a specifically oriented guanidinium side chain. This residue combines the ability of Pro to induce turn-like conformations with the Arg side-chain functionality. The conformational profile of the CREKA analogue incorporating this Arg substitute has been investigated by a combination of simulated annealing and molecular dynamics. Comparison of the results with those previously obtained for the natural CREKA shows that (Pro)hArg significantly reduces the conformational flexibility of the peptide. Although some changes are observed in the backbone...backbone and side-chain...side-chain interactions, the modified peptide exhibits a strong tendency to accommodate turn conformations centered at the (Pro)hArg residue and the overall shape of the molecule in the lowest energy conformations characterized for the natural and the modified peptides exhibit a high degree of similarity. In particular, the turn orients the backbone such that the Arg, Glu, and Lys side chains face the same side of the molecule, which is considered important for bioactivity. These results suggest that replacement of Arg by (Pro)hArg in CREKA may be useful in providing resistance against proteolytic enzymes while retaining conformational features which are essential for tumor-homing activity.

ABSINTH: a new continuum solvation model for simulations of polypeptides in aqueous solutions

A new implicit solvation model for use in Monte Carlo simulations of polypeptides is introduced. The model is termed ABSINTH for self-Assembly of Biomolecules Studied by an Implicit, Novel, and Tunable Hamiltonian. It is designed primarily for simulating conformational equilibria and oligomerization reactions of intrinsically disordered proteins in aqueous solutions. The paradigm for ABSINTH is conceptually similar to the EEF1 model of Lazaridis and Karplus (Proteins 1999, 35, 133). In ABSINTH, the transfer of a polypeptide solute from the gas phase into a continuum solvent is the sum of a direct mean field interaction (DMFI), and a term to model the screening of polar interactions. Polypeptide solutes are decomposed into a set of distinct solvation groups. The DMFI is a sum of contributions from each of the solvation groups, which are analogs of model compounds. Continuum-mediated screening of electrostatic interactions is achieved using a framework similar to the one used for the DMFI. Promising results are shown for a set of test cases. These include the calculation of NMR coupling constants for short peptides, the assessment of the thermal stability of two small proteins, reversible folding of both an alpha-helix and a beta-hairpin forming peptide, and the polymeric properties of intrinsically disordered polyglutamine peptides of varying lengths. The tests reveal that the computational expense for simulations with the ABSINTH implicit solvation model increase by a factor that is in the range of 2.5-5.0 with respect to gas-phase calculations.

Replica exchange simulations of the thermodynamics of Abeta fibril growth

Replica exchange molecular dynamics and an all-atom implicit solvent model are used to probe the thermodynamics of deposition of Alzheimer's Abeta monomers on preformed amyloid fibrils. Consistent with the experiments, two deposition stages have been identified. The docking stage occurs over a wide temperature range, starting with the formation of the first peptide-fibril interactions at 500 K. Docking is completed when a peptide fully adsorbs on the fibril edge at the temperature of 380 K. The docking transition appears to be continuous, and occurs without free energy barriers or intermediates. During docking, incoming Abeta monomer adopts a disordered structure on the fibril edge. The locking stage occurs at the temperature of approximately 360 K and is characterized by the rugged free energy landscape. Locking takes place when incoming Abeta peptide forms a parallel beta-sheet structure on the fibril edge. Because the beta-sheets formed by locked Abeta peptides are typically off-registry, the structure of the locked phase differs from the structure of the fibril interior. The study also reports that binding affinities of two distinct fibril edges with respect to incoming Abeta peptides are different. The peptides bound to the concave edge have significantly lower free energy compared to those bound on the convex edge. Comparison with the available experimental data is discussed.

Ab initio folding of proteins with all-atom discrete molecular dynamics

Discrete molecular dynamics (DMD) is a rapid sampling method used in protein folding and aggregation studies. Until now, DMD was used to perform simulations of simplified protein models in conjunction with structure-based force fields. Here, we develop an all-atom protein model and a transferable force field featuring packing, solvation, and environment-dependent hydrogen bond interactions. We performed folding simulations of six small proteins (20-60 residues) with distinct native structures by the replica exchange method. In all cases, native or near-native states were reached in simulations. For three small proteins, multiple folding transitions are observed, and the computationally characterized thermodynamics are in qualitative agreement with experiments. The predictive power of all-atom DMD highlights the importance of environment-dependent hydrogen bond interactions in modeling protein folding. The developed approach can be used for accurate and rapid sampling of conformational spaces of proteins and protein-protein complexes and applied to protein engineering and design of protein-protein interactions.

The energy landscape of a selective tumor-homing pentapeptide

Recently, a potentially powerful strategy based on phage-display libraries has been presented to target tumors via homing peptides attached to nanoparticles. The Cys-Arg-Glu-Lys-Ala (CREKA) peptide sequence has been identified as a tumor-homing peptide that binds to clotted plasmas proteins present in tumor vessels and interstitium. The aim of this work consists of mapping the conformational profile of CREKA to identify the bioactive conformation. For this purpose, a conformational search procedure based on modified simulated annealing combined with molecular dynamics was applied to three systems that mimic the experimentally used conditions: (i) the free peptide; (ii) the peptide attached to a nanoparticle; and (iii) the peptide inserted in a phage display protein. In addition, the free peptide was simulated in an ionized aqueous solution environment, which mimics the ionic strength of the physiological medium. Accessible minima of all simulated systems reveal a multiple interaction pattern involving the ionized side chains of Arg, Glu, and Lys, which induces a beta-turn motif in the backbone observed in all simulated CREKA systems.

Folding pathways and rates for the three-stranded beta-sheet peptide Beta3s using discrete path sampling

The discrete path sampling method was used to investigate the folding of a three-stranded antiparallel beta-sheet peptide, Beta3s, described by an empirical potential and implicit solvent model. After application of a coarse-graining scheme that groups together sets of minima in local equilibrium, the calculated folding time was in reasonable agreement with other simulations and consistent with the experimental upper bound. The folding mechanism exhibited by the most significant discrete paths involves early formation of the C-terminal hairpin followed by docking of the N-terminal strand.

Temperature-induced dissociation of Abeta monomers from amyloid fibril

Using all-atom molecular dynamics, we study the temperature-induced dissociation of Abeta monomers from the fibril protofilament. To accelerate conformational sampling, simulations are performed at elevated temperatures and peptide concentrations. By computing free energy disconnectivity graphs we mapped the free energy landscape of monomers on the surface of Abeta fibril. We found that Abeta monomers sample diverse sets of low free energy states with different degrees of association with the fibril. Some of these states have residual amounts of fibril interactions, whereas others lack fibril-like content. Generally, Abeta monomers with partially formed fibril-like interactions have the lowest free energies, but their backbone conformations may differ considerably from those in the fibril interior. Overall, Abeta amyloid protofilaments seem to be highly resistant to thermal dissociation. Monomer dissociation from the fibril edge proceeds via multiple stages and pathways. Our simulation findings are discussed in the context of recent experimental results.

Computational study of the Trp-cage miniprotein based on the ECEPP/3 force field

Using a newly developed Monte Carlo global optimization method called basin paving, we have performed an ab initio computation for the structure of Trp-cage based on the ECEPP/3 force field in vacuo. The lowest energy minimum has been located. Its corresponding configuration is comparable to the native structure of Trp-cage (PDB code 1L2Y) with a backbone root mean square deviation of 2.24 A.

Replica exchange simulation of reversible folding/unfolding of the Trp-cage miniprotein in explicit solvent: On the structure and possible role of internal water

We simulate the folding/unfolding equilibrium of the 20-residue miniprotein Trp-cage. We use replica exchange molecular dynamics simulations of the AMBER94 atomic detail model of the protein explicitly solvated by water, starting from a completely unfolded configuration. We employ a total of 40 replicas, covering the temperature range between 280 and 538 K. Individual simulation lengths of 100 ns sum up to a total simulation time of about 4 micros. Without any bias, we observe the folding of the protein into the native state with an unfolding-transition temperature of about 440 K. The native state is characterized by a distribution of root mean square distances (RMSD) from the NMR data that peaks at 1.8A, and is as low as 0.4A. We show that equilibration times of about 40 ns are required to yield convergence. A folded configuration in the entire extended ensemble is found to have a lifetime of about 31 ns. In a clamp-like motion, the Trp-cage opens up during thermal denaturation. In line with fluorescence quenching experiments, the Trp-residue sidechain gets hydrated when the protein opens up, roughly doubling the number of water molecules in the first solvation shell. We find the helical propensity of the helical domain of Trp-cage rather well preserved even at very high temperatures. In the folded state, we can identify states with one and two buried internal water molecules interconnecting parts of the Trp-cage molecule by hydrogen bonds. The loss of hydrogen bonds of these buried water molecules in the folded state with increasing temperature is likely to destabilize the folded state at elevated temperatures.

Comparison of solvation-effect methods for the simulation of peptide interactions with a hydrophobic surface

In this study we investigated the interaction behavior between thirteen different small peptides and a hydrophobic surface using three progressively more complex methods of representing solvation effects: a united-atom implicit solvation method [CHARMM 19 force field (C19) with Analytical Continuum Electrostatics (ACE)], an all-atom implicit solvation method (C22 with GBMV), and an all-atom explicit solvation method (C22 with TIP3P). The adsorption behavior of each peptide was characterized by the calculation of the potential of mean force as a function of peptide-surface separation distance. The results from the C22/TIP3P model suggest that hydrophobic peptides exhibit relatively strong adsorption behavior, polar and positively-charged peptides exhibit negligible to relatively weak favorable interactions with the surface, and negatively-charged peptides strongly resist adsorption. Compared to the TIP3P model, the ACE and GBMV implicit solvent models predict much stronger attractions for the hydrophobic peptides as well as stronger repulsions for the negatively-charged peptides on the CH(3)-SAM surface. These comparisons provide a basis from which each of these implicit solvation methods may be reparameterized to provide closer agreement with explicitly represented solvation in simulations of peptide and protein adsorption to functionalized surfaces.

Balancing solvation and intramolecular interactions: toward a consistent generalized Born force field

The efficient and accurate characterization of solvent effects is a key element in the theoretical and computational study of biological problems. Implicit solvent models, particularly generalized Born (GB) continuum electrostatics, have emerged as an attractive tool to study the structure and dynamics of biomolecules in various environments. Despite recent advances in this methodology, there remain limitations in the parametrization of many of these models. In the present work, we demonstrate that it is possible to achieve a balanced implicit solvent force field by further optimizing the input atomic radii in combination with adjusting the protein backbone torsional energetics. This parameter optimization is guided by the potentials of mean force (PMFs) between amino acid polar groups, calculated from explicit solvent free energy simulations, and by conformational equilibria of short peptides, obtained from extensive folding and unfolding replica exchange molecular dynamics (REX-MD) simulations. Through the application of this protocol, the delicate balance between the competing solvation forces and intramolecular forces appears to be better captured, and correct conformational equilibria for a range of both helical and beta-hairpin peptides are obtained. The same optimized force field also successfully folds both beta-hairpin trpzip2 and mini-protein Trp-Cage, indicating that it is quite robust. Such a balanced, physics-based force field will be highly applicable to a range of biological problems including protein folding and protein structural dynamics.

Long dynamics simulations of proteins using atomistic force fields and a continuum representation of solvent effects: calculation of structural and dynamic properties

Long dynamics simulations were carried out on the B1 immunoglobulin-binding domain of streptococcal protein G (ProtG) and bovine pancreatic trypsin inhibitor (BPTI) using atomistic descriptions of the proteins and a continuum representation of solvent effects. To mimic frictional and random collision effects, Langevin dynamics (LD) were used. The main goal of the calculations was to explore the stability of tens-of-nanosecond trajectories as generated by this molecular mechanics approximation and to analyze in detail structural and dynamical properties. Conformational fluctuations, order parameters, cross correlation matrices, residue solvent accessibilities, pKa values of titratable groups, and hydrogen-bonding (HB) patterns were calculated from all of the trajectories and compared with available experimental data. The simulations comprised over 40 ns per trajectory for ProtG and over 30 ns per trajectory for BPTI. For comparison, explicit water molecular dynamics simulations (EW/MD) of 3 ns and 4 ns, respectively, were also carried out. Two continuum simulations were performed on each protein using the CHARMM program, one with the all-atom PAR22 representation of the protein force field (here referred to as PAR22/LD simulations) and the other with the modifications introduced by the recently developed CMAP potential (CMAP/LD simulations). The explicit solvent simulations were performed with PAR22 only. Solvent effects are described by a continuum model based on screened Coulomb potentials (SCP) reported earlier, i.e., the SCP-based implicit solvent model (SCP-ISM). For ProtG, both the PAR22/LD and the CMAP/LD 40-ns trajectories were stable, yielding C(alpha) root mean square deviations (RMSD) of about 1.0 and 0.8 A respectively along the entire simulation time, compared to 0.8 A for the EW/MD simulation. For BPTI, only the CMAP/LD trajectory was stable for the entire 30-ns simulation, with a C(alpha) RMSD of approximately 1.4 A, while the PAR22/LD trajectory became unstable early in the simulation, reaching a C(alpha) RMSD of about 2.7 A and remaining at this value until the end of the simulation; the C(alpha) RMSD of the EW/MD simulation was about 1.5 A. The source of the instabilities of the BPTI trajectories in the PAR22/LD simulations was explored by an analysis of the backbone torsion angles. To further validate the findings from this analysis of BPTI, a 35-ns SCP-ISM simulation of Ubiquitin (Ubq) was carried out. For this protein, the CMAP/LD simulation was stable for the entire simulation time (C(alpha) RMSD of approximately 1.0 A), while the PAR22/LD trajectory showed a trend similar to that in BPTI, reaching a C(alpha) RMSD of approximately 1.5 A at 7 ns. All the calculated properties were found to be in agreement with the corresponding experimental values, although local deviations were also observed. HB patterns were also well reproduced by all the continuum solvent simulations with the exception of solvent-exposed side chain-side chain (sc-sc) HB in ProtG, where several of the HB interactions observed in the crystal structure and in the EW/MD simulation were lost. The overall analysis reported in this work suggests that the combination of an atomistic representation of a protein with a CMAP/CHARMM force field and a continuum representation of solvent effects such as the SCP-ISM provides a good description of structural and dynamic properties obtained from long computer simulations. Although the SCP-ISM simulations (CMAP/LD) reported here were shown to be stable and the properties well reproduced, further refinement is needed to attain a level of accuracy suitable for more challenging biological applications, particularly the study of protein-protein interactions.

Polyproline and the "spectroscopic ruler" revisited with single-molecule fluorescence

To determine whether Forster resonance energy transfer (FRET) measurements can provide quantitative distance information in single-molecule fluorescence experiments on polypeptides, we measured FRET efficiency distributions for donor and acceptor dyes attached to the ends of freely diffusing polyproline molecules of various lengths. The observed mean FRET efficiencies agree with those determined from ensemble lifetime measurements but differ considerably from the values expected from Forster theory, with polyproline treated as a rigid rod. At donor-acceptor distances much less than the Forster radius R(0), the observed efficiencies are lower than predicted, whereas at distances comparable to and greater than R(0), they are much higher. Two possible contributions to the former are incomplete orientational averaging during the donor lifetime and, because of the large size of the dyes, breakdown of the point-dipole approximation assumed in Forster theory. End-to-end distance distributions and correlation times obtained from Langevin molecular dynamics simulations suggest that the differences for the longer polyproline peptides can be explained by chain bending, which considerably shortens the donor-acceptor distances.

Peptide and protein folding and conformational equilibria: theoretical treatment of electrostatics and hydrogen bonding with implicit solvent models

Since biomolecules exist in aqueous and membrane environments, the accurate modeling of solvation, and hydrogen bonding interactions in particular, is essential for the exploration of structure and function in theoretical and computational studies. In this chapter, we focus on alternatives to explicit solvent models and discuss recent advances in generalized Born (GB) implicit solvent theories. We present a brief review of the successes and shortcomings of the application of these theories to biomolecular problems that are strongly linked to backbone H-bonding and electrostatics. This discussion naturally leads us to explore existing areas for improvement in current GB theories and our approach towards addressing a number of the key issues that remain in the refinement of these models. Specifically, the critical importance of balancing solvation forces and intramolecular forces in GB models is illustrated by examining the influence of backbone hydrogen bond strength and backbone dihedral energetics on conformational equilibria of small peptids.