Free energy landscape for protein folding kinetics: Intermediates, traps, and multiple pathways in theory and lattice model simulations

Heterogeneity in Protein Folding and Unfolding Reactions

Proteins have dynamic structures that undergo chain motions on time scales spanning from picoseconds to seconds. Resolving the resultant conformational heterogeneity is essential for gaining accurate insight into fundamental mechanistic aspects of the protein folding reaction. The use of high-resolution structural probes, sensitive to population distributions, has begun to enable the resolution of site-specific conformational heterogeneity at different stages of the folding reaction. Different states populated during protein folding, including the unfolded state, collapsed intermediate states, and even the native state, are found to possess significant conformational heterogeneity. Heterogeneity in protein folding and unfolding reactions originates from the reduced cooperativity of various kinds of physicochemical interactions between various structural elements of a protein, and between a protein and solvent. Heterogeneity may arise because of functional or evolutionary constraints. Conformational substates within the unfolded state and the collapsed intermediates that exchange at rates slower than the subsequent folding steps give rise to heterogeneity on the protein folding pathways. Multiple folding pathways are likely to represent distinct sequences of structure formation. Insight into the nature of the energy barriers separating different conformational states populated during (un)folding can also be obtained by resolving heterogeneity.

Stochastic diffusion framework determines the free-energy landscape and rate from single-molecule trajectory

A theoretical stochastic diffusion framework is developed that characterizes the position-dependent diffusion coefficient [D(Q)] and drift velocity [ v (Q)] by analysing single-molecule time traces [Q(t)]. The free-energy landscape [F(Q)] that governs the dynamics is reconstructed with the calculated D and v . There are many computational tools that perform this task in which some are computationaly demanding, difficult to run, and, most of the time, not directly available to the community. This is a first attempt to implement the simplified stochastic diffusion framework that is fast, easy to run in a Python environment, and available to be extended as needed. It does not require adjustable parameters, inference methods, or sampling bias such as Monte Carlo Bayesian estimators or umbrella samplings. The stochastic framework was applied in the protein-like lattice model with Monte Carlo simulations, which accurately predicted the folding rates with the coordinate-dependent D and F plugged into Kramers' theory. The results were compared with two other independently developed methodologies (the Bayesian analysis and fep1D algorithm) presenting a good match, which confirms its validity. This theoretical framework might be useful in determining the free-energy and rates by providing time series only from biological or condensed-phase systems. The code is freely available at https://github.com/ronaldolab/stochastic_diffusion.

Accurate Protein-Folding Transition-Path Statistics from a Simple Free-Energy Landscape

A central goal of protein-folding theory is to predict the stochastic dynamics of transition paths-the rare trajectories that transit between the folded and unfolded ensembles-using only thermodynamic information, such as a low-dimensional equilibrium free-energy landscape. However, commonly used one-dimensional landscapes typically fall short of this aim, because an empirical coordinate-dependent diffusion coefficient has to be fit to transition-path trajectory data in order to reproduce the transition-path dynamics. We show that an alternative, first-principles free-energy landscape predicts transition-path statistics that agree well with simulations and single-molecule experiments without requiring dynamical data as an input. This "topological configuration" model assumes that distinct, native-like substructures assemble on a time scale that is slower than native-contact formation but faster than the folding of the entire protein. Using only equilibrium simulation data to determine the free energies of these coarse-grained intermediate states, we predict a broad distribution of transition-path transit times that agrees well with the transition-path durations observed in simulations. We further show that both the distribution of finite-time displacements on a one-dimensional order parameter and the ensemble of transition-path trajectories generated by the model are consistent with the simulated transition paths. These results indicate that a landscape based on transient folding intermediates, which are often hidden by one-dimensional projections, can form the basis of a predictive model of protein-folding transition-path dynamics.

Principles of Protein Stability and Their Application in Computational Design

Proteins are increasingly used in basic and applied biomedical research. Many proteins, however, are only marginally stable and can be expressed in limited amounts, thus hampering research and applications. Research has revealed the thermodynamic, cellular, and evolutionary principles and mechanisms that underlie marginal stability. With this growing understanding, computational stability design methods have advanced over the past two decades starting from methods that selectively addressed only some aspects of marginal stability. Current methods are more general and, by combining phylogenetic analysis with atomistic design, have shown drastic improvements in solubility, thermal stability, and aggregation resistance while maintaining the protein's primary molecular activity. Stability design is opening the way to rational engineering of improved enzymes, therapeutics, and vaccines and to the application of protein design methodology to large proteins and molecular activities that have proven challenging in the past.

Quantifying Nonnative Interactions in the Protein-Folding Free-Energy Landscape

Protein folding is a central problem in biological physics. Energetic roughness is an important aspect that controls protein-folding stability and kinetics. The roughness is associated with conflicting interactions in the protein and is also known as frustration. Recent studies indicate that an addition of a small amount of energetic frustration may enhance folding speed for certain proteins. In this study, we have investigated the conditions under which frustration increases the folding rate. We used a Cα structure-based model to simulate a group of proteins. We found that the free-energy barrier at the transition state (ΔF) correlates with nonnative-contact variation (ΔA), and the simulated proteins are clustered according to their fold motifs. These findings are corroborated by the Clementi-Plotkin analytical model. As a consequence, the optimum frustration regime for protein folding can be predicted analytically.

Communication: Entropic measure to prevent energy over-minimization in molecular dynamics simulations

This work examines the impact of energy over-minimization on an ensemble of biological molecules subjected to the potential energy minimization procedure in vacuum. In the studied structures, long potential energy minimization stage leads to an increase of the main- and side-chain entropies in proteins. We show that such over-minimization may diverge the protein structures from the near-native attraction basin which possesses a minimum of free energy. We propose a measure based on the Pareto front of total entropy for quality assessment of minimized protein conformation. This measure may help in selection of adequate number of energy minimization steps in protein modelling and, thus, in preservation of the near-native protein conformation.

Visualization of protein folding funnels in lattice models

