The effect of integration time on fluctuation measurements: calibrating an optical trap in the presence of motion blur

Dynamical instrument limitations, such as finite detection bandwidth, do not simply add statistical errors to fluctuation measurements, but can create significant systematic biases that affect the measurement of steady-state properties. Such effects must be considered when calibrating ultra-sensitive force probes by analyzing the observed Brownian fluctuations. In this article, we present a novel method for extracting the true spring constant and diffusion coefficient of a harmonically confined Brownian particle that extends the standard equipartition and power spectrum techniques to account for video-image motion blur. These results are confirmed both numerically with a Brownian dynamics simulation, and experimentally with laser optical tweezers.


INTRODUCTION
Investigations of micro-to nano-scale phenomena at finite temperature (e.g. single-molecule measurements, microrheology) require a detailed treatment of the Brownian fluctuations that mediate weak interactions and kinetics [1,2,3,4]. Experimental quantification of such fluctuations are affected by instrument limitations, which can introduce errors in surprising ways. Dynamical limitations, such as finite detection bandwidth, do not simply add statistical errors to fluctuation measurements, but can create significant systematic biases that affect the measurement of steady-state properties such as fluctuation amplitudes and probability densities (e.g. position histograms).
Motion blur, which results from time-averaging a signal over a finite integration time, can create significant problems when imaging fast moving objects. It is particularly relevant when measuring the position fluctuations of a Brownian particle, where even fast detection methods can have long integration times with respect to the relevant time scale, as we will demonstrate in this paper. Instrument bandwidth limitations that arise from motion blur affect a variety of fluctuation-based measurement techniques, including the quantification of forces with magnetic tweezers using lateral fluctuations [5], and microrheology measurements based on the video-tracking of small particles [6]. The issue of videoimage motion blur has recently been addressed in the single-molecule literature [7] and in the field of microrheology, where the static and dynamic errors resulting from video-tracking have been carefully analyzed [8,9]. However, discussion has been notably absent in the area of ultra-sensitive force-probes, despite the significant effect that it can have on quantitative measurements. This paper focuses on the practical problem of calibration an optical trap by analyzing the confined Brownian motion of a trapped particle [10,11,12,13,14] in the presence of video-image motion blur.
In this article, we present a novel method for extracting the true spring constant and diffusion coefficient of a harmonically confined Brownian particle that extends the standard equipartition and power spectrum techniques to account for motion blur. In section (2) we describe how the measured variance of the position of a harmonically trapped Brownian particle depends on the integration time of the detection apparatus, the diffusion coefficient of the bead and the trap stiffness. Next, this theoretical relationship is compared with both simulated data (section (3)) and experimental data using an optical trap (section (4)). Practical strategies for trap calibration are given in the discussion section (5), where we show that motion blur is not a liability once it is understood, but rather provides valuable information about the dynamics of bead motion. In particular, we show how both the spring constant and the diffusion coefficient can be determined by measuring position fluctuations while varying either the shutter speed of the acquisition system or the confinement strength of the trap.

BIAS IN THE MEASURED VARIANCE OF A HARMONICALLY TRAPPED BROWNIAN PARTICLE
Detection systems, such as video cameras and photodiodes, do not measure the instantaneous position of a particle. Rather, the measured position X m is an average of the true position X taken over a finite time interval, which we call the integration time W . In the simplest model, where both the measured and true positions of the particle have been expressed as functions of time t. More complex situations can be treated by multiplying X(t ′ ) by an instrument-dependent function within the integral, i.e. by using a non-rectangular integration kernel. We consider the case of a particle undergoing Brownian motion within a harmonic potential, U (x) = 1 2 kx 2 . In equilibrium, we expect the probability density of the particle position to be established by the Boltzmann weight exp(−U (x)/k B T ), where k B is the Boltzmann constant and T is the absolute temperature: The variance of the position should then satisfy the equipartition theorem, However, these equations do not hold for the measured position X m . In particular, motion blur introduces a systematic bias in the measured variance, Following standard techniques (e.g. [15]), the necessary correction can be calculated precisely as a function of the spring constant k, the friction factor of the particle γ, and the integration time of the imaging device W [7,8,9]. First, we define the dimensionless parameter α by expressing the exposure time W in units of the trap relaxation time τ = γ/k, i.e.
Note that α can also be expressed in terms of the diffusion coefficient D by using the Einstein relation γ = k B T /D, i.e. α = W Dk/(k B T ). Then as presented in appendix (A), the measured variance is given by: var(X m ) = var(X)S(α) (6) where S(α) is the motion blur correction function

