Energy landscape and multiroute folding of topologically complex proteins adenylate kinase and 2ouf-knot

Accurate yet Concise Modelling of Folding.

Previously, we developed an AICG using a multiscale protocol (30). AICG uses one-bead-per-residue representation, but its interresidue nonlocal contact energies were derived from atomic interactions. Here, we extended AICG by introducing a flexible local potential (AICG2). The local interactions are described by combination of generic statistical potentials of the residue types and local contact potentials that are based on the native structure. The relative weight and average strength of the local contact potentials were derived by a multiscale protocol, the same as for the nonlocal contact interactions. Incorporating sequence-specificity and the locally flexible nature of chains, AICG2 can reproduce folding energy landscape better than the simple Go model. As an illustration, the asymmetric folding pathways of the topologically symmetric proteins, protein G and protein L, were well reproduced (SI Appendix, Figs. S1 and S2).

High Free Energy Barrier and Substates in the AKE Folding.

AKE is a 214-residue kinase that plays important roles in maintaining the balance of ATP and ADP in the cell and, more importantly here, it has been used as a model multidomain protein (31). AKE consists of three domains, i.e., LID domain, NMP domain, and CORE domain (31) (Fig. 1A). AKE folds relatively slowly and a recent single molecule experiment revealed complex multiroute folding (24).

We first tried to realize reversible folding and unfolding transitions using AICG2 model and unbiased molecular dynamics (MD) simulations, which turned out to be difficult due to its intrinsically slow folding. Then, we used the modified multicanonical method (32) for better conformational sampling. The converged one-dimensional free energy profiles near the folding temperature T_F (Fig. 1B) show relatively high and broad free energy barrier, consistent with slow folding in vitro (33) and in silico. Albeit no apparent intermediate states in Fig. 1B, the two shoulders labelled by the arrows suggest that the folding/unfolding is not a simple two-state transition (see SI Appendix, Fig. S3). Some substates, virtually hidden in Fig. 1B, may exist between denatured and native states.

Because the high temperature simulations tend to hide the possible substates, we next conducted folding simulations starting from fully unfolded structures at the three temperatures, 0.90T_F, 0.92T_F, and 0.94T_F. A representative time course of formation of native contacts, i.e., so-called Q-score, at 0.90T_F (Fig. 1C) shows successive transitions among different substates. Particularly, the CORE domain encounters at least five substates beside the fully unfolded and native states. Such transitions among multiple substates closely resemble smFRET trajectories observed (24), although directly mapping the substates by the experiment with those by the simulation is difficult due to the lack of structural information in the experiment.

To identify these substates robustly and statistically, we plotted the kinetic free energy landscape on the reaction coordinates Q_tot and Q_CORE in Fig. 1D using 200 folding trajectories at 0.90T_F. Such a nonequilibrium free energy landscape has been used in characterizing the folding/unfolding pathways of proteins and RNA (34–36). In addition to the fully unfolded and native states, seven other substates can be identified. Qualitatively similar results can be obtained when a different reaction coordinate was used (SI Appendix, Fig. S4). In comparison, the corresponding kinetic free energy landscape obtained by the simple Go model was much smoother (SI Appendix, Fig. S5), suggesting that sequence specificity and realistic local flexibility included in the AICG2 model are responsible for making these substates in AKE.

Structural Features of the Substates in the AKE Folding.

An advantage of the MD simulation is the ability to provide detailed structural information for each of the substates. Fig. 2A shows, for each substate, the average contact scores at every residue by color; continuously from blue (native-like) to red (unfolded). The corresponding contact probability maps are shown in SI Appendix, Fig. S6. Each substate has a well-defined structural feature, which further rationalizes the identified substates. For example, in substates 1, 2, 3, and 4, the LID domain is folded. These substates differ in the order of the CORE domain or the NMP domain. Because of its depth, the state 1 corresponds to the denatured state. Whereas, in substates 5, 6, and 7, the LID domain is unfolded. Interestingly, among the seven substates, the CORE domain has five distinct structural features, consistent with the observations in Fig. 1D. We note that the number of the substates found is less than the product of the number of states of individual domains, suggesting coupling between domains. We also calculated the distributions of the distance between residues 73 and 203, which correspond to the label sites in the single molecule experiment, for each of the substates. Both the peak positions and the distribution widths largely vary among the substates, and resemble the results of the single molecule experiment (SI Appendix, Fig. S7).