Protein folding occurs in a very high dimensional phase space with an exponentially large number of states, and according to the energy landscape theory it exhibits a topology resembling a funnel. In this statistical approach, the folding mechanism is unveiled by describing the local minima in an effective one-dimensional representation. Other approaches based on potential energy landscapes address the hierarchical structure of local energy minima through disconnectivity graphs. In this paper, we introduce a metric to describe the distance between any two conformations, which also allows us to go beyond the one-dimensional representation and visualize the folding funnel in 2D and 3D. In this way it is possible to assess the folding process in detail, e.g., by identifying the connectivity between conformations and establishing the paths to reach the native state, in addition to regions where trapping may occur. Unlike the disconnectivity maps method, which is based on the kinetic connections between states, our methodology is based on structural similarities inferred from the new metric. The method was developed in a 27-mer protein lattice model, folded into a 3×3×3 cube. Five sequences were studied and distinct funnels were generated in an analysis restricted to conformations from the transition-state to the native configuration. Consistent with the expected results from the energy landscape theory, folding routes can be visualized to probe different regions of the phase space, as well as determine the difficulty in folding of the distinct sequences. Changes in the landscape due to mutations were visualized, with the comparison between wild and mutated local minima in a single map, which serves to identify different trapping regions. The extension of this approach to more realistic models and its use in combination with other approaches are discussed.

Topography of funneled landscapes determines the thermodynamics and kinetics of protein folding

The energy landscape approach has played a fundamental role in advancing our understanding of protein folding. Here, we quantify protein folding energy landscapes by exploring the underlying density of states. We identify three quantities essential for characterizing landscape topography: the stabilizing energy gap between the native and nonnative ensembles δE, the energetic roughness ΔE, and the scale of landscape measured by the entropy S. We show that the dimensionless ratio between the gap, roughness, and entropy of the system Λ=δE/(ΔE√(2S)) accurately predicts the thermodynamics, as well as the kinetics of folding. Large Λ implies that the energy gap (or landscape slope towards the native state) is dominant, leading to more funneled landscapes. We investigate the role of topological and energetic roughness for proteins of different sizes and for proteins of the same size, but with different structural topologies. The landscape topography ratio Λ is shown to be monotonically correlated with the thermodynamic stability against trapping, as characterized by the ratio of folding temperature versus trapping temperature. Furthermore, Λ also monotonically correlates with the folding kinetic rates. These results provide the quantitative bridge between the landscape topography and experimental folding measurements.

Why do protein folding rates correlate with metrics of native topology?

For almost 15 years, the experimental correlation between protein folding rates and the contact order parameter has been under scrutiny. Here, we use a simple simulation model combined with a native-centric interaction potential to investigate the physical roots of this empirical observation. We simulate a large set of circular permutants, thus eliminating dependencies of the folding rate on other protein properties (e.g. stability). We show that the rate-contact order correlation is a consequence of the fact that, in high contact order structures, the contact order of the transition state ensemble closely mirrors the contact order of the native state. This happens because, in these structures, the native topology is represented in the transition state through the formation of a network of tertiary interactions that are distinctively long-ranged.

Analyzing protein folding cooperativity by differential scanning calorimetry and NMR spectroscopy

Some marginally stable proteins undergo microsecond time scale folding reactions that involve significant populations of partly ordered forms, making it difficult to discern individual steps in their folding pathways. It has been suggested that many of these proteins fold non-cooperatively, with no significant barriers to separate the energy landscape into distinct thermodynamic states. Here we present an approach for studying the cooperativity of rapid protein folding with a combination of differential scanning calorimetry (DSC), nuclear magnetic resonance (NMR) relaxation dispersion experiments, and an analysis of the temperature dependence of amide (1)H and (15)N chemical shifts. We applied this method to the PBX homeodomain (PBX-HD), which folds on the microsecond time scale and produces a broad DSC thermogram with an elevated and steeply sloping native-state heat capacity baseline, making it a candidate for barrierless folding. However, by globally fitting the NMR thermal melt and DSC data, and by comparing these results to those obtained from the NMR relaxation dispersion experiments, we show that the native form of the protein undergoes two-state exchange with a small population of the thermally denatured form, well below the melting temperature. This result directly demonstrates the coexistence of distinct folded and unfolded forms and firmly establishes that folding of PBX-HD is cooperative. Further, we see evidence of large-scale structural and dynamical changes within the native state by NMR, which helps to explain the broad and shallow DSC profile. This study illustrates the potential of combining calorimetry with NMR dynamics experiments to dissect mechanisms of protein folding.

Comparison of the thermodynamic landscapes of unfolding and formation of the energy dissipative state in the isolated light harvesting complex II

In biochemistry and cell biology, understanding the molecular mechanisms by which physiological processes are regulated is regarded as an ultimate goal. In higher plants, one of the most widely investigated regulatory processes occurs in the light harvesting complexes (LHCII) of the chloroplast thylakoid membranes. Under limiting photon flux densities, LHCII harvests sunlight with high efficiency. When the intensity of incident radiation reaches levels close to the saturation of the photosynthesis, the efficiency of light harvesting is decreased by a process referred to as nonphotochemical quenching (NPQ), which enhances the singlet-excited state deactivation via nonradiative dissipative processes. Conformational rearrangements in LHCII are known to be crucial in promoting and controlling NPQ in vitro and in vivo. In this article, we address the thermodynamic nature of the conformational rearrangements promoting and controlling NPQ in isolated LHCII. A combined, linear reaction scheme in which the folded, quenched state represents a stable intermediate on the unfolding pathway was employed to describe the temperature dependence of the spectroscopic signatures associated with the chlorophyll fluorescence quenching and the loss of secondary structure motifs in LHCII. The thermodynamic model requires considering the temperature dependence of Gibbs free energy difference between the quenched and the unquenched states, as well as the unfolded and quenched states, of LHCII. Even though the same reaction scheme is adequate to describe the quenching and the unfolding processes in LHCII monomers and trimers, their thermodynamic characteristics were found to be markedly different. The results of the thermodynamic analysis shed light on the physiological importance of the trimeric state of LHCII in stabilizing the efficient light harvesting mode as well as preventing the quenched conformation of the protein from unfolding. Moreover, the transition to the quenched conformation in trimers reveals a larger degree of cooperativity than in monomers, explained by a small characteristic entropy (DeltaH(q) = 85 +/- 3 kJ mol(-1) compared to 125 +/- 5 kJ mol(-1) in monomers), which enables the fine-tuning of nonphotochemical quenching in vivo.

