The first phase is akin to a romantic stage of research, where the final goal of studies may not be directly applicable to a broader understanding, or exploitable in a relevant science. The final goal is to determine the basic principles that relate protein sequence and structure. The second phase is a more pragmatic stage of research, where the applications drive research in the field and the rational manipulation of derived knowledge allows engineering of tools for advancement of a relevant science. For example, understanding the functional intermediates that accompany the transition of a protein en route to its native state may allow rational manipulation of protein structure via protein design. This example not only relates protein sequence, structure and function, but also demonstrates the engineering aspect of the modern protein folding field.

The discovery of a class of simple, single-domain proteins which fold via two-state kinetics without any detectable intermediates in the early 1990s [6, 7], the development of experimental techniques with improved spatial/temporal resolution[8-13], and the application of computer simulations using simplified lattice and off-lattice models[14, 15] greatly enhanced our understanding of various aspects of the protein folding problem. Based on the nucleation theory [16-18], one of the early proposed mechanisms for protein folding, the nucleation-condensation model was formulated[19-21]. In this scenario, a small number of residues (folding nucleus) need to form their native contacts in order for the folding reaction to proceed fast into the native state. The cooperativity of the protein folding process is analogous to that exhibited in first-order phase transitions, which proceed via a nucleation and growth mechanism [22]. Because of these similarities, terminology used in studies of phase transitions, such as energy landscapes and nucleation, was introduced into the discussion of protein folding. The concepts of the nucleation and the free energy landscape have promoted much of the recent progress in understanding the process of protein folding. Proteins are generally thought to have evolved to exhibit globally funneled energy landscapes [23-25] which allow proteins to fold to their native states through a stochastic process in which the free energy decreases spontaneously. The unfolded state, transition state, native state and possible intermediates correspond to local minima or saddle points in the free-energy landscape.

Advances in experimental techniques such as protein engineering, nuclear magnetic resonance (NMR), mass spectrometry, hydrogen exchange, fluorescence resonance energy transfer (FRET), atomic force microscopy (AFM), have made it possible to obtain detailed information about the different conformations occurring in the folding process[26, 27]. At the same time, computational methods have been developed to better interpret experimental data by using simulations to obtain structural information about the states which are populated during the folding process. In Table 1, we list several advances in experimental and computational methodologies used for investigating the folding of model proteins.

All-atom protein models with explicit or implicit solvents were developed to study the folding thermodynamics and the unfolding dynamics of specific proteins. Technological advances in computation allowed folding simulations of small proteins and peptides at atomic detail [28-30]. However, due to the complexity and vast dimensionality of protein conformational space, all-atom MD simulations have severe limitations on the time and length scales that can be studied. Novel simulation protocols have been proposed to improve conformational sampling efficiency, including biased sampling of the free-energy surface and non-equilibrium unfolding simulations [24]. In addition, world-wide parallel computing (e.g. Folding/at/Home [31]) and generalized ensemble sampling techniques that involve parallel simulations of molecular systems coupled with a Monte Carlo (MC) protocol [32, 33] have been successfully applied to protein folding [25, 34-36].

Multi-scale modeling approaches have also been used to combine efficient conformational sampling of coarse-grained models and accuracy of all-atom models to study protein folding pathways. In this approach, iterative simulations and inter-conversion between high and low-resolution protein models are performed. Feig et al. developed a multi-scale modeling tool set, MMTSB [37], which integrates a simplified protein model with the MC simulation engine, MONSSTER [38], and the all-atom MD packages AMBER [39] or CHARMM[40]. Using a combination of CHARMM and discrete molecular dynamics (DMD) [41-46], Ding et al. reconstructed the transition state ensemble of the src-SH3 protein domain through multi-scale simulations [47]. The protein folding studies can also be facilitated by sampling protein conformations near the native state. Several native-state sampling algorithms [48, 49] have been successfully utilized to study plasticity [50], cooperative interactions [51], and allostery [52] in proteins. Considering native-state ensemble naturally takes into account protein flexibility, which is shown to be crucial in structure based drug designs.

