Calibrating optical tweezers with Bayesian inference

: We present a new method for calibrating an optical-tweezer setup that does not depend on input parameters and is less affected by systematic errors like drift of the setup. It is based on an inference approach that uses Bayesian probability to infer the diffusion coefficient and the potential felt by a bead trapped in an optical or magnetic trap. It exploits a much larger amount of the information stored in the recorded bead trajectory than standard calibration approaches. We demonstrate that this method outperforms the equipartition method and the power-spectrum method in input information required (bead radius and trajectory length) and in output accuracy.


Introduction
Optical tweezers, first constructed and used by Arthur Ashkin during the 1970s and 1980s [1], are widely used for single molecule manipulation [2]. The movement of a micron-sized dielectric bead in the potential well produced by a focused laser is used to detect forces on the pico-Newton level and displacements down to nanometer resolution. Applications range from studying protein folding [3][4][5], and elastic properties of DNA [6][7][8] to molecular motor function in vitro [9] or even in vivo [10][11][12]. As an experimental tool, it has become highly valuable to biophysicists and biologists working at the single-molecule level. More recently, magnetic tweezers allowing multiplexed measurements have been introduced and applied to a variety of biological questions [13,14], and further innovation in the field of optical tweezers is continuously taking place [15,16].
Before quantitative results can be extracted from an optical or magnetic tweezer experiment, however, the trap must be calibrated. Its effective spring constant needs to be determined to relate the bead displacement within the trap to the corresponding force experienced by the bead. To this end, a variety of methods have been developed using the same beam as the optical trapping beam to record the bead trajectory inside the trap and extract the trap stiffness from it. The most commonly employed techniques are the equipartition method [17] and the power-spectrum method [18,19]. Alternative methods have been proposed, such as the drag-force method [20], the escape-force method [21], and the step-response method [22]. Moreover, two-beam approaches have been introduced using one laser beam to trap the bead and a second low-power one to measure its displacement inside the trap [23]. This configuration is more complicated to implement and to align but is interesting for cases where the absolute bead location has to be known or for multiple-trap experiments.
The Brownian motion of the bead inside the confining potential can be recorded either with a quadrant photodiode (QPD), with position sensitive detectors (PSD) or with a camera [17]. For stiff traps, large bead displacements take place very fast and the acquisition rate must be fast enough to correctly capture the bead motion. Although acquisitions with a camera allow simultaneous recording of multiple trapped beads, they are limited by the image readout rate to the calibration of relatively shallow traps. We here show that single-beam setups can be more rapidly and more accurately calibrated by better exploiting the information stored in the recorded bead trajectories with Bayesian inference. We will demonstrate this for bead displacements recorded with a QPD, the results however remain valid independently of how the bead displacement is recorded.
The equipartition method calculates the mean squared displacement of the bead in one dimension, 2 x , and uses the equipartition theorem to obtain the corresponding spring coefficient in this dimension: 2

22
Bx kT= k x (1) where x k is the trap stiffness in the x-direction, T the temperature and B k the Boltzmann constant. However, this method suffers strongly from drift during acquisition, which leads to an increase of the apparent mean squared displacement and thus to an underestimation of the effective spring constant of the trap.
The power-spectrum method, on the other hand, exploits the stochastic motion, x(t), of the bead in one dimension versus time and, in particular, its Fourier transform, FT[x(t)]. The resulting data is fitted with the following Lorentzian: Here, 6 is the drag coefficient, the liquid viscosity, and a the bead radius. c f represents the cutoff frequency related to the spring constant as follows: . 2 x c k f (3) Displacement in stiffer traps will contain more high-frequency components that will raise the cutoff frequency. One can readily see that, in this case, calibration depends on precise knowledge of several parameters, in particular particle size and viscosity of the surrounding liquid whose effective value depends on the presence of nearby surfaces [24].
The escape force and drag force method suffer from similar constraints. A force is applied by displacing the liquid surrounding the trapped bead (with a piezo translation stage, for instance) producing a drag force according to Stokes' law, 6 drag F av where v is the fluid displacement speed [20,21,23]. This force is equated to the restoring force, i.e. the product of the beam displacement from the trap center and the effective spring constant. As in the power-spectrum method, the viscosity and the bead radius a are required in addition to the velocity v of the fluid relative to the bead. Furthermore, lateral displacement during drag also induces axial displacement of the bead [20] perturbing the effective spring constant and further complicating the calibration attempt.
The principle of the step-response method is to abruptly displace the trap and observe the response as the bead moves back into equilibrium towards the center of the trap [22]. The restoring motion is described by a characteristic time constant that can be related to the trap stiffness k by 6 a k . Again, precise knowledge of additional parameters i.e. the bead size and the viscosity of the surrounding liquid is required to apply this method. We here propose a new method to calibrate an optical trap based on a statistical Bayesian inference approach applied to the bead trajectory inside the trap that does not suffer from the above drawbacks. We have previously shown that this technique is ideally suited to extract the confining potential felt by membrane receptors diffusing inside cell membrane microdomains [25][26][27][28]. The stiffness of the potential in that case is very low (3 4 pN/nm). Moreover, Bayesian inference was used to extract the viscoelastic properties of the medium between two traps in a double-trap setup [29]. We here show that this approach is applicable to the much stiffer potentials encountered in optical-tweezer experiments with stiffnesses of 10 2 -10 1 pN/nm. This inference approach circumvents the drawbacks of single-beam calibration techniques: it is less influenced by minor drift during acquisition, does not require any prior knowledge of the experimental parameters, can detect deviations from the assumed 2nd order spring potential, and can obtain calibration results much faster (i.e. requires much less trajectory points). All that is required for the calibration is a recorded trajectory of the bead within the trap a few hundred points long with an appropriate time step between consecutive points.

