Decoding the mechanical fingerprints of biomolecules

Abstract

The capacity of biological macromolecules to act as exceedingly sophisticated and highly efficient cellular machines – switches, assembly factors, pumps, or motors – is realized through their conformational transitions, that is, their folding into distinct shapes and selective binding to other molecules. Conformational transitions can be induced, monitored, and manipulated by pulling individual macromolecules apart with an applied force. Pulling experiments reveal, for a given biomolecule, the relationship between applied force and molecular extension. Distinct signatures in the force–extension relationship identify a given biomolecule and thus serve as the molecule's ‘mechanical fingerprints’. But, how can these fingerprints be decoded to uncover the energy barriers crossed by the molecule in the course of its conformational transition, as well as the associated timescales? This review summarizes a powerful class of approaches to interpreting single-molecule force spectroscopy measurements – namely, analytically tractable approaches. On the fundamental side, analytical theories have the power to reveal the unifying principles underneath the bewildering diversity of biomolecules and their behaviors. On the practical side, analytical expressions that result from these theories are particularly well suited for a direct fit to experimental data, yielding the important parameters that govern biological processes at the molecular level.

1. Introduction: ‘Fingerprinting’ biological macromolecules

Conformational transitions in biological macromolecules – such as folding, unfolding, binding to and dissociating from other molecules – serve as the mechanism that brings macromolecules into a working shape, enabling them to act as switches, assembly factors, pumps, or force- and displacement-generating motors (Alberts et al. 2014). Let us consider the following conceptual approach to probe these conformational transitions: we hold one end of the macromolecule fixed and pull on the other end with an applied force, for example, through a Hookean spring. As the macromolecule undergoes conformational rearrangements, e.g. unfolds, refolds, binds, unbinds or twists, we use the instantaneous extension of the spring as a reporter of these conformational dynamics. This simple concept is played out in a class of rather sophisticated biophysical experiments – single-molecule pulling experiments – in which the role of the abstract ‘pulling spring’ is fulfilled by the cantilever of an atomic force microscope, a focused laser beam in optical tweezers, or a magnetic field in magnetic tweezers. The force can also be exerted with an electric field that ‘threads’ the macromolecule through a nanopore, or with a hydrodynamic flow that pulls the macromolecule apart by acting on a bead to which the macromolecule is tethered (De Souza, 2012; Dudko et al. 2010; Greenleaf et al. 2007). The biomolecule is often tethered to the pulling device through molecular linkers, typically double-stranded DNA or polypeptides; however, the behavior of the linkers under force is usually not simple enough (i.e. it is non-linear) for them to be idealized as Hookean springs. These pulling experiments have made it possible to monitor and manipulate conformational transitions in individual DNA, RNA, and proteins, as well as assemblies of different macromolecules.

In single-molecule force spectroscopy experiments, the applied force is increased (in the force ramp protocol) or decreased (in the force-relaxation protocol) at a given loading rate. The response of a biomolecule to force can be captured in the form of a relationship between the instantaneous force, experienced by the biomolecule, and its end-to-end extension. Interestingly, force–extension curves of nucleic acids, proteins, and molecular complexes often exhibit abrupt drops and increases, or ‘rips’ (Fig. 1). These rips are signs of conformational transitions. An abrupt drop in the force is associated with an unfolding or unbinding event: the rupture of a structural element within the molecule and the subsequent release of the unstructured segment leads to a decrease in the extension of the pulling spring and hence in the measured force. Likewise, an abrupt increase in the force is associated with a folding or binding event: a previously unstructured segment within the molecule forms bonds with other parts of the molecule, decreasing the molecule's extension and thereby increasing the extension of the pulling spring. Specific patterns of the rips in the force–extension curves identify a given biomolecule and thus serve as its ‘mechanical fingerprints’ (Marszalek et al. 2001). The possibility of revealing the fine details of the conformational transitions in biomolecules by taking their mechanical fingerprints opens an unprecedented opportunity to gain new insights into biomolecular interactions. However, this opportunity comes with a challenge: how can we interpret the mechanical fingerprints in quantitative terms to reveal the kinetic barriers and timescales of the associated biological process?

Fig. 1. Decoding the mechanical fingerprints of a complex macromolecule. Upper: A free-energy profile with a cascade of barriers separating the folded and unfolded states of the macromolecule. The macromolecule (green) is depicted in three of its many conformational states. Red arrows indicate the barriers and timescales that are sought to be reconstructed. Lower: The force–extension curve featuring ‘rips’, or ‘mechanical fingerprints’, from a force-ramp experiment on a single macromolecule. Such an experiment can be performed, for example, using optical tweezers as shown schematically above the force–extension curve. The force–extension curve is generated via Brownian dynamics simulations (Zhang & Dudko, 2013). Each rip, that is, abrupt drop or increase, in the curve signifies a conformational transition over the corresponding barrier. From multiple repeats of this experiment, histograms (shown in insets) of the transition forces for different types of transitions are collected. The histograms of transition forces contain a wealth of information about the conformational dynamics and constitute the key output of force spectroscopy experiments. How can we decode these mechanical fingerprints so that they reveal the kinetic barriers and timescales This review summarizes recent analytical theories that address this question.

Conformational transitions in most biological macromolecules are governed by complex energy landscapes (Fig. 1) (Anthony et al. 2012; Kaiser et al. 2011; Neupane et al. 2011; Stigler et al. 2011), often with multiple barriers and even multiple reaction pathways. This complexity is reflected in the elaborate features of the resulting mechanical fingerprints. The problem of interpreting the mechanical fingerprints is further complicated by the fact that the force that induces and probes these conformational transitions is time dependent, which makes the rates of the transitions over individual barriers also time dependent. One of the central challenges in single-molecule biophysics is to develop theoretical foundations that capture this complex behavior in sufficient generality but, at the same time, are sufficiently simple to be analytically tractable. Indeed, while computational approaches have proven to be a valuable tool for understanding and predicting the mechanisms of biomolecular transitions under mechanical force (Chen & Garcia, 2012; Dudko et al. 2011; Izrailev et al. 1997; Heyon & Thirumalai, 2011; Konda et al. 2014; Sulkowska et al. 2009), there is a special value to analytical approaches: analytical expressions for the experimentally measurable quantities (i) have the generality that unifies the behavior of a broad range of systems and experimental techniques; and (ii) allow one to extract mechanistic insights through a direct fit to experimental data.

