1. Introduction

Optical traps have become an important tool in the areas of biology and biophysics since their invention [1]. The change in the momentum of light as it is diffracted and reflected through a dielectric particle results in a force that is proportional to the power of the beam. The relatively small magnitudes (~1–100 pN) and noninvasive nature of the forces that an optical trap can apply have made them especially useful for the manipulation of cells and motor proteins. Optical traps have been used to measure forces associated with the motor proteins myosin [2] and kinesin [3] as well as the compliance of DNA [4]. Typically biological specimen are tethered to a microsphere, which acts as a handle. The optical forces on these handles is described by Hooke’s law (F=kx), where F is the force, k is the trapping stiffness, and x is the displacement from the center of the trap. Many methods exist for measuring the particle’s displacement [5].

A few common methods for calibrating the trapping stiffness, k, are the drag or escape method, the equipartition theorem method, and the power spectral method [6]. The first method calibrates the stiffness by creating a fluid flow around the particle and observing the displacement. The stiffness can be calculated from the velocity, ν, and Stoke’s drag coefficient γ=6πηr (η is fluid viscosity and r is the particle diameter) according to k=γν/x. This method requires an actuator to move the particle with respect to the fluid, e.g., a moving stage, pump, or the trap itself. A second method uses the equipartition theorem in which the mean-squared displacement and thermal energy of the system is related to stiffness according to k=k_BT/〈x ²〉, where k_B is Boltzmann’s constant and T is the absolute temperature. This requires a calibrated position sensor to measure the mean-square displacement. The power spectral method curve fits a plot of the auto-spectrum, which is the Fourier transform of the autocorrelation function, of the displacement signal. The corner frequency of the auto-spectrum, Ω, is related to the stiffness according to Ω=k/γ [7]. By using the power spectral method to calibrate the stiffness, the equipartition method can then be used to calibrate the position senor, since the integral of the auto-spectrum is the mean-squared displacement.

The basis for the power spectral method lies in the modeling of an optically trapped particle as a spring-mass-damper system. The equation of motion is as follows: mẍ+γẋ+kx=kd̃, where m is the object’s mass, x is the position of the object with respect to the trapping center, and d̃ is the fluctuating Brownian disturbance, which, in this case, is white, zero mean, with a variance of ${\sigma }_d^2=\frac{2\gamma k_{B}T}{k^2}$ . In water, the system is highly over damped due to the dominance of viscous forces over inertial forces. In this case, the system can be simplified to a first order system according to: γẋ+kx=kd̃. A continuous-time frequency-domain representation of this system can be made by Laplace transforming this equation (assuming zero initial conditions). A frequency-domain transfer function is the ratio of the output of the system to the input. The transfer function for this system is: G(s)=X(s)/D̃(s)=k/(γs+k), where s is the Laplace variable. The auto-spectrum of this system is

where ω is the frequency in radians/second. The corner frequency, Ω, dictates the speed of which a trapped particle can be manipulated.

Every particle that is trapped is different, which results in different corner frequencies and, therefore, trapping stiffnesses. In addition, significant changes occur when switching between particles of different sizes or indices of refraction, and, even if the same particle is used, changing the laser power changes the trapping stiffness. Figure 1 illustrates this by showing a single particle trapped at multiple power levels as well as different particles trapped at the same power setting; the resulting corner frequencies vary dramatically. This is especially critical considering that force measurements use feedback control to regulate the forces applied to trapped particles, and changes in the system, such as the stiffness, alter the system dynamics and the performance of the force controller. Furthermore, controller design is a time consuming process that requires knowledge of the system for optimal controller design. Potentially changing plant dynamics motivate a method to automate the process of system identification and stiffness calibration.

 

Fig. 1. Examples of changing plant dynamics depending on trapping power and particle size. (a) A 1 µmpolystyrene sphere trapped at power levels varying from20 mWto 80 mW. (b) Multiple particle sizes (0.5,1,5 µm) and indices of refraction (polystyrene and silica) trapped with the same power.

Download Full Size | PDF