Bayesian inference
We assume that the movement is dictated by the overdamped Langevin equation: (4) where r(t) is the position at time t, V(r) the potential at r, the term representing stochastic Brownian motion, and D the diffusion coefficient, with the Einstein-Stokes fluctuation-dissipation theorem relating D and by 6 BB kT kT D a . In our previous work, we showed that it is possible to assume and infer a position-dependent friction and diffusion coefficient [25,28]. In the present case of a bead in water, however, there is no reason to introduce such a position dependence and we therefore use a constant and D.

P( ,t | ,t ) rr
is the probability of going from space-time point 00 () ,t r to space-time point () ,t r . Splitting the trajectory into subdomains S ij in which the potential gradient () ij V r is constant, allows solving Eq. (5) and obtaining the probability of going from one space-time point to another given the diffusion coefficient D and force () ij V r : Dt t tt rr F rr F (6) F ij is the force vector in the subdomain with row value i and column value j, and is the positioning noise of the acquired data. The size of the subdomains subdomain L depends on the domain size and the average bead displacement during the time interval between trajectory points (L step 0.006 µm) [25] and was chosen such that 1 25 We can then calculate the probability for the entire trajectory T, where r µ and r µ + 1 represent two consecutive points within subdomain S ij . Then, by using Bayes' theorem we can find the most likely motion parameters given the observed trajectory, i. e. the parameter values that maximize the posteriori probability We therefore assume a second-order potential and infer its coefficients C x , C y , C xx , C xy and C yy (C is not inferred and is assumed to be 0), as in [27], rather than inferring the individual force values in each subdomain, as in [26]: 2 and 2 .
xx x yy y kC kC (10) Note that, in this case, the algorithm still optimizes the potential derivative (force) values in each subdomain, only these values are not independent but governed by Eq. (9). Typically, only C xx and C yy take on values of significant magnitude, confirming the spring-like potential of the trap and all the other terms can be neglected. However, assuming a full second-order potential allows the detection of any notable deviations from the ideal spring potential [30].
x k and y k represent the spring constants related to the restoring forces in the x and y directions, respectively. Note that all other calibration methods typically neglect the linear and xy cross terms and extract directly x k and y k [18]. Figure 1 shows an experimental trajectory obtained for a laser power (measured at the microscope entrance) of 377 mW and the inferred confinement potential. We define the dimensionless confinement factor tt t is the delay time between two trajectory points. The confinement factor u represents the fraction of the domain size covered by the displacement during one trajectory step. Ideally, u should be as small as possible and, in all cases, much smaller than 1. Depending on the value of u, we have shown that the inference algorithm may somewhat overestimate the spring constant k and underestimate the diffusion coefficient D [25]. This bias is due to the fact that, in Eq. (4), the potential term and the diffusion term can compensate for slight deviations of each other. Moreover, when the displacement between two trajectory points starts becoming comparable to the size of the domain, some displacements will be underestimated because the delay between successive recordings of the position is not small enough and the particle will bounce off the potential barrier during this delay t . Thus, the full displacement may not be captured from one recording to the next. This leads to an underestimation of the diffusion coefficient and, consequently, to an overestimation of the confining potential. This bias is deterministic and can be corrected for by performing a set of simulations for parameter values in the range of the experimental conditions to yield a bias curve that gives the relevant bias value for each confinement factor value. Figure 2 shows the bias of the motion parameters, D and k x , as a function of the confinement parameter u. The obtained bias curves can be fitted with Eq. (11) and (12) and have been used in the following to correct the inferred parameter values.
The subscript 'inf' denotes that this is the value inferred by the Bayesian inference algorithm.