NUMERICAL STUDIES
To verify Eq. (6) numerically, we use a simple Brownian dynamics [17] simulation of a bead fluctuating in a harmonic potential. For each time step ∆t, the change in the bead position ∆x is given by a discretization of the overdamped Langevin equation: where δx(∆t) is a Gaussian random variable with δx = 0 and (δx) 2 = 2D∆t, and the deterministic force f det = −kx corresponds to a harmonic potential as in our calculation. Motion blur is simulated by time-averaging the simulated bead positions over a finite integration time W . To minimize errors due to discretization, the simulation sampling time is much smaller than both W and Γ/m. Fig. 1(a) shows the simulation results for 3 different bead and spring constant settings as described in the caption. Agreement with the motion blur correction function S(α) of Eq. (7) is within the fractional standard error of the variance, ∼ 2/N [19]. We observe in Fig. 1(b) that the distribution of measured positions is a Gaussian random variable with variance var(X m ) < k B T /k, and is therefore fully characterized by the mean (trap center) and the measured variance calculated in Eq. (6). It is a Gaussian distribution as expected [16]. The normal curve with the predicted variance is superimposed showing excellent agreement. The expected distribution for an ideal "blur-free" measurement system is superimposed as a dotted line.

Instrument description
The optical trap is formed by focusing 1064 nm near-IR laser light (Coherent Compass 1064-4000M Nd:YVO 4 laser) through a high numerical aperture oil immersion objective (Zeiss Plan Neofluar 100x/1.3) into a closed, water filled chamber. Laser power is varied with a liquid-crystal power controller (Brockton Electro-Optics). This optical tweezers system is integrated into an inverted light microscope (Zeiss Axiovert S100).
The trapped bead is imaged with transmitted bright field illumination provided by a 100 W halogen lamp (Zeiss HAL 100). The image is observed with a high-speed cooled CCD camera with adjustable exposure time (Cooke high performance SensiCam) connected to a computer running custom data acquisition software [20]. Each video frame is processed in real-time to determine the position of the trapped bead. Fast one-dimensional position detection is accomplished by analyzing the intensity profile of a single line passing through the bead center. To increase the signal to noise ratio and the frame rate, the camera bins (i.e. spatially integrates) several lines about the bead center (32 in this experiment) to form the single line used in analysis. A third order polynomial is fit to the two minima corresponding to the one-dimensional "edges" of the bead, giving sub-pixel position detection with a measured accuracy of about 2 nm.
Tracking errors can be included in the measured variance by adding the parameter ε 2 to Eq.
ε 2 can be interpreted as the measured variance of a stationary particle, provided there are no correlations between the tracking error and the measured positions. A detailed treatment of such particle tracking errors is provided in reference [8], where the authors also stress the importance of using identical conditions of noise and signal quality when comparing the ε 2 parameter between different experimental runs. Care was taken to achieve these conditions as is described below.

Experimental conditions
The sample chamber was prepared with pure water and polystyrene beads (Duke Scientific certified size standards 4203A, 3.063 µm ± 0.027 µm). Experiments were performed by holding a bead in the optical trap and varying the power and the exposure time. The bead was held 30 µm from the closest surface, and the lamp intensity was varied with exposure time to ensure a similar intensity profile for each test. This ensured that the noise and signal quality between experimental runs was very similar, as validated by the results. For each test, both edges of the bead in one dimension were recorded and averaged to estimate the center position.

