Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Appl Phys B. Author manuscript; available in PMC 2011 April 1.
Published in final edited form as:
Appl Phys B. 2010 April 1; 99(1-2): 23–30.
doi:  10.1007/s00340-009-3843-y
PMCID: PMC2850131

Optimal strategy for trapping single fluorescent molecules in solution using the ABEL trap


Trapping of 10-nm-sized single fluorescent bio-molecules in solution has been achieved using high-speed position sensing and electrokinetic feedback forces in the Anti-Brownian ELectrokinetic (ABEL) trap. The high diffusion coefficient of small objects in solution requires very fast, real-time sensing of position, and this has been previously achieved using a simple rotating beam, but improved strategies are needed for the smallest objects, such as single nanometer-sized fluorescent molecules. At the same time, single molecules are limited in photon emission rate and total number of photons, so each emitted photon must be used as efficiently as possible. We describe a new controller design for the ABEL trap which features fast, knight’s tour scanning of an excitation beam on a 2D square lattice and a Kalman filter-based estimator for optimal position sensing. This strategy leads directly to a maximum-likelihood-based method to extract the diffusion coefficient of the object held in the trap. The effectiveness of the algorithms are demonstrated and compared to the simple rotating beam design through Monte Carlo simulations. Our new approach yields tighter trapping and a much improved ability to extract diffusion coefficients.

1 Introduction

Recently, researchers have made significant progress in feedback control and tracking of single Brownian objects in aqueous solution using optical microscopy, typically fluorescence [1]. Observing a single molecule in the liquid phase for an extended amount of time provides the capability to measure intrinsic properties of the single molecule over timescales from nanoseconds to seconds [2], without perturbations from surface attachment or other methods of immobilization. The general approach has been to combine real-time position sensing and closed-loop feedback actuation to compensate for thermally driven diffusion [24]. Both tasks are technically challenging, especially for single biomolecules or single small fluorophores, which are both fast and dim. High diffusion coefficients are unavoidable for nanoscale objects, placing a strong requirement for fast updating of position sensing and feedback. In addition, single molecules are limited in photon emission rate and total number of emitted photons before photobleaching occurs, so the trapping design must use each emitted photon as efficiently as possible.

The Anti-Brownian Electrokinetic (ABEL) trap, invented by Cohen and Moerner [46], is one elegant approach. The ABEL trap uses two-dimensional electrokinetic forces to push the object back to the target position once position deviation caused by Brownian motion is detected. Electrokinetic forces provide an effective handle on colloidal objects because they scale more favorably for small objects, both in terms of strength (compared to dielectrophoretic forces in optical tweezers) and frequency response (compared to mechanical stage motion). At the position sensing end, the first implementation of the ABEL trap used a CCD camera and centroid-based image processing to extract position in-formation of a fluorescently labeled Brownian object [4, 7] and achieved trapping of objects as small as 20 nm in size, limited only by the software update time (~4 ms). To trap smaller objects, an all-hardware version of the trap [6, 8] was developed that used a rotating beam/lock-in detection method [912] to sense position. These improvements reduced the update time and enabled successful trapping of 10 nm fluorescent objects in solution. However, it is still difficult to trap dimmer (<20 kHz count rate) and/or faster (D>40 µm2 s−1) objects; furthermore, extracting diffusion information from the object held in the trap directly from the feedback voltages is less straightforward compared to the camera-based design [5, 6].

Many feedback controller designs have been described [10, 13, 14]; however, they all require relatively high photon count rate in order to average among the many photon detection events within the integration window. Efficient feedback strategies in the low-count-rate regime are lacking, mainly due to the limitation of the highest achievable bandwidth set by the mechanical stage used in the tracking experiments [2, 3, 15]. However, with much faster electrokinetic feedback, we are in the good position to explore optimal strategies in the low-count-rate regime for the first time. In this paper, we propose a feedback controller specifically tailored to trap fast and dim objects in the ABEL trap platform, as well as a maximum likelihood estimator (MLE) to extract the diffusion coefficient of the trapped object. The rest of the paper is organized as follows. In Sect. 2, we introduce a new beam-scanning scheme and describe the structure of a real-time feedback algorithm. In Sect. 3, we develop the MLE for estimating diffusion coefficient. Simulation analyses of the proposed schemes are provided in Sect. 4. We discuss practical implementation issues and conclude in Sect. 5.

