Interpreting the widespread nonlinear force spectra of intermolecular bonds

Single Bond.

In a DFS experiment the force transducer [atomic force microscope (AFM) cantilever, biomembrane force probe, or an optical trap] typically exhibits a Hookean response to displacement, and thus if a single bond is characterized by a single barrier, the combined bond/transducer system comprises two states: the bound state, b, when the bond is formed, and the unbound state, u, when the bond is broken and the system fluctuates in the parabolic well of the transducer. If we assume that the time scale of intra-well relaxation is faster than any other time scales of the system, then the probability of finding the system in the bound state, p_b(t), is described by the master equation of a two-state Markov process (8), [1]where the instantaneous unbinding, k_u(t), and binding, k_b(t), rates are time-dependent to reflect their dependence on force.

The transition rates used in Eq. 1 must adequately describe the dynamics, yet they do not need to be complicated for modeling typical laboratory experiments. Several recent publications (20⇓⇓–23) showed that it is possible to produce nonlinearities in the force spectrum at very fast loading rates under an assumption of a force-dependent transition state (which can result from certain model-based potentials). However, the loading rate required to explore this regime of nonlinearity is 4 to 7 orders of magnitude greater than the upper end of loading rates reached in the laboratory (see SI Text for details). Therefore, for this work we chose to use the simpler unbinding rate given by the Bell model (7), which estimates the effect of force on the barrier to be linear and depends on a minimal number of parameters [Eq. 2]. In this model, the rebinding rate is approximately the natural frequency of the transducer, , scaled by the Boltzmann-weighted energy of a spring extended between the barrier location x_t and the relative displacement of the spring minimum, f/k_c (see Fig. 1A). Thus the unbinding and binding transition rates are given by, [2][3]where is the intrinsic unbinding rate of the molecular complex (where the probe-molecule connection may contribute additional fluctuations and a modified damping), β^-1 = k_BT is Boltzmann’s constant times the temperature, and defines the equilibrium free energy ΔG_bu = G_u - G_b of the bond G_b, relative to the free cantilever G_u.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Potential energy profile, two-state kinetics, and comparison of the single-bond model force spectrum to the two-state master equation. (A) 1D energy profile of a single bond (blue curve) loaded by a harmonic potential (dotted green curve). The total potential of the system is represented by the thick black curve. At zero load (I), the harmonic and bond potential minima coincide and a single stable minimum exists. As the harmonic potential is pulled beyond the bond barrier (II) a second minimum is created and the system becomes bistable; however, the barrier impeding rebinding is negligible. Pulling further (III) the re-entry barrier increases and significantly reduces the chance of rebinding. The black dashed line represents the pathway that is commonly assumed beyond the barrier. That is, once the system overcomes the barrier it is carried away irreversibly. (B) Transition rates for unbinding [Eq. 2] and rebinding [Eq. 3] using typical parameters of , x_t = 1 Å, ΔG_bu = 10 k_BT, and k_c = 100 pN/nm. Rates are in units of and force is in units of k_BT/x_t. The regimes corresponding to stages in the upper panels are partitioned for comparison. (C) The differential equation defining the two-state master equation [Eq. 1] is solved by numerical integration (symbols) with the Bulirsch-Stoer method using Richardson extrapolation (Igor Pro, Wavemetrics). Results are shown for ΔG_bu = 10 k_BT, and three different confinement force values f_k = k_cx_t, which govern the effects of the probing spring constant on the resulting spectrum. The single-bond model in Eq. 6 is calculated with identical parameters for comparison (solid curves).

The two rates in Eqs. 2 and 3 cross at a unique force, f_eq, given by (1, 24) [4]which defines the equilibrium force for the bond/transducer system. It can also be shown that for βΔG_bu≫1 the mean force in the limit of vanishing loading rate is also equal to f_eq. The existence of this limiting force—the lowest force required to break a bond for a given force transducer stiffness—has been predicted by a number of theoretical studies (8, 9, 15), and we have recently presented an experimental demonstration of this phenomenon (24). On a more fundamental level, the existence of this force reflects major differences between the dynamics of free bonds and bonds attached to a force probe that adds a confining potential to the system and creates a defined unbound state. Note that even polymer linkers between the molecule and the probe behave like an entropic spring over large fractional extensions when the contour length is significantly larger than the persistence length (25). This effect results in an effectively spring-like potential for small loads where near-equilibrium effects are important. The major underlying assumption of the linear first-passage model is that a free noncovalent bond will break on its own in absence of external force if left alone for a sufficiently long time. However, when the same bond is attached to a force transducer, the presence of the transducer potential precludes the bond from breaking until a minimal force is applied, at which point a second state is created. Given the existence of the two states, another profound consequence is that the overall dynamics of the system should now be dependent on the shape of the second well and thus on the transducer stiffness k_c in all regimes of bond rupture (26⇓–28).

