Optimal tracking of a Brownian particle

Alexander P. Fields; Adam E. Cohen

doi:10.1364/OE.20.022585

1. Introduction

The motion of a molecule in solution encodes information about its size, interactions, and surroundings. Single-particle tracking experiments [1] have yielded information about the dynamics of molecules in the plasma membrane [2,3], the walking mechanisms of motor proteins [4], and intracellular transport of nucleic acids [5] and protein complexes [6]. Anti-Brownian traps combine optical tracking with active feedback, applied via stage motion or electrokinetic drift, to confine a single particle within an observation volume for an extended time [7]. These devices have yielded new insights into the dynamics of biomolecules [8–10] and nanoparticles [11] but face the additional challenge that the tracking must be performed in close to real time.

In choosing an imaging strategy, one typically faces a tradeoff between photon shot noise and Brownian motion. During a brief exposure a particle remains relatively localized in space, but the precision of the measurement is degraded by photon shot noise. Diffractive blurring causes each photon to carry imperfect information about the location of the particle, and background photons can contribute spurious information. During a long exposure, the photon statistics are improved but the precision of the measurement is degraded by Brownian motion. Specialized hardware can reach exposure times shorter than the mean interval between photons, so that each “image” consists of at most a few photons [12], while video cameras can achieve arbitrarily long exposure times.

Given a series of noisy position measurements, the fundamental challenge of particle tracking is to link these data into a trajectory. Naively, one might simply concatenate successive measurements. However, recent measurements often contain useful information about the present location, and thus one can combine data taken at different times to compute a more likely trajectory than is obtained by simple concatenation. Here we derive a hierarchy of algorithms for calculating likely trajectories. These algorithms allow tradeoffs between accuracy and computational efficiency.

For post-processing of tracking data we present an Assumed Density Filter (ADF) algorithm that achieves near-optimal performance in particle tracking by applying an assumed model of particle motion. This algorithm is useful to obtain maximum likelihood estimates of molecular transport coefficients (electrokinetic mobility and diffusion coefficient). Anti-Brownian traps require real-time information on particle location. We develop a Kalman filter algorithm which is a simplified and computationally efficient form of the ADF. The Kalman filter presumes knowledge, or an estimate, of the transport coefficients, and seeks to estimate the position of the particle. These theoretical results proved crucial to our efforts to trap individual fluorophores and other molecules in solution and to estimate the transport coefficients of trapped molecules [13].

Our ADF algorithm is a discrete-time recursive Bayesian estimator [14,15]. This approach uses a stochastic model of particle motion to combine a series of noisy measurements into a maximum likelihood estimate of the present position. We model the motion as Brownian diffusion plus electrokinetic drift, and we explicitly include Poisson-distributed photon shot noise and background photons. A related model [16–18] previously treated the noise and background as Gaussian-distributed white noise. The key innovations in the present work are to use a nonlinear filter, the ADF, to enhance tracking accuracy, and to account for the discreteness of photon arrivals in both the Kalman filter and the ADF. We present our model in the context of the anti-Brownian electrokinetic trap (ABEL trap) [19], but the treatment can be adapted to other uses (e.g. for other feedback modalities or in the absence of feedback).

The ABEL trap uses fluorescence microscopy to track a particle diffusing in two dimensions, and applies electrokinetic forces to oppose the observed motion. Out-of-plane motion is blocked by the walls of the sample chamber. The earliest ABEL trap [19] imaged particles in the trap on a camera and was consequently limited by frame readout time. Recent incarnations [13,20] scan a focused laser beam among a set of discrete spots in the sample plane, simultaneously tallying photon arrivals on a single-point detector, and determine particle location from the position of the laser at the instant each photon is detected. Our implementation scans the laser among 27 spots in a hexagonal lattice, with a dwell time of 3.1 μs at each spot [Fig. 1 ]. The challenge is to construct a likelihood distribution for the particle’s position given the history of detected photons. We advise caution in applying the mathematics described here when the underlying motion is not Brownian (e.g. active transport); such cases would be better treated by constructing a variety of models and selecting the one best fit by the data.

Fig. 1 Tracking and feedback in an anti-Brownian electrokinetic trap. A fluorescent particle diffuses in a quasi-two-dimensional sample chamber at the center of four platinum electrodes. A focused laser is scanned through a set of discrete spots in the sample plane, triggering the particle to emit fluorescence when the laser overlaps with the particle. Detected photons are tallied and assigned to laser spots according to their detection time. A Kalman filter combines this data with a prediction for the location of the particle based on previous data to generate an updated position estimate. The position estimate is used to calculate feedback voltages which are applied to counteract the particle’s motion.

Download Full Size | PDF

2. General treatment

We consider the case of a single fluorescent particle diffusing in free solution, where the particle diameter is much less than the optical wavelength. This particle might have complex internal structure, but the structure is not optically resolvable, so we treat the particle as a point-like emitter and refer to it generically as a “fluorophore”. Although our instrument tracks particles in two dimensions, the mathematical description is similar for one- or three-dimensional geometries.

Assuming laser intensity below saturation, the rate γ at which photons are detected is well approximated as a Gaussian function of the displacement of the fluorophore from the center of the laser spot, plus a constant background rate [21,22]. Mathematically,

γ (x) = s \exp (- \frac{1}{2} {(x - c)}^{T} W^{- 1} (x - c)) + b,

where s is the maximum “signal” photon rate (in photons per second) from the fluorophore, b is the “background” photon rate (also in photons per second), x is the particle position (written as a column vector), c is the laser spot center, W is the spatial covariance of the spot, ^T denotes transposition, and ⁻¹ denotes matrix inversion. The probability of detecting n photons during a time bin of length ∆t during which the laser is stationary is the Poisson distribution

