Equilibrium folding and unfolding pathways for a model protein

The folding–unfolding process of reduced bovine pancreatic trypsin inhibitor was investigated with an idealized model employing approximate free energies. The protein is regarded to consist of only C^α and C^β atoms. The backbone dihedral angles are the only conformational variables and are permitted to take discrete values at every 10°. Intraresidue energies consist of two terms: an empirical part taken from the observed frequency distributions of (?,ψ) and an additional favorable energy assigned to the native conformation of each residue. Interresidue interactions are simplified by assuming that there is an attractive energy operative only between residue pairs in close contact in the native structure. A total of 230,000 molecular conformations, with no atomic overlaps, ranging from the native state to the denatured state, are randomly generated by changing the sampling bias. Each conformation is classified according to its conformational energy, F; a conformational entropy, S(F) is estimated for each value of F from the number of samples. The dependence of S(F) on energy reveals that the folding–unfolding transition for this idealized model is an “all-or-none” type; this is attributable to the specific long-range interactions. Interresidue contact probabilities, averaged over samples representing various stages of folding, serve to characterize folding intermediates. Most probable equilibrium pathways for the folding–unfolding transition are constructed by connecting conformationally similar intermediates. The specific details obtained for bovine pancreatic trypsin inhibitor are as follows: (1) Folding begins with the appearance of nativelike medium-range contacts at a β-turn and at the α-helix. (2) These grow to include the native pair of interacting β-strands. This state includes intact regular secondary conformations, as well as the interstrand sheet contacts, and corresponds to an activated state with the highest free energy on the pathway. (3) Additional native long-range contacts are completely formed either toward the amino terminus or toward the carboxyl terminus. (4) In a final step, the missing contacts appear. Although these folding pathways for this model are not consistent with experimental reports, it does indicate multiple folding pathways. The method is general and can be applied to any set of calculated conformational energies and furthermore permits investigation of gross folding features.     

Identifying folding nucleus based on residue contact networks of proteins

In the native structure of a protein, all the residues are tightly parked together in a specific order following its folding and every residue contacts with some spatially neighbor residues. A residue contact network can be constructed by defining the residues as nodes and the native contacts as edges. During the folding of small single-domain proteins, there is a set of contacts (or bonds), defined as the folding nucleus (FN), which is formed around the transition state, i.e., a rate-limiting barrier located at about the middle between the unfolded states and the native state on the free energy landscape. Such a FN plays an essential role in the folding dynamics and the residues, which form the related contacts called as folding nucleus residues (FNRs). In this work, the FNRs in proteins are identified by using quantities which characterize the topology of residue contact networks of proteins. By comparing the specificities of residues with the network quantities K(R), L(R), and D(R), up to 90% FNRs of six typical proteins found experimentally are identified. It is found that the FNRs behave the full-closeness centrals rather than degree or closeness centers in the residue contact network, implying that they are important to the folding cooperativity of proteins. Our study shows that the FNRs can be identified solely from the native structures of proteins based on the analysis of residue contact network without any knowledge of the transition state ensemble.     

An empirical energy function for threading protein sequence through the folding motif

In this paper we present a new residue contact potantial derived by statistical analysis of protein crystal structures. This gives mean hydrophobic and pairwise contact energies as a function of residue type and distance interval. To test the accuracy of this potential we generate model structures by “threading” different sequences through backbone folding motifs found in the structural data base. We find that conformational energies calculated by summing contact potentials show perfect specificity in matching the correct sequences with each globular folding motif in a 161-protcin data set. They also identify correct models with the core folding motifs of heme-rythrin and immunoglobulin McPC603 V₁-do- main, among millions of alternatives possible when we align subsequences with α-helices and β-strands, and allow for variation in the lengths of intervening loops. We suggest that contact potentials reflect important constraints on nonbonded interaction in native proteins, and that “threading” may be useful for structure prediction by recognition of folding motif. © 1993 Wiley-Liss, Inc.     

