Effective energy function for proteins in solution

Forces and energetics of hapten-antibody dissociation: a biased molecular dynamics simulation study

The unbinding of fluorescein from the single-chain Fv fragment of the 4D5Flu antibody is investigated by biased molecular dynamics with an implicit solvation model. To obtain statistically meaningful results, a large number of unbinding trajectories are calculated; they involve a total simulation time of more than 200 ns. Simulations are carried out with a time-dependent perturbation and in the presence of a constant force. The two techniques, which provide complementary information, induce unbinding by favoring an increase in the distance between the ligand and the antibody. This distance is an appropriate progress variable for the dissociation reaction and permits direct comparison of the unbinding forces in the simulations with data from atomic force microscopy (AFM). The time-dependent perturbation generates unfolding pathways that are close to equilibrium and can be used to reconstruct the mean force; i.e. the derivative of the potential of mean force, along the reaction coordinate. This is supported by an analysis of the overall unbinding profile and the magnitude of the mean force, which are similar to those of the unbinding force (i.e. the external force due to the time-dependent perturbation) averaged over several unbinding events. The multiple simulations show that unbinding proceeds along a rather well-defined pathway for a broad range of effective pulling speeds. Initially, there is a distortion of the protein localized in the C-terminal region followed by the fluorescein exit from the binding site. This occurs in steps that involve breaking of specific electrostatic and van der Waals interactions. It appears that the simulations do not explore the same barriers as those measured in the AFM experiments because of the much higher unfolding speed in the former. The dependence of the force on the logarithm of the loading rate is linear and the slope is higher than in the AFM, in agreement with experiment in other systems, where different slopes were observed for different regimes. Based on the unbinding events, mutations in the 4D5Flu antigen binding site are predicted to result in significant changes in the unbinding force.

Structure of the human 3alpha-hydroxysteroid dehydrogenase type 3 in complex with testosterone and NADP at 1.25-A resolution

The first crystallographic structure of human type 3 3alpha-hydroxysteroid dehydrogenase (3alpha-HSD3, AKR1C2), an enzyme playing a critical role in steroid hormone metabolism, has been determined in complex with testosterone and NADP at 1.25-A resolution. The enzyme's 17beta-HSD activity was studied in comparison with its 3alpha-HSD activity. The enzyme catalyzes the inactivation of dihydrotestosterone into 5alpha-androstane-3alpha,17beta-diol (3alpha-diol) as well as the transformation of androstenedione into testosterone. Using our homogeneous and highly active enzyme preparation, we have obtained 150-fold higher 3alpha-HSD specificity as compared with the former reports in the literature. Although the rat and the human 3alpha-HSDs share 81% sequence homology, our structure reveals significantly different geometries of the active sites. Substitution of the Ser(222) by a histidine in the human enzyme may compel the steroid to adopt a different binding to that previously described for the rat (Bennett, M. J., Albert, R. H., Jez, J. M., Ma, H., Penning, T. M., and Lewis, M. (1997) Structure 5, 799-T812). Furthermore, we showed that the affinity for the cofactor is higher in the human 3alpha-HSD3 than the rat enzyme due to the presence of additional hydrogen bonds on the adenine moiety and that the cofactor is present under its reduced form in the active site in our preparation.

Protein molecular dynamics with the generalized Born/ACE solvent model

Implicit solvent models are increasingly important for the study of proteins in aqueous solution. Here, the generalized Born (GB) solvent polarization model as implemented in the analytical ACE potential [Schaefer and Karplus (1996) J Phys Chem 100:1578] is used to perform molecular dynamics simulations of two small, homologous proteins: the immunoglobulin-binding domain of streptococcal protein G and the Ras binding domain of Raf. Several model parameterizations are compared through more than 60 ns of simulation. Results are compared with two simpler solvent models-an accessible surface area model and a distant-dependent dielectric model, with finite-difference Poisson calculations, with existing explicit solvent simulations, and with experimental data. The simpler models yield stable but distorted structures. The best GB/ACE implementation uses a set of atomic Voronoi volumes reported recently, obtained by averaging over a large database of crystallographic protein structures. A 20% reduction is applied to the volumes, compensating in an average sense for an excessive de-screening of individual charges inherent in the ACE self-energy and for an undersolvation of dipolar groups inherent in the GB screening function. This GB/ACE parameterization yields stable trajectories on the 0.5-1-ns time scale that deviate moderately (approximately 1.5-2.5 A) from the X-ray structure, reproduce approximately the surface distribution of charged, polar, and hydrophobic groups, and reproduce accurately backbone flexibility as measured by amide NMR-order parameters. Over longer time scales (1.5-3 ns), some of the protein G runs escape from the native energy basin and deviate strongly (3 A) from the native structure. The conformations sampled during the transition out of the native energy basin are overstabilized by the GB/ACE solvation model, as compared with a numerical treatment of the full dielectric continuum model.

