This Article Abstract Full Text (PDF) All Versions of this Article: biophysj.106.087684v1 92/2/547    most recent Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via Google Scholar Google Scholar Articles by Li, M. S. Articles by Hu, C.-K. Search for Related Content PubMed PubMed Citation Articles by Li, M. S. Articles by Hu, C.-K.

Refolding upon Force Quench and Pathways of Mechanical and Thermal Unfolding of Ubiquitin

Mai Suan Li ^*, Maksim Kouza ^* and Chin-Kun Hu 

^* Institute of Physics, Polish Academy of Sciences, Warsaw, Poland; and Institute of Physics, Academia Sinica, Nankang, Taipei, Taiwan

Correspondence: Address reprint requests to Prof. Mai Suan Li, E-mail: masli{at}ifpan.edu.pl; or Prof. Chin-Kun Hu, E-mail: huck{at}phys.sinica.edu.tw.

	ABSTRACT

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
The refolding from stretched initial conformations of ubiquitin (PDB ID: 1ubq) under the quenched force is studied using the C-G model and the Langevin dynamics. It is shown that the refolding decouples the collapse and folding kinetics. The force-quench refolding-times scale as _F exp(f_qx_F/k_BT), where f_q is the quench force and x_F 0.96 nm is the location of the average transition state along the reaction coordinate given by the end-to-end distance. This value is close to x_F 0.8 nm obtained from the force-clamp experiments. The mechanical and thermal unfolding pathways are studied and compared with the experimental and all-atom simulation results in detail. The sequencing of thermal unfolding was found to be markedly different from the mechanical one. It is found that fixing the N-terminus of ubiquitin changes its mechanical unfolding pathways much more drastically compared to the case when the C-end is anchored. We obtained the distance between the native state and the transition state x_UF 0.24 nm, which is in reasonable agreement with the experimental data.

	INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
Deciphering the folding and unfolding pathways and free energy landscape of biomolecules remains a challenge in molecular biology. Traditionally, folding and unfolding are monitored by changing temperature or concentration of chemical denaturants. In these experiments, due to thermal fluctuations of initial unfolded conformations it is difficult to describe the folding mechanisms in an unambiguous way. With the help of the atomic force microscopy, mechanical force has been used to prepare well-defined initial states of proteins (1,2). Using the initial force, f_I, which is higher than the equilibrium critical force, f_c, to unfold the tandem of poly ubiquitin (Ub), Fernandez and Li (2) have shown that the refolding can be initiated starting from stretched conformations or force-denaturated ensemble (FDE) and quenching the force to a low constant value, f_q (f_q < f_c). Monitoring folding events as a function of the end-to-end distance (R), they have made the following important observations:

1. Contrary to the standard folding from the thermal denaturated ensemble (TDE), the refolding under the quenched force is a multiple stepwise process.
   
2. The force-quench refolding time obeys the Bell formula (3), , where is the folding time in the absence of the quench force and x_F is the average location of the transition state (TS).
   

Motivated by the experiments of Fernandez and Li (2), we have studied (4) the refolding of the domain I27 of the human muscle protein using the C-G model (5) and the four-strand ß-barrel model sequence S1 (6) (for this sequence the nonnative interactions are also taken into account). Basically, we have reproduced qualitatively the major experimental findings listed above. In addition, we have shown that the refolding is a two-state process in which the folding to the native basin attractor (NBA) follows the quick collapse from initial stretched conformations with low entropy. The corresponding kinetics can be described by the biexponential time dependence, contrary to the single exponential behavior of the folding from the TDE with high entropy.

To make the direct comparison with the experiments of Fernandez and Li (2), in this article we performed simulations for a single domain Ub using the C-G model (see Materials and Methods for more details). Because the study of refolding of 76-residue Ub (Fig. 1 a) by all-atom simulations is beyond present computational facilities, the G modeling is an appropriate choice. Most of the simulations have been carried out at T = 0.85 T_F = 285 K. Our present results for refolding upon the force quench are in qualitative agreement with the experimental findings of Fernandez and Li, and with those obtained for I27 and S1 theoretically (4). A number of quantitative differences between I27 and Ub will be also discussed. For Ub, we have found the average location of the transition state x_F 0.96 nm, which is in reasonable agreement with the experimental value 0.8 nm (2).

View larger version (36K): [in this window] [in a new window]   FIGURE 1  (a) Native state conformation of ubiquitin taken from the PDB (PDB ID: 1ubq). There are five ß-strands: S1 (2–6), S2 (12–16), S3 (41–45), S4 (48–49), and S5 (65–71), and one helix A (23–34). (b) Structures B, C, D, and E consist of pairs of strands (S1,S2), (S1,S5), (S3,S5), and (S3,S4), respectively. In the text we also refer to helix A as structure A.