Uncovering the properties of energy-weighted conformation space networks with a hydrophobic-hydrophilic model

The conformation spaces generated by short hydrophobic-hydrophilic (HP) lattice chains are mapped to conformation space networks (CSNs). The vertices (nodes) of the network are the conformations and the links are the transitions between them. It has been found that these networks have "small-world" properties without considering the interaction energy of the monomers in the chain, i. e. the hydrophobic or hydrophilic amino acids inside the chain. When the weight based on the interaction energy of the monomers in the chain is added to the CSNs, it is found that the weighted networks show the "scale-free" characteristic. In addition, it reveals that there is a connection between the scale-free property of the weighted CSN and the folding dynamics of the chain by investigating the relationship between the scale-free structure of the weighted CSN and the noted parameter Z score. Moreover, the modular (community) structure of weighted CSNs is also studied. These results are helpful to understand the topological properties of the CSN and the underlying free-energy landscapes.

Mutational robustness and geometrical form in protein structures

Theoretical studies of RNA and lattice protein models suggest that mutationally robust or the so-called designable phenotypes tend to have special geometric features such as being more compact and more geometrically regular. Such geometrical forms have been also linked to speed of folding and stability properties that may also assist in promoting mutational robustness. Here we test these theoretical predictions on a non-redundant collection of 2,660 experimentally determined structures from the PDB (Protein Data Bank) and CATH (Class Architecture Topology Homologous superfamily) database. We first developed an index summarizing the geometrical regularity of the structures and then used this index to show that the statistical pattern of empirical data is consistent with the theoretical predictions relating geometry to mutational robustness. Mutationally robust proteins tend to be more symmetric and compact. But, the relationship between compactness and robustness cannot be explained simply by the geometrical packing of individual amino acids in proteins; rather, it is the property of the whole system that is related to the statistical characteristics of the folding landscape. Finally, we hypothesize that a triplet relationship between mutational robustness, stability and form is a general properties of objects that optimize real-valued relationships between sequences and discrete structures.

Kinetic barriers and the role of topology in protein and RNA folding

This review compares the folding behavior of proteins and RNAs. Topics covered include the role of topology in the determination of folding rates, major folding events including collapse, properties of denatured states, pathway heterogeneity, and the influence of the mode of initiation on the folding pathway.

Pathway to copolymer collapse in dilute solution: uniform versus random distribution of comonomers

Monte Carlo simulations show that copolymers with uniformly (or periodically) distributed sticky comonomers collapse "cooperatively," abruptly forming a compact intermediate comprising a monomer shell surrounding a core of the aggregated comonomers. In comparison, random copolymers collapse through a relatively less-compact intermediate comprising a comonomer core surrounded by a fluffy monomer shell that densifies over a wide temperature range. This difference between the collapse pathways for random and uniform copolymers persists to higher chain lengths, where uniform copolymers tend to form multiple comonomer cores. In this paper, we describe the formation of such an intermediate state, and the subsequent collapse, by recognizing that these arise from the expected balance between comonomer aggregation enthalpy and loop formation entropy dictated by the chain microstructure.

Configuration-dependent diffusion can shift the kinetic transition state and barrier height of protein folding

We show that diffusion can play an important role in protein-folding kinetics. We explicitly calculate the diffusion coefficient of protein folding in a lattice model. We found that diffusion typically is configuration- or reaction coordinate-dependent. The diffusion coefficient is found to be decreasing with respect to the progression of folding toward the native state, which is caused by the collapse to a compact state constraining the configurational space for exploration. The configuration- or position-dependent diffusion coefficient has a significant contribution to the kinetics in addition to the thermodynamic free-energy barrier. It effectively changes (increases in this case) the kinetic barrier height as well as the position of the corresponding transition state and therefore modifies the folding kinetic rates as well as the kinetic routes. The resulting folding time, by considering both kinetic diffusion and the thermodynamic folding free-energy profile, thus is slower than the estimation from the thermodynamic free-energy barrier with constant diffusion but is consistent with the results from kinetic simulations. The configuration- or coordinate-dependent diffusion is especially important with respect to fast folding, when there is a small or no free-energy barrier and kinetics is controlled by diffusion. Including the configurational dependence will challenge the transition state theory of protein folding. The classical transition state theory will have to be modified to be consistent. The more detailed folding mechanistic studies involving phi value analysis based on the classical transition state theory also will have to be modified quantitatively.

Sampling of states for estimating the folding funnel entropy and energy landscape of a model alpha-helical hairpin peptide

Protein folding times are many orders of magnitude shorter than would occur if the peptide chain randomly sampled possible configurations, which implies that protein folding is a directed process. The detailed shape of protein's energy landscape determines the rate and reliability of folding to the native state, but the large number of structural degrees of freedom generates an energy landscape that is hard to visualize because of its high dimensionality. A commonly used picture is that of an energy funnel leading from high energy random coil state down to the low energy native state. As lattice computer models of protein dynamics become more realistic, the number of possible configurations becomes too large to count directly. Statistical mechanic and thermodynamic approaches allow us to count states in an approximate manner to quantify the entropy and energy of the energy landscape within a folding funnel for an alpha-helical protein. We also discuss the problems that arise in attempting to count the huge number of individual states of the random coil at the top of the funnel.

Folding time distributions as an approach to protein folding kinetics