An alternative derivation of the equations of motion in torsion space for a branched linear chain

Torsion space molecular dynamics may be more efficiently encoded if the global motions are separated from the internal motions. The equations of motion for single, non-cyclic chains are shown to be first order in the backbone angle parameters when the global frame of reference is ignored and second order otherwise. Adding a simple heuristic substitute for the global motions enables the encoding of dynamics for mixed constrained/unconstrained model systems.

Protein docking along smooth association pathways

We propose a docking method that mimics the way proteins bind. The method accounts for the dominant driving forces at the different length scales of the protein binding process, allowing for an efficient selection of a downhill path on the evolving receptor-ligand-free energy landscape. Starting from encounter complexes with as much as 10 A rms deviation from the native conformation, the method locally samples the six dimensional space of rigid-body receptor-ligand structures subject to a van der Waals constraint. The sampling is initially biased only by the desolvation and electrostatic components of the free energy, which capture the partial affinity of unbound structures that are more than 4 A away from the native state. Below this threshold, improved discrimination is attained by adding an increasing fraction of the van der Waals energy to the force field. The method, with no free parameters, was tested in eight different sets of independently crystallized receptor-ligand structures consistently predicting bound conformations with the lowest free energies and appropriate stability gap around 2 A from the native complex. This multistage approach is consistent with the underlying kinetics and internal structure of the free energy funnel to the bound state. Implications for the nature of the protein binding pathways are also discussed.

Structured disorder and conformational selection

Traditionally, molecular disorder has been viewed as local or global instability. Molecules or regions displaying disorder have been considered inherently unstructured. The term has been routinely applied to cases for which no atomic coordinates can be derived from crystallized molecules. Yet, even when it appears that the molecules are disordered, prevailing conformations exist, with population times higher than those of all alternate conformations. Disordered molecules are the outcome of rugged energy landscapes away from the native state around the bottom of the funnel. Ruggedness has a biological function, creating a distribution of structured conformers that bind via conformational selection, driving association and multimolecular complex formation, whether chain-linked in folding or unlinked in binding. We classify disordered molecules into two types. The first type possesses a hydrophobic core. Here, even if the native conformation is unstable, it still has a large enough population time, enabling its experimental detection. In the second type, no such hydrophobic core exists. Hence, the native conformations of molecules belonging to this category have shorter population times, hindering their experimental detection. Although there is a continuum of distribution of hydrophobic cores in proteins, an empirical, statistically based hydrophobicity function may be used as a guideline for distinguishing the two disordered molecule types. Furthermore, the two types relate to steps in the protein folding reaction. With respect to protein design, this leads us to propose that engineering-optimized specific electrostatic interactions to avoid electrostatic repulsion would reduce the type I disordered state, driving the molten globule (MG) --> native (N) state. In contrast, for overcoming the type II disordered state, in addition to specific interactions, a stronger hydrophobic core is also indicated, leading to the denatured --> MG --> N state.

Assessment of the quality of energy functions for protein folding by using a criterion derived with the help of the noisy go model

We propose a simple criterion based on the Z-scoreto assess the quality of energy functions for protein folding: one should obtain Z>-10 for the equilibrium ensembleat about native conditions. We derive this criterionby studying a Go model with random errors added to the native interactions. The dependence of the Z-score on the thermodynamic parameters,including the noise, can be precisely obtained in this case,as the ground state of the model is known exactly.We apply this criterion to rapidly rule out two otherwise promisingpairwise energy approximations.The advantage of adopting the present criterionis that it is not necessary to know the ground state of an energy function to assess its quality. It is sufficient to compute the Z-scorefrom a single equilibrium simulation at around the folding temperature.

Optimization of solvation models for predicting the structure of surface loops in proteins

