Magnetic tweezers force calibration for molecules that exhibit conformational switching

Single-molecule magnetic tweezers (MT) provide a powerful means for measuring the elastic and conformational properties of individual biomolecules.¹ In a typical MT setup (Fig. 1(a)),^1,2 a single molecule is chemically modified so that it can be attached, at one end, to the surface of a glass flow cell and, at the other end, to a micron-scale paramagnetic bead. When a magnetic field gradient is applied across the sample, a magnetic force acts on the bead and stretches the molecule. Unlike other force spectroscopy techniques, MT instruments are able to apply inherently constant forces without the use of feedback routines.

The constant-force capability of MT is particularly well-suited to studying structured biomolecules. When the force applied to certain RNAs or proteins is maintained at a particular value, the molecules exhibit conformational switching between their folded and unfolded states, as illustrated in Fig. 1(b). Through measurements of equilibrium constants and state lifetimes, such studies give thermodynamic and kinetic information about the folding process.^3–6 Such data can be further analyzed to reveal changes in the number of ligands bound to a biomolecule as it folds/unfolds.^7–10 Additionally, in combination with the application of torque, studies can be made of DNA and RNA supercoiling transitions.^11,12

Conformational switching manifests as a telegraphic pattern in the time-series of molecular extension, z (Fig. 1(b)), and introduces two parameters into the system: the change in molecular extension between the two states, Δz, and the kinetic rate constant for switching between the states at coexistence, k_rate. Coexistence occurs at a particular applied force, F_1/2. Away from coexistence there will be multiple rate constants, one corresponding to the lifetime in each conformational state. Values of F_1/2, Δz, and k_rate for several simple biomolecular systems that exhibit equilibrium conformational switching are listed in Table I.

One of the main tasks in analyzing MT data is to calibrate the mechanical force experienced by the biomolecule. This is typically done via analysis of the Brownian fluctuations of the bead orthogonal to the direction of the pulling force (⁠$xˆ$ in Fig. 1(a)) by computing the variance,¹³ Allan variance,¹⁴ or power spectral density^14–18 of those fluctuations. In all of these cases, the force is obtained by treating the bead, tethered by the biomolecule, as a damped pendulum subject to Brownian fluctuations;¹³ the “spring constant,” κ, of the system for fluctuations in the $xˆ$ direction is obtained from the analysis of the horizontal bead fluctuations and the force, F, is found via

where is the time-average vertical extension of the molecule. Associated with these fluctuations is a characteristic corner frequency that varies inversely with extension. Accurate direct force calibrations for short molecules can therefore only be obtained when sampled at a correspondingly high frame rate,^14,15,17 e.g., using modern high-speed MT instruments.^19–21 Traditional MT instruments, with 60–100 Hz frame rates, will give significantly biased force calibrations for such molecules. In practice, short molecules are often handled by pre-calibrating the MT using a longer tether such as double-stranded DNA; however, this introduces issues of bead-to-bead variability for which careful corrections must be made.^22,23

Conformational switching may bias direct force calibration, since the Brownian fluctuations that are analyzed to obtain the force could be contaminated by additional fluctuations arising from the stochastic switching. These additional fluctuations are seen in Fig. 1(c), where conformational switching—clearly monitored in the z trajectory (Fig. 1(b))—leads to a non-constant amplitude of the x fluctuations. In Sec. II of this paper, we use theoretical arguments to show that conformational switching-induced bias can arise via two different mechanisms. Then, in Sec. III, we use simulations to verify these effects and to identify the portion of the experimental parameter space in which they are important.

MT force calibration can be accomplished by measuring the Brownian fluctuations of the bead/tether system. The fluctuations can be quantified by a single value, the variance,¹³ or by computing quantities that more completely use the available data: the Allan variance (AV) or power spectral density (PSD).¹⁴ The conceptually simple variance-based method assigns an average energy to each degree of freedom of the system using the equipartition theorem. The expression for the calibrated force, F_cal, of a single-state system is

where k_B is the Boltzmann constant, T is the temperature, and is the variance in the bead position orthogonal to the pulling direction.¹³ In the presence of conformational switching, each state has a different extension and variance ⁠. Assuming that the equilibration kinetics are sufficiently fast (see Sec. II B), Eq. (2) is shown also to hold for the switching molecule under the substitutions and (bar denotes an average over conformational states weighted by the time spent in each state).

A more robust calibration scheme utilizes the AV or PSD (in terms of a time difference, τ, or frequency, respectively), since variance-based calibration is particularly sensitive to noise. The AV and PSD are fit by functions having two free parameters: κ and α, the dissipation due to viscous drag on the bead. For example, the AV is well fit by¹⁴ 

A similar function exists for the PSD. Because the PSD and AV contain the same force calibration information, we will confine our discussion in this paper to the AV; however, our findings hold for the PSD as well.

