Protein folded states are kinetic hubs

Many proteins appear to fold via a single cooperative transition from a rapidly mixing ensemble of unfolded conformations to a well defined native structure (16, 18). However, based on chemical intuition, one would expect to find many more metastable states, corresponding to the numerous favorable interactions that could form in the absence of the full native structure as well as dynamics within the native state. To reconcile these points, one typically assumes a single dominate free energy barrier that serves as the rate limiting step for folding. Other barriers are often assumed to be small relative to the thermal energy (or at least to the dominant barrier) and the equilibrium probability of any intermediate is assumed to be too small to detect.

However, in some cases modeling experimental data requires the use of at least three states (19–21) and simple toy models have shown that even three-state systems can yield apparent two-state behavior (22). Thus, it is natural to hypothesize that many systems may have more complex arrangements of metastable states (5, 23, 24) yet still exhibit apparent two-state behavior.

To test this hypothesis, we first turn to an MSM for the villin headpiece. This MSM was recently built from atomistic simulations and, by assuming stochastic jumps between its states, was shown to give quantitative agreement with experimental structures and folding rates in addition to recapitulating the raw simulation data (6). Thus, the presence of numerous metastable states in this model would be strong evidence for their actual existence and the stochastic nature of transitions between them. Indeed, with a lagtime on the order of 10ns (Fig. S1), analysis of this MSM reveals the existence of at least 500 metastable states. At least 2,000 are found for NTL9 (13). The free energy barriers between our villin states have an average height of about 5.9(+/-2.5) kT (SI Text), indicating that they are nontrivial and potentially detectable. Moreover, no single dominant barrier is apparent.

To better understand the system specific results from our all-atom models, we now consider three simple models for dynamics capable of providing insight into protein folding in general. Each of these networks has six metastable states and is depicted in Fig. 1. These models have a single folding pathway (S), parallel folding pathways (P), and a heterogeneous unfolded state (H, with multiple unfolded basins that each interconvert more rapidly with the native state than with one another) as discussed in Materials and Methods.

One may be tempted to associate the states in these models with folding nuclei (25), preorganized secondary structure (26), foldons (27), or the elements of some other model of protein folding (28). However, we simply require that they all be metastable. That is, a system within one state is more likely to stay there than to transition to a different state. Moreover, we propose that the concept of metastability unifies many of the previously proposed folding mechanisms, each of which describes some systems better than others, as all consist of basic units that are stable on some timescale.

We can now imagine monitoring stochastic transitions within each of these representative systems (or ensembles thereof) with a device that can only detect the native state. This hypothetical setup is equivalent to experiments wherein unfolded molecules are allowed to relax to an observable folded state where they are trapped to prevent unfolding and refolding. Fig. 2 shows that such an experiment yields the exponential behavior typical of an ideal two-state system. In fact, exponential fits to the data after the initial lag phase only give slight underestimates of the true mean first passage times (MFPTs) between the unfolded and folded states (Table S1). Thus, even these simple systems are qualitatively consistent with both stochastic jumps between numerous metastable states and apparent two-state behavior. This is particularly surprising for model H because it cannot be divided into a single, rapidly mixing unfolded basin separated from the native state by one dominant barrier (i.e., it is not two-state).

A kinetic perspective on our simple networks helps to explain why two-state behavior is often observed even when there are many large barriers. As discussed previously, when there is a single dominant rate then faster transitions will tend to be lost in the noise. Multiple slow rates will also be lost in the noise if they are too similar. Moreover, this same logic applies even when there are multiple folding routes from different starting points (and thus no kinetically relevant two-state decomposition). Thus, observing anything other than mainly single exponential kinetics requires a delicate balance wherein the slowest rates differ sufficiently to distinguish them but not so much that one dominates the rest, not to mention extremely precise measurements.

Fortunately, there is ample evidence that achieving this balance and the precision necessary to detect it are possible. Multiexponential behavior is often consistent with the experimental data, but fit to stretched exponentials (29, 30). Increasing the temporal resolution of single molecule pulling experiments has also steadily revealed more metastable states and kinetic measurements can be probe dependent (31, 32). We propose that the ability to simultaneously monitor multiple degrees of freedom (such as extension and FRET) in single molecule experiments would reveal even more metastable states, particularly if MSMs were used to choose the number of probes employed and their placement.