Non-native interactions play an effective role in protein folding dynamics

Systematic Monte Carlo simulations of simple lattice models show that the final stage of protein folding is an ordered process where native contacts get locked (i.e., the residues come into contact and remain in contact for the duration of the folding process) in a well‐defined order. The detailed study of the folding dynamics of protein‐like sequences designed as to exhibit different contact energy distributions, as well as different degrees of sequence optimization (i.e., participation of non‐native interactions in the folding process), reveals significant differences in the corresponding locking scenarios—the collection of native contacts and their average locking times, which are largely ascribable to the dynamics of non‐native contacts. Furthermore, strong evidence for a positive role played by non‐native contacts at an early folding stage was also found. Interestingly, for topologically simple target structures, a positive interplay between native and non‐native contacts is observed also toward the end of the folding process, suggesting that non‐native contacts may indeed affect the overall folding process. For target models exhibiting clear two‐state kinetics, the relation between the nucleation mechanism of folding and the locking scenario is investigated. Our results suggest that the stabilization of the folding transition state can be achieved through the establishment of a very small network of native contacts that are the first to lock during the folding process.     

Validity of Gō models: comparison with a solvent-shielded empirical energy decomposition

Do Gō-type model potentials provide a valid approach for studying protein folding? They have been widely used for this purpose because of their simplicity and the speed of simulations based on their use. The essential assumption in such models is that only contact interactions existing in the native state determine the energy surface of a polypeptide chain, even for non-native configurations sampled along folding trajectories. Here we use an all-atom molecular mechanics energy function to investigate the adequacy of Gō-type potentials. We show that, although the contact approximation is accurate, non-native contributions to the energy can be significant. The assumed relation between residue-residue interaction energies and the number of contacts between them is found to be only approximate. By contrast, individual residue energies correlate very well with the number of contacts. The results demonstrate that models based on the latter should give meaningful results (e.g., as used to interpret phi values), whereas those that depend on the former are only qualitative, at best.     

Structure-based design of model proteins

A structure-based, sequence-design procedure is proposed in which one considers a set of decoy structures that compete significantly with the target structure in being low energy conformations. The decoy structures are chosen to have strong overlaps in contacts with the putative native state. The procedure allows the design of sequences with large and small stability gaps in a random-bond heteropolymer model in both two and three dimensions by an appropriate assignment of the contact energies to both the native and nonnative contacts. The design procedure is also successfully applied to the two-dimensional HP model. Proteins 31:10–20, 1998. © 1998 Wiley-Liss, Inc.     

Comparing folding codes for proteins and polymers

Proteins fold to unique compact native structures. Perhaps other polymers could be designed to fold in similar ways. The chemical nature of the monomer “alphabet” determines the “energy matrix” of monomer interactions—which defines the folding code, the relationship between sequence and structure. We study two properties of energy matrices using two-dimensional lattice models: uniqueness, the number of sequences that fold to only one structure, and encodability, the number of folds that are unique lowest-energy structures of certain monomer sequences. For the simplest model folding code, involving binary sequences of H (hydrophobic) and P (polar) monomers, only a small fraction of sequences fold uniquely, and not all structures can be encoded. Adding strong repulsive interactions results in a folding code with more sequences folding uniquely and more designable folds. Some theories suggest that the quality of a folding code depends only on the number of letters in the monomer alphabet, but we find that the energy matrix itself can be at least as important as the size of the alphabet. Certain multi-letter codes, including some with 20 letters, may be less physical or protein-like than codes with smaller numbers of letters because they neglect correlations among inter-residue interactions, treat only maximally compact conformations, or add arbitrary energies to the energy matrix.     

Statistical mechanics of protein folding by exhaustive enumeration