1. The barrier to the mechanical unfolding has not been computed.
   
2. Experiments of Schlierf et al. (7) have suggested that Cluster 1 (i.e., strands S1, S2, and the helix A) unfolds after Cluster 2 (strands S3, S4, and S5). However, this observation has not yet been studied theoretically.
   
3. Since the structure C, which consists of the strands S1 and S5, unzips first, Irbäck et al. (9) pointed out that strand S5 unfolds before S2 (the terminal strands follows the unfolding pathway S1 S5 S2). This conclusion may be incorrect, because it has been obtained from the breaking of the contacts within the structure C.
   
4. In pulling and force-clamp experiments, the external force is applied to one end of the proteins, while the other end is kept fixed. Therefore, one important question emerges as to how fixing one terminus affects the unfolding sequencing of Ub. This issue has not been addressed by Irbäck et al. (9).
   
5. Using a simplified all-atom model, it was shown (9) that mechanical intermediates occur more frequently than in experiments (7). It is relevant to ask if a C-G model can capture similar intermediates as this may shed light on the role of nonnative interactions.
   

In this article, from the force dependence of mechanical unfolding times we estimated the distance between the native state and the transition state to be x_UF 0.24 nm, which is close to the experimental results of Carrion-Vazquez et al. (10) and Schlierf et al. (7). In agreement with the experiments (7), Cluster 1 was found to unfold after Cluster 2 in our simulations. Applying the force to the both termini, we studied the mechanical unfolding pathways of the terminal strands in detail and obtained the sequencing S1 S2 S5, which is different from the result of Irbäck et al. (9). When the N-terminus is fixed and the C-terminus is pulled by a constant force, the unfolding sequencing was found to be very different from the previous case. The unzipping initiates, for example, from the C-terminus but not from the N-terminus. Anchoring the C-end is shown to have a little effect on unfolding pathways. We have demonstrated that the present C-G model does not capture rare mechanical intermediates, presumably due to the lack of nonnative interactions. Nevertheless, it can correctly describe the two-state unfolding of Ub (7).

It is well known that thermal unfolding pathways may be very different from the mechanical ones, as has been shown for the domain I27 (11). This is because the force is applied locally to the termini while thermal fluctuations have the global effect on the entire protein. In the force case, unzipping should propagate from the termini whereas under thermal fluctuations the most unstable part of a polypeptide chain unfolds first.

The unfolding of Ub under thermal fluctuations was investigated experimentally by Cordier and Grzesiek (12) and by Chung et al. (13). If one assumes that unfolding is the reverse of the refolding process then one can infer information about the unfolding pathways from the experimentally determined -values (14) and -values (15,16). The most comprehensive -value analysis is that of Went and Jackson. They found that the C-terminal region, which has very low -values, unfolds first and then the strand S1 breaks before full unfolding of the -helix fragment A occurs. However, the detailed unfolding sequencing of the other strands remains unknown.

Theoretically, the thermal unfolding of Ub at high temperatures has been studied by all-atom molecular dynamics (MD) simulations by Alonso and Daggett (17) and Larios et al. (18). In the latter work, the unfolding pathways were not explored. Alonso and Daggett have found that the -helix fragment A is the most resilient toward temperature but the structure B breaks as early as the structure C. The fact that B unfolds early contradicts not only the results for the -values obtained experimentally by Went and Jackson (14) but also findings from a high resolution NMR (12). Moreover, the sequencing of unfolding events for the structures D and E was not studied.

What information about the thermal unfolding pathways of Ub can be inferred from the folding simulations of various coarse-grained models? Using a semiempirical approach, Fernandez predicted (19) that the nucleation site involves the ß-strands S1 and S5. This suggests that thermal fluctuations break these strands last, but what happens to the other parts of the protein remain unknown. Furthermore, the late breaking of S5 contradicts the unfolding (12) and folding (14) experiments. From later folding simulations of Fernandez et al. (20,21), one can infer that the structures A, B, and C unzip late. Since this information is gained from -values, it is difficult to determine the sequencing of unfolding events even for these fragments. Using the results of Gilis and Rooman (22), we can only expect that the structures A and B unfold last. In addition, with the help of a three-bead model it was found (23) that the C-terminal loop structure is the last to fold in the folding process and most likely plays a spectator role in the folding kinetics. This implies that strands S4, S5, and the second helix (residues 38–40) would unzip first but again the full unfolding sequencing cannot be inferred from this study.