Here we present an adaptive system identification method for the modeling of optical trap dynamics. The identified system dynamics can be used to locate the system’s corner frequency to calibrate trapping stiffness using the power spectral method. The method is based on an adaptive least-mean-square (LMS) algorithm, which replicates a discrete time model of the system dynamics using a gradient descent method. The identified model is a high order finite impulse response (FIR) filter. The filter order is then reduced using balanced model reduction to determine the corner frequency. This method has an advantage over other techniques in that it is quick, does not explicitly require operator interaction, and can be acquired in real time. It is also a necessary step for the implementation of an adaptive controller which can automate the controller design process and update the controller for changes in the trapping dynamics (e.g., power fluctuations) and for particles of different sizes and refractive indices.

2. Adaptive least-mean-square algorithm

The system identification method uses the LMS algorithm [8] to adapt the tap weights of an FIR filter to mimic the system dynamics. The FIR filter is a discrete time filter in which the output is the sum of the incrementally time step delayed inputs where each time step is multiplied by a tap weight (see Fig. 2). This is also known as a transversal or tapped delay line filter. Discrete time transfer functions are in the z-domain where z ^-1 represents a single time step delay while z ^-m represents an mth time step delay. For an mth order filter, the output of the filter, y, at a discrete time step n is defined as: y(n)=w′(n)u(n), where w(n)=[w ₀(n),w ₁(n), …,w_m(n)]′ is the vector of tap weights, and u(n)=[u(n),u(n-1), …, u(n-m)]′ is the vector of tap inputs.

 

Fig. 2. Transversal or tapped delay line filter where u(n) is the input, y(n) is the output, z ^-1 is a time step delay, and w is the tap weights.

Download Full Size | PDF

These filters are easily produced in Matlab’s Simulink environment using integer delay blocks.

The FIR filter’s tap delay weights are updated using the LMS algorithm. This is the simplest stochastic gradient-based algorithm since it does not require measurement of correlation functions or matrix inversion, unlike other adaptive filtering algorithms [9], though the approach presented here would be possible with other adaptive algorithms. The LMS algorithm is designed to minimize a user defined error function, e, which will be defined later. The LMS algorithm has a mean-squared error cost function: $J=\frac{1}{2}ℰ\left\{e{\left(n\right)}^2\right\}$ , where ℰ{·} denotes the statistical expectation operator. The cost function is quadratic in w, creating a surface on which a global minimum can be found using a gradient descent method. To minimize the cost function, the tap weights are adjusted in the negative gradient direction of the cost function. The minimum is achieved when the gradient is zero, ∇J=-ℰ{u(n)e(n)}=0; that is, when the filter error and tap input are orthogonal. Since determining the expected value of the gradient is difficult, an instantaneous estimate of the gradient is used in the LMS algorithm. The resultant weight update equation is

where µ is the step size.

The LMS algorithm is used to identify the system dynamics in a series arrangement. The relative position of the trapped particle is the FIR filter input, or to put this in the context of the notation used previously, u(n)=x(n). The resulting FIR filter is an inverse of the system’s transfer function. In the z-domain, the resulting transfer function model of the system dynamics is

where the zeros (or numerator) of the FIR filter end up being the poles (or denominator) of the system. There is no guarantee that the resulting model Ĝ(z) is stable since the zeros of the FIR filter could be anywhere in the z-plane; but, in this situation this is not important since the objective is to match the frequency response of the system dynamics with the identified model, and ultimately determine the stiffness. It is important to note that this configuration only finds the corner frequency of the system and not the magnitude of the response. The magnitude is fit using the equipartition theorem.

In order to prevent all of the weights from being made arbitrarily small to reduce the error, the zeroth weight is constrained to unity, w ₀=1, and the remaining m weights are allowed to adapt. This is done by defining the error function as: e(n)=y(n)+x(n). Intuitively, this is can be thought of as comparing a Laplace transform to a z transform, which for the system presented here is:

The first value of the denominator of the z transform is unity and the second value is dependent on the corner frequency, a, and the sample period, T. By allowing the rest of the weights to adapt, the discrete time transfer function mimics the actual continuous time system.

The system dynamics are a first order system. If a single tap delay is used to find the corner frequency of the system, noise in the system leads to an underestimation of the corner frequency. A third order filter is the smallest order necessary in order to match these dynamics, as determined experimentally. In order to determine the corner frequency of the higher order FIR model, the identified model in Eqn. 3 is reduced. In this work, a balanced model reduction is used [10] to capture the important characteristics of the system’s frequency response by retaining only those states that contribute to signal throughput as measured by the Hankel singular values.

