Computer Simulations of Peptides from the p53 DNA Binding Domain

We have studied the dynamics and thermodynamics of two of the four evolutionarily conserved segments from the p53 DNA binding domain, using molecular dynamics and replica exchange simulations. These two regions contain well-defined elements of secondary structure (a β hairpin for region II and an α helix for region V) and bind to DNA in the intact protein. They are also mutational hot spots. The goal of our study was to determine the stability and folding propensity of these peptides in isolation. We used three force fields and solvent models (CHARMM19 with EEF1, CHARMM27 with GBMV, GROMOS96 with SPC). The predicted stability, folding propensity, and secondary structures depend upon the potential. Secondary structure predictors identify helical propensity for region II, in agreement with one of the force fields (CHARMM/GBMV). However, the other two potentials favor β structure for this peptide, although the conformations may differ from the crystal. For region V secondary structure predictions are unclear. Only one force field (CHARMM/GBMV) produces low-lying free energy minima that retain some of the α helical structure from the crystal structure. The other two potentials appear to favor β structure for this peptide.

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.

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.

Conformational sampling with implicit solvent models: application to the PHF6 peptide in tau protein

Implicit solvent models approximate the effects of solvent through a potential of mean force and therefore make solvated simulations computationally efficient. Yet despite their computational efficiency, the inherent approximations made by implicit solvent models can sometimes lead to inaccurate results. To test the accuracy of a number of popular implicit solvent models, we determined whether implicit solvent simulations can reproduce the set of potential energy minima obtained from explicit solvent simulations. For these studies, we focus on a six-residue amino-acid sequence, referred to as the paired helical filament 6 (PHF6), which may play an important role in the formation of intracellular aggregates in patients with Alzheimer's disease. Several implicit solvent models form the basis of this work--two based on the generalized Born formalism, and one based on a Gaussian solvent-exclusion model. All three implicit solvent models generate minima that are in good agreement with minima obtained from simulations with explicit solvent. Moreover, free-energy profiles generated with each implicit solvent model agree with free-energy profiles obtained with explicit solvent. For the Gaussian solvent-exclusion model, we demonstrate that a straightforward ranking of the relative stability of each minimum suggests that the most stable structure is extended, a result in excellent agreement with the free-energy profiles. Overall, our data demonstrate that for some peptides like PHF6, implicit solvent can accurately reproduce the set of local energy minimum arising from quenched dynamics simulations with explicit solvent. More importantly, all solvent models predict that PHF6 forms extended beta-structures in solution, a finding consistent with the notion that PHF6 initiates neurofibrillary tangle formation in patients with Alzheimer's disease.

Thermal unfolding simulations of a multimeric protein--transition state and unfolding pathways

The folding of an oligomeric protein poses an extra challenge to the folding problem because the protein not only has to fold correctly; it has to avoid nonproductive aggregation. We have carried out over 100 molecular dynamics simulations using an implicit solvation model at different temperatures to study the unfolding of one of the smallest known tetramers, p53 tetramerization domain (p53tet). We found that unfolding started with disruption of the native tetrameric hydrophobic core. The transition state for the tetramer to dimer transition was characterized as a diverse ensemble of different structures using Phi value analysis in quantitative agreement with experimental data. Despite the diversity, the ensemble was still native-like with common features such as partially exposed tetramer hydrophobic core and shifts in the dimer-dimer arrangements. After passing the transition state, the secondary and tertiary structures continued to unfold until the primary dimers broke free. The free dimer had little secondary structure left and the final free monomers were random-coil like. Both the transition states and the unfolding pathways from these trajectories were very diverse, in agreement with the new view of protein folding. The multiple simulations showed that the folding of p53tet is a mixture of the framework and nucleation-condensation mechanisms and the folding is coupled to the complex formation. We have also calculated the entropy and effective energy for the different states along the unfolding pathway and found that the tetramerization is stabilized by hydrophobic interactions.

Triggering loops and enzyme function: identification of loops that trigger and modulate movements

