Exploring the Folding Mechanism of Small Proteins GB1 and LB1

Article pubs.acs.org/JCTC

Cite This: J. Chem. Theory Comput. 2019, 15, 3432−3449

Exploring the Folding Mechanism of Small Proteins GB1 and LB1 Qianyi Cheng,†,∥ InSuk Joung,‡,∥ Juyong Lee,‡ Kunihiro Kuwajima,¶,∥ and Jooyoung Lee*,§,∥ †

Department of Chemistry, University of Memphis, Memphis, Tennessee 38152, United States Department of Chemistry, Kangwon National University, Chuncheon 24341, South Korea ¶ Department of Physics, University of Tokyo, Tokyo 113-0033, Japan § Center for In Silico Protein Science, Korea Institute for Advanced Study, Seoul 02455, South Korea ∥ School of Computational Sciences, Korea Institute for Advanced Study, Seoul 02455, South Korea Downloaded via GUILFORD COLG on August 1, 2019 at 05:49:26 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

‡

S Supporting Information *

ABSTRACT: The computational atomistic description of the folding reactions of the B1 domains, GB1 and LB1, of protein G and protein L, respectively, is an important challenge in current protein folding studies. Although the two proteins have overall very similar backbone structures (β-hairpin−α-helix−β-hairpin), their apparent folding behaviors observed experimentally were remarkably different. LB1 folds in a two-state manner with the single-exponential kinetics, whereas GB1 folds in a more complex manner with an early stage intermediate that may exist on the folding pathway. Here, we used a new method of all-atom molecular dynamics simulations to investigate the folding mechanisms of GB1 and LB1. With the Lorentzian energy term derived from the native structure, we successfully observed frequent folding and unfolding events in the simulations at a high temperature (414 K for GB1 or 393 K for LB1) for both the proteins. Three and two transition-state structures were predicted for the GB1 and LB1 folding, respectively, at the high temperature. Two of the three transition-state structures of GB1 have a better formed second β-hairpin. One of the LB1 transition states has a better formed first hairpin, and the other has both hairpins equally formed. The structural features of these transition states are in good agreement with experimental transition-state analysis. At 300 K, more complex folding processes were observed in the simulations for both the proteins. Several intermediate structures were predicted for the two proteins, which led to the conclusion that both the proteins folded through similar mechanisms. However, the intermediate state accumulated in a sufficient amount only in the GB1 folding, which led to the double-exponential feature of its folding kinetics. On the other hand, the LB1 folding kinetics were well fitted by a single-exponential function. These results are fully consistent with those previously observed experimentally.

■

INTRODUCTION The B1 domains of protein G (GB1) and protein L (LB1), as classic models of folding, have been subjected to numerous experimental and computational studies.1−20 They both are small in size, bind to immunoglobulin G, and possess similar backbone topologies, composed of a hairpin−helix−hairpin structure. In GB1, residues 2−19 form the N-terminal β-hairpin, residues 23−36 form an α-helix, and residues 42−55 form the Cterminal β-hairpin, while in LB1, residues 4−24 form the Nterminal β-hairpin, residues 26−44 form an α-helix, and residues 47−62 form the C-terminal β-hairpin (Figure 1). However, their sequence identity is low (only about 15%), and their folding behaviors are quite different.21,22 Therefore, the computational atomistic description of the folding behaviors of the two proteins is an important challenge in current protein folding studies. Early kinetic studies of the GB1 folding, using stopped-flow fluorescence methods under various solvent conditions, showed that the folding of GB1 was described by a sequential three-state model U ⇌ I ⇌ N, where the intermediate state appeared at an early stage of folding with a largely solvent-shielded tryptophan.23 Later, the combined continuous- and stoppedflow refolding experiments of GB1 showed that its folding © 2019 American Chemical Society

process could be better described by a sum of two exponentials, rather than a single exponential,24 further indicating the presence of the intermediate state. The rate-limiting step was predicted to be between the intermediate and native states. Very recently, the folding of GB1 was reexamined via a variety of probes,3 with each probe yielding different kinetics. Thus, the GB1 folding is a complex process. On the other hand, both kinetics and thermodynamics of the LB1 folding were well fitted to a simple two-state mechanism between the native and the unfolded states without accumulation of intermediates.22 It should also be noted that LB1 folds considerably slower than GB1.22 The mutational ϕ-value analyses also predicted a difference in the folding pathway between GB1 and LB1. The second βhairpin of GB1 was formed early in the folding transition state ensemble,25 whereas the first β-hairpin formed first in LB1.26,27 It has been speculated that this difference is due to differences in side-chain interactions between the two proteins; the side-chain interactions are stronger in the N-terminal than in the CReceived: November 16, 2018 Published: May 9, 2019 3432

DOI: 10.1021/acs.jctc.8b01163 J. Chem. Theory Comput. 2019, 15, 3432−3449

Article

Journal of Chemical Theory and Computation

Figure 1. Residue−residue contact maps of the native structures of GB1 and LB1, where the upper and lower triangles represent the contacts in GB1 and LB1, respectively. In parts a and b, two residues are defined as contacting each other when the distance between their Cα atoms is less than 6 Å. (a) Whole Cα contact map between different residues and (b) the two residues are separated by at least four residues along the primary sequence. In parts c and d, two residues are defined as contacting each other when two heavy atoms each from a different residue are within 4 Å. (c) Whole heavy-atom contact map between different residues and (d) the two residues are separated by at least four residues along the primary sequence.

terminal region of LB1, whereas the C-terminal hairpin is slightly shorter and has stronger interactions than the N-terminal hairpin in GB1.17 Computational studies have also indicated that GB1 folds in a much more complicated manner than LB1,17,19 and a complex pathway in folding of GB1 has been explored by simulation and experiment.3 From an all-atom Monte Carlo simulation, multiple intermediate states have been identified in the folding process of GB1,14 indicating that multiple pathways may coexist in a parallel manner. However, the atomistic details of the folding mechanisms of GB1 and LB1 have not yet been fully understood, and it is a challenging task to reproduce the experimentally observed folding kinetics of the two proteins by an all-atom molecular dynamics (MD) simulation. During the last 2 decades, the MD simulation has often been used in computational studies of protein folding.28,29 However, conventional MD simulation methods for protein folding need to be rather long (on the order of microseconds) to observe a single folding event even when the force field is properly chosen. Therefore, it has been practically impossible to observe multiple folding events by the simulations and to compare the accumulated trajectories of folding with the real experimental folding kinetics of bulk protein molecules. Because of this difficulty, various methods have been developed to enhance the capability of the MD simulation, and these methods may include all-atom unfolding simulations at a high temperature14,15,30−32

and the use of a biased potential energy for conformational sampling.33,34 Only very recently, advanced computation resources as well as more suitable force field development enabled atomic-level simulations on protein folding on the milliseconds time scale.35,36 To investigate the folding mechanisms of GB1 and LB1 by the MD simulations and to compare the simulation results with the real experimental folding kinetics of the proteins previously studied, here, we apply a new all-atom structure-based method combined with the simulation. This approach uses Lorentzian attractive interactions derived from the native structures of the proteins and has been successfully applied for investigating the conformational transitions of adenylate kinase (ADK).37 A big advantage of this approach is that the Lorentzian attractive interactions are intrinsically bound, and hence, the approach better utilizes the native structure information in facilitating the simulations. In the case of ADK, using the structure information on the open and closed states, more than 1000 conformational transitions between the two states were realized in a total of a 6μs MD simulation. In this study, the native structures of GB1 and LB1 have been used for determining the Lorentzian energy term, and this guided the proteins to fold much faster than in the conventional MD simulations. Furthermore, this method can efficiently and adequately illustrate the folding processes and help us identify key structures in the processes as well as profile their potential-energy surfaces. We observed accumulated 3433


Article