A 27-residue lattice heteropolymer subject to Monte Carlo dynamics on a simple cubic lattice is studied over a range of temperatures. Folding time distributions are used to obtain information concerning the details of folding kinetics. The results are compared with those from methods based on mean force surfaces expressed in terms of a reduced set of variables and on a disconnectivity graph for the same system. A detailed analysis of the folding trajectories is given, and the importance of dead-end traps in determining the folding time is demonstrated. We show that the calculated folding kinetics can be modeled by a system of kinetic equations, with the essential rate constants determined from the Monte Carlo simulations and the resulting folding time distributions. The kinetic equations make possible an analysis of the variation of the importance of different channels with temperature. In particular, we show that the presence of intermediates may be masked in the folding time distributions, with the mean folding time being independent of the height of the barrier between the intermediates and collapsed globule state of the system. This and other results demonstrate that care has to be used in interpreting experimental folding data in terms of the underlying kinetics. Correspondingly, simulations are shown to have to satisfy certain requirements to obtain proper sampling of the dead-end traps.

A ternary nucleation model for the nucleation pathway of protein folding

Recently [Y. S. Djikaev and E. Ruckenstein, J. Phys. Chem. B 111, 886 (2007)], the authors proposed a kinetic model for the nucleation mechanism of protein folding where a protein was modeled as a heteropolymer consisting of hydrophobic and hydrophilic beads and the composition of the growing cluster of protein residues was assumed to be constant and equal to the overall protein composition. Here, they further develop the model by considering a protein as a three-component heteropolymer and by allowing the composition of the growing cluster of protein residues to vary independently of the overall one. All the bonds in the heteropolymer (now consisting of hydrophobic, hydrophilic, and neutral beads) have the same constant length, and all the bond angles are equal and fixed. As a crucial idea of the model, an overall potential around the cluster wherein a residue performs a chaotic motion is considered to be a combination of the average dihedral and average pairwise potentials assigned to the bead. The overall potential as a function of the distance from the cluster center has a double well shape which allows one to determine its emission and absorption rates by using a first passage time analysis. Knowing these rates as functions of three independent variables of a ternary cluster, one can develop a self-consistent kinetic theory for the nucleation mechanism of folding of a protein using a ternary nucleation formalism and evaluate the size and composition of the nucleus and the protein folding time. As an illustration, the model is applied to the folding of bovine pancreatic ribonuclease consisting of 124 amino acids whereof 40 are hydrophobic, 81 hydrophilic, and 3 neutral. With a reasonable choice of diffusion coefficients of the residues in the native state and potential parameters, the model predicts folding times in the range of 1-100 s.

A lattice protein with an amyloidogenic latent state: stability and folding kinetics

We have designed a model lattice protein that has two stable folded states, the lower free energy native state and a latent state of somewhat higher energy. The two states have a sizable part of their structures in common (two "alpha-helices") and differ in the content of "alpha-helices" and "beta-strands" in the rest of their structures; i.e. for the native state, this part is alpha-helical, and for the latent state it is composed of beta-strands. Thus, the lattice protein free energy surface mimics that of amyloidogenic proteins that form well organized fibrils under appropriate conditions. A Go-like potential was used and the folding process was simulated with a Monte Carlo method. To gain insight into the equilibrium free energy surface and the folding kinetics, we have combined standard approaches (reduced free energy surfaces, contact maps, time-dependent populations of the characteristic states, and folding time distributions) with a new approach. The latter is based on a principal coordinate analysis of the entire set of contacts, which makes possible the introduction of unbiased reaction coordinates and the construction of a kinetic network for the folding process. The system is found to have four characteristic basins, namely a semicompact globule, an on-pathway intermediate (the bifurcation basin), and the native and latent states. The bifurcation basin is shallow and consists of the structure common to the native and latent states, with the rest disorganized. On the basis of the simulation results, a simple kinetic model describing the transitions between the characteristic states was developed, and the rate constants for the essential transitions were estimated. During the folding process the system dwells in the bifurcation basin for a relatively short time before it proceeds to the native or latent state. We suggest that such a bifurcation may occur generally for proteins in which native and latent states have a sizable part of their structures in common. Moreover, there is the possibility of introducing changes in the system (e.g., mutations), which guide the system toward the native or misfolded state.

Model for the Nucleation Mechanism of Protein Folding

A nucleation-like pathway of protein folding involves the formation of a cluster containing native residues that grows by including residues from the unfolded part of the protein. This pathway is examined by using a heteropolymer as a protein model. The model heteropolymer consists of hydrophobic and hydrophilic beads with fixed bond lengths and bond angles. The total energy of the heteropolymer is determined by the pairwise repulsive/attractive interactions between nonlinked beads and by the contribution from the dihedral angles involved. The parameters of these interactions can be rigorously defined, unlike the ill-defined surface tension of a cluster of protein residues that constitutes the basis of a previous nucleation model. The main idea underlying the new model consists of averaging the dihedral potential of a selected residue over all possible configurations of all neighboring residues along the protein chain. The resulting average dihedral potential depends on the distance between the selected residue and the cluster center. Its combination with the average pairwise potential of the selected residue and with a confining potential caused by the bonds between the residues leads to an overall potential around the cluster that has a double-well shape. Residues in the inner (closer to the cluster) well are considered as belonging to the folded cluster, whereas those in the outer well are treated as belonging to the unfolded part of the protein. Transitions of residues from the inner well into the outer one and vice versa are considered as elementary emission and absorption events, respectively. The double-well character of the potential well around the cluster allows one to determine the rates of both emission and absorption of residues by the cluster using a first passage time analysis. Once these rates are found as functions of the cluster size, one can develop a self-consistent kinetic theory for the nucleation mechanism of folding of a protein. The model allows one to evaluate the size of the nucleus and the protein folding time. The latter is evaluated as the sum of the times necessary for the first nucleation event to occur and for the nucleus to grow to the maximum size (of the folded protein). Depending on the diffusion coefficients of the native residues in the range from 10(-6) to 10(-8) cm2/s, numerical calculations for a protein of 2500 residues suggest that the folding time ranges from several seconds to several hundreds of seconds.

Criteria for downhill protein folding: Calorimetry, chevron plot, kinetic relaxation, and single-molecule radius of gyration in chain models with subdued degrees of cooperativity