Inspection of Eqs. 3 and 4 reveals that the rebinding rate can be expressed as a function of the unbinding rate and the equilibrium force f_eq. Surprisingly, in this case the complete two-state process described by Eq. 1 includes only one additional free parameter, f_eq, to the usual unknowns , and x_t. To understand the significance of f_eq consider the evolution of the transition rates with force (Fig. 1B). For forces less than f_eq the rebinding rate dominates the dynamics, acting to keep the system in the bound state. At forces beyond f_eq the rebinding rate falls off rapidly due to the quadratic term -f²/2k_c in the exponent, leaving the unbinding rate as the only significant contributor to the dynamics. Note that, although the distribution of rupture forces extends below f_eq, it is impossible to observe a mean force below f_eq. This important feature of the bond rupture process allows us to simplify the two-state master equation (Eq. 1) by starting the bond loading process at f = f_eq, instead of starting at f = 0 [Seifert (15) used the same approach for analyzing the two-state problem through a mean-field approach]. With this simplification we truncate Eq. 1 to, [5]and inserting Eq. 2 in Eq. 5 we find the mean rupture force , [6]where f_β = k_BT/x_t is the thermal force scale, and is the exponential integral. We can also produce an analytical approximation to Eq. 6 by using the relationship e^zE₁(z) ≅ ln(1 + e^-γ/z), where γ = 0.577... is Euler’s constant. The most significant feature of Eq. 6 is that it extends over the two primary regimes of the single-bond dynamic force spectrum: a linear regime in the slow loading limit, , and a nonlinear regime in the fast loading limit, , leading to the characteristic ln r dependence expected for the irreversible rupture of single bonds. Indeed Eq. 6 is remarkably accurate: If we compare its predictions with the numerical solution of Eq. 1 we see negligible differences over a variety of parameters and loading rates (Fig. 1C), with the largest error of roughly 5% occurring at the “elbow” of the spectrum where R(f_eq) ∼ 1.

Multiple Bonds.

We now turn to the system that comprises N bonds arranged in parallel between the probe and substrate. As we discussed in the introduction, this is arguably the common state of affairs in the laboratory (8). In a continuum approach, the number of formed parallel bonds N_b, of which N_t total bonds are allowed to independently unbind and rebind, is represented by the differential equation (7, 16⇓–18) [7]where k_u(f) is given by Eq. 2 and k_a is the constant rebinding rate of an individual molecule acting over a short distance from its retracted position while the overall cluster remains bound. This multiple bond model has been extensively examined numerically and by asymptotic analysis (7, 16⇓–18, 29⇓–31). Again, just as with the single bond case, two primary regimes emerge. At slow loading rates, fluctuations in the number of closed bonds is relatively fast (Fig. 2), and the system establishes a meta-stable number of bonds, N(f) = N_t/(1 + k_u(f/N)/k_a). From here we define N as the equilibrated number of formed bonds at zero force, . Stability of the cluster is completely lost when the force exceeds the equilibrium force (7, 17, 30) [8]where W(x) is the Lambert W-function, which is defined as the solution for W to the equation We^W = x, and defines the equilibrium free energy ΔG_ba = G_a - G_b of the closed bond G_b relative to an open bond G_a within the bound cluster. In systems where fluctuations in N_b are large, rupture can occur long before stability in N is lost, leading to a more gradual decrease in rupture force with decreasing loading rate than predicted here. The second primary regime appears at fast loading rates where the rebinding rate k_a, and hence fluctuations in N_b, become less important. Here the force reaches higher loads and the bonds rupture rapidly after the first bond fails, leading to the irreversible rupture force, (16, 17). Note that the force spectrum in this regime follows the same classic f ∼ ln r trend, therefore the resulting force spectrum will look qualitatively very similar to the single bond case. We obtain an interpolation between these two regimes (see SI Text for details): [9]Again, this equation is remarkably accurate in describing the force spectrum over the complete range of loading rates. When we compared the spectrum predicted by Eq. 9 to Monte Carlo simulation of the complete model in Eq. 7 over seven decades in loading rate (Fig. 2) using the Gillespie algorithm (32), we found that this model accurately captured the two primary trends, as well as the crossover region between the two loading rate regimes.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