Experimental setup and methods
The optical set-up has been described elsewhere [22]. Briefly, it is built using an inverted microscope (Olympus IX70) and an oil-immersion objective (Olympus PlanApo 60X, NA = 1.45). A Nd:YAG laser (Quantum Laser, model Forte 1064, TEM00, 1W cw) was used both for trapping and bead detection. The laser beam is diffracted by an acousto-optic deflector (AOD) (Intra Action Corp. DTD-274HA6) conjugated with the back focal plane (BFP) of the objective. The first-order diffracted beam by the AOD was enlarged to fill the pupil of the objective. A feedback system was used to keep the trapping laser power constant at the entrance of the objective during the experiments. Transmitted light was collected with a high numerical aperture condenser (Olympus Aplanat Achromat, NA = 1.4) and directed to a quadrant photodiode (QPD) (SPOT-9DMI, OSI Optoelectronics) located in a plane conjugated with the BFP of the microscope objective. The cutoff frequency of the QPD is on the order of a few kHz (5-10 kHz depending on the incident laser power). The 4 signals of the QPD were digitized simultaneously (sampling rate of 65536 Hz) using a Delta Sigma DAC (National Instrument, PCI 4474) and further processed using LabView 8.2. Displacement signals were normalized by the sum of all quadrant values before trap calibration. Silica beads The QPD calibration was performed using the so-called step-response method [22]. As the QPD is located in a plane conjugated with the condenser BFP, it is only sensitive to motions of the bead relative to the trap center and not to absolute motions of the trapping laser. Thus, a rapid displacement of the trap using the AOD corresponds to a spike on the QPD signal due to the fact that the trapped bead cannot follow that motion instantaneously. After AODmotion calibration, one can easily and accurately calibrate the QPD V to µm conversion.
Two series of experimental trajectories were recorded with 45 000 points for each trajectory. In the first series, the laser power was varied and the trap was located sufficiently far from the coverslip so that the effective viscosity seen by the bead is that of bulk water. In the second series, the laser power was kept constant and the distance from the coverslip was varied by displacing the microscope objective.
Numerical simulations and inference analysis were conducted in largely the same manner as in [25]. Each Brownian step in x and y directions was constructed using a Gaussian distribution of standard deviation 2Dt with D t chosen in accordance to t, each particle takes 10 000 substeps that are not averaged. The displacement due to the force corresponding to the confining potential, Only every sixth point of the experimental trajectories was used for the inference to remove the non-instantaneous QPD response. Inference analysis was conducted on a PC (dual-core 3.4 GHz, 4 GB RAM) using C language. Inferring the potential and diffusion coefficient from a trajectory with 10 00 points takes approximately a few tens of seconds.
The power-spectrum analysis was done with the MATLAB algorithm described in [33,34]. For the short trajectories of Fig. 6, in order to obtain a cutoff frequency, the number of data points per block had to be chosen below the values recommended in [33,34].