It is hard to construct theories for the folding of globular proteins because they are large and complicated molecules having enormous numbers of nonnative conformations and having native states that are complicated to describe. Statistical mechanical theories of protein folding are constructed around major simplifying assumptions about the energy as a function of conformation and/or simplifications of the representation of the polypeptide chain, such as one point per residue on a cubic lattice. It is not clear how the results of these theories are affected by their various simplifications. Here we take a very different simplification approach where the chain is accurately represented and the energy of each conformation is calculated by a not unreasonable empirical function. However, the set of amino acid sequences and allowed conformations is so restricted that it becomes computationally feasible to examine them all. Hence we are able to calculate melting curves for thermal denaturation as well as the detailed kinetic pathway of refolding. Such calculations are based on a novel representation of the conformations as points in an abstract 12-dimensional Euclidean conformation space. Fast folding sequences have relatively high melting temperatures, native structures with relatively low energies, small kinetic barriers between local minima, and relatively many conformations in the global energy minimum's watershed. In contrast to other folding theories, these models show no necessary relationship between fast folding and an overall funnel shape to the energy surface, or a large energy gap between the native and the lowest nonnative structure, or the depth of the native energy minimum compared to the roughness of the energy landscape. Proteins 32:425–437, 1998. © 1998 Wiley-Liss, Inc.     

Folding protein alpha-carbon chains into compact forms by Monte Carlo methods.

A method is presented for generating folded chains of specific amino acid sequences on a simple cubic lattice. Monte Carlo simulations are used to transform extended geometries of simplified alpha-carbon chains for eight small monomeric globular proteins into folded states. Permitted chain transitions are limited to a few types of moves, all restricted to occur on the lattice. Crude residue-residue potentials derived from statistical structure data are used to describe the energies for each conformer. The low resolution structures obtained by this procedure contain many of the correct gross features of the native folded architectures with respect to average residue energy per nonbonded contact, segment density, and location of surface loops and disulfide pairs. Rms deviations between these and the native X-ray structures and percentage of native long-range contacts found in these final folded structures are 7.6 +/- 0.7 A and 48 +/- 3%, respectively. This procedure can be useful for predicting approximate tertiary interactions from amino acid sequence.     

Sequence determinants of protein folding rates: positive correlation between contact energy and contact range indicates selection for fast folding

In comparison with intense investigation of the structural determinants of protein folding rates, the sequence features favoring fast folding have received little attention. Here, we investigate this subject using simple models of protein folding and a statistical analysis of the Protein Data Bank (PDB). The mean-field model by Plotkin and coworkers predicts that the folding rate is accelerated by stronger-than-average interactions at short distance along the sequence. We confirmed this prediction using the Finkelstein model of protein folding, which accounts for realistic features of polymer entropy. We then tested this prediction on the PDB. We found that native interactions are strongest at contact range l = 8. However, since short range contacts tend to be exposed and they are frequently formed in misfolded structures, selection for folding stability tends to make them less attractive, that is, stability and kinetics may have contrasting requirements. Using a recently proposed model, we predicted the relationship between contact range and contact energy based on buriedness and contact frequency. Deviations from this prediction induce a positive correlation between contact range and contact energy, that is, short range contacts are stronger than expected, for 2/3 of the proteins. This correlation increases with the absolute contact order (ACO), as expected if proteins that tend to fold slowly due to large ACO are subject to stronger selection for sequence features favoring fast folding. Our results suggest that the selective pressure for fast folding is detectable only for one third of the proteins in the PDB, in particular those with large contact order.     

On the design and analysis of protein folding potentials

Pairwise interaction models to recognize native folds are designed and analyzed. Different sets of parameters are considered but the focus was on 20 x 20 contact matrices. Simultaneous solution of inequalities and minimization of the variance of the energy find matrices that recognize exactly the native folds of 572 sequences and structures from the protein data bank (PDB). The set includes many homologous pairs, which present a difficult recognition problem. Significant recognition ability is recovered with a small number of parameters (e.g., the H/P model). However, full recognition requires a complete set of amino acids. In addition to structures from the PDB, a folding program (MONSSTER) was used to generate decoy structures for 75 proteins. It is impossible to recognize all the native structures of the extended set by contact potentials. We therefore searched for a new functional form. An energy function U, which is based on a sum of general pairwise interactions limited to a resolution of 1 angstrom, is considered. This set was infeasible too. We therefore conjecture that it is not possible to find a folding potential, resolved to 1 angstrom, which is a sum of pair interactions.     

