We present a comprehensive theory of closed-loop particle tracking for calculating the statistics of a diffusing fluorescent particle’s motion relative to the tracking lock point. A detailed comparison is made between the theory and experimental results, with excellent quantitative agreement found in all cases. A generalization of the theory of (open-loop) fluorescence correlation spectroscopy is developed, and the relationship to previous results is discussed. Two applications of the statistical techniques are given: a method for determining a tracked particle’s localization and an algorithm for rapid particle classification based on real-time analysis of the tracking control signal.
©2007 Optical Society of America
Various methods for controlling the Brownian diffusion of individual fluorescent particles have been developed over the past several years. Experimental demonstrations of closed-loop, single-particle tracking have included two-dimensional [1, 2, 3, 4] and three-dimensional [5, 6] tracking with piezo-electric actuation and two-dimensional control of individual [7, 8, 9] and multiple [10, 11] particles with electrophoretic actuation. Concurrently, theoretical results have been developed for understanding the statistics and performance limits of these new control systems [3, 4, 12, 13, 14, 15, 16, 10, 11, 17]. However, a comprehensive theory has not been previously available.
In this paper, we give a full theory of these control systems together with a detailed comparison to our own experimental results. In Sect. 2, we consider the theoretical description of a feedback-controlled, single-particle tracking experiment as a linear stochastic control system in state space (i.e., in the time domain). A detailed model is presented for the general case of a Brownian particle tracked by an arbitrary linear control law and subject to Gaussian position-sensor noise. We give a prescription for calculating statistical quantities such as the tracking error and various correlation functions and give explicit results for first-and second-order systems. For a fluorescent particle tracked (or trapped) in the focus of a spatially-modulated Gaussian excitation laser, fluorescence fluctuations depend on the competition between diffusion and control through the statistics of the tracking error. In Sect. 3, we exploit this fact to calculate generalized equations of closed-loop fluorescence correlation spectroscopy (FCS), which reduce to familiar results in the appropriate limits. Our model completely characterizes any linear fluorescent-particle tracking control system that uses a Gaussian excitation laser. In Sect. 4, we give an extensive comparison between the theory presented here and our own experimental results. In particular, we use the theoretical model to infer tracking control parameters. We then give two applications of such an analysis. First, we consider the problem of determining a particle’s localization, i.e. the standard deviation in the tracking error, and show how this parameter can be determined from typical tracking data in at least three complementary ways. Second, we use the statistics of the tracking control signal to develop a particle classification procedure for distinguishing between species in a binary mixture based on very little information collected over small spatial and temporal scales. The results presented in this paper are applicable not only to our own two-dimensional tracking results, but to all optical, linear-feedback particle tracking and trapping control systems, which differ (theoretically) only in the choice of reference frame.
2. Linear control system model
In this section, the basic model of a particle-tracking control system is presented in the language of linear stochastic control theory. The statistics of the resulting model will be calculated by solving an appropriate Fokker-Planck equation. The analytical methods and terminology used here can be found in many standard textbooks such as Refs. [18, 19, 20, 21].
Consider the feedback network shown in Fig. 1, where the tracking controller is represented by the transfer function C(s) and the control actuator (for example, the piezoelectric stage in our experiment) is represented by P(s). We assume that the control system does not couple the different Cartesian coordinates of a particle’s diffusion, so we may consider each coordinate separately. The tracked particle moves by Brownian motion with “velocity” U(t):
where D is the diffusion coefficient and dWp(t) is a stochastic Wiener increment. (As discussed in most textbook treatments of stochastic processes , the Wiener process Wp(t) is continuous but nowhere differentiable, and strictly speaking, the “velocity” U(t) is undefined. In our analysis, however, this term will always drive a linear differential equation, where it serves as convenient shorthand notation for a more carefully written, and well-defined, stochastic differential equation.) The time-integral of U(t) is the position of the particle, Xp(t), at time t:
The position of the sample stage is denoted by X(t), and the error signal
is the difference between the stage position X(t) and the particle’s position Xp(t). The control objective is to lock this error signal to zero, i.e., to achieve perfect tracking. In the experimental scenarios considered here, the tracked particle’s position Xp(t) is sensed by optical methods. In order to capture the noisy nature of such a sensor, Gaussian white noise
with a spectral density n (with dimensions of, for example, nm/) is added to the error signal; Wn(t) is another Wiener process statistically independent of Wp(t). For optical detection methods, the noise density n arises from photon counting statistics and no amount of technical sophistication or signal processing can suppress it. For otherwise optimal tracking, this sensor noise places a fundamental limit on the performance of the tracking controller [4, 16]. Furthermore, all of the optical parameters of a typical experiment, including a particle’s brightness, the excitation laser beam waist, and the point spread function of an imaging system, are absorbed into the noise density n and need not be considered in detail until we return to a discussion of fluorescence fluctuations.
2.1. Specification of transfer functions
In the previous section, we specified our tracking control system along one Cartesian coordinate as a linear control system in which the inputs, U(t) and N(t), drive the outputs, X(t) and E(t). The two output functions X(t), the sample stage position, and E(t), the deviation of the particle from the laser centroid, play central roles in analyzing a tracking experiment. In particular, during tracking, we cannot access the particle’s position Xp(t) directly; rather, we can only measure the sample stage position X(t), which tracks Xp(t) but is not identically equal to it. Similarly, we measure the particle’s fluorescence, which is a function of E(t), the deviation of the particle’s position from the excitation laser centroid. Therefore, we are particularly interested in the statistics of these two signals. In this section, we will calculate the joint, two-time probability distributions of the processes X(t) and E(t) for a generic stable control system of arbitrary order, and we will explicitly record the results for particular parameterizations of first-and second-order systems.
Let ~ represent the Laplace transform of a time-domain function, so that, for example, X̃(s) is the Laplace transform of X(t). In the Laplace domain, we now find a linear algebraic system
A transfer function T(s) is strictly proper if it is the ratio of two polynomials in s with the order of the denominator greater than the order of the numerator. It is stable if all of its poles, i.e., the zeros of its denominator polynomial, lie in the left half of the complex s-plane. Note that TEU(s) and TXU(s) as defined in Eq. (5) have poles at the origin s = 0, and it is at least possible (if this pole is not canceled by a zero in the numerator) that these are unstable. However, a pole at the origin represents a particularly innocuous form of instability, marginal stability. Because the transfer function representing the time derivative of such a process is stable, we will be able to handle these marginally stable cases with little difficulty.
Because our system is linear, each output is given by the sum of the individual response to each input. Letting X̃U(s) = TXU(s)Ũ(s) and similarly for X̃N(s), X̃U(s), and ẼN(s), we have
We may now consider each of the four quantities X̃U(s), X̃N(s), ẼU(s), and ẼN(s) separately and sum them to find the desired output statistics.
2.2. State-space realizations and the Fokker-Planck equation
Let us now take a step back to calculate the statistics of a generic, stable linear system. The following analysis will apply to any of the stable input-output pairs in our model; marginally stable cases will be treated in the next section. Consider a system with transfer function T(s) driven by an input U(t) ↔ Ũ(s) with output X(t) ↔ X̃(s). It is a straightforward task to find a state-space realization of the system  consisting of three matrices A, B, C, satisfying
Table 1. Table of useful first- and second-order statistics of the generic process X(t) defined by Eq. (11) of Sect. 2.2, for τ ≥ 0. ΔX Δt(t) = X(t + Δt) -X(t) is the discrete time-derivative of X(t) when data is sampled at time intervals Δt and E[∙] denotes an expectation value.
(A general state-space realization includes a fourth matrix D representing the direct feedthrough of input signal to output signal. We assume that the system under consideration is strictly proper, so that we may take D = 0.) If T(s) is of order m, then A, B, and C have sizes m×m,m×1 and 1×m respectively. The specification of these matrices is not unique, but as long as they satisfy Eq. (10), they represent a valid realization. Letting q(t) be an m-component internal state vector, the system’s dynamics can now be written in the form
Now consider a stochastic input U(t)dt = √αdW(t). The equation of motion for q(t) then becomes the multivariate Ornstein-Uhlenbeck process
q(t) is a vector-valued random process, whose statistics are given by the time-dependent probability distribution p(q, t) satisfying a linear Fokker-Planck equation:
where j and k represent the components of their corresponding vectors or matrices. When the system begins in state q 0 at time t = 0, the full time-dependent solution to time t = 0, the full time-dependent solution to Eq. (14) is given by the conditional probability distribution pt(q∣q 0):
represents an m-dimensional multivariate Gaussian distribution with mean vector m and (symmetric) covariance matrix ∑. The covariance matrix ∑t in Eq. (15) satisfies
For a stable system in which the eigenvalues of A all have negative real part [corresponding to left-half-plane poles of T(s)], Eq. (17) admits a finite stationary solution ∑∞ defined algebraically by the Lyapunov equation
The full solution to Eq. (17) can be written
The integral solution, Eq. (19), holds for all A even when ∑t becomes unbounded; for example, in the simple uncontrolled Brownian motion case where and A = 0 we find ∑t = 2Dt.
The two-time probability distribution for the state vector q(t) is easily found with Eqs. (15) and (20) and a little manipulation. Denoting the joint probability that q(t + t) = q 2 and q(t) = q 1 by p τ(q 2,q 1), we find for τ > 0,
We can now marginalize Eq. (22) to find the joint probability that X(t + τ) X(t) = X 1:
The statistics of the Gaussian process X(t) can be read off of the mean and covariance of Eq. (23). It will be useful to record the statistics of a few functions of X(t) as well. The results for various first- and second-order moments are recorded in Table 1.
For any strictly proper, stable transfer function T(s) driven by a stochastic input signal of the form dU(t) = √αdW(t), the results summarized in Table 1 give the statistics of the resulting output signal in terms of the state-space realization of T(s) and the solution of the Lyapunov equation, Eq. (18). For low-order systems, we may calculate these quantities explicitly (see the first- and second-order examples below). However, much of the analysis here is included in standard numerical analysis software. With these tools, it is a straightforward task to investigate quite complicated systems, including (for example) multiple resonances and time delays using polynomial (Padé) approximations.
2.3. Marginally stable systems
Some of the systems defined by the transfer functions of Eqs. (5) may be only marginally stable because they have a pole at the origin s = 0. In this section, we will modify the analysis of section 2.2 to account for this slight technical complication. To begin, consider the system X̃(s) = T(s)Ũ(s). If T(s) is marginally stable, then we can rewrite this system as
where T̄(s) = sT(s) is the closed-loop transfer function of a stable system. For this case, we can simply consider the time derivative
Table 2. Table of first- and second-order statistics of the marginally stable process X(t) defined by the state-space realization A̅, B̅, and C̅ as in Sect. 2.3.
This new system has stable dynamics represented by
and we can solve for the statistics of X˙(t) as before.
To begin, we solve for the statistics of X˙(t) in terms of a state-space realization A̅,B̅, C̅ and ∞ of T̅(s), with U(t)dt = √αdW(t) as in Sect. 2.2. In this case, the statistics of X(t) can be found by integration:
Note that we have suppressed an initial condition X(0). We quickly see that the mean is given by
Second-order moments can be constructed by explicit integration. For example, for τ ≥ 0, we consult Table 1 to find
After performing this and similar integrals, we find the results listed in Table 2.
2.4. First- and second-order systems
Having finished all the necessary calculations, we may now proceed to find the statistics for a few low-order systems of interest. In particular, consider the sample stage position X(t) and tracking error E(t) driven by the particle’s motion U(t) and measurement noise N(t) as described in Sect. 2. We consider two systems, specified by the transfer functions C(s) and P(s):
In both cases, we consider an integral control law C(s) with unity-gain frequency γc/2π. However, in the first case given by Eq. (30), the plant has a flat transfer function that can be driven arbitrarily hard with no amplitude or phase rolloff. This first-order system corresponds to the ideal tracking case, in which the bandwidth is set by the level of aggression of the control law, via γc. The second case given by Eq. (31) has the same control law, but the plant transfer function is now a low-pass filter, exhibiting both amplitude and phase rolloff at frequencies above γp/2π. We will use the second-order model to analyze our own experimental results in later sections.
In order to calculate statistical quantities, we must first find state-space realizations of the various transfer functions in Eq. (5). These are given in the Appendix for both the first- and second-order models of C(s) and P(s) and can be checked for consistency with Eq. (10). Using Tables 1 and 2 and the realizations given in the Appendix, we can explicitly calculate expectation values.
As an example, let us calculate the statistics of the tracking error, E(t), in the first-order model. For TEU(s), we find A = -γc, C = 1, and ∑∞ = D/γc so that, for example,
To construct the full statistics of E(t), we must consider the responses EU(t) and EN(t) due to both inputs. EU(t) and EN(t) are uncorrelated because they are driven by uncorrelated processes, so the statistics of the full error signal are just given by summing the means and variances calculated separately:
with . Because it can be driven arbitrarily hard with no degradation in amplitude or phase response, the first-order system is optimal and Eq. (33) represents the upper bound on tracking performance derived in Ref. .
For a richer example, consider the tracking error E(t) for the second-order system. Following the same procedure as for the first-order case, we find
where . Similarly, we find
These autocorrelation functions exhibit damped oscillatory behavior for γc > γp/4, where ν becomes imaginary, but they remain stable for all γp,γc > 0. The first-order tracking results are reproduced in the limit γp → ∞, that is, in the limit that the plant rolloff becomes much larger than the controller bandwidth.
3. Closed-loop fluorescence correlation spectroscopy
In Sect. 2, we found a prescription for calculating the statistics of the position of the sample stage X(t) and tracking error E(t). The former signal can be measured directly during a tracking experiment. The latter cannot be measured directly, but it can be sensed through the statistics of the fluorescence photon count rate. To see this, consider a one-dimensional system with a particle at position Xp(t) when the sample stage is at position X(t). For a fluorescence detectivity profile Φ(x) centered at the sample stage position, the rate of photon arrivals from this particle is given by
(We use the term “detectivity profile” to refer to the fluorescence photon detection rate taking into account the laser excitation profile, instrument response, and detection efficiency.)
Just as in the open-loop case of traditional fluorescence correlation spectroscopy (FCS) [22, 23, 24, 25], where fluorescence fluctuations arise from the particle’s free Brownian motion, the fluorescence autocorrelation function measured in closed-loop particle tracking probes the particle’s motion through the tracking error E(t) where the fluctuations arise from competition between free diffusion and feedback-assisted damping. In this section, we will calculate these fluorescence autocorrelation functions for Gaussian Φ(x), including the contribution from a deterministic spatial modulation pattern of the laser. These closed-loop calculations are a generalization of the open-loop FCS case, and we will show that they reduce to those results in the appropriate (weak feedback) limit. We will only calculate the expectation values of the fluorescence correlation functions; the variances could in principle be calculated using the same methods, because we have already derived the full Gaussian distribution of the relevant statistical quantities. As an example, see Ref.  for calculations of the variance in open-loop FCS.
3.1. Calculation of the fluorescence autocorrelation function
To begin, consider tracking a particle in the x-direction only. In many tracking scenarios, the laser moves in a deterministic (often circular) modulation pattern around a centroid position X(t). Denote the time-dependent offset of the laser from the beam centroid by xL(t), and let the detectivity function be Gaussian with beam waist w,
The (stochastic) fluorescence photon detection rate is
and the expectation value of the fluorescence correlation function is given by 
For convenience, let us define a shorthand vector notation, suppressing the time dependence of the error signal and the spatial path of the excitation laser:
x L is deterministic, and E is stochastic with Gaussian probability distribution p(E) characterized by its mean and covariance
where σ 0 2 and σ τ 2 were calculated in Sect. 2:
Using these expressions, we find the simple result
Finally, letting 0 2 = σ 0 2 + w 2/4, we find for the normalized fluorescence correlation function
The laser may also move along any time-dependent path described by x L as long as the lock point of the tracking control moves with it. For tracked diffusion in higher dimensions, even with asymmetric tracking and beam profiles, the full fluorescence autocorrelation function is just a product of terms of the form of Eq. (43), calculated along each Cartesian coordinate,
3.2. Recovery of open-loop results in the weak-tracking limit
In order to recover the standard results for open-loop FCS, we want to find the limit of G(t;τ) when the tracking is very weak. Because the system damping is represented by the matrix A, the weak-tracking limit can be found by letting this matrix approach 0. This limit is a bit tricky, however, because we must simultaneously let the standard deviation in the particle’s position σ 0 2 → ∞. That is, the damping becomes infinitesimally small while the particle’s confinement becomes correspondingly large. If we just blindly try to take the limit A → 0, we can easily find nonsensical results in which the unnormalized correlation function goes to 0 and the normalized version diverges, as the mean and variance of the fluorescence signal both go to zero at different rates. Furthermore, the closed-loop case includes only a single tracked particle, whereas the open-loop case is formulated for an ensemble of particles characterized by their concentration.
In order to take this open loop limit, we must simultaneously let A tend to zero and σ 0 2 approach infinity. To accomplish this, note that for very small A, we have
where we used the fact that ∑∞ is symmetric, and in the last line we took CB = , which can be seen by considering the case A = 0 for the system in Eq. (11).
We may suppress the distinction between the single-particle closed-loop case and the ensemble of particles in the open-loop case by fixing the prefactor in G(t;τ) at τ = 0. Using Eq. (47) and taking the large σ 0 2 limit, we find:
with . Apart from the prefactor, the final expression in Eq. (48) is the usual open-loop FCS result for one-dimensional motion in a time-dependent laser intensity.
3.3. Two-dimensional tracking in a rotating laser
In this section, we will find the specific form of the fluorescence autocorrelation function relevant to our own experiment. We track isotropically along the x and y directions, with the excitation laser rotating at angular frequency ω 0 and radius r, so that
The fluorescence autocorrelation function does not depend on t for this case:
which we can normalize to give
Finally, the deterministic oscillatory factor cos ω 0 τ may be distracting in the measured value of g(τ), but we can suppress it by averaging g(τ) over the rotation period. Denoting this averaged correlation function by g̅(τ), we have
where I 0 is the zeroth-order modified Bessel function, and the approximation holds when ω 0 is much larger than the feedback tracking bandwidth (the largest eigenvalue of A).
3.4. Behavior of g(τ)for τ ≈ 0
In a traditional open-loop FCS measurement, the value of the correlation function g(τ) at τ = 0 is a measure of the fraction of the beam which is filled with particles, or equivalently, it is a measure of the overlap between the beam profile and the distribution of particles in the sample. That is, it is a measure of the sample concentration. In that situation, a lower concentration leads to greater fluctuations, relative to the mean intensity. The one-dimensional fluorescence correlation function takes the form 
where N̅ is the average number of particles in the laser focus and the sample concentration can be determined from the relation g(0) = 1/N̅.
In closed-loop tracking, the value of g(τ) near τ = 0 is also a measure of the overlap of the trapped particle’s position distribution with the beam profile. However, in closed-loop tracking, there is only one particle in the laser focus at any time, and the concentration becomes difficult to define. Furthermore, as the tracking becomes better, the fluctuations decrease and g(τ) tends to 0. However, for the two-dimensional rotating laser case in our experiment, the value of the correlation function near τ = 0 still gives a measure of the particle’s confinement, i.e., the steady-state variance of the tracking error, σ 0 2.
To see this, first consider the simple case in which a particle is tracked but the laser is stationary, r = 0. From Eq. (51), we may define g 0 = g(τ = 0) to find
Although the expression is slightly more complicated for this closed-loop scenario, the basic physical principle is the same in the open- and closed-loop cases: the overlap between the distribution of particles (whether confined or free) and the detectivity profile determines g(0), the variance of the fluorescence fluctuations.
Now let us consider the experimentally relevant case of two-dimensional tracking in a circularly rotating laser, where g(τ) takes the form of Eq. 45 as discussed in Sect. 3.3. We suppose that the laser rotation frequency ω0 is much larger than any diffusion or control timescale (i.e. ω 0/2π is much larger than the largest eigenvalue of A). We may then assume that at τ = π/ω 0 (i.e., at the first minimum of cos ω 0 τ) we have σ τ=π/ω0 2≈ σ 0 2. Now define two quantities g 0 ± by
Using the preceding approximation for or σ τ=π/ω0 2 we find
Using Eq. (57), together with the value of r/w, we can determine a tracked particle’s confinement from the fluorescence correlation function through the value of g(τ) near τ = 0.
3.5. Relation to other literature results
Our analytical results are a generalization of the theory of open-loop correlation spectroscopy, and a number of models from the literature are contained in the general form of Eqs. (43),(45), and (46). For x L constant and A → 0, Eq. (46) reproduces the standard open-loop FCS result as shown above. For A → 0 and x L describing a two-dimensional circular orbit with the radius of rotation r much larger than the beam waist w, Eq. (46) reproduces the “fluorescence particle counting” results of Ref. . Under the same conditions on x L and A, but with an arbitrary radius of rotation, we find the recent results of Ref.  for the temporal autocorrelation in a laser scanning configuration. Note that all of these are open loop (A = 0) models.
The most closely related theoretical work in the literature is that of Enderlein , who studied the case of a modulated laser intensity tracing a circular orbit for use in tracking control. In that work, numerical simulations of closed-loop fluorescence correlation spectroscopy are presented including simple kinetic state transitions; particle escape probabilities, or tracking failures, are investigated as well. In this work, we have analytically solved the linearized version of that model, with no kinetic transitions. Note that a linearized model will be valid whenever tracking (or trapping) control is sufficiently good that the particle does not deviate far from the lock point . However, the inclusion of kinetic state transitions in the closed-loop FCS model discussed here is much more difficult than for open-loop FCS. This difficulty arises because a transition that affects a particle’s fluorescent brightness may consequently affect the tracking control dynamics, and when the generators of the tracking control dynamics for different internal states no longer commute, the model becomes largely intractable. Some results can be derived in certain limits , but no general solution has been developed. Analytical calculation of escape statistics requires a truly nonlinear model beyond the scope of this work.
4. Experimental results
In order to demonstrate the validity of our theoretical model, we tracked 60 and 210 nm diameter fluorescently labeled polystyrene nanoparticles in aqueous solution. Our experimental apparatus is depicted in Fig. 2, and is described in considerable detail in Ref. . It consists of a home-built fluorescence microscope in an epifluorescence configuration. We monitor fluorescent particles diffusing in a thin (~ 1 μm) liquid layer between two microscope coverslips. When a fluorescent nanoparticle diffuses away from the laser focus, its deviation from the laser centroid is detected, filtered by analog controller circuits and used to drive a piezoelectric stage in order to translate the sample and bring the particle back to the laser centroid. The 532 nm excitation laser beam is deflected in a circular pattern at angular frequency ω 0 = 2π×8 kHz, and a fluorescent particle’s x and y positions are detected in real time by phase-sensitive demodulation of the fluorescence signal at the rotation frequency. A CCD camera detects elastically scattered excitation light; the diffraction patterns caused by interference between the tightly focused incident and reflected beams  are used as a visual aid for focusing the microscope optics. An electronic servo automatically adjusts the laser power by feeding back to the radio frequency drive power of an acousto-optic modulator in order to stabilize either the fluorescent count rate or the excitation intensity, detected by a photodiode monitoring the leakage of excitation light through a dichroic filter. In the fluorescence stabilization mode, whenever a bright (dim) particle enters the laser focus, the servo system reduces (increases) the laser intensity to maintain a constant count rate at the detectors. Our background count rate is approximately 1000 s-1 and we typically lock to total photon count rates between 5 and 15 kHz (fluorescent count rates between 4 and 14 kHz). Because of this servo, nanoparticles with different fluorescence characteristics exhibit no difference in brightness during a tracking trajectory.
4.1. Tracking error and fluorescence fluctuations
Three tracking trajectories recorded from a single 60 nm fluorescent nanoparticle at three different fluorescent set points are shown in Fig. 3. These three trajectories were recorded under exactly the same experimental conditions, except that the intensity servo lock point was varied in order to track at three different fluorescence photon count rates. Because the tracking controller is linear, the overall loop gain of the tracking control system is proportional to the photon count rate. We may therefore use these three trajectories as examples of tracking control with different values of the overall loop gain. If we denote the x position of the sample stage during a tracking trajectory by X(t), then the mean-square deviation calculated over time intervals Δt provides an estimate of the diffusion coefficient:
where Var[∙] denotes the variance. For perfect tracking of a Brownian particle with diffusion coefficient D, we expect D̂ (Δt) = D, independent of Δt; however, for a realistic tracking control system, D̂(Δt) will not be equal to D for time intervals Δt that are small compared to the tracking bandwidth.
In order to compare the theory of Sect. 2 with our experimental results, we plotted D̂ (Δt) for each of the three trajectories displayed in Fig. 3. These plots are shown in Fig. 4, together with least-squares fits to the second-order model of Sect. 2.4, including a noise density n. We also calculated the fluorescence autocorrelation function g(τ) for each trajectory . These are displayed in Fig. 5 together with least-squares fits to the theory of Sect. 3. As the system gain increases along with the brightness of a tracked particle, we clearly see the increase in control bandwidth followed by the onset of oscillation at the highest gain (brightness) setting; the oscillation appears as a peaking behavior near Δt = 10 ms in D̂ (Δt) and as a revival near τ = 10 ms in g(τ). We find outstanding agreement between the theory and experiment, an indication that our model successfully accounts for fluctuations in both the tracking error and fluorescence count rate. The fit parameters for each set of curves are shown in Table 3, where we find good agreement between parameters determined through these two separate data channels.
4.2. Particle localization
As an application of the theory developed previously, we can use the parameters determined from fitting the values of D̂ (Δt) and g(τ) to infer the standard deviation in the tracking error, that is, the particle’s localization due to feedback control. Recall that the tracking error E(t) is the difference between the particle’s position and the laser centroid position. We define the localization L (along one axis) to be
We considered the calculation of this quantity from D̂ (Δt) in Ref. ; here, we discuss the estimation of L from both spatial information, through D̂ (Δt), and fluorescence fluctuations, through g(τ). For our investigation of the tracking localization L, we return to the same data set presented in Ref.  in which tracking trajectories were recorded from a binary mixture of roughly equal concentration 60 nm and 210 nm diameter fluorescent nanoparticles. For each trajectory, we calculate D̂ (Δt) and g(τ), fit these curves to the second-order-plus-noise model of Sects. 2-3 then calculate the localization using the theoretical expression of Sects. 2-3 then calculate the localization using the theoretical expression
which can be found from Eqs. (34–35) at τ= 0. We can determine the parameters D, γc, γp, and n in two independent ways, using fits to either D̂ (Δt) or g(τ). However, the value of D can be found much more reliably from the spatial information in D̂ (Δt), so we constrain the value of D in g(τ) according to the value found from D̂ (Δt). Finally, we can also determine L from the value of g(t) near t = 0 using Eq. (57) together with our calibrated values of the beam waist w = 1.0 μm and rotation radius r = 0.6 μm. In Fig. 6, we compare the value of L determined from these three methods, with the results summarized in Table 4. Note the unacceptably large spread in the measurement noise n determined from D̂ (Δt) for the smaller particles. We find the most reliable localization values to be those given by fitting g(t) to find γc, γp, and n with the diffusion coefficient D constrained by the asymptotic value of D̂ (Δt) at large Δt; this method combines the sensitivity of D̂ (Δt) at long times with the high resolution of g(τ) at short times.
4.3. Fast classification through hypothesis testing
As a final application of our results, in this section we will investigate the use of our statistical theory to classify particles as large (210 nm) or small (60 nm) using as little data as possible. In the microscopy literature, the phrase “single-particle tracking” typically refers to ex post facto reconstruction of individual particle trajectories via off-line analysis of a sequence of images captured with a fluorescence microscope . Such techniques are attractive for many purposes because they can be applied to wide-field images of complex biological systems. However, the amount and complexity of data processing required to extract individual particle trajectories make them inappropriate if one wants to react promptly to the detection of a particle trajectory satisfying some criterion of interest [33, 34]. For example, one might want to trigger an optical/mechanical procedure for extracting, catalyzing or destroying particles of interest as soon as they are identified; alternatively, one could consider adaptive reconfiguration of the microscope in order to focus light-collection and/or data acquisition resources on specific particles of interest. While fluorescence characteristics such as spectral or brightness information may be used to rapidly distinguish between particles in such a scenario, these are not always available and rapid characterization based on real-time motional statistics may extend the utility of these applications.
In contrast to these video tracking techniques, in our closed-loop tracking system, we have access to real-time information about a particle’s diffusional motion through the feedback tracking signal. We may therefore seek to classify particles based on the tracking feedback signal rapidly and with high fidelity. We may form an estimate of the diffusion coefficient along each direction by calculating the variance in the trajectory step size over N time intervals of length Δt, where the choice of N and Δt will determine the statistical accuracy of the estimate as well as the total estimation time T = NΔt. Over each time interval T, we may form an unbiased estimator of the diffusion coefficient D from D̂ (Δt) as defined in Eq. (58). Other estimators of D may be defined, emphasizing Bayesian analysis  or detection of dynamical changes in D , but we use D̂ (Δt) specifically for its ease of implementation and potential for real-time applications.
For Brownian motion, the increments in a particle’s position (in one Cartesian direction) obey Gaussian statistics with zero mean and variance given by 2DΔt, while the estimate D̂(Δt) obeys χ 2 statistics with mean value D . Because our feedback control system is linear, D̂(Δt) still obeys χ2 statistics even for small Δt, albeit with a mean value that deviates from the asymptotic (underlying) value of the particle diffusion coefficient D. This fact is confirmed in Fig. 7 for each type of particle (60 and 210 nm), where we show the distribution of D̂(Δt) for Δt = 10 ms together with the expected χ 2 distributions for varying sample numbers N.
Now consider a binary mixture, such as the one used here, consisting of a fraction λ 1 of particles of type 1 (diffusion coefficient D 1) and a fraction λ 2 = 1 - λ 1 of particles of type 2 (diffusion coefficient D 2 ≥ D 1). We wish to find a threshold value Dth such that we may assign a particle to class 1 if D̂(Δt) Dth and class 2 for D̂(Δt) ≥ Dth. Let Pcorr denote the probability that a classification is correct under this thresholding algorithm. A straightforward calculation shows that, for the expected χ 2 statistics of D̂(Δt), the value of Pcorr is maximized by choosing
D * th is a weak function of λ 1 and N for even moderately large N, and if λ 1 = λ 2, i.e. if the particles occur with equal likelihood (or we have no prior knowledge of their distribution), then D * th does not depend on N at all. D * th given by Eq. (61) may become negative or unbounded, but these limits simply indicate regimes in which a measurement is too noisy to warrant any correction to the prior distribution.
In a more general scenario, particle classification based on a measurement record may be formulated as a problem of hypothesis testing , which we briefly review. In the binary form of this problem, an m-component measurement is made with result θ, which may represent a single measurement or a sequence of measurements. The experimenter knows that this measurement result was drawn from one of two distributions, p 1(θ) or p 2(θ), with corresponding probabilities λ 1 and λ 2 = 1 - λ 1, and wishes to decide which of these distributions was most likely to have produced the observed value. Let H 1 represent the hypothesis that the underlying was p 1(λ) and similarly for H 2. Then we define a test procedure by the following threshold criterion: we accept hypothesis H 1 if λ 1 p 1(θ)p 1(θ) > λ 2 p 2(θ), and reject it otherwise. A standard theorem states that this test procedure minimizes the probability of an incorrect classification; furthermore, the Neyman-Pearson Lemma states that any other test procedure that decreases the probability of incorrectly accepting H 1 necessarily increases the probability of incorrectly accepting H 2 .
For our specific case of particle classification based on the scalar estimate D̂(Δt), the test described above divides the positive real line into regions corresponding to particles of type 1 (D̂(Δt) < D * th) and type 2 (D̂(Δt) ≥ D * th) where D * th is simply the point where
exactly as defined in Eq. (61). This general formulation can be applied to the case of a two-dimensional estimate of diffusion coefficients Dx and Dy along two Cartesian directions, even for the case that these are not identically distributed, i.e., the diffusion is not isotropic. In that case, the measurement vector θ = [D̂x(Δt), D̂y(Δt)] lies in a Dx-Dy plane, and the threshold criterion is a line dividing the plane into regions corresponding to each type of particle. For higher-dimensional measurement vectors, the hypothesis testing criterion is a surface partitioning the measurement space into regions corresponding to H 1 and H 2.
We implemented the above classification procedure on our data set consisting of 48 individual tracking trajectories, for various values of T, Δt and N. Because we form our estimate using short segments of very long trajectories, we may confirm whether a particular sample correctly identified a particle by comparing it to a high-fidelity identification based on the entire trajectory. In this way, we retain the ability to calculate the success probability, Pcorr. These results are shown in Fig. 8. At fixed estimation time T, Pcorr strictly increases with decreasing Δt down to about 10 ms. Beyond this point, the nanoparticles used here move (on average) only 250–500 nm and we are no longer able to make fast determinations of position with the accuracy required to form a good estimate D̂(Δt).
Note that at T = 60 ms we collect only (on average) 275 fluorescence photons but still identify particles with > 75% success; over this interval the larger (smaller) particles move on average only 600 (1150) nm. When the observation time is doubled to T = 120 ms, the success rate reaches 90%. At T = 1 s, the success rate is > 99% and even the faster particles move less than 5 μm. This method exhibits impressive fidelity based on information from very few photons collected over small spatial volumes and should be applicable to experiments involving single quantum dots or fluorescent biomolecules. Furthermore, these success rates are limited by the tracking error and feedback bandwidth, which are in turn nearly limited by photon counting shot-noise in our experiment . Note that since individual particles are tracked for times much longer than are required for accurate characterization, our method could be used to detect real-time dynamical changes in the diffusive behavior of an individual nanoparticle caused, e.g., by binding events or conformational switching (see also Ref. ).
In this work, we developed a detailed model of fluctuations in closed-loop fluorescent particle tracking using the language of linear stochastic control theory. Our results can be used to predict and analyze the statistics of both the control signal and the fluorescence fluctuations arising from competition between free diffusion and tracking control. We gave a detailed comparison between the theory and measurements made using our two-dimensional tracking apparatus, and found outstanding agreement in all cases. Two applications of the statistical theory, determination of a tracked particle’s localization and fast particle classification by hypothesis testing, were presented. Because tracking and trapping differ only in the choice of reference frame, the same results can be applied to both cases.
A. State-space realizations for the first- and second-order models
In this Appendix, we give explicit state-space realizations for the first- and second-order tracking models defined in Eqs. (30) and (31), respectively. The matrices A, B, and C, (or A̅, B̅, C̅ for the marginally stable TXU(s) case), together with Tables 1 and 2, give a prescription for calculating the statistics in both models. In each case, the transfer input-output transfer function is given by T(s)=C(s I - A)-1 B.
K. McHale acknowledges a National Defense Science and Engineering Graduate fellowship. This work was supported by the National Science Foundation under grant CCF-0323542 and by the Institute for Collaborative Biotechnologies through grant DAAD19-03-D-0004 from the Army Research Office. Any opinions, findings, and conclusions or recommendations expressed in this material are those of the authors, and do not necessarily reflect the views of either funding agency. A. Berglund’s present address: Center for Nanoscale Science and Technology, National Institute of Standards and Technology, Gaithersburg, MD 20899.
References and links
02. V. Levi, Q. Ruan, M. Plutz, A. S. Belmont, and E. Gratton, “Chromatin dynamics in interphase cells revealed by tracking in a two-photon excitation microscope,” Biophys. J. 89, 4275–4285 (2005). [CrossRef] [PubMed]
04. A. J. Berglund, K. McHale, and H. Mabuchi, “Feedback localization of freely diffusing fluorescent particles near the optical shot-noise limit,” Opt. Lett. 32, 145–147 (2007). [CrossRef]
05. V. Levi, Q. Ruan, and E. Gratton, “3-D particle tracking in a two-photon microscope. Application to the study of molecular dynamics in cells,” Biophys. J. 88, 2919–2928 (2005). [CrossRef] [PubMed]
06. H. Cang, C. M. Wong, C. S. Xu, A. H. Rizvi, and H. Yang, “Confocal three dimensional tracking of a single nanoparticle with concurrent spectroscopic readout,” Appl. Phys. Lett. 88, 223901 (2006). [CrossRef]
07. A. E. Cohen and W. E. Moerner, “Method for trapping and manipulating nanoscale objects in solution,” Appl. Phys. Lett. 86, 093109 (2005). [CrossRef]
10. S. Chaudhary and B. Shapiro, “Arbitrary steering of multiple particles independently in an electro-osmotically driven microfluidic system,” IEEE Trans. Contr. Syst. Technol. 14, 669–680 (2006). [CrossRef]
11. M. D. Armani, S. V. Chaudhary, R. Probst, and B. Shapiro, “Using feedback control of microflows to independently steer multiple particles,” IEEE J. Microelectromech. Syst. 15, 945–956 (2006). [CrossRef]
12. J. Enderlein, “Tracking of fluorescent molecules diffusing within membranes,” Appl. Phys. B 71, 773–777 (2000). [CrossRef]
13. J. Enderlein, “Positional and temporal accuracy of single molecule tracking,” Sing. Mol. 1 , 225–230 (2000).
14. A. J. Berglund and H. Mabuchi, “Feedback Controller design for tracking a single fluorescent molecule,” Appl. Phys. B 78, 653–659 (2004). [CrossRef]
15. S. B. Andersson, “Tracking a single fluorescent molecule in a confocal microscope,” Appl. Phys. B 80, 809–816 (2005). [CrossRef]
16. A. J. Berglund and H. Mabuchi, “Performance bounds on single-particle tracking by fluorescence modulation,” Appl. Phys. B 83, 127–133 (2006). [CrossRef]
17. D. Montiel, H. Cang, and H. Yang, “Quantitative characterization of changes in dynamical behavior for single-particle tracking studies,” J. Phys. Chem. B 110, 19763–19770 (2006). [CrossRef] [PubMed]
18. O. L. R. Jacobs, Introduction to Control Theory, 2nd ed. (Oxford University Press, 1996).
19. N. G. Van Kampen, Stochastic processes in physics and chemistry (Elsevier Science Pub. Co., North-Holland, Amsterdam, 2001).
20. C. W. Gardiner, Handook of Stochastic Methods for Physics, Chemistry and the Natural Sciences, 2nd ed. (Springer-Verlag, 1985).
21. H. Risken, The Fokker-Planck Equation: Methods of Solution and Applications, 2nd ed. (Springer, 1959).
22. D. Magde, E. L. Elson, and W. W. Webb, “Thermodynamic fluctuations in a reacting system - measurement by fluorescence correlation spectroscopy,” Phys. Rev. Lett. 29, 705–708 (1972). [CrossRef]
23. E. L. Elson and D. Magde, “Fluorescence correlation spectroscopy. 1. Conceptual basis and theory,” Biopolymers 13, 1–27 (1974). [CrossRef]
25. O. Krichevsky and G. Bonnett, “Fluorescence correlation spectroscopy: the technique and its applications,” Rep. Prog. Phys. 65, 251–297 (2002). [CrossRef]
28. M. A. Digman, C. M. Brown, P. Sengupta, P. W. Wiseman, A. R. Horwitz, and E. Gratton, “Measuring fast dynamics in solutions and cells with a laser scanning microscope,” Biophys. J. 90, 1317–1327 (2005). [CrossRef]
29. A. J. Berglund, “Feedback Control of Brownian Motion for Single-Particle Fluorescence Spectroscopy,” Ph.D. thesis, California Institute of Technology (2006), http://etd.caltech.edu/etd/available/etd-10092006-165831/.
30. L. Novotny, R. D. Grover, and K. Karrai, “Reflected image of a strongly focused spot,” Opt. Lett. 26, 789–791 (2001). [CrossRef]
32. M. J. Saxton and K. Jacobson, “Single-particle tracking: applications to membrane dynamics,” Annu. Rev. Bio-phys. Biomolec. Struct. 26, 373–399 (1997). [CrossRef]
33. S. Bonneau, M. Dahan, and L. D. Cohen, “Single quantum dot tracking based on perceptual grouping using minimal paths in a spatiotemporal volume,” IEEE Trans. Image Process. 14, 1384–1395 (2005). [CrossRef] [PubMed]
34. E. Meijering, I. Smal, and G. Danuser, “Tracking in Molecular Bioimaging,” IEEE Signal Processing Mag. 23, 46–53 (2006). [CrossRef]
36. M. H. DeGroot, Probability and Statistics (Addison-Wesley, Reading, MA, 1986).