Bayesian inference analysis of experimental trajectories
Before applying our approach to experimental trajectories in optical traps, we will take into account the fact that, in most cases, four-quadrant photodiodes with a finite cut-off frequency are used to acquire the bead position. Thus, the recorded data suffer from a "memory effect" due to the non-instantaneous response time of the photodiode [18]. The response function of the photodiode, g(t), can be described by an instantaneous response fraction () diode and a noninstantaneous term with a characteristic decay time .
gt t e (13) The detected signal is thus given by: where S(t) is the actual signal and S (det) the recorded signal. The experimental data were recorded every 15.3 µs. Given that the cut-off frequency of the photodiode is a few kHz, we consider that by taking every sixth point of the trajectory we are effectively eliminating the memory effect as the non-instantaneous contribution due to the previously recorded trajectory point will have decayed sufficiently. We verified this with simulated trajectories where the non-instantaneous response was added. Taking every 6th point of these trajectories, effectively eliminated the non-instantaneous response and lead to correctly inferred stiffness values (data not shown). Obviously, the delay time between successive trajectory points will then be six times larger. This is the reason why the bias calculation [ Fig. 2 Experimental trajectories were recorded for six different laser powers and six different heights from the glass coverslip surface. The inferred confinement potentials are shown in Fig. 3. As expected, the parabolic potential becomes stiffer as the laser power increases. On the other hand, as the distance of the bead to the coverslip surface becomes comparable to its radius, the effective viscosity felt by the bead increases [18,24,35]. This leads to a decrease of the diffusion coefficient and, therefore, the trap stiffness determined is overestimated, as can be seen from the two terms of Eq. (4). This is indeed what is observed in Fig. 3(b). Note that, for large distances from the coverslip surface, the trap stiffness may also be affected by spherical aberrations which set in because of the refractive index change at the coverslipwater interface [36].

Comparison with the equipartition and the power-spectrum methods
The stiffness values obtained in the x-direction (k x ) with Bayesian inference, the equipartition method and the power-spectrum method are shown in Fig. 4 for varying laser power values (A) and for varying distances from the coverslip (B). Both the inference and the equipartition method yield similar stiffness values. The power-spectrum method, on the other hand, finds consistently higher stiffness values than the inference approach and the equipartition method. This difference is probably caused by the fact that the power-spectrum method requires input parameters, in particular the radius of the bead which was taken equal to the nominal value of the supplier, 500 nm. To confirm this, we measured the hydrodynamic radius of the beads using dynamic light scattering. The light scattering curve yields a maximum for a bead radius of 436 nm, which is indeed lower that the nominal value, and has a full width at halfmaximum of 184 nm. Based on Eq. (3), an overestimated bead radius yields an overestimated stiffness value and explains the higher stiffness values found by the power-spectrum approach. Fig. 4. Spring constants extracted from experimental bead trajectories with the Bayesian inference (black data points), the power-spectrum (red) and the equipartition method (blue) as a function of laser power (A) and distance from the coverslip (B). The solid lines are a guide to the eye. The dashed black line in (A) is a linear fit to the spring constants extracted with Bayesian inference (BI). Error bars for the power spectrum data are derived from the error in the Lorentzian fit; for the Bayesian inference data, the error bars indicate the width of the posteriori distribution and, for the equipartition data, the error bars are the error on the mean of <x 2 >. (C) Diffusion coefficient as a function of the distance from the coverslip inferred from the data set in (B) (orange), effective viscosity normalized to the bulk water viscosity 0 (purple) and diffusion coefficient after correcting for the effective viscosity effect (green). The green dashed line represents the resulting diffusion coefficient: 0.63 ± 0.02 µm 2 /s (the error is that of the fit with a constant value).
In contrast, the Bayesian inference approach does not require any input values. Moreover, the Bayesian inference method also yields the bead diffusion coefficient which can then be used to obtain the bead radius. The average diffusion coefficient found by the inference approach for the experimental data in Fig. 4(a) is 0.62 ± 0.01 µm 2 /s. Using the Einstein-Stokes relation and the water viscosity at 20°C, 3 10 Pa s , we find a radius of 358 ± 4 nm. Given the polydispersity of the bead solution and the fact that dynamic light scattering tends to overestimate the contribution of larger particles [37], the bead radius found with the inference approach is perfectly compatible with the dynamic light scattering measurement.
The bead radius may also be extracted from the y-intercept of the power spectrum curve for zero frequencies, from which the diffusion coefficient is calculated and, hence, the radius. However, the data at small frequency values are usually noisy and do not allow a precise diffusion coefficient and radius determination. Indeed, the radius obtained from the power spectrum in this manner was 154 ± 5 nm, which is much lower than both the dynamic light scattering and the inference values. The use of the bead radius determined from the dynamic light scattering data is not very helpful for a precise stiffness determination either. Indeed, the presence of polydispersity does not allow determination of the radius of the bead used in the actual experiment. Gosse and Croquette have proposed an alternative way to determine the bead radius based on the fact that the power spectrum of velocity fluctuations asymptotically reaches a value proportional to the diffusion coefficient at high frequencies [13]. However, this asymptotic value can only be measured if the cutoff frequency of the acquisition system is much higher than the one of the trap. In practice, for a trap stiffness on the order of 10 4 pN/nm (f c 2 kHz for a 1-µm bead), typical for optical tweezers, the cutoff frequency of the most commonly used QPDs is too low to allow determination of the diffusion coefficient and the bead radius in this manner.
When fitting the inferred values, we indeed obtain a linear relation between the laser power and the spring constant k x , as expected [ Fig. 4(a)] [1,21,38]. As discussed above, when the distance from the coverslip decreases, the bead sees an increased effective viscosity due to the proximity of the coverslip surface and, therefore, a decreased diffusion coefficient.  (15) where R is the bead radius and h the distance of its center from the solid surface. We indeed observe a decrease in the diffusion coefficient inferred with the Bayesian approach [ Fig. 4(c), orange]. If we correct the inferred diffusion coefficient by dividing with the ratio between effective and bulk viscosity (purple), we suppress the diffusion coefficient variation and obtain a fairly constant D value of 0.63 ± 0.02 µm 2 /s [green curve in Fig. 4(c)], as expected.
To examine the effect of drift on the equipartition method results, we generated numerical simulations with an increasing linear drift during the acquisition time [ Fig. 5]. We observe that experimental drift causes a much larger drop in the determined spring constant when using the equipartition method than when applying the inference method. Indeed, the inference analysis finds the correct spring constant with negligible bias (less than 5%) for drifts up to 10 nm at least, whereas the equipartition approach becomes increasingly biased already for drifts above 5 nm. An important additional advantage of the Bayesian inference approach is the fact that it needs a lot less points to obtain a meaningful result than the other methods. This may be useful because typically the trap calibration is performed after the end of the experiment. If the bead inadvertently escapes the trap, the ability to calibrate with only very few points may be crucial to exploit the experiment. Figure 6 shows the k x values determined with Bayesian inference and with the power-spectrum approach normalized to the input values for simulated trajectories with different numbers of points, N. As the trajectory lengths become smaller, the power-spectrum approach yields values with increasing uncertainty which become highly biased for very short trajectories, whereas the inference method can still accurately determine the spring constant even for trajectories with as few as 100 points. Furthermore, the algorithm for the power-spectrum approach we used, described in [34], increasingly fails to determine a trap stiffness value as the length of input trajectories decreases. By N = 500, the algorithm fails for about one quarter of the input trajectories. For the Bayesian inference approach, the posteriori distributions for k x [ Fig. 6(b)] remain narrow and unbiased for trajectories of lengths down to 1000 points and only become broader and somewhat more biased for trajectories of 600 points. In the case of the shortest trajectories of 600 points or less, to avoid wasting information and using only 100 points, as is the case for the posteriori distribution shown in black, we still used every 6th point to remove the non-instantaneous QPD response, but also used all sets of intermediate points sampled every 6th point and combined them consecutively to form a trajectory of the full length [blue curve in Fig. 6(b)]. This trajectory then has the original number of points but individual values remain decoupled from the previous values. Posteriori probability distributions of kx for two numerical trajectories generated with D = 0.3 µm 2 /s and kx = 0.06 pN/nm with N = 600 and 3000 points. For N = 600, if only every 6th point is taken into account, the posterior distribution is given by the black curve which is broad and more biased. If all 600 points are taken into account (see text), the blue posteriori distribution is obtained which is narrow and peaks at the same value as that corresponding to the longer N = 3000 trajectory (red). The small bias observed with respect to the input value of kx = 0.06 pN/nm (vertical dotted line) is corrected by using the bias curves of Fig. 2. Finally, the Bayesian inference approach may also be used to determine deviations from the expected parabolic profile. In addition to detecting the presence of linear or cross terms, the inference method can determine if a fourth-or sixth-order potential is more appropriate to describe the confinement, that is detect the presence of higher order terms, which may prove useful to correct aberrations in holographic optical traps [39] or when the bead explores areas away from the central part of the trap where the potential was found to be anharmonic [30]. In the case of our experimental trajectories, this was not the case: no matter the order of the potential assumed for the Bayesian inference, the second order coefficients do not change and additional higher order coefficients are found to be negligible [ Fig. 7]. However, this property of the Bayesian inference approach may also be helpful during the alignment of an optical-or magnetic-tweezer setup. Fig. 7. 2nd (red), 4th (dark blue) and 6th (light blue) order potentials inferred from an experimental trajectory obtained for a laser power of 251 mW (same trajectory as in Fig. 4(a)).

Summary
To summarize, the inference method is more accurate and reliable than both the other two most commonly used single-beam methods, the equipartition and the power-spectrum approach. Bayesian inference analysis of the recorded bead trajectories does not depend on external drift and does not require any input parameters. In contrast, the precision of the equipartition method is quite susceptible to drift during acquisition, while the power-spectrum approach requires the bead radius as input information, a parameter that is difficult to determine precisely with this approach. Furthermore, Bayesian inference needs the least amount of points to obtain a precise calibration. In addition, Bayesian inference yields the diffusion coefficient and therefore the bead radius. We here used an independent calibration of the QPD using the step-response method [22]. Alternatively, if the bead radius value is known, for example in the case of monodisperse bead solutions, we can use the inferred diffusion coefficient to calibrate the QPD which further facilitates the experimental procedure. In this work, we used experimental data from an optical-tweezer setup, however, the approach is directly applicable to magnetic tweezers and more generally to single-particle Brownian trajectories exploring all kinds of confining potentials [26][27][28].