Recent investigations of possible downhill folding of small proteins such as BBL have focused on the thermodynamics of non-two-state, "barrierless" folding/denaturation transitions. Downhill folding is noncooperative and thermodynamically "one-state," a phenomenon underpinned by a unimodal conformational distribution over chain properties such as enthalpy, hydrophobic exposure, and conformational dimension. In contrast, corresponding distributions for cooperative two-state folding are bimodal with well-separated population peaks. Using simplified atomic modeling of a three-helix bundle-in a scheme that accounts for hydrophobic interactions and hydrogen bonding-and coarse-grained C(alpha) models of four real proteins with various degrees of cooperativity, we evaluate the effectiveness of several observables at defining the underlying distribution. Bimodal distributions generally lead to sharper transitions, with a higher heat capacity peak at the transition midpoint, compared with unimodal distributions. However, the observation of a sigmoidal transition is not a reliable criterion for two-state behavior, and the heat capacity baselines, used to determine the van't Hoff and calorimetric enthalpies of the transition, can introduce ambiguity. Interestingly we find that, if the distribution of the single-molecule radius of gyration were available, it would permit discrimination between unimodal and bimodal underlying distributions. We investigate kinetic implications of thermodynamic noncooperativity using Langevin dynamics. Despite substantial chevron rollovers, the relaxation of the models considered is essentially single-exponential over an extended range of native stabilities. Consistent with experiments, significant deviations from single-exponential behavior occur only under strongly folding conditions.

Quantifying the kinetic paths of flexible biomolecular recognition

Biomolecular recognition often involves large conformational changes, sometimes even local unfolding. The identification of kinetic pathways has become a central issue in understanding the nature of binding. A new approach is proposed here to study the dynamics of this binding-folding process through the establishment of a path-integral framework on the underlying energy landscape. The dominant kinetic paths of binding and folding can be determined and quantified. The significant coupling between the binding and folding of biomolecules often exists in many important cellular processes. In this case, the corresponding kinetic paths of binding are shown to be intimately correlated with those of folding and the dynamics becomes quite cooperative. This implies that binding and folding happen concurrently. When the coupling between binding and folding is weak (strong), the kinetic process usually starts with significant folding (binding) first, with the binding (folding) later proceeding to the end. The kinetic rate can be obtained through the contributions from the dominant paths. The rate is shown to have a bell-shaped dependence on temperature in the concentration-saturated regime consistent with experiment. The changes of the kinetics that occur upon changing the parameters of the underlying binding-folding energy landscape are studied.

A protein evolution model with independent sites that reproduces site-specific amino acid distributions from the Protein Data Bank

BACKGROUND

Since thermodynamic stability is a global property of proteins that has to be conserved during evolution, the selective pressure at a given site of a protein sequence depends on the amino acids present at other sites. However, models of molecular evolution that aim at reconstructing the evolutionary history of macromolecules become computationally intractable if such correlations between sites are explicitly taken into account.

RESULTS

We introduce an evolutionary model with sites evolving independently under a global constraint on the conservation of structural stability. This model consists of a selection process, which depends on two hydrophobicity parameters that can be computed from protein sequences without any fit, and a mutation process for which we consider various models. It reproduces quantitatively the results of Structurally Constrained Neutral (SCN) simulations of protein evolution in which the stability of the native state is explicitly computed and conserved. We then compare the predicted site-specific amino acid distributions with those sampled from the Protein Data Bank (PDB). The parameters of the mutation model, whose number varies between zero and five, are fitted from the data. The mean correlation coefficient between predicted and observed site-specific amino acid distributions is larger than <r> = 0.70 for a mutation model with no free parameters and no genetic code. In contrast, considering only the mutation process with no selection yields a mean correlation coefficient of <r> = 0.56 with three fitted parameters. The mutation model that best fits the data takes into account increased mutation rate at CpG dinucleotides, yielding <r> = 0.90 with five parameters.

CONCLUSION

The effective selection process that we propose reproduces well amino acid distributions as observed in the protein sequences in the PDB. Its simplicity makes it very promising for likelihood calculations in phylogenetic studies. Interestingly, in this approach the mutation process influences the effective selection process, i.e. selection and mutation must be entangled in order to obtain effectively independent sites. This interdependence between mutation and selection reflects the deep influence that mutation has on the evolutionary process: The bias in the mutation influences the thermodynamic properties of the evolving proteins, in agreement with comparative studies of bacterial proteomes, and it also influences the rate of accepted mutations.

Measurements of protein sequence-structure correlations

Correlations between protein structures and amino acid sequences are widely used for protein structure prediction. For example, secondary structure predictors generally use correlations between a secondary structure sequence and corresponding primary structure sequence, whereas threading algorithms and similar tertiary structure predictors typically incorporate interresidue contact potentials. To investigate the relative importance of these sequence-structure interactions, we measured the mutual information among the primary structure, secondary structure and side-chain surface exposure, both for adjacent residues along the amino acid sequence and for tertiary structure contacts between residues distantly separated along the backbone. We found that local interactions along the amino acid chain are far more important than non-local contacts and that correlations between proximate amino acids are essentially uninformative. This suggests that knowledge-based contact potentials may be less important for structure predication than is generally believed.

Multiple folding mechanisms of protein ubiquitin