Thus, neither the direct MD (17) nor indirect folding simulations (19–23) provide a complete picture of the thermal unfolding pathways for Ub. One of our aims is to decipher the complete thermal unfolding sequencing and compare it with the mechanical one. The mechanical and thermal routes to the denaturated states have been found to be very different from each other. Under the force, e.g., the ß-strand S1, unfolds first, while thermal fluctuations detach strand S5 first. The later observation is in good agreement with NMR data of Cordier and Grzesiek (12). A detailed comparison with available experimental and simulation data on the unfolding sequencing will be presented. The free energy barrier to thermal unfolding was also calculated.

To summarize, in this article we have obtained the following novel results. We have shown that the refolding of Ub is a two-stage process in which the "burst" phase exists on very short timescales. The construction of the T – f phase diagram allows us to determine the equilibrium critical force f_c separating the folded and unfolded regions. Using the exponential dependence of the refolding and unfolding times on f, x_F and x_UF were computed. Our results for f_c, x_F and x_UF are in acceptable agreement with the experiments. It has been demonstrated that fixing the N-terminus of Ub has much stronger effect on mechanical unfolding pathways compared to the case when the C-end is anchored. In comparison with previous studies, we provide a more complete picture for thermal unfolding pathways, which are very different from the mechanical ones.

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
C-G model for Ub
We use coarse-grained continuum representation for Ub in which only the positions of C-carbons are retained. The interactions between residues are assumed to be G-like and the energy of such a model is (5)

(1)

Here _i = _i – _0i, r_{i, i+1} is the distance between beads i and i + 1, _i is the bond angle between bonds (i – 1) and i, and _i is the dihedral angle around the i^th bond and r_ij is the distance between the i^th and j^th residues. Subscripts 0, NC, and NNC refer to the native conformation, native contacts, and nonnative contacts, respectively. Residues i and j are in native contact if r_0ij is less than a cutoff distance d_c taken to be d_c = 6.5 Å, where r_0ij is the distance between the residues in the native conformation. With this choice of d_c and the native conformation from the PDB (Fig. 1 a), we have the total number of native contacts Q_max = 99.

The first harmonic term in Eq. 1 accounts for chain connectivity and the second term represents the bond-angle potential. The potential for the dihedral angle degrees of freedom is given by the third term in Eq. 1. The interaction energy between residues that are separated by at least three beads is given by 10–12 Lennard-Jones potential. A soft-sphere repulsive potential (the fourth term in Eq. 1) disfavors the formation of nonnative contacts. The last term accounts for the force applied to C- and N-termini along the end-to-end vector . We choose K_r = 100 _H/Ų, K = 20 _H/rad², , and K⁽³⁾ = 0.5 _H, where _H is the characteristic hydrogen bond energy and C = 4 Å. Since T_F = 0.675 _H (see below) and T_F = 332.5 K (24), we have _H = 4.1 kJ/mol = 0.98 kcal/mol. Then the force unit [f] = /Å = 68.0 pN.

We assume the dynamics of the polypeptide chain obeys the Langevin equation. The equations of motion (see (25) for details) were integrated using the velocity form of the Verlet algorithm (26) with the time step t = 0.005_L, where _L = (ma²/_H)^1/2 3 ps.

Simulations
To obtain the T – f phase diagram we use the fraction of native contacts or the overlap function (27)

(2)

where _ij is equal to 1 if residues i and j form a native contact and 0 otherwise, and (x) is the Heaviside function. The argument of this function guarantees that a native contact between i and j is classified as formed when r_ij is shorter than 1.2 r_0ij. The probability of being in the native state, f_N, which can be measured by various experimental techniques, is defined as f_N = , where ... stands for a thermal average. The T – f phase diagram (a plot of 1 – f_N as a function of f and T) and thermodynamic quantities were obtained by the multiple histogram method (28), extended to the case when the external force is applied to the termini (29,30). In this case, the reweighting is carried out not only for temperature but also for force. We collected data for six values of T at f = 0 and for five values of f at a fixed value of T. The duration of MD runs for each trajectory was chosen to be long enough to get the system fully equilibrated (9 x 10⁵_L from which 1.5 x 10⁵_L were spent on equilibration). For a given value of T and f, we have generated 40 independent trajectories for thermal averaging.

For the mechanical unfolding we have considered two cases. In the first case, the external force is applied via both termini N and C. In the second case it is applied to either N- or C-terminus.