\Pr [n | x] = \frac{{(γ (x) Δ t)}^{n} e^{- γ (x) Δ t}}{n!} .

Throughout the experiment, the fluorophore diffuses with diffusion coefficient D. In the case of the ABEL trap, a series of electric fields E_k are applied to the trap, imparting motion to the fluorophore in proportion to its electrokinetic mobility μ. The probability distribution of the molecule’s motion ∆x during time step k is

g_{k} (Δ x) = \frac{\exp (- \frac{{| Δ x - μ Δ t E_{k} |}^{2}}{4 D Δ t})}{4 π D Δ t} .

This equation is indifferent to the method by which the applied electric field E_k is selected. In the case of the ABEL trap, the field strength is modulated so as to counteract the observed diffusion,

E = - \frac{\hat{x}}{μ Δ t},

where

\hat{x}

is the most recent estimate of the molecule’s position, but Eq. (3) and the following derivation remain valid if E is held constant or modulated independently of the molecular motion.

The transport coefficients μ and D reflect the charge, size, and interactions of the particle, so extracting these (possibly time-dependent) parameters is a goal of particle tracking and feedback trapping. The molecular brightness s, which is sensitive to the electrochemical environment of the fluorophore, may also be of interest.

Let {n_j, j ≤ k} denote the complete history of photon counts in discrete time bins up to and including bin k, where each bin corresponds to a new laser position. We seek to construct the likelihood distribution Pr[x_k|{n_j, j ≤ k}] for the position of the fluorophore during the k^th bin. The position likelihood distribution f_a_|_b ≡ Pr[x_a|{n_j, j ≤ b}] is estimated recursively from:

f_{k | k} = L_{k}^{- 1} P r [n_{k} | x_{k}] f_{k | k - 1}

f_{k | k - 1} = \int g_{k - 1} (x_{k} - x_{k - 1}) f_{k - 1 | k - 1} d x_{k - 1}

where the integral is over the entire two-dimensional plane and L_k ≡ Pr[n_k|{n_j, j ≤k − 1}]. Equations (5) and (6) are known as “update” and “prediction” equations, respectively, and together constitute a “recursive Bayesian estimator” of the fluorophore position. To be applicable, these equations require that the molecule’s unperturbed motion be truly Brownian (e.g., no significant momentum), that the effect of the electric field E be independent of the molecule’s position x, and that the motion of the fluorophore during each bin be negligible.

3. Gaussian assumed-density filter (ADF)

In principle, Eqs. (1-6) enable exact calculation of the likelihood distribution of particle position given a set of photon counts and corresponding laser positions. In practice, it is unfeasible to perform this calculation in real-time at the speed required to track fast-moving fluorophores – indeed, in our experience the full calculation is too slow even for post-processing (see section 6). The difficulty stems from the complicated shape of the measurement distribution [Eq. (2)], which renders the estimate distribution analytically intractable if no approximations are made. On the other hand, treating the full probability distributions numerically is impractical because of the slowness of two-dimensional numerical convolution [Eq. (6)], which needs to be computed every time the laser moves to a new spot (every 3.1 μs in our implementation). Propagating the distribution on a 128 × 128 grid, for instance, requires >10⁵ multiplication operations for each convolution. It is therefore necessary to make approximations. Myriad schemes have been devised to handle nonlinear measurements (several of which are described in [23]). The Gaussian assumed-density filter (ADF) [24,25] is well suited to particle tracking in the presence of shot noise, Brownian motion, and background fluorescence. An ADF is a type of recursive Bayesian estimator in which the complicated shape of the likelihood densities f_a_|_b are approximated by simpler, mathematically tractable distributions.

Let us assume that the likelihood functions f_a_|_b are normally distributed with mean ${\hat{x}}_{a | b}$ and covariance ${\hat{Σ}}_{a | b}$ ; that is,

f_{a | b} \approx Ν ({\hat{x}}_{a | b}, {\hat{Σ}}_{a | b}) \equiv {(2 π {| {\hat{Σ}}_{a | b} |}^{\frac{1}{2}})}^{- 1} \exp (- \frac{1}{2} {(x - {\hat{x}}_{a | b})}^{T} {\hat{Σ}}_{a | b}^{- 1} (x - {\hat{x}}_{a | b})) .

This shape of distribution vastly simplifies the prediction step [Eq. (6)], because the convolution of two normal distributions is another normal distribution. Therefore, we can immediately reduce Eq. (6) to

\begin{matrix} {\hat{x}}_{k | k - 1} = {\hat{x}}_{k - 1 | k - 1} + μ Δ t E_{k - 1} \\ {\hat{Σ}}_{k | k - 1} = {\hat{Σ}}_{k - 1 | k - 1} + 2 D Δ t I . \end{matrix}

On the other hand, the update step [Eq. (5)] does not preserve the Gaussian shape of the likelihood distribution. Instead, combining Eqs. (1,2,5,7) produces

\begin{matrix} f_{k | k} = L_{k}^{- 1} {(n_{k}!)}^{- 1} {(S \exp (- \frac{1}{2} {(x - c)}^{T} W^{- 1} (x - c)) + B)}^{n_{k}} \\ \times \exp (- (S \exp (- \frac{1}{2} {(x - c)}^{T} W^{- 1} (x - c)) + B)) Ν ({\hat{x}}_{k | k - 1}, {\hat{Σ}}_{k | k - 1}), \end{matrix}

where S ≡ s∆t and B ≡ b∆t. This distribution is no longer normally distributed, but we can approximate it as such by calculating its mean and covariance and dropping higher moments. An alternative strategy would be to project the posterior distribution onto a finite sum of Gaussian distributions, similar to a Gaussian sum filter [26]. In principle, this approach could increase the accuracy of the estimator, at the expense of greater computational complexity; we have not attempted such a strategy.

