Mechanical Unfolding of Single Polyubiquitin Molecules Reveals Evidence of Dynamic Disorder.

Mechanical unfolding of single polyubiquitin molecules subjected to a constant stretching force showed nonexponentiality in the measured probability density of unfolding (waiting time distribution) and the survival probability of the folded state during the course of the measurements. These observations explored the relevance of disorder present in the system under study with implications for a static disorder approach to rationalize the experimental results. Here, an approach for dynamic disorder is presented based on Zwanzig's fluctuating bottleneck (FB) model, in which the rate of the reaction is controlled by the passage through the cross-sectional area of the bottleneck. The radius of the latter undergoes stochastic fluctuations that in turn is described in terms of the end-to-end distance fluctuations of the Rouse-like dynamics using a non-Markovian generalized Langevin equation with a memory kernel and Gaussian colored noise. Our results are comprised of analytical expressions for the survival probability and waiting time distribution, which show excellent agreement with the experimental data throughout the range of the applied forces. In addition, by fitting the survival probabilities at different stretching forces, we quantify two system parameters, namely, the average free energy ΔG av and the average distance to the transition state Δx av, both perfectly recovered the experimental estimates. These agreements validate the present model of polymer dynamics, which captures the very essence of dynamic disorder in single-molecule pulling experiments.

Introduction

Mechanical flexibility of protein molecules is a key component for their biological functioning.[1] This intrinsic property can be made accessible experimentally using single-molecule force spectroscopy, which shed much light on the extraction of inherent dynamical information about the system under study.[2−8] When a mechanical force is applied at one end of a single protein molecule, its length is increased, but the height of the activation energy barrier for the unfolding process is lowered by an amount equal to the magnitude of the external force times the distance between the native conformation and the transition state conformation. Interestingly, such a system manifests the significant effect of inherent disorder, which results in deviation from the Arrhenius behavior, typically showing nonexponentiality in the measured survival probability of unfolding. Here, in the present work, we consider one such example of single-molecule pulling experiment by Kuo et al.[6] where using force-clamp spectroscopy, a single polyubiquitin molecule was stretched by a constant force. We develop a theory and show the consequences of dynamic disorder on the unfolding kinetics otherwise treated originally with a static disorder approach.

During the course of the measurements, Kuo et al. observed the sequential unfolding of nine individual repeat units of the polyubiquitin molecules with a stepwise increment of ∼20 nm in the length while varying the magnitude of the external force in the region between 90 and 190 pN. The durations of the discrete unfolding events were not fixed, but spanned a range of timescales from a fraction of a second to several seconds. A histogram of 2799 unfolding events (probability density of unfolding) at a force of 110 pN clearly showed deviation from the single exponential decay at short dwell times. The justification for the experimental results was made on the basis of the heterogeneity of the population such that the different polyubiquitin molecules that adopt a particular conformation, distinct from each other, during the experiments have to surmount different barrier heights for the unfolding. The combined results from a set of large number of experiments do not follow the simple Arrhenius kinetics, and the measured survival probabilities exhibited nonexponential behavior. To analyze the experimental data, Kuo et al. presented a static disorder model based on Zwanzig’s work[9] and derived an expression for the ensemble-averaged survival probability with a force-modified rate constant. The latter depends on a random variable r that modulates the height of the barrier during the unfolding processes and follows a Gaussian distribution. Although the expression for the survival probability lacks a closed form, it could very well reproduce the experimental data numerically adjusting the mean and variance of the barrier height for the best fit. By fitting the logarithm of rate of crossing the barrier height as a function of applied force, the system parameters Δxav and ΔGav were estimated as 0.23 nm and 51 kJ mol–1 (85.1 pN nm), respectively. The authors deduced a closed-form expression for the survival probability in the limit of lower r. Although the resulting expression retains the nonexponential nature, at short times, it is well approximated by an exponential function,[10] not consistent with the experimental results. This therefore suggests that the emergence of the concept of dynamics disorder in the unfolding event would be quite natural. In this context, it is to be noted that both static disorder and dynamic disorder approaches lead to the nonexponential survival probability and therefore reproduce the experimental data when the relevant fit parameters are adjusted suitably.