We now investigate which of the simple network topologies is most representative of real protein free energy landscapes. As a first step, we have calculated that every state can reach the native basin of our villin model in one or two steps. This eliminates the possibility of a single pathway because states with that topology could require up to 499 steps to reach the native basin.

Determining whether the parallel pathway model (17, 33, 34) or the heterogeneous unfolded state model is more representative of villin requires a definition of the unfolded state(s). Because every nonnative state can reach the native basin in one or two steps it is natural to label every state that is not directly connected to the native state (332 in all) as unfolded and all other nonnative states (167 in all) as intermediates.

Taking this definition, we can now examine the distribution of MFPTs from each unfolded state to the native state as well as the distribution of MFPTs between all pairs of unfolded states. Doing so reveals that the average MFPT to the native state is 880(+/-270) ns, in reasonable agreement with the experimentally predicted folding time of 720 ns (14). Moreover, this value is much lower than the average MFPT between pairs of unfolded states (approximately 370 ns), as shown in Fig. 3 A and B. Considering every nonnative state as part of the unfolded ensemble also gives similar distributions (Fig. S2), implying that these results are robust to the exact definition of the unfolded state. Similar results are found for NTL9 as well (Fig. S3). Thus, we can conclude that the heterogeneous unfolded state model is most representative of our villin and NTL9 models and possibly proteins in general. This result is in contrast to existing theories of protein folding, which assume rapid equilibration within the unfolded ensemble (17, 35, 36).

Examination of representative structures suggests that nonnative interactions (often in the context of relatively compact conformations) and the enormity of conformational space are responsible for slow transitions between unfolded basins (Fig. S4). Nonnative contacts can easily have free energies on the order of native contacts, making nonnative states reasonably metastable. Once a set of nonnative contacts is broken, the probability of forming a particular set of other nonnative contacts is quite small due to the large number of other possibilities. This small probability is equivalent to a slow rate. In contrast, evolutionary pressure to fold makes transitioning to the native state reasonably probable, which equates to fast folding relative to slower transitions between unfolded basins.

The tight distribution of MFPTs to the native state is also consistent with our explanation of apparent two-state behavior. Due to experimental noise, it is difficult to justify using more than one or two exponentials to fit the relaxation of our coarse-grained villin model with 500 states, as shown in Fig. 3C. Only with an extremely high signal to noise ratio can one accurately identify the deviations from single exponential relaxation shown in Fig. S5. We also note that more fine-grained models for villin can capture the burst phase in its relaxation (Fig. S6) but here we emphasize the ability of our coarse-grained model to capture the apparent two-state behavior that dominates this system’s relaxation (14).

Our ability to reconcile our model with existing experimental data on the nature of the unfolded ensemble (specifically under native conditions, as opposed to the more rapidly mixing denatured state) indicates that more experiments will be required to definitively falsify or support our conclusions. For example, Nettels et al. have reported a 50 ns global relaxation time within the unfolded ensemble (37). Our model, however, would suggest that this may be due to relaxation within individual unfolded basins, not between them. This hypothesis is consistent with recent measurements of slow dynamics in the unfolded ensemble from the Lapidus lab (38, 39). Therefore, we suggest that this may be an interesting direction for future experimental work. In addition to existing methodologies for probing the unfolded ensemble, single molecule experiments monitoring multiple degrees of freedom could help to falsify or support our conclusions.

If our heterogeneous unfolded state model is indeed generally true then protein folding kinetics cannot be accurately described by two-states separated by a single barrier. Instead, folding must be understood in terms of multiple pathways starting from a number of distinct states. Mixing between pathways adds another layer of complexity to the folding process. Modeling the effects of mutations will thus require considering changes in the relative free energies of numerous states and barrier heights. Understanding the global effects of small changes on networks will likely also be important for protein design.