2 Beam scanning on a 2D grid and photon-by-photon position estimation

In feedback control of a Brownian object, position information is carried by the fluorescent (or scattered) photons.1 When designing the controller, one trades off localization precision for feedback bandwidth. Berglund et al. [11] suggest that the optimal feedback bandwidth for lockin detection should be chosen so that the rms distance of a Brownian displacement during the integration window (Δt) balances the localization error (1/Δt). In the low-count-rate/small-object regime, the photon detection rate ([Gamma with macron]) can easily be on the order of or less than the characteristic feedback bandwidth needed to trap (or track) a fast diffusing object (Γ¯4D/a02, a0 is the maximum displacement that can be tolerated before loss of the object). In this regime, it might be more advantageous to design a controller that provides feedback on every detected photon, as was described for the rotating beam case by Cohen et al. [6].

We propose a new beam scanning pattern on a 2D square lattice as illustrated in Fig. 1, where the Gaussian excitation spot dwells on each reference position (yj, j = 1–40) for a short period of time (~1 µs) before making a “knight’s tour” jump to the subsequent position. Such a scan pattern aims to encode more position information on each detected photon. After a complete traversal over the 40 points (defined as one frame), each spot will be covered exactly once, resulting in a flat excitation profile in the time-averaged sense (Fig. 1b). Such a “flat-top” excitation profile uniformly covers the trapping region, and would be beneficial to study fluorescence intensity fluctuations on time scales much longer than the beam scanning period [8].

Fig. 1
Two-dimensional beam scanning on a grid with full traverse defined as one frame. (a) Scan trajectory. (b) Frame-averaged intensity produced by the scan pattern in (a)

The tracking algorithm works as follows. For each detected photon, the position of the scan beam (one of the 40 possible positions in Fig. 1) at the time of the photon-detection event is recorded to be the unfiltered position estimate (yj ) of the object. We now examine the statistical meaning of the unfiltered measurement by seeking an expression for P(yj |x), namely, the probability of detecting a photon when the beam is at yj , given that the real position of the object is at x. If intrinsic intensity fluctuations that may arise from photophysics or conformational dynamics can be ignored, we write the instantaneous photon emission rate of the fluorescently labeled object as


where Γ0 is the count rate when the object overlaps with the center of the Gaussian excitation profile, x is the position vector of the object, yj is the position vector of the scanning beam at index j , and w is the 1/e2 waist (radius) of the excitation Gaussian profile. In the low-count-rate regime, the average number of photon detections per frame is smaller than unity. To make progress, we make two additional assumptions concerning the object’s motion:

  1. The movement (diffusion and translation) of the object is much smaller than the size of the scanning region, allowing x to be approximated as constant during the time of one frame
  2. The object remains close to the center of the trapping region, i.e., we consider situations when trapping is tight

The first assumption places an upper limit on the diffusion coefficient, and the second will be verified by direct simulation (vide infra): Under these assumptions, when a single photon is detected, the probability that the photon originates from scan point yj while the object is at x can be expressed as


where the Poisson probability of exactly one photon detection has been included.

We can now use (2) to develop a model for the measurement process, which is driven by the actual photon detection events at times τk. The photon detection selects one of the beam positions, which we call yk, indexed by k, because this set of positions is critical to the time evolution. Note that since our set of yj (Fig. 1) discretizes the 2D plane with uniform density, (2) means that yk samples the real position of the object (x) with Gaussian measurement noise nk ~ N(xk, Q), where the covariance matrix resulting from finite beam profile is given by Q = (w/2)2 I (I is the 2×2 identity matrix). The unfiltered measurement can thus be modeled by


Between photon detection events, the dynamics of the Brownian particle in the trap can be modeled in a fine-grained state space as [8]