To calculate the mean and covariance of the distribution [Eq. (9)], we first rearrange it as

f_{k | k} = L_{k}^{- 1} \sum_{m = 0}^{\infty} l_{m} Ν ({\hat{χ}}_{m}, {\hat{Ψ}}_{m})

where

\begin{matrix} l_{m} \equiv e^{- B} \sum_{i = \max (n_{k} - m, 0)}^{n_{k}} \frac{{(- 1)}^{m + i - n_{k}} B^{i} S^{m}}{(n_{k} - i)! i! (m + i - n_{k})!} \frac{{| {\hat{Ψ}}_{m} |}^{\frac{1}{2}}}{{| {\hat{Σ}}_{k | k - 1} |}^{\frac{1}{2}}} \\ \times \exp (- \frac{1}{2} m {({\hat{x}}_{k | k - 1} - c_{k})}^{T} {(m {\hat{Σ}}_{k | k - 1} + W)}^{- 1} ({\hat{x}}_{k | k - 1} - c_{k})) \\ {\hat{χ}}_{m} \equiv {\hat{Ψ}}_{m} ({\hat{Σ}}_{k | k - 1}^{- 1} {\hat{x}}_{k | k - 1} + m W^{- 1} c_{k}) \\ {\hat{Ψ}}_{m} \equiv {({\hat{Σ}}_{k | k - 1}^{- 1} + m W^{- 1})}^{- 1} . \end{matrix}

Equation (10) rewrites Eq. (9) as an infinite sum of normal distributions. The overall normalization factor L_k is simply

L_{k} = \sum_{m = 0}^{\infty} l_{m} .

The mean and covariance are

{\hat{x}}_{k | k} = \int x f_{k | k} d x = L_{k}^{- 1} \sum_{m = 0}^{\infty} l_{m} {\hat{χ}}_{m}

{\hat{Σ}}_{k | k} = (\int x x^{T} f_{k | k} d x) - {\hat{x}}_{k | k} {\hat{x}}_{k | k}^{T} = (L_{k}^{- 1} \sum_{m = 0}^{\infty} l_{m} ({\hat{χ}}_{m} {\hat{χ}}_{m}^{T} + {\hat{Ψ}}_{m})) - {\hat{x}}_{k | k} {\hat{x}}_{k | k}^{T} .

These series are convergent by the alternating series test and so are well approximated by partial sums. We typically truncate each summation when subsequent terms contribute less than some fraction ε << 1 (typically 10⁻⁶) to the total.

The normalization factor L_k is the likelihood of observing n_k photons given all previous observations and the system parameters; that is, L_k ≡ Pr[n_k|{n_j, j ≤ k − 1}; D, μ] (parameter conditioning previously suppressed for brevity), so the log-likelihood of the entire trajectory is

\ln (Λ) = \sum_{k} \ln (L_{k}) .

Maximum likelihood parameter estimates are found by gradient ascent of this function. Simulations testing the performance of this maximum likelihood estimator are in section 6, below. Application to experimental data is in section 7 and in [13].

4. Kalman filter

Equations (12-14) can be evaluated quickly enough for maximum likelihood parameter estimation in post-processing, but are still too complicated to be computed in real-time (i.e. during a 3.1 μs sampling bin) on current hardware. To simplify further the calculation, we map the ADF onto a Kalman filter. A Kalman filter can be thought of as a special type of assumed-density filter, in which the measurement distribution [Eq. (2)], the prediction equation [Eq. (3)], and the estimate distributions f_a_|_b are all taken to be Gaussian.

To map the ADF onto a Kalman filter, we return to Eq. (10) and restrict ourselves to scenarios in which the signal-to-background ratio is high (B << S) and the scan rate is fast relative to the rate of photon detection (S << 1). In this case, the dominant term in the series is the one with m = n_k and i = 0. So, the update step [Eqs. (13,14)] becomes

\begin{matrix} {\hat{x}}_{k | k} = {({\hat{Σ}}_{k | k - 1}^{- 1} + n_{k} W^{- 1})}^{- 1} ({\hat{Σ}}_{k | k - 1}^{- 1} {\hat{x}}_{k | k - 1} + n_{k} W^{- 1} c_{k}) \\ {\hat{Σ}}_{k | k} = {({\hat{Σ}}_{k | k - 1}^{- 1} + n_{k} W^{- 1})}^{- 1} . \end{matrix}

These equations, together with Eq. (8) allow propagation of a Gaussian likelihood function for the particle’s position using five variables (two for the position, three for the covariance), and are an instance of the well-known Kalman filter [27,28], with associated likelihood function

L_{k}^{k a l} = \frac{S^{n_{k}}}{n_{k}!} \frac{{| {\hat{Ψ}}_{n_{k}} |}^{\frac{1}{2}}}{{| {\hat{Σ}}_{k | k - 1} |}^{\frac{1}{2}}} \exp (- \frac{1}{2} n_{k} {({\hat{x}}_{k | k - 1} - c_{k})}^{T} {(n_{k} {\hat{Σ}}_{k | k - 1} + W)}^{- 1} ({\hat{x}}_{k | k - 1} - c_{k})) .

In the common case that W is diagonal (any beam asymmetry is along the same axes that the feedback is applied), the Kalman covariance is also diagonal and the cross-covariance term can be dropped. If the beam profile is circularly symmetric, then W = w²I for some w, and the variance in the x- and y-dimensions is always equal and can be denoted ${\hat{p}}_{a | b}$ . This yields the simple expressions

\begin{matrix} {\hat{x}}_{k | k} = \frac{w^{2} {\hat{x}}_{k | k - 1} + n_{k} {\hat{p}}_{k | k - 1} c_{k}}{w^{2} + n_{k} {\hat{p}}_{k | k - 1}} \\ {\hat{p}}_{k | k} = \frac{w^{2} {\hat{p}}_{k | k - 1}}{w^{2} + n_{k} {\hat{p}}_{k | k - 1}} . \end{matrix}