This review summarizes recent analytical theories that address the challenge of interpreting single-molecule force spectroscopy measurements. It will be shown that the difficulty of deriving an analytical form of the key experimental output of force spectroscopy – the distribution of transition forces – for complex biomolecules with multiple kinetic barriers can be bypassed by applying a general, model-free transformation to these forces (Section 2). At the same time, despite the abovementioned complication arising from the time dependence of the transition rates (which is the consequence of the time-dependent applied force), analytical expressions for the distribution of transition forces can nevertheless be derived for several important cases (Sections 3.1 and 3.2). The closed-form expressions, which result from the analytical theories, can be used as a fitting tool, yielding the important timescales and energy barriers that govern the underlying conformational transitions.

2. Decoding the complex mechanical fingerprints

2.1 Transforming the fingerprints into a decodable form

Complex biomolecules, i.e. biomolecules that unfold or dissociate via intermediate states rather than in a simple two-state fashion, are characterized by equally complex mechanical fingerprints. Individual rips in the fingerprints reflect the molecule crossing the corresponding barriers on its multi-barrier energy landscape. However, extracting the parameters of the barriers, as well as the timescales associated with these barriers, from the fingerprints of a multi-barrier system is not a trivial task. The challenge originates from the fact that the transition force at a given rip (say, rip number i) is not determined by the barrier number i alone, but also by all the preceding barriers that the molecule had to cross on the way to the ith barrier. This dependence of measured forces on the history of the multi-barrier transition arises from the time-dependent nature of the applied force, or the fact that, while the molecule wanders among the intermediate barriers on the way to the ith barrier, the external force keeps changing. Consider, for example, the ramping protocol, in which the applied force is increased with time. If the molecule is significantly delayed in the intermediate states because of a high barrier and/or simply because of ‘bad luck’ (i.e. lack of favorable fluctuations), by the time the molecule arrives to the ith barrier, the force has already been ramped up significantly. As the result, the measured force upon the escape over the ith barrier will be higher than it would have been in the absence of preceding barriers.

The consequence of this interplay between the changing force that perturbs the energy landscape and the intrinsic timescales set by the barriers on this landscape is that the measured transition forces for individual barriers contain entangled contributions from other barriers. Thus, in order to interpret the measured transition forces quantitatively, it would be helpful to find a way to disentangle the contributions from individual barriers in the measured forces. A simple solution to this problem can be found by exploiting the quasi-adiabatic approximation (detailed in Section 3.1.1) of the force-ramp experiments (Dudko et al. 2006, 2008; Evans et al. 2009; Raible et al. 2004). The solution has the form of a mathematical transformation (Zhang & Dudko, 2013):

(1a) $$\matrix{ {k_{ij} \left( F \right) = \left \vert {\dot F_i (F)} \right \vert \displaystyle{{P_{ij} (F)} \over {{\rm {\cal N}}_i (F)}},} & {\forall \; i \to j} \cr}. $$

Here, P _ij (F) is the force histogram for the transition from state i to state j, ${\rm {\cal N}}_i (F)$ is the population and $\vert {\dot F_i (F)} \vert$ is the absolute value of the force-loading rate in state i at force F (Fig. 2). With the characteristics of the mechanical fingerprints as the input on the right-hand side, this transformation yields, as an output on the left-hand side, the force-dependent transition rates for all the transitions that can possibly occur in the system. The resulting full set of force-dependent rates constitutes ‘the rate map’ of the system (Fig. 3) (Zhang & Dudko, 2013).

Fig. 2. Applying the transformation in Eqs. (1a , b ) to the mechanical fingerprints: step-by-step illustrations. (a) Free-energy profile of a system with two sequential barriers. Indicated are the intrinsic transition rates (each rate is the inverse of the corresponding timescale) and the barrier heights and locations, which are the parameters sought to be reconstructed. (b) Two representative force-extension curves from a stretching and relaxation cycle. Indicated are the transition forces for different types of transitions. (c) Force histograms for the different types of transitions. The histograms are collected from the force–extension curves from stretching (left column) and relaxation (right column) cycles at the nominal loading rates indicated (in piconewtons per second) next to each histogram. Dividing the raw number of counts in the bin by the corresponding bin width yields P _ij (F) in Eq. (1a ) for a transition from state i to state j. (d) Finding ${\rm {\cal N}}_i (F)$ , the number of molecules (trajectories) in the ith state at force F. (e) Determining ${\dot F}_i (F)$ , the force-loading rate in different states, from the slopes of the force trajectory in these states.

Fig. 3. The transformation at work: from mechanical fingerprints to the rate map to activation barriers. Left: A system with two sequential barriers (top) and the corresponding rate map (bottom) obtained via Eqs. (1a , b ). The four branches on the rate map correspond to the four possible transitions on the two-barrier potential. The colors of the arrows on the potential graph correspond to different types of transitions and match the colors of the corresponding branches on the rate map. Middle: A system with a cascade of barriers and eight possible transitions. Right: A system with two competing pathways, one of which features an intermediate. To test the robustness of the transformation, the effect of an anharmonic linker was incorporated in the simulations of the systems in the left and middle panels. Error bars are calculated with Eq. (2). Lines are the fit to Eq. (3). As an example, Table 1 lists the heights and locations of the barriers and the associated rates extracted from the fit for the system in the left panel.

