Near real-time measurement of forces applied by an optical trap to a rigid cylindrical object

Abstract. An automated data acquisition and processing system is established to measure the force applied by an optical trap to an object of unknown composition in real time. Optical traps have been in use for the past 40 years to manipulate microscopic particles, but the magnitude of applied force is often unknown and requires extensive instrument characterization. Measuring or calculating the force applied by an optical trap to nonspherical particles presents additional difficulties which are also overcome with our system. Extensive experiments and measurements using well-characterized objects were performed to verify the system performance.

Near real-time measurement of forces applied by an optical trap to a rigid cylindrical object 1 Introduction Our laboratory is interested in measuring the mechanical properties of a biological object, the primary cilium. 1 Optical traps provide a localized noncontact method to apply a well-controlled force while simultaneously permitting the observation of the deformation (bending) of this organelle. Here, we present our progress toward this goal by demonstrating our ability to accurately measure the spring constant of our optical trap applied to a cylindrical object, in this case an Escherichia coli bacterium.
Although use of optical traps 2 to apply forces to spherical and near-spherical objects is well understood in the context of Mie scattering and its generalizations, [3][4][5][6][7][8][9] trapping cylindrical objects such as bacteria and certain virus particles has largely consisted of qualitative experiments, [10][11][12] and use of scattering models to quantify the optical trapping of slender, cylindrical objects is greatly complicated by geometry. [13][14][15][16][17] Use of approximate methods 18 including the discrete-dipole approximation 19,20 the T-matrix, 21 or the finite-difference time-domain 22 to calculate the applied force is one common approach, but accurate calculation of the applied force (equivalently, the trap stiffness) requires detailed trapping beam specifications; knowledge of the refractive index of both the trapped object and the solvent; viscosity of solvent; and size, shape, and composition of the trapped particle, and typically assumes that the trapped object is homogeneous. For our experimental case of interest, none of this information is readily available and may require additional complex measurements to obtain. Similarly, optical trap calibration methods involving measurements of the power spectrum 23,24 or scattering efficiencies [25][26][27] have only been experimentally verified using spherical objects. Consequently, alternative methods of calculating the applied force based on particle tracking have been developed. [28][29][30] The primary cilium 31 is a singular protruding hair-like organelle possessed by most mammalian differentiated cells, genetically related to flagella, but in contrast to motile cilia present on some specialized tissue (e.g., airway epithelia), it does not actively move. Cilia are slender, protruding cylinders, approximately 0.2 μm in diameter, ranging in length from several to tens of microns and are anchored at one end to the centrosome within the cell body. A large and still growing body of evidence has demonstrated that the primary cilium coordinates an ever-lengthening list of cell signaling pathways including Hedgehog, Wnt, platelet-derived growth factor receptor alpha, and the polycystins PC1 and PC2. [32][33][34][35][36][37][38][39][40][41][42][43][44][45][46][47] The primary cilium has recently emerged as an organizing center for a wide range of cellular functions including environmental sensing of light, 48 odorants 49 and osmolarity, 50 mechanosensation, 51,52 cell development, 32 migration 53 and differentiation, 54 and planar cell polarity. 46 The mechanotransduction response has been shown to significantly modify many cellular functions and yet is currently the least-understood ciliary function due to a lack of information about the mechanical properties of the sensor. In contrast to measurements that apply a known disturbance with fluid flow, 55,56 optical trapping of cilia and flagella could provide improved information about the mechanical properties of these important organelles. Specifically, it is possible to measure the bending modulus 57,58 and generated force, 59,60 both of which modulate the biological role of motile cilia and flagella (mucociliary transport and microorganism propulsion) and sensory cilia (flow sensing). 51,52 Finally, the capability of near-simultaneous trapping and observation of a biological response (for example, intracellular Ca 2þ ratiometric imaging) 1,37 could provide improved understanding of the biological role of this organelle.
In the present study, we have constructed a calibrated optical trap apparatus that provides near real-time measurements of the transverse spring constant without requiring knowledge of optical properties of the trapped object, suspending solvent, and trapping beam geometry. Because the primary cilium is anchored at one end to the cell and thus a more complicated system than a free particle, we demonstrate here a first step-measuring the trapping force applied to an E. coli bacterium, a rod-shaped bacterium 1 μm in diameter and 3 μm long.