A novel procedure for optimizing the atomic solvation parameters (ASPs) sigma(i) developed recently for cyclic peptides is extended to surface loops in proteins. The loop is free to move, whereas the protein template is held fixed in its X-ray structure. The energy is E(tot) = E(FF)(epsilon = nr) + summation operator sigma(i)A(i), where E(FF)(epsilon = nr) is the force-field energy of the loop-loop and loop-template interactions, epsilon = nr is a distance-dependent dielectric constant, and n is an additional parameter to be optimized. A(i) is the solvent-accessible surface area of atom i. The optimal sigma(i) and n are those for which the loop structure with the global minimum of E(tot)(n, sigma(i)) becomes the experimental X-ray structure. Thus, the ASPs depend on the force field and are optimized in the protein environment, unlike commonly used ASPs such as those of Wesson and Eisenberg (Protein Sci 1992;1:227-235). The latter are based on the free energy of transfer of small molecules from the gas phase to water and have been traditionally combined with various force fields without further calibration. We found that for loops the all-atom AMBER force field performed better than OPLS and CHARMM22. Two sets of ASPs [based on AMBER (n = 2)], optimized independently for loops 64-71 and 89-97 of ribonuclease A, were similar and thus enabled the definition of a best-fit set. All these ASPs were negative (hydrophilic), including those for carbon. Very good (i.e., small) root-mean-square-deviation values from the X-ray loop structure were obtained with the three sets of ASPs, suggesting that the best-fit set would be transferable to loops in other proteins as well. The structure of loop 13-24 is relatively stretched and was insensitive to the effect of the ASPs.

Protein folding and function: the N-terminal fragment in adenylate kinase

Three-dimensional protein folds range from simple to highly complex architectures. In complex folds, some building block fragments are more important for correct protein folding than others. Such fragments are typically buried in the protein core and mediate interactions between other fragments. Here we present an automated, surface area-based algorithm that is able to indicate which, among all local elements of the structure, is critical for the formation of the native fold, and apply it to structurally well-characterized proteins. In particular, we focus on adenylate kinase. The fragment containing the phosphate binding, P-loop (the "giant anion hole") flanked by a beta-strand and an alpha-helix near the N-terminus, is identified as a critical building block. This building block shows a high degree of sequence and structural conservation in all adenylate kinases. The results of our molecular dynamics simulations are consistent with this identification. In its absence, the protein flips to a stable, non-native state. In this misfolded conformation, the other local elements of the structure are in their native-like conformations; however, their association is non-native. Furthermore, this element is critically important for the function of the enzyme, coupling folding, and function.

Review: protein design--where we were, where we are, where we're going

Protein design has become a powerful approach for understanding the relationship between amino acid sequence and 3-dimensional structure. In the past 5 years, there have been many breakthroughs in the development of computational methods that allow the selection of novel sequences given the structure of a protein backbone. Successful design of protein scaffolds has now paved the way for new endeavors to design function. The ability to design sequences compatible with a fold may also be useful in structural and functional genomics by expanding the range of proteins used for fold recognition and for the identification of functionally important domains from multiple sequence alignments.

Solvation energy density occlusion approximation for evaluation of desolvation penalties in biomolecular interactions

In calculations involving many displacements of an interacting pair of biomolecules, such as brownian dynamics, the docking of a substrate/ligand to an enzyme/receptor, or the screening of a large number of ligands as prospective inhibitors for a particular receptor site, there is a need for rapid evaluation of the desolvation penalties of the interacting pair. Although continuum electrostatic treatments with distinct dielectric constants for solute and solvent provide an account of the electrostatics of solvation and desolvation, it is necessary to re-solve the Poisson equation, at considerable computational cost, for each displacement of the interacting pair. We present a new method that uses a formulation of continuum electrostatic solvation in terms of the solvation energy density and approximates desolvation in terms of the occlusion of this density. We call it the SEDO approximation. It avoids the need to re-solve the Poisson equation, as desolvation is now estimated by an integral over the occluded volume. Test calculations are presented for some simple model systems and for some real systems that have previously been studied using the Poisson equation approach: MHC class I protein-peptide complexes and a congeneric series of human immunodeficiency virus type 1 (HIV-1) protease--ligand complexes. For most of the systems considered, the trends and magnitudes of the desolvation component of interaction energies obtained using the SEDO approximation are in reasonable correlation with those obtained by re-solving the Poisson equation. In most cases, the error introduced by the SEDO approximation is much less than that of the often-used test-charge approximation for the charge-charge components of intermolecular interactions. Proteins 2001;43:12-27.