Let us summarize the key aspects of the transformation in Eq. (1a ), both of fundamental and of practical nature:

1. (1) The transformation is model-free: it makes no model assumptions regarding the shape of the free-energy barriers or the functional form of the force-dependent rates. The model-free nature of the transformation stems from its ability to implement, with no model assumptions, the statistical–mechanical ‘flux-over-population’ formulation of the transition rate (van Kampen, 2007; Zhou, 2010).

2. (2) The transformation is quite general, as it is applicable to transitions that involve arbitrarily large number of barriers and reaction pathways, as long as the individual transitions can be resolved as rips in the force–extension curves.

3. (3) The transformation is easy to implement in practice, as it only involves the quantities that can be directly accessed in force spectroscopy experiments. Figure 2 provides a step-by-step illustration of the implementation. To highlight the ease of the implementation, we can rewrite the transformation in the form of a word equation:

   (1b) $${\rm rate}_{{\rm state}\; i \to {\rm state}\; j} \left( {{\rm force}\,F} \right) = \displaystyle{{\left \vert {{\rm loading}\;{\rm rate}\;{\rm in}\;{\rm state}\;\; i\,\; {\rm at}\,{\rm force}\,F} \right \vert} \over {{\rm number} \,\,{\rm of}\;{\rm trajectories\;} {\rm in\;} {\rm state}\,\; i\,\; {\rm at}\,\; {\rm force}\,\; F}} \times \displaystyle{{{\rm number}\,\,{\rm of}\;{\rm counts}\,{\rm in}\,{\rm bin}\,F} \over {{\rm width\;} {\rm of}\,{\rm bin}\,\; F}}.$$

4. (4) The transformation is applicable to the fingerprints obtained both in a ramping (stretching) and in a relaxation force protocol.

5. (5) The transformation naturally accounts for the reversible behavior, or back-and-forth ‘hopping’ between states. Hopping is often ignored in the analytical treatments of force–extension relationships because it tends to complicate the equations for the experimental observables. In contrast, hopping behavior in the mechanical fingerprints works to one's advantage when using the transformation in Eqs. (1a , b ): hopping is equivalent to repeatedly sampling a given transition, which provides additional data points for the transition forces that can be used as the input in the transformation. For example, hopping yields two values of F _NI in a single force–extension curve in Fig. 2b .

6. (6) The transformation easily accounts for the dependence of the force-loading rate on the force itself. Such dependence, which typically complicates an analytical treatment, can result from the nonlinear force–extension relationship of the linker that tethers the molecule under study to the pulling device. In the transformation, however, the force-loading rate ${\dot F}(F)$ appears as an independent factor. This factor can be readily determined from the slope of the force trajectory in the state of interest, as illustrated in Fig. 2e .

7. (7) Even though force-dependent rates can, in principle, be measured in a force clamp experiment, the immediate outcome of the transformation – the rate map – yields the rates of the entire spectrum of transitions over a broad range of forces (Fig. 3). The covered force range is typically beyond what can be accessed directly in a force clamp. The coverage of the broad range of forces is achieved through the use of data from force ramp experiments that are, in turn, performed in a broad range of loading rates, with each loading rate effectively ‘probing’ the corresponding range of forces.

In practice, the number of measurements is, of course, finite, which contributes to the statistical errors in the force histograms P _ij (F) and in the population ${\rm {\cal N}}_i (F)$ in a given state. How do the errors in the histograms and populations propagate into the rates that are obtained by transforming the histograms with Eqs. (1a , b )? The expression for the standard deviation of the rates on the rate map can be found by computing the standard deviations of each component in Eqs. (1a , b ) and subsequently applying the standard error propagation rule. The resulting standard deviation in the logarithmic rate k _ij (F _m ) for transitions from state i to state j at the force F _m is (Zhang & Dudko, 2013)

(2) $$\sigma _{{\rm ln}k_{ij} (F_m )} \approx \left[ {\displaystyle{1 \over {P_{ij} (F_m )\Delta F_m}} + \displaystyle{1 \over {{\rm {\cal N}}_i (F_m )}}} \right]^{1/2}, $$

where P _ij (F _m )ΔF _m is the number of counts in the mth bin with the bin width ΔF _m , and ${\rm {\cal N}}_i (F_m )$ is the number of trajectories found in the state i at the force F _m (Fig. 2d ). The vertical error bars seen in the rate maps in Fig. 3 have been determined with Eq. (2).

2.2 Extracting rates and rate-limiting barriers

The effect of the transformation in Eqs. (1a , b ) is to decompose the response of a multi-barrier system to force into contributions from the individual barriers of the system. Specifically, Eqs. (1a , b ) convert complex mechanical fingerprints into a form (the rate map) that readily reveals the barriers and timescales of the system. Indeed, each branch on the rate map (Fig. 3) represents the force-dependent transition rate k _ij (F) for a particular barrier, namely, the barrier separating states i and j. Consequently, each branch on the rate map can be analyzed as a single-barrier problem. To carry out this analysis, the following expression for the force-dependent rate can be used to fit the individual branches on the rate map (Zhang & Dudko, 2013):

(3) $$\eqalign{k(F) = & k_0 \left[ {1 + \displaystyle{{\nu \kappa x^{{\ddag} 2}} \over {2 \Delta G^{{\ddag}}}} \mp \displaystyle{{\nu Fx^{{\ddag}}} \over {{\Delta} G^{{\ddag}}}} \left( {1 + \displaystyle{{(1 - \nu )\kappa x^{{\ddag} 2}} \over {2 \Delta G^{{\ddag}}}}}\right) }\right]^{1/\nu - 1} \cr & \times \exp \left\{ {\beta {\Delta} G^{{\ddag}} \left[ {1 - \left[ {1 + \displaystyle{{\nu \kappa x^{{\ddag} 2}} \over {2 \Delta G^{{\ddag}}}} \mp \displaystyle{{\nu Fx^{{\ddag}}} \over {{\Delta} G^{{\ddag}}}} \left( {1 + \displaystyle{{(1 - \nu )\kappa x^{{\ddag} 2}} \over {2 \Delta G^{{\ddag}}}}} \right) }\right]^{1/\nu}} \right]} \right\},}$$