end, both the proteins folded to the nativelike structures. Detailed information about how to determine the above weight parameters for the van der Waals and the Lorentzain restraint energy terms as well as the temperatures for the hightemperature transition state analysis are given in the Supporting Information (Figures S1−S6). All simulations, started from fully stretched structures, were performed at 300 and 414 K for GB1 and 300 and 393 K for LB1. The Berendsen thermostat43 and the Beeman integrator44 were used to perform constant temperature simulations. The integration time step was 2 fs, and snapshots were recorded at each 10 ps. With different initial velocities, 48 independent MD simulations were performed for 400 ns, which amounted to a total 19.2 μs of MD simulations for both the proteins at 300 K. Similarly, 24 independent MD simulations were carried out for 300 ns at the higher temperature (414 K for GB1 and 393 K for LB1), which resulted in a 7.2-μs MD simulation in total. It should be noted that the timescale reported here (in units of nanoseconds) is different from the real experimental timescale. Clustering Analysis of the Transition-State Ensembles at the High Temperature. Using the TM-score,45 which ranged from 0 to 1 during folding from the unfolded to the native structure, as a reaction coordinate, we constructed the profiles of the potential of the mean force at 414 K for GB1 and 393 K for LB1. On each profile, there were two apparent basins with a large population of the structures located at TM-scores of 0.25 and 0.70. The transition state (TS) connecting the two basins was defined as the structure that has a TM-score higher than 0.40 and lower than 0.55. From the 24 trajectories, 5% of the frames were randomly selected, and among the structures in the selected frames, 338 and 178 structures in GB1 and LB1, respectively, fell in the TS TM-score range of 0.40−0.55. A weighted network was then generated for each protein from the conformational ensembles. Each frame was a node, and two nodes were considered as connected if the RMSD between them was less than 4.4 Å, with the inverse values of their RMSDs taken as the edge weights. The RMSD is considered to be a good metric, as it is more sensitive than the TM-score in capturing the local structure difference, such as β-strand orientations. Using our in-house community/module detection method ModCSA,46 the structures in the GB1 TS network were partitioned into three groups, while the structures in the LB1 TS network were partitioned into two groups. Clustering Analysis of All Possible States at 300 K. Simulations at 300 K were quite different from those at the high temperatures. Although the TM-score oscillated frequently between the folded state and the unfolded state at 414 K for GB1 and 393 K for LB1, all trajectories soon reached a state with a high TM-score (close to 1.0) and stayed in the nativelike state for both GB1 and LB1 at 300 K. For each simulation, if 1000 continuous steps with a 10 ps interval reached a TM-score higher than 0.65 and an RMSD lower than 2.75 Å, then those 1000 frames along with all of the frames preceding these 1000 frames will be kept for structural analysis (i.e., folding phase). There were 9265 frames for GB1 and 5060 frames for LB1 for the folding phase, and they covered a TM-score range of 0.10− 0.90. We randomly selected 10% of these frames for the clustering analysis. Weighted networks were then constructed with the edges defined as the inverse RMSDs between two frames only when the RMSD was smaller than 4.0 Å for GB1 and 5.0 Å for LB1, respectively. Exponential Fitting of the Simulation Trajectories. For both the proteins, the folding kinetics were obtained by

misfolded structures in the folding process of GB1 but not in the process of LB1, demonstrating that GB1 is a non-two-state folder but that LB1 is apparently a two-state folder.

■

METHODS All-Atom MD Simulations with the Lorentzian Restraint Energy Term. Using the Tinker molecular modeling software,38 supplemented with Lorentzian structure-based potentials and the CHARMM22/CMAP force field39 without hydrogen atoms, all-atom MD simulations were conducted for both GB1 and LB1. As shown in our previous study,37 where we applied this approach to investigate the conformational transitions of ADK, the potential energy was defined as E = Estereochemistry + wvdwEvdw + wrestraintErestraint (1) where Estereochemistry denotes the stereochemistry term which we adopted from Modeler,40 and Evdw denotes the van der Waals potential.41 Here, the structure-based Lorentzian attractive (restraint) potential was formulated as

∑

Erestraint =

l K (rij)

restraints ∈ K

∑

=

restraints ∈ K

(rij − rijK )2 1 σijK (rij − rijK )2 + (σijK )2

(2)

where the atom−pair distances (rKij ) between atoms i and j in the native structures (K) were extracted as the reference distance from the PDB structures of 1PGB (residues 1−56) and 1HZ6 (residues 4−64) for GB1 and LB1, respectively. The intrinsic variability of the distance restraint (σKij ) were generated by a statistical prediction model, Sigma-RF,42 which is based on a set of known sequence−structure alignment information. Sigma-RF predicts the intrinsic variability of inter-residue distances of proteins. To validate that the use of Sigma-RF values leads to a more accurate description of dynamics of the proteins, we compared the correlation coefficients (CCs) values between the B-factors obtained from the simulated folding trajectories and those experimentally from X-ray crystallographic analyses (PDB codes, 1PGB and 1HZ6 for GB1 and LB1, respectively). The correlation coefficients were calculated in two ways: (a) using the predicted valued by Sigma-RF and (b) using a constant value. For the constant sigma value, we used the average of the predicted sigma values in this study. For both proteins, using the predicted sigma value yielded significantly better correlations of the B-factors. For GB1, the CC value obtained with the predicted values was 0.49. While the constant sigma value led to a much worse CC value of 0.26. For LB1, a similar result was obtained; the predicted sigma values yielded a CC of 0.42 while a constant sigma value resulted in a CC of 0.14. We also observed similar results in our previous study with ADK.37 The arithmetic sum, S, of prefactors in eq 2 given by S=

∑ restraints ∈ K

1 σijK

(3)

This parameter S helps us identify which part of the protein contributes to the energy term. The values of wvdw and wrestraint in eq 1 were tuned by a careful grid search by short MD simulations performed with different combinations of wvdw and wrestraint. The final values of wvdw and wrestraint were determined to be 0.30 and 0.10, respectively, for GB1, and 0.25 and 0.09, respectively, for LB1. In each simulation, both proteins did not fold too fast at the beginning so that details of folding were observed, and at the 3434


Article