Long- and short-range interactions in native protein structures are consistent/minimally frustrated in sequence space

We show that long- and short-range interactions in almost all protein native structures are actually consistent with each other for coarse-grained energy scales; specifically we mean the long-range inter-residue contact energies and the short-range secondary structure energies based on peptide dihedral angles, which are potentials of mean force evaluated from residue distributions observed in protein native structures. This consistency is observed at equilibrium in sequence space rather than in conformational space. Statistical ensembles of sequences are generated by exchanging residues for each of 797 protein native structures with the Metropolis method. It is shown that adding the other category of interaction to either the short- or long-range interactions decreases the means and variances of those energies for essentially all protein native structures, indicating that both interactions consistently work by more-or-less restricting sequence spaces available to one of the interactions. In addition to this consistency, independence by these interaction classes is also indicated by the fact that there are almost no correlations between them when equilibrated using both interactions and significant but small, positive correlations at equilibrium using only one of the interactions. Evidence is provided that protein native sequences can be regarded approximately as samples from the statistical ensembles of sequences with these energy scales and that all proteins have the same effective conformational temperature. Designing protein structures and sequences to be consistent and minimally frustrated among the various interactions is a most effective way to increase protein stability and foldability.     

No evidence that mRNAs have lower folding free energies than random sequences with the same dinucleotide distribution

This work investigates whether mRNA has a lower estimated folding free energy than random sequences. The free energy estimates are calculated by the mfold program for prediction of RNA secondary structures. For a set of 46 mRNAs it is shown that the predicted free energy is not significantly different from random sequences with the same dinucleotide distribution. For random sequences with the same mononucleotide distribution it has previously been shown that the native mRNA sequences have a lower predicted free energy, which indicates a more stable structure than random sequences. However, dinucleotide content is important when assessing the significance of predicted free energy as the physical stability of RNA secondary structure is known to depend on dinucleotide base stacking energies. Even known RNA secondary structures, like tRNAs, can be shown to have predicted free energies indistinguishable from randomized sequences. This suggests that the predicted free energy is not always a good determinant for RNA folding.     

Predicted residue–residue contacts can help the scoring of 3D models

During the 7th Critical Assessment of Protein Structure Prediction (CASP7) experiment, it was suggested that the real value of predicted residue–residue contacts might lie in the scoring of 3D model structures. Here, we have carried out a detailed reassessment of the contact predictions made during the recent CASP8 experiment to determine whether predicted contacts might aid in the selection of close‐to‐native structures or be a useful tool for scoring 3D structural models. We used the contacts predicted by the CASP8 residue–residue contact prediction groups to select models for each target domain submitted to the experiment. We found that the information contained in the predicted residue–residue contacts would probably have helped in the selection of 3D models in the free modeling regime and over the harder comparative modeling targets. Indeed, in many cases, the models selected using just the predicted contacts had better GDT‐TS scores than all but the best 3D prediction groups. Despite the well‐known low accuracy of residue–residue contact predictions, it is clear that the predictive power of contacts can be useful in 3D model prediction strategies. Proteins 2010. © 2010 Wiley‐Liss, Inc.     

The Limited Role of Nonnative Contacts in the Folding Pathways of a Lattice Protein