where β = (k _B T)⁻¹ is the inverse of the product of the Boltzmann constant and the absolute temperature, ∓ applies to forward/backward transitions in both the stretching and relaxation cycles, and κ is the spring constant of the pulling device. Equation (3) has been derived by solving the Kramers problem (Kramers, 1940) on a force-dependent one-dimensional (1D) free-energy profile. The unifying nature of the solution in Eq. (3) is evident from the scaling factor ν, which specifies different models of the free-energy profile: a smooth barrier that is described by a linear–cubic polynomial corresponds to ν = 2/3, while a cusp-like barrier with a harmonic well corresponds to ν = 1/2. A simple least-squares fit of each branch on the rate map to Eq. (3) yields the intrinsic parameters for the height ΔG ^‡ and location x ^‡ of the corresponding activation barrier and the transition rate k ₀ associated with this barrier. The fit (with the value of ν fixed at 1/2) is shown as lines in Fig. 3. Table 1 lists the parameters extracted from the fit for the system with two sequential barriers.

Table 1. Intrinsic rates (in s⁻¹), barrier heights (in k_BT) and barrier distances (in nm) extracted from the fit of the rate map to Eq. (3). The underlying free-energy profile with two sequential barriers is shown in Fig. 3 (top left) along with the corresponding rate map (bottom left, fit is shown as lines).

Equation (3) is accurate beyond the approximation (F ≈ κVt) typically used for the applied force F = κ(Vt − x). This higher degree of accuracy makes Eq. (3) applicable not only to transitions that originate in ‘stiff’ states (i.e. states characterized by narrow minima on the free-energy landscape) but also transitions that originate in ‘soft’ states, such as the unfolded state, for which the abovementioned approximation may not be sufficiently accurate. Detailed analysis shows that Eq. (3) is valid for the biomolecular stiffness values as low as ~5κ, where κ is, as before, the pulling spring constant. Given the typical spring constant of an AFM ~5 pN nm⁻¹ and an optical trap ~0.3 pN nm⁻¹, and with an additional softening effect of the linker, Eq. (3) applies to most experimental situations in single-molecule force spectroscopy of biomolecules (Zhang & Dudko, 2013).

The force-dependent rates obtained via Eqs. (1a , b ) from the force ramp should, of course, agree with the rates measured directly under constant force (i.e. in the force clamp). However, when comparing the rates obtained in these two pulling modes, it has to be taken into account that the force clamp and the force ramp have different effects on the molecular potential. The bias imposed by the force clamp on the molecular potential is described by the term −Fx and is thus linear in the extension x, while the corresponding bias in the force ramp is described by 1/2κ(x − Vt)² and is thus non-linear in x. To account for this difference, the data point for the rate k(F) measured at constant force F should be transferred to the rate map as k(F′), where $F^{\prime} = (F \pm \kappa x^{\ddag} /2)/(1 + (1 - \nu )\kappa x^{\ddag 2} /(2 \Delta G^{\ddag} ))$ . Here ‘+’ applies to forward transitions (i.e. those toward larger values of x, such as the unfolding) and ‘–’ to backward transitions (such as the refolding).

3. Special but important cases

3.1 Irreversible transitions over a single barrier

The case when conformational transitions involve a single rate-limiting barrier is worth special consideration. In this case, the problem of force-induced molecular transitions can be solved analytically for the major outputs of both the force clamp and force ramp experiments. The derivation can be carried out within the 1D picture where the molecular extension x serves as a reaction coordinate. This case covers two-state systems, in which the folded and unfolded states are separated by a single barrier. In multi-state systems, this case applies to transitions over the first barrier.

3.1.1 Unfolding and unbinding

Force-induced conformational transitions can be treated as diffusive escape over a free-energy barrier in a bias field (Evans & Ritchie, 1997; Izrailev et al. 1997; Shapiro & Qian, 1997). This treatment is an extension of the view of conformational transitions as Brownian motion over a potential barrier, which is the idea underlying Kramers theory of chemical reactions (Kramers, 1940; van Kampen, 2007; Zwanzig, 2001). Posing and solving Kramers problem – that of the rate at which a Brownian particle escapes from a potential well – on a class of model potentials that feature a single barrier, the following analytical solution has been derived for the rate k(F) at constant force F (Dudko et al. 2006):

(4) $$k(F) = {K_0} \left({1 - \displaystyle{{vF{x^{{\ddag}}}} \over {\Delta {G^{{\ddag}}}}}}\right)^{1/v - 1} \exp \left\{{ \beta \Delta {G^{{\ddag}}}\left[{1 - \left({1 - \displaystyle{{vF{x^{{\ddag}}}} \over {\Delta {G^{{\ddag}}}}}}\right)}\right]}\right\}.$$

The predictive value of Eq. (4) lies in its capacity to relate a measurable quantity – the rate k(F) at force F – to the intrinsic (zero-force) characteristics of the system: the height ΔG ^‡ and location x ^‡ of the activation barrier and the unfolding (or unbinding) rate k ₀. These intrinsic parameters can be readily extracted with Eq. (4) through a least-squares fit of the rates measured in a force clamp. As in the case of Eq. (3), the unifying nature of Eq. (4) is evident from the scaling factor ν, which specifies different models of the 1D free-energy profile. As a test of relative insensitivity of the extracted parameters to the details of the shape of the barrier, the force-dependent rates can be fit to Eq. (4) with the values of ν fixed at 1/2 and 2/3, and the resulting values of ΔG ^‡, x ^‡ and k ₀ compared. A notable property of the solution in Eq. (4) (as well as of its generalized version, Eq. (3)) at physically meaningful values of ν (1/2 and 2/3) is a characteristic concave down curvature in the rate at intermediate-to-high forces on a semi-log scale (Fig. 4a ). This curvature reflects the intrinsic property of the potential barrier and well to move toward each other with an increase in the external force. In contrast, this property is ignored in the Bell's postulate for the force-dependent rate (Bell, 1978). Bell's expression is recovered in Eq. (4) with ν = 1.