In the paper by Chatterjee et al.,[10] the authors introduced the time-dependent nature of the random variable r whose stochastic evolution was governed by a generalized Langevin equation (GLE) subjected to fractional Gaussian noise with a power-law memory function.[11] The detailed calculations, invoking the Wilemski–Fixman approximation,[12,13] yielded analytical expressions for the probability density in short-time and long-time limits, which satisfied the experimental results separately in those limits with different fit parameters. However, the force dependence of the survival probability was not discussed and the system parameters were not estimated. These were accomplished by Zheng et al.[14] applying the Kramers’ rate theory[15] to the polyubiquitin unfolding. The authors considered a force-modified free energy surface, and the unfolding process was described by the passage of a single particle over it along the reaction coordinate. The latter was associated with the distance fluctuations in the protein molecule rather than the phenomenological variable r related to the fluctuations in the barrier height for the unfolding. The time evolution of the reaction coordinate was described by a similar generalized Langevin equation (GLE) presented in ref (10). The analytical results for the survival probability and waiting time distribution provided satisfactory comparisons with the experimental data. An estimate of the system parameter Δxav, average extension, as 0.26 nm was reasonable relative to the estimate of 0.23 nm by Kuo et al. However, ΔGav was calculated as 24.8 kJ mol–1, which deviates from the estimated value of 51 kJ mol–1. Moreover, a knowledge of the harmonic frequencies in well and barrier regions in the free energy profile was a prerequisite without which ΔGav cannot be estimated directly, even by fitting the experimental data.

In another study, Hyeon et al.[16] adopted the fluctuating bottleneck (FB) model, first introduced by Zwanzig to account for dynamic disorder.[17] In this model, the rate of the reaction was controlled by the passage through the cross-sectional area of the bottleneck and the reaction sink was taken proportional to the latter. The fluctuating radius of the bottleneck characterized the intrinsic dynamics of the protein, and the force dependence of the rate constant (fluctuation-independent rate) was assumed to obey the Bell approximation.[18] The resulting expression for the survival probability was essentially the same as provided by Zwanzig (eq of ref (17)). The comparison of the analytical results with the experimental data showed excellent agreement between the two. Also, in accordance with the Bell model, the quantification of Δxav gives a near-perfect result with Kuo et al. However, the average free energy was not estimated from their analysis and the microscopic nature of the conformational fluctuations is obscured. To this end, we would like to point out that by fitting the different experimental data, Hyeon et al. found a rough exponential relationship between the frequency of conformational transitions (λ), governing the internal dynamics, and the applied tension (f). As discussed by the authors, the rate of change of the cross-sectional area increases with higher λ as a consequence of larger f. Therefore, one may expect the possible influences of the applied forces on the protein dynamics, which will in turn tune the timescales of relaxation. In contrast, the relaxation time was assumed to be constant under the action of different stretching forces in the work by Zheng et al.

In this work, we revisit the nonexponential unfolding kinetics with a motivation of exploring the microscopic nature of the intrinsic dynamics of protein as well as the influence of an applied force on the dynamic characteristics. We account for the general platform of Zwanzig’s theory of fluctuating bottleneck and borrow the new perspective into the old problem of unfolding describing the time-varying radius of the bottleneck, R(t), by the end-to-end distance fluctuations of a Rouse chain that characterize the intrinsic dynamical nature of protein’s conformations, which in turn modulate the rate of the reaction.[19] The previous dynamic disorder studies are devoid of such an idea. In Zwanzig’s theory, the escape rate is the equilibrium flux through the bottleneck, which is proportional to the cross-sectional area of the bottleneck undergoing stochastic fluctuations due to the fluctuations in protein’s conformations. The application of a polymer dynamics model is a novel idea from the perspective of the simplest description of a protein, sequence of amino acids, which can effectively be described by a coarse-grained bead spring model. Previous studies showed that the protein molecules have a high degree of conformational flexibility similar to simple polymers.[21,22] Also, the Rouse model describes the dynamics of a flexible polymer chain. Therefore, the present model offers a more realistic approach and is much less coarse-grained compared to the traditional works where the description of the polyubiquitin molecule goes to the reduced level of a single particle. The immediate consequence of this difference will be discussed later. The stochastic evolution of R(t) occurs according to a non-Markovian generalized Langevin equation with a memory kernel and Gaussian colored noise.[23,24] During the course of the reaction, the folded polyubiquitin molecules unfold marked by a larger separation of the two ends. A notable point in this context is that R(t) in our model does not coincide with the physical distance between the two ends of the polyubiquitin molecule, rather the microscopic origin of the nonexponential kinetics due to the presence of dynamic disorder in the reaction pathway finds its association with the dynamics of R(t). We derive analytical expressions for the force-dependent survival probability and waiting time distribution from the solution of a non-Markovian reaction−diffusion equation using the Wilemski–Fixman approximation. Our results provide the most accurate reproduction of the experimental results along with perfect estimates of the system parameters Δxav and ΔGav. Also, we discuss at length the possible influence of the applied forces on the intrinsic conformational dynamics.