where uk is the feedback voltage vector updated at photon detection events, μ is the mobility matrix, and wi is the Brownian increment term with wi ~ N(0, 2D dt I). In (3) and (4), index k denotes photon detection events, while index i represents small fixed sampling intervals (dt) where Brownian motion is updated. Under the conditions of Gaussian process noise and Gaussian measurement noise, the optimal state estimation problem can be solved analytically, and a recursive implementation is the celebrated Kalman filter [16].

Following the standard algorithm, we thus construct a state estimator comprised of a prediction update (or prior estimate, superscript –)


which describes how the mean of the estimate propagates according to the deterministic feedback forces, and a filtering update (or posterior estimate, superscript +)


which describes how the mean of the estimate gets refined given the newest measurement yk . Kk is the Kalman gain that weighs the importance of the most recent information. The feedback voltages are taken to be proportional the posterior estimate


where g is the voltage gain in units of V/µm. The duration of each electrokinetic kick is set to be


where Δtmax is a predetermined upper bound of feedback voltage duration, which might be useful in practice to prevent escape of the object during rarely occurring long photon waiting times. The estimator ((5) and (6)) and controller (7) are updated at photon detection events. We ignore the computational delay and electrokinetic response time in this model.

We now consider position tracking accuracy. Because the two Cartesian coordinates are independent and statistically equivalent (provided that Assumption 2 is satisfied), we restrict the following discussion to 1D. The variances of estimation errors, defined as Pk±=E[(xkx^k±)2], can be shown to propagate according to [16]



Equation (9) expresses the broadening of the estimation error variance due to diffusion, and (10) describes how the observer refines the estimation uncertainty after new measurement information becomes available. If the diffusion coefficient is known a priori, given the known photon arrival times, Pk can be updated at each measurement instant and used to determine the optimal Kalman gain that minimizes the posterior estimation error variance. We find


The optimal Kalman gain dynamically adjusts the importance of previous knowledge and new measurement information in forming the optimal posterior estimate (note that when the photon waiting time is large, 2D·(τk+1τk)(w/2)2,Kk+1opt1, which means the previous estimates become obsolete so that only the newest measurement is used in forming the posterior estimate in (6)). Practically speaking, we often face the problem in reverse: the diffusion coefficient is usually unknown, and it is desirable to extract information about the object’s diffusion statistics from experimentally available data. We show in Sect. 3 that such a parameter estimation problem can be solved by a maximum likelihood approach.

We note that the statistical properties of the estimation error and the expression for the optimal Kalman gain are only valid when there are no background photons. In case of low signal-to-background, the estimator and controller in (5), (6), and (7) can still be used, but the interpretation of each measurement becomes obscured. Due to the fact that D is unknown and the presence of background photons, a fixed Kk should be used to provide sub-optimal position tracking in practice. The value of Kk might be chosen by simulation or trial and error.

3 Maximum likelihood estimation of diffusion coefficient

Extracting the diffusion coefficient of the trapped (tracked) object can be easily done when the position can be sampled with considerable precision during the experiment [3, 5, 15]. Here we develop an algorithm to extract D from our photon-by-photon position estimation data.

In a real trapping experiment as described in Sect. 2, the available data include photon arrival times ({τk}), unfiltered measurement vectors ({yk}), prior and posterior estimates2 ({x^k±}), feedback voltage records ({uk}), excitation beam waist w, and Kalman gain ({Kk}) (or a fixed Kk ). We estimate D by maximizing the likelihood of observing N photons in sequence during the trapping experiment. The likelihood function is given by


where [Gamma]k is the (estimated) local rate of photon emission at time τk. How [Gamma] k is related to the diffusion coefficient D can be seen from the following argument. We start with the background-free case. At time τk−1, our knowledge of the object’s position is described by a Gaussian distribution centered at x^k1+, with a variance of Pk1+ along both Cartesian coordinates (assuming isotropic diffusion). After time τk − τk−1, right before the detection of the kth photon, the mean of the distribution translates to x^k by the applied feed-back force while the variance broadens to Pk due to diffusion (Fig. 2). Meanwhile, the excitation Gaussian profile is located at yk, with a variance of (w/2)2. The likelihood of emitting the kth photon is maximized when the prior estimation distribution and the excitation beam profile share the biggest overlap. The instantaneous rate of photon emission can therefore be expressed as the overlap integral,