To simulate the mechanical unfolding the computation has been performed at T = 285 K and mainly at the constant force f = 70, 100, 140, and 200 pN. This allows us to compare our results with the mechanical unfolding experiments (7) and to see if the unfolding pathways change at low forces. Starting from the native conformation but with different random number seeds the unfolding sequencing of helix A and five ß-stands is studied by monitoring fraction of native contacts as a function of the end-to-end extension. In the case of structures A, B, C, D, and E we consider not only the evolution of the number of intrastructure contacts as has been done by Irbäck et al. (9), but also the evolution of all contacts (intrastructure contacts and the contacts formed by a given structure with the rest of a protein).

In the thermal unfolding case the simulation is also started from the native conformations and it is terminated when all of the native contacts are broken. Due to thermal fluctuations there is no one-to-one correspondence between R and time. Therefore R ceases to be a good reaction coordinate for describing unfolding sequencing. To rescue this, for each i^th trajectory we introduce the progressive variable , where is the unfolding time. Then we can average the fraction of native contacts over a unique window 0 _i 1 and monitor the unfolding sequencing with the help of the progressive variable .

	RESULTS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS DISCUSSION ACKNOWLEDGEMENTS REFERENCES

 
Temperature-force phase diagram and thermodynamic quantities
The T – f phase diagram, obtained by the extended histogram method (see Materials and Methods), is shown in Fig. 2 a. The folding-unfolding transition, defined by the yellow region, is sharp in the low temperature region but it becomes less cooperative (the fuzzy transition region is wider) as T increases. The weak reentrancy (the critical force slightly increases with T) occurs at low temperatures. This seemingly strange phenomenon occurs as a result of competition between the energy gain and the entropy loss upon stretching. The similar cold unzipping transition was also observed in a number of models for heteropolymers (31) and proteins (29) including the C-G model for I27 (M. S. Li, unpublished results). As follows from the phase diagram, at T = 285 K, the critical force f_c 30 pN, which is close to f_c 25 pN, is estimated from the experimental pulling data (to estimate f_c from experimental pulling data, we use f_max f_cln(v/v_min) (32), where f_max is the maximal force needed to unfold a protein at the pulling speed v. From the raw data in Fig. 3 b of Carrion-Vasquez et al. (10), we obtain f_c 25 pN). Given the simplicity of the model this agreement can be considered satisfactory and it validates the use of the G model.

View larger version (36K): [in this window] [in a new window]   FIGURE 2  (a) The T – f phase diagram obtained by the extended histogram method. The force is applied to termini N and C. The color code for 1 – (T, f) is given on the right. The blue color corresponds to the state in the NBA, while the red color indicates the unfolded states. The vertical dashed line refers to T = 0.85, TF 285 K, at which most of simulations have been performed. (b) The temperature dependence of fN (open circles) defined as the renormalized number of native contacts (see Material and Methods). The solid line refers to the two-state fit to the simulation data. The dashed line represents the experimental two-state curve with Hm = 48.96 kcal/mol and Tm = 332.5 K (24).

View larger version (12K): [in this window] [in a new window]   FIGURE 3  (a) The dependence of the free energy on Q for selected values of f at T = TF. D and N refer to the denaturated and native states, respectively. (b) The dependence of folding and unfolding barriers, obtained from the free energy profiles, on f. The linear fits y = 0.36 + 0.218x and y = 0.54 – 0.029x correspond to FF and FUF, respectively. From these fits we obtain xF 1 nm and xUF 0.13 nm.

The cooperativity of the denaturation transition may be characterized by the cooperativity index, _c (see (34) and (35) for definition). From simulation data for f_N presented in Fig. 2 b we have _c 57, which is considerably lower than the experimental value _c 384 obtained with the help of H_m = 48.96 kcal/mol and T_m = 332.5 K (24). The underestimation of _c in our simulations is not only a shortcoming of the off-lattice G model (36) but also a common problem of much more sophisticated force fields in all-atom models (33).

Another measure of the cooperativity is the ratio between the van' t Hoff and the calorimetric enthalpy ₂ (37). For the G Ub we obtained ₂ 0.19. Applying the base line subtraction (38) gives ₂ 0.42, which is still much below ₂ 1 for the truly one-or-none transition. Since ₂ is an extensive parameter, its low value is due to the shortcomings of the off-lattice G models but not due to the finite size effects. More rigid lattice models give better results for the calorimetric cooperativity ₂ (39).