Principles of Optical Tweezer Operation
As constructed, 61,62 optical tweezers create a single-beam three-dimensional (3-D) potential well due to a spatial gradient of the electric field, created by focusing a laser beam with a high-numerical-aperture microscope objective. Objects possessing a refractive index greater than the surrounding fluid medium experience a spring-like restoring force as they move away from the trap center. The trapped particles thus execute a modified form of Brownian motion due to the confining potential well. For a particle undergoing free one-dimensional Brownian motion, the mean-squared displacement (MSD) relationship is given by the wellknown diffusion expression 63 where the D is the diffusion constant for the particle, the angle brackets indicate an average over all time, and τ is the lag-time. A particle confined within a Gaussian potential well (corresponding to the optical trapping by a focused Gaussian beam) presents a modified MSD relationship 64,65 where k B is the Boltzmann constant, T is the temperature of the fluid medium, ξ is the Stokes (viscous) drag coefficient of the trapped particle (ξ ¼ 6πμR for a sphere with radius R suspended in a fluid with viscosity μ), and κ is the spring constant of the optical trap. Note that D ¼ k B T∕ξ. While the spring constant and Stokes drag coefficient can be calculated analytically for the idealized case of a homogeneous spherical particle, a homogeneous suspending fluid, and aberration-free Gaussian beams, in general it is not possible to analytically calculate the spring constant of nonspherical particles of unknown composition trapped within an aberrated beam. Because our eventual goal is to calculate the force applied to a slender cylinder, we developed an algorithm based on particle tracking and fitted the calculated MSD curve using two free parameters corresponding to κ and ξ. We first track the position of a trapped object with a quadrant photodiode (QPD), calculate the MSD of the particle's trajectory, and then use a fitting algorithm in which the spring constant and viscous drag are free parameters to calculate the stiffness of the optical trap. That is, measuring the trajectory of a trapped particle allows the calculation of the gradient forces generated by the optical trap without a priori knowledge of the physical properties of the trapped object.
3 System Hardware Configuration and Software Architecture

Hardware Configuration
The laser tweezer apparatus shown in Fig. 1 is largely the same as discussed previously 61 and operates as a single-beam 3-D trap. The source was a 0.5-W diode-pumped Nd:YAG continuous-wave single-mode laser. The optical tweezer was aligned to the optical axis of the microscope (Leica DM6000B, Wetzlar, Germany) using a 5-degree-of-freedom (x-axis, y-axis, z-axis, pitch, and yaw) mount. The laser beam was expanded to fill the back aperture of the objective lens. The objective lens used was a Leica 63X NA 0.9 U-V-I HCX long working distance plan apochromat dipping objective. The tweezer couples into the microscope through a lateral port, providing direct optical access to the fluorescence turret. A side-looking 1064-nm dichroic mirror (Chroma Technology, Bellows Falls, Vermont) mounted within the fluorescence turret provides the ability to perform the normal (visible) transillumination microscope viewing while the tweezers are operating. A KG-1 IR cutoff filter (Newport, Irvine, California) inserted above the fluorescence turret prevents the observation of the tweezer spot by a camera during operation. The trapping beam applies a force to objects by moving the microscope sample stage (Prior Proscan II Motorized H101/2 Stage); the trapping beam does not move.
To record the movement of a trapped object, a QPD (first sensor model QP50-6-SD2) records the location of the trapping beam centroid. The QPD is located downstream from the condenser lens. Use of a dichroic mirror (Qoptic Microbench, Fairport, New York) and laser line filter (1064 nm, Edmund Optics, Barrington, New Jersey) ensures that only the Nd: YAG light is incident on the QPD. Intensities were converted to voltage readouts within the QPD housing and then transferred to a National Instruments data acquisition device. This device recorded voltages at a user-specified frequency and exported a digital record of particle positions to the host computer to be saved and analyzed later.

Tracking protocol
We created a virtual instrument (VI) environment in LabView 2010 with a front panel that allows for customization of sampling frequency and duration of data acquisition. Our tracking VI recorded voltage outputs from the QPD, normalized the data via a reference voltage to remove the circuit noise, and displayed the particle information in real time on the front panel. Figure 2 shows the flow of data through the VI, and Fig. 3 presents an example of the VI's voltage/position readout.
The recorded datasets are then converted and stored as .TDMS binary files, which are used for their low dataloss rate at high sampling frequency and low memory footprint. This file can then be directly imported into our MATLAB analysis algorithm.
Stock data analysis methods are not well suited for the large data files produced by high sampling rates. Therefore, an algorithm was developed to manage the large sets of acquired experimental data. All the analysis procedures are called using the function QPDanalysisBulk (TDMSFolder,TDMSfileID,f). This function has user input variables for the root directory of the .TDMS file to be analyzed and the filename of said .TDMS.