Fig. 2
Maximum likelihood estimation of diffusion coefficient from a one step prediction process.N(x^k1+,Pk1+) and N(x^k,Pk) are 2D Gaussian distributions that characterize previous step posterior estimates and current step ...

For a given diffusion coefficient and the complete record of photon arrival times, the prior and posterior estimates as well as their error variances can be recovered at each photon detection event according to (5), (6), (9), and (10). With the reconstructed trajectory containing N prior estimates and their variances, the most likely value of D is the one that maximizes the overlap integral between two N × 2 dimensional Gaussians, constructed by all the prior distributions and the corresponding beam profiles, which can be seen mathematically by plugging (13) into (12).

We now consider the case with background photons. We model background photons as being generated by a Poisson process with constant rate Γbg. In calculating the instantaneous fluorescence rate as in (13), we take into account the possibility that photon detected at τk−1 might be a background photon. If that is the case, no position information is gained, so we heuristically replace the prior estimate distribution with a 2D uniform distribution (defined within the square-shaped trapping region) in the integrand. The approximation here is to ignore the fact that the error variance Pk in the presence of background photons should be modified. Specifically, let the signal-to-background ratio be SBR [equivalent] [Gamma with macron]bg; then


where the second integral is defined on the square-shaped trapping region with an area of (2xtrap)2. The likelihood function in (12) is modified consequently to be


The value of D maximizing (15) is determined numerically.

4 Simulation results

In this section, we validate the proposed algorithms by numerical simulations. The object’s two-dimensional dynamics is updated according to (4) with a time resolution of dt = 150 ns. The beam traversal pattern is updated simultaneously as illustrated in Fig. 1 with a dwell time of 1.05 µs on each reference position. The 1/e2 beam waist radius3 is set to be w = 0.4 µm, and the spacing of the 2D grid is 0.3 µm. The trapping region is approximated to be a square with edge length 2xtrap of 0.8 µm (central constant-intensity region of Fig. 1b). Photon detection is simulated as a Poisson process with a time-dependent rate determined by the distance between the excitation beam and the object (1). The observer is implemented as described by (5) and (6), while the feedback law is implemented according to (7) and (8). In the rare cases when more than one photon is detected within one dwell time of the scan beam, only the first photon is retained.

Figure 3 shows an example trapping trace, where a fixed Kalman gain is used, with simulation parameters chosen to mimic trapping of a small protein labeled with a single bright fluorophore. In Fig. 3a, the brightness fluctuations at one millisecond appear to arise from standard Poisson photon emission statistics at constant intensity, as opposed to the large fluctuations that would occur for a molecule diffusing through a simple Gaussian-shaped focal spot. (Of course, there are emission fluctuations near the frame rate frequency.) In Fig. 3b, we confirm that the object remains close to the center of the trapping region; moreover, the object would have diffused a root mean square distance of 4DΔt ~ 13 µm without the trap.

Fig. 3
Simulated performance of the ABEL trap controller as described in Sect. 2. The diffusion coefficient of the fluorescent object is 100 µm2/s. Other parameters include: g · μ = 7000 s−1 (overall feedback gain), Kk = 0.4, ...