The organization of the paper is as follows. We discuss our theoretical model of dynamic disorder in Dynamic Disorder Model of Single Polyubiquitin Unfolding Kinetics section, where the key steps involved in the calculations lead to the development of a non-Markovian reaction−diffusion equation. The evaluations of the survival probability of the folded state and the probability density of unfolding from the solution of the above equation are also presented in the same section. Next, in Comparison with Experiment section, we show the comparisons of our theoretical results with the experimental data of protein unfolding from ref (6), followed by the Conclusions section.

Dynamic Disorder Model of Single Polyubiquitin Unfolding Kinetics

The ubiquitous presence of nonexponential kinetics is a characteristic feature of biological reactions. Numerous examples are available in the literature, including ligand binding of heme proteins,[25−28] electron transfer kinetics,[29,30] enzyme catalysis,[31−33] escape kinetics through nanopores,[34−36] etc. The origin of such behavior is rooted in the inherent disorders of the systems that are divided into static disorder and dynamic disorder according to Zwanzig. The protein molecules possess a large number of conformational substates. If the rate of interconversions between them is much slower than the rate of the reaction, a situation for the static disorder prevails. On the other hand, dynamic disorder refers to a circumstance when the rate of conformational transitions is comparable to or slower than the rate of the reaction. Explicit evidences of dynamic disorder have been reported in the literature.[37,38] In the present study, we reveal the effects of dynamic disorder on the reaction kinetics taking into account the representative example of polyubiquitin unfolding under the action of a stretching force. Though the polyubiquitin molecules reside in the folded state before unfolding, they may undergo conformational fluctuations during the measurements. However, there are a number of proteins that exhibit conformational fluctuations in the folded state.[39] Therefore, it is not surprising that polyubiquitin is also subjected to conformational dynamics being in the same state. The coarse-grained Rouse description serves as a natural model of conformational dynamics if there exists considerable structural flexibility within the molecule. However, the conformational flexibility is associated with the intersegment distance fluctuations, most often described by the end-to-end distance fluctuations. In the following, we implement this coarse-grained description to depict the effects of dynamic disorder associated with the pulling of single polyubiquitin molecules in the following subsections.

Rate of Force-Induced Unfolding

Assuming the unfolding process of a single polyubiquitin molecule to be irreversible, a first-order rate equation was proposed in refs (6) and (10), where the rate coefficients were coupled with a random variable associated with the fluctuations of barrier height. The fluctuation-independent rate was proportional to the force-modified free energy, given by kF ∝ exp[−β(ΔGav – F Δxav)], where ΔGav is the average barrier height for the unfolding process in the absence of the stretching force F, Δxav is the average distance to the transition state, and β = 1/kBT, where kB is the Boltzmann constant and T is the temperature. Following the same, we start our calculations with a similar first-order rate expression for the evolution of the survival probability, given bywhere R(t) represents the fluctuating radius of the bottleneck described by the time-varying end-to-end distance of a Rouse chain. The rate coefficient depends on R(t), which refers to the nonexistence of a well-defined rate constant of the reaction, characteristic of a system with dynamic disorder. Also, the notation of the survival probability signifies its dependent nature on the random variable and the external force. In conventional Bell’s formula,[18] the effect of external perturbation is included in the Arrhenius law, but there exists no coupling with the inherent dynamics. The latter was incorporated by Hyeon et al. borrowing the idea of fluctuating bottleneck from Zwanzig’s work in which the time-dependent fluctuations of the area of the bottleneck account for inherent conformational dynamics.