Fig. 4. Unfolding/unbinding transitions over a single barrier under applied force, from simulations (symbols and histograms) and analytical theory (black lines). (a) The force-dependent rate from the force-clamp pulling mode. Note the characteristic nonlinearity in the rate on a semi-log plot. Line is Eq. (4) with ν = 2/3. (b) A representative force–extension curve from the force–ramp pulling mode. The ‘rip’ signals an unfolding/unbinding event, the corresponding unfolding/unbinding force value is circled. From multiple repeats of this experiment, histograms of the unfolding forces (see inset) are collected. (c) Distributions of unfolding/unbinding forces from the force–ramp mode at three values of the loading rate. Note the characteristic negative skew in the distributions. Lines are Eq. (5) with ν = 2/3. (d) Average force as a function of the loading rate, or ‘the force spectrum’. Note the characteristic upward curvature. Line is Eq. (6) with ν = 2/3.

The functional form of the rate versus force, established in Eq. (4), can be used to derive an analytical expression for the distribution of rupture forces, the key output of a force ramp experiment (Fig. 4b, c ). Under the assumption that the external force is ramped up with time slowly enough such that a barrier of at least a few k _B T is present at all rupture events (the quasi-adiabatic approximation), which reflects a typical experimental situation, the following expression has been derived (Dudko et al. 2006):

(5) $$p(F) = \displaystyle{{k(F)} \over {\kappa V}}{\rm exp}\left( {\displaystyle{{k_0} \over {\beta x^{\ddag} \kappa V}}} \right){\rm exp}\left[ { - \displaystyle{{k(F)} \over {\beta x^{\ddag} \kappa V}}\left( {1 - \displaystyle{{\nu Fx^{\ddag}} \over {{\Delta} G^{\ddag}}}} \right)^{1 - 1/\nu}} \right].$$

An analytical form of the solution in Eq. (5) makes it convenient for a fit of transition-force histograms obtained in force-ramp experiments. A fit (Fig. 4c , lines) readily yields the values of the height ΔG ^‡ and location x ^‡ of the barrier and the rate k ₀ in the absence of a force. The accuracy of the extracted parameters is enhanced in a global fit. A global fit is fit of multiple histograms at once with a single set of fitting parameters {k ₀, ΔG ^‡, x ^‡}, where the histograms are measured over a broad range of values of the pulling speed.

The force distribution in Eq. (5) has a characteristic negative skew (Fig. 4c ). The longer low-force tail of the distributions reflects the rare unfolding events that occur at short times – and hence at low forces – when the barrier is still high. Increasing the loading rate shifts the distribution to higher forces (Fig. 4c ). This trend in the shift is quite intuitive: when the applied force increases rapidly (which happens at high loading rates), there is less time for thermal fluctuations to kick the molecule over a high barrier, causing transitions to occur later in time – which means, at higher forces (hence the shift to higher forces) – when the barrier is lower.

The probability distribution of rupture forces in Eq. (5) has been used to derive an analytical expression for the mean rupture force as a function of the loading rate κV, or ‘the force spectrum’ (Dudko et al. 2006):

(6) $$\left\langle F \right\rangle \cong \displaystyle{{{\Delta} G^{\ddag}} \over {\nu x^{\ddag}}} \left[ {1 - \left( {\displaystyle{1 \over {\beta {\Delta} G^{\ddag}}} {\rm ln}\displaystyle{{k_0 \,{\rm e}^{\beta {\Delta} G^{\ddag} + \gamma}} \over {\beta x^{\ddag} \kappa V}}} \right)^\nu} \right],$$

where γ = 0.577… is the Euler–Mascheroni constant. The expression in Eq. (6) predicts an increase in the mean rupture force with the loading rate as well as a convex upward curvature in the force spectrum (F versus ln(κV)) at intermediate-to-high loading rates on the semi-log scale at 1/2 > ν > 2/3 (Fig. 4d ). This curvature is a manifestation of the force-induced motion of the transition barrier toward the potential well on the energy landscape, as mentioned above following Eq. (4).

3.1.2 Folding and binding

Equation (5) describes the distribution of unfolding or unbinding forces that are measured in force–ramp experiments, in which the force is increased with time until an unfolding event is detected. Likewise, in force-relaxation experiments, the force is decreased with time, which allows measurements of the reverse processes, the folding or binding (Fig. 5b ). The expression for the folding or binding force distribution, which is analogous to Eq. (5) for unfolding or unbinding, has been derived as (Pierse & Dudko, 2013)

(7) $$p(F) = \displaystyle{{k_ \leftarrow (F)} \over {\left \vert {\dot F(F)} \right \vert}} {\rm exp}\left[ { - \displaystyle{{k_ \leftarrow (F)} \over {\beta \left \vert {\dot F(F)} \right \vert x^{\ddag}}} \left( {1 + \displaystyle{{\kappa _S} \over {\kappa _U (F)}}} \right)^{\nu /(\nu - 1)} \left( {1 + \displaystyle{{\nu Fx^{\ddag}} \over {{\Delta} G^{\ddag}}}} \right)^{1 - 1/\nu}} \right].$$