The accessibility of villin’s native state implies the hub-like connectivity characteristic of small-world and scale-free networks (40, 41). We can test this hypothesis by counting the number of connections observed between states because only those transitions with probabilities above some threshold are observed with our finite sampling (all transitions would be observed with infinite sampling). Examining subsets of the states independently, one finds that the average degree (or number of connections) increases as one moves from the unfolded states to the native basin. The unfolded states have an average degree of 12 whereas the intermediate states have an average degree of 25. The native state acts as a hub, connected to 167 other states. Similar results are found for a small β-sheet peptide (42) and NTL9.

Reduced connectivity between nonnative states results in slow dynamics within the unfolded ensemble. This connectivity contradicts other models, which predict bottlenecks close to the native state and high connectivity in nonnative regions (17, 34, 36, 43), as depicted in Fig. 4A. A more thorough discussion of the similarities and differences between our model and those proposed previously is given in the next section.

The native hub explains how villin folds so quickly. Just as there are only about six degrees of separation between people in the United States (44), it is possible to reach villin’s native state in one or two jumps (each 15 ns). Therefore, it is possible to fold from anywhere in the landscape in 30 ns or less. This result is consistent with recent experimental work showing that the transition path time between the unfolded and native ensembles can be as much as four orders of magnitude faster than the average folding time (45) and likely results from evolutionary pressure to fold quickly.

Due to the kinetic proximity of the native state with a 15 ns lagtime, we see that villin can fold in just 30 ns; however, such trajectories are rare because the metastability and connectivity of nonnative states makes taking a direct route to the native state improbable. Instead, villin will often spend considerable time in a web of nonnative states before finally folding, resulting in an average folding time on the microsecond timescale. In the future, it will be interesting to test whether slower folding proteins have unfolded states further from the native one or just more strongly metastable states, which equates to higher barriers and slower transitions between states. Preliminary analysis of NTL9 suggests every basin can reach the native state in 5 steps (approximately 100 ns) or less.

We have also found a rough correlation between the connectivity of states and their equilibrium probabilities. The average probabilities of unfolded and intermediate states are approximately 0.0005 and approximately 0.004, respectively. The native state has an equilibrium probability of approximately 0.2. Fig. 4B shows a schematic of a protein folding network that attempts to capture all of these observations in a humanly comprehendible manner. All of these observations are in qualitative agreement regardless of the degree of lumping; that is, whether one uses smaller and more numerous states to capture more local minima in the landscape or fewer and more voluminous states to obtain an even more coarse-grained model. Whereas one may be tempted to consider Fig. 4B merely an alternative depiction of a funnel, we emphasize that the kinetic connectivity of the native state and lack of connectivity within the unfolded ensemble are important qualitative deviations from traditional funnel theory (17).

An important methodological consequence of the network topology found here is that many short, parallel simulations (or experiments) started from arbitrary initial points are an excellent way of exploring the entire free energy landscape. In the extreme case of using a single starting point, one could still reach every free energy basin despite the presence of numerous metastable states so long as each simulation was longer than the diameter of the network (the minimal time that allows one to reach any state from an arbitrary starting point). However, reaching every state would be impossible with simulations that were shorter than the diameter of the network. Thus, our network theory provides an alternate explanation for the previously noted need to have simulations longer than some minimal lag phase, which was then attributed to the need to equilibrate within the unfolded state before folding in two-state systems (46).

Another simple but more efficient strategy would be to start simulations from multiple conformations dispersed throughout phase space and run them long enough to ensure mixing between them and coverage of the entire space. In fact, Fig. 5 and Fig. S7 show that such a scheme is actually more valuable than a few long trajectories, using a relative entropy metric for MSMs from ref. 47 to measure the information content of different datasets relative to our validated villin model. However, this trend can be seen to break down for simulations that are insufficiently long or too few as they are unlikely to reach every state or traverse every possible pathway between pairs of states. The simulation length at which this breakdown occurs decreases as the number of simulations increases though. Even better performance can be obtained using adaptive sampling algorithms (47, 48), which direct sampling to where it is needed most to improve a model.