Multiroute Folding of AKE and Its Modulation by Temperature.

To characterize folding routes, we analyzed the transitions among the above defined substates, identifying five most probable folding routes (arrows of each color in Fig. 2 represent one route). In the figure, the breadth of the arrows represents the probabilities of the routes. (see also SI Appendix, Table S2 and Fig. S8 for more details). At 0.90T_F (Fig. 2A), the dominant folding route is U → 1(D) → 3 → 4 → N (green arrows). In this route, the LID domain folds first (or it is folded in the denatured state). Then the central β sheet of the CORE domain folds (substate 3), which is followed by folding of the NMP domain (substate 4). At last, N- and C-terminal α helices of the CORE domain dock to make the native state. There was a slightly less probable alternative route, U → 1(D) → 3 → 6 → 7 → N (black arrows), where the LID domain, which was folded in 1(D) and substate 3, transiently unfolded in substate 6, and 7.

This second route became dominant at higher temperatures, 0.92T_F and 0.94T_F (Fig. 2 B and C). Thus, the AKE folding took heterogeneous routes, and the populations of the folding routes are sensitive to temperatures. These results are in agreement with the single molecule observations in which a few folding routes were identified and were found sensitive to the denaturant concentration. We note that although we can approximately model the effect of denaturant by high temperatures, the actual temperature dependence of folding can be very complicated as shown in a theoretical work (37). Both the equilibrium simulations around the folding temperature and the folding simulations at low temperatures showed that the LID domain frequently transits between the folded and unfolded states during the whole folding trajectory (Fig. 1C, SI Appendix, Fig. S9). As higher temperature tends to destabilize the LID domain (SI Appendix, Fig. S10), increasing the temperature promotes the folding route U → 1(D) → 3 → 6 → 7 → N as observed in Fig. 2. Interestingly, such multiroute folding was also observed in some simple proteins (38), and has been generally linked to topological complexity of the folded state and the structure symmetry (39).

Slow Folding of Knotted Protein 2ouf-knot.

We move to a protein with another type of topological complexity, namely, knot formation. Here we focus on a designed knotted protein, 2ouf-knot, and its unknotted counterpart, 2ouf-ds, which have almost identical sequences and tertiary structures, thus being ideal model systems to study the role of a knot in folding (Fig. 3 A and B). These proteins were designed based on a natural dimeric protein, in which each subunit has four helices. Connecting two subunits with a nine-residue loop (L1 in Fig. 3A) led to 2ouf-knot, with left-handed trefoil knot of 28-residue depth in N terminus and 57-residue depth in C terminus. In contrast, the two subunits were connected with a disulfide bond to make 2ouf-ds that lacks a knot. Experimentally, both proteins reversibly folded, but 2ouf-knot folded an order of magnitude slower than 2ouf-ds (26).

Here, using AICG2, we conducted extensive unbiased MD simulations at a wide temperature range for both 2ouf-knot and 2ouf-ds. In simulations, both proteins exhibited reversible folding and unfolding transitions (Fig. 3 C and D) and we confirmed that the knot is indeed formed for 2ouf-knot. Consistent with the experiment, the folding of the knotted one was much slower than the unknotted counterpart (Its folding rate is nearly seven times slower than that of the 2ouf-ds). The higher free energy barrier of 2ouf-knot in the one-dimensional free energy profile (Fig. 3E) also indicates slower folding of the knotted protein than that of 2ouf-ds. Apparently, such slowing of 2ouf-knot folding purely results from the knot because it is the only major difference between 2ouf-knot and 2ouf-ds.

Interaction Heterogeneity and Chain Flexibility Promote the Knot Formation.