Enzyme function often involves a conformational change. There is a general agreement that loops play a vital role in correctly positioning the catalytically important residues. Nevertheless, predicting the functional loops and most importantly their role in enzyme function remains a difficult task. A major reason for this difficulty is that loops that undergo conformational change are frequently not well conserved in their primary sequence. beta1,4-Galactosyltransferase is one such enzyme. There, the amino acid sequence of a long loop that undergoes a large conformational change upon substrate binding is not well conserved. Our molecular dynamics simulations show that the large conformational change in the long loop is brought about by a second, interacting loop. Interestingly, while the structural change of the second loop is much smaller than that of the long loop, its sequence (particularly glycine residues) is highly conserved. We further examine the generality of the proposition that there are loops that trigger movements but nevertheless show little or no structural changes in crystals. We focus on two other enzymes, enolase and lipase. We chose these enzymes, since they too undergo conformational change upon ligand binding, however, they have different folds and different functions. Through multiple sets of simulations we show that the conformational change of the functional loop(s) is brought about through communication of flexibility by triggering loops that have several glycine residues. We further propose that similar to the conservation of common favorable fold types and structural motifs, evolution has also conserved common "skillful" mechanisms. Mechanisms may be conserved across different folds, sequences and functions, with adaptation to specific enzymatic roles.

Effective energy function for proteins in lipid membranes

A simple extension of the EEF1 energy function to heterogeneous membrane-aqueous media is proposed. The extension consists of (a) development of solvation parameters for a nonpolar phase using experimental data for the transfer of amino acid side-chains from water to cyclohexane, (b) introduction of a heterogeneous membrane-aqueous system by making the reference solvation free energy of each atom dependent on the vertical coordinate, (c) a modification of the distance-dependent dielectric model to account for reduced screening of electrostatic interactions in the membrane, and (d) an adjustment of the EEF1 aqueous model in light of recent calculations of the potential of mean force between amino acid side-chains in water. The electrostatic model is adjusted to match experimental observations for polyalanine, polyleucine, and the glycophorin A dimer. The resulting energy function (IMM1) reproduces the preference of Trp and Tyr for the membrane interface, gives reasonable energies of insertion into or adsorption onto a membrane, and allows stable 1-ns MD simulations of the glycophorin A dimer. We find that the lowest-energy orientation of melittin in bilayers varies, depending on the thickness of the hydrocarbon layer.

Protein unfolding transitions in an intrinsically unstable annexin domain: molecular dynamics simulation and comparison with nuclear magnetic resonance data

Unfolding transitions of an intrinsically unstable annexin domain and the unfolded state structure have been examined using multiple approximately 10-ns molecular dynamics simulations. Three main basins are observed in the configurational space: native-like state, compact partially unfolded or intermediate compact state, and the unfolded state. In the native-like state fluctuations are observed that are nonproductive for unfolding. During these fluctuations, after an initial loss of approximately 20% of the core residue native contacts, the core of the protein transiently completely refolds to the native state. The transition from the native-like basin to the partially unfolded compact state involves approximately 75% loss of native contacts but little change in the radius of gyration or core hydration properties. The intermediate state adopts for part of the time in one of the trajectories a novel highly compact salt-bridge stabilized structure that can be identified as a conformational trap. The intermediate-to-unfolded state transition is characterized by a large increase in the radius of gyration. After an initial relaxation the unfolded state recovers a native-like topology of the domain. The simulated unfolded state ensemble reproduces in detail experimental nuclear magnetic resonance data and leads to a convincing complete picture of the unfolded domain.

Anti-cooperativity and cooperativity in hydrophobic interactions: Three-body free energy landscapes and comparison with implicit-solvent potential functions for proteins