The salient idea of our work is to implement a polymer dynamics model that describes the conformational fluctuations of protein molecule during the reaction. The characteristic of dynamic disorder, having similar timescales for these simultaneous processes, requires that the rate expression should depend on the parameter describing the dynamical process. However, unlike electron transfer and energy transfer rate expressions, which carry the explicit distance dependence, there is no simple way to consider the latter for the unfolding rate expression. In this circumstance, Zwanzig’s idea of fluctuating bottleneck model gives a way to couple the conformational dynamics with the rate of the reaction. Therefore, in our work, Zwanzig’s FB model represents the general mechanism for dynamic disorder and we describe the rate coefficient given in eq with the following expressionwhere the fluctuation-independent rate can be expressed asIn this context, our approach is similar to Hyeon et al. Nevertheless, we assert a microscopically realized model as the bottleneck radius represents a distance between the two segments of the polymer chain, which adequately captures the simplest description of a protein. In eq , k′ is the preexponential factor, which will be combined later on with the other parameters, shown in Survival Probability of Unfolding subsection, to give an effective parameter, which is to be adjusted for the best fit. Also, ΔGav and Δxav in eq have their usual meaning. The rate expression of eq with kF given by eq is physically consistent as (i) an increase of F lowers the barrier height for the unfolding process, thereby increasing kF and (ii) the rate of the unfolding reaction increases with the increase of accessibility of the area, proportional to R(t)2, by the solvent molecules. At a given temperature when the molecule is stretched with a constant force, the modified average free energy (expression inside the exponential term) becomes a constant quantity, given by A = β(ΔGav – FΔxav). Kuo et al. used a preexponential factor of 106 s–1 to determine the system parameters from the Arrhenius equation between kF and F. Assigning an appropriate estimate for the effective parameter, a plot of A vs F would be linear, from which ΔGav and Δxav can be estimated. However, to accomplish this, we require a set of A values in the range of the applied forces. The single-molecule pulling experiments measured the survival probability S(t) of the folded state during an unfolding event under the action of a constant stretching force. Therefore, by fitting S(t) at different F’s, one can extract the dimensionless energy parameter A. Hence, it is of primary interest to determine the force-modified survival probability. The rest of the present section is devoted for the same.

Conformational Fluctuations of Single Polyubiquitin Molecules

Dynamic disorder in the reaction pathway lies in the conformational fluctuations of the protein molecules. The rate of fluctuations for a biomolecule like protein depends on its domain size. Larger the domain size, smaller is the rate of fluctuations and vice versa. However, it is to be noted that the timescales for those fluctuations which are comparable to the timescales of the reaction contribute significantly when dynamic disorder prevails. In the present work, we consider the distance fluctuations of an entire polymer chain, characterized by the end-to-end distance fluctuations. We parametrized the latter by R(t) = R0 + u(t), where R0 is the position vector of the nth monomer of the polymer chain and u(t) is the deviation from it. Therefore, the conformational dynamics can effectively be described by the stochastic dynamics of the end-to-end vector R(t) = R(t) – R0(t) = R0 + u(t) – u0(t).[23] The scalar component of R(t), given by , represents the radius of the bottleneck. The polyubiquitin molecules undergo conformational fluctuations during the unfolding process. The time evolution of R(t) is described by an overdamped non-Markovian generalized Langevin equation with a memory kernel and Gaussian colored noise, given by[23]The first term on the right-hand side is due to the elastic force originating from the chain connectivity. In the coarse-grained description, ⟨R2⟩ = Nb2, with N being the number of monomers of size b.[40]f(t) is the Gaussian colored noise characterized by ⟨f(t)⟩ = 0 and ⟨f(t)f(0)⟩ = kBTK(t), where K(t) is the friction kernel. Multiplying both sides of eq by R(0), averaging followed by the Laplace transform, a closed-form expression for the memory kernel can be derived, given bywhere χ(s) is the Laplace transform of the function χ(t), defined as χ(s) = ∫0∞ dt exp(−st)χ(t), where χ(t) is the normalized end-to-end distance autocorrelation function, given by . The Rouse dynamics yield the following expression for χ(t)where τ is the longest relaxation time associated with the first normal mode (p = 1). Interestingly, from eqs and 6, one can ensure the power-law behavior of the memory kernel, K(t) ∝ t–1/2, observed experimentally.[20] The details of the steps are given elsewhere.[23]