Based on the C(alpha) Go-type model, the folding kinetics and mechanisms of protein ubiquitin with mixed alpha/beta topology are studied by molecular dynamics simulations. The relaxation kinetics shows that there are three phases, namely the major phase, the intermediate phase and the slowest minor phase. The existence of these three phases are relevant to the phenomenon found in experiments. According to our simulations, the folding at high temperatures around the folding transition temperature T(f) is of a two-state process, and the folding nucleus is consisted of contacts between the front end of alpha-helix and the turn(4). The folding at low temperature (approximately T = 0.8) is also studied, where an A-state like structure is found lying on the major folding pathway. The appearance of this structure is related to the stability of the first part (residue 1-51) of protein ubiquitin. As the temperature decreases, the formation of secondary structures, tertiary structures and collapse of the protein are found to be decoupled gradually and the folding mechanism changes from the nucleation-condensation to the diffusion-collision. This feature indicates a unifying common folding mechanism for proteins. The intermediate phase is also studied and is found to represent a folding process via a long-lived intermediate state which is stabilized by strong interactions between the beta(1) and the beta(5) strand. These strong interactions are important for the function of protein ubiquitin as a molecular chaperone. Thus the intermediate phase is assumed as a byproduct of the requirement of protein function. In addition, the validity of the current Go-model is also investigated, and a lower limited temperature for protein ubiquitin T(limit) = 0.8 is proposed. At temperatures higher than this value, the kinetic traps due to glass dynamics cannot be significantly populated and the intermediate states can be reliably identified although there is slight chevron rollover in the folding rates. At temperature lower than T(limit), however, the traps due to glass dynamics become dominant and may be mistaken for real intermediate states. This limitation of valid temperature range prevents us to reveal the burst phase intermediate in the major folding phase since it might only be stabilized at temperatures lower than T(limit), according to experiments. Our works show that caution must be taken when studying low-temperature intermediate states by using the C(alpha) Go-models.

Principal eigenvector of contact matrices and hydrophobicity profiles in proteins

With the aim of studying the relationship between protein sequences and their native structures, we adopted vectorial representations for both sequence and structure. The structural representation was based on the principal eigenvector of the fold's contact matrix (PE). As has been recently shown, the latter encodes sufficient information for reconstructing the whole contact matrix. The sequence was represented through a hydrophobicity profile (HP), using a generalized hydrophobicity scale that we obtained from the principal eigenvector of a residue-residue interaction matrix, and denoted as interactivity scale. Using this novel scale, we defined the optimal HP of a protein fold, and, by means of stability arguments, predicted to be strongly correlated with the PE of the fold's contact matrix. This prediction was confirmed through an evolutionary analysis, which showed that the PE correlates with the HP of each individual sequence adopting the same fold and, even more strongly, with the average HP of this set of sequences. Thus, protein sequences evolve in such a way that their average HP is close to the optimal one, implying that neutral evolution can be viewed as a kind of motion in sequence space around the optimal HP. Our results indicate that the correlation coefficient between N-dimensional vectors constitutes a natural metric in the vectorial space in which we represent both protein sequences and protein structures, which we call vectorial protein space. In this way, we define a unified framework for sequence-to-sequence, sequence-to-structure and structure-to-structure alignments. We show that the interactivity scale is nearly optimal both for the comparison of sequences to sequences and sequences to structures.

Reconstruction of the src-SH3 Protein Domain Transition State Ensemble using Multiscale Molecular Dynamics Simulations

We use an integrated computational approach to reconstruct accurately the transition state ensemble (TSE) for folding of the src-SH3 protein domain. We first identify putative TSE conformations from free energy surfaces generated by importance sampling molecular dynamics for a fully atomic, solvated model of the src-SH3 protein domain. These putative TSE conformations are then subjected to a folding analysis using a coarse-grained representation of the protein and rapid discrete molecular dynamics simulations. Those conformations that fold to the native conformation with a probability (P(fold)) of approximately 0.5, constitute the true transition state. Approximately 20% of the putative TSE structures were found to have a P(fold) near 0.5, indicating that, although correct TSE conformations are populated at the free energy barrier, there is a critical need to refine this ensemble. Our simulations indicate that the true TSE conformations are compact, with a well-defined central beta sheet, in good agreement with previous experimental and theoretical studies. A structured central beta sheet was found to be present in a number of pre-TSE conformations, however, indicating that this element, although required in the transition state, does not define it uniquely. An additional tight cluster of contacts between highly conserved residues belonging to the diverging turn and second beta-sheet of the protein emerged as being critical elements of the folding nucleus. A number of commonly used order parameters to identify the transition state for folding were investigated, with the number of native Cbeta contacts displaying the most satisfactory correlation with P(fold) values.

Chevron behavior and isostable enthalpic barriers in protein folding: successes and limitations of simple Gō-like modeling

It has been demonstrated that a "near-Levinthal" cooperative mechanism, whereby the common Gō interaction scheme is augmented by an extra favorability for the native state as a whole, can lead to apparent two-state folding/unfolding kinetics over a broad range of native stabilities in lattice models of proteins. Here such a mechanism is shown to be generalizable to a simplified continuum (off-lattice) Langevin dynamics model with a Calpha protein chain representation, with the resulting chevron plots exhibiting an extended quasilinear regime reminiscent of that of apparent two-state real proteins. Similarly high degrees of cooperativity are possible in Gō-like continuum models with rudimentary pairwise desolvation barriers as well. In these models, cooperativity increases with increasing desolvation barrier height, suggesting strongly that two-state-like folding/unfolding kinetics would be achievable when the pairwise desolvation barrier becomes sufficiently high. Besides cooperativity, another generic folding property of interest that has emerged from published experiments on several apparent two-state proteins is that their folding relaxation under constant native stability (isostability) conditions is essentially Arrhenius, entailing high intrinsic enthalpic folding barriers of approximately 17-30 kcal/mol. Based on a new analysis of published data on barnase, here we propose that a similar property should also apply to a certain class of non-two-state proteins that fold with chevron rollovers. However, several continuum Gō-like constructs considered here fail to predict any significant intrinsic enthalpic folding barrier under isostability conditions; thus the physical origin of such barriers in real proteins remains to be elucidated.

The complex kinetics of protein folding in wide temperature ranges

The complex protein folding kinetics in wide temperature ranges is studied through diffusive dynamics on the underlying energy landscape. The well-known kinetic chevron rollover behavior is recovered from the mean first passage time, with the U-shape dependence on temperature. The fastest folding temperature T0 is found to be smaller than the folding transition temperature Tf. We found that the fluctuations of the kinetics through the distribution of first passage time show rather universal behavior, from high-temperature exponential Poissonian kinetics to the relatively low-temperature highly non-exponential kinetics. The transition temperature is at Tk and T0 < Tk < Tf. In certain low-temperature regimes, a power law behavior at long time emerges. At very low temperatures (lower than trapping transition temperature T < T0/(4 approximately 6)), the kinetics is an exponential Poissonian process again.