Native topology or specific interactions: what is more important for protein folding?

Fifty-five molecular dynamics runs of two three-stranded antiparallel beta-sheet peptides were performed to investigate the relative importance of amino acid sequence and native topology. The two peptides consist of 20 residues each and have a sequence identity of 15 %. One peptide has Gly-Ser (GS) at both turns, while the other has d-Pro-Gly ((D)PG). The simulations successfully reproduce the NMR solution conformations, irrespective of the starting structure. The large number of folding events sampled along the trajectories at 360 K (total simulation time of about 5 micros) yield a projection of the free-energy landscape onto two significant progress variables. The two peptides have compact denatured states, similar free-energy surfaces, and folding pathways that involve the formation of a beta-hairpin followed by consolidation of the unstructured strand. For the GS peptide, there are 33 folding events that start by the formation of the 2-3 beta-hairpin and 17 with first the 1-2 beta-hairpin. For the (D)PG peptide, the statistical predominance is opposite, 16 and 47 folding events start from the 2-3 beta-hairpin and the 1-2 beta-hairpin, respectively. These simulation results indicate that the overall shape of the free-energy surface is defined primarily by the native-state topology, in agreement with an ever-increasing amount of experimental and theoretical evidence, while the amino acid sequence determines the statistically predominant order of the events.

Solvent effects on the energy landscapes and folding kinetics of polyalanine

The effect of a solvation on the thermodynamics and kinetics of polyalanine (Ala(12)) is explored on the basis of its energy landscapes in vacuum and in an aqueous solution. Both energy landscapes are characterized by two basins, one associated with alpha-helical structures and the other with coil and beta-structures of the peptide. In both environments, the basin that corresponds to the alpha-helical structure is considerably narrower than the basin corresponding to the beta-state, reflecting their different contributions to the entropy of the peptide. In vacuum, the alpha-helical state of Ala(12) constitutes the native state, in agreement with common helical propensity scales, whereas in the aqueous medium, the alpha-helical state is destabilized, and the beta-state becomes the native state. Thus solvation has a dramatic effect on the energy landscape of this peptide, resulting in an inverted stability of the two states. Different folding and unfolding time scales for Ala(12) in hydrophilic and hydrophobic chemical environments are caused by the higher entropy of the native state in water relative to vacuum. The concept of a helical propensity has to be extended to incorporate environmental solvent effects.

Exploration of partially unfolded states of human alpha-lactalbumin by molecular dynamics simulation

Molecular dynamics simulations are used to probe the properties of non-native states of the protein human alpha-lactalbumin (human alpha-LA) with a detailed atomistic model in an implicit aqueous solvent environment. To sample the conformational space, a biasing force is introduced that increases the radius of gyration relative to the native state and generates a large number of low-energy conformers that differ in terms of their root-mean-square deviation, for a given radius of gyration. The resulting structures are relaxed by unbiased simulations and used as models of the molten globule and partly denatured states of human alpha-LA, based on measured radii of gyration obtained from nuclear magnetic resonance experiments. The ensembles of structures agree in their overall properties with experimental data available for the human alpha-LA molten globule and its more denatured states. In particular, the simulation results show that the native-like fold of the alpha-domain is preserved in the molten globule. Further, a considerable proportion of the antiparallel beta-strand in the beta-domain are present. This indicates that the lack of hydrogen exchange protection found experimentally for the beta-domain is due to rearrangement of the beta-sheet involving transient populations of non-native beta-structures. The simulations also provide details concerning the ensemble of structures that contribute as the molten globule unfolds and shows, in accord with experimental data, that unfolding is not cooperative; i.e. the various structural elements do not unfold simultaneously.

Molecular dynamics simulation of Escherichia coli dihydrofolate reductase and its protein fragments: relative stabilities in experiment and simulations