Models of protein energetics that neglect interactions between amino acids that are not adjacent in the native state, such as the Gō model, encode or underlie many influential ideas on protein folding. Implicit in this simplification is a crucial assumption that has never been critically evaluated in a broad context: Detailed mechanisms of protein folding are not biased by nonnative contacts, typically argued to be a consequence of sequence design and/or topology. Here we present, using computer simulations of a well-studied lattice heteropolymer model, the first systematic test of this oft-assumed correspondence over the statistically significant range of hundreds of thousands of amino acid sequences that fold to the same native structure. Contrary to previous conjectures, we find a multiplicity of folding mechanisms, suggesting that Gō-like models cannot be justified by considerations of topology alone. Instead, we find that the crucial factor in discriminating among topological pathways is the heterogeneity of native contact energies: The order in which native contacts accumulate is profoundly insensitive to omission of nonnative interactions, provided that native contact heterogeneity is retained. This robustness holds over a surprisingly wide range of folding rates for our designed sequences. Mirroring predictions based on the principle of minimum frustration, fast-folding sequences match their Gō-like counterparts in both topological mechanism and transit times. Less optimized sequences dwell much longer in the unfolded state and/or off-pathway intermediates than do Gō-like models. For dynamics that bridge unfolded and unfolded states, however, even slow folders exhibit topological mechanisms and transit times nearly identical with those of their Gō-like counterparts. Our results do not imply a direct correspondence between folding trajectories of Gō-like models and those of real proteins, but they do help to clarify key topological and energetic assumptions that are commonly used to justify such caricatures.     

Structured pathway across the transition state for peptide folding revealed by molecular dynamics simulations

Small globular proteins and peptides commonly exhibit two-state folding kinetics in which the rate limiting step of folding is the surmounting of a single free energy barrier at the transition state (TS) separating the folded and the unfolded states. An intriguing question is whether the polypeptide chain reaches, and leaves, the TS by completely random fluctuations, or whether there is a directed, stepwise process. Here, the folding TS of a 15-residue β-hairpin peptide, Peptide 1, is characterized using independent 2.5 μs-long unbiased atomistic molecular dynamics (MD) simulations (a total of 15 μs). The trajectories were started from fully unfolded structures. Multiple (spontaneous) folding events to the NMR-derived conformation are observed, allowing both structural and dynamical characterization of the folding TS. A common loop-like topology is observed in all the TS structures with native end-to-end and turn contacts, while the central segments of the strands are not in contact. Non-native sidechain contacts are present in the TS between the only tryptophan (W11) and the turn region (P7-G9). Prior to the TS the turn is found to be already locked by the W11 sidechain, while the ends are apart. Once the ends have also come into contact, the TS is reached. Finally, along the reactive folding paths the cooperative loss of the W11 non-native contacts and the formation of the central inter-strand native contacts lead to the peptide rapidly proceeding from the TS to the native state. The present results indicate a directed stepwise process to folding the peptide.     

154 Protein folding and binding funnels: a common driving force and a common mechanism