Downhill kinetics of biomolecular interface binding: globally connected scenario

We study the kinetics of the biomolecular binding process at the interface using energy landscape theory. The global kinetic connectivity case is considered for a downhill funneled energy landscape. By solving the kinetic master equation, the kinetic time for binding is obtained and shown to have a U-shape curve-dependence on the temperature. The kinetic minimum of the binding time monotonically decreases when the ratio of the underlying energy gap between native state and average non-native states versus the roughness or the fluctuations of the landscape increases. At intermediate temperatures, fluctuations measured by the higher moments of the binding time lead to non-Poissonian, non-exponential kinetics. At both high and very low temperatures, the kinetics is nearly Poissonian and exponential.

Looking at structure, stability, and evolution of proteins through the principal eigenvector of contact matrices and hydrophobicity profiles

We review and further develop an analytical model that describes how thermodynamic constraints on the stability of the native state influence protein evolution in a site-specific manner. To this end, we represent both protein sequences and protein structures as vectors: structures are represented by the principal eigenvector (PE) of the protein contact matrix, a quantity that resembles closely the effective connectivity of each site; sequences are represented through the "interactivity" of each amino acid type, using novel parameters that are correlated with hydropathy scales. These interactivity parameters are more strongly correlated than the other hydropathy scales that we examine with: (1) the change upon mutations of the unfolding free energy of proteins with two-states thermodynamics; (2) genomic properties as the genome-size and the genome-wide GC content; (3) the main eigenvectors of the substitution matrices. The evolutionary average of the interactivity vector correlates very strongly with the PE of a protein structure. Using this result, we derive an analytic expression for site-specific distributions of amino acids across protein families in the form of Boltzmann distributions whose "inverse temperature" is a function of the PE component. We show that our predictions are in agreement with site-specific amino acid distributions obtained from the Protein Data Bank, and we determine the mutational model that best fits the observed site-specific amino acid distributions. Interestingly, the optimal model almost minimizes the rate at which deleterious mutations are eliminated by natural selection.

Prediction of site-specific amino acid distributions and limits of divergent evolutionary changes in protein sequences

We derive an analytic expression for site-specific stationary distributions of amino acids from the structurally constrained neutral (SCN) model of protein evolution with conservation of folding stability. The stationary distributions that we obtain have a Boltzmann-like shape, and their effective temperature parameter, measuring the limit of divergent evolutionary changes at a given site, can be predicted from a site-specific topological property, the principal eigenvector of the contact matrix of the native conformation of the protein. These analytic results, obtained without free parameters, are compared with simulations of the SCN model and with the site-specific amino acid distributions obtained from the Protein Data Bank. These results also provide new insights into how the topology of a protein fold influences its designability, i.e., the number of sequences compatible with that fold. The dependence of the effective temperature on the principal eigenvector decreases for longer proteins, as a possible consequence of the fact that selection for thermodynamic stability becomes weaker in this case.

Genomic determinants of protein folding thermodynamics in prokaryotic organisms

Here we investigate how thermodynamic properties of orthologous proteins are influenced by the genomic environment in which they evolve. We performed a comparative computational study of 21 protein families in 73 prokaryotic species and obtained the following main results. (i) Protein stability with respect to the unfolded state and with respect to misfolding are anticorrelated. There appears to be a trade-off between these two properties, which cannot be optimized simultaneously. (ii) Folding thermodynamic parameters are strongly correlated with two genomic features, genome size and G+C composition. In particular, the normalized energy gap, an indicator of folding efficiency in statistical mechanical models of protein folding, is smaller in proteins of organisms with a small genome size and a compositional bias towards A+T. Such genomic features are characteristic for bacteria with an intracellular lifestyle. We interpret these correlations in light of mutation pressure and natural selection. A mutational bias toward A+T at the DNA level translates into a mutational bias toward more hydrophobic (and in general more interactive) proteins, a consequence of the structure of the genetic code. Increased hydrophobicity renders proteins more stable against unfolding but less stable against misfolding. Proteins with high hydrophobicity and low stability against misfolding occur in organisms with reduced genomes, like obligate intracellular bacteria. We argue that they are fixed because these organisms experience weaker purifying selection due to their small effective population sizes. This interpretation is supported by the observation of a high expression level of chaperones in these bacteria. Our results indicate that the mutational spectrum of a genome and the strength of selection significantly influence protein folding thermodynamics.

Probing the kinetics of single molecule protein folding

We propose an approach to integrate the theory, simulations, and experiments in protein-folding kinetics. This is realized by measuring the mean and high-order moments of the first-passage time and its associated distribution. The full kinetics is revealed in the current theoretical framework through these measurements. In the experiments, information about the statistical properties of first-passage times can be obtained from the kinetic folding trajectories of single molecule experiments (for example, fluorescence). Theoretical/simulation and experimental approaches can be directly related. We study in particular the temperature-varying kinetics to probe the underlying structure of the folding energy landscape. At high temperatures, exponential kinetics is observed; there are multiple parallel kinetic paths leading to the native state. At intermediate temperatures, nonexponential kinetics appears, revealing the nature of the distribution of local traps on the landscape and, as a result, discrete kinetic paths emerge. At very low temperatures, exponential kinetics is again observed; the dynamics on the underlying landscape is dominated by a single barrier. The ratio between first-passage-time moments is proposed to be a good variable to quantitatively probe these kinetic changes. The temperature-dependent kinetics is consistent with the strange kinetics found in folding dynamics experiments. The potential applications of the current results to single-molecule protein folding are discussed.

3D-Lattice Monte Carlo simulations of model proteins. Size effects on folding thermodynamics and kinetics