The prediction step [Eq. (8)] becomes

\begin{matrix} {\hat{x}}_{k | k - 1} = {\hat{x}}_{k - 1 | k - 1} + μ Δ t E_{k - 1} \\ {\hat{p}}_{k | k - 1} = {\hat{p}}_{k - 1 | k - 1} + 2 D Δ t . \end{matrix}

These four equations enable high-speed two-dimensional tracking of a particle using just three variables.

One concern with the tracking scheme outlined here is its requirement for prior knowledge of parameter values. Four parameters are required: the laser spot locations (c_k), the laser spot covariance (W or w²), the diffusion constant (D), and the electrokinetic mobility (μ). The first two can be measured by scanning a point source such as a surface-immobilized fluorescent bead through the laser during a scan. The population-average diffusion constant and electrokinetic mobility can be measured independently using multi-spot fluorescence correlation spectroscopy (FCS) in the presence of an applied electric field [29]. Alternatively, guesses can be made (or trapping can be attempted using a variety of parameter values), and maximum likelihood estimates can be calculated from trapping data and used to inform later experiments. Recently, a method for adaptive real-time parameter tuning was experimentally realized [20], in which trapping parameters were iteratively adjusted to “whiten” the autocorrelation of the innovation sequence, taking advantage of the fact that parameter errors give rise to correlations in the innovation (the displacement of each newly measured position from the prior estimate).

The ADF and the Kalman filter differ only in the update step. The Kalman update equation [Eq. (16)] is a simple linear combination of the prior estimate mean and the new measurement. The relative weighting is tuned by the relative confidences in the new measurement and in the prior estimate, but does not itself depend on the innovation. The ADF update also adjusts the prior estimate mean along the direction of each new measurement, but its relative weighting is itself a complicated nonlinear function of the innovation. This allows the ADF to better account for the true shape of the measurement distribution.

The assumption that B << S means that the Kalman filter does not distinguish between signal and background photons, so spurious counts significantly degrade the estimate. The assumption that S << 1 introduces a subtle bias, most noticeable when no photons are observed during a bin (n_k = 0). Absence of photons suggests that the particle is not in the laser spot, and so is more likely to be elsewhere. Analogously, if one were to peer through a telescope and see darkness, one could declare that the moon was not where the telescope was pointed. The ADF update step [Eq. (13)] accounts for this effect, while the Kalman update step [Eq. (16)] does not. Consequently, the Kalman estimator is slightly biased towards the most recent bin, leading to a modulation of the estimated particle position at the frequency of the laser scan.

5. Fundamental constraints

The Kalman filter is simpler and faster than the ADF, but has greater tracking error. In section 6, we present simulations quantifying the difference in performance between the two filters; here, we derive analytical estimates of these differences.

In each update step, the ADF incorporates new information about the fluorophore position. The information gain is quantified by the Fisher information matrix of the updated estimate, the (a,b)^th element of which is

- 〈 \frac{\partial^{2} \ln (L_{k})}{\partial x^{(a)} \partial x^{(b)}} 〉 = - 〈 \frac{\partial^{2} \ln (P r [n_{k} | x_{k}])}{\partial x^{(a)} \partial x^{(b)}} 〉 - 〈 \frac{\partial^{2} \ln (f_{k | k - 1})}{\partial x^{(a)} \partial x^{(b)}} 〉

where angle brackets indicate the expectation (integral) over the prior probability distribution for the molecule’s position (f_k_|_k_–1). The last term is the Fisher information of the prior distribution, which, under the Gaussian assumption of the ADF, is simply

{\hat{Σ}}_{k | k - 1}^{- 1}

. Of greater interest is the first term, which is the average information gain during each time bin. Substitution of Eqs. (1,2) into this term and dropping the time bin subscripts k for clarity gives

F_{A D F}^{(a, b)} \equiv - 〈 \frac{\partial^{2} \ln (P r [n | x])}{\partial x^{(a)} \partial x^{(b)}} 〉 = 〈 \frac{1}{Γ (x)} \frac{\partial Γ (x)}{\partial x^{(a)}} \frac{\partial Γ (x)}{\partial x^{(b)}} 〉

where Γ ≡ γΔt.

The overall average information gain for an arbitrary bin is this expression averaged over the laser scan positions c_k. If the trapping region is well sampled by the scan pattern (the spots are spaced regularly, the distance between them is less than or comparable to the width of the beam, and the overall size of the pattern is much larger than the spread of the beam width around the particle), then the sum over the laser scan positions is well approximated by an integral. Under these conditions, the Fisher information matrix for a two-dimensional scan is approximately

\bar{F_{A D F}} \approx W^{- 1} (1 + \frac{B}{S} L i_{2} (- \frac{S}{B})) \bar{Γ_{S}}

where over-bars indicate averaging over the bins k of the scan pattern, Li₂ is the polylogarithm function of order 2, and

\bar{Γ_{S}}

is the scan-averaged expected number of signal photons per time bin. (Under the given assumptions,

\bar{Γ_{S}} \approx \frac{2 π S {| W |}^{\frac{1}{2}}}{L Δ c^{2}}

, where L is the number of spots in the scan pattern and Δc is the spacing between scan positions.) The inverse of the Fisher information is the Cramér-Rao lower bound on the measurement error covariance associated with a particular scan pattern and signal-to-background ratio. In the absence of background counts, therefore, the measurement covariance is simply the beam covariance divided by the number of detected photons, akin to the localization precision of super-resolution techniques based on Gaussian fitting [4]. Unsurprisingly, the measurement covariance increases for B > 0.

