Robustness of atomistic Go models in predicting native-like folding intermediates.

Go models are exceedingly popular tools in computer simulations of protein folding. These models are native-centric, i.e., they are directly constructed from the protein's native structure. Therefore, it is important to understand up to which extent the atomistic details of the native structure dictate the folding behavior exhibited by Go models. Here we address this challenge by performing exhaustive discrete molecular dynamics simulations of a Go potential combined with a full atomistic protein representation. In particular, we investigate the robustness of this particular type of Go models in predicting the existence of intermediate states in protein folding. We focus on the N47G mutational form of the Spc-SH3 folding domain (x-ray structure) and compare its folding pathway with that of alternative native structures produced in silico. Our methodological strategy comprises equilibrium folding simulations, structural clustering, and principal component analysis.

Identification of a conserved aggregation-prone intermediate state in the folding pathways of Spc-SH3 amyloidogenic variants.

We compared the folding pathways of selected mutational variants of the α-spectrin SH3 domain (Spc-SH3) by using a continuum model that combines a full atomistic protein representation with the Go potential. Experimental data show that the N47G mutant shows very little tendency to aggregate while the N47A and triple mutant D48G(2Y) are both amyloidogenic, with the latter being clearly more aggregation prone. We identified a strikingly similar native-like folding intermediate across the three mutants, in which strand ß(1) is totally unstructured and more than half of the major hydrophobic core residues are highly solvent exposed. Results from extensive docking simulations show that the ability of the intermediates to dimerize is largely driven by strand ß(1) and is consistent with the in vitro aggregation behavior reported for the corresponding mutants. They further suggest that residues 44 and 53, which are key players in the nucleation-condensation mechanism of folding, are also important triggers of the aggregation process.

Identifying critical residues in protein folding: Insights from phi-value and P(fold) analysis.

We apply a simulational proxy of the phi-value analysis and perform extensive mutagenesis experiments to identify the nucleating residues in the folding "reactions" of two small lattice Go polymers with different native geometries. Our findings show that for the more complex native fold (i.e., the one that is rich in nonlocal, long-range bonds), mutation of the residues that form the folding nucleus leads to a considerably larger increase in the folding time than the corresponding mutations in the geometry that is predominantly local. These results are compared to data obtained from an accurate analysis based on the reaction coordinate folding probability P(fold) and on structural clustering methods. Our study reveals a complex picture of the transition state ensemble. For both protein models, the transition state ensemble is rather heterogeneous and splits up into structurally different populations. For the more complex geometry the identified subpopulations are actually structurally disjoint. For the less complex native geometry we found a broad transition state with microscopic heterogeneity. These findings suggest that the existence of multiple transition state structures may be linked to the geometric complexity of the native fold. For both geometries, the identification of the folding nucleus via the P(fold) analysis agrees with the identification of the folding nucleus carried out with the phi-value analysis. For the most complex geometry, however, the applied methodologies give more consistent results than for the more local geometry. The study of the transition state structure reveals that the nucleus residues are not necessarily fully native in the transition state. Indeed, it is only for the more complex geometry that two of the five critical residues show a considerably high probability of having all its native bonds formed in the transition state. Therefore, one concludes that, in general, the phi-value correlates with the acceleration/deceleration of folding induced by mutation, rather than with the degree of nativeness of the transition state, and that the "traditional" interpretation of phi-values may provide a more realistic picture of the structure of the transition state only for more complex native geometries.

Statistically enhanced promiscuity of structurally correlated patterns.

We predict that patterns with correlated surface density of atoms have statistically higher promiscuity (ability to bind stronger to an arbitrary pattern) as compared with noncorrelated patterns with the same average surface density. We suggest that this constitutes a generic design principle for highly connected proteins (hubs) in protein interaction networks. We develop an analytical theory for this effect. We show that our key predictions are generic and independent, qualitatively, on the specific form of the interatomic interaction potential, provided it has a finite range.

Structural similarity enhances interaction propensity of proteins.