We have carried out molecular dynamics simulations of the native dihydrofolate reductase from Escherichia coli and several of its folded protein fragments at standard temperature. The simulations have shown fragments 1--36, 37--88, and 89--159 to be unstable, with a C(alpha)RMSD (C(alpha) root mean squared deviation) >5 A after 3.0 nsec of simulation. The unfolding of fragment 1--36 was immediate, whereas fragments 37--88 and 89--159 gradually unfolded because of the presence of the beta-sheet core structure. In the absence of residues 1--36, the two distinct domains comprising fragment 39--159 associated with each other, resulting in a stable conformation. This conformation retained most of its native structural elements. We have further simulated fragments derived from computational protein cutting. These were also found to be unstable, with the exception of fragment 104--159. In the absence of alpha(4), the loose loop region of residues 120--127 exhibited a beta-strand-like behavior, associating itself with the beta-sheet core of the protein fragment. The current study suggests that the folding of dihydrofolate reductase involves cooperative folding of distinct domains which otherwise would have been unstable as independent folded units in solution. Finally, the critical role of residues 1--36 in allowing the two distinct domains of fragment 104--159 to fold into the final native conformation is discussed.

On the transferability of atomic solvation parameters: Ab initio structural prediction of cyclic heptapeptides in DMSO

A statistical mechanics methodology for predicting the solution structures and populations of peptides developed recently is based on a novel method for optimizing implicit solvation models, which was applied initially to a cyclic hexapeptide in DMSO (C. Baysal and H. Meirovitch, Journal of American Chemical Society, 1998, vol. 120, pp. 800-812). Thus, the molecule has been described by the simplified energy function E(tot) = E(GRO) + summation operator(k) sigma(k)A(k), where E(GRO) is the GROMOS force-field energy, sigma(k) and A(k) are the atomic solvation parameter (ASP) and the solvent accessible surface area of atom k, respectively. In a more recent study, these ASPs have been found to be transferable to the cyclic pentapeptide cyclo(D-Pro(1)-Ala(2)-Ala(3)-Ala(4)-Ala(5)) in DMSO (C. Baysal and H. Meirovitch, Biopolymers, 2000, vol. 53, pp. 423-433). In the present paper, our methodology is applied to the cyclic heptapeptides axinastatin 2 [cyclo(Asn(1)-Pro(2)-Phe(3)-Val(4)-Leu(5)-Pro(6)-Val(7))] and axinastatin 3 [cyclo(Asn(1)-Pro(2)-Phe(3)-Ile(4)-Leu(5)-Pro(6)-Val(7))], in DMSO, which were studied by nmr by Mechnich et al. (Helvetica Chimica Acta, 1997, vol. 80, pp. 1338-1354). The calculations for axinastatin 2 show that special ASPs should be optimized for the partially charged side-chain atoms of Asn while the rest of the atoms take their values derived in our previous work; this suggests that similar optimization might be needed for other side chains as well. The solution structures of these peptides are obtained ab initio (i.e., without using experimental restraints) by an extensive conformational search based on E(GRO) alone and E(*)(tot), which consists of the new set of ASPs. For E(*)(tot), the theoretical values of proton-proton distances, (3)J coupling constants, and other properties are found to agree very well with the nmr results, and they are always better than those based on E(GRO).

Folding simulations of a three-stranded antiparallel beta -sheet peptide

Protein folding is a grand challenge of the postgenomic era. In this paper, 58 folding events sampled during 47 molecular dynamics trajectories for a total simulation time of more than 4 micros provide an atomic detail picture of the folding of a 20-residue synthetic peptide with a stable three-stranded antiparallel beta-sheet fold. The simulations successfully reproduce the NMR solution conformation, irrespective of the starting structure. The sampling of the conformational space is sufficient to determine the free energy surface and localize the minima and transition states. The statistically predominant folding pathway involves the formation of contacts between strands 2 and 3, starting with the side chains close to the turn, followed by association of the N-terminal strand onto the preformed 2-3 beta-hairpin. The folding mechanism presented here, formation of a beta-hairpin followed by consolidation, is in agreement with a computational study of the free energy surface of another synthetic three-stranded antiparallel beta-sheet by Bursulaya and Brooks [(1999) J. Am. Chem. Soc. 121, 9947-9951]. Hence, it might hold in general for antiparallel beta-sheets with short turns.

Native protein sequences are close to optimal for their structures