Assuming that the ADF estimator is efficient (i.e., it makes optimal use of all available information), we approximate the steady-state value of the estimate covariance matrix (following the update step) according to

\bar{Σ^{- 1}} \approx \bar{F_{A D F}} + {(\bar{Σ} + 2 D Δ t I)}^{- 1}

where over-bars indicate averaging over the bins k of the scan pattern. All of the matrices commute, so we arrive at a final value for the average ADF tracking covariance of

\bar{Σ} \approx D Δ t ({(\frac{2}{D Δ t} \bar{F_{A D F}^{- 1}} + I)}^{\frac{1}{2}} - I)

where the square root is understood to indicate the unique positive definite symmetric matrix square root.

To determine the Kalman filter tracking error, we start by combining the update and prediction equations for its variance [Eqs. (8, 16)]

{\hat{Σ}}_{k | k} = {({({\hat{Σ}}_{k - 1 | k - 1} + 2 D Δ t I)}^{- 1} + n_{k} W^{- 1})}^{- 1} .

Neglecting fluctuations due to stochastic photon arrivals, this equation achieves the steady-state value

\bar{{\hat{Σ}}_{k | k}} \approx D Δ t ({(\frac{2}{D Δ t} \frac{W}{\bar{Γ}} + I)}^{\frac{1}{2}} - I)

where

\bar{Γ}

is the average number of photons detected per bin. Because the derivation of the Kalman filter assumed the absence of background, this equation is only accurate for B = 0. In fact, Eq. (26) reproduces the background-free value of the ADF tracking error [Eq. (24)], indicating that the filters’ tracking performance differs appreciably only when the background is significant.

Introduction of background photons broadens the spread of the detected photons around the molecule’s true position, increasing the effective beam waist and degrading the Kalman filter’s performance. The magnitude of the broadening is difficult to calculate because it depends on the distribution of the particle in the trap. The ADF tracking error estimate [Eq. (24)] is a lower bound on the Kalman tracking error in the presence of background; we use simulations (section 6) to quantify performance more precisely.

In the limit of fast scanning (∆t << (2w²/Dγ)^½), the background-free tracking error of the ADF and the Kalman filter [Eq. (26)] simplifies to

\bar{Σ} \approx \sqrt{\frac{2 D W}{\bar{γ_{s}}}} .

This result can be interpreted as the product of the width of the laser spot and the mean distance the particle diffuses between photon detection events, and reproduces the tracking error formula for a continuous-time linear filtration scheme based on lock-in detection [17].

During the update step of either the ADF or Kalman filter, the prior estimate mean is multiplied by a weighting factor κ that reflects the confidence that it retains validity. In the case of the Kalman filter with no background and circular symmetry [Eq. (18)], the weighting factor is

κ = \frac{w^{2}}{w^{2} + n_{k} {\hat{p}}_{k | k - 1}} .

Using the approximate average tracking error from Eq. (26), the average weighting factor is

\bar{κ} \approx \frac{w^{2}}{w^{2} + \bar{Γ} D Δ t ({(\frac{2}{D Δ t} \frac{w^{2}}{\bar{Γ}} + 1)}^{\frac{1}{2}} + 1)} .

The estimate constructed j bins ago will typically be weighted by κ^j. Recasting in terms of elapsed time T ≡ j∆t and setting ∆t → 0 gives

\bar{κ_{T}} \approx \exp (- \frac{\sqrt{2 D \bar{γ}}}{w} T) .

Thus, the timescale during which information from each photon retains utility is

τ \approx \frac{w}{\sqrt{2 D \bar{γ}}} .

The denominator in Eq. (31) represents the average velocity of the particle between photon detection events. This formula has also been derived in the continuous-time case [17].

To illustrate, Alexa 647 has a diffusion constant of 330 μm²/s [30], and we typically have $\bar{γ_{s}}$ ≈50 kHz and $\bar{γ_{b}}$ ≈9 kHz (corresponding to S ≈1.5 and B ≈0.03 for our scan pattern). Our 27-point scan pattern has spot spacing of 0.65 μm, beam width w ≈0.4 μm, and scan rate of 12 kHz. Equation (24) predicts a minimum tracking error of 315 nm. The timescale over which measurements retain utility [Eq. (31)] is approximately 70 μs, corresponding to a scan rate of 14 kHz. Our choice of 12 kHz approaches this limit, and trap performance shows little change when the scan rate is increased.

6. Simulation: trapping single fluorophores in solution

To demonstrate application to trapping single fluorophores in solution, we simulated trapping experiments at a variety of parameter values and calculated maximum likelihood parameter estimates from the simulated data. All calculations were performed in MATLAB. The default transport coefficients were D = 100 μm²/s, μ = −1000 μm²/(V s). The feedback electric field strength was capped at 3 V/μm. The expected photon count rate for a molecule in the center of the trap was set to 100 kHz, with a background count rate of 25 kHz. The laser scan matched our experimental setup, with 27 spots of width w = 0.4 μm on a hexagonal lattice of 0.65 μm spacing, scanned with a per-spot dwell time ∆t = 3.1 μs. The Kalman filter tracking and feedback algorithm also matched the experimental implementation as closely as possible. All parameters were held at these values unless otherwise noted. For each parameter set, we simulated one hundred individual molecular traces, each of length 0.1 s. Molecules were considered successfully “trapped” if their maximum displacement from the trap center never exceeded 3.4 μm.

We compared the tracking accuracy of the Kalman filter, ADF, and a full recursive Bayesian estimator [Eqs. (5,6)] on each simulated molecular trajectories. The model parameters were set to match the true simulated values. The full Bayesian estimator, which was propagated on a 128 × 128 grid with 38 nm spacing (encompassing the entire scan pattern), was indistinguishable from the ADF algorithm, except for a few cases in which its performance was hindered by the grid spacing. All of the algorithms were programmed in C using the MATLAB mex interface; the ADF typically ran more than 1500 times faster than the full grid estimator, and the Kalman filter ran 9 times faster than the ADF.