Fig. 3 a shows the free energy as a function of Q for several values of force at T = T_F. Since there are only two minima, our results support the two-state picture of Ub (7,13). As expected, the external force increases the folding barrier, F_F (F_F = F_TS – F_D) and it lowers the unfolding barrier, F_UF (F_UF = F_TS – F_N). From the linear fits in Fig. 3 b we obtain the average distance between the TS and D states, x_F = F_F/f 1 nm, and the distance between TS and the native state, x_UF = F_UF/f 0.13 nm. Note that x_F is very close to x_F 0.96 nm obtained from refolding times at a bit lower temperature T = 285 K (see Fig. 6 below). However, x_UF is lower than value 0.24 nm followed from mechanical unfolding data at f > f_c (Fig. 8). This difference may be caused by either sensitivity of x_UF to the temperature, or the determination of x_UF from the approximate free energy landscape, since a function of a single coordinate Q is not sufficiently accurate.

View larger version (6K): [in this window] [in a new window]   FIGURE 6  The dependence of folding times on the quench force at T = 285 K. The value F was computed as the average of the first passage times (F is the same as extracted from the biexponential fit for Q(t)). The result is averaged over 30–50 trajectories depending on fq. From the linear fit, y = 3.951 + 0.267x. With correlation level equal to –0.96, we obtain xF 0.96 nm.

View larger version (11K): [in this window] [in a new window]   FIGURE 8  Dependence of mechanical unfolding time on the force. Circles refer to the process when the force is applied to both N- and C-termini. Squares signify the case when the N-end is fixed and the C-end is pulled. For the first case the linear fit (y = 9.247 – 0.067x) gives the distance between the native state and TS xUF 0.24 nm. In the second case, from the linear fit (y = 9.510 – 0.062x) we obtained xUF 0.22 nm. Thus, within error bars, fixing one end does not affect the value of xUF. The inset shows the linear dependence of UF on f in the high force regime.

Refolding under quenched force
Our protocol for studying the refolding of Ub is identical to that used in the experiments of Fernandez and Li (2). We first apply the force f_I 70 pN to prepare initial conformations (the protein is stretched if R 0.8 L, where the contour length L = 28.7 nm). Starting from the FDE we quenched the force to f_q < f_c and then monitored the refolding process by following the time dependence of the number of native contacts Q(t), R(t), and the radius of gyration R_g(t) for typically 50 independent trajectories.

Fig. 4 shows considerable diversity of refolding pathways. In accord with experiments (2) and simulations for I27 (4), the reduction of R occurs in a stepwise manner. In the f_q = 0 case (Fig. 4 a), R decreases continuously from 18 nm to 7.5 nm (stage 1) and fluctuates around this value for 3 ns (stage 2). The further reduction to R 4.5 nm (stage 3) until a transition to the NBA. The stepwise nature of variation of Q(t) is also clearly shown up but it is more masked for R_g(t). Although we can interpret another trajectory for f_q = 0 (Fig. 4 b) in the same way, the timescales are different. Thus, the refolding routes are highly heterogeneous.

View larger version (43K): [in this window] [in a new window]   FIGURE 4  (a,b) The time dependence of Q, R, and Rg for two typical trajectories starting from FDE (fq = 0 and T = 285 K). Arrows 1, 2, and 3 in panel a correspond to time 3.1 (R = 10.9 nm), 9.3 (R = 7.9 nm), and 17.5 ns (R = 5 nm). Arrow 4 marks the folding time F = 62 ns (R = 2.87 nm) when all 99 native contacts are formed. Panels c and d are the same as in panels a and b, but for fq = 6.25 pN. The corresponding arrows refer to t = 7.5 (R = 11.2 nm), 32 (R = 9.4 nm), 95 ns (R = 4.8 nm), and F = 175 ns (R = 3.65 nm).

Fig. 5 shows the time dependence of R(t), Q(t), and R_g(t), where ... stands for averaging over 50 trajectories. The left and right panels correspond to the long and short time windows, respectively. For the TDE case (Fig. 5, a and b), the single exponential fit works pretty well for R(t) for the whole time interval. A little departure from this behavior is seen for Q(t) and R_g(t) for t < 2 ns (Fig. 5 b). Contrary to the TDE case, even for f_q = 0 (Fig. 5, c and d) the difference between the single and biexponential fits is evident not only for Q(t) and R_g(t) but also for R(t). The timescales, above which two fits become eventually identical, are slightly different for three quantities (Fig. 5 d). The failure of the single exponential behavior becomes more and more evident with the increase of f_q, as demonstrated in Fig. 5, e and f, for the FDE case with f_q = 6.25 pN.