In Fig. 4, we compare the algorithm developed in this work to the circular rotation/lock-in demodulation scheme [6, 15]. The comparison is made at constant signal to background ([Gamma with macron] = 14 kHz and Γbg = 2 kHz). Parameters for the circular rotation scheme are chosen to be rcirc = 0.85 µm, wcirc=2×rcirc and ωcirc = 2π×26 kHz (Fig. 4 inset left). This set of parameters guarantees that the time-averaged intensity is flat near the center [11, 13], and the size of the “top-hat” region (as well as the frame rate of scan) approximately matches that of the 2D grid scan scheme (Fig. 4 inset right). The lock-in algorithm generates a quadrature output by multiplying a rotating phasor reference by the pulse train of photon detection events, here averaged over one cycle of rotation. The resulting vector output of the lockin is used to generate the feedback voltages. Previous work [11] has shown that at high count rates, the quadrature outputs approach a Gaussian distribution in x and y centered at the origin. However, in the low-count-rate regime, the circular rotation algorithm degrades to a phase-only scheme, as evidenced by the cumulative vector plot of the lock-in outputs and radius histogram in Fig. 4b and c. Here, phase-only means that most of the feedback vectors are located on the circle. As a comparison, the outputs of our Kalman-filtered position estimator and a histogram of the estimates along one dimension are plotted in Fig. 4d and e. The estimated positions follow a Gaussian distribution, a direct result of the Gaussian measurement model (3) used in the controller design. Our new algorithm is thus more capable of estimating positions near the trapping center, presumably due to more positional information being encoded on every detected photon. We chose the standard deviation of position to gauge the performance of the two controllers in Fig. 4a. Our new algorithm is capable of trapping the objects more tightly, with up to 30% suppression of position fluctuation in the trap.

Fig. 4
Comparison of trapping performance between the circular rotation/lock-in demodulation scheme [6, 15] (black) and the algorithm developed in this work (red). The setup and time-averaged intensity of the circular rotation scheme are depicted in the inset ...

In Fig. 5, we show the performance of the maximum likelihood D estimator developed in Sect. 3. In this part of the simulation, objects are trapped at 14 kHz average signal, and the first 150 ms (~2400 detected photons) of trapping data are used to estimate the object’s diffusion coefficient using the likelihood function in (14) and (15). A constant Kalman gain of K = 0.4 is used. We tested the MLE for different levels of background photons, and it performs well at a signal-to-background ratio as low as 4.7 (Fig. 5b), although the uncertainty in the estimates increased significantly compared to the low-background case (Fig. 5a). The MLE estimator has also been tested with different values of trapping parameters (feedback gain, Kalman gain, beam waist, etc.) and proved to be quite robust except at the very lowest D (~1 µm2/s, where the object does not move much between photon detection events) and very highest D (>120 µm2/s, where the detected count rate is insufficient for stable trapping, and the assumptions of the model fail).

Fig. 5
Maximum likelihood estimation of diffusion coefficients using the algorithm developed in Sect. 3. The parameters used in the simulation are K = 0.4, g · μ = 8000 s−1, Δtmax = 76.8 µs. Only the first 150 ms of trapping ...

5 Conclusions

Using electrokinetic forces for feedback actuation enables controller designs that can operate at several tens of kHz and will hopefully make trapping of smaller and faster-moving objects possible. In this paper we have developed a controller that estimates the position and applies feedback control for every detected photon. We show that by scanning the excitation beam on a 2D square lattice, every detected photon yields more information about the object’s position, with an accuracy limited by the size of the scanning beam. In our approach, the dominant Gaussian measurement noise is a result of the TEM00 mode profile of the laser excitation spot, even at low count rate. Our position estimator is based on the well-established Kalman filter that takes advantage of the Gaussian statistics in providing (tractable) optimal state estimation. We also propose an MLE for estimating D, and the accuracy of this MLE is a direct validation of the measurement model in (3).

We expect the algorithm to work best in the regime where it is designed (i.e., fast diffusion with low signal count rate). To trap slow or bright objects, it is preferable to integrate the measurements to obtain sub-diffraction-limit localization accuracy before providing feedback [11, 13].

We envision the algorithm to be easily implemented on a Field-Programmable-Gate-Array platform [18]. Two-dimensional beam steering can be achieved by a pair of acousto-optical beam deflectors (AOBD), and fast switching can be accomplished by direct digital synthesizers [19], where the switching speed would be limited only by the aperture time of the AOBDs. In principle, one would like to scan a smaller beam using more points on a 2D grid for better position estimation accuracy. However, too small a scan beam would produce a large intensity variation in the z direction, and too many scan points would slow down the frame rate. The beam scanning scheme described here is thus a trade off in consideration of practical implementation issues. Moreover, we note that in calculating the MLE for the diffusion coefficient, (14) can be expressed in analytical form involving Gaussian and error functions (Erf(x)), so an on-line implementation is also possible. Real-time information about the object’s diffusion coefficient could be used to fine-tune the feedback loop as in (11) when the signal-to-background ratio is high. Such D estimates should also have a number of biophysical applications such as sorting biomolecules and detecting binding/unbinding events of different species.