Figure 2 summarizes the simulation results for varying diffusion coefficient or electrokinetic mobility; Fig. 3 displays the same information as a function of fluorophore brightness and signal-to-background ratio. Stable trapping is achieved for a wide range of simulated parameters, but is hindered primarily when the diffusion coefficient is high or the photon count rate and signal-to-background ratio are low. As expected, the RMS tracking error of the Kalman filter is uniformly greater than that of the ADF, but the difference is typically small unless the background count rate is significant. Equation (24) offers a good approximation of the ADF tracking error, with significant deviations only when molecules are not stably trapped. The derivation of Eq. (24) assumes that the region of confinement is well sampled by the scan pattern; a significant probability of escape implies that the molecules spend significant time at the edges of the scan pattern, violating this assumption and increasing the tracking error above the analytical estimate. The Kalman filter performs well for tracking when S > 1 (such as when the molecule’s count rate is 100 kHz), in accordance with Eq. (26), despite the assumption S << 1 in the derivation. The primary effect of violation of the S << 1 approximation is to bias the Kalman filter towards the most recent bin (see section 4), and although this bias worsens as S increases, it is more than compensated by the additional information provided by the larger number of photons.

Fig. 2 Simulated trapping and parameter fitting as a function of diffusion coefficient (left) or electrokinetic mobility (right). Top, fraction of molecules remaining within 3.4 μm of the trap center for the 100 ms simulation. Center, RMS tracking error of the Kalman filter and ADF algorithms. The dashed curve is the analytical prediction for trapping performance [Eq. (24)]. Bottom, maximum likelihood parameter estimation using the ADF algorithm. Parameter estimation using the Kalman filter did not converge for any of these cases. Error bars are ± one standard deviation about the mean of the relative difference between the estimated and true parameter values.

Download Full Size | PDF

Fig. 3 Simulated trapping and parameter fitting as a function of photon count rate (not including background photons) and signal-to-background ratio (SBR). Plots are as in Fig. 2, except that the bottom plots all show diffusion coefficient fitting, and fits were attempted using both the ADF algorithm and the Kalman filter. The parameter fitting results for the two algorithms are offset horizontally for clarity.

Download Full Size | PDF

Maximum likelihood estimates of diffusion coefficient and electrokinetic mobility are plotted in the bottom of Fig. 2. For high values of the diffusion coefficient, the relative precision of the fit is improved by the increased magnitude of the diffusive motion. Similarly, the relative precision of the fit of the mobility is improved for higher electrokinetic mobility because the motion due to the applied feedback becomes greater relative to diffusion.

Figure 3 (bottom) shows maximum likelihood estimates of the diffusion coefficient as a function of the molecular photon count rate, at three values of the signal-to-background ratio (SBR, $\bar{γ_{s}} / \bar{γ_{b}}$ ). Unsurprisingly, the estimate improves as either the signal rate or the SBR increase. The strong precision and accuracy of the estimator even when the SBR is 1 highlights the background rejection afforded by the ADF.

We compared our parameter estimation results using the ADF to those derived using the Kalman filter and its associated likelihood function [Eq. (17)]. For many choices of parameters, the Kalman filter failed to find a reasonable fit (defined by maximum likelihood diffusion coefficient ≥1500 μm²/s); the successful cases are plotted in Fig. 3. (None of the parameter sets plotted in Fig. 2 resulted in a successful Kalman fit.) Consistent with the derivation in section 4, the Kalman filter performs well in cases where B << S << 1; indeed, its performance is nearly indistinguishable from that of the ADF in the case of zero background and $\bar{γ_{s}}$ < 30 kHz (S < 0.9). Because the Kalman filter treats all photons as real, it vastly overestimates molecular diffusion in the presence of background and hence is extremely sensitive to SBR. Unlike the ADF, the Kalman filter fit worsens as the signal strength increases, to the point that it fails when the photon count rate is 100 kHz (S ≈3), even in the absence of background.

Violations of the assumptions B << S << 1 cause the Kalman filter to report that the particle is moving more than it actually is, due to the treatment of background photons as real when B > 0, and due to the bias towards the most recent bin when S > 0. When the Kalman filter is supplied with a diffusion coefficient not too far from the true value (e.g. when it is used for tracking), this extraneous motion is largely suppressed by the filter and the tracking error remains low. However, when the Kalman filter is used for parameter fitting, the model inflates the diffusion coefficient to explain the excess motion. The ADF does not encounter this problem because it can increase the parameters B or S to account for the detected photons.

Finally, we calculated the tracking error associated with both filters if the model parameters were inaccurate (Fig. 4 ), at three values of the SBR. Even fairly large discrepancies in the diffusion coefficient and electrokinetic mobility result in only modest increases in the tracking error. When B > 0, the Kalman filter’s tracking error is minimized when the model diffusion coefficient is set below its true value (or, equivalently, when the beam width is set to a larger value); this is because the introduction of background broadens the effective beam width, an effect that is compensated in the ADF but not in the Kalman filter.

Fig. 4 Tracking error as a function of filter parameter mismatch, for three values of signal-to-background ratio (SBR) and both algorithms. Fluorophore photon emission rate is 100 kHz for all plots. The ADF performs best with parameter values matching the true values. For nonzero background count rates, the Kalman filter performs better when its assumed diffusion coefficient is less than the true value.

Download Full Size | PDF

7. Experiment: trapping single molecules of Alexa 647