Recently, we devised an energy scale to vary systematically amino-acid residue-solvent interactions for Monte Carlo simulations of lattice-model proteins in water. For 27-mer proteins, the folding behavior varies appreciably with the choice of interaction parameters. We now perform similar simulations with 64-mers to study the size dependence of the optimal energy parameter set for representing realistic behavior typical of many real proteins (i.e. fast folding and high cooperativity for single chains). We find that 64-mers are considerably more stable and more cooperative compared to 27-mers. The optimal interfacial-interaction-energy parameter set, however, is relatively size independent.

Folding probabilities: A novel approach to folding transitions and the two-dimensional Ising-model

The theoretical concept of folding probability, p(fold), has proven to be a useful means to characterize the kinetics of protein folding. Here, we illustrate the practical importance of p(fold) and demonstrate how it can be determined theoretically. We derive a general analytical expression for p(fold) and show how it can be estimated from simulations for systems where the transition rates between the relevant microstates are not known. By analyzing the Ising model we are able to determine the scaling behavior of the numerical error in the p(fold) estimate as function of the number of analyzed Monte Carlo runs. We apply our method to a simple, newly developed protein folding model for the formation of alpha helices. It is demonstrated that our technique highly parallelizes the calculation of p(fold) and that it is orders of magnitude more efficient than conventional approaches.

Single-molecule studies highlight conformational heterogeneity in the early folding steps of a large ribozyme

The equilibrium folding of the catalytic domain of Bacillus subtilis RNase P RNA is investigated by single-molecule fluorescence resonance energy transfer (FRET). Previous ensemble studies of this 255-nucleotide ribozyme described the equilibrium folding with two transitions, U-to-I(eq)-to-N, and focused on the I(eq)-to-N transition. The present study focuses on the U-to-I(eq) transition. Comparative ensemble measurements of the ribozyme construct labeled with fluorescein at the 5' end and Cy3 at the 3' end show that modifications required for labeling do not interfere with folding and help to define the Mg(2+) concentration range for the U-to-I(eq) transition. Histogram analysis of the Mg(2+)-dependent single-molecule FRET efficiency reveals two previously undetermined folding intermediates. The single-molecule FRET trajectories exhibit non-two-state and nonergodic behaviors at intermediate Mg(2+) concentrations on the time scale of seconds. The trajectories at intermediate Mg(2+) concentrations are classified into five classes based on three FRET levels and their dynamics of interconversion within the measured time range. This heterogeneity, together with the observation of "nonsudden jump" FRET transitions, indicates that the early folding steps of this ribozyme involve a series of intermediates with different degrees of kinetic isolation and that folding occurs under kinetic control and involves many "local" conformational switches. A free energy contour is constructed to illustrate the complex folding surface.

Simple two-state protein folding kinetics requires near-levinthal thermodynamic cooperativity

Simple two-state folding kinetics of many small single-domain proteins are characterized by chevron plots with linear folding and unfolding arms consistent with an apparent two-state description of equilibrium thermodynamics. This phenomenon is hereby recognized as a nontrivial heteropolymer property capable of providing fundamental insight into protein energetics. Many current protein chain models, including common lattice and continuum Gō models with explicit native biases, fail to reproduce this generic protein property. Here we show that simple two-state kinetics is obtainable from models with a cooperative interplay between core burial and local conformational propensities or an extra strongly favorable energy for the native structure. These predictions suggest that intramolecular recognition in real two-state proteins is more specific than that envisioned by common Gō-like constructs with pairwise additive energies. The many-body interactions in the present kinetically two-state models lead to high thermodynamic cooperativity as measured by their van't Hoff to calorimetric enthalpy ratios, implying that the native and denatured conformational populations are well separated in enthalpy by a high free-energy barrier. It has been observed experimentally that deviations from Arrhenius behavior are often more severe for folding than for unfolding. This asymmetry may be rationalized by one of the present modeling scenarios if the effective many-body cooperative interactions stabilizing the native structure against unfolding is less dependent on temperature than the interactions that drive the folding kinetics.

Contact order dependent protein folding rates: kinetic consequences of a cooperative interplay between favorable nonlocal interactions and local conformational preferences

Physical mechanisms underlying the empirical correlation between relative contact order (CO) and folding rate among naturally occurring small single-domain proteins are investigated by evaluating postulated interaction schemes for a set of three-dimensional 27mer lattice protein models with 97 different CO values. Many-body interactions are constructed such that contact energies become more favorable when short chain segments sequentially adjacent to the contacting residues adopt native-like conformations. At a given interaction strength, this scheme leads to folding rates that are logarithmically well correlated with CO (correlation coefficient r = 0.914) and span more than 2.5 orders of magnitude, whereas folding rates of the corresponding Gō models with additive contact energies have much less logarithmic correlation with CO and span only approximately one order of magnitude. The present protein chain models also exhibit calorimetric cooperativity and linear chevron plots similar to that observed experimentally for proteins with apparent simple two-state folding/unfolding kinetics. Thus, our findings suggest that CO-dependent folding rates of real proteins may arise partly from a significant positive coupling between nonlocal contact favorabilities and local conformational preferences.

Origins of chevron rollovers in non-two-state protein folding kinetics

Chevron rollovers of some proteins imply that their logarithmic folding rates are nonlinear in native stability. This is predicted by lattice and continuum Gō models to arise from diminished accessibilities of the ground state from transiently populated compact conformations under strongly native conditions. Despite these models' native-centric interactions, the slowdown is due partly to kinetic trapping caused by some of the folding intermediates' non-native topologies. Notably, simple two-state folding kinetics of small single-domain proteins are not reproduced by common Gō-like schemes.

Structural determinant of protein designability

Here we present an approximate analytical theory for the relationship between a protein structure's contact matrix and the shape of its energy spectrum in amino acid sequence space. We demonstrate a dependence of the number of sequences of low energy in a structure on the eigenvalues of the structure's contact matrix, and then use a Monte Carlo simulation to test the applicability of this analytical result to cubic lattice proteins. We find that the lattice structures with the most low-energy sequences are the same as those predicted by the theory. We argue that, given sufficiently strict requirements for foldability, these structures are the most designable, and we propose a simple means to test whether the results in this paper hold true for real proteins.