TDMS conversion to .MAT file type
In order to begin the data analysis, our algorithm first converts the data from the .TDMS file into a .MAT file type for MATLAB to quickly perform the rest of the data analysis. The called function, convertTDMS(true, TDMSfileID), has been developed within a community of developers since late 2009. The function can be found via the MATLAB Central File Exchange and is offered open source to those who wish to make their own additions to the development of the code. In our toolbox, we utilized version 9 of the software and made no major changes to the function.

Calculate the mean-square-displacement of trajectory
With our time-series voltage data loaded into MATLAB, the algorithm proceeds to calculate the MSD for the trajectory of the trapped particle. For any given trajectory rðtÞ discretized into N T þ 1 time steps consisting of N A ¼ N T -τ þ 1 overlapping time intervals of duration τ, the position autocorrelation function is defined as 63 The MSD(PosArray) function iteratively calculates the MSD for each lag-time, τ. PosArray is an array where each row is a position vector for a specific time.  utilizes the fact that MATLAB can efficiently vectorize calculations on an entire array rather than computing them element-by-element. For example, the deltaCoords variable, which is the difference in positions between two row-vectors separated by dt, is computed by deltaCoords = PosArray[(1 +dt:end),(1:2)]-PosArray[(1:end-dt),(1:2)]. This allows us to use the sum(deltaCoords^2, 2) function to sum each of the squares of the deltaCoord variable's row elements and to finally take the average over all times for each dt.

Remove linear drift of MSD
The MSD curve can be split into roughly two sections. The initial (short-time) section of the curve represents an inertial timescale, whereas the latter section represents a long time limit corresponding to the restoring force of the optical trap. However, the long time limit is often not reached and as a result, the calculated MSD contains a linear drift. This experimental artifact is removed through the use of the function RemoveLinarDrift(MSDArray,maxlagtime,freq). The function first applies a linear regression to the long-time calculated MSD data. The linear regression is then subtracted from the experimentally acquired MSD data. The corrected MSD data allow the fitting of the idealized curve equation to the experimental data. An example of this curve can be seen in Fig. 4, which shows the removal of linear drift and resultant, corrected MSD.

Fit primary dataset to idealized relationship and store fitting parameters
Once the experimental data have been corrected with the linear drift of the long-term limit removed, the algorithm fits the data to the idealized MSD function. Since our data have an arbitrary length and time units, we must introduce two conversion constants, α (normalized volts per meter) and β (seconds per desired time unit). The new model can be expressed as Thus, our model can be expressed in terms of a and b To enhance our fitting routine's accuracy, we first fit the long-time portion of our data with a linear model to obtain the value of a. Next, we fit our data using a robust nonlinear least squares fit, supported by the MATLAB Curve Fitting toolbox. During this process, the stored value of a is substituted into the model so that b is the only free parameter. A graph of this fit is then exported as a MATLAB .FIG file and a high-resolution .PNG to the root folder. Figure 5 details an example of the quality of the fit achieved by the algorithm. Utilizing the relationships of the constants a and b, the software then solves for the spring and diffusion constants for the dataset. Finally, these values are stored in their own array to be appended with more values down the line.

Subdivide primary dataset into data blocks
During QPD data acquisition, it is sometimes observed that a second object either passes through the trapping beam or is even trapped as well, resulting in either two trapped objects or the initial object being ejected from the trap. This results in QPD data containing position jumps not characteristic to confined Brownian motion, requiring sections of the QPD data to be excised. To ease the process in manually excluding these datasets, the main data series is broken up into data blocks consisting of position data over 1-s time intervals.

Repeat analysis on data blocks
With the primary dataset subdivided into smaller data blocks consisting of 50 to 100 k data points, our analysis process is then iteratively applied to every data block. By repeating this analysis, we can produce results for the linear correction and fitting process as well as return a two-dimensional trajectory plot over the time interval. This trajectory plot will enable easy detection for which datasets should be excluded from analysis due to particle-particle interactions within the trap.