During the last five years, several tools for performing web-based analyses of protein folding dynamics have been developed. The Fold-Rate server (http://psfs.cbrc.jp/fold-rate/) [53] predicts rates of protein folding using the amino-acid sequence. The Parasol folding server (http://parasol.tamu.edu/groups/amatogroup/foldingserver) [54] predicts protein folding pathways using “probabilistic roadmaps”-based motion planning techniques. The iFold server (http://ifold.dokhlab.org) [55] allows discrete molecular dynamics (DMD) simulations of protein dynamics using simplified two-bead per residue protein models. These tools facilitate the second phase of protein folding research, whereby targeted simulations may be performed for probing the dynamics of protein folding and unfolding under controlled conditions.

DMD approaches [43-46] with simplified structural models of proteins have been extensively used for investigating general principles of protein folding and unfolding [56-60]. Dokholyan et al.[61] have highlighted the differences between molecular dynamics and DMD approaches. As opposed to the traditional MD approach of iteratively solving Newtonian equations of motion for evolving protein folding trajectory, DMD simulations solve ballistic equations of motion with square-well approximation to inter-particle interaction potentials. DMD algorithm gains efficiency over traditional MD simulations in multiple ways. First, due to ballistic modeling of particle dynamics, a larger time step can be used in DMD simulations on average, which corresponds to the time interval between fastest ballistic interactions; secondly a faster inter-particle collision detection and velocity updating algorithm is used, since only the coordinates of colliding atoms need to be updated at each collision. Additionally, faster simulation speeds are attainable with the DMD approach through simplification of protein models. Overall, an increase in simulation speed of 5–10 orders of magnitude is attainable using DMD [62]. Jang et al. [56] used DMD and simplified protein models with Gō interactions [63] to probe protein folding kinetics. Protein folding kinetics studies using DMD simulations are reviewed in [42]. Recently, DMD simulations we used in uncovering the structural mechanisms of protein aggregation [64-66]. Among the fundamental challenges in studying protein folding using computer simulations are the time-scales and length-scales that can be investigated. DMD simulations have been shown to be useful for investigating long-timescale folding dynamics of complex biological systems such as poly-alanine aggregation [66, 67] and the nucleosome core particle [68].

In addition to the extensive in silico and in vitro studies of protein folding, significant progress has been made in understanding protein folding in vivo. There are two major differences between protein folding in vivo and in vitro. First, protein folding in vivo is usually assisted by molecular machinery, such as chaperones (in an ATP-dependent manner), and often involves small molecule cofactors. Molecular chaperones such as the heat shock protein Hsp70 and chaperonin proteins facilitate protein folding, in part, by isolating the proteins from bulk cytosol[69, 70]. Hartl and Horwich pioneered the research of chaperone-mediated protein folding[71], highlighting the differences between in vivo and in vitro folding mechanisms[69]. The mechanism of chaperonin GroEL mediated folding, including in vivo folding intermediates, has been extensively studied by Horwich and Gierasch [72, 73]. Work by Landry et al. [74] showed that chaperon binding promotes α-helix formation in partially folded polypeptide chains. Horowitz et al. have investigated the role of chaperonin Cpn60-mediated hydrophobic exposure in protein folding[75, 76]. Nearly one third of all proteins in living cells are coordinated to small molecule cofactors. The pioneering work of Wittung-Stafshede and coworkers on the role of cofactors in in vivo protein folding[77, 78] demonstrated that bound metals stabilize the native fold, suggesting cofactor binding to unfolded polypeptides dramatically accelerates folding timescales[77].

A second notable difference between in vivo and in vitro protein folding is the fact that the concentrations of macromolecular solutes in cells can reach hundreds of grams per liter in cells[79], but most in vitro studies are performed in buffered solution with <1% of the cellular macromolecule concentration. The crowding environment in vivo can have a significant impact on protein stability and the native structure by changing the energy landscape of protein folding[80, 81]. Dedmon et al.[82] showed that FlgM, a 97-residue protein from Salmonella typhimurium is unstructured in dilute solution, but in E.coli cells its C-terminal half is structured. McPhie et al.[83] found that a molten globular state of apomyoglobin at low pH is stabilized by high concentration of the inert polymer, dextran, compared to the unfolded state. Moreover, it was found that aggregate formation from human apolipoprotein C-II is significantly accelerated by the addition of dextran[84], suggesting a direct effect of molecular crowding on protein aggregation.

Experimentally, ΔG values can be obtained from denaturing experiments [6, 87-90] where the protein unfolds by increasing temperature or by adding denaturing agents such as urea and guanidinium HCl (GdHCl). Theoretically, given the interactions, free energy can be obtained from statistical mechanics using the partition function, Z, as G = -RT lnZ. Analytical calculations of partition functions require integration over all degrees of freedom in the protein’s conformational space, which is impossible in practice, except for simple models. Advances in computational biology have made possible the direct calculation of ΔΔG by MD or Monte Carlo simulations[28, 91-94], a comprehensive review of which can be found in Ref. [95]. Using MD simulations, a ΔΔG value has been calculated for T157V mutant of T4 lysozyme which is in close agreement with experimental measurements [92].

The computational cost of direct ΔΔG estimation is still too high, however, preventing it from being applied to a large number of mutations for protein engineering. More heuristic approaches have been adapted which try to describe the free energy using empirical or effective functions taking advantage of the vast amount of known protein structures and stability measurements. Such simplifications significantly decrease the computational overhead and allow ΔΔG calculations of large numbers of mutants that can be compared with experimental results. In recent years, various methods have been proposed for large-scale ΔΔG predictions with reasonable prediction accuracies (which are commonly assessed by the linear correlation coefficient between the predicted and measured ΔΔG values). A comparison between these methods is listed in Table 2.

One approach utilizes statistical potentials that are developed using information from known protein structures, where the frequency distribution of amino acids conformations (such as pairwise distances and torsion angles) is used to extract effective potentials for free energy evaluations[96]. Gilis and Rooman first applied database-derived backbone dihedral potentials to study the change of thermodynamic stability upon point mutations[97-99]. They found that torsion-angle potentials predict ΔΔG accurately for mutations of solvent-exposed residues and that distance-dependent statistical potentials are more accurate for predicting the ΔΔG of buried residues. They obtained correlation coefficients of 0.55 to 0.87 for a dataset of 238 mutations. Zhou et al. [100] developed a knowledge-based potential using the distance-scaled finite ideal-gas reference state (DFIRE) approach and calculated ΔΔG for 895 mutants, which have a correlation of 0.67 with experimental measurements. Similarly, statistical potentials utilizing side-chain rotamer libraries[101], direction- and distance-dependent distributions[102], and four-body interactions[103] were also adapted for ΔΔG predictions and significant agreement with experimental measurements was achieved.

Another approach for large-scale ΔΔG predictions uses empirical functions to describe free energy changes induced by mutations and trains the parameters to recapitulate the experimental results. Guerois et al. [104] developed the FOLD-X energy function to study the stabilities of 1088 mutants. They used a comprehensive set of parameters to describe the van der Waals, solvation, hydrogen-bonding, electrostatic and entropic contribution to the protein stability, and obtained a correlation of 0.64 for the blind test set after training their parameters on 339 mutants[104]. Khatun et al. [105] utilized contact potentials to predict ΔΔG of three sets of 303, 658 and 1356 mutants and their prediction correlations varied between 0.45 to 0.78. Bordner and Abagyan [106] used a combination of physical energy terms, statistical energy terms and a structural descriptor with weight factors scaled to experimental data for ΔΔG predictions, and found a correlation of 0.59 on 908 test mutants. Saraboji et al. classified the available thermal denaturing data on mutations according to substitution types, secondary structures and the area of solvent accessibilities, and used the average value from each category for the prediction and obtained a correlation of 0.64 [107].

Taking advantage of the vast amount of experimental ΔΔG data now available, machine-learning techniques have been introduced for ΔΔG estimation. Capriotti et al. [108, 109] trained a support vector machine using temperature, pH, mutations, nearby residues and relative solvent accessible area as input vectors. The support vector machine, when applied to a test set, gives a prediction correlation of 0.71. Cheng et al. [110] improved the support vector machine model to directly include the sequence information and obtained higher correlation prediction accuracy.

There are two significant drawbacks with training-based studies. First, improvement of the prediction accuracy relies on the available experimental stability data for parameter trainings. It is questionable whether parameters obtained from these trainings are transferable to other studies [105] since the experimentally-available mutation data may be biased (e.g., towards substitutions of large residues for small ones). Second, some mutations introduce strains in the protein backbone. To properly estimate the ΔΔG values, it is necessary to estimate the structural rearrangement that a protein undergoes to release the strain. To our knowledge, protein dynamics and flexibility have not been explicitly modeled in previous methods. Ignoring protein flexibility prohibits the application of current prediction methods to a wide range of mutations [100, 104].

To address both these caveats, Yin et al. [111] developed a novel method, Eris (http://eris.dokhlab.org), for accurate and rapid evaluation of the ΔΔG values using the recently-developed Medusa redesign suite[112]. Eris features an all-atom force field developed from x-ray crystal structures, a fast side-chain packing algorithm, and a backbone relaxation method. The ΔΔG values of 595 mutants from five proteins were calculated and compared with the experimental data from the Protherm database and other sources [104, 105, 113, 114]. Significant correlations of ≈0.75 were found between predicted and experimental ΔΔG values. Eris identifies and efficiently relaxes strains in the backbone, especially when clashes and backbone strains are introduced by a small-to-large amino acid substitution. Interestingly, when high-resolution structures are not available, Eris allows refinement of the backbone structure, which yields better prediction accuracy. Compared with other ΔΔG prediction methods, the Eris method is a unique approach that combines physical energies with efficient atomic modeling, resulting in fast and unbiased ΔΔG predictions.

Probing the structural properties of intermediate states and determining folding pathways have been major experimental challenges. Nonetheless, there have been significant advances in experimental techniques. Experimental methods with different temporal resolutions allow a more detailed dissection of the folding process (Table 3). Other methods allow for investigation of protein structure (e.g. NMR, ultraviolet/visual light CD, site-directed mutagenesis, Φ-value analysis, isotope labeling) and global properties (e.g., mass spectroscopy, quasi elastic light scattering, ultracentrifugation) [118].

As experimental data on folding rates of various proteins accumulated over the years, people sought the determinants of folding rates. Plaxco and co-workers [119] showed that there is a high correlation between the folding rate and the structural properties of proteins, as defined by contact order CO, , where L is the sequence length, N is the total number of inter-residue atomic contacts within a cutoff distance, ΔL_ij is the sequence separation of contacting residues i and j. Interestingly, assessment of other protein folding rate determinants such as local and long-range contacts were found to perform equally well as the contact order[53]. Also, contact order is a geometric property which does not take into account the distribution and strength of the interactions on the rate [120]. Fersht and coworkers found that specific interactions in the folding nucleus are equally important determinants of folding rates. For example, mutations in the folding nucleus of CI2 did not change its contact order, but result in a three order magnitude increase in the folding rate [90, 121, 122]. Sequence-based prediction of folding rates has also been proposed and was found to be of comparable performance to that of contact order [123, 124]. Thus, there is still a debate as to whether structure-based or sequence-based prediction is a more reliable predictor of folding rates [123, 124].

Extensive studies have likewise been made on the prediction of folding rates using molecular dynamics simulations. The primary bottleneck in this approach was sampling the time scales where the folding transitions are observable. Thus, early studies pioneered by Caflisch et al. [125-127] employed continuum solvent with low viscosity to observe multiple folding transitions. However, there is a nonlinear relationship between the folding time and viscosity [128], hence, the precise effect of very low viscosity on the protein folding kinetics of various systems remains unclear. To circumvent this problem, Pande et al. used “coupled ensemble dynamics” to simulate the folding of a β-hairpin from protein G using continuum solvent model [129] and united atom force-field [130] with water-like viscosity. In this and other subsequent simulations of other β-hairpins, the calculated folding rate was in close agreement with experimental measurements. A remarkable simulation in this class is a 1 μs folding simulation on the villin headpiece by Duan and Kollman [28]. Folding rate predictions of this type have been limited to small two-state proteins. As protein size increases, it is difficult to computationally study folding kinetics. Rate predictions have been likewise performed using molecular dynamics simulations using explicit water models such as the TIP3P [131] and SPC [132] to gain additional insight into folding kinetics. Examples of MD simulations using explicit solvent which yielded experimentally consistent rates were performed by Pande et al., who observed helix-coil transitions [133] and protein folding [34]. However, a potential drawback in the use of water models is that they are parameterized to a single temperature (~298 K), and thus may bias the dynamics in non-native temperature simulations. Overall, theoretical rate determinations and their subsequent comparisons with experiments provide a test of our understanding of protein folding kinetics.

Structural investigation of the transition state ensemble (TSE) is extremely challenging experimentally, since the TSE is an unstable state whose experimental detection is very difficult. Computationally, transition state conformations may be identified using unfolding simulations (cf. section on unfolded protein states), projection into one or two reaction coordinates, validation of putative transition states through calculation of the probability to fold (P_fold), and path sampling. The intuitive appeal of low-dimensional energy landscapes in explaining simple chemical reactions inspired people to develop a similar formulation for protein folding. However, unlike simple molecules, the high degrees of freedom of a protein make the analysis formidable. Thus, several groups proposed dimensional reduction by projecting the multi-dimensional energy landscape into few relevant coordinates. The proposed reaction coordinates could be the volume of the molecule [134], the fraction of amino acids in their native conformation [135], the number of contacts between amino acids [136, 137], and the fraction of native contacts in a conformation [136, 137]. Some others directly tackled the transition state by developing rigorous path sampling techniques [138]. They constructed a large ensemble of transition paths, and through statistical analysis, they determined conformations whose P_fold=0.5. Although this method is computationally expensive, it is advantageous since there is no presupposed reaction coordinate.

Characteristics of the transition state ensemble have mainly been investigated by Φ-analysis[139], which involves measuring the folding kinetics and equilibrium thermodynamics of mutants containing amino acid substitutions throughout a protein. This method provides means to identify interactions mediated by specific amino acid side chains that stabilize the folding transition state[140]. Φ-analysis has been applied to a large number of proteins (such as BPTI, myoglobin, protein A, ubiquitin, SH3 domain, and the WW domain) and, recently, even to amyloidogenesis[141, 142]. However, despite the prevalence of the methodology, there is debate regarding the validity and conventional interpretation of Φ-analysis, especially when the ΔΔG between wild type and mutant is less than 1.7 kcal/mol [2, 42, 143, 144]. However, Fersht and coworkers argued that reliable ϕ-values can be derived from mutations in suitable proteins with 0.6 < ΔΔG < 1.7 kcal/mol[145]. Plaxco and coworkers disproved the assumed independence of the changes in free energy of transition and folded states when calculating error estimates in Φ-values [143]. They proposed a new method of error estimation that accounts for the interdependence of changes in free energy of transition and folded states.

Using simplified protein models and rapid sampling DMD, Dokholyan et al. directly observed and characterized the transition state ensemble of Src homology 3 (SH3) (Fig. 2) [47, 59, 146, 147]. To probe the contribution of each amino acid residue to the transition state ensemble, they calculated the Φ-values, and found high correlation between simulation and experimental Φ-values. Moreover, they also predicted that the two most kinetically important residues in folding are L24 and G64. Both L24 and G64 are experimentally-verified to be important kinetically [148].

Experiments suggest that proteins may be kinetically trapped en route to the native state[149-151], but how do proteins avoid kinetic traps? By using a Go-model scaled to include sequence-specific interactions, Khare et al. [152] found that the residues which contribute most to Cu, Zn SOD1 stability also function as “gatekeepers” that avoid kinetic traps and protein misfolding. “Gatekeeper” residues were also identified in later computational studies of the ribosomal protein S6. Stoycheva et al. [153] and Matysiak et al.[154] show that the mutations of gatekeeper residues can alter the folding landscape of S6 and shift the balance between its folding and aggregation. All these computational studies are fully consistent with experimental observations and the existence of gatekeeper residues suggests a selective pressure on avoiding misfolding in natural proteins.

Among the traditional experimental methods, hydrogen exchange (HX)[169] is a unique tool that allows the detection and characterization of not only kinetic, but also weakly populated thermodynamic intermediates. In a typical equilibrium HX experiment, the rate at which an individual main-chain amide hydrogen exchanges with solvent deuterons is measured by NMR or mass spectrometry [170]. The exchange rates of different amide hydrogens are sensitive to local and global structural changes of proteins, and thus contain useful structural information about different protein conformations. From the amide hydrogen exchange rates measured by HX it is possible to detect and characterize weakly-populated intermediates which are inaccessible by bulk methods. Over the years, HX methods have provided considerable insight into the coarse features of intermediate state conformations for a wide variety of proteins[169]. However, HX is limited in its ability to describe the detailed structures of intermediates, which severely restricts its applications. In contrast, computational methods can reveal structural information at much higher resolutions that cannot be accessed by experiments. Also, unlike bulk experimental methods like HX, computational methods are able to study protein motions at the single-molecule level, which enables the detection of heterogeneous conformations in an ensemble of molecules.

In recent studies, Gsponer et al.[171] and Dixon et al.[172] have developed new computational approaches to incorporate HX protection factors from NMR experiments as constraints into MD simulations to detect and/or characterize the conformations of intermediate states. Using their new approach, Gsponer et al.[171] were able to define the thermodynamic folding intermediates of the bacterial immunity protein Ig7. A structural comparison between this thermodynamic intermediate and kinetic intermediate determined from other methods indicates that the kinetic and thermodynamic intermediates of Ig7 are similar. In a parallel study, Dixon et al.[172] computationally predicted and characterized a folding intermediate of the focal adhesion targeting (FAT) domain of focal adhesion kinase (FAK), which plays critical roles in cell proliferation and migration. The detected intermediate was hypothesized to expose a cryptic phosphorylation site in the FAT domain and regulate the localization of FAK. The existence of this intermediate and the predicted structural features were later experimentally confirmed [165]. Besides HX methods, the site-specific structural information obtained from other NMR-based techniques such as relaxation dispersion spectroscopy [173], has also been incorporated into MD simulations to study weakly-populated folding intermediates [174].

A major limitation of traditional experimental methods such as HX is that they can only measure the average behavior of an ensemble of molecules and often cannot distinguish individual folding and unfolding routes or intermediates in the ensemble. By contrast, single-molecule methods enable the observation of the folding and unfolding of individual molecules and play increasingly important roles in studying intermediates[175]. Two single-molecule techniques: fluorescence resonance energy transfer [176] and force spectroscopy including optical tweezers and AFM[177], have been utilized to probe intermediate states of proteins [178-180]. Direct comparisons of results from single-molecule stretching experiments by AFM and computational simulations shed light on the mechanical unfolding of proteins and how intermediates contribute to protein function.

Using AFM, Marszalek et al. [179] uncovered a force-induced unfolding intermediate of Ig domain I27 from titin, a modular protein which is responsible for muscle elasticity. This experimentally-discovered intermediate was also predicted by steered-molecular dynamics (SMD) simulations [179]. Based on the results of SMD, the authors further predicted that the hydrogen-bonding network between strand A’ and G in I27 mainly determines its mechanical stability. This prediction was verified in a mutagenesis study where Li et al. [181] showed that the point mutations on the residues participating in this hydrogen-bonding network dramatically altered the mechanical stability of I27. In a similar study of another protein, domain 10 of type III fibronectin module, which plays a pivotal role in mechanical coupling between cell surface and extracellular matrix (ECM), Gao et al. [182] predicted force-induced unfolding intermediates using SMD simulations. The discovered intermediates were considered to expose binding sites which are necessary for the assembly of ECM fibronectin fibrils. One of the computationally-predicted unfolding intermediates is in excellent agreement with the one observed in a more recent AFM stretching experiment [178].

It has long been postulated that the denatured state of proteins is composed of an ensemble of featureless random coil-like conformations. According to Flory’s random coil theory[188], the size of a random coil polymer, characterized by the radius of gyration R_g, follows a power law dependence on the length of the polymer chain, n, Rg = R₀n^ν. Here, R₀ is the scaling constant, which is a function of persistence length, and ν represents the power law scaling exponent. Flory predicted the exponent to be 0.6 and later a more accurate renormalization calculation obtained ν=0.588 [189]. Tanford et al. first confirmed this random coil scaling behavior for denatured proteins [190]. Using intrinsic viscosity measurements for 12 proteins denatured by 5-6 M GuHCl, the authors obtained a scaling exponent ν=0.67±0.09. Wilkins et al. [191] showed that the hydrodynamic radii of sets of 8 highly denatured, disulfide-free proteins follow a power law scaling with ν=0.58±0.11. Recently, Kohn et al. [192] reassessed the scaling behavior of denatured proteins using small angle x-ray scattering for 17 proteins of lengths varying from 8 to 549 residues. They found the scaling exponent to be ν=0.598±0.029. All these experimental results confirm the random coil scaling of denatured proteins. In a random coil model of the denatured state, proteins are believed to lack persistent structure both locally and globally. The distribution of end-to-end distances or radii of gyration can be fit by a Gaussian distribution. A recent computational study [193] confirmed this behavior by generating an ensemble of protein conformations whereby only steric interactions between amino acids were considered for four different proteins. A scaling exponent of 0.58±0.02 was obtained.

The random coil scaling behavior [188, 194] was originally derived for homopolymers. However, proteins are heteropolymers for which specific interactions between amino acids play an important role and determine a unique native structure. Hence, while the scaling of the sizes of denatured proteins follows the random coil scaling as shown in experiments, it does not necessarily exclude the possibility that the denatured proteins can have residual native-like structures. Under denaturing conditions, the protein can still exhibit native conformational bias and retain a certain amount of residual native structures. Mounting experimental evidence [168, 195-198] supports residual native-like structural elements in the denatured state for a variety of proteins. Using residual dipolar coupling from NMR measurements, Shortle and Ackerman [195] showed that native-like topology persists under strong denaturing conditions as high as 8 M Urea for a truncated staphylococcal nuclease. By applying quasi-elastic neutron-scattering on α-lactalbumin, Bu et al. [196] demonstrated residual helical structure and tertiary-like interactions even in the absence of disulfide bonds and under highly denaturing conditions. Similarly, by using triple-resonance NMR, native-like topology has also been observed in protein L [197].

Several theoretical and computational studies [199-202] have addressed the role of specific interactions in conformational biasing toward the native state in the denatured states. Using a simple force field with only steric and hydrogen bond interactions, Pappu et al. [200] demonstrated that denatured protein states have a strong preference for the native structure. It was suggested [199-202] that the conformational bias of native structures in the denatured state is a possible explanation of the Levinthal’s paradox [3].

To reconcile the seemingly controversial properties of denatured proteins — the random coil scaling of their sizes and the presence of residual native structures — several computational works have recently been reported [203-207]. Fitzkee and Rose [206] reproduced the random coil scaling exponent using a denatured protein model with fixed secondary structure elements for a set of proteins. Tran et al. [204, 207] constructed the protein denatured state ensemble at atomic resolution in the excluded volume limit. Jha et al. [205] built the denatured state using a statistical coil library. Both the Pappu and Sosnick groups found that the putative denatured state ensemble features transient local structures such as turns, strand, and helices. In the mean time, the dimension of the denatured state follows the experimentally-observed random-coil scaling exponent. Ding et al. [203] developed a computational method to model denatured proteins using a structure-based potential [63]. This interaction model is commonly used in thermodynamic and kinetics studies of protein folding [58, 208-210] to model amino acid interactions. This study [203] suggested that denatured proteins follow the random coil scaling sizes and retain residual secondary structures akin to those observed in native protein states. Hence, these computational works provide a conceptual reconciliation between two seemingly mutually-exclusive views of protein unfolded states.

What is the physical origin of the random coil scaling of protein size along with the seemingly contradictory persistence of local structures in denatured proteins? Previous computational studies [203-205, 207] suggest that the residual structures in the denatured state are limited to short-range elements, which only extend approximately seven to ten residues. The correlated fluctuation of residual structures diminishes quickly along the sequence and long-range contact formation is purely governed by the random diffusion of peptide chains. Hence, a coarse-graining process, which groups locally-interacting amino acids along the polypeptide chain into renormalized structural units (Fig. 4) reduces a denatured protein to a renormalized homopolymer. The renormalization process will result in an effective homopolymer which forms contacts due to chain diffusion. Protein sizes follow the renormalized power law scaling as proposed by Flory [188]: Rg=R₀(N/L)^ν=(R₀L^-ν)N^ν. Thus, we expect that the scaling exponent of denatured proteins ν is the same as for homopolymers, whose structural units are locally interacting amino acids.

Extant protein sequences are the result of a long history of evolutionary refinement establishing a set of interactions defining the native state. However, the same inherent recognition that occurs between sequences within a protein is the basis for a type of self-association termed 3-dimensional domain swapping ([214], extensively catalogued in 2002 by Liu and Eisenberg [215]). Domain swapping is an important phenomenon, taking part in both normal and disease-related processes, and is intimately tied to protein folding. Domain swapping may be viewed as a natural mechanism for dealing with instability due to evolutionary changes in the amino acid sequence [216]. For example, a mutation that rigidifies a loop connecting two parts of a protein induces strain, which can be relieved without the loss of function by “swapping” the portion of the protein on one side of the loop with the corresponding part of a similar protein. Dimerization by this mechanism has advantages including a high local concentration of enzymatic activity, since two functional proteins are joined together.

Many diseases are now associated with protein aggregation and particularly with a form of ordered aggregate called the amyloid fibrils, which, regardless of the native sequence and structure of the precursor proteins, share distinct structural characteristics. From a protein folding standpoint, the inherent properties of the polypeptide chain that allow proteins with little or no sequence or structural similarity to misfold and assemble into similar high-order structures are of vital interest. Studies of aggregate structure reveal defined characteristics such as extensive hydrogen bond networks perpendicular to the fiber axis, called a cross-β conformation[217], and an α to β transition known to occur during the oligomerization of amyloid-forming proteins with significant helical content. Evidence for domain swapping as an early step in the aggregation process has been reported for several proteins [218-220]. In several aggregation-associated diseases including familial Amyotrophic Lateral Sclerosis and the transthyretin amyloidoses, the dissociation of a protein from its multimeric native state is known to be the rate limiting step for aggregation, suggesting a method for preventing aggregation by stabilizing the native interfaces of these assemblies [221-223].

Amyloid fibrils have alternately been deemed responsible for the pathologies of their various associated diseases and, more recently, credited with delaying or counteracting the observed pathologies by acting as a sink for highly cytotoxic soluble oligomers [224]. The viewpoint that soluble oligomers act as cytotoxic species has garnered widespread attention since at least 1999 when it was noted that the abundance of soluble Aβ 1-42 oligomers is inversely correlated with neuronal degeneration in Alzheimer’s disease whereas amyloid levels do not correlate [225, 226]. In 2003, Glabe and co-workers discovered that soluble oligomeric species from several disease-related proteins shared a common structural epitope to which an antibody was developed [227]. Later studies showed that soluble oligomers are able to disrupt the polarity of cellular membranes [144, 228, 229], one possible basis for disease-associated toxicity.

Various cellular protective mechanisms have evolved to ensure the proper folding of proteins. Molecular chaperones, for example, recognize misfolded proteins and provide an environment conducive to the formation of the appropriate native contacts [230]. Vast proteosomal machinery clears proteinaceous debris from the cell using ubiquitin ligases to tag misfolded proteins for degradation and removal [231, 232]. One hypothesis formulated to explain the prevalence of protein aggregation in neurodegenerative diseases is that the ubiquitin proteasome loses efficiency over time causing a buildup of protein aggregates and debris in post-mitotic cells such as neurons [233]. Also, parkin, an E3 ubiquitin ligase was found to be mutated in at least half of autosomal recessive juvenile parkinsonism patients, suggesting that a deficit in the clearance of its target protein leads to an early onset of symptoms [234].

Computational studies of protein aggregation have traditionally been inadequate due to the massive complexity of the system in both time and length scale. Several approaches have been used to overcome this complexity (Fig. 5). Traditional all-atom molecular dynamics simulations have been carried out to model the aggregation of disease-related peptides such as Aβ and polyglutamine [235-237]. Other techniques seek to identify which parts within larger proteins are responsible for their aggregation behavior, resulting in the identification of sequence stretches in various proteins that are “amyloidogenic,” or “hot spots” for aggregation [238-240]. Underlying this work is the idea that evolution acts to prevent aggregation by burying aggregation-prone protein sequences or otherwise prohibiting their apposition in protein structures and during folding¹. To study the nature of subunit assembly and extension during aggregation, which are out of reach for all-atom molecular dynamics in both the size scale and time scale, simplified protein models are being utilized [65]. By accessing aggregation events that are out of the reach of experimentalists, computational studies of aggregation are an essential compliment to the experimental findings regarding aggregate structure and formation mechanisms.