Fig. 5. Folding/binding transition over a single barrier against applied force, from simulations (symbols and histograms) and analytical theory (lines). (a) The force-dependent folding rate from the force–clamp pulling mode. Note the characteristic nonlinearity on a semi-log plot. Line is Eq. (4) (with x ^‡ < 0 for the folding or binding process) evaluated at ν = 2/3. (b) A representative force–extension curve from a force-relaxation experiment with the value of the folding force circled. From multiple repeats of this experiment, histograms of the folding forces are collected (inset). (c) Distributions of folding/binding forces from the force-relaxation regime at three values of the relaxation rate. Note the characteristic positive skew in the distributions. Lines are Eq. (7) with ν = 2/3.

The force-dependent folding rate that enters this expression is given by

$$k_\leftarrow (F) = k_0 \left( {1 + \displaystyle{{\kappa _S} \over {\kappa _U (F)}}} \right)^{1/\nu - 1/2} \left( {1 + \displaystyle{{\nu Fx^{\ddag}} \over {{\Delta} G^{\ddag}}}} \right)^{1/\nu - 1} {\rm exp}\left\{ {\beta {\Delta} G^{\ddag} \left[ {1 - \left( {1 + \displaystyle{{\kappa _S} \over {\kappa _U (F)}}} \right)^{2\nu /(1 - \nu ) - 1} \left( {1 + \displaystyle{{\nu Fx^{\ddag}} \over {{\Delta} G^{\ddag}}}} \right)^{1/\nu}} \right]} \right\},$$

where the stiffness of the unfolded biomolecule at force F and the force-relaxation rate are, respectively,

$$\kappa_U (F) = \displaystyle{\nu \over {(1 - \nu)^2}} \displaystyle{{{\Delta} G^{\ddag}} \over {x^{{\ddag}2}}} \left({1 + \displaystyle{{\nu Fx^{\ddag}} \over {{\Delta} G^{\ddag}}}} \right)^{2 - 1/\nu}, \quad {\dot F}(F) = \displaystyle{{-\kappa_S V} \over {1 + \kappa_S /\kappa_U (F)}}.$$

The solution in Eq. (7) reveals properties of the distribution of folding forces that distinguish this distribution from its unfolding counterpart, Eq. (5). The distinguishing properties of these two distributions can also be seen from the comparison between Figs 4c and 5c . First, the distribution of folding/binding forces has a positive skew (Fig. 5c ). The longer tail at higher forces reflects the rare folding events that occur at short times – and hence at high forces – when the barrier is still high. Second, increasing the absolute value of the relaxation rate $\vert {\dot F} \vert$ shifts the distribution to lower forces (Fig. 5c ). This trend can be understood with simple physical arguments: when the applied force decreases rapidly (which is the case at high relaxation speeds), there is less time for thermal fluctuations to kick the molecule over a high barrier, causing transitions to occur later – and hence at lower forces – when the barrier is lower. Third, the stiffness κ _s of the pulling device and speed V enter Eq. (7) separately, providing independent routes to control the range of folding forces and thus enhance the robustness of a fit of experimental distributions to Eq. (7). A global fit (the meaning of a global fit in the present context is discussed in Section 3.1.1) yields the intrinsic height ΔG ^‡ and location x ^‡ of the barrier to folding as well as the on-rate k ₀. The application of the Eq. (7) to binding experiments on a ligand and receptor connected by a tether involves an additional step associated with the decoupling of the effect of the tether and allows the reconstruction of the parameters of ligand–receptor binding (Pierse & Dudko, 2013).

3.2 Barely so, but still analytically tractable: two sequential barriers

Section 3.1 described a direct analytical approach to the analysis and interpretation of the data from pulling experiments. This approach takes advantage of the availability of analytical expressions that are suitable for fitting the information-rich characteristics of the conformational dynamics – the distribution of transition forces p(F). The fit, in turn, yields the activation barrier and timescale of the conformational transition. Balancing the generality with simplicity, this direct approach is particularly efficient when we deal with relatively simple systems for which an analytical expression for p(F) can be derived. However, even modest complications of the system, such as the presence of more than one energy barrier, rapidly complicate the corresponding formulae and require additional parameters, increasing the risk of overfitting. Mathematically, the problem is challenging due to the time-varying perturbation caused by the external force F(t). The force introduces the time-dependence in the rate coefficients k _ij (t), which, on the one hand, enriches the behavior of the system but, on the other hand, makes it challenging to solve the corresponding set of rate equations analytically.

Let us consider a case of relatively complex dynamics, for which the direct analytical approach nevertheless leads to expressions simple enough to enable a meaningful fit of the force–ramp data. This is the case of a system in which the native and unfolded states are separated by one intermediate state (Fig. 6). For a broad class of topologies of energy landscapes with an intermediate state, the expression for the distribution of forces at which the system transits into the final (unfolded) state has been derived as (Garai et al. 2014)

(8) $$p(F) = \displaystyle{1 \over {\kappa V}}\left\{ {k_{IU} \left( F \right){\rm e}^{ - g_{IU} (F)} \left[ {1 - {\rm e}^{ - g_{NI} (F^{\prime})}} \right] + k_{NI} (F){\rm e}^{ - g_{NI} (F^{\prime} )} \left[ {1 - {\rm e}^{ - g_{IU} (F)} \left( {1 + g_{IU} (F)} \right)} \right]\left[ {1 + g_{IU}^{ - 1} (F)} \right]} \right\}.$$