The reversible folding/unfolding transitions in Fig. 3C surprised us because previous folding simulations of another knotted protein with a Go model showed quite a low ratio of successful knot formation (29). To quantify the ratio of successful folding, we conducted extensive independent folding simulations in a range of temperatures below T_F. The successful folding ratio as a function of temperature showed a bell-like shape with the peak of 96% near 0.95T_F (Fig. 3F). In comparison with the previous study, is this high ratio due to interaction features in AICG2 or due to difference in topology? To this end, we conducted additional folding simulations with three different models, namely, a simple Go model and two intermediate models (Fig. 3F; see SI Appendix for the details of each model). Consistent with the previous study (29), the simple Go model led to quite low successful folding ratios (approximately 30%). For majority of the unsuccessful events, the protein stayed at highly trapped structures (SI Appendix, Fig. S11), in which the segment H6-H7 is misplaced around the H3, without forming the knot. Apparently, knotting from this trapped structure requires the detachment of H8 from the H2 and H5, and the threading of the longer C-terminal segment H7-H8 through the H2-H3-H4-L1 loop, which is difficult. On the contrary, for all the successful folding events by simple Go model and AICG2, the protein formed a knot by threading the segment H1 through the H4-L1-H5-H6 loop, which is much easier due to the short length of H1. A chimera model where local interactions are simple Go and nonlocal interactions are AICG (designated AICG) gave an intermediate ratio of successful folding. Similarly, another chimera model where local interactions are of AICG2 type but with homogeneous strength and nonlocal interactions are from simple Go model (designated FLP) also showed intermediate ratio. Thus, we can conclude that both the sequence specificity of contact interactions and the local flexibility in AICG2 contributed to the ability of avoiding the highly trapped structure during folding, although it is not easy to trace such ability to certain specific contacts (SI Appendix, Fig. S12). Actually, a previous study by Onuchic et al. showed that by selectively enhancing certain contact interactions, the successful ratio can be improved to a large extent (29), indicating the importance of the sequence specific interactions. Yet, it is of some surprise because 2ouf-knot is a designed protein and has not evolved to achieve folding ability.

Folding of the 2ouf-knot Is More Topologically Frustrated Than Its Unknotted Counterpart.

To compare folding mechanisms of knotted and unknotted proteins, we next drew a number of two-dimensional free energy surfaces (Fig. 4, SI Appendix, Figs. S13 and S14). In Fig. 4, although the overall shapes of the free energy surfaces of two proteins show some similarity, the knotted one possesses some unique features. In early stages of folding (Q_tot ∼ 0.2 and Q_tot ∼ 0.5), the H1-H4 contacts can be formed partly (Fig. 4A). However, because the formation of such H1-H4 contacts hinders the knot formation, they need to be transiently broken up (Q_tot ∼ 0.35 and Q_tot ∼ 0.6). Thus, the knot formation introduces topological frustrations, which was also observed in the folding of other proteins with complex topology (32). In contrast, the free energy landscape of 2ouf-ds is much smoother (Fig. 4B). The folding order of the H4-H6 contacts also has large differences. For the 2ouf-knot, the H4-H6 contacts are predominantly formed early on (Fig. 4C), whereas for the 2ouf-ds, the H4-H6 contacts can be formed in a later stage with certain probability (Fig. 4D). In addition, for the 2ouf-ds, the folding of the two halves is symmetric due to the symmetry in sequences and the tertiary structures (SI Appendix, Fig. S14). In comparison, for the 2ouf-knot, the pathway with the H1-H4 contacts form before the H5-H8 contact formation is much more probable. Such asymmetry in folding routes is crucial for the efficient folding of the knot, which will be discussed in more detail below.

Structural Features of the Substates in 2ouf-knot Folding.

Similar to the case of AKE, we conducted extensive folding simulations of 2ouf-knot starting from fully unfolded structures at the three temperatures, 0.90T_F, 0.925T_F, and 0.95T_F. Fig. 5 A and B shows two representative trajectories. The stepwise folding of different parts of the 2ouf-knot implies the existence of several substates. We plotted the kinetic free energy landscape projected onto the reaction coordinate Q_tot and rmsd in Fig. 5C. Six substates besides the unfolded state and the native state were observed clearly. The structural features of these substates are given in Fig. 6A by the average contact scores at every residue. See SI Appendix, Fig. S15 for more details. In substate 1, the knotting loop consisting of H4, L1, H5, and H6 is formed by H4-H6 or H3-H6 contacts. In substate 2, the contacts between H2 and H6 are also formed. From substate 2, H1 may pack with H4 without correctly threading the knotting loop, leading to the substate 3 (the knot not yet formed completely). The knot formation is complete at substate 4 when H1 correctly threads the knotting loop, as featured by the high contact probability between H2 and L1. In this substate, H7 and H8 are not yet packed to the final positions. Opposite to substate 4, in the substate 6, H7 and H8 are well packed, while H1 is not correctly threaded through L1. The substate 5 corresponds to the highly trapped structure shown in SI Appendix, Fig. S11 C and D.