Potentials of mean force (PMFs) of three-body hydrophobic association are investigated to gain insight into similar processes in protein folding. Free energy landscapes obtained from explicit simulations of three methanes in water are compared with that predicted by popular implicit-solvent effective potentials for the study of proteins. Explicit-water simulations show that for an extended range of three-methane configurations, hydrophobic association at 25 degrees C under atmospheric pressure is mostly anti-cooperative, that is, less favorable than if the interaction free energies were pairwise additive. Effects of free energy nonadditivity on the kinetic path of association and the temperature dependence of additivity are explored by using a three-methane system and simplified chain models. The prevalence of anti-cooperativity under ambient conditions suggests that driving forces other than hydrophobicity also play critical roles in protein thermodynamic cooperativity. We evaluate the effectiveness of several implicit-solvent potentials in mimicking explicit water simulated three-body PMFs. The favorability of the contact free energy minimum is found to be drastically overestimated by solvent accessible surface area (SASA). Both the SASA and a volume-based Gaussian solvent exclusion model fail to predict the desolvation barrier. However, this barrier is qualitatively captured by the molecular surface area model and a recent "hydrophobic force field." None of the implicit-solvent models tested are accurate for the entire range of three-methane configurations and several other thermodynamic signatures considered.

Conformational polymorphism of wild-type and mutant prion proteins: Energy landscape analysis

Conformational transitions are thought to be the prime mechanism of prion diseases. In this study, the energy landscapes of a wild-type prion protein (PrP) and the D178N and E200K mutant proteins were mapped, enabling the characterization of the normal isoforms (PrP(C)) and partially unfolded isoforms (PrP(PU)) of the three prion protein analogs. It was found that the three energy landscapes differ in three respects: (i) the relative stability of the PrP(C) and the PrP(PU) states, (ii) the transition pathways from PrP(C) to PrP(PU), and (iii) the relative stability of the three helices in the PrP(C) state. In particular, it was found that although helix 1 (residues 144-156) is the most stable helix in wild-type PrP, its stability is dramatically reduced by both mutations. This destabilization is due to changes in the charge distribution that affects the internal salt bridges responsible for the greater stability of this helix in wild-type PrP. Although both mutations result in similar destabilization of helix 1, they a have different effect on the overall stability of PrP(C) and of PrP(PU) isoforms and on structural properties. The destabilization of helix 1 by mutations provides additional evidences to the role of this helix in the pathogenic transition from the PrP(C) to the pathogenic isoform PrP(SC).

Contributions to the binding free energy of ligands to avidin and streptavidin

The free energy of binding of a ligand to a macromolecule is here formally decomposed into the (effective) energy of interaction, reorganization energy of the ligand and the macromolecule, conformational entropy change of the ligand and the macromolecule, and translational and rotational entropy loss of the ligand. Molecular dynamics simulations with implicit solvation are used to evaluate these contributions in the binding of biotin, biotin analogs, and two peptides to avidin and streptavidin. We find that the largest contribution opposing binding is the protein reorganization energy, which is calculated to be from 10 to 30 kcal/mol for the ligands considered here. The ligand reorganization energy is also significant for flexible ligands. The translational/rotational entropy is 4.5-6 kcal/mol at 1 M standard state and room temperature. The calculated binding free energies are in the correct range, but the large statistical uncertainty in the protein reorganization energy precludes precise predictions. For some complexes, the simulations show multiple binding modes, different from the one observed in the crystal structure. This finding is probably due to deficiencies in the force field but may also reflect considerable ligand flexibility.

Thermal unfolding molecular dynamics simulation of Escherichia coli dihydrofolate reductase: thermal stability of protein domains and unfolding pathway

Temperature induced unfolding of Escherichia coli dihydrofolate reductase was carried out by using molecular dynamic simulations. The simulations show that the unfolding generally involves an initial end-to-end collapse of the adenine binding domain into partially extended loops, followed by a gradual breakdown of the remaining beta sheet core structure. The core, which consists of beta strands 5-7, was observed to be the most resistant to thermal unfolding. This region, which is made up of part of the N terminus domain and part of the large domain of the E. coli dihydrofolate reductase, may constitute the nucleation site for protein folding and may be important for the eventual formation of both domains. The unfolding of different domains at different stages of the unfolding process suggests that protein domains vary in stability and that the rate at which they unfold can affect the overall outcome of the unfolding pathway. This observation is compared with the recently proposed hierarchical folding model. Finally, the results of the simulation were found to be consistent with a previous experimental study (Frieden, Proc Natl Acad Sci USA 1990;87:4413-4416) which showed that the folding process of E. coli dihydrofolate reductase involves sequential formation of the substrate binding sites.