Before proceeding further, we would like to point out that from eq , one can clearly see the physical basis for the power-law behavior of the memory kernel, which lies in the intrinsic conformational fluctuations. Equation reveals that a microscopic polymer dynamics model suggests the time correlation of distance fluctuations, given by χ(t) in eq , is non-Markovian in nature (sum over the exponential terms). This arises because the end-to-end vector in Rouse-like dynamics is related to the sum of the normal modes.[23] The K(t) ∝ t–1/2 behavior is a natural consequence of this model, which does not require any beforehand input from the experimental results or an ad hoc power-law assumption. On the other hand, it might also be interesting to see the direct comparison of eq with the data from ref (20). However, in the present experimental context, this is not a necessary one.

Consideration of the protein molecule as a Rouse chain requires much less coarse graining relative to the previous traditional description.[10,14] Proteins are complex molecules formed by folded strands of amino acids. However, having complicated quaternary structure, a protein molecule also possesses a high degree of conformational flexibility similar to simple polymers. The simplest description of a protein is the sequence of amino acids connected by covalent peptide bonds. These two key features, conformational flexibility and sequence of monomers, suggest that a reduced description of the protein molecule may rely on the coarse-grained Rouse chain, which in turn describes the dynamics of a flexible polymer. However, it is to be noted that the complete description of a protein molecule and its detailed dynamics are beyond the scope of the coarse-grained Rouse chain, which may effectively be studied using the tools of biomolecular simulations. Nevertheless, it is important to see that our simplest description will adequately capture the very essence of the experimental observations, and thus the purpose of an analytical model is well accounted for by the Rouse description. Within this model, the natural emergence of the power-law behavior manifests the immediate consequence of choosing the Rouse-like dynamics over the one-dimensional GLE approach presented in refs (10) and (14). The latter considered a power-law memory kernel on ad hoc basis, which depends on the numerical choice of the Hurst index[10] or the related exponent γ.[14] Also, as shown in our previous work, ref (23), the Rouse-like dynamics is able to recover the expected δ-correlated Markovian behavior at long times. The said GLE approach however fails to find the same. Therefore, the generality of the present friction memory kernel is far-reaching.

The formal solution of eq can be written asEquation shows that the survival probability is a functional of R (also R is a functional of f(t)), and therefore it is required to average it over the distribution of f(t). This is usually done following a method known as Zwanzig’s indirect approach of noise averaging.[17] For this purpose, an equation for the joint probability density of S(t) and R(t) is needed, in which f(t) has been averaged out. If P(S,R,t) = ⟨δ(S – S(t))δ(R – R(t))⟩ denotes the probability density that at time t, S(t) and R(t) are given by S and R, respectively, the non-Markovian generalized Langevin equation, eq , can be transformed into the following Smoluchowski equation. The details of the steps are outlined in our earlier work[23] and ref (41)where in which χ(t) is given by eq .

When eq is multiplied by S and integrated over all S from 0 to 1, one gets the equation for the noise average survival probability, S̅(R,t) ≡ ∫01 dSSP(S,R,t), the same is given bywhere the operator D is defined asIt is to be noted that the structure of eq , reaction-diffusion equation, is essentially the same as the result of Zwanzig.[17] However, the difference is accompanied by the explicit nature of the fluctuation-independent rate coefficient (eq ) and the identity of the control parameter with the end-to-end distance of a Rouse chain. In ref (17), the corresponding equation was solved using an exponential form of the noise-averaged concentration of the ligand, while our method is based on an approximation technique described below. On the other hand, a recent work suggests the exact solution to the problem of a quadratic sink for a subdiffusive Brownian oscillator.[42] Although the result is exact, the method relies on the similar ad hoc power-law memory function, the limitations of which are discussed above and therefore away from an improved realistic description.

Survival Probability of Unfolding