How large is the volume of sequence space that is compatible with a given protein structure? Starting from random sequences, low free energy sequences were generated for 108 protein backbone structures by using a Monte Carlo optimization procedure and a free energy function based primarily on Lennard-Jones packing interactions and the Lazaridis-Karplus implicit solvation model. Remarkably, in the designed sequences 51% of the core residues and 27% of all residues were identical to the amino acids in the corresponding positions in the native sequences. The lowest free energy sequences obtained for ensembles of native-like backbone structures were also similar to the native sequence. Furthermore, both the individual residue frequencies and the covariances between pairs of positions observed in the very large SH3 domain family were recapitulated in core sequences designed for SH3 domain structures. Taken together, these results suggest that the volume of sequence space optimal for a protein structure is surprisingly restricted to a region around the native sequence.

Modeling protein density of states: additive hydrophobic effects are insufficient for calorimetric two-state cooperativity

A well-established experimental criterion for two-state thermodynamic cooperativity in protein folding is that the van't Hoff enthalpy DeltaH(vH) around the transition midpoint is equal, or very nearly so, to the calorimetric enthalpy DeltaH(cal) of the entire transition. This condition is satisfied by many small proteins. We use simple lattice models to provide a statistical mechanical framework to elucidate how this calorimetric two-state picture may be reconciled with the hierarchical multistate scenario emerging from recent hydrogen exchange experiments. We investigate the feasibility of using inverse Laplace transforms to recover the underlying density of states (i.e., enthalpy distribution) from calorimetric data. We find that the constraint imposed by DeltaH(vH)/DeltaH(cal) approximately 1 on densities of states of proteins is often more stringent than other "two-state" criteria proposed in recent theoretical studies. In conjunction with reasonable assumptions, the calorimetric two-state condition implies a narrow distribution of denatured-state enthalpies relative to the overall enthalpy difference between the native and the denatured conformations. This requirement does not always correlate with simple definitions of "sharpness" of a transition and has important ramifications for theoretical modeling. We find that protein models that assume capillarity cooperativity can exhibit overall calorimetric two-state-like behaviors. However, common heteropolymer models based on additive hydrophobic-like interactions, including highly specific two-dimensional Gō models, fail to produce proteinlike DeltaH(vH)/DeltaH(cal) approximately 1. A simple model is constructed to illustrate a proposed scenario in which physically plausible local and nonlocal cooperative terms, which mimic helical cooperativity and environment-dependent hydrogen bonding strength, can lead to thermodynamic behaviors closer to experiment. Our results suggest that proteinlike thermodynamic cooperativity may require a cooperative interplay between local and nonlocal interactions. The prospect of using calorimetric data to constrain Z-scores of knowledge-based potentials is discussed.

Molecular dynamics estimates of ion diffusion in model hydrophobic and KcsA potassium channels

Molecular dynamics simulations are carried out to obtain estimates of diffusion coefficients of biologically important Na+, K+, Ca2+ and Cl- ions in hydrophobic cylindrical channels with varying radii and large reservoirs. Calculations for the cylindrical channels are compared to those for the KcsA potassium channel, for which the protein structure has recently been determined from X-ray diffraction experiments. Our results show that ion diffusion is maintained at reasonably high levels even within narrow channels, and does not support the very small diffusion coefficients used in some continuum models in order to fit experimental data. The present estimates of ion diffusion coefficients are useful in the calculation of channel conductance using the Poisson-Nernst-Planck theory, or Brownian dynamics.

Ab initio construction of protein tertiary structures using a hierarchical approach

We present a hierarchical method to predict protein tertiary structure models from sequence. We start with complete enumeration of conformations using a simple tetrahedral lattice model. We then build conformations with increasing detail, and at each step select a subset of conformations using empirical energy functions with increasing complexity. After enumeration on lattice, we select a subset of low energy conformations using a statistical residue-residue contact energy function, and generate all-atom models using predicted secondary structure. A combined knowledge-based atomic level energy function is then used to select subsets of the all-atom models. The final predictions are generated using a consensus distance geometry procedure. We test the feasibility of the procedure on a set of 12 small proteins covering a wide range of protein topologies. A rigorous double-blind test of our method was made under the auspices of the CASP3 experiment, where we did ab initio structure predictions for 12 proteins using this approach. The performance of our methodology at CASP3 is reasonably good and completely consistent with our initial tests.

Unfolding proteins by external forces and temperature: the importance of topology and energetics