Experimental Results
The one-dimensional variance of a single bead in an optical trap was measured at various laser powers and exposure times. Low frequency instrument drift was filtered out as described in Appendix B. For each power, measured variance vs. exposure time data was fit with Eq. (9) to yield values for the spring constant k, friction factor γ, and tracking error ε 2 . Error estimates in the variance were calculated from the standard error due to the finite sample size, and variations due to vertical drift.
Error in the fitting parameters indicate that the best estimates for γ and ε 2 occur at the lowest and highest powers, respectively. These estimates both agree within 2% of the error weighted average for all powers. For the nominal bead size, γ agrees with the Stokes' formula calculation to within 11%, indicating a slightly smaller bead or lower water viscosity than expected. Additionally, the estimate of tracking error ε determined from the fit compares favorably with the standard deviation in position of a stuck bead, differing by about half a nanometer.
For a single bead observed under identical measurement conditions, γ and ε 2 are expected to remain essentially constant as the laser power is varied. Good consistency was found between the determined values of ε 2 from different experimental runs, due to the protocol of matching the signal strength between tests. While laser heating could cause γ to decrease with increasing power, this effect should be small for the < 500 mW powers used here [21,22], so this effect was neglected.
Holding γ and ε 2 constant for all powers, the raw data was re-fit with Eq. (9) to yield k. The data for all 4 powers was error-corrected by subtracting ε 2 and was rescaled according to Eq. (6) and Eq. (7). This non-dimensionalized data is plotted alongside the motion blur correction function in Fig. 2, showing near-perfect quantitative agreement. This exceptional agreement further validates our treatment of γ and ε 2 . A plot of spring constant vs. dimensionless power is shown in Fig. 3, demonstrating the discrepancy between the blur-corrected spring constant and naïve spring constant for different integration times. Even for a modest spring constant of 0.03 pN/nm and a reasonably fast exposure time of 1 ms, the expected error is roughly 50%. We also note that the blur-corrected spring constant increases linearly with laser power as expected from optical-trapping theory. Once confirmed for a given system, this linearity can be exploited to determine not only the spring constant as a function of power but also the diffusion coefficient of the bead. This is discussed in subsection (5.3), and presented in Fig. 3.
When the data acquisition rate is sufficiently high relative to 1/τ = k/γ, it is feasible to calibrate the trap using the bead position power spectrum, allowing comparisons to the previous results at low laser power. Power spectrum fitting with the blur-corrected and aliased expression (Eq. (A14)) at the lowest power yielded both a spring constant and friction factor that agree with the blur-corrected equipartition values to within 1%. Fits of the same data using the naïve expression (Eq. (A4), not corrected for exposure time or aliasing) provided slightly worse results, overestimating the spring constant by 3% and the friction factor by 7%. (See Appendix C for procedural details.) For an additional check that does not rely on fluctuations, a purely mechanical test was performed and compared with the corrected power spectrum fit. This test consisted of a bead drop experiment to determine the bead radius and friction factor, and a trap recoil experiment to determine the spring constant. The bead drop was performed by releasing a bead and recording its average velocity over a known distance. The trap recoil experiment was performed by measuring the exponential decay of the same bead as it returned to the trap center after deviation in one dimension. This mechanical test agrees with the corrected power spectrum fit to within 5% for the determination of both the spring constant and the friction factor.

DISCUSSION: PRACTICAL SUGGESTIONS FOR CALIBRATING AN OPTICAL TRAP
In this section, we present some practical techniques for measuring the spring constant k and diffusion coefficient D of a harmonically confined Brownian particle. We will assume that the temperature T is known. The approaches here are generic, and can be used even if the confining potential is not an optical trap (e.g. beads embedded in a gel, etc.) We will continue to treat the measured position as an unweighted time average of the true position over the If the diffusion coefficient D of the confined particle and the integration time W of the instrument is known, the correct spring constant k can be directly obtained from the measured variance var(X m ) by using equation 6. If the tracking error ε is significant, ε 2 should first be subtracted from the measured variance as in equation 9. While we cannot in general isolate k in this transcendental equation, it can easily be found numerically by utilizing a standard root-finding method. Alternatively, an approximate closed form solution for k is derived in Appendix D.