Multiroute Folding of 2ouf-knot.

In the same way as in AKE, we analyzed the folding routes of 2ouf-knot at 0.90T_F (Fig. 6A), 0.925T_F (Fig. 6B), and 0.95T_F (Fig. 6C). More detailed analysis is given in SI Appendix, Table S3 and Fig. S16. Not surprisingly, multiple routes can populate with high probabilities, suggesting the heterogeneity of the 2ouf-knot folding routes. At 0.90T_F, the dominant folding route U → 1 → 2 → 3 → 4 → N is the most sequential one. The major feature of this route is that the knot forms before the packing of the H8 to the H2 and H5 (Fig. 5A, Movie S1). In this case, the L1 and H5, parts of the knotting loop, can have high flexibility during the knot formation. Such high flexibility of the L1 and H5 tends to expand and tilt the knotting loop, which could make the threading of H1 much easier. Thus, the knot formation is not just the threading of H1 itself, but it is accompanied with high flexibility of L1 and H5. We thus call it “concerted threading.” Detailed analysis shows that in all the trajectories with this folding route at the temperature of 0.90T_F, the knot forms by the concerted threading mechanism. In comparison, if the pathway involves the substate 6, the knot mostly forms via a slipknot intermediate (Fig. 5B, Movie S2). For example, among the trajectories with the folding route U → 1 → 2 → 6 → N, which is the most probable route via substate 6, approximately 80% of the trajectories follow the “slipknot intermediate” knotting mechanism. In this case, the H8 is well packed to H5 before the knot formation (substate 6), which makes the L1 and H5 well fixed. The rigid knotting loop makes the direct threading of H1 difficult. Consequently, the formation of the knot by the slipknot intermediate is preferred (Fig. 5B, Movie S2). Such results suggest that the knotting mechanism depends on the folding order of other parts of the protein. Notably, no route leading to the folded state via the substate 5, which requires the threading of the much longer C terminus (H8 and H9) through a different knotting loop (H2-H3-H4-L1), is observed. We found that most routes involving substate 5 terminate at the substate 5, which is a dead-end trap.

In Fig. 6, we also observe some shortcuts of the folding routes, e.g., U → 1 → 2 → N in which the placements of the N-terminal segment (H1) and the C-terminal segments (H7 and H8) are cooperative. These shortcut routes become more probable at higher temperatures, perhaps because higher temperature usually tends to make the folding transition more cooperative.

Folding of Two Other Knotted Proteins.

Because folding mechanisms of natural proteins and designed proteins could be very different as clarified in ref. 40, we also conducted the folding simulations for two natural proteins with trefoil knot, i.e., MJ0366 and YibK (28, 29, 41), which we briefly describe here. The trefoil knot of YibK is much deeper than those of the 2ouf-knot and MJ0366 (SI Appendix, Figs. S17–S19). For the MJ0366, we found that the maximal ratio of successful folding with AICG2 was 99%, a much higher ratio than that with the simple Go model (50%). However, the maximal ratios of successful folding for YibK by both simple Go model and AICG2 were very low (0.1–0.2%), which may suggest that in addition to the interaction ingredients included in the current AICG2, more realistic interaction features, e.g., nonnative interactions (28, 40), need to be considered in order to achieve efficient folding of YibK.

Reaction Coordinates and Data Analysis.

In describing the folding of AKE, the Q_tot, Q_CORE, Q_LID, and Q_NMP, defined as the fraction of native contacts for the full-length protein, the CORE domain, the LID domain, and the NMP domain, respectively, were used as reaction coordinates. Similarly, in describing the folding of 2ouf-knot and its unknotted counterpart, besides Q_tot, we used Q_H1-H4, Q_H2-L1, Q_H2-H6, Q_H2-H8, Q_H3-H6, and Q_H4-H6, defined as the fraction of native contacts between the corresponding segments. We also used the rmsd of C_α atoms. The Koniaris–Muthukumar–Taylor algorithm was used in detecting the knot (42). In constructing the folding routes, the snapshots of the folding trajectories were assigned to the substates according to the values of the reaction coordinates. See SI Appendix for details.

Simulation Details.

We used CafeMol (43) developed in our laboratory for general CG simulations (44). The all-atom MD simulations were conducted by AMBER 11 (45). The protein structures were visualized by VMD software (46). See SI Appendix for details.