Return fitting parameters, graphs, and raw data for user
After analysis is performed on all data blocks, the software finally deals with the problem of length-scale conversion. As detailed elsewhere, 23 for an object with a known diffusion constant, the conversion constant α is found by taking the ratio of the known diffusion constant with the fitted counterpart. We solved for α given the fitting constants a and b. The software then converts all the spring constants to physical units [N∕m]. Finally, the algorithm stores all the retuned constant arrays, graphs, and raw time-series data into a folder hierarchy within the root directory. This allows for a simple organization scheme for the user's convenience. The total data processing time varies with the data acquisition time, but as a guideline, 10 s of QPD data is fully analyzed in 5 min using an Intel Core 2 Quad Q9400 CPU running at 2.66 GHz with 3.5 GB of RAM.

Algorithm Calibration
Simulation of a trapped particle executing confined Brownian motion was used as an internal consistency check of our data analysis algorithm. The Langevin equation for a Brownian particle moving within a 3-D harmonic potential (neglecting the inertial term) is given as whererðtÞ is the 3-D coordinates of the particle at some time t,r · ðtÞ is the 3-D velocity vector,k is the strength of the trap along each of the axes, andWðtÞ is the white noise characteristic of Brownian motion. This differential equation can be solved for the position and gives the finite difference equatioñ wherer i ¼ ½x i ; y i ; z i represents the position of the particle at time t i , andw i ¼ ½wx i ; wy i ; wz i is a vector of Gaussian random numbers with zero mean and a variance of 1. For our simulation,w i is generated bỹ whereR is a 11 × 3 matrix of random numbers that are uniformly distributed in (0,1) such that the sum producesw i . The MATLAB code can be found freely at the link provided in Ref. 66 and can simulate 3-D trajectories for any number of particles given a sampling rate, time span, trap stiffness constant, and particles' "radii." Following Refs. 64 and 65, we independently generated a trajectory for a spherical particle in an optical trap using our system specifications and analyzed the simulated data with our toolkit. Our toolkit calculated the MSD fitting parameters for several sets of simulation data with varying sampling frequencies, time spans, trap strengths, and initial positions.
The returned constants fit the provided values within a range of 5% from the expected values.

Instrument Calibration
We calibrated the optical trap fitting parameters using homogeneous spherical particles, and we will show that once calibrated, the fitting algorithm correctly returns the spring constant for (1) spherical particles trapped near a glass surface, (2) spherical particles immersed in glycerol, and (3) cylindrical objects (E. coli bacteria). Each test was chosen to verify the fitting algorithm using uncertainties in viscous drag (tests 1 and 2), changes in trap beam geometry (test 2), and changes in object geometry and optical properties (test 3).
The optical trap was calibrated using 0.5-, 1-, 1.98-, 2.5-, and 5.46-μm diameter polystyrene microspheres (Bangs Laboratories, Fishers, Indiana), each stock solution diluted (by volume) 1∶10 9 in phosphate-buffered saline (PBS, Cellgro, Pittsburgh, Pennsylvania) to prevent the aggregation. Viscosity of PBS was measured with a Cannon-Fenske Routine viscometer (Induchem Lab Glass, Roselle, New Jersey) and pycnometer (Cole Parmer, Vernon Hills, Illinois), resulting in a value of 0.8816 AE 0.0002 cP at 25°C . The dilute suspension of (monodisperse) microspheres was placed within a hanging drop-type microscope slide, the laser trap was turned on, and the microscope stage was manipulated to bring a microsphere into the trap far from any solid surface using a camera for visual confirmation that a single microsphere was trapped. The height of the trapped particle relative to the microscope slide was recorded using the microscope's internal z-axis encoder readout. During lateral movements of the stage, the distance between trapped particle and microscope slide remained constant.
Comparing QPD output to Stokes drag measurements, shown in Fig. 7, allowed us to obtain a calibrated relation between QPD data and applied force for spherical objects, because the only variable is the particle diameter. The QPD recorded the positions of spheres before the stage was moved. After the QPD data was acquired, the stage was moved laterally at a known speed from its home position to second position approximately 2 mm away and then back to its home position. The test was performed at progressively higher stage speeds until the trapped particle fell out of the laser trap due to viscous drag. The stage speed at which the particle fell out of the laser trap was used to determine the drag force from the PBS solution on the microsphere using Stokes' law.