Determining k and D by varying W
Even if the data acquisition rate of the system is not fast enough to permit a blur-corrected power spectrum fit (as described in Appendix C), k and D can still be determined by measuring the variance at different shutter speeds and fitting to the blur-corrected variance function. This technique is demonstrated in the experimental results section, and yields accurate measurements provided the integration time is not too much larger than the trap relaxation time (α is not much larger than 1). Practically speaking, this is a useful technique, as the maximum shutter speed of a camera is often much faster than the maximum data acquisition speed (e.g. it is much easier to obtain a video camera with a 0.1 ms shutter speed than a camera with a frame rate of 10 kHz). Furthermore, this approach for quantifying the power spectrum from the blur is quite general, and could be used in other systems. As long as the form of the power spectrum is known, the model parameters could be determined by measuring the total variance over a suitable spectrum of shutter speeds.

Determining k and D by varying k
Other approaches are possible if the confinement of the particles can be varied in a controlled way, i.e. by varying the laser power of the optical trap. If the spring constant varies linearly with laser power, (which is typically true and was confirmed for our system in subsection (4.3)), the first observation is that the spring constant only needs to be measured at a single laser power, as it can be extrapolated to other laser powers. Typically calibration should be done at a low power, as this usually increases the accuracy of both the power spectrum fit and the blur correction technique.
Linearity between the spring constant and laser power can be further exploited to determine both k and D by measuring the variance of a trapped bead at different laser powers but with the same shutter speed. Such data can be fit to the blur model (equation 6, recalling that α = W Dk/(k B T )) by introducing an additional fitting parameter c that relates the laser power P to the spring constant, i.e. we make the substitution k = cP , and perform a non-linear fit to var(X m ) vs. power data in order to determine c and D. Equivalently, we can express the naïve spring constant k m = k B T /var(X m ) as a function of c and P , and perform a fit to k m vs. P data as shown in Fig. 3 of the experimental results subsection (4.3), where the viability of this method is demonstrated.

Design strategies for using the blur technique
When using these motion blur techniques to characterize the dynamics of confined particles, we reiterate that it is the shutter speed and not the data acquisition speed that limits the dynamic range of a measurement. Thus, even inexpensive cameras with fast shutter speeds can make dynamical measurements without requiring the investment of a fast video camera. Alternative methods for controlling the exposure time are the use of optical shutters or strobe lights.

CONCLUSIONS
We have experimentally verified a relationship between the measured variance of a harmonically confined particle and the integration time of the detection device. This yields a practical prescription for calibrating an optical trap that corrects and extends both the standard equipartition and power spectrum methods. By measuring the variance at different shutter speeds or different laser powers, the true spring constant can be determined by application of the motion blur correction function of Eq. (7). Additionally, this provides a new technique for determining the diffusion coefficient of a confined particle from time-averaged fluctuations.
The dramatic results from our experiment indicate that integration time of the detection device cannot be overlooked, especially with video detection. Furthermore, we have shown that motion blur need not be a detriment if it is well understood, as it provides useful information about the dynamics of the system being studied. In this appendix, a derivation of the motion blur correction function Eq. (7) is presented. This quantifies how the measured variance depends upon the spring constant k, the diffusion coefficient of the particle D, and the integration time of the imaging device W (notation is as introduced in section (2)). The derivation follows standard techniques (e.g. [15]) and is similar to calculations presented in references [7,8,9,16]. For completeness, we present the calculation in two different ways: (1) a frequency-space calculation that convolves the true particle trajectory with the appropriate moving-average filter, and (3) a real-space calculation using Green's functions. An expression for the modified power spectrum of the harmonically confined bead that accounts for the effects of filtering and aliasing is included in the frequency-space calculation (see subsection (2))