View larger version (35K): [in this window] [in a new window]   FIGURE 5  (a) The time dependence of Q(t), R(t), and Rg(t) when the refolding starts from TDE. (b) The same as in panel a, but for the short timescale. (c,d) The same as in panels a and b, but for FDE with fq = 0. (e,f) The same as in panels c and d, but for fq = 6.25 pN.

Fig. 6 shows the dependence of the folding times on f_q. Using the Bell-type formula (3) and the linear fit in Fig. 6, we obtain x_F 0.96 nm, which is in acceptable agreement with the experimental value x_F 0.8 nm (2). The linear growth of the free energy barrier to folding with f_q is due to the stabilization of the random coil states under the force. Our estimate for Ub is higher than x_F 0.6 nm obtained for I27 (4). One of possible reasons for such a pronounced difference is that we used the cutoff distance d_c = 0.65 and 0.6 nm in the G model (1) for Ub and I27, respectively. The larger value of d_c would make a protein more stable (more native contacts) and it may change the free energy landscape leading to enhancement of x_F. This problem requires further investigation.

Absence of mechanical unfolding intermediates in C-G model
To study the unfolding dynamics of Ub, Schlierf et al. (7) have performed the AFM experiments at a constant force f = 100, 140, and 200 pN. The unfolding intermediates were recorded in 5% of 800 events at different forces. The typical distance between the initial and intermediate states is R = 8.1 ± 0.7 nm (7). However, the intermediates do not affect the two-state behavior of the polypeptide chain. Using the all-atom models, Irbäck et al. (9) have also observed the intermediates in the region 6.7 nm < R < 18.5 nm. Although the percentage of intermediates is higher than in the experiments, the two-state unfolding events remain dominating. To check the existence of force-induced intermediates in our model, we have performed the unfolding simulations for f = 70, 100, 140, and 200 pN. Because the results are qualitatively similar for all values of force, we present the f = 100 pN case only.

Fig. 7 shows the time dependence of R(t) for 15 runs starting from the native value R_N 3.9 nm. For all trajectories the plateau occurs at R 4.4 nm. As seen below, passing this plateau corresponds to breaking of intrastructure native contacts of structure C. At this stage, the chain ends are almost stretched out, but the rest of the polypeptide chain remains nativelike. The plateau is washed out when we average over many trajectories and R(t) is well fitted by a single exponential (Fig. 7), in accord with the two-state behavior of Ub (7).

View larger version (20K): [in this window] [in a new window]   FIGURE 7  Time dependence of the end-to-end distance for f = 100 pN. The thin curves refer to 15 representative trajectories. The average over 200 trajectory R(t) values is represented by the thick line. The dashed curve is the single exponential fit R(t) = 21.08 – 16.81 exp(–x/UF), where UF 11.8 ns.

Mechanical unfolding barrier
We now try to determine the barrier to the mechanical unfolding from the dependence of the unfolding times _UF on f. It should be noted that this way of determination of the unfolding barrier is exact and it would give a more reliable estimate compared to the free energy landscape approach in which the free energy profile is approximated as a function of only one order parameter.

We first consider the case when the force is applied via both termini N and C. Since the force lowers the unfolding barrier, _UF should decrease as f increases (Fig. 8). The present G model gives _UF smaller than the experimental values by approximately eight orders of magnitude. E.g., for f = 100 pN, _UF 12 ns whereas the experiments gives _UF 2.77 s (7). As seen from Fig. 8, for f < 140 pN _UF depends on f exponentially. In this regime, , where x_UF is the average distance between the N and TS states. From the linear fit in Fig. 8 we obtained x_UF 0.24 nm. Using different fitting procedures, Schlierf et al. (7) obtained x_UF 0.14 nm and 0.17 nm. The larger value x_UF 0.25 nm was reported in the earlier experiments (10). Thus, given experimental uncertainty, the C-G model provides a reasonable estimate of x_UF for the two-state Ub.

In the high force regime (f > 140 pN), instead of the exponential dependence, _UF scales with f linearly (inset in Fig. 8). The crossover from the exponential to the linear behavior is in full agreement with the earlier theoretical prediction (32). Similar crossover has been also observed (41) for the another G-like model of Ub but x_UF has not been estimated. At very high forces, _UF is expected to be asymptotically independent of f.

One can show that fixing one terminus of a protein has the same effect on unfolding times no matter whether the N- or C-terminus is fixed. Therefore, we show the results obtained for the case when the N-end is anchored. As seen from Fig. 8, the unfolding process is slowed down nearly by a factor of 2. It may imply that diffusion-collision processes (42) play an important role in the Ub unfolding. Namely, as follows from the diffusion-collision model, the time, required for formation (breaking) contacts, is inversely proportional to the diffusion coefficient, D, of a pair of spherical units. If one of them is idle, D is halved and the time needed to break contacts increases accordingly. Although fixing one end increases the unfolding times, it does not change the distance between the TS and the native state, x_UF (Fig. 8).