Unlike in the simple case of Eq. (2), we find that the AV and PSD fitting functions are formally modified by conformational switching. This occurs because the power of the fluctuations at a particular frequency, corresponding to a particular τ, is the time-weighted average of the power in that channel across states. Therefore, the overall AV is the time-weighted average of the AVs of the constituent states (indexed by i) plus the power from the transitions between the states,

So, to fit the AV of an N-state system, it is necessary to fit the average of N copies of Eq. (3)—each with a different value of κ and weighted by the fraction of the total time spent in that state—and to include an additional term from the switching signal. Note that this switching refers to the subtle changes in the fluctuations in x (Fig. 1(c)), not the telegraphic pattern in z (Figure 1(b)). The simulations of Sec. III suggest that the switching contribution can be neglected under experimentally relevant conditions. If the extensions and occupancy times for each state are known, the multi-state Eq. (4) can be fit to the data without increasing the number of fitting parameters.

Our analysis thus far has neglected any effects that may arise from the kinetics of the system, either from the rate of switching between the conformational states or the rate at which the bead diffuses through the solution. In Sec. II A, we assumed that the fluctuations within each state were independent. However, if the lifetime within a particular state is too short, this assumption can break down. As an example, suppose that the system begins in a compact state and then transitions to an extended state: the bead must diffuse to larger distances from the origin in order to probe the larger fluctuations expected. If the system then transitions back to the compact state before this diffusion can take place, due to a sufficiently short lifetime of the extended state, these larger fluctuations will not be adequately sampled.

We can estimate when this diffusion limitation will occur by considering the interplay between the diffusion constant, D = k_BT/α, and the effective spring constant, /F. The spring constant sets the length scale of the fluctuations, δ_x, through Eq. (2). A simple model of diffusion can be used to obtain the expected displacement, Δx, in a particular time, Δt, given the diffusion constant: $Δx=2DΔt$⁠. This Δt is related to the kinetic rate constant for switching between the two states (assumed to be the same for both states, i.e., coexistence): Δt = 1/k_rate. We can expect to obtain biased force calibrations if the bead cannot diffuse far enough to sample the expected amplitude of the fluctuations in the time set by k_rate. Thus we can identify a critical rate constant, k_crit, above which the force calibration will become inaccurate,

This result is a function of F and ⁠, which are properties of the biomolecule under study (see values of F_1/2 and Δz in Table I), and of α, which depends on the magnetic bead size and solution viscosity. Eq. (5) assumes that the value of in the extended state is much larger than that in the compact state.

We simulated two-state bead trajectories using a Brownian dynamics approach based upon the over-damped Langevin equation governing this system,

where F_L is the stochastic Langevin force. The simulations were carried out iteratively at 300 kHz and then down-sampled, by averaging over bins, to f_s = 3 kHz to approximate the effects of video tracking with a finite frame rate. In all of our simulations, we set κ based on Eq. (1)—with F = 12 pN and depending on the state occupied—and set α = 1.5 × 10⁻⁵ pN s nm⁻¹: reasonable values for an MT experiment using micron-diameter magnetic beads. Transitions between states are modeled as stochastic events with k_rate and actualized by switching the spring constant between its two possible values. We model both the folding and unfolding transitions as having the same kinetic rate. This is the condition under which the effect of conformational switching is most pronounced: for unequal values of k_rate, the system tends toward one-state behavior. An example MATLAB script for carrying out these simulations and the fitting required for force calibration, and for computing estimated force calibration biases, may be freely downloaded at http://www.engr.ucsb.edu/~saleh/#Code.

To verify that the AV of a system with conformational switching behaves as the time-weighted average of the AVs of the individual states, we simulated a two-state system with extensions of 300 and 3000 nm and with k_rate = 5 s⁻¹. The computed AV from this simulation is shown in Fig. 2(a), with error bars representing the 68% confidence interval (akin to the one-sigma level) derived from the gamma distribution of the AV.¹⁴ Also plotted are the predicted curves from modeling either a one-state system with an effective, averaged spring constant (Eq. (3)) or a two-state system in which the two AVs have been averaged (Eq. (4)). The models are evaluated, without fitting, using the known values of the state extensions and of F and α given above. Comparing the two models, we see that the two-state model is in good agreement with the data, whereas the one-state model disagrees across a broad domain of τ.

To confirm the switching-rate-dependent effect of bead diffusion, we performed simulations over k_rate = {10, 100, 1000} s⁻¹, a range encompassing the predicted critical rate for the parameters of this system: k_crit = 533 s⁻¹ (from Eq. (5)). The AV obtained from these simulations is plotted in Fig. 2(b). There is little variation between the two curves below k_crit; however, above k_crit, the maximum in the AV shifts to lower τ, commensurate with a force calibration biased towards higher forces.