We built an ABEL trap that used the Kalman filter derived in section 4 for real-time particle tracking and feedback [13]. The improved tracking precision of the Kalman filter (see section 6) enabled the trap to hold individual molecules of Alexa 647 fluorophore for an average of 1 s per trapping event. Additionally, the ADF algorithm derived in section 3 allowed us to find the maximum likelihood diffusion coefficient and electrokinetic mobility of each molecule [Fig. 5 ]. The diffusion coefficient averaged over n = 137 trapping events was D = 348 ± 2 μm²/s (s.e.m.), in reasonable agreement with the value 325 μm²/s we obtained from a different data set [13] and the value 330 μm²/s measured by FCS [30].

Fig. 5 Trapping of Alexa 647 molecules in an anti-Brownian electrokinetic trap. (a) Molecular entry or exit induced large step-like changes in fluorescence intensity. Intensity trace binned to 10 ms. (b) The diffusion coefficient D and electrokinetic mobility μ of each trapped molecule was extracted using the ADF algorithm. The population-average diffusion coefficient was 348 ± 2 μm²/s (s.e.m.); the electrokinetic mobility revealed the presence of two subpopulations, most likely differing in charge.

Download Full Size | PDF

The mean electrokinetic mobility was μ = −5.0 × 10³ μm²/(V s). For a free particle of size much smaller than the Debye screening length, the Einstein-Smoluchowski relation predicts an electrophoretic mobility μ = qD/k_BT, where q is the charge of the particle. From our measurement of D, this relation predicts μ = −1.4 × 10⁴ μm²/(V s) assuming each molecule bears a single negative charge, significantly different from our data. However, the Einstein-Smoluchowski formula does not take into account electroosmotic flow, the fluid motion induced by application of an electric field to water in a small capillary. This flow velocity depends on details of chemical state of the channel walls, and is expected to be counter to the direction of the electrophoretic motion in our experiments.

Unexpectedly, the distribution of electrokinetic mobilities was bimodal; the clusters most likely represent distinct charge states of the dye molecule, possibly corresponding to structural isomers. These results demonstrate that the tracking and feedback strategy derived here enable trapping of single fluorescent dyes in water at room temperature. The algorithm reports molecule-by-molecule transport coefficients – information that would be difficult or impossible to measure by other means. Additional experimental results are presented in [13].

8. Conclusion

The mathematical framework presented here constitutes a near-optimal approach to track a fluorescent particle undergoing Brownian diffusion plus linear drift in the face of background photons and shot noise. The successive approximations we derive enable real-time particle tracking and offline maximum likelihood estimation of transport parameters, with well-controlled estimates of the uncertainties. This work was motivated by the ABEL trap; our simulation and experimental results demonstrate the algorithms’ efficacy in this context. The same treatment also applies to other anti-Brownian traps, such as those functioning by stage motion.

The recursive Bayesian framework derived here is adaptable to other measurement modalities or models of particle motion via modification of the update or prediction equations, respectively. A simple adjustment to the prediction step accounts for the case that the applied feedback is itself stochastic [13]. The case of camera tracking can be handled by treating each pixel as a “laser spot” and applying the update steps corresponding to each of the pixels in a frame simultaneously (instead of applying a prediction step in between, as in the case of scanned laser tracking). Non-Brownian motion is less simple to handle; the treatment offered here could be used as-is, with the understanding that it represents an inexact approximation of the motion as memoryless and Gaussian, or a more specific model could be constructed.

The traditional approach of treating every frame in a particle tracking experiment independently is mathematically equivalent to a model of molecular motion in which the particle is equally likely to move anywhere regardless of its previous position, so information from recent frames is ignored during the analysis of the current frame. This approach is statistically inefficient if a model of the motion is known. The algorithms we have derived demonstrate how to apply an explicit model of particle motion to get the most from the data collected in a tracking experiment.

Acknowledgements

This work was supported in part by National Science Foundation (NSF) Grant CHE-0910824; the Materials Research Science and Engineering Center (MRSEC) of Harvard University under NSF grant DMR-0820484; and the Harvard Faculty of Arts and Sciences (FAS) Center for Nanoscale Systems, a member of the National Nanotechnology Infrastructure Network, supported by NSF award ECS-0335765. Computations were run on the Odyssey cluster supported by the FAS Research Computing Group.

References and links

1. M. J. Saxton and K. Jacobson, “Single-particle tracking: applications to membrane dynamics,” Annu. Rev. Biophys. Biomol. Struct. 26(1), 373–399 (1997). [CrossRef] [PubMed]

2. A. D. Douglass and R. D. Vale, “Single-molecule microscopy reveals plasma membrane microdomains created by protein-protein networks that exclude or trap signaling molecules in T cells,” Cell 121(6), 937–950 (2005). [CrossRef] [PubMed]

3. I. Chung, R. Akita, R. Vandlen, D. Toomre, J. Schlessinger, and I. Mellman, “Spatial control of EGF receptor activation by reversible dimerization on living cells,” Nature 464(7289), 783–787 (2010). [CrossRef] [PubMed]

4. A. Yildiz, J. N. Forkey, S. A. McKinney, T. Ha, Y. E. Goldman, and P. R. Selvin, “Myosin V walks hand-over-hand: single fluorophore imaging with 1.5-nm localization,” Science 300(5628), 2061–2065 (2003). [CrossRef] [PubMed]

5. M. A. Thompson, J. M. Casolari, M. Badieirostami, P. O. Brown, and W. E. Moerner, “Three-dimensional tracking of single mRNA particles in Saccharomyces cerevisiae using a double-helix point spread function,” Proc. Natl. Acad. Sci. U.S.A. 107(42), 17864–17871 (2010). [CrossRef] [PubMed]

