## 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 ${\mathrm{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.

## 2.

## 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 well-known diffusion expression^{63}

^{64}

^{,}

^{65}where ${k}_{B}$ is the Boltzmann constant, $T$ is the temperature of the fluid medium, $\xi $ is the Stokes (viscous) drag coefficient of the trapped particle ($\xi =6\pi \mu R$ for a sphere with radius $R$ suspended in a fluid with viscosity $\mu $), and $\kappa $ is the spring constant of the optical trap. Note that $D={k}_{B}T/\xi $.

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 $\kappa $ and $\xi $. 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

## 3.1.

### 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.

## 3.2.

### Software Architecture

## 3.2.1.

#### 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 data-loss 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.

## 3.2.2.

#### 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.

## 3.2.3.

#### 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}\u2013\tau +1$ overlapping time intervals of duration $\tau $, the position autocorrelation function is defined as^{63}

## (3)

$$[\mathrm{\Delta}{r}^{2}(\tau )]=\frac{1}{{N}_{A}}\sum _{t=0}^{{N}_{T}-\tau}{[r(t+\tau )-r(t)]}^{2}.$$The *MSD(PosArray)* function iteratively calculates the MSD for each lag-time, $\tau $. *PosArray* is an array where each row is a position vector for a specific time. The function 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$.

## 3.2.4.

#### 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.

## 3.2.5.

#### 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, $\alpha $ (normalized volts per meter) and $\beta $ (seconds per desired time unit). The new model can be expressed as

## (4)

$$[\mathrm{\Delta}{r}^{2}(\tau )]=\frac{2{k}_{B}T{\alpha}^{2}}{\kappa}(1-{e}^{-\kappa \beta \tau /\xi})\equiv a(1-{e}^{-b\tau}).$$Thus, our model can be expressed in terms of $a$ and $b$

## (5)

$$a=\frac{2{k}_{B}T{\alpha}^{2}}{\kappa};\phantom{\rule[-0.0ex]{1em}{0.0ex}}b=\frac{\kappa \beta}{\xi}.$$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.

## 3.2.6.

#### 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.

## 3.2.7.

#### 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.

## 3.2.8.

#### 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 $\alpha $ is found by taking the ratio of the known diffusion constant with the fitted counterpart. We solved for $\alpha $ given the fitting constants $a$ and $b$. The software then converts all the spring constants to physical units [$\mathrm{N}/\mathrm{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.

## 4.

## Experiments and Results

## 4.1.

### 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

## (6)

$$\stackrel{\xb7}{\overrightarrow{r}}(t)=-\frac{1}{\xi}\overrightarrow{k}\xb7\overrightarrow{r}(t)+\sqrt{2D}\overrightarrow{W}(t),$$## (7)

$${\overrightarrow{r}}_{i}=(1-\frac{\mathrm{\Delta}t}{\xi}\overrightarrow{k})\xb7{\overrightarrow{r}}_{i-1}+\sqrt{2D\mathrm{\Delta}t}\xb7{\overrightarrow{w}}_{i},$$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.

## 4.2.

### 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\pm 0.0002\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{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.

## 4.3.

### 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.

## 4.4.

### 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 ($\mu =1.416\mathrm{P}$), but because the MSD fit curve allows separate fitting of trap stiffness $\kappa $ and viscous drag $\xi $, 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 $\lambda =1064\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mathrm{nm}$, ${n}_{\mathrm{gly}}=1.46$, while ${n}_{\mathrm{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 $\rho $ and viscosity $\mu $) 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 $\gamma =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.

## 5.

## Conclusion and Future Work

In conclusion, we present a data acquisition and analysis toolkit based on particle tracking that allows near real-time 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-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.

## Acknowledgments

Funding for this work, provided by the National Institutes of Health under Grant R15DK092716, is gratefully acknowledged.

## References

*orpk*mice,” Am. J. Physiol. 289(5), F978–F988 (2005).AJPHAP0002-9513 http://dx.doi.org/10.1152/ajprenal.00260.2004 Google Scholar

## Biography

**Joseph Glaser** received his BS degree in physics from Cleveland State University in 2014 and is currently enrolled in the physics graduate program at Drexel University.

**David Hoeprich** received his MS degree in medical physics from Cleveland State University in 2012. He is currently employed as a board-certified medical physicist at Applied Medical Physics In Radiology, Inc., located in Kirtland Hills, Ohio.

**Andrew Resnick** is an assistant professor at Cleveland State University. He received his BS degree in physics from Rensselaer Polytechnic Institute in 1991 and his PhD degree in physics from the University of Alabama in Huntsville in 1996. His current research interests include cellular mechanotransduction and microscopy. He is a member of the American Physical Society, American Physiological Society, and Optical Society of America.