If eq is solved and integrated over R, one obtains the expression for the average survival probability ⟨S(t)⟩ ≡ ∫dRS̅(R,t), which is compared with the experimental results. Equation is a nonlinear integral equation which cannot be solved in closed form. Therefore, we use an approximate method, known as the Wilemski–Fixman approximation,[12,13] to get the required expression. The basic idea of this approximation is given by S̅(R,t) – S̅eq(R)B(t), where S̅eq(R) = S̅(R,0), which describes the system in thermal equilibrium at t = 0. With this initial condition, the formal solution of the reaction-diffusion equation (eq ) can be written as[10,12−14,43,44]where G(R, t – t′|R′, 0) is Green’s function, which satisfies the following equation, given byThe solution of eq is known aswhich in the limit of t → ∞ defines S̅eq(R) given belowTo determine the survival probability, we substitute the Wilemski–Fixman approximation into the right-hand side of eq , which on subsequent integration over R yieldsConsidering and in the above equation, the subsequent calculations lead towhere k1 = k′Nb2 and A = β(ΔGav – F Δxav). On the other hand, if both sides of eq are multiplied by k(R) and then integrated over R, we get the following expressionwhereandThe Laplace transform of eq results inB(s) in eq can further be substituted with an expression obtained from the Laplace transform of eq , given byThis yields the following expression for the survival probability in the Laplace domainThe explicit expression for ⟨S(s)⟩ can be deduced when ka and C(t) are evaluated from their definitions. Substituting eqs , 3, and 14 into the definition of ka, the result of the integration givesOn the other hand, using eqs , 3, 13, and 14, C(t) can be evaluated asInserting eq and the Laplace transform of eq into eq , the latter appears asHere, we note that explicit closed-form expressions for ⟨S(t)⟩ cannot be obtained because an exact form of Laplace transform for Green’s function is not available owing to the presence of the correlation function χ(t). Therefore, performing the inverse Laplace transform of eq numerically, we obtain ⟨S(t)⟩, which gives the probability of the polyubiquitin molecule to be in the folded state that survives up to time t under the action of an applied tension. This quantity is directly compared with the measured data from ref (6) by adjusting the fit parameters to have the best agreement between our theory and the experiment. The nonexponentiality in the survival probability results from the appearance of the correlation function, the latter being a multiexponential one.

The waiting time distribution can be calculated from the expression for the survival probability by the following relationWe will compare f(t) numerically with the experimental results at a force of 110 pN using the same fit parameters required to satisfy the S(t) data at 110 pN. Additionally, we determine the mean time of unfolding of the polyubiquitin molecules from the definition

Comparison with Experiment

When a single polyubiquitin molecule is stretched at a constant force in atomic force microscopy (AFM), the individual unfolding events are characterized by a dwell time. For an applied force of 110 pN, the histogram of unfolding events showed the nonexponential nature of the waiting time distribution characterized by the failure of the single exponential fit at short dwell times. The authors determined the force-modified survival probability from the dwell time histogram, which also exhibited the pronounced nonexponential behavior. Similar results were obtained at different magnitudes of the stretching forces. We are now in the position to use our theoretical results to interpret the sets of experimental data from the force-clamp experiments of Kuo et al.[6]

We first compare f(t), given by eq , with the experimental results at 110 pN in Figure a, where the symbols are the data points and the solid line is the theoretical fit. One can see the good agreement between our theory and experiment in the entire timescale of the measurements. Unlike ref (10), we use a single expression for f(t) to satisfy the observed behavior. The corresponding fit parameters used in Figure a are listed in Table . The quality of the fit carries a signature of the presence of dynamic disorder during the unfolding events. It is interesting to note that the distribution becomes exponential at long times (Figure C of Kuo et al.). In this limit, the different Rouse modes, characterized by the mode index p, are relaxed and the nature of the correlation functions becomes exponential. We also calculate the waiting time distribution at other magnitudes of the stretching forces. However, histograms at other constant forces, 90 and 130–190 pN, unlike survival probability, were not available from ref (6) that would allow us to make a direct comparison of our results with the experimental data. Therefore, we determine f(t) using the parameters given in Table , which provides the best agreement with the measured data for survival probability, as discussed in the next paragraph. The time variation of the normalized waiting time distribution, f(t)/f(0), is shown in Figure b. We recover the expected trend for the same as a function of the applied force over the timescale of the measurements.

Figure 1