Unfolding of proteins by forced stretching with atomic force microscopy or laser tweezer experiments complements more classical techniques using chemical denaturants or temperature. Forced unfolding is of particular interest for proteins that are under mechanical stress in their biological function. For beta-sandwich proteins (a fibronectin type III and an immunoglobulin domain), both of which appear in the muscle protein titin, the results of stretching simulations show important differences from temperature-induced unfolding, but there are common features that point to the existence of folding cores. Intermediates detected by comparing unfolding with a biasing perturbation and a constant pulling force are not evident in temperature-induced unfolding. For an alpha-helical domain (alpha-spectrin), which forms part of the cytoskeleton, there is little commonality in the pathways from unfolding induced by stretching and temperature. Comparison of the forced unfolding of the two beta-sandwich proteins and two alpha-helical proteins (the alpha-spectrin domain and an acyl-coenzyme A-binding protein) highlights important differences within and between protein classes that are related to the folding topologies and the relative stability of the various structural elements.

Computer simulations of protein folding by targeted molecular dynamics

We have performed 128 folding and 45 unfolding molecular dynamics runs of chymotrypsin inhibitor 2 (CI2) with an implicit solvation model for a total simulation time of 0.4 microseconds. Folding requires that the three-dimensional structure of the native state is known. It was simulated at 300 K by supplementing the force field with a harmonic restraint which acts on the root-mean-square deviation and allows to decrease the distance to the target conformation. High temperature and/or the harmonic restraint were used to induce unfolding. Of the 62 folding simulations started from random conformations, 31 reached the native structure, while the success rate was 83% for the 66 trajectories which began from conformations unfolded by high-temperature dynamics. A funnel-like energy landscape is observed for unfolding at 475 K, while the unfolding runs at 300 K and 375 K as well as most of the folding trajectories have an almost flat energy landscape for conformations with less than about 50% of native contacts formed. The sequence of events, i.e., secondary and tertiary structure formation, is similar in all folding and unfolding simulations, despite the diversity of the pathways. Previous unfolding simulations of CI2 performed with different force fields showed a similar sequence of events. These results suggest that the topology of the native state plays an important role in the folding process.

Ab initio prediction of the solution structures and populations of a cyclic pentapeptide in DMSO based on an implicit solvation model

Using a recently developed statistical mechanics methodology, the solution structures and populations of the cyclic pentapeptide cyclo(D-Pro(1)-Ala(2)-Ala(3)-Ala(4)-Ala(5)) in DMSO are obtained ab initio, i.e., without using experimental restraints. An important ingredient of this methodology is a novel optimization of implicit solvation parameters, which in our previous publication [Baysal, C.; Meirovitch, H. J Am Chem Soc 1998, 120, 800-812] has been applied to a cyclic hexapeptide in DMSO. The molecule has been described by the simplified energy function E(tot) = E(GRO) + summation operator(k) sigma(k)A(k), where E(GRO) is the GROMOS force-field energy, sigma(k) and A(k) are the atomic solvation parameter (ASP) and the solvent accessible surface area of atom k. This methodology, which relies on an extensive conformational search, Monte Carlo simulations, and free energy calculations, is applied here with E(tot) based on the ASPs derived in our previous work, and for comparison also with E(GRO) alone. For both models, entropy effects are found to be significant. For E(tot), the theoretical values of proton-proton distances and (3)J coupling constants agree very well with the NMR results [Mierke, D. F.; Kurz, M.; Kessler, H. J Am Chem Soc 1994, 116, 1042-1049], while the results for E(GRO) are significantly worse. This suggests that our ASPs might be transferrable to other cyclic peptides in DMSO as well, making our methodology a reliable tool for an ab initio structure prediction; obviously, if necessary, parts of this methodology can also be incorporated in a best-fit analysis where experimental restraints are used.

Leap-dynamics: efficient sampling of conformational space of proteins and peptides in solution

A molecular simulation scheme, called Leap-dynamics, that provides efficient sampling of protein conformational space in solution is presented. The scheme is a combined approach using a fast sampling method, imposing conformational 'leaps' to force the system over energy barriers, and molecular dynamics (MD) for refinement. The presence of solvent is approximated by a potential of mean force depending on the solvent accessible surface area. The method has been successfully applied to N-acetyl-L-alanine-N-methylamide (alanine dipeptide), sampling experimentally observed conformations inaccessible to MD alone under the chosen conditions. The method predicts correctly the increased partial flexibility of the mutant Y35G compared to native bovine pancreatic trypsin inhibitor. In particular, the improvement over MD consists of the detection of conformational flexibility that corresponds closely to slow motions identified by nuclear magnetic resonance techniques.