6. N. P. Wells, G. A. Lessard, P. M. Goodwin, M. E. Phipps, P. J. Cutler, D. S. Lidke, B. S. Wilson, and J. H. Werner, “Time-resolved three-dimensional molecular tracking in live cells,” Nano Lett. 10(11), 4732–4737 (2010). [CrossRef] [PubMed]

7. A. P. Fields and A. E. Cohen, “Anti-Brownian traps for studies on single molecules,” Methods Enzymol. 475, 149–174 (2010). [CrossRef] [PubMed]

8. A. E. Cohen and W. E. Moerner, “Principal-components analysis of shape fluctuations of single DNA molecules,” Proc. Natl. Acad. Sci. U.S.A. 104(31), 12622–12627 (2007). [CrossRef] [PubMed]

9. Y. Jiang, Q. Wang, A. E. Cohen, N. Douglas, J. Frydman, and W. E. Moerner, “Hardware-based anti-Brownian electrokinetic trap (ABEL trap) for single molecules: control loop simulations and application to ATP binding stoichiometry in multi-subunit enzymes,” Proc. Soc. Photo Opt. Instrum. Eng. 7038, 1–12 (2008). [PubMed]

10. R. H. Goldsmith and W. E. Moerner, “Watching conformational- and photodynamics of single fluorescent proteins in solution,” Nat. Chem. 2(3), 179–186 (2010). [CrossRef] [PubMed]

11. H. Cang, D. Montiel, C. S. Xu, and H. Yang, “Observation of spectral anisotropy of gold nanoparticles,” J. Chem. Phys. 129(4), 044503 (2008). [CrossRef] [PubMed]

12. J. Enderlein, “Tracking of fluorescent molecules diffusing within membranes,” Appl. Phys. B 71(5), 773–777 (2000). [CrossRef]

13. A. P. Fields and A. E. Cohen, “Electrokinetic trapping at the one nanometer limit,” Proc. Natl. Acad. Sci. U.S.A. 108(22), 8937–8942 (2011). [CrossRef] [PubMed]

14. A. H. Jazwinski, Stochastic Processes and Filtering Theory (Academic Press, 1970).

15. K. McHale, A. J. Berglund, and H. Mabuchi, “Bayesian estimation for species identification in single-molecule fluorescence microscopy,” Biophys. J. 86(6), 3409–3422 (2004). [CrossRef] [PubMed]

16. A. J. Berglund, K. McHale, and H. Mabuchi, “Fluctuations in closed-loop fluorescent particle tracking,” Opt. Express 15(12), 7752–7773 (2007). [CrossRef] [PubMed]

17. A. J. Berglund and H. Mabuchi, “Performance bounds on single-particle tracking by fluorescence modulation,” Appl. Phys. B 83(1), 127–133 (2006). [CrossRef]

18. A. J. Berglund and H. Mabuchi, “Tracking-FCS: Fluorescence correlation spectroscopy of individual particles,” Opt. Express 13(20), 8069–8082 (2005). [CrossRef] [PubMed]

19. A. E. Cohen and W. E. Moerner, “Method for trapping and manipulating nanoscale objects in solution,” Appl. Phys. Lett. 86(9), 093109 (2005). [CrossRef]

20. Q. Wang and W. E. Moerner, “An adaptive anti-Brownian electrokinetic trap with real-time information on single-molecule diffusivity and mobility,” ACS Nano 5(7), 5792–5799 (2011). [CrossRef] [PubMed]

21. K. I. Mortensen, L. S. Churchman, J. A. Spudich, and H. Flyvbjerg, “Optimized localization analysis for single-molecule tracking and super-resolution microscopy,” Nat. Methods 7(5), 377–381 (2010). [CrossRef] [PubMed]

22. B. Zhang, J. Zerubia, and J. C. Olivo-Marin, “Gaussian approximations of fluorescence microscope point-spread function models,” Appl. Opt. 46(10), 1819–1829 (2007). [CrossRef] [PubMed]

23. M. S. Arulampalam, S. Maskell, N. Gordon, and T. Clapp, “A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking,” IEEE Trans. Signal Process. 50(2), 174–188 (2002). [CrossRef]

24. T. P. Minka, “A family of algorithms for approximate Bayesian inference,” Ph.D. thesis, Massachusetts Institute of Technology (2001). http://research.microsoft.com/en-us/um/people/minka/papers/ep/minka-thesis.pdf.

25. P. S. Maybeck, Stochastic models, estimation and control (Academic press, 1979).

26. H. W. Sorenson and D. L. Alspach, “Recursive Bayesian estimation using Gaussian sums,” Automatica 7(4), 465–479 (1971). [CrossRef]

27. R. E. Kalman, “A new approach to linear filtering and prediction problems,” J. Basic Eng. Trans. ASME 82(1), 35–45 (1960). [CrossRef]

28. G. Welch and G. Bishop, “An introduction to the Kalman filter,” University of North Carolina at Chapel Hill technical report TR 95–041 (2006). http://www.cs.unc.edu/~welch/media/pdf/kalman_intro.pdf.

29. M. Brinkmeier, K. Dörre, J. Stephan, and M. Eigen, “Two-beam cross-correlation: a method to characterize transport phenomena in micrometer-sized structures,” Anal. Chem. 71(3), 609–616 (1999). [CrossRef] [PubMed]

30. P. Kapusta, “Absolute diffusion coefficients: compilation of reference data for FCS calibration,” http://www.picoquant.com/technotes/appnote_diffusion_coefficients.pdf.

Optimal tracking of a Brownian particle

Abstract

1. Introduction

2. General treatment

3. Gaussian assumed-density filter (ADF)

4. Kalman filter

5. Fundamental constraints

6. Simulation: trapping single fluorophores in solution

7. Experiment: trapping single molecules of Alexa 647

8. Conclusion

Acknowledgements

References and links

Cited By

Figures (5)

Equations (31)

Optics Express