Interpretation of the experimental data of polyubiquitin unfolding from ref (6). (a) Measured data of the probability density of unfolding at 110 pN (circles), digitized from Figure 1C of ref (6), are compared with the waiting time distribution (eq ) using a logarithmic scale. The fit parameters are given in Table . The inset shows the corresponding linear plot. (b) Normalized waiting time distribution plotted as a function of time at different forces in the range of 90–190 pN using the parameters given in Table .

Table 1

Force-Dependent Parameters A and τ for the Best Fits of the Experimental Results from ref (6)a

force (pN)	A	τ (s)
90	15.65	31.0
110	14.01	4.8
130	13.40	4.0
150	12.33	1.9
170	10.74	0.43
190	9.96	0.12

The effective parameter k1 is set to be 106 s–1.

The survival probability of the folded polyubiquitin molecules, given by the inverse Laplace transform of eq , is compared with the experimental estimates, shown in Figure , measured at six different forces in the range of 90–190 pN. As before, the circles are the data points and the solid lines represent our theoretical results. We see excellent fits to the data at each of the stretching forces throughout the timescale of the experiments. It is important to mention here that an approach of dynamic disorder is a better alternative to explain the nonexponential kinetics as the expression for the survival probability depends on the correlation function χ(t), which essentially prevents S(t) to behave exponentially at the early times, a limitation of the static disorder model discussed earlier. The extracted parameters to satisfy the best agreement can be found in the tabulated results. In particular, we need to vary the dimensionless energy parameter A and the relaxation time τ keeping the effective parameter k1 constant. The analogue of k1, the preexponential factor, was assigned a value of 106 s–1 to achieve the best results in ref (6). We also set the value of k1 to be 106 s–1 to reproduce the data. In reality, a number of noncovalent interactions play a dominant role in a folded protein. However, the details of these interactions and how they guide the unfolding of a folded protein is beyond the scope of a simple coarse-grained analytical theory. Importantly, this does not create any serious problem in explaining the experimental results, as observed in the present work. Nevertheless, we note that studies based on molecular dynamics simulations are rather appropriate to address those molecular details.

Figure 2

Survival probability data (circles) at different forces, 90 pN (left) and 110–190 pN (right), digitized from Figure 2 of ref (6), are compared with the theoretical results (solid lines) given by the inverse Laplace transform of eq . The corresponding fit parameters are listed in Table .

Given the results of our work, we would like to point out that the use of general platforms like Zwanzig’s idea of fluctuating bottleneck and Wilemski–Fixman approximation significantly improved the quality of the work over the previous studies.[10,14] In particular, although ref (10) applied the Wilemski–Fixman approximation, the use of dynamic disorder approach based on single particle dynamics seems to be very limited in explaining the experimental results of force-dependent survival probability and system parameters. On the other hand, the nonlinear fit suggested in Figure of ref (14) is rather poor to satisfy the extracted data of the force-dependent parameter C(F). Even though one cannot estimate the average barrier height directly unless an approximation regarding the frequencies in the well and barrier regions is employed. Our work overcomes all these limitations and produces most accurate estimates.

Figure 4

Variation of the dimensionless energy parameter A, given in Table , with the applied force F. The linear relationship yields an estimate of the average barrier height in the absence of stretching force (ΔGav) as well as the average distance to the transition state (Δxav).

The variation in the values of A as a function of F is physically consistent because with increasing the magnitude of the applied force, the corresponding barrier height for the unfolding process decreases, resulting in an easy extension of the protein molecule. The qualitative nature of unfolding remains the same in the range of the applied forces as can be seen from the shapes of the different decay curves. On the other hand, in contrast to the work by Zheng et al., it was essential in the present study to alter the relaxation times at different stretching forces to achieve the best comparisons. The observed trend is that with increasing F, the relaxation time becomes smaller. We justify the same with the following discussion. During an unfolding event, the polyubiquitin molecule undergoes simultaneous conformational fluctuations described by the end-to-end dynamics of a Rouse chain. The relaxation time for the end-to-end fluctuation can be defined as the time required for the polymer to diffuse through a distance of its own size. In the Rouse model, both the diffusion and relaxation have the same origin and the diffusion coefficient is inversely related to the relaxation time.[45] As a result of applied tension, the molecule experiences enhanced diffusion, which lowers the relaxation time. The higher the magnitude of the stretching force, the lower is the relaxation time. Our results are supported by the work of Hyeon et al. in a qualitative manner, where the frequency of conformational transitions was enhanced with the applied force and a rough exponential relationship was established between them. Interestingly, the authors also revealed similar observations of enhancement of the transition frequency with the pulling speed from the force-ramp experiments.