To quantify how conformational switching biases the force calibrations obtained in an experiment, we performed simulations over a range of k_rate (10–1000 s⁻¹) and Δz (0–1000 nm). Mirroring experimental practice, we treated the conformationally switching molecule as being attached to a handle of fixed length; here nm. We analyzed the resulting traces using both the effective one-state (Eq. (3)) and two-state averaged (Eq. (4)) AV fitting models. To fit the AV, we use a previously described maximum likelihood approach.¹⁴ Both models are fit with two free parameters (κ and α), since, for the two-state model, we take as known the fraction of time spent in each state; in experimental practice, these data can be obtained from telegraphic z-traces (Fig. 1(b)). At each pair of k_rate and Δz, ten simulations, each 10 s long, were carried out. The results of this analysis are shown in Fig. 3. The bias values reported are the percentage deviations of the obtained forces from the expected force, F = 12 pN; the uncertainties reported are the standard errors of the mean of these biases across the ten trials.

Both methods give an acceptably unbiased calibration of the force over most of the parameter space. There is no discernible difference in the results from the single-state and two-state models up to Δz = 1000 nm at small k_rate. Large biases are seen in the data with Δz = 1000 nm and k_rate = 1000 s⁻¹ for both the single- and two-state AV analyses, as expected based on the bead diffusion effect (compare k_rate with k_crit in Fig. 3). This effect is confirmed in simulations with smaller values of Δz and F, but having the same value of k_crit, in which appreciable bias is still seen. The bias does decrease for smaller values of Δz, presumably due to the breakdown of Eq. (5) when one state is no longer much more extended than the other.

While the results of these simulations are given in Fig. 3 in terms of Δz, an accurate interpretation of the underlying effects must also account for the handle length, ⁠. Because the bead diffusion effect (Sec. II B) depends on the time it takes the bead to diffuse to the limits of its range of fluctuations in the extended state, its magnitude will depend predominately on Δz, provided, as mentioned above, that ⁠. In contrast, the non-kinetic averaging effect requiring the use of the two-state AV model (Sec. II A) depends on there being a significant difference in the AV of the two individual states and will thus depend predominately on ⁠.

To assess the practical experimental implications of conformational switching, we monitored the folding and unfolding of a 14 base pair RNA hairpin using a high-speed MT instrument (f_s = 11 kHz)¹⁹ and then analyzed the resulting traces to obtain force calibrations. We used a molecular construct, akin to that of Ref. 10, consisting of a hairpin separated from the flow cell surface by an 828 base pair DNA spacer and from the magnetic bead by a 20 nucleotide length of single-stranded DNA. At a given bulk salt concentration, an appropriate force is applied to cause the system to enter folding/unfolding equilibrium between the two states, as in Fig. 1(b).

Figure 4(a) shows the AV computed from a representative trace near 50:50 folding/unfolding equilibrium and a fit of the effective single-state model (Eq. (3)). Because the AV and the fitting function are in good agreement, we conclude that the two-state correction is not required for such small values of Δz. This agrees with our simulation results (Fig. 3), which predict force bias values of <1% and standard deviations of ∼1% for both one- and two-state fittings.

Force calibrations using the one- and two-state AV models are shown in Fig. 4(b) for a representative molecule over a range of force and bulk salt concentrations. The calibrated forces, as a function of the position of the magnets, z_mag, relative to a reference position are fit by a function that approximates the force calibration behavior: F = A/(B − z_mag)⁴. The standard deviation of the calibrated forces about this curve is 5.6%, larger than the 1% deviation expected from the simulations. Because the single-state and multi-state force calibrations are in such good agreement, and because the model functions fit the data so well, we conclude that this scatter arises from sources other than the determination of κ.

We have carried out theoretical and simulation studies to determine whether conformational switching biases the calibration of force spectroscopy experiments. We have shown that perturbations to the expected behavior can occur due to the way the AV/PSD of the individual states must be combined (by averaging) and due to the competition between bead diffusion and conformational switching kinetics; the former effect can be corrected by using the appropriately averaged fitting function. We see from simulated (Fig. 3) and actual (Fig. 4) experiments that, in practice, conformational switching is not a dominant source of noise over the parameter space of contemporary biophysical studies, where Δz values tend to be less than 100 nm (Table I) and ⁠. In these cases, error in force calibration is dominated by other sources. Biases arising from conformational switching will become more pronounced in systems with faster switching kinetics, in solutions with higher viscosity, and in experiments employing larger beads or shorter biomolecular handles.

We thank Bob Lansdorp for his assistance with the high-speed magnetic tweezers experiments and Sarah Innes-Gold for helpful discussions regarding error bars in plots of the Allan variance. This work was supported by the National Science Foundation (Grant Nos. DMR-1006737 and DMR-1309414). D.R.J. was supported by an NSF graduate research fellowship (Grant No. DGE-1144085).