There is a long history of theoretical models for protein folding (28) so it is important to put our work in the context of these previous theoretical approaches. In particular, folding funnel models (17, 36, 43) have dominated much of how the field currently conceptualizes protein folding and hence it is natural to compare our model to such theories. One of the most similar funnel categories is type0B, which is characterized by overall downhill folding interrupted by a glass transition along the reaction coordinate (17). Whereas this regime does include slow dynamics between compact states, it also results in a small number of folding pathways relative to higher connectivity in the unfolded ensemble. In addition, this and other previous funnel-based models have explicitly described rapidly interconverting unfolded states, as reflected in the “bottleneck” discussed in previous works (34, 35), as well as the choice of structurally based reaction coordinates like the number of native contacts (Q) (17, 35), which directly requires that dynamics along orthogonal degrees of freedom, such as interconversion between unfolded conformations, is rapid compared to folding. In contrast, we find a large number of folding pathways, slow dynamics between unfolded states relative to folding, and no glass transition. Our folding rates are also quite similar, rather than the different rates characteristic of the folding pathways in type0B folding.

Other funnel models have recognized the possibility of a large number of folding pathways (17, 33, 43), but still in the context of fast dynamics within the unfolded basin relative to slower transitions to the folded state. Some have even gone so far as to assume global connectivity (49, 50); however, even these emphasize that local connectivity would dominate in the full dimensional conformational space and global connectivity only arises when projecting onto a few order parameters. Furthermore, they argue global connectivity will not give an activation barrier and, therefore, these models are primarily intended for studies of downhill folding or the early activationless stages of folding. Our model, on the other hand, has a native hub and slow dynamics in the unfolded state relative to faster folding regardless of the degree of coarse-graining one employs. We also demonstrate that this can result in apparent two-state folding (i.e., activated kinetics) and that this occurs in nondownhill folding proteins, such as the millisecond folding NTL9.

Ref. 6 describes the use of the MSMBuilder package (https://simtk.org/home/msmbuilder/) (23) to construct an MSM with 10,000 microstates for the villin headpiece (HP-35 NleNle). This model was based on approximately 450 all-atom, explicit solvent simulations, each up to 2 μs in length, for a total simulation time of 354 μs (51). Whereas the longest timescale transitions in the model from ref. 6 were found to be Markovian, implying memoryless transitions between metastable states, not every state was metastable. We used MSMBuilder to lump kinetically related microstates into 500 metastable macrostates to ensure a direct correlation between states in the MSM and free energy basins, as described in SI Text. This is equivalent to common experimental analyses in which the potential is smoothed and the friction is rescaled. We note, however, that the free energy landscape for this system is actually a hierarchy of basins so it is possible to build many valid MSMs with different numbers of states. As a result, one would not expect there to be exactly 500 experimentally detectable states. Regardless of the resolution at which one examines this hierarchy, however, requiring that each state is metastable ensures that they are directly related to a free energy basin. Thus, our networks of metastable states are an important step beyond previously described networks, which often used simpler approximations to define state boundaries and the transition rates between states (17, 34, 42, 52, 53). An additional 40,000 simulations, each up to 400 ns in length (for a total simulation time of 14 milliseconds), were also assigned to this MSM to explore the effect of using more simulations.

Preliminary results for a 39-residue fragment of NTL9 are based on an MSM built from approximately 1.5 milliseconds of simulation in implicit solvent with a different force field (13). Similarities between these two systems thus suggest our results are not a force field artifact.

We have designed three simple networks, depicted in Fig. 1, that capture the essence of various protein folding mechanisms. Each of these models has six metastable states with approximately the same equilibrium and transition probabilities so that differences between their behaviors may be attributed to differences in their topologies (SI Text).

The first model (S) has a single folding pathway. This model is a natural extension of the common U↔I↔N model (19, 54) and is often used to justify the expense of running long simulations as shorter ones could fail to reach every state.

The second model (P) has parallel folding pathways. Parallel folding pathways have been proposed for a number of systems (20, 21, 33, 51). In addition, this model emphasizes the need to observe numerous folding and unfolding transitions to obtain sufficient statistics on the entire process. The increased connectivity relative to S also results in faster timescales.

The third model (H) has a heterogeneous unfolded state—multiple unfolded basins that each interconvert more rapidly with the native state than with one another. Thus, there is no kinetic decomposition of this model into two states, one folded and one unfolded. This model was inspired by a growing body of work on the presence of deep minima and gutters in unfolded regions of conformational space (38, 39, 55–57).