An extension of the preceding discussion emphasizes the finding of a possible relationship that might exist between the relaxation times and the applied forces. In an attempt for the same, we show the variations of the relaxation times normalized by τ, τ/τ (symbols), with the relative force difference, given by ΔF = Fapplied – 90 pN, in the range of 90–190 pN in Figure . The solid line is the single exponential fit by the expression 0.9 exp[−0.052 ΔF], which shows a near-exponential decreasing trend in τ/τ with ΔF. Therefore, our result is qualitatively equivalent to the observation made by Hyeon et al., from which the possible influence of an external perturbation on the intrinsic dynamics is rather clear. In this regard, the use of a constant relaxation time in ref (14) at each of the experimental stretching forces seems not to be very clear.

Figure 3

Variation of the normalized relaxation times, τ/τ (symbols), with the relative force difference, ΔF. The decreasing behavior follows a near-exponential relationship (solid line), given by 0.9 exp[−0.052 ΔF].

Having the best agreement with the experimental observations, we focus our attention to estimate the different system parameters. As mentioned previously, from the linear nature of the A vs F plot, it is possible to quantify Δxav and ΔGav for the unfolding process. In Figure , we use the extracted values of A from Table at different stretching forces and achieve the linear fit for which coefficient of determination is about 0.99. From the slope of the A vs F plot, we estimate Δxav as 0.23 nm, which perfectly reproduces the result of Kuo et al. The corresponding results from refs (14) and (16) are 0.26 and 0.24 nm, respectively. On the other hand, the average barrier height ΔGav was found to be 50.9 kJ mol–1 (24.8 kJ mol–1 from ref (14)), which also provides the perfect estimate of the experimental result, given by 51 kJ mol–1 (85.1 pN nm).[6] Here, it is to be mentioned that a nonlinear fit was used by Zheng et al. to determine the system parameters. However, the direct estimation of ΔGav from the fit was not possible unless an approximation for the harmonic frequencies in the well and barrier regions was taken into account.

The estimates of the mean time of unfolding, ⟨τunfold⟩, in the range of the applied forces are shown in Table . With increasing F from 90 to 190 pN, ⟨τunfold⟩ decreases, signifying a faster unfolding under the action of a higher force. Interestingly, the given results allow us to make a comparison with the relaxation times used for fitting at different forces from which the evidence of dynamic disorder can be inferred. It is easily noticeable that the relaxation timescales are either comparable to or slower than the mean timescales of the reaction, which satisfy the criteria for dynamic disorder to be present in the mechanical unfolding of single polyubiquitin molecules.

Table 2

Relaxation Time τ for the Best Fits of the Survival Probability Data from Ref (6) and the Mean Time of Unfolding ⟨τunfold⟩ at Different Forces Calculated from eq a

force (pN)	τ (s)	⟨τunfold⟩ (s)
90	31.0	13.94
110	4.8	2.40
130	4.0	1.65
150	1.9	0.70
170	0.43	0.15
190	0.12	0.051

Evidence of dynamics disorder in the unfolding kinetics is revealed by the comparisons between τ and ⟨τunfold⟩.

Conclusions

In summary, the microscopic origin of the nonexponential unfolding kinetics has been well accounted for by our work based on the dynamic disorder in the reaction pathway due to the conformational fluctuations of the polyubiquitin molecules. Considering the Rouse-like dynamics for the fluctuating radius of the bottleneck that characterize the intrinsic dynamics, our theoretical results provide excellent fits to the experimental data throughout the range of the stretching forces. Using the extracted parameters from the fits, we quantify the system parameters perfectly. Therefore, the present study suggests a plausible alternative approach for interpreting the results from single-molecule pulling experiments, which may also be found applicable to rationalize the effects of other perturbations such as membrane potentials in the studies of DNA escape kinetics from biological nanopores.