Validation of Trap Strength Algorithm
As internal consistency checks, QPD data were collected for 60 s at sampling rates of 10, 50, and 100 kHz to determine potential effects from sampling rate. Other datasets were acquired for spot centroids located at various positions on the QPD to determine effects from misalignment. The aforementioned analysis was performed on each dataset, and the calculated spring constant was unchanged, indicating our system is quite robust. Next, a single large dataset was divided into smaller sets of 20,000 to 30,000 data points. The same analysis was performed on these subsets of data. The resulting spring constant of the entire set was compared with the spring constant results of the subsets of data, and the algorithm returned consistent spring constants.
Taken together, our internal consistency checks demonstrate internal validation of our trap strength algorithm.

Validation of Trap Strength Calculation
Three experimental checks of the fitting algorithm were performed. First, we trapped microspheres near the glass surface as an experimental test of Faxen's law. 67 Next, we trapped microspheres suspended in glycerol rather than PBS. Finally, we trapped an E. coli bacterium, a rod-shaped object for which a closed-form solution to Stokes drag can be obtained, and for all tests, we compared the calibrated fitting algorithm output to the (maximum) viscous drag force.
The first validation step was performed by comparing the QPD data to Stokes drag data for particles trapped at various heights above the slide (Fig. 6). This check was performed because while the fluid drag varies with height, approximately given by Faxen's law, 68 the trapping force does not (the small amount of light back-reflected off the fluid-glass interface can be neglected). The calculated spring constant did not show a variation in height, as expected.
Next, we trapped a 1.98-μm-diameter microsphere in glycerol (Fisher Scientific, Pittsburgh, Pennsylvania). Glycerol is approximately 1000 times as viscous as PBS (μ ¼ 1.416P), but because the MSD fit curve allows separate fitting of trap stiffness κ and viscous drag ξ, the trap force can be distinguished from viscous drag. As seen in Fig. 7, the calibrated algorithm correctly measured the spring constant. Thus, we demonstrated that our calibration curve is not sensitive to variations in solvent viscosity.
Additionally, because the objective lens is uncorrected for glycerol (at λ ¼ 1064 nm, n gly ¼ 1.46, while n PBS ¼ 1.33), the objective lens introduced wavefront aberrations (primarily spherical aberration) to the trapping beam. Thus, this measurement confirmed the results of test 1 and also demonstrated that our fitting algorithm does not require detailed information about the trapping beam wavefront.
We next optically trapped a cylindrical object, a single E. coli bacterium. E. coli is a rod-shaped bacterium 1 μm in diameter and 3 μm in length. E. coli normally swim using a flagellum, so we first killed the bacteria with a brief exposure to intense ultraviolet light to avoid experimental artifacts. Once trapped, the bacterium oriented along the optical axis, and both QPD and Stokes drag measurements were performed. Critically important, the bacterium does not deform during the Stokes drag measurement. Modeling the bacterium as a cylinder of length L and radius R capped with hemispherical ends, the viscous drag at velocity U (fluid density ρ and viscosity μ) is given by The first term results from the two hemispherical caps, and the second term is the drag force per unit length for a cylindrical object 51 and contains Euler's constant γ ¼ 0.577. Again, as seen in Fig. 7, the calibrated QPD algorithm correctly computed the spring constant, demonstrating that our algorithm results are insensitive to object shape.

Conclusion and Future Work
In conclusion, we present a data acquisition and analysis toolkit based on particle tracking that allows near realtime measurements of the force applied by an optical trap to an object of unknown shape and composition. In addition, this analysis toolkit allows for a convenient calibration method for optical traps via the use of standard microspheres. This calibration step accurately calculates length-  Fig. 7 This graph shows the experimental results for the QPD algorithm calibration using microspheres (N ¼ 5 for each data point) overlaid on results of our algorithm output for N ¼ 5 microsphere suspended in glycerol and our algorithm output for a trapped (N ¼ 5) Escherichia coli bacterium. The error bars correspond to measured uncertainty in (vertical) Stokes drag and (horizontal) results of our algorithm.
scale conversions, via diffusion constant ratios, as well as the effective strength of the trapping potential. Future work will be focused on applying the optical trap to a primary cilium, in an effort to measure the mechanical properties of this important organelle.