We study statistical properties of interacting protein-like surfaces and predict two strong, related effects: (i) statistically enhanced self-attraction of proteins; (ii) statistically enhanced attraction of proteins with similar structures. The effects originate in the fact that the probability to find a pattern self-match between two identical, even randomly organized interacting protein surfaces is always higher compared with the probability for a pattern match between two different, promiscuous protein surfaces. This theoretical finding explains statistical prevalence of homodimers in protein-protein interaction networks reported earlier. Further, our findings are confirmed by the analysis of curated database of protein complexes that showed highly statistically significant overrepresentation of dimers formed by structurally similar proteins with highly divergent sequences ("superfamily heterodimers"). We suggest that promiscuous homodimeric interactions pose strong competitive interactions for heterodimers evolved from homodimers. Such evolutionary bottleneck is overcome using the negative design evolutionary pressure applied against promiscuous homodimer formation. This is achieved through the formation of highly specific contacts formed by charged residues as demonstrated both in model and real superfamily heterodimers.

Statistically enhanced self-attraction of random patterns.

In this work we develop a theory of interaction of randomly patterned surfaces as a generic prototype model of protein-protein interactions. The theory predicts that pairs of randomly superimposed identical (homodimeric) random patterns have always twice as large magnitude of the energy fluctuations with respect to their mutual orientation, as compared with pairs of different (heterodimeric) random patterns. The amplitude of the energy fluctuations is proportional to the square of the average pattern density, to the square of the amplitude of the potential and its characteristic length, and scales linearly with the area of surfaces. The greater dispersion of interaction energies in the ensemble of homodimers implies that strongly attractive complexes of random surfaces are much more likely to be homodimers, rather than heterodimers. Our findings suggest a plausible physical reason for the anomalously high fraction of homodimers observed in real protein interaction networks.

Constructing, verifying, and dissecting the folding transition state of chymotrypsin inhibitor 2 with all-atom simulations.

Experimentally, protein engineering and phi-value analysis is the method of choice to characterize the structure in folding transition state ensemble (TSE) of any protein. Combining experimental phi values and computer simulations has led to a deeper understanding of how proteins fold. In this report, we construct the TSE of chymotrypsin inhibitor 2 from published phi values. Importantly, we verify, by means of multiple independent simulations, that the conformations in the TSE have a probability of approximately 0.5 to reach the native state rapidly, so the TSE consists of true transition states. This finding validates the use of transition state theory underlying all phi-value analyses. Also, we present a method to dissect and study the TSE by generating conformations that have a disrupted alpha-helix (alpha-disrupted states) or disordered beta-strands 3 and 4 (beta-disrupted states). Surprisingly, the alpha-disrupted states have a stronger tendency to fold than the beta-disrupted states, despite the higher phi values for the alpha-helix in the TSE. We give a plausible explanation for this result and discuss its implications on protein folding and design. Our study shows that, by using both experiments and computer simulations, we can gain many insights into protein folding.

Understanding hierarchical protein evolution from first principles.

We propose a model that explains the hierarchical organization of proteins in fold families. The model, which is based on the evolutionary selection of proteins by their native state stability, reproduces patterns of amino acids conserved across protein families. Due to its dynamic nature, the model sheds light on the evolutionary time-scales. By studying the relaxation of the correlation function between consecutive mutations at a given position in proteins, we observe separation of the evolutionary time-scales: at short time intervals families of proteins with similar sequences and structures are formed, while at long time intervals the families of structurally similar proteins that have low sequence similarity are formed. We discuss the evolutionary implications of our model. We provide a "profile" solution to our model and find agreement between predicted patterns of conserved amino acids and those actually observed in nature.

Characterization of the folding kinetics of a three-helix bundle protein via a minimalist Langevin model.