Frequency-space calculation
The measured trajectory of a particle in the presence of motion blur X m (t) can be calculated by convolving the true trajectory X(t) with a rectangular function, where H(t) is defined by: The integral is taken over the full range of values (i.e. t ′ is integrated from −∞ to +∞), which is our convention whenever limits are not explicitly written. The width of the rectangle W is simply the integration time as previously defined. This convolution acts as an ideal moving average filter in time, and is consistent with the integral expression for X m (t) given in Eq. (1). Taking the power spectrum of Eq. (A1) yields: Where the Fourier transform is denoted by a tilde, e.g.X(ω) = X(t) exp(iωt)[ . t], and ω is the frequency in radians/second. The theoretical power spectrum P (ω) is given by: where γ is the friction factor of the particle, and is related to the diffusion coefficient by the Einstein relation γ = k B T /D. This power spectrum has been well-described previously [11,12,16], and is derived in section (3). The power spectrum of the moving average filter can be expressed as a squared sinc function: Using Parseval's Theorem and integrating the power spectrum P (ω) yields the true variance of X(t), which is in agreement with the equipartition theorem. Similarly, we calculate the measured variance var(X m ) as a function of the exposure time W and the friction factor γ by integrating the power spectrum of the measured position (Eq. (A3)): var(X m ) = 1 2π where τ = γ/k, the trap relaxation time. Writing this formula in terms of the dimensionless exposure time, and the variance of the true bead position var(X) = k B T /k yields: var(X m ) = var(X)S(α) (A10) where S(α) is the motion blur correction function:

Blur-corrected filtered power spectrum
Often, trap calibration is performed by fitting the power spectrum of a confined particle. Here we provide a modification to the standard functional form P (ω) that accounts for both exposure time effects and aliasing. Combining expressions A3, A4 and A5, we can see the effect of exposure time on the measured power spectrum: Additionally, the effect of aliasing can be accounted for: where ω s is the angular sampling frequency (i.e. the data acquisition rate times 2π). Aliasing changes the shape of the power spectrum, so neglecting it when fitting can cause errors. The sum in Eq. (A14) can be calculated numerically and fit to experimental data. It is typically sufficient to calculate only the first few terms. It is important to note that aliasing does not affect our result for the measured variance, Eq. (A8). Aliasing shifts power into the wrong frequencies, but does not change the integral of the power. Hence, var(X m ) is unchanged. A detailed discussion of power spectrum calibration with an emphasis on photodiode detection systems is given in reference [14].

Real-space calculation
Since a Brownian particle follows a random trajectory X(t), the measured position X m is a random function of the true position of the particle at the start of the integration time, i.e.
where X(t | x 0 ) is the actual position of the bead at time t given that it is at position x 0 at time zero, and W is the integration time as defined previously. In other words, even with knowledge of the initial particle position, it is not possible to predict what the measured position will be. However, the distribution of X m is well-defined, and one can determine its moments. The variance of the measured position is given by Notice that to calculate the ensemble average . . . , we must average over both the random initial position X, and the measured position for a given initial position X m (x). For the harmonic potential U (x) = 1 2 kx 2 , X m (X) = 0 by symmetry, so the variance reduces to where ρ X (x 0 ) is the probability density of the initial position, and the integral is taken over all space (consistent with our previously stated convention). In equilibrium, ρ X (x 0 ) is simply the Boltzmann distribution given in Eq. (2). Using Eq. (A15), we express X m (x 0 ) 2 as the double integral In the second step, the ensemble average is brought into the integral, and the averaging condition t 2 > t 1 is added, which changes the limits of integration. The time-ordered auto-correlation function X(t 1 )X(t 2 ) t2>t1 can be calculated using the Green's function of the diffusion equation for a harmonic potential, ρ (x, t | x 0 , t 0 ). The Green's function represents the probability density for finding the particle at position x at time t given that it is at x 0 at time t 0 . It can be found by solving the diffusion equation with the initial conditions ρ(x, t 0 ) = δ(x − x 0 ). The solution to this problem is well-known [16,18] and is given by: where we have defined the dimensionless function: As before τ = k B T /(kD) = γ/k. Notice that this is simply a spreading Gaussian distribution with the mean given by the deterministic (non-Brownian) position of a particle connected to a spring in an overdamped environment, and with a variance that looks like free diffusion at short time scales (i.e. initially increasing as 2D(t − t 0 )), but exponentially approaching the equilibrium value of k B T /k on longer time scales. The time-ordered auto-correlation function can be written as follows: Putting in the Green's function of Eq. (A21) and evaluating the integrals gives the result: where V (t) and τ are as defined above. Carrying out the double time integral in Eq. (A18), followed by the integral over the initial position x 0 of Eq. (A17) we obtain the final result for the measured variance: This reproduces the result of the frequency-space calculation presented in Eq. (A8). The ideal power spectrum can be obtained from the position auto-correlation function of Eq. (A24). We determine the long-time limit of the auto-correlation function by letting t 1 ≫ τ , which yields the simplified equation: Next, by taking the Fourier transform of this equation with respect to (t 2 − t 1 ) we obtain the standard result of Eq. (A4).