Mechanical unfolding pathways: force is applied to both termini
Here we focus on the mechanical unfolding pathways by monitoring the number of native contacts as a function of the end-to-end extension R R – R_eq, where R_eq is the equilibrium value of R. For T = 285 K, R_eq 3.4 nm. Following Schlierf et al. (7), we first divide Ub into two clusters. Cluster 1 consists of strands S1, S2, and the helix A (42 native contacts) and cluster 2, strands S3, S4, and S5 (35 native contacts). The dependence of fraction of intracluster native contacts is shown in Fig. 9 for f = 70 and 200 pN (similar results for f = 100 and 140 pN are not shown). In agreement with the experiments (7), Cluster 2 unfolds first. The unfolding of these clusters becomes more and more synchronous upon decreasing f. At f = 70 pN the competition with thermal fluctuations becomes so important that two clusters may unzip almost simultaneously. Experiments at low forces are needed to verify this observation.

View larger version (16K): [in this window] [in a new window]   FIGURE 9  The dependence of fraction of the native contacts on the end-to-end extension for Cluster 1 (solid lines) and Cluster 2 (dashed lines) at f = 70 pN and 200 pN. The results are averaged over 200 independent trajectories. The arrow points to the extension R = 8.1 nm.

To compare the mechanical unfolding pathways of Ub with the all-atom simulation results (9), we discuss the sequencing of helix A and structures B, C, D, and E in more detail. We monitor the intrastructure native contacts and all contacts separately. The later include not only the contacts within a given structure but also the contacts between it and the rest of the protein. It should be noted that Irbäck et al. have studied the unfolding pathways based on the evolution of the intrastructure contacts. Fig. 10 a shows the dependence of the fraction of intrastructure contacts on R at f = 100 pN. At R 1 nm, which corresponds to the plateau in Fig. 7, most of the contacts of C are broken. In agreement with the all-atom simulations (9), the unzipping follows C B D E A. Since C consists of the terminal strands S1 and S5, it was suggested that these fragments unfold first. However, this scenario may be no longer valid if one considers not only intrastructure contacts but also other possible ones (Fig. 10 b). In this case the statistically preferred sequencing is B C D E A, which holds not only for f = 100 pN but also for other values of f. If it is true, then S2 will unfold even before S5. To make this point more transparent, we plot the fraction of contacts for S1, S2, and S5 as a function of R (Fig. 11 a) for a typical trajectory. Clearly, S5 detaches from the core part of a protein after S2 (see also the snapshot in Fig. 11 b). So, instead of the sequencing S1 S5 S2 proposed by Irbäck et al., we obtain S1 S2 S5.

View larger version (16K): [in this window] [in a new window]   FIGURE 10  (a) The dependence of fraction of the intrastructure native contacts on R for structures A, B, C, D, and E at f = 100 pN. (b) The same as in panel a, but for all native contacts. The results are averaged over 200 independent trajectories.

View larger version (23K): [in this window] [in a new window]   FIGURE 11  (a) The dependence of fraction of the native contacts on R for strands S1, S2, and S5 (f = 200 pN). The vertical dashed line marks the position of the plateau at R 1 nm. (b) The snapshot, chosen at the extension marked by the arrow in a, shows that S2 unfolds before S5. At this point, all native contacts of S1 and S2 have already broken, while 50% of the native contacts of S5 are still present.

View larger version (19K): [in this window] [in a new window]   FIGURE 12  (a) The dependence of fraction of the native contacts on extension for A and all ß-strands at f = 70 pN. (b) The same as in panel a, but for f = 200 pN. The arrow points to R = 8.1 nm where the intermediates are recorded on the experiments (7). The results are averaged over 200 trajectories.

(3)

It should be noted that the sequencing (3) is valid in the statistical sense. In some trajectories, S5 unfolds even before S1 and S2 or the native contacts of S1, S2, and S5 may be broken at the same timescale (Table 1). From Table 1, it follows that the probability of having S1 unfolded first decreases with lowering f but the main trend (3) remains unchanged. One has to stress again that the sequencing of the terminal strands S1, S2, and S5 given by Eq. 3 is different from that proposed by Irbäck et al. (9), based on the breaking of the intrastructure contacts of C. Unfortunately, there are no experimental data available for comparison with our theoretical prediction.