3. Experimental methods and results

The experiments were conducted using a custom built optical trap based on an inverted microscope (Zeiss Axiovert 200) and a Nd:YAG laser at 1064 nm (Coherent Compass 1064-500). The laser enters the microscope though the epi-fluorescence port and is introduced to the microscope’s optical path using a dichroic mirror located in the microscope’s filter cube turret. The objective is a Zeiss Plan-Apochromat 63x/1.4 NA. Displacements of the trapped particle are measured using the forward scattered laser,which is collected by a high NA condenser (Zeiss Achromatic-aplanatic condenser 1.4 NA). The laser beam is separated from the illumination with a dichroic mirror located above the condenser and directed through a collimating lens and onto a quadrant photodiode (QPD) [11]. The high NA condenser allows the collection of most of the scattered wave resulting in a high precision displacement sensor.

The displacement data from the QPD was processed using a dSPACE D/A board and software. A Simulink model was built in Matlab and compiled onto the dSPACE D/A board to produce a normalized displacement signal as well as to update the filter weights at each time step according to Eqn. 2. The filter output was calculated by taking the dot product of the tap delayed displacement signal and the updated filter weights. The error was calculated by summing the displacement signal with the output of the FIR filter. The system model was reduced using a balanced model reduction after the weights were allowed to converge. (An example of the Simulink model is available upon request from the authors.) The particle displacement was sampled at 12.8 kHz.

In order to test the capabilities of the adaptive system identification method, multiple particles were trapped in varying power levels. Four different particles were used to test the algorithm: 0.5 µm, 1 µm, 5 µm polystyrene and 1 µm silica microspheres. The multiple sizes of polystyrene tests the system with different sized particles while the silica microspheres allows the comparison of different refractive indexes. The particles also give a wide range of corner frequencies for which to compare.

The experimental procedure was to trap each particle 5 µm above the coverslip (10 µm for the 5 µm spheres) in order to minimize the wall effects and to gain an accurate representation of the particle. Adaptive system identification used three tapped delays (m=3). Multiple delays, from 1 to 10 were tested, with 3 being the smallest size that accurately captured the system dynamics. The step-size, µ, was set to 0.08. The adaptive system identification program was then executed and allowed 1 min for the weights to converge. Particles with higher corner frequencies converged in less time, but 1 min was used for consistency. The movement of the particle was tracked on the QPD and recorded, as well as the final value of the weight vector. The resulting identified plant transfer function was then reduced to a first order system and plotted along with the auto-spectrum of the data. The magnitude of the reduced order model was fit using the equipartition theorem. As can be seen in Fig. 3, the modeled system accurately tracks the corner frequency as well as the magnitude. The corner frequencies vary from 17 Hz for the 5 µm to 4213 Hz for the 0.5 µm microspheres.

 

Fig. 3. Examples of changing plant dynamics depending on trapping power and particle size. (a) 1 µm silica sphere trapped in 20 mW, Ω=282 Hz. (b) 0.5 µm polystyrene sphere trapped in 20 mW, Ω=740 Hz

Download Full Size | PDF

Table 1. Resultant corner frequencies and stiffnesses for various particle diameters.

View Table

Adaptive system identification successfully identified the corner frequency of the spherical particles confined in an optical trap. The corner frequency was then used to compute the trapping stiffness according to k=Ωγ. The results for various trapping conditions are presented in Table 1. These results are comparable to those achieved previously [12]. While this presents an automated method for characterizing the stiffness of the optical trap, the real benefit is the ability to monitor the stiffness of the trap in real time.

4. Conclusions

We have successfully demonstrated an adaptive LMS approach for determining the trapping stiffness of an optical trap in real time. The method is capable of characterizing the corner frequency, which is the important system characteristic for determining trap stiffness. Multiple particle sizes with varying indices of refraction were experimentally verified under different power levels, all of which result in different system dynamics, and all with results that validate our technique. This work has the potential to improve the sensitivity of stiffness measurements as well as act as a starting point for the creation of an adaptive control of optical traps. The work presented here will ultimately lead to the design of optical trap controllers through adaptive techniques for servo control of trapped particles as well as force controller.