Journal of Chemical Theory and Computation averaging multiple simulation trajectories at 300 K, and the number of the trajectories used was 24 or 48. The kinetics were represented by the TM-score and also by the fractions of individual conformational states (the unfolded state, the native state, and also intermediates if present) as a function of time in nanoseconds (t). The kinetic folding curves thus obtained were fitted to an exponential function as given by i t yz zz zz k τi {

∑ ai expjjjjj− n

f (t ) = a 0 +

i=1

(4)

where f(t) is the TM-score or the fraction of each conformational state, a0 is the f(t) at infinite time, and ai and τi are the amplitude and the time constant, respectively, of the ith kinetic phase. The fitting parameters were a0, ai, and τi, and their best-fit values were evaluated by the nonlinear least-squares method with the Levenberg−Marquardt algorithm by using gnuplot 5.2 and Igor Pro, version 6.37 (Wave Metrics).

Figure 2. S-values map derived from the Lorentzian terms in eq 3.

Nine trajectories belonged to type B (Figure 3d−f): starting from a fully stretched structure, the protein folded through misfolded structures, stayed in an N-terminal misfolded form for a long time, and then reached the native structure. Six trajectories belonged to type C (Figure 3g−i): the beginning of the folding process for these was similar to that for type B, but the folding stalled and then passed through a C-terminal misfolded structure instead of the N-terminal misfolded structure before it reached the native state. Finally, nine trajectories belonged to type D (Figure 3j−l): all the trajectories got trapped in a misfolded structure, where transitions between the initial unfolded structure and the misfolded structure were very frequent. Some of type D trajectories then passed through a less nativelike structure before reaching the nativelike structure. In 45 out of the 48 simulations, the folded native state was reached at the end of 400 ns and hence the remaining 3 simulations of type B did not reach the native state. The final structures of these 3 simulations all had a misfolded N-terminus. For these 3 simulations, an additional 400 ns simulation was conducted starting from their final structures and MDsimulation conditions and all trajectories reached the native state. Using the 48 400 ns trajectories of GB1, the average TM-score was calculated as a function of folding time (t). The kinetic folding curve thus obtained may be equivalent to a folding reaction curve observed experimentally in bulk protein solution, and the curve shown could not be fitted by a single-exponential function (the RMSD of residuals of 0.015 and the variance of 0.0002) (Figure 4a). However, a double-exponential function fitted the folding curve of the averaged TM-score reasonably well (the RMSD of residuals of 0.0086 and the variance of 7.5 × 10 −5), apparently suggesting the presence of a stable kinetic folding intermediate between the initial unfolded and final native states. Here, τ1 and τ2 were estimated at 55.7 and 5.4 ns. The folding rates of the two phases thus differed by an order of magnitude. The detailed single- and double-exponential fitting results of the averaged TM-score of 48 trajectories of GB1 as well as the fitting residual plots are displayed in Figure S7. For LB1, all 48 simulations rapidly folded into the nativelike state with or without passing through a misfolded structure. The typical folding trajectories are shown in Figure 3m−o. The final TM-scores in all the 48 trajectories fluctuated in a slightly lower range than those of GB1, and the kinetic folding curve of the averaged TM-score could be fitted adequately well by a singleexponential function with a τ estimated at 11.4 ns (Figure 4b) (the RMSD of residuals of 0.0068 and the variance of 4.6 × 10 −5 ). The detailed single- and double-exponential fitting results of

■

RESULTS Native Structures of GB1 and LB1. Figure 1 represents the native structures of GB1 and LB1 by the residue−residue contact maps; the upper and lower triangles show the GB1 and LB1 structures, respectively. In Figure 1a, two residues are defined as contacting each other when the distance between their Cα atoms is less than 6 Å. If the two residues are in contact and also separated by four or more residues along the primary sequence, the contact formed by this pair of the residues was defined as a “Cα tertiary contact” or a long-range contact (Figure 1b). In Figure 1c, two residues are defined as contacting each other when two heavy atoms each from the different residues are within 4 Å, and the number of the atom−atom contacts between the two residues is indicated in the contact map by a color given by the color palette beside Figure 1d. If the two residues are in contact and the residue separation along the primary sequence is four or more, the contact is defined as a “tertiary contact” (Figure 1d). As shown by the Cα atom contacts (Figure 1a), GB1 and LB1 are quite similar in their backbone topologies, although their sequence identity is only 15%. However, when considering all heavy atom contacts (Figure 1c), there are significantly more atom−atom contacts in the C-terminal hairpin than in the Nterminal hairpin in GB1, while there are only slightly more contacts in the N-terminal hairpin than in the C-terminal hairpin in LB1. The same trend can be found in the heavy-atom longrange contacts as well (Figure 1d). The Lorentzian energy terms of eq 2 were generated from the reference native structure of each protein. Figure 2 shows the S values of the Lorentzian terms (eq 3). We observe higher values of this term in the C-terminal hairpin than in the N-terminal hairpin for GB1 (Figure 2). In LB1, however, the two hairpins have similar values. Multiple MD Simulations and the Folding Kinetics Monitored by the TM-Score. A total of 48 independent MD simulations, each for 400 ns at 300 K, were carried out for both GB1 and LB1 with different initial velocities. Here, we first describe the results of GB1. For GB1, the 48 trajectories followed quite different paths, and they were classified into four types. The representative trajectories are shown in Figure 3. A total of 24 trajectories belonged to type A (Figure 3a−c): folding to the native state occurred, and the TM-score jumped above 0.5 at an early stage. Some of these trajectories passed through a less nativelike structure before they reached the nativelike folded structure. 3435


Article


Figure 3. 400 ns folding trajectories for GB1 and LB1 in terms of the TM-score at 300 K: (a−l) typical ones for GB1 and (m−o) typical ones for LB1. All 48 trajectories of GB1 can be grouped into four types: (a−c) type A, in which the unfold structure folds to the nativelike structure immediately; (d− f) type B, in which the unfold structure folds through the N-terminal misfolded structure; (g−i) type C, in which the unfold structure folds through the C-terminal misfolded structure; (j−l) type D, in which the unfolded structure folds through the misfolded structure and the nativelike structure. For LB1 (m−o), folding in all trajectories happens at a very early stage and either passes the N-terminal misfolded structure or the nativelike structure.

the averaged TM-score of 48 trajectories of LB1 as well as the fitting residual plots are displayed in Figure S8. Interestingly, when averaging over only the type A trajectories of GB1, the folding kinetics could be fitted quite well by the single exponential with the fitting-parameter values that were close to those for LB1 (Figure 4c). When averaging over both

the type A and the type B trajectories, the single exponential could still fit the folding kinetics but not as good as the case of using only the type A trajectories (Figure 4d). When averaging over both the type A and the type C trajectories, the single exponential was not sufficient to characterize the kinetics, where a1 increases from almost 0 (Figure 4c) to 0.18 (Figure 4e) when 3436


Article


Figure 4. Averaged folding trajectories (red color) at 300 K in terms of the TM-score vs time (ns) with single- (green color) and double- (blue color) exponential fitting for (a) all 48 trajectories of GB1, (b) all 48 trajectories of LB1, (c) 24 type A trajectories of GB1, (d) 24 type A and 9 type B trajectories of GB1, (e) 24 type A and 6 type C trajectories of GB1, and (f) 24 type A and 9 type D trajectories of GB1.

the kinetics were fitted by the double-exponential function. The type D trajectories added even further the double-exponential feature to the folding kinetics (Figure 4f), where a2 decreases to 0.28 and a1 increases to 0.29. Clearly, the type C and D trajectories were mainly responsible for the folding kinetics of the multiexponential feature. Typical MD-simulation trajectories of GB1 and LB1 at the high temperature (414 K for GB1 and 393 K for LB1) are shown in Figure 5a−c and d−f, respectively. For both the proteins, all

24 trajectories oscillated frequently between the two TM-score ranges of 0.2−0.4 and 0.6−0.8. Reaction Path and the Density of States. Figure 6 shows potentials of mean force using the TM-score as reaction coordinate for GB1 at 300 and 414 K and for LB1 at 300 and 393 K. At the high temperature, two basins with almost equal populations are observed, and they may correspond to the unfolded (U) and the native (N) states (Figure 6c,d). The 3437


Article


Figure 5. Typical folding trajectories for GB1 (a−c) at 414 K and LB1 (d−f) at 393 K in terms of the TM-score.

There are only two basins, one basin is located in the RMSD range of 7−15 Å and the TM-score range is 0.15−0.35 and the other in the RMSD range of 2−5 Å and the TM-score range is 0.50−0.85 for both GB1 and LB1 (Figure 7b,d). Transition-State Analysis at High Temperature. Because of the advantage of using the Lorentzian energy function, we could observe numerous fast folding and unfolding events during the simulations at high temperature, 414 or 393 K, for GB1 and LB1, respectively. As shown in the free-energy profiles in Figure 6c,d, the TM-score alone is sufficient for picturing the folding process of GB1 and LB1 at the high temperature. As described in the Methods section, the structures located in the TM-scores range of 0.40−0.55 were identified as the transition-state structures, and 338 frames for GB1 and 186 frames for LB1 were randomly selected in this TM-score range for clustering. Using the in-house community/module detection method, Mod-CSA, 338 frames of GB1 were grouped into three clusters (Figure 8): TS1 (144 frames), TS2 (131 frames), and TS3 (63 frames). The center structure of TS1 (Figure 8a) had a TMscore of 0.529, and when we traced the 144 structures back to the original trajectories, 124 eventually folded to the N state. The center structure of TS2 (Figure 8b) had a TM-score of 0.474, and its C-terminus was misfolded. Among the 131 structures, 40 structures folded to the N state. The center structure of TS3 had a TM-score of 0.494 (Figure 8c), among the 63 structures, 37 structures folded to the N state. The central α-helix was largely folded in the three center structures of each community, and the C-terminal hairpins in both the TS1 and TS3 center structures had more atom−atom contacts than the N-terminal hairpin. Large contacts were also formed between strands 1 and 4. These features were not observed in the TS2 center structure. For LB1, 186 frames were grouped into 2 clusters: TS1 (99 frames) and TS2 (87 frames). The cluster center structure of TS1 (Figure 9a) had a TM-score of 0.492, in which the Nterminal hairpin was largely formed, and when we traced these 99 frames back into the original trajectories, only 19% folded to the N state. For TS2, however, 89% of the 87 structures folded to the N state, and the center structure of TS2 had a TM-score of

barrier between the two basins was estimated be be less than 2kBT for both GB1 and LB1. The potential of mean force at 300 K was quite different from that at the high temperature for both GB1 and LB1, and its profile at 300 K was more complicated. For GB1, two basins are located at TM-scores around 0.30 and 0.85 and the higher TMscore basin has higher population (Figure 6a,b). These two basins may correspond to the initial U state and the final N state. In addition to these two basins, a small third basin is detected around a TM-score of 0.65 and this may correspond to an intermediate (I) state. However, the barrier from U to I is higher than that from I to N, suggesting that the I state may not accumulate to an experimentally detectable extent, and hence, the folding kinetics should be fitted by a single exponential function if the folding can be fully described by the chosen reaction coordinate. This question will be addressed later in this manuscript (section Exponential Fitting Using Fraction of Different Conformers Formed on GB1 Folding Path). For LB1, only the batch of nativelike folded states around a TM-score of 0.8 became largely accumulated. Obviously, using only the TM-score as the reaction coordinate appears not sufficient for describing the folding process. We thus decided to use RMSD as well as the TM-score for our analysis. Figure 7 shows the free-energy profiles for GB1 and LB1 with the TM-score and the RMSD as the twodimensional (2D) reaction coordinates. Now the potentials of mean force for the two proteins are better revealed, and at 300 K, the basin of the U state, ranging from 9 to 12 Å of RMSD, covers a wider range of the TM-score in GB1 than it does in LB1. In the potential profile of LB1, the U basin did not show up clearly when the TM-score was used as the reaction coordinate (Figure 6b). The existence of this U-state basin became more evident when the RMSD instead of TM-score was used as the reaction coordinate (Figure 7c). The I basin was observed at a TM-score of 0.6−0.7 and at an RMSD of 3−4 Å for GB1, but this feature did not exist for LB1. The barrier height from U to I was about 5kBT for GB1, and this was higher than that from I to N. For LB1, the barrier from U to N at 300 K was much smaller, i.e., less than 2kBT. The free energy profile at the high temperature (414 K for GB1 and 393 K for LB1) was simpler for the two proteins. 3438


Article


Figure 6. Potentials of mean force using the TM-score as a reaction coordinate for GB1 and LB1: (a) using all 48 trajectories of GB1 at 300 K, (b) using all 48 trajectories of LB1 at 300 K, (c) all 24 trajectories of GB1 at 414 K, and (d) all 24 trajectories of LB1 at 393 K.

0.534 (Figure 9b), in which both of the hairpins were largely formed. Detailed Analysis of Populating States during Folding at 300 K. Because the folding at 300 K was much more complicated than that at the high temperature, especially for GB1, a more careful and detailed analysis was required for fully understanding the folding events. We again applied the ModCSA clustering algorithm for the analysis but on a group of structures which covered a much larger TM-score range (the TM-score of 0.2−0.9 and the RMSD of 0.5−14.5 Å). We found 9265 structures in this range from the folding phases of 300 K trajectories. We then randomly selected 10% of these frames, calculated all pairwise RMSDs, and used a cutoff value of 4.0 Å to construct a weighted network of the size of 662 frames. ModCSA identified 24 clusters, and the center structures were determined. Next, we grouped the whole 9265 structures into these 24 clusters based on the similarity of each frame to the 24 clustercenter structures. The total number of clusters was thus still 24, but the center of each cluster shifted slightly after including more structures. Among the 24 clusters, we found 7 biggest clusters with the size of 1860 to 307 frames; the rest of the 17 clusters

covered less than 4% of the 9265 structures and a TM-score range of 0.26−0.52. The center structures of these 7 clusters are shown in Figure 10, where all the structures are aligned in the way that the helix goes from left (N-terminus) to right (Cterminus); in the following, the seven clusters are denoted by Roman numerals, I, II, III, IV, V, VI, and VII. Most of the frames in the first biggest cluster (cluster I) had the closely native structure, with the center frame having a TMscore of 0.85 and an RMSD of 1.78 Å (Figure 10a). Cluster I may thus be regarded as the N state. The second biggest cluster (cluster II) included the structures which had the misfolded Nterminus, and the cluster center structure had a TM-score of 0.64 and an RMSD of 3.78 Å (Figure 10b). The third biggest cluster (cluster III) was composed of the misfolded C-terminus with the center structures having a TM-score of 0.70 and an RMSD of 3.47 Å (Figure 10c). Most frames in the fourth biggest cluster (cluster IV) were unfolded (Figure 10d) with the center structure having a much lower TM-score (0.26) and a much higher RMSD (9.89 Å), and this cluster was considered as the initial unfolded state. The center structures of the following two clusters (clusters V and VI) were both nativelike, and their TMscores were 0.71 (cluster V) (Figure 10e) and 0.59 (cluster VI) 3439


Article


Figure 7. 2D free-energy profiles using the TM-score and RMSD as two-dimensional reaction coordinates for GB1 at (a) 300 K and (b) 414 K and for LB1 at (c) 300 K and (d) 393 K.

misfolded state (cluster II) and the C-terminal misfolded state (cluster III), and the protein stayed in these structures for a long period of time and then finally folded into the N state (cluster I). We also performed the Mod-CSA analysis for the LB1 trajectories at 300 K using a total of 925 frames, which covered a TM-score range of 0.16−0.88 and an RMSD range of 41.78− 1.78 Å with a cutoff value of 5.0 Å for pairwise RMSDs. We identified eight clusters, and whole 5060 structures were grouped into the 8 clusters. Among these eight clusters, the structures in the four biggest clusters (clusters I, II, III, and IV) having higher TM-scores (TM-scores above 0.74) than the other clusters were very close to the native structure (Figure 11). The fifth cluster (cluster V) was much smaller in size; its center structure had a TM-score of 0.54, and the protein had the misfolded N-terminus. Structures in the rest of the clusters (clusters VI, VII, and VIII) were all misfolded or unfolded with much smaller populations. The atom−atom contact maps of the center structures of these eight clusters are illustrated in Figure S10. To investigate the roles of the different conformational states, represented by the different clusters, in the folding reactions, we

(Figure 10f). The last cluster (cluster VII) had the center structure similar to that of cluster IV, but the C-terminal hairpin already folded (Figure 10g); however, the orientation of the third and fourth strands was different from that in cluster IV. Cluster VII may be regarded as an initial misfolded state. The atom−atom contact maps of the center structures of these 7 clusters are shown in Figure S9, where the C-terminal hairpin is slightly better formed than the N-terminal hairpin in parts b, c, and f (clusters II, III, and VI), while the N-terminal hairpin is better formed in part e (cluster V). Because clusters II, III, and VI may represent the kinetic folding intermediates, this result may be in accordance with the previous result on the transitionstate structure of GB1 at 414 K, where we have observed the Cterminal hairpin more structured than the N-terminal one. When we mapped these structures of the seven clusters back to the original trajectories, we found different folding paths connecting these structures: (i) there were fast exchanges between clusters IV, VII (the initial unfolded state (part d) and the initial misfolded state (part g)) and the nativelike I states (cluster V and VI) then the protein folded into the N state (cluster I); (ii) cluster IV can also form the N-terminal 3440


Article


Figure 8. Predicted 3 communities of the 338 frames of transition states for GB1: (a) TS1 is the largest community, and the RMSD center structure has the highest TM-score of 0.53; (b) TS2 is the second largest community, and the RMSD center structure has the lowest TM-score of 0.47 and a misfolded C-terminus; (c) TS3 is the smallest community, and the RMSD center structure has the a TM-score of 0.49. For each TS, the atom−atom contacts in two residues less than 4 Å are shown in the upper triangles, while those in the native structure are shown in the lower triangles as comparison.

calculated the 2D potentials of mean force using the RMSDs from two different structures as the X−Y reaction coordinates. We employed the RMSDs from the native structure (1PGB for GB1 and 1HZ6 for LB1) and from the N-terminal misfolded structure (Figure 10b (cluster II) for GB1 and Figure 11e (cluster V) for LB1). The 2D free-energy profiles thus generated from the whole 48 trajectories of GB1 and LB1 are shown in Figure 12a,b, respectively. The biggest difference of the two proteins is observed around the following RMSD values: 3 Å from the N-terminal misfolded structure and 3.5−4 Å from the native structure. We observe the third basin with remarkable deepness in this location of the free-energy profile of GB1 (Figure 12a), but this is not observed in the LB1 free-energy profile (Figure 12b). The states populated in this location of GB1 were mainly composed of the N-terminal misfolded and the

C-terminal misfolded species (clusters II and III), and there is a significant free-energy barrier between the third basin and the final N state (Figure 12a). The presence of the deep third basin in GB1 and the presence of the free-energy barrier between the third basin and the N state reasonably explain the doubleexponential folding kinetics of GB1, as the misfolded intermediate, accumulated in the third basin at an early stage of folding, may slowly fold into the N state. For LB1, the third basin is much shallower (Figure 12b), so that the accumulation of the misfolded species is insignificant, leading to the apparent two-state single-exponential kinetics of folding. It is important to note that when we used just the TM-score as a reaction coordinate (Figures 6 and 7), we could not clearly distinguish between the misfolded species and the N state, and hence, the double-exponential character of the GB1 folding could not 3441


Article


Figure 9. Predicted 2 communities of the 186 frames of transition states for LB1: (a) TS1 is the bigger community, and the RMSD center structure has a TM-score 0.49; (b) TS2 is the smaller community, and the RMSD center structure has a higher TM-score of 0.53. For each TS, the atom−atom contacts in two residues less than 4 Å are shown in the upper triangles, while those in the native structure are shown in the lower triangles as comparison.

When all the conformations which were similar to Figure 10b−g were assumed as the unfolded state, both the formation of the native state (Figure 10a) and the disappearing of the unfolded state of GB1 could be well fitted by two exponentials with reasonably small residuals (Figure 13). Here, τ1 and τ2 were determined to be 3.9 and 59.0 ns, and hence, they were 1 order of magnitude different. The same analysis was also performed for LB1 using the center structures of the 8 clusters (Figure 11). All the conformations similar to the structures in Figure 11e−h were considered as unfolded, and all the conformations similar to the structures in Figure 11a−d were considered as native. Both the formation of the native state and the disappearance of the unfolded state could be well fitted by a single exponential function with small residuals (Figure 14). Here, τ1 was predicted to be 8.6 ns. Further analysis was performed for GB1 by separating the intermediate conformers, belonging to clusters II, III, V, and VI (Figure 10b,c,e,f) from the unfolded conformers, belonging to clusters IV and VII (Figure 10d,g). Global fitting was then carried out for the time courses of the fractions (f U, f I, and f N) of the three states (U, I, and N) (Figure 15). The best results to fit all the three time courses were obtained by three exponentials. The two major phases were found to have τ values (τ1 = 5.8 ns and τ2 = 69.2 ns) similar to those in the earlier doubleexponential fitting. However, we found a lag phase at an early stage before 1 ns for the formation of the N state. The lag phase is represented by the fastest phase of the three exponentials for the N state (a3 = 0.02 and τ3 = 0.52 ns), and this is a strong indication that the intermediate (I) is an on-pathway folding

reasonably be interpreted. It is also noticeable that when we used only the 24 type A trajectories of GB1, the free-energy profile became very similar to that of LB1 (Figure 12c), leading to the apparent single-exponential kinetics of folding (Figure 4c). Another difference between GB1 and LB1 was also found in the U state in a high RMSD region with both RMSDX and RMSDY values of around 11 Å. We observed more diverse conformations in this region on the GB1 potential energy surface than on the LB1. The conformations distributed in this region were composed of the initial U state (cluster IV) and the initial misfolded state (cluster VII) in GB1, whereas only the U state was present in this region in LB1. The initial U state of GB1 can thus go either directly to the nativelike I states (clusters V and VI) or go through the misfolded I states (clusters II and III) via the initial misfolded state (cluster VII) then to the N state. For LB1, however, the majority of the U state go directly to the N state. Folding Kinetics of GB1 Monitored by the Fractions of Different Conformational States. The folding kinetics of GB1 and LB1 were also monitored by the fractions of different conformational states. To this end, we used the whole 48 trajectories at 300 K for each protein and calculated the RMSDs for each single frame with respect to the center structures of the 7 (GB1) or 8 (LB1) different clusters (Figures 10 and 11). Each frame belonged to one of the 7 (or 8) conformations if the RMSD from its center structure was the smallest. The fractions of the different conformational states thus obtained at different folding times were then calculated for further analysis. 3442


Article


Figure 11. Center structures of the predicted 8 total communities of the 5060 frames of LB1: (a) RMSD center structure of the largest community (cluster I), which has the lowest RMSD of 1.81 Å and highest TM-score of 0.87; (b) RMSD center structure of the second largest community (cluster II), which has an RMSD of 2.43 Å and a TM-score of 0.83; (c) RMSD center structure of community 3 (cluster III), which has a TM-score of 0.74 and an RMSD of 2.89 Å; (d) RMSD center structure of the fourth largest group (cluster IV), which has a TM-score of 0.84 and an RMSD of 2.41 Å; (e) RMSD center structure for community 5 (cluster V), which has a TM-score of 0.54 and an RMSD of 4.90 Å, and all the structures have a misfolded N-terminus. The RMSD center structures for communities 6 and 7 (clusters VI, VII, and VIII) are all misfolded with a TM-score ≤0.54.

Figure 10. Center structures of the predicted biggest 7 communities of the 9265 frames of GB1: (a) RMSD center structure of the largest community (cluster I), which has the lowest RMSD of 1.78 Å and highest TM-score of 0.85. (b) RMSD center structure of the second largest community (cluster II), which has an RMSD of 3.78 Å and a TM-score of 0.64; all structures have a misfolded N-terminus. (c) RMSD center structure of community 3 (cluster III), which has a TMscore of 0.70 and an RMSD of 3.47 Å; all structures contain the misfolded C-terminus. (d) RMSD center structure of the fourth largest group (cluster IV), which has the lowest TM-score of 0.26 and the highest RMSD of 9.89 Å; all structures are misfolded. (e and f) RMSD center structures for communities 5 and 6 (clusters V and VI), and the TM-scores are 0.71 and 0.59, respectively; all the structures in both groups are folded correctly, but cluster V is more nativelike. (g) RMSD center structure in community 7 (cluster VII), which is similar to the fourth community center structure but has a better formed second hairpin; TM-score is 0.31.

belonging to U in the above analysis, there exist a number of intermediates similar to the structures in parts b, c, e, and f in Figure 10. Among the four, the structures in parts e and f can be grouped as the I1 state which are more nativelike, while the structures in parts b and c, both having a misfolded conformation, can be grouped as the I2 state. Before these four I structures, there also exist such low TM-score and high RMSD I structures (grouped as Ib), which have relatively low population and show a rapid pre-equilibrium with the U state. The folding reaction is thus represented by the transition from the U plus I b to the N state via I1 or I2. The lag phase was observed in the fitting of the GB1 folding curve (Figure 15), and this suggests the direct folding pathway from the U to N state may not exist.

intermediate.47 In fact, when we focused on the region around 0.5 ns, we observed not only a decreasing phase of U but also an increasing phase of I, indicating that I was accumulating while consuming U and that the formation of N was delayed during accumulation of I. It should also be noted that the f I value extrapolated back to zero time was significantly positive (f I = 0.26 at t = 0 ns). This is probably due to very fast exchanges occurring between some of low TM-score and high RMSD conformers (Ib) that may have structures between I and U. Clearly, more than one intermediate conformational state contributed to the folding reaction of GB1. The proposed mechanism of the GB1 folding is shown in Figure 16. Besides the initial unfolded and misfolded states, both

■

DISCUSSION GB1 and LB1 have been used extensively as a model system in computational as well as experimental studies of protein folding 3443


Article


Figure 12. 2D free-energy profiles using RMSDs of each frame to the native structure and the N-terminal misfolded structure for GB1 and LB1 at 300 K. (a) Analysis on all 48 trajectories of GB1, U and high RMSD I (cluster IV and VII) structures occupy the high RMSD region, the nativelike I (clusters V and VI) and misfolded N-terminal I (cluster II) as well as the misfolded C-terminal I (cluster III) at RMSDs (3,4) Å, N structures (cluster I) at RMSDs (4,2) Å; (b) analysis on all 48 trajectories of LB1, there are mostly the two wells composed of U and N structures; (c) analysis on 24 type A trajectories of GB1; (d) analysis on 24 type A and 9 type B trajectories of GB1; (e) analysis on 24 type A and 6 type C trajectories of GB1; and (f) analysis on 24 type A and 9 type D trajectories of GB1.

Different from the conventional computational approaches, the Lorentzian restraint energy terms were included in addition to the stereochemistry and van der Waals energy terms. By using this method, we successfully observed numerous folding− unfolding events at high temperature (414 K for GB1 and 393 K for LB1) and many folding processes at 300 K for both the proteins.

over the years. Because of the complexity of the molecular mechanisms of protein folding themselves, various computational methods have so far been developed, including all-atom unfolding simulations at a high temperature,6,14,15,48−51 the use of a biased potential12−14,52−55 (such as a Go̅-like model56−61), or the use of simplified models of a protein.17,19,20,62−65 Here, we used the Lorentzian structure-based potential in the all-atom MD simulations to study the folding kinetics of GB1 and LB1. 3444


Article


Figure 14. Single-exponential fitting for LB1 native and unfold conformers accumulated in the 48 trajectories at 300 K. The native conformers all have similar structures to cluster center structures in Figure 11a−d, the unfolded conformers all have similar structures to cluster center structures in Figure 11e−h. The fitting function is f(x) = 1.00 − 0.92 exp(−x/8.63), g(x) = 0.92 exp(−x/8.63).

Figure 13. Double-exponential fitting for GB1 native and unfolded conformers accumulated in the 48 trajectories at 300 K. The native conformers all have similar structures to cluster center structures in Figure 10a, and the unfolded conformers all have similar structures to the cluster center structures in Figure 10b−g. The fitting function is f(x) = 0.96 − 0.57 exp(−x/58.95) − 0.39 exp(−x/3.90), g(x) = 0.04 + 0.57 exp(−x/58.95) + 0.39 exp(−x/3.90).

terminal hairpin (Figure 8). This is in good agreement with experimental results by McCallister and co-workers.25 In their experiment, mutational ϕ-value analysis observed higher ϕ values in the second β-hairpin and the second β-turn than that in the first β-hairpin and helix, which indicated that the second hairpin was largely formed in the folding transition state, but the first hairpin and helix was relatively disordered. While for LB1, mutational ϕ-value analysis observed only higher ϕ values in the first β-hairpin in the transition state.27 However, ψ-value analysis with engineered metal ion binding sites found that β1−β4, β3−β4, and β1−β2 all have a ψ value higher than 0.75 and a much lower value in the helix, which indicated that the entire four-stranded β-sheet was formed.54 In our study, Mod-CSA identified two TSs for LB1 at 393 K, with one that has a better formed first β-hairpin and the other has both hairpins equally formed (Figure 9). Our findings again are in good agreement with the experimental results. Folding Mechanisms at 300 K as Revealed by the Simulations. For GB1, we identified the multiple intermediate

The kinetic times for folding obtained in this study may not be directly connected with the real physical time of folding. However, the significant differences of folding behavior were observed between the two proteins, and this gives a strong indication that the two proteins fold apparently in a different manner. The folding simulations using the Go̅-like potential, also known as a structure-based potential, have been widely used in previous studies for the purpose to obtain the information on the kinetics of protein folding.59,60 Our model keeps the all-atom details that may be important for describing the dynamics of biomolecules.61 Structures of the Transition-State Ensemble. Using Mod-CSA, three transition states have been identified for GB1 at 414 K. Among the three, we found that the C-terminal hairpin of both TS1 and TS3 have more atom−atom contacts than the N3445


Article


terminal and the C-terminal misfolded I states. While a large portion of the molecule may already be quite nativelike with a similar high TM-score, they either lack correct side-chain interactions or possess a misfolded conformation in a certain region. This means that it will take a long time for the protein to fix such “wrong” interactions. Because of complexity of the GB1 folding, neither the 1D (TM-score) potential (Figure 6) nor the conventional 2D (TM-score/RMSD) potential (Figure 6) of the mean force could explain the non-two-state folding behavior of the protein. Because the I states had these parameter values close to those of the N state, we could not easily distinguish I from N. Using the RMSD to the N-terminal misfolded I state, we could clearly distinguish between I and N and better describe the folding processes of GB1 at 300 K. This has an important implication that the conventional 2D potential of mean force, often employed in the MD trajectory analysis, may not be useful for the folding reactions in which the misfolded intermediate is significantly accumulated. As shown in the folding model of Figure 16, there are multiple parallel pathways of the GB1 folding at 300 K, and on each of the folding pathways, a rapidly formed I state is accumulated. The presence of such multiple parallel pathways has also been supported by recent computational as well as experimental studies of the folding kinetics of GB1.3,6,14 The global fitting analysis of the time courses of the fractional U, I, and N states indicated that about 50% of the molecules exhibited the fast phase of folding with a time constant of 5.8 ns after the formation of the I state (Figure 15). This fast phase of folding may correspond to the Type A trajectories depicted in Figure 4c, because the average of the Type A trajectories resulted in the single-exponential kinetics with a time constant very close to that of the fast phase. The I state accumulated in the fast phase may correspond to the nativelike I state, because the nativelike structures were observed in the Type A trajectories. The global fitting analysis also indicated that the remaining 50% of the molecules folded in the slow phase with a time constant of 69.2 ns (Figure 15). The I states accumulated in this slow phase may correspond to various misfolded I states, including the Nterminal and the C-terminal misfolded I states, and hence the misfolding events significantly decelerated the folding kinetics of GB1. For LB1, the folding kinetics were much simpler and exhibited the single-exponential behavior at 300 K (Figures 4b and 14). However, our detailed clustering analysis has also revealed the presence of the N-terminal misfolded and the nativelike intermediates, similar to the I states of GB1, during the LB1 folding (Figures 11 and 12b). Therefore, the basic molecular mechanisms of folding are similar between the two homologous proteins, GB1 and LB1, and this may be not inconsistent with a recent report by Baxa et al.,55 who have proposed that homologous proteins have similar folding mechanisms even when non-native interactions are present in the transition state. However, the accumulation of the misfolded intermediate was much smaller in LB1 than in GB1, so that the LB1 folding was well approximated as the three-state, U ⇌ I ⇌ N, mechanism with the nativelike intermediate or the even simpler two-state, U ⇌ N, mechanism. The N-terminal and C-terminal hairpins are similar in length and similar in the number of atom−atom contacts, and this is in contrast with the hairpin structures in native GB1 (Figure 1). Furthermore, the interactions of the Nterminal and C-terminal hairpins with the central α helix distribute more widely in LB1 than in GB1 (Figure 1c). These differences in the hairpin structures between the two proteins

Figure 15. Three exponential fitting for GB1 native, intermediate, and unfolded conformers accumulated in the 48 trajectories at 300 K. The native conformers all have similar structures to cluster center structures in Figure 10a; the intermediate conformers have similar structures to cluster center structures in Figure 10b,c,e,f; the unfolded conformers all have similar structures to the cluster center structures in Figure 10d,g. The fitting function is f(x) = 0.97 − 0.50 exp(−x/69.19) − 0.46 exp(−x/5.83) + 0.02 exp(− x/0.52), g(x) = 0.00 + 0.13 exp(−x/69.19) + 0.39 exp(−x/5.83) + 0.19 exp(−x/0.52), h(x) = 0.03 + 0.37 exp(−x/ 69.19) + 0.07 exp(−x/5.83) − 0.21 exp(−x/0.52).

states in the simulation trajectories at 300 K, including the Nterminal misfolded I state (cluster II), the C-terminal misfolded I state (cluster III), and the nativelike I state (clusters V and VI). Although both the TM-score 1D potential and the TM-score/ RMSD 2D potential of the mean force indicated that GB1 folded in a more complicated manner than LB1 (Figures 6 and 7), these potentials could not well describe the accumulation of the above I states during the folding of GB1. The complexity of the GB1 folding was caused by misfolding events. There are more and stronger contacts in the C-terminal hairpin than in the Nterminal hairpin in the native structure of GB1 (Figure 1), and this might lead to faster folding of the C-terminal region without correct folding of the rest of the protein molecule. Therefore, this may increase the possibility that the protein forms the N3446


Article


Figure 16. Possible folding mechanism of GB1. There is a rapid pre-equilibrium between between U (d and g) and Ib within the burst phase, and the exchange between Ib and I1 (e and f) is fast and exchange between Ib and I2 (b and c) is slower, but both I1 and I2 leaded to N (a).

required for analyzing their potential energy surfaces. For GB1, we observed the U state, the multiple I states (N-terminal misfolded, C-terminal misfolded, and nativelike) and the N state during the kinetic folding process, and the I states were accumulated to a significant extent during the process. For LB1, we also identified the U, the I (nativelike and N-terminal misfolded) and the N states, but the I states were not significantly accumulated, leading to the apparent two-state behavior of LB1. Therefore, the folding kinetics of GB1 were better fitted by a double-exponential function, while the LB1 folding kinetics were fitted well by a single-exponential function.

might give us a reasonable interpretation for the difference in their folding behaviors, and the folding of LB1 may occur more cooperatively without significant accumulation of the Nterminal misfolded intermediate. Comparison with Previous Simulation Studies on the Kinetics of GB1 and LB1. Brown and Head-Gordon studied the folding kinetics of GB1 and LB1 by constant-temperature Langevin dynamics.19 They used a minimalist protein model represented as a sequence of hydrophilic, hydrophobic, and neutral beads. Using such a very simplified model, they also observed double-exponential kinetics for GB1 while the LB1 folding was single-exponential. However, the difference in the time constant between the fast and slow phases of the GB1 folding was only 3.4-fold, much smaller than that observed in the present study. In the present study, we observed 1 order of magnitude difference between the two phases of the GB1 folding (Figures 4a and 13), being consistent with the experimentally observed difference between the two phases.24 Therefore, the all-atom presentation of the protein model may be required for correct description of the folding kinetics. Shimada and Shachnovich studied the GB1 folding by an all-atom Monte Carlo simulation with the use of a square-well Go̅ potential for native atom−atom contacts.14 They also observed the presence of the multiple parallel pathways in the GB1 folding, but the folding kinetics after the burst phase were single-exponential, in contrast to the double-exponential kinetics observed in the present study as well as in the real folding experiments,24 suggesting that the use of the Lorentzian structure-based potential may be superior to the use of the square-well potential.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.8b01163. Figures S1−S4, MD simulations for van der Waals and restraint weights selection for GB1 and LB1 at 300 K; Figures S5 and S6, replica-exchange simulation for GB1 and LB1 for high-temperature selection; Figures S7 and S8, fitting details of the averaged folding curve of all 48 trajectories of GB1 and LB1; Figures S9 and S10, atom− atom contact maps of the 7 GB1 and 8 LB1 center structures shown in Figure 10 and 11 (ZIP)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

■

Qianyi Cheng: 0000-0002-4640-2238 Jooyoung Lee: 0000-0002-4432-6163

CONCLUSIONS In this work, all-atom MD simulation with structure-based Lorentzian restraint energy terms was used to study the GB1 and LB1 folding. Frequent folding and unfolding events were observed in the simulations of GB1 and LB1 at 414 and 393 K, respectively. At the high temperature, three and two transitionstate structures for GB1 and LB1 were predicted, respectively. The structural features of these transition states were in good agreement with the experimental transition-state analysis. At 300 K, the folding processes were more complicated for both the proteins, and the 2D reaction coordinates of the RMSDs from the N state and from the N-terminal misfolded I state were

Funding

This research was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science and ICT (NRF2017R1E1A1A01077717) and JSPS KAKENHI Grant Number JP16K07314. Authors also acknowledge Center for Advanced Computation (CAC) in the Korea Institute for Advanced Study (KIAS) for providing computing resources (Linux Cluster System). Notes

The authors declare no competing financial interest. 3447


Article


■ ■

(23) Park, S.-H.; O’Neil, K. T.; Roder, H. An early intermediate in the folding reaction of the B1 domain of protein G contains a native-like core. Biochemistry 1997, 36, 14277−14283. (24) Park, S.-H.; Shastry, M.; Roder, H. Folding dynamics of the B1 domain of protein G explored by ultrarapid mixing. Nat. Struct. Biol. 1999, 6, 943−947. (25) McCallister, E. L.; Alm, E.; Baker, D. Critical role of β-hairpin formation in protein G folding. Nat. Struct. Biol. 2000, 7, 669−673. (26) Gu, H.; Kim, D.; Baker, D. Contrasting Roles for Symmetrically Disposed b-Turns in the Folding of a Small Protein. J. Mol. Biol. 1997, 274, 588−596. (27) Kim, D. E.; Fisher, C.; Baker, D. A breakdown of symmetry in the folding transition state of protein L. J. Mol. Biol. 2000, 298, 971−984. (28) Dror, R. O.; Dirks, R. M.; Grossman, J.; Xu, H.; Shaw, D. E. Biomolecular Simulation: A Computational Microscope for Molecular Biology. Annu. Rev. Biophys. 2012, 41, 429−452. (29) Perez, A.; Morrone, J. A.; Simmerling, C.; Dill, K. A. Advances in free-energy-based simulations of protein folding and ligand binding. Curr. Opin. Struct. Biol. 2016, 36, 25−31. (30) Mayor, U.; Johnson, C. M.; Daggett, V.; Fersht, A. R. Protein folding and unfolding in microseconds to nanoseconds by experiment and simulation. Proc. Natl. Acad. Sci. U. S. A. 2000, 97, 13518−13522. (31) Alonso, D. O. V.; Daggett, V. Staphylococcal protein A: Unfolding pathways, unfolded states, and differences between the B and E domains. Proc. Natl. Acad. Sci. U. S. A. 2000, 97, 133−138. (32) Pan, Y.; Daggett, V. Direct Comparison of Experimental and Calculated Folding Free Energies for Hydrophobic Deletion Mutants of Chymotrypsin Inhibitor 2: Free Energy Perturbation Calculations Using Transition and Denatured States from Molecular Dynamics Simulations of Unfoldi. Biochemistry 2001, 40, 2723−2731. (33) Sheinerman, F. B.; Brooks, C. L. Calculations on folding of segment B1 of streptococcal protein G 1 1Edited by F. Cohen. J. Mol. Biol. 1998, 278, 439−456. (34) Shea, J.-E.; Brooks, C. L., III From folding theories to folding proteins: a review and assessment of simulation studies of protein folding and unfolding. Annu. Rev. Phys. Chem. 2001, 52, 499−535. (35) Nguyen, H.; Maier, J.; Huang, H.; Perrone, V.; Simmerling, C. Folding simulations for proteins with diverse topologies are accessible in days with a physics-based force field and implicit solvent. J. Am. Chem. Soc. 2014, 136, 13959−13962. (36) Lindorff-Larsen, K.; Piana, S.; Dror, R. O.; Shaw, D. E. How fastfolding proteins fold. Science 2011, 334, 517−520. (37) Lee, J.; Joo, K.; Brooks, B. R.; Lee, J. The Atomistic Mechanism of Conformational Transition of Adenylate Kinase Investigated by Lorentzian Structure-Based Potential. J. Chem. Theory Comput. 2015, 11, 3211−3224. (38) Kundrot, C. E.; Ponder, J. W.; Richards, F. M. Algorithms for calculating excluded volume and its derivatives as a function of molecular conformation and their use in energy minimization. J. Comput. Chem. 1991, 12, 402−409. (39) Mackerell, A. D., Jr.; Feig, M.; Brooks, C. L., III Extending the treatment of backbone energetics in protein force fields: Limitations of gas-phase quantum mechanics in reproducing protein conformational distributions in molecular dynamics simulations. J. Comput. Chem. 2004, 25, 1400−1415. (40) Fiser, A.; Sali, A. Modeller: Generation and Refinement of Homology-Based Protein Structure Models. Methods Enzymol. 2003, 374, 461−491. (41) Joo, K.; Joung, I.; Lee, J.; Lee, J.; Lee, W.; Brooks, B.; et al. Protein structure determination by conformational space annealing using NMR geometric restraints. Proteins: Struct., Funct., Genet. 2015, 83, 2251− 2262. (42) Lee, J.; Lee, K.; Joung, I.; Joo, K.; Brooks, B. R.; Lee, J. Sigma-RF: prediction of the variability of spatial restraints in template-based modeling by random forest. BMC Bioinf. 2015, 16, 94. (43) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 1984, 81, 3684−8.

ACKNOWLEDGMENTS We thank Professor Ayori Mitsutake for insightful discussion. REFERENCES

(1) Maity, H.; Reddy, G. Folding of Protein L with Implications for Collapse in the Denatured State Ensemble. J. Am. Chem. Soc. 2016, 138, 2609−2616. (2) Skinner, J. J.; Yu, W.; Gichana, E. K.; Baxa, M. C.; Hinshaw, J. R.; Freed, K. F.; et al. Benchmarking all-atom simulations using hydrogen exchange. Proc. Natl. Acad. Sci. U. S. A. 2014, 111, 15975−15980. (3) Lapidus, L. J.; Acharya, S.; Schwantes, C. R.; Wu, L.; Shukla, D.; King, M.; et al. Complex Pathways in Folding of Protein G Explored by Simulation and Experiment. Biophys. J. 2014, 107, 947−955. (4) Granata, D.; Camilloni, C.; Vendruscolo, M.; Laio, A. Characterization of the free-energy landscapes of proteins by NMR-guided metadynamics. Proc. Natl. Acad. Sci. U. S. A. 2013, 110, 6817−6822. (5) Shenker, S.; O’Donnell, C. W.; Devadas, S.; Berger, B.; Waldispuhl, J. Efficient traversal of beta-sheet protein folding pathways using ensemble models. J. Comput. Biol. 2011, 18, 1635−1647. (6) Camilloni, C.; Broglia, R. A.; Tiana, G. Hierarchy of folding and unfolding events of protein G, CI2, and ACBP from explicit-solvent simulations. J. Chem. Phys. 2011, 134, No. 045105. (7) Guarnera, E.; Pellarin, R.; Caflisch, A. How does a simplifiedsequence protein fold? Biophys. J. 2009, 97, 1737−1746. (8) Narzi, D.; Daidone, I.; Amadei, A.; Di Nola, A. Protein Folding Pathways Revealed by Essential Dynamics Sampling. J. Chem. Theory Comput. 2008, 4, 1940−1948. (9) Banushkina, P.; Meuwly, M. Diffusive dynamics on multidimensional rough free energy surfaces. J. Chem. Phys. 2007, 127, 135101. (10) Ekonomiuk, D.; Kielbasinski, M.; Kolinski, A. Protein modeling with reduced representation: statistical potentials and protein folding mechanism. Acta Biochim. Polym. 2005, 52, 741−748. (11) Hubner, I. A.; Shimada, J.; Shakhnovich, E. I. Commitment and Nucleation in the Protein G Transition State. J. Mol. Biol. 2004, 336, 745−761. (12) Yup Lee, S.; Fujitsuka, Y.; Kim, D. H.; Takada, S. Roles of physical interactions in determining protein-folding mechanisms: molecular simulation of protein G and α spectrin SH3. Proteins: Struct., Funct., Genet. 2004, 55, 128−138. (13) Song, G.; Thomas, S.; Dill, K. A.; Scholtz, J. M.; Amato, N. M. A path planning-based study of protein folding with a case study of hairpin formation in protein G and L. Pac. Symp. Biocomput. 2003, 240−251. (14) Shimada, J.; Shakhnovich, E. I. The ensemble folding kinetics of protein G from an all-atom Monte Carlo simulation. Proc. Natl. Acad. Sci. U. S. A. 2002, 99, 11175−11180. (15) Sheinerman, F. B.; Brooks, C. L., III Calculations on folding of segment B1 of streptococcal protein G. J. Mol. Biol. 1998, 278, 439− 456. (16) Sheinerman, F. B.; Brooks, C. L., III Molecular picture of folding of a small alpha/beta protein. Proc. Natl. Acad. Sci. U. S. A. 1998, 95, 1562−1567. (17) Karanicolas, J.; Brooks, C. L., III The origins of asymmetry in the folding transition states of protein L and protein G. Protein Sci. 2002, 11, 2351−2361. (18) Weikl, T. R.; Dill, K. A. Folding rates and low-entropy-loss routes of two-state proteins. J. Mol. Biol. 2003, 329, 585−598. (19) Brown, S.; Head-Gordon, T. Intermediates and the folding of proteins L and G. Protein Sci. 2004, 13, 958−970. (20) Kmiecik, S.; Kolinski, A. Folding Pathway of the B1 Domain of Protein G Explored by Multiscale Modeling. Biophys. J. 2008, 94, 726− 736. (21) Kamagata, K.; Arai, M.; Kuwajima, K. Unification of the Folding Mechanisms of Non-two-state and Two-state Proteins. J. Mol. Biol. 2004, 339, 951−965. (22) Scalley, M. L.; Yi, Q.; Gu, H.; McCormack, A.; Yates, J. R.; Baker, D. Kinetics of Folding of the IgG Binding Domain of Peptostreptoccocal Protein L. Biochemistry 1997, 36, 3373−3382. 3448


Article

Journal of Chemical Theory and Computation (44) Beeman, D. Some Multistep Methods for Use in Molecular Dynamics Calculations. J. Comput. Phys. 1976, 20, 130−139. (45) Zhang, Y.; Skolnick, J. Scoring function for automated assessment of protein structure template quality. Proteins: Struct., Funct., Genet. 2004, 57, 702−10. (46) Lee, J.; Gross, S. P.; Lee, J. Modularity optimization by conformational space annealing. Phys. Rev. E 2012, 85, No. 056702. (47) Morrone, A.; Giri, R.; Toofanny, R. D.; Travaglini-Allocatelli, C.; Brunori, M.; Daggett, V.; et al. GB1 Is Not a Two-State Folder: Identification and Characterization of anOn-Pathway Intermediate. Biophys. J. 2011, 101, 2053−2060. (48) Rocco, A. G.; Mollica, L.; Ricchiuto, P.; Baptista, A. M.; Gianazza, E.; Eberini, I. Characterization of the Protein Unfolding Processes Induced by Urea and Temperature. Biophys. J. 2008, 94, 2241−2251. (49) Zhao, L.; Wang, J.; Dou, X.; Cao, Z. Studying the unfolding process of protein G and protein L under physical property space. BMC Bioinf. 2009, 10, S44. (50) Arad-Haase, G.; Chuartzman, S. G.; Dagan, S.; Nevo, R.; Kouza, M.; Mai, B. K.; et al. Mechanical Unfolding of Acylphosphatase Studied by Single-Molecule Force Spectroscopy and MD Simulations. Biophys. J. 2010, 99, 238−247. (51) Kouza, M.; Hu, C.-K.; Li, M. S.; Kolinski, A. A structure-based model fails to probe the mechanical unfolding pathways of the titin I27 domain. J. Chem. Phys. 2013, 139, No. 065103. (52) Koga, N.; Takada, S. Roles of native topology and chain-length scaling in protein folding: A simulation study with a Go̅-like model. J. Mol. Biol. 2001, 313, 171−180. (53) Prieto, L.; de Sancho, D.; Rey, A. Thermodynamics of Go̅-type models for protein folding. J. Chem. Phys. 2005, 123, 154903. (54) Yoo, T. Y.; Adhikari, A.; Xia, Z.; Huynh, T.; Freed, K. F.; Zhou, R.; et al. The folding transition state of protein L is extensive with nonnative interactions (and not small and polarized). J. Mol. Biol. 2012, 420, 220−234. (55) Baxa, M. C.; Yu, W.; Adhikari, A. N.; Ge, L.; Xia, Z.; Zhou, R.; et al. Even with nonnative interactions, the updated folding transition states of the homologs Proteins G & L are extensive and similar. Proc. Natl. Acad. Sci. U. S. A. 2015, 112, 8302−8307. (56) Taketomi, H.; Ueda, Y.; Go̅, N. Studies on Protein Folding, Unfolding and Fluctuations by Computer Simulation. Int. J. Pept. Protein Res. 1975, 7, 445−459. (57) Abe, H.; Go̅, N. Noninteracting local-structure model of folding and unfolding transition in globular proteins. II. Application to twodimensional lattice proteins. Biopolymers 1981, 20, 1013−1031. (58) Takada, S. Go-ing for the prediction of protein folding mechanisms. Proc. Natl. Acad. Sci. U. S. A. 1999, 96, 11698−11700. (59) Schonbrun, J.; Dill, K. A. Fast protein folding kinetics. Proc. Natl. Acad. Sci. U. S. A. 2003, 100, 12678−12682. (60) Linhananta, A.; Boer, J.; MacKay, I. The equilibrium properties and folding kinetics of an all-atom Go model of the Trp-cage. J. Chem. Phys. 2005, 122, 114901. (61) Baweja, L.; Roche, J. Pushing the Limits of Structure-Based Models: Prediction of Nonglobular Protein Folding and Fibrils Formation with Go-Model Simulations. J. Phys. Chem. B 2018, 122, 2525−2535. (62) Sali, A.; Shakhnovich, E.; Karplus, M. How does a protein fold. Nature 1994, 369, 248−251. (63) Dinner, A. R.; Sali, A.; Karplus, M. The folding mechanism of larger model proteins: role of native structure. Proc. Natl. Acad. Sci. U. S. A. 1996, 93, 8356−8361. (64) Dinner, A. R.; Abkevich, V.; Shakhnovich, E.; Karplus, M. Factors that affect the folding ability of proteins. Proteins: Struct., Funct., Genet. 1999, 35, 34−40. (65) Munoz, V.; Eaton, W. A. A simple model for calculating the kinetics of protein folding from three-dimensional structures. Proc. Natl. Acad. Sci. U. S. A. 1999, 96, 11311−11316.

3449


Exploring the Folding Mechanism of Small Proteins GB1 and LB1

Recommend Documents