Fig. 6. Signatures of sequential barriers on the free-energy profile. (a) Free-energy profile featuring two sequential barriers: a stiff (narrow) and low barrier is followed by a soft (broad) and high barrier. (b) Gray line: force-dependent effective rate from a force-clamp pulling mode. Also shown are the individual rates for the forward transitions over each of the two barriers (orange and red) and the backward transition over the first barrier (green). (c) Representative force–extension curves from a force-ramp pulling mode at low, intermediate and high loading rates. Arrows indicate the transition events into the unfolded state; the corresponding transition force values are collected in histograms that are shown in (d). The ‘trapping’ effects of the native and intermediate states at different timescales are illustrated in the insets with the populations of particles accumulated in the corresponding ‘trap’. (d) The distributions of unfolding forces at three values of the loading rate. A transient bimodality in the distribution can be seen at intermediate loading rates. Lines are Eq. (8) with ν = 2/3. (e) A visual representation of the origin of bimodality in the force distribution. Shown on the left are the snapshots of the free-energy profile at different instants during a force-ramp experiment. The distribution, P(F), is equivalent to the flux into the final state times the inverse of the force-loading rate ${\dot F}$ . The flux, in turn, is given by the transition rate (governed by the barrier) times the population. The width of the arrows over the barriers reflects the relative value of the flux at different instants. The number of particles in each potential well reflects the relative population in this well at the given instant. Shown on the right are the distributions of forces at the transition into the final state with early transition events contributing to the first peak and late events contributing to the second peak (highlighted in red).

In this expression, the force-dependent rates for transitions between various states (denoted as i and j) are given by

$$k_{ij} (F) = k_{ij}^0 \left( {1 - \displaystyle{{\nu Fx_{ij}^{\ddag}} \over {{\Delta} G_{ij}^{\ddag}}}} \right)^{1/\nu - 1} {\rm exp}\left\{ {\beta \Delta G_{ij}^{\ddag} \left[ {1 - \left( {1 - \displaystyle{{\nu Fx_{ij}^{\ddag}} \over {\Delta G_{ij}^{\ddag}}}} \right)^{1/\nu}} \right]} \right\}$$

and the functions g _ij (F) are

$$g_{ij} (F) = \displaystyle{{k_{ij}^0} \over {\beta \kappa Vx_{ij}^{\ddag}}} \left\{ {{\rm exp}\left[ {\beta \Delta G_{ij}^{\ddag} \left[ {1 - \left( {1 - \displaystyle{{\nu Fx_{ij}^{\ddag}} \over {\Delta G_{ij}^{\ddag}}}} \right)^{1/\nu}} \right]} \right] - 1} \right\}.$$

The above expressions account for a generally non-negligible drop in the force from F′ to F occurring in the course of the transition, upon escape from N to I, due to partial unfolding. These expressions have been derived under the assumption that backward transitions are negligible, which is justified beyond the range of very small forces.

An unusual feature in the distribution of transition forces – namely, bimodality (Fig. 6d ) – is predicted by Eq. (8) if the landscape has a relatively low and stiff barrier that is followed by a high and soft barrier (Fig. 6a ). The bimodality originates from the ‘trapping’ effect of the native and intermediate states at different stages of a force ramp experiment (Fig. 6e ). One of these ‘traps’, the intermediate state, releases its accumulated population at low forces when the second barrier is lowered sufficiently to allow fast escape. Another ‘trap’, the native state, releases the remaining population at high forces when the first, stiff, barrier has been lowered sufficiently.

The expression in Eq. (8) is suitable for a global fit (global fits are discussed in Section 3.1.1) of transition force distributions in a broad range of pulling speeds. The fit yields the locations and heights of each of the two sequential barriers, the location of the intermediate state, and the rate constants for all the transitions on the energy landscape with an intermediate (Garai et al. 2014).

4. When a single coordinate does not suffice

The conformational dynamics of biological macromolecules are inherently multidimensional, as they involve a multitude of degrees of freedom of the constituent atoms. Indeed, a small protein consisting of a few thousands of atoms is characterized by at least the order of 10⁴ degrees of freedom, to which the degrees of freedom of its surrounding environment should be added. Despite the hyper-dimensionality of the problem, the analytical expressions presented in Sections 2.2 and 3 have been derived within a 1D picture. In this picture, the conformational dynamics were captured by a single coordinate – the molecular extension x. What makes such a drastic simplification possible? The reduction of multidimensional dynamics to 1D dynamics along x is justified when x is a slow variable: changes in x are accompanied by rapid adjustments of all the other, fast, degrees of freedom. The fast degrees of freedom can then be eliminated from the explicit description. The result of this ‘adiabatic elimination’ is 1D dynamics along x, in which the effect of the fast degrees of freedom is treated as part of the thermal bath.

But, what if the system possesses a slow degree of freedom – let us call it Q – which is different from x? In this case, the elimination of Q from the explicit consideration is no longer justified. A minimal free-energy landscape in the space of two coordinates – x and Q – is shown in Fig. 7 (left). From the figure, it is apparent that the system faces a barrier in the direction of Q while equilibrating rapidly in a harmonic potential along the direction of x. Mathematically, this two-dimensional (2D) free-energy landscape is described by Dudko et al. (2008) and Suzuki and Dudko (2010) as

$$G(x,Q) = G_0 (Q) + \displaystyle{1 \over 2}K\left( Q \right)\left[ {x - x_0 (Q)} \right]^2 - Fx.$$

Fig. 7. Effect of the multidimensionality of the energy landscape on the response of a biomolecule to force. Left: A minimalist 2D model of the energy landscape in the space of two coordinates: the pulling coordinate x and a slow coordinate Q. Right: In contrast to the traditional 1D description, the 2D model gives rise to a rich spectrum of scenarios in the response of a biomolecule to force. One of these scenarios is a ‘rollover’ in the lifetime τ(F). A rollover is realized in this model when the pathway that is aligned unfavorably at low force (as seen in F = 0 snapshot, where the transition state (■) is seen to have a smaller value of x than the native state (•)) becomes aligned favorably at high force (as seen in F = 110 pN snapshot, where the transition state has moved to larger values of x, i.e. in the direction of the force application). Cartoons of the macromolecule at each snapshot illustrate how the molecular extension at the transition state relative to that of the folded state changes as the force is increased. The rollover is the consequence of the basic property of the transition state to move with respect to the folded state as the force is increased, resulting in a deformation of the reaction pathway (solid line connecting • and ■). Inset in the lifetime versus force plot shows how the barrier height increases with force at low forces and subsequently decreases at high forces as the result of the gradual alignment of the reaction pathway along the direction of the force.