The multi-bond model and comparison to Monte Carlo simulation of N-bonded systems. (A) Schematics of the two limiting pathways when several bonds are allowed to independently unbind and rebind, which are dictated by loading rate. At slow loading, the motion along the N_b coordinate is relatively slow, and rupture occurs once N_b diffusively crosses a threshold defined by a maximum in the free energy (18, 30). At fast loading, the force is quickly raised and thus the drop in N_b will be fast after the initial bond ruptures. In this case the spatial diffusion of one bond over its energy barrier is the slow process which dictates rupture of the cluster. (B) The mean rupture force for the indicated number of N_t bonds, as described by Eq. 7, is solved by the Gillespie algorithm for . Rupture forces for a given trajectory are determined when N_b = 0, and hundreds of trajectories are used to produce the mean rupture force (symbols) plotted here against the dimensionless loading rate. For comparison, the multi-bond model in Eq. 9 is calculated (solid curves) for the same ratio of transition rates, and with .

Interestingly, the multi-bond model of Eq. 9 shows that, if we combine the number of bonds N into an apparent force scale, [10]along with a generic definition of an equilibrium force f_eq, then the multi-bond force spectrum takes the same mathematical form as the single-bond model in Eq. 6. Eq. 10 shows that the multiple-bond model modifies the single-bond model simply through the parameter N, which factors inversely to the transition state x_t everywhere in the function. Significantly, this observation explains the unreasonably small transition state distances reported in many force spectroscopy experiments. When an N-bonded system is analyzed in the kinetic regime using the common-practice single-bond model, the fitted apparent transition state, , will be N times smaller than the true distance x_t.

Comparison to Experimental Data.