We use a simple off-lattice Langevin model of protein folding to characterize the folding and unfolding of a fast-folding, 46 residue three-helix bundle. Under conditions at which the C-terminal helix is 30 % stable, we observe a clear three-state folding mechanism. In the on-pathway intermediate state, the middle and C-terminal helices are folded and in contact with each other, while the N-terminal region remains disordered. Nevertheless, under these conditions this intermediate is thermodynamically unstable relative to its unfolded state. The first and highest folding barrier corresponds to the organization of the hinge between the middle and C-terminal helices. A subsequent major barrier corresponds to the organization of the hinge between the middle and N-terminal helices. Hyperstabilizing the hinge regions leads to twice the folding rate that is obtained from hyperstabilizing the helices, even though much fewer contacts are involved in hinge hyperstabilization than in helix hyperstabilization. Unfolding follows single-exponential kinetics, even at temperatures only slightly above the folding transition temperature.

Excluded volume in protein side-chain packing.

The excluded volume occupied by protein side-chains and the requirement of high packing density in the protein interior should severely limit the number of side-chain conformations compatible with a given native backbone. To examine the relationship between side-chain geometry and side-chain packing, we use an all-atom Monte Carlo simulation to sample the large space of side-chain conformations. We study three models of excluded volume and use umbrella sampling to effectively explore the entire space. We find that while excluded volume constraints reduce the size of conformational space by many orders of magnitude, the number of allowed conformations is still large. An average repacked conformation has 20 % of its chi angles in a non-native state, a marked reduction from the expected 67 % in the absence of excluded volume. Interestingly, well-packed conformations with up to 50 % non-native chi angles exist. The repacked conformations have native packing density as measured by a standard Voronoi procedure. Entropy is distributed non-uniformly over positions, and we partially explain the observed distribution using rotamer probabilities derived from the Protein Data Bank database. In several cases, native rotamers that occur infrequently in the database are seen with high probability in our simulation, indicating that sequence-specific excluded volume interactions can stabilize rotamers that are rare for a given backbone. In spite of our finding that 65 % of the native rotamers and 85 % of chi(1) angles can be predicted correctly on the basis of excluded volume only, 95 % of positions can accommodate more than one rotamer in simulation. We estimate that, in order to quench the side-chain entropy observed in the presence of excluded volume interactions, other interactions (hydrophobic, polar, electrostatic) must provide an additional stabilization of at least 0.6 kT per residue in order to single out the native state.

Mode-coupling theory for heteropolymers.

We study the Langevin dynamics of a heteropolymer by means of a mode-coupling approximation scheme, giving rise to a set of coupled integro-differential equations relating the response and correlation functions. The analysis shows that there is a regime at low temperature characterized by out-of-equilibrium dynamics, with violation of time-translational invariance and of the fluctuation-dissipation theorem. The onset of aging dynamics at low temperatures gives insight into the nature of the slow dynamics of a disordered polymer. We also introduce a renormalization-group treatment of our mode-coupling equations, which supports our analysis, and might be applicable to other systems.

The folding thermodynamics and kinetics of crambin using an all-atom Monte Carlo simulation.

We present a novel Monte Carlo simulation of protein folding, in which all heavy atoms are represented as interacting hard spheres. This model includes all degrees of freedom relevant to folding, all side-chain and backbone torsions, and uses a Go potential. In this study, we focus on the 46 residue alpha/beta protein crambin and two of its structural components, the helix and helix hairpin. For a wide range of temperatures, we recorded multiple folding events of these three structures from random coils to native conformations that differ by less than 1 A C(alpha) dRMS from their crystal structure coordinates. The thermodynamics and kinetic mechanism of the helix-coil transition obtained from our simulation shows excellent agreement with currently available experimental and molecular dynamics data. Based on insights obtained from folding its smaller structural components, a possible folding mechanism for crambin is proposed. We observed that the folding occurs via a cooperative, first order-like process, and that many folding pathways to the native state exist. One particular sequence of events constitutes a "fast-folding" pathway where kinetic traps are avoided. At very low temperatures, a kinetic trap arising from the incorrect packing of side-chains was observed. These results demonstrate that folding to the native state can be observed in a reasonable amount of time on desktop computers even when an all-atom representation is used, provided the energetics sufficiently stabilize the native state.