Our goal is to optimally trap and study single bio-molecules in aqueous buffer. 100 kDa and 20 kDa proteins are estimated to have diffusion coefficientsz of ~60 µm2/s and ~100 µm2/s (in water), respectively [20]. The simulation results in this paper show that both proteins can be trapped tightly with only a few fluorescent labels and their diffusion statistics be recovered by our new approach.


We thank Sam Bockenhauer, Yan Jiang, and Dr. Randall Goldsmith for helpful discussions and Dr. Alexandre Fürstenberg for suggestions on an early version of the manuscript. This work is supported in part by the U. S. Department of Energy Grant No. DE-FG02-07ER15892 and by Grant No. 1R21-RR023149 from the National Center for Research Resources of the National Institutes of Health.


1In this paper, the focus will be on emitted fluorescence, since nano-sized objects are poor scatterers.

2We assume that electrokinetic mobility of the object can be characterized separately, for example, by an FCS approach [17].

3The corresponding confocal distance of the beam would be 2zR = 2πw2/λ = 1.5232 µm (for λ = 0.66 µm), still much larger than the z confinement depth of the ABEL trap (0.5–0.8 µm) [7]. Thus we expect trapping performance to be approximately the same as the object freely diffuses in the trapping channel in the z dimension.

Contributor Information

Q. Wang, Departments of Electrical Engineering and Chemistry, Stanford University, Stanford, CA 94305, USA.

W.E. Moerner, Departments of Chemistry and Applied Physics, Stanford University, Stanford, CA 94305, USA ; ude.drofnats@renreomw.


1. Cang H, Xu CS, Yang H. Chem. Phys. Lett. 2008 457;:285.
2. McHale K, Berglund AJ, Mabuchi H. Nano Lett. 2007;7(11):3535. [PubMed]
3. Cang H, Xu CS, Montiel D, Yang H. Opt. Lett. 2007;32(18):2729. [PubMed]
4. Cohen AE, Moerner WE. Appl. Phys. Lett. 2005;86:093109.
5. Cohen AE, Moerner WE. Proc. Natl. Acad. Sci. USA. 2006;103(12):4362. [PubMed]
6. Cohen AE, Moerner WE. Opt. Exp. 2008;16(10):6941. [PMC free article] [PubMed]
7. Cohen AE, Moerner WE. Proc. SPIE. 2005;5699:296.
8. Jiang Y, Wang Q, Cohen AE, Douglas N, Frydman J, Moerner WE. Proc. SPIE. 2008;7038(703807) [PMC free article] [PubMed]
9. Enderlein J. Appl. Phys. B. 2000;71:773.
10. Berglund AJ, Mabuchi H. Appl. Phys. B. 2004;78:653.
11. Berglund AJ, Mabuchi H. Appl. Phys. B. 2006;83(1):127.
12. Levi V, Ruan Q, Gratton E. Biophys. J. 2005;88:2919. [PubMed]
13. Andersson SB. Appl. Phys. B. 2005;80:809.
14. Cang H, Wong CM, Xu CS, Rizvi AH, Yang H. Appl. Phys. Lett. 2006;88:223901.
15. Berglund AJ, McHale K, Mabuchi H. Opt. Exp. 2007;15(12):7752. [PubMed]
16. Simon D. Optimal State Estimation: Kalman, H Infinity, and Non-linear Approaches. New York: Wiley; 2006.
17. Magde D, Webb WW, Elson EL. Biopolymers. 1978;17(2):361.
18. Stockton J, Armen M, Mabuchi H. J. Opt. Soc. Am. B. 2002;19(12):3019.
19. Vučinić D, Sejnowski TJ. PLoS ONE. 2007;2(8):e699. [PMC free article] [PubMed]
20. Tyn MT, Gusek TW. Biotech. Bioeng. 1989;35(4):327. [PubMed]