These models provide a theoretical basis for extracting a consistent set of interaction parameters of a bond over the broader context of loading conditions, not just those that satisfy the far-from-equilibrium condition. To test the validity of our model experimentally, we measured rupture forces for the nickel-mediated bond between nitrilotriacetic acid (NTA) and poly-Histidine (His₆) over four decades of loading rates (see SI Text). As Fig. 3A demonstrates, the mean force with log-loading-rate displays the expected two-regime behavior. The fit of Eq. 6 to the data is excellent, and we obtain reasonable values for all parameters: x_t = 0.92 ± 0.15 Å, and . The kinetic unbinding rate of 4.3 s^-1 is of the same order as 1.8 ± 0.40 s^-1 found for a single Ni-NTA unbinding from His₆ by time-resolved fluorescence (33). The same system was previously investigated by another group (34), with the only difference that the probe-surface linkages of His₆ and NTA functionalities were inverted. Those measurements (Fig. 3A) appear to diverge from ours. Attempts to corroborate the two data sets using the multi-bond model, through varying N, were unsuccessful. However, when we replace the value for the spring constant in our experiments (38 pN/nm) with that used in the previous study (10 pN/nm) and use exactly the same parameters obtained from our data, the force spectrum predicted by the single-bond model of Eq. 6 provides an excellent fit to these data points. This result highlights the significance of the transducer stiffness on the resulting force spectrum (9, 24, 26, 27) and demonstrates that proper analysis can unify the force spectra collected under different experimental conditions.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Model fits to the nonlinear force spectra of intermolecular bonds. (A) Force spectrum of the Ni-NTA/His₆ bond measured in this work along with the data of Verbelen et al. (34). Measurements made without Ni²⁺ demonstrate the specificity of the interaction (open circles). Solid lines represent fits to Eq. 6 using identical parameters for both data sets except for the respective spring constants. (B) Force spectra of 10 data sets taken from the literature are fit to Eq. 9 assuming a generic equilibrium force and apparent transition state distance. Data are exploded along the loading rate axis for clarity. Inset: The same force spectra in their raw form illustrating the theoretical range of equilibrium force and cross-over loading rate span (shaded regions). The upward inflection at high forces in the Biotin/Avidin data of ref. 2 may be due to limited sampling rate effects (41) and enhanced force error at very fast loading rates (42). (C) The same data in B plotted in the natural coordinates of Eq. 9 (see Eq. S3) show that all spectra collapse onto a single line. (D) Biotin-avidin bond rupture data of Teulon et al. (35) are globally fitted to Eq. 9 assuming N = 1, 2, and 3 parallel bonds. Only f_eq is independently fit, while x_t and are shared. Fitted values are x_t = 0.78 Å, , and f_eq = 24.6 pN (N = 1), 58.5 pN (N = 2), 142.3 pN (N = 3). Legend in (B) refers to references and corresponding bonds as follows: Biotin/Avidin (2); LFA-1/ICAM-1 [rest 3A9] (10); Aβ-40/Aβ-40 (11); N,C,N-pincer/pyridine (12); Si₃N₄/Mica in Ethanol (14); peptide/steel (43); Integrin/Fibronectin (44; Lysozyme/Anti-Lysozyme (45); Dig/Anti-Dig (46); Actomyosin/ADP (47).

The deeper implications of the model in Eq. 9 become clearer when we apply it to DFS data for a wide range of systems. In addition to our data, we have compiled force spectra from 10 different studies, each probing a different intermolecular bond and extending over a sufficient span of loading rates to exhibit the two dynamic regimes. Eq. 9 fits the full range of all data sets very well, even for those datasets in which the lowest value of the pulling rate is still too large to fully reach the plateau of the near-equilibrium regime (Fig. 3; see also Table 1 for the fitted parameter values). Moreover, all data collapse onto a single universal line when plotted using the natural coordinates of that equation (Fig. 3C; and see Eq. S3). We also compare our model to the outstanding analysis of DFS measurements for multiples of biotin-avidin bonds by Teulon et al. (35). To be valid, the parameters predicted by Eq. 9 should remain constant regardless of the number of bonds, N. In Fig. 3D we present global fits to the three curves of Teulon et al. (assuming N = 1, 2, and 3, and keeping x_t and unchanged for all datasets) and find good agreement with the data for all three curves.

As predicted by the multiple bond analysis, a number of spectra we examined produce apparent transition state distances that are unrealistically small, for example 6.3 × 10^-4 Å, 1.3 × 10^-3 Å, 1.5 × 10^-2 Å, and 1.4 × 10^-2 Å from refs. 13 and 36⇓–38, respectively, with several on the order of 10^-2 Å (see Table S1 for complete list of parameters). Note that these values are very similar to those given in the original publications using the linear f ∼ ln r model on the high-force regime. Potential of mean force calculations generally find the distance from the minimum to the primary barrier to be on the order of 1 Å (39). With this rough estimate the apparent distances found on the order of reflect bond numbers around N = 100, which is the same number of bonds predicted when assuming thiol monolayers on Au-coated probes with 100 nm tip radii using JKR contact mechanics (40). We note that if our model is to be used to extract true transition state distances from experimental data from a multi-bond system, the number of bonds will need to be determined independently. Alternatively, independent methods of ensuring that measurements probe only single bond interactions (8) can be used prior to applying the single-bond model to the data. However, the intrinsic unbinding rate of a single bond , and the single-bond free energy between the bound and closely unbound states ΔG_ba, can both be determined without knowledge of N. Further discussion of the correlation between equilibrium force and transition state as it applies to bond valency is provided in SI Text (see Fig. S1).

What other arguments support or refute our interpretation? It is logical to assume that if rebinding plays a significant role in the bond rupture at slow loading rates then we should be able to observe these events in the raw data. However, this is not a trivial phenomenon to observe, because the rare sojourns in the unbound state may not be long enough to be observed with typical sampling rates on the order of 1 kHz. Hence fast sampling rates should be used at slow loading rates to confirm or discard near-equilibrium effects. In the multi-bonded scenario the question regarding rebinding events becomes irrelevant, because macroscopic rebinding from the completely unbound state is not required to establish the quasi-steady-state number of bonds and therefore the equilibrium force.

Finally, we consider the conditions required to observe nonlinearity within this model—that is, the loading rate that marks the crossover between near-equilibrium to kinetic regimes. Inspection of Figs. 1 and 2 shows—as a consequence of the exponential integral appearing in Eq. 9—a useful relationship for benchmarking where this crossover occurs is when R(f_eq/N)/N = 1, or, [11]Using the extremes of experimentally observed values for f_eq ∼ 10–200 pN, , with nominal x_t ∼ 1 Å and N ∼ 1–50, we find crossover loading rates encompass a range r_x ∼ 1–10⁵ pN/s. Remarkably, this range brackets the typical loading rates explored in a laboratory, as shown by the shaded area in the inset to Fig. 3B. Likewise, we show in SI Text that the theoretical equilibrium unbinding force for a single-bond falls in the range of f_eq ∼ 40 to 108 pN. Not only is this range above the noise limits of widely used AFM systems, but these forces can be quite high, and are remarkably consistent with the experimentally measured low-force regime data (see inset to Fig. 3B).