Different circular permutations produced different folding nuclei in proteins: a computational study.

There have been many studies about the effect of circular permutation on the transition state/folding nucleus of proteins, with sometimes conflicting conclusions from different proteins and permutations. To clarify this important issue, we have studied two circular permutations of a lattice protein model with side-chains. Both permuted sequences have essentially the same native state as the original (wild-type) sequence. Circular permutant 1 cuts at the folding nucleus of the wild-type sequence. As a result, the permutant has a drastically different nucleus and folds more slowly than wild-type. In contrast, circular permutant 2 involves an incision at a site unstructured in the wild-type transition state, and the wild-type nucleus is largely retained in the permutant. In addition, permutant 2 displays both two-state and multi-state folding, with a native-like intermediate state occasionally populated. Neither the wild-type nor permutant 1 has a similar intermediate, and both fold in an apparently two-state manner. Surprisingly, permutant 2 folds at a rate identical with that of the wild-type. The intermediate in permutant 2 is stabilised by native and non-native interactions, and cannot be classified simply as on or off-pathway. So we advise caution in attributing experimental data to on or off-pathway intermediates. Finally, our work illuminates the results on alpha-spectrin SH3, chymotrypsin inhibitor 2 and beta-lactoglobulin, and supports a key assumption in the experimental efforts to locate potential nucleation sites of real proteins via circular permutations.

Energy profile of the space of model protein sequences.

A numerical study of the energy landscape of the space of model proteinsequences is carried out. As a consequence of the heterogeneity of thecontact energies among amino acids, the energy landscape displays a veryrough profile, a behaviour typical of frustrated systems. This givesraise to a hierarchical clustering of low-energy sequences and can have evolutionary consequences.

Proteins with selected sequences: a heteropolymeric study.

Protein sequences are expected not to be random but selected in order to form a stable native structure that is kinetically accessible. Therefore our model contains a selective temperature in sequence space (see [S. Ramanathan and E. Shakhnovich, Phys. Rev. E 50, 1303 (1994)] ) to optimize the sequence for the target conformation statistically. Replica calculations, which go beyond quadratic approximations in the field-theoretical Hamiltonian, are presented. A phase diagram indicating the temperatures and selective temperatures at which transitions to a frozen globule, i.e., the native state, occur is obtained. It is shown that going beyond the quadratic approximation in the field Hamiltonian is very important, since it results in a significant change of the phase diagram. Moreover, we suggest that a one-step replica permutation symmetry scheme is sufficient to solve the model. In addition to this we present a result for the sequence correlation function along the chain in the case of a short-ranged potential between the monomers. A correlation function between monomers that form a contact in the native state is given depending on the temperature and the interaction parameter.

Differential stabilization of two hydrophobic cores in the transition state of the villin 14T folding reaction.

We report the distribution of hydrophobic core contacts during the folding reaction transition state for villin 14T, a small 126-residue protein domain. The solution structure of villin 14T contains a central beta-sheet with two flanking hydrophobic cores; transition states for this protein topology have not been previously studied. Villin 14T has no disulfide bonds or cis-proline residues in its native state; it folds reversibly, and in an apparently two-state manner under some conditions. To map the hydrophobic core contacts in the transition state, 27 point mutations were generated at positions spread throughout the two hydrophobic cores. After each point mutation, comparison of the change in folding kinetics with the equilibrium destabilization indicates whether the site of mutation is stabilized in the transition state. The results show that the folding nucleus, or the sub-region with the strongest transition state contacts, is located in one of the two hydrophobic cores (the predominantly aliphatic core). The other hydrophobic core, which is mostly aromatic, makes much weaker contacts in the transition state. This work is the first transition state mapping for a protein with multiple major hydrophobic cores in a single folding unit; the hydrophobic cores cannot be separated into individual folding subdomains. The stabilization of only one hydrophobic core in the transition state illustrates that hydrophobic core formation is not intrinsically capable of nucleating folding, but must also involve the right specific interactions or topological factors in order to be kinetically important.