Under the free energy landscape theory, both the protein-folding and protein–ligand binding processes are driven by the decrease in total Gibbs free energy of the protein-solvent or protein–ligand-solvent system, which involves the non-complementary changes between the entropy and enthalpy, ultimately leading to a global free energy minimization of these thermodynamic systems (Ji & Liu, 2011; Liu et al., 2012; Yang, Ji & Liu, 2012). In the case of protein folding, the lowering of the system free energy coupled with the gradual reduction in conformational degree of freedom of the folding intermediates determines that the shape of the free energy landscape for protein folding must be funnel-like (Dill & Chan, 1997), rather than non-funneled shapes (Ben-Naim, 2012). In the funnel-like free energy landscape, protein folding can be viewed as going down the hill via multiple parallel routes from a vast majority of individual non-native states on surface outside the funnel to the native states located around the bottom of the funnel. The first stage of folding, i.e. the rapid hydrophobic collapse process, is driven by the solvent entropy maximization. Concretely, the water molecules squeeze and sequestrate the hydrophobic amino acid side chains within the interior of the folding intermediates while exposing the polar and electrostatically charged side chains on the intermediate surface so as to minimize the solvent-accessible surface area of the solute and thus, the minimal contacts between the folding intermediates and the water molecules. This will maximize the entropy of the solvent, thus contributing substantially to lowering of the system free energy due to an absolute advantage of the solvent in both quantity and mass (Yang, Ji & Liu, 2012). The resulting molten globule states (Ohgushi & Wada, 1983), within which a few transient secondary structural components and tertiary contacts have been formed but many native contacts or close residue–residue interactions has yet to form, need to be further sculptured into the native states. This is a relatively slow “bottleneck” process because the competitive interactions between protein residues within the folding intermediates and between residues and water molecules may repeat many rounds to accumulate a large enough number of stable noncovalent bonds capable of counteracting the conformational entropy loss of the intermediates, thus putting this bottleneck stage under the enthalpy control (i.e. negative enthalpy change), contributing further to the lowering of the system free energy. Although the protein–ligand association occurs around the rugged bottom of the free energy landscape, the exclusion of water from the binding interfaces and the formation of noncovalent bonds between the two partners can still lower the system free energy. In conjunction with the loss of the rotational and translational degrees of freedom of the two partners as well as the loss of the conformational entropy of the protein, these processes could merge, downwards expand, and further narrow the free energy wells within which the protein–ligand binding process takes place, thereby making them look like a funnel, which we term the binding funnel. In this funnel, the free energy downhill process follows a similar paradigm to the protein-folding process. For example, if the initial collisions/contacts occur between the properly complementary interfaces of the protein and ligand, a large amount of water molecules (which usually form a water network around the solute surface) will be displaced to suit the need for maximizing the solvent entropy. This process is similar to that of the hydrophobic collapse during protein folding, resulting in a loosely associated protein–ligand complex that needs also to be further adapted into a tight complex, i.e. the second step which is mainly driven by the negative enthalpy change through intermolecular competitive interactions to gradually accumulate the noncovalent bonds and ultimately, to stabilize the complex at a tightly bound state. Taken together, we conclude that whether in the protein-folding or in the protein–ligand binding process, both the entropy-driven first step and the enthalpy-driven second step contribute to the lowering of the system free energy, resulting in the funnel-like folding or binding free energy landscape.     

Use of a quantitative structure-property relationship to design larger model proteins that fold rapidly

A quantitative structure-property relationship (QSPR) was used to design model protein sequences that fold repeatedly and relatively rapidly to stable target structures. The specific model was a 125-residue heteropolymer chain subject to Monte Carlo dynamics on a simple cubic lattice. The QSPR was derived from an analysis of a database of 200 sequences by a statistical method that uses a genetic algorithm to select the sequence attributes that are most important for folding and a neural network to determine the corresponding functional dependence of folding ability on the chosen attributes. The QSPR depends on the number of anti-parallel sheet contacts, the energy gap between the native state and quasi-continuous part of the spectrum and the total energy of the contacts between surface residues. Two Monte Carlo procedures were used in series to optimize both the target structures and the sequences. We generated 20 fully optimized sequences and 60 partially optimized control sequences and tested each for its ability to fold in dynamic MC simulations. Although sequences in which either the number of anti-parallel sheet contacts or the energy of the surface residues is non-optimal are capable of folding almost as well as fully optimized ones, sequences in which only the energy gap is optimized fold markedly more slowly. Implications of the results for the design of proteins are discussed.     

mRNAs have greater negative folding free energies than shuffled or codon choice randomized sequences

An examination of 51 mRNA sequences in GenBank has revealed that calculated mRNA folding is more stable than expected by chance. Free energy minimization calculations of native mRNA sequences are more negative than randomized mRNA sequences with the same base composition and length. Randomization of the coding region of genes yields folding free energies of less negative magnitude than the original native mRNA sequence. Randomization of codon choice, while still preserving original base composition, also results in less stable mRNAs. This suggests that a bias in the selection of codons favors the potential formation of mRNA structures which contribute to folding stability.