APPENDIX B: HIGH-PASS FILTERING IN VARIANCE MEASUREMENTS
Calculation of the variance requires special attention, since low frequency noise or drift can inflate the variance dramatically, causing an underestimation of the spring constant. A high pass filter can be used to remove low frequency noise, but the use of any ideal filter lowers the variance by neglecting the contribution from the removed frequencies (note Eq. (A6)).
To reliably estimate the variance while accounting for low frequency drift, we first progressively high-pass filter the data over a range of increasing cut-off frequencies. A plot of measured variance vs. cutoff frequency (Fig. 4) clearly shows a linear trend at frequencies below the corner frequency (f c =k/2πγ). However, as the filtering frequencies approach zero, drift causes the measured variance to increase beyond its expected value. By applying a linear fit and extrapolating to the 0 Hz cutoff, we can reliably estimate the "drift-free" variance of bead position. Power spectrum calibrations were performed by fitting the one-sided power spectrum with Eq. (A14). The original 65536 data points taken at ∼1500 samples per second were blocked into 128 non-overlapping segments. The power spectrum of the blocks were calculated separately and averaged to produce the data in Fig. 5. This procedure is well described in the literature [12,14]. This data was fit with the blur-corrected and aliased model of Eq. (A14) and compared with the commonly used non-corrected power spectrum of Eq. (A4), both with and without aliasing. The quality of the fit to Eq. (A14) was further investigated by examining the fractional deviation in the power (the measured data divided by the model fit) as in reference [14]. A scatter plot and histogram of the fractional deviation is presented in Fig. 6. The histogram agrees well with a Gaussian distribution with a standard deviation of 1/ √ 128 (see reference [14] for a thorough discussion of power-spectrum fitting, including the expected scatter from unity).
Accounting for tracking error in the power spectrum fit is more difficult than in the equipartition case, requiring knowledge of the frequency dependence of the error. To investigate this in the current study, the power spectrum of a stationary bead was subtracted from the calibration power spectrum. We found that the fit parameters remained practically unchanged (within 2%), allowing us to neglect tracking error in our power spectrum fits at low power. It should be noted that in other situations (e.g. different bead size or power), modifications to the power spectrum due to tracking error could be significant.  [14]. Right: Histogram of the fractional deviation data overlaid with a Gaussian distribution with a standard deviation of 1/ √ 128 (solid red line).

APPENDIX D: APPROXIMATE ANALYTICAL EXPRESSION FOR k
When W is not significantly larger than the trap relaxation time, i.e. α = W/τ = W k/γ is not much larger than 1, an approximate version of equation 6 can be inverted to give a closed form solution for k. First, we use a Padé approximation to express the motion blur correction function as: S(α) ≈ 1 − 2α/15 + α 2 /60 1 + α/5 (D1) Substituting this expression into equation 6 yields a quadratic equation that is easily solved for k. This results in the following approximation for the true spring constant: The Padé approximation is good to within 3% for α < 3, which corresponds to a blur correction factor of S(3) . = 0.46. In other words, if the uncorrected equipartition method gives a spring constant which is within a factor of 2 of the true value, this approximation formula should be accurate to within 3%, as we have tested numerically.