Statistical significance of protein structure prediction by threading.

In this study, we estimate the statistical significance of structure prediction by threading. We introduce a single parameter epsilon that serves as a universal measure determining the probability that the best alignment is indeed a native-like analog. Parameter epsilon takes into account both length and composition of the query sequence and the number of decoys in threading simulation. It can be computed directly from the query sequence and potential of interactions, eliminating the need for sequence reshuffling and realignment. Although our theoretical analysis is general, here we compare its predictions with the results of gapless threading. Finally we estimate the number of decoys from which the native structure can be found by existing potentials of interactions. We discuss how this analysis can be extended to determine the optimal gap penalties for any sequence-structure alignment (threading) method, thus optimizing it to maximum possible performance.

Comparison of two optimization methods to derive energy parameters for protein folding: perceptron and Z score.

Two methods were proposed recently to derive energy parameters from known native protein conformations and corresponding sets of decoys. One is based on finding, by means of a perceptron learning scheme, energy parameters such that the native conformations have lower energies than the decoys. The second method maximizes the difference between the native energy and the average energy of the decoys, measured in terms of the width of the decoys' energy distribution (Z-score). Whereas the perceptron method is sensitive mainly to "outlier" (i.e., extremal) decoys, the Z-score optimization is governed by the high density regions in decoy-space. We compare the two methods by deriving contact energies for two very different sets of decoys: the first obtained for model lattice proteins and the second by threading. We find that the potentials derived by the two methods are of similar quality and fairly closely related. This finding indicates that standard, naturally occurring sets of decoys are distributed in a way that yields robust energy parameters (that are quite insensitive to the particular method used to derive them). The main practical implication of this finding is that it is not necessary to fine-tune the potential search method to the particular set of decoys used.

What can disulfide bonds tell us about protein energetics, function and folding: simulations and bioninformatics analysis.

We study the impact of disulfide bonds on protein stability and folding. Using lattice model simulations, we show that formation of a disulfide bond stabilizes a protein to an extent that depends on the distance along the chain between linked cysteine residues. However, the impact of disulfide bonds on folding kinetics varies broadly, from acceleration when disulfides are introduced in or close to the folding nucleus, to slowing when disulfides are introduced outside the nucleus. Having established the effect of disulfide bonds on stability, we study the correlation between the number of disulfide bonds and the composition of certain amino acid classes with the goal to use it as a statistical probe into factors that contribute to stability of proteins. We find that the number of disulfides is negatively correlated with aliphatic hydrophobic but not aromatic content. It is surprising that we observe a strong correlation of disulfide content with polar (Q,S,T,N) amino acid content and a strong negative correlation with charged (E,D,K,R) content. These findings provide insights into factors that determine protein stability and principles of protein design as well as possible relations of disulfide bonds and protein function.

Analysis of knowledge-based protein-ligand potentials using a self-consistent method.

We propose a self-consistent approach to analyze knowledge-based atom-atom potentials used to calculate protein-ligand binding energies. Ligands complexed to actual protein structures were first built using the SMoG growth procedure (DeWitte & Shakhnovich, 1996) with a chosen input potential. These model protein-ligand complexes were used to construct databases from which knowledge-based protein-ligand potentials were derived. We then tested several different modifications to such potentials and evaluated their performance on their ability to reconstruct the input potential using the statistical information available from a database composed of model complexes. Our data indicate that the most significant improvement resulted from properly accounting for the following key issues when estimating the reference state: (1) the presence of significant nonenergetic effects that influence the contact frequencies and (2) the presence of correlations in contact patterns due to chemical structure. The most successful procedure was applied to derive an atom-atom potential for real protein-ligand complexes. Despite the simplicity of the model (pairwise contact potential with a single interaction distance), the derived binding free energies showed a statistically significant correlation (approximately 0.65) with experimental binding scores for a diverse set of complexes.