From Fig. 13, we obtain the unfolding pathways

(4a)

(4b)

which are also valid for the other values of force (f = 70, 100, and 140 pN). Similar to the case when the force is applied to both ends, the structure C unravels first and the helix A remains the most stable. However, the sequencing of B, D, and E changes markedly compared to the result obtained by Irbäck et al. (9) (Fig. 10 a).

View larger version (36K): [in this window] [in a new window]   FIGURE 13  (a) The dependence of fraction of the intrastructure native contacts on extension for all structures at f = 200 pN. The N-terminus is fixed and the external force is applied via the C-terminus. (b) The same as in panel a, but for the native contacts of all individual ß-strands and helix A. The results are averaged over 200 trajectories. (c) A typical snapshot to show that S5 is fully detached from the core while S1 and S2 still have 50% and 100% contacts, respectively. (d) The same as in panel b, but the C-end is anchored and N-end is pulled. The strong drop in the fraction of native contacts of S4 at R 7.5 nm does not correspond to the substantial change of structure as it has only three native contacts in total.

C-terminus is fixed
One can show that unfolding pathways of structures A, B, C, D, and E remain exactly the same as in the case when Ub has been pulled from both termini (see Fig. 10). Concerning the individual strands, a slight difference is observed for S₅ (compare Fig. 13 d and Fig. 12). Most of the native contacts of this domain break before S₃ and S₄, except the long tail at extension nm due to high mechanical stability of only one contact between residues 61 and 65 (the highest resistance of this pair is probably due to the fact that among 25 possible contacts of S₅ it has the shortest distance |61 – 65| = 4 in sequence). This scenario holds in 90% of trajectories, whereas S₅ unravels completely earlier than S₃ and S₄ in the remaining trajectories. Thus, anchoring C-terminus has much less effect on unfolding pathways than in the case when the N-terminus is immobile.

It is worthwhile to note that, experimentally, one has studied the effect of extension geometry on the mechanical stability of Ub fixing its C-terminus (10). The greatest mechanical strength (the longest unfolding time) occurs when the protein is extended between N- and C-termini. This result has been supported by Monte Carlo (10) as well as MD (8) simulations. However, the mechanical unfolding sequencing has not been studied yet. It would be interesting to check our results on the effect of fixing one end on Ub mechanical unfolding pathways by experiments.

Thermal unfolding pathways
To study the thermal unfolding we follow the protocol described in Materials and Methods. Two-hundred trajectories were generated starting from the native conformation with different random seed numbers. The fractions of native contacts of helix A and five ß-strands are averaged over all trajectories for the time window 0 1. The unfolding routes are studied by monitoring these fractions as a function of . Above T 500 K, the strong thermal fluctuations (entropy-driven regime) make all strands and helix A unfold almost simultaneously. Below this temperature, the statistical preference for the unfolding sequencing is observed. We focus on T = 370 and 425 K. As in the case of the mechanical unfolding, Cluster 2 unfolds before Cluster 1 (results not shown). However, the main departure from the mechanical behavior is that the strong resistance to thermal fluctuations of Cluster 1 is mainly due to the stability of strand S2 but not of helix A (compare Fig. 14, c and d, with Fig. 12). The unfolding of Cluster 2 before Cluster 1 is qualitatively consistent with the experimental observation that the C-terminal fragment (residues 36–76) is largely unstructured while nativelike structure persists in the N-terminal fragment (residues 1–35) (45–47). This is also consistent with the data from the folding simulations (23) as well as with the experiments of Went and Jackson (14), who have shown that the -values 0 in the C-terminal region. However, our finding is at odds with the high -values obtained for several residues in this region by all-atom simulations (48) and by a semiempirical approach (19). One possible reason for high -values in the C-terminal region is the force fields. For example, Marianayagam and Jackson have employed the GROMOS 96 force field (49) within the GROMACS software package (50). It would be useful to check whether the other force fields give the same result.

View larger version (48K): [in this window] [in a new window]   FIGURE 14  (a) The dependence of fraction of intrastructure native contacts on the progressive variable for all structures at T = 425 K. (b) The same as in panel a, but for T = 370 K. (c) The dependence of the all native contacts of the ß-strands and helix A at T = 425 K. (d) The same as in panel c, but for T = 370 K.

We now consider the thermal unstability of individual ß-strands and helix A. At T = 370 K (Fig. 14 d), the trend that S2 unfolds after S4 is more evident compared to the T = 425 K case (Fig. 14 c). Overall, the simple G model leads to the sequencing given by:

(5a)

(5b)

From Eqs. 3, 4b, and 5b