Effects of turn residues on beta-hairpin folding--a molecular dynamics study

Folding of beta-hairpin structures of synthetic peptides has been simulated using the molecular dynamics method with a solvent-referenced potential. Two similar sequences, Ac-MQIFVKS(D)PGKTITLKV-NH(2) and Ac-MQIFVKS(L)PGKTITLKV-NH(2), derived from the N-terminal beta-hairpin of ubiquitin, were used to study the effects of turn residues in beta-hairpin folding. The simulations were carried out for 80 ns at 297 K. With extended initial conformation, the (D)P-containing peptide folded into a stable 2:2 beta-hairpin conformation with a type II' beta-turn at (D)PG. The overall beta-hairpin ratio, calculated by the DSSP algorithm, was 32.6%. With randomly generated initial conformations, the peptide also formed the stable 2:2 beta-hairpin conformation. The interactions among the side chains in the 2:2 beta-hairpin were almost identical to those in the native protein. These interactions reduced the solvation energy upon folding and stabilized the beta-hairpin conformation. Without the solvent effect, the peptide did not fold into stable beta-hairpin structures. The solvent effect is crucial for the formation of the beta-hairpin conformation. The effect of the temperature has also been studied. The (L)P-containing peptide did not fold into a stable beta-hairpin conformation and had a much lower beta-hairpin ratio (16.6%). The( L)P-containing peptide has similar favorable side-chain interactions, but the turn formed by (L)PG does not connect well with the right-handed twist of the beta-strands. For comparison, the isolated N-terminal peptide of ubiquitin, Ac-MQIFVKTLTGKTITLEV-NH(2), was also simulated and its beta-hairpin ratio was low, indicating that the beta-hairpin in the native structure is stabilized by the interaction with the protein environment. These simulation results agreed qualitatively with the available experimental findings.

Understanding beta-hairpin formation

The kinetics of formation of protein structural motifs (e.g., alpha-helices and beta-hairpins) can provide information about the early events in protein folding. A recent study has used fluorescence measurements to monitor the folding thermodynamics and kinetics of a 16-residue beta-hairpin. In the present paper, we obtain the free energy surface and conformations involved in the folding of an atomistic model for the beta-hairpin from multicanonical Monte Carlo simulations. The results suggest that folding proceeds by a collapse that is downhill in free energy, followed by rearrangement to form a structure with part of the hydrophobic cluster; the hairpin hydrogen bonds propagate outwards in both directions from the partial cluster. Such a folding mechanism differs from the published interpretation of the experimental results, which is based on a helix-coil-type phenomenological model.

Discrimination of the native from misfolded protein models with an energy function including implicit solvation

An essential requirement for theoretical protein structure prediction is an energy function that can discriminate the native from non-native protein conformations. To date most of the energy functions used for this purpose have been extracted from a statistical analysis of the protein structure database, without explicit reference to the physical interactions responsible for protein stability. The use of the statistical functions has been supported by the widespread belief that they are superior for such discrimination to physics-based energy functions. An effective energy function which combined the CHARMM vacuum potential with a Gaussian model for the solvation free energy is tested for its ability to discriminate the native structure of a protein from misfolded conformations; the results are compared with those obtained with the vacuum CHARMM potential. The test is performed on several sets of misfolded structures prepared by others, including sets of about 650 good decoys for six proteins, as well as on misfolded structures of chymotrypsin inhibitor 2. The vacuum CHARMM potential is successful in most cases when energy minimized conformations are considered, but fails when applied to structures relaxed by molecular dynamics. With the effective energy function the native state is always more stable than grossly misfolded conformations both in energy minimized and molecular dynamics-relaxed structures. The present results suggest that molecular mechanics (physics-based) energy functions, complemented by a simple model for the solvation free energy, should be tested for use in the inverse folding problem, and supports their use in studies of the effective energy surface of proteins in solution. Moreover, the study suggests that the belief in the superiority of statistical functions for these purposes may be ill founded.