In the above expression, the last term reflects the effect of force acting along the pulling coordinate x.

The force-dependent rate on this 2D landscape can be derived analytically from the Kramers theory generalized to many dimensions, resulting in (Suzuki & Dudko, 2010)

(9) $$k_{{\rm 2D}} (F) = k_0 \Delta q(F)\,{\rm exp}\left[ {\beta {\Delta} G^{\ddag} \left( {1 - \Delta q^3 (F)} \right)} \right],$$

where Δq(F) is the distance in the Q-direction between the native well and the transition barrier, normalized by its zero-force value ΔQ ^‡:

$$\Delta q^2 (F) = 1 - \displaystyle{{2{\Delta}x^{\ddag}} \over {3{\Delta}G^{\ddag}}} F - \left( {\displaystyle{{{\Delta} K^{\ddag} /(3{\Delta} G^{\ddag} K_0^2 )} \over {1 - {\Delta} K^{\ddag}/K_0}} - \displaystyle{{{\varepsilon}^2 {\Delta}Q^{\ddag 4}} \over {9{\Delta} G^{\ddag 2}}}} \right)F^2.$$

Here, K ₀ is the curvature of the native well in the x-direction at zero force, ΔK ^‡ is the difference of this curvature and the curvature of the transition state, and ε defines the shape of the dissociation pathway.

Interestingly, despite a simple functional form of the force-dependent rate in Eq. (9) and in striking contrast with the 1D picture, the minimalist 2D model, featuring only a single barrier, predicts a rich spectrum of responses of the molecule to force. Each response is characterized by the corresponding dependence of the rate on force. In one such response (Fig. 7, right), the lifetime (equal to the inverse of the rate) is predicted to depend non-monotonically on force, with the conformational transition – quite counter-intuitively – being slowed down by low force before it becomes accelerated by higher force.

We can identify the mathematical origin of this rich behavior by inspecting the expression in the exponent of Eq. (9): this expression can change its sign as the value of the force F is changed. In physical terms, the expression in the exponent of Eq. (9) describes the change in the barrier height at force F relative to the barrier height at zero force. The non-monotonic dependence of the barrier height (and hence the lifetime) on force is a manifestation of the natural property of the transition state to move with respect to the native well under force. Snapshots of the 2D energy landscape in Fig. 7 (top right) at three values of the external force illustrate how the force-induced motion of the transition state relative to the native well causes an increasingly favorable alignment of the dissociation pathway along the direction of the force. This property of the energy landscapes under an external force is discussed in Section 3 within a 1D picture and is now generalized to two dimensions.

5. Conclusions and outlook

This review summarized some of the recent analytical theories that address the challenge of interpreting the mechanical fingerprints, or signatures in the force–extension curves, of biological macromolecules. The inherent limitation of analytical approaches is that many details of the biomolecular system have to be overlooked to make the problem analytically tractable. However, when analytical treatment is indeed possible, such approaches have the power to reveal the general principles that unify seemingly unrelated biomolecules and their behaviors in different experimental situations. Furthermore, analytical expressions for experimentally measurable quantities are particularly well suited for a direct fit of experimental data. Such a fit results in concrete predictions of specific mechanisms that are instantiations of the general principles captured by the theory.

Single-molecule pulling experiments are rapidly expanding their scope to more complex biomolecules and multi-component biomolecular assemblies. More complex systems typically entail more complex behavior, thereby presenting new fascinating physics problems of biological significance. On the other hand, complexity usually hampers analytical tractability. While full analytical expressions for experimental outputs – for example, the probability distribution of transition forces – may no longer be possible, one may be able to analytically describe distinct signatures of these outputs that can be observed in the experiment. Examples of such distinct observable signatures include a multi-modality (multiple peaks) in the force distribution or a non-monotonic force profile of the transition rate. These distinct signatures, captured in a form of analytical expressions that are applicable to a broad range of systems, can serve as diagnostic tools to help identify the underlying mechanism of biomolecular behavior among several candidate mechanisms.

Although single-molecule force experiments are a powerful tool for probing biomolecular behavior with an unprecedented resolution, they suffer of an inherent limitation: they typically only monitor molecular movements in the direction of the force application. Advancements in instrument development are making it possible to overcome this limitation by the incorporation of a second measurement mode, fluorescence, so that multiple variables of motion can be monitored simultaneously. Hybrid force–fluorescence spectroscopy experiments have the unique capacity to exert picoNewton forces via single-molecule manipulation while reporting concurrent structural rearrangements by fluorescence, thus combining the advantages of two principally different single-molecule approaches in a single assay (Hohng et al. 2007; Hugel et al. 2007; Lang et al. 2003; Shroff et al. 2005; Zhou et al. 2011). These new experimental capabilities require a theoretical foundation to quantitatively correlate the dynamics of multiple simultaneously measured variables. Analytical frameworks of hybrid force–fluorescence spectroscopy will enable the multidimensional reconstruction of the energy landscapes and the important timescales of complex biomolecular assemblies.

Ultimately, we want to understand single-molecule behavior in the native environment of these molecules – in the living organism with all its complexity. The development of analytical approaches to more complex and less well-understood instances of biomolecular behavior will build a theoretical foundation for extending single-molecule studies to the domain of utmost complexity – molecular behavior in vivo.

A major strength of analytical theories lies in the balance between their abstraction, which makes a complex behavior simple enough to help develop intuition about it, and their power to generate concrete, experimentally testable predictions. It is my hope that this review will inspire the reader to further explore the intuitively appealing nature of analytical–theoretical approaches and their power to reveal the simplicity of the physical principles underneath the complexity of biological phenomena.