1. Introduction

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 dW_p(t) is a stochastic Wiener increment. (As discussed in most textbook treatments of stochastic processes [20], the Wiener process W_p(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, X_p(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 X_p(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 X_p(t) is sensed by optical methods. In order to capture the noisy nature of such a sensor, Gaussian white noise

 

Fig. 1. Block diagram of the particle tracking control system. The system is driven by the particle’s velocity U(t) and the measurement noise N(t); the system outputs are the tracking stage position X(t) and the error signal E(t). For perfect tracking control, E(t) = 0.

Download Full Size | PDF

with a spectral density n (with dimensions of, for example, nm/$\sqrt{\text{Hz}}$) is added to the error signal; W_n(t) is another Wiener process statistically independent of W_p(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 X_p(t) directly; rather, we can only measure the sample stage position X(t), which tracks X_p(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

The four transfer functions in Eq. (5) are given by inspection of Fig. 1 and standard block-diagram algebra [18]. Letting L(s) = C(s)P(s) denote the loop transfer function, we find

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 T_EU(s) and T_XU(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) = T_XU(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 [18] 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 ₀ 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 p_t(q∣q ₀):

where

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 $\sqrt{\alpha }\mathbf{B}=\sqrt{2D}$ 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 ₂ and q(t) = q ₁ by p _τ(q ₂,q ₁), we find for τ > 0,

We can now marginalize Eq. (22) to find the joint probability that X(t + τ) X(t) = X ₁:

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 $\overline{\sum }$ _∞ 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 T_EU(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 E_U(t) and E_N(t) due to both inputs. E_U(t) and E_N(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 $\overline{D}=D+\frac{n^2{\gamma }_c^3}{2}$. 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. [16].

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 $\nu =\sqrt{{\gamma }_p^2-4{\gamma }_p{\gamma }_c}$. 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 X_p(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. [26] 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 x_L(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 [3]

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 σ ₀ ² and σ _τ ² were calculated in Sect. 2:

Using these expressions, we find the simple result

with

Finally, letting $\overline{\sigma }$ ₀ ² = σ ₀ ² + w ²/4, we find for the normalized fluorescence correlation function

Eqs. (43) and (45) represent the fluorescence autocorrelation function for a particle tracked by an arbitrary linear feedback law, represented through the covariances σ ₀ ² and σ _τ ² of Eq. (42).

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 σ ₀ ² → ∞. 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 σ ₀ ² 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 = $\sqrt{2D}$, 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 σ ₀ ² limit, we find:

with ${\tau }_D=\frac{w^2}{4D}$. 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 ω ₀ 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 ω ₀ τ 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 ₀ is the zeroth-order modified Bessel function, and the approximation holds when ω ₀ 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 [25]

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, σ ₀ ².

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 ₀ = 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 ω₀ is much larger than any diffusion or control timescale (i.e. ω ₀/2π is much larger than the largest eigenvalue of A). We may then assume that at τ = π/ω ₀ (i.e., at the first minimum of cos ω ₀ τ) we have σ _τ=π/ω0 ²≈ σ ₀ ². Now define two quantities g ₀ ^± by

and let

Using the preceding approximation for or σ _τ=π/ω0 ² 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. [27]. Under the same conditions on x _L and A, but with an arbitrary radius of rotation, we find the recent results of Ref. [28] 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 [13], 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 [16]. 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 [29], but no general solution has been developed. Analytical calculation of escape statistics requires a truly nonlinear model beyond the scope of this work.

 

Fig. 2. Schematic diagram of the experimental apparatus for tracking freely diffusing fluorescent nanoparticles in two dimensions as described in the main text. Key: AOM, acousto-optic modulator; HWP, half-wave plate; PBS, polarizing beam splitter; QWP, quarter-wave plate; CCD, charge-coupled device camera; PD, photodiode; APD, avalanche photodiode single-photon counter; HV, high voltage.

Download Full Size | PDF

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. [29]. 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 ω ₀ = 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 [30] 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.

 

Fig. 3. Tracking trajectories for the same particle at three different excitation intensities. The rate of fluorescent photon arrivals during each trajectory is shown in the upper plots, while the x and y positions of the sample stage are shown in the lower plots. A single particle was tracked for around 100 s. After each 20 s interval, data collection was paused for a few seconds while the intensity servo set point was changed manually; these pauses are indicated by vertical black lines, but are not shown on the time (t) axis.

Download Full Size | PDF

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 [31]. 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.

 

Fig. 4. Mean-square deviations D̂ (Δt)/2Δt for each of the trajectories in Fig. 3, plotted together with fits to the theory developed in Sect. 2.

Download Full Size | PDF

 

Fig. 5. Fluorescence autocorrelation functions g(τ) for each of the trajectories in Fig.3, plotted together with fits to the theory developed in Sect. 3. The curves have been offset along the vertical axis for clarity.

Download Full Size | PDF

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

 

Fig. 6. Three-dimensional scatter plot of the localization L determined from D̂ (Δt), g(τ), and g ₀ as described in the text. Two-dimensional scatter plots comparing each pair of methods are shown in lighter shades projected onto their respective axes. Localization values determined for the 60 nm beads are shown in blue and 210 nm beads are shown in red. Dashed lines indicate where values would lie on the projection planes if all estimates were identically equal.

Download Full Size | PDF

We considered the calculation of this quantity from D̂ (Δt) in Ref. [4]; 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. [4] 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

Table 3. Table of fit parameters for the mean-square deviation curves of Fig. 4 and fluorescence autocorrelation functions of Fig. 5. Each set of fit parameters is labeled above by the photon count rate in kHz. For the fits to g(τ), the beam waist and diffusion coefficient were constrained to w = 1 μm (determined by scanning the excitation laser over an immobilized fluorescent nanoparticle) and D = 5.1 μm²/s (the average value determined from the mean-square deviation curves); constrained parameters are indicated by ^*. Aside from the constraint on D in g(τ), the two parameter sets are otherwise independent.

View Table | View all tables in this article

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 [32]. 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.

Table 4. Table comparing the localization L and measurement noise n determined from D̂ (Δt), g(τ), and g ₀ for both 60 and 210 nm diameter particles. Shot-noise limited, theoretical optimum values are denoted by “Opt.” The error values are the observed standard deviation.

View Table | View all tables in this article

 

Fig. 7. Histogram of estimates D̂ (Δt = 10ms) calculated for a single trajectory of each type of particle (60 and 210 nm) versus the number of samples N used for the estimate. The solid black curves show the expected distribution with mean value set equal to the mean value from each data set.

Download Full Size | PDF

 

Fig. 8. Measured probability of correct classification P_corr as the estimation time T and the sample time Δt are varied. The dashed curves show the expected success probability calculated from χ² statistics.

Download Full Size | PDF

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 [35] or detection of dynamical changes in D [17], 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 χ ² statistics with mean value D [36]. Because our feedback control system is linear, D̂(Δt) still obeys χ² 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 χ ² distributions for varying sample numbers N.

Now consider a binary mixture, such as the one used here, consisting of a fraction λ ₁ of particles of type 1 (diffusion coefficient D ₁) and a fraction λ ₂ = 1 - λ ₁ of particles of type 2 (diffusion coefficient D ₂ ≥ D ₁). We wish to find a threshold value D_th such that we may assign a particle to class 1 if D̂(Δt) D_th and class 2 for D̂(Δt) ≥ D_th. Let P_corr denote the probability that a classification is correct under this thresholding algorithm. A straightforward calculation shows that, for the expected χ ² statistics of D̂(Δt), the value of P_corr is maximized by choosing

D ^* _th is a weak function of λ ₁ and N for even moderately large N, and if λ ₁ = λ ₂, 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 [36], 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 ₁(θ) or p ₂(θ), with corresponding probabilities λ ₁ and λ ₂ = 1 - λ ₁, and wishes to decide which of these distributions was most likely to have produced the observed value. Let H ₁ represent the hypothesis that the underlying was p ₁(λ) and similarly for H ₂. Then we define a test procedure by the following threshold criterion: we accept hypothesis H ₁ if λ ₁ p ₁(θ)p ₁(θ) > λ ₂ p ₂(θ), 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 ₁ necessarily increases the probability of incorrectly accepting H ₂ [36].

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 D_x and D_y 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 D_x-D_y 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 ₁ and H ₂.

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, P_corr. These results are shown in Fig. 8. At fixed estimation time T, P_corr 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 [4]. 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. [17]).

5. Conclusion

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.