1. Introduction

Due to the invention of specialized instruments with nanometer-scale resolution, such as the atomic force microscope and the scanning electron microscope, our knowledge of biological features at the nanometer scale has greatly improved in recent times. These instruments, however, are limited by their relatively slow imaging speeds and by other restrictions in terms of applicability to imaging live cells. Because of the manner by which they acquire their images (via point-by-point measurements through raster scanning) their temporal resolutions are far slower than inherently parallel instruments, such as widefield light microscopes. The spatial resolution of light, unfortunately, remains limited by diffraction; even cutting edge instruments, like stimulated emission depletion (STED) microscopes, have yet to reach single nanometer spatial resolution in three-dimensions with millisecond temporal resolution [1]. The class of techniques broadly labeled as single particle tracking, however, has shown its potential for achieving excellent temporal and spatial resolution.

In a typical particle tracking experiment, one or more sub-diffraction limit sized objects of interest, whether they be viruses, proteins, or other biomolecular macromolecules, are labeled by individually binding small fluorescent particles to them. Commonly used labels include proteins, such as green fluorescent protein (GFP), dyes, such as rhodamine, or crystals, such as quantum dots. By selectively labeling a feature with a fluorophore, its dynamic characteristics may be inferred by measuring the photon emission of the label over some temporal period and spatial area. For example, by labelling an individual leg of a myosin V motor with a fluorophore and observing the fluorophore’s emission as the motor traversed an actin filament, both qualitative and quantitative information regarding its walking behavior were derived [2]. The effectiveness of this method was also demonstrated in the context of tracking phospholipids in rat kidney fibroblasts [3] and of the influenza virus during infection [4] while its utility and impact is illustrated by several recent review articles [5–7].

Due to its simplicity, the prevailing method for tracking particles involves passively measuring the emitted fluorescence as the particles move during the experiment; this approach, however, is not without challenges. Since the particles are not actively contained within the observer’s field of view, the sensor of choice must consist of many pixels to compensate for the uncertainty of the particle’s potential range of motion. One ramification for this method is the amount of time required to read each individual pixel value, resulting in limited frame rates and causing motion blur when the dynamics of the feature are particularly fast relative to the imaging rate. This motion blur obscures the image and contributes additional localization uncertainty. In addition, because of the axial symmetry of the point spread function (PSF) these methods cannot natively track particles in three dimensions without augmenting the instrument with specialized optics that distort the image with a deterministic axial dependence [8–11] or with multiple cameras to image multiple planes simultaneously [12]. Many of these methods have limited axial range or limit the measured intensity on any single detector by spreading the photons out over multiple detectors.

There is an alternative class of techniques which make a different set of tradeoffs, choosing to limit the field of view to a single particle of interest and then to follow it actively as it moves throughout its environment. To accomplish this, most active methods localize the particle in real-time and use the resulting estimated position in a feedback algorithm that reduces the tracking error. These methods are typically implemented using confocal or two-photon imaging modalities and can be broadly categorized as either orbital, in which the excitation path is modified to allow for rapid scanning of the excitation volume along a given a given closed path (typically circular), or multi-point, in which the detection path is modified to incorporate multiple detectors at different points and different planes. Under the orbital scheme, fluctuations in the measured intensity along the closed path are used both to estimate the particle location and to center the scanning circle on the (estimated) particle position. In one implementation, axial information was provided by sequentially scanning circles in two different planes [13, 14]. A similar technique was introduced by one of the authors but rather than rapid scanning of the excitation volume, it was moved through a finite constellation of points. The technique was designed to be optimal in the sense that the constellation provided a minimum variance unbiased estimate of the particle location when localized with the fluoroBancroft algorithm [15–17]. Tracking of particles diffusing at rates on the order of 0.1 μm²/s have been reported [14]. In general, orbital methods are limited in their ability to track quickly moving particles due to the time required to collect enough information to localize the particle in real-time as part of the control loop. While the former method was limited by the maximum angular speed of the circle, the time to switch between axial planes, and the need to average measurements over multiple traversals of the circular path to improve the signal-to-noise ratio, the latter was limited by the time it took to travel to each constellation point. Multi-point approaches avoid the need for scanning beams entirely by introducing multiple pinholes and detectors into the output path [19,20]. This allows the PSF to be sampled at multiple points simultaneously and tracking of particles diffusing at rates on the order of 1 μm²/s have been reported. Clearly, if more focal volumes were added so that more measurements could be performed in parallel, then the quality of tracking could be improved. Consequently, a high-speed dual-beam instrument which combined both the orbital approach on the excitation path and multi-point detection on the output path was built which demonstrated the effectiveness of the active paradigm by tracking particles with diffusion coefficients up to 20 μm²/s [18]. However, the spatial accuracy of this system on such fast particles was reported to be 352 nm in the lateral directions and 272 nm along the optical axis.

While both orbital and multi-point techniques are quite promising in their ability to track quickly moving particles, they are significantly more complicated and expensive than standard or home-built confocal or two-photon instruments. A more optically straightforward approach has been developed that uses multi focal-plane imaging in a wide-field setting combined with active laser positioning for tracking of a single particle [22]. Since it is based on a wide-field modality, however, it still requires on-line position estimation and is subject to issues of motion blur and similar challenges common to this setup.

In this work, we present an active method for tracking fluorescent particles that can be implemented on conventional confocal and multi-photon microscopes. Such instruments are both relatively common and easy to set up in the lab [21]. The only requirements of the proposed method are that the instrument take point-by-point measurements of intensity, that the intensity pattern (i.e., the PSF) has an isolated maximum, and that the distance between the measurement point and the particle must have three independent degrees of freedom to allow for three-dimensional Cartesian actuation. A sample-actuated epifluorescent confocal microscope is one potential instrument since it moves the sample relative to the excitation-detection volume and measures the intensity at a specified point in space. A laser-scanning confocal microscope may also be used for the same reasons except that it moves the excitation-detection volume instead of the sample. The method is able to track fluorescent particles in real-time by employing an algorithm which “seeks” the extremum of the PSF as intensity measurements are acquired. In essence, the algorithm explores the local structure of the PSF and implicitly ascends its gradient to move towards the peak intensity. When locked on to a particle, the algorithm shares some similarities with the orbital approach in that the steady state motion is spherical-like corresponding to three distinct angular velocities. The radius of the exploratory sphere R and the convergence gain K_p are specified by the experimenter; both of these parameters affect the tracking performance and allow the algorithm to be experimentally tuned for a given microscope and experiment configuration. As the speed of motion of the particle increases, the algorithm naturally transitions away from the orbital pattern to a chase mode where the trajectory of the detection volume begins to resemble the characteristics of the motion of the particle being tracked. The algorithm is simple to implement and requires no modification to the optical setup other than the ability to control the actuators moving the relative position of the measurement point and the particle.

The algorithm is described in Sec. 2. It is a modification of an earlier algorithm of the first author for tracking single fluorescent particles [24]. This previous work was studied through simulation and shown to be able to track particles diffusing at rates on the order of 1 μm²/s at moderate fluorescence intensities with relatively slow piezo-actuated state-scanning, and at rates on the order of 90 μm²/s at extremely high fluorescence intensities and very high bandwidth, acousto-optic actuated beam scanning. The modification was presented in an earlier work by the authors [25] where it was described in the context of robotic potential field exploration. The modification to the original algorithm allows for a simpler description of the steady state behavior of the algorithm as well as better tracking performance in terms of stability. In this article, we solely focus on how the algorithm is used for tracking fluorescent particles. In this context, the role of the potential field is played by the PSF and the role of the robot is played by the actuation of the distance between the focal volume and the particle center. Proof of the algorithm’s convergence to a spherical orbit about a stationary extremum is given in [25].

The method is reactive, that is the focal volume is moved in response only to the measured intensity. There is no explicit estimation or measurement of the position of the tracked particle and thus a potential benefit of using the proposed tracking method is that it yields raw measurements. The process by which the particle is localized may be performed offline without any real-time constraints. In this work, we employ the SMC-EM algorithm which uses the expectation maximization (EM) algorithm in conjunction with sequential Monte Carlo (SMC) methods to simultaneously localize the particle and infer model-based motion parameters (e.g., diffusion coefficients). The SMC-EM method yields an approximate maximum likelihood estimate of the motion parameters in addition to discrete approximations of the time-dependent posterior probability densities of the particle position given all the measurements. The estimation technique is applicable to any single particle tracking method, including both passive (wide-field) and active schemes, through appropriate modeling of the measurement. In Sec. 3 we describe its application in the context of the confocal measurements used in our reactive tracker; further details on SMC-EM, including validation of estimation accuracy can be found in [23].

2. Tracking algorithm

We assume the location of the particle, relative to some fixed reference frame and denoted by x_p(t) = (x_p(t), y_p(t), z_p(t)), moves in time with unknown dynamics, and that the focal volume location, relative to the same reference frame and denoted by x_s(t) = (x_s(t), y_s(t), z_s(t)), moves according to the continuous-time extremum-seeking dynamics

In practice, the intensity measurement is formed by accumulating photons over a fixed time interval. Thus, the controller in Eq. (1) must be implemented in discrete time. We assume the accumulation period is Δt, yielding the discrete-time intensity measurement I_k. Additionally, we assume the commanded focal volume position is updated synchronously with the intensity measurement. By numerically integrating with a first-order Euler method, the commanded focal volume position is governed by

Assuming a relatively fast update frequency relative to the largest of the three orbital frequencies (i.e., ω₁ + ω₂), this discretization yields an accurate representation of the continuous time dynamics. The four selectable parameters maintain their original characteristics in this case.

Convergence to the steady state deviation from the particle center, given by Eq. (2), is only true for a stationary particle and for a radial PSF. We note that any deviation from these two assumptions will result in the addition of a time-varying disturbance to Eq. (2). Despite the disturbances, the algorithm will still move to the neighborhood of the maximum intensity and is thus robust against aberations in the PSF. When the source particle is moving, the algorithm attempts to reject any particle motion so that the deviation Eq.(2) is maintained; in fact, this is precisely why the algorithm is able to track moving particles. In general, the frequencies ω₁ and ω₂ should be taken as large as possible, subject to the bandwidth of the actuators, so that the sensing trajectory is (approximately) in its steady state for as broad a range of particle motion as possible. However, even when not in steady state, tracking can occur and the algorithm naturally converts to a “chase” mode for faster particles.

Because the algorithm as expressed in Eq. (3) and Eq. (4) is very simple to implement in a digital controller, it can be performed at extremely high sample rates (easily on the order of 10–100 kHz). As with any active scheme, its tracking capability is ultimately limited by the bandwidth of the actuators used and by the available photon rate in the experiment. In essence, so long as the rate of motion of the particle is such that it cannot move completely out of the detection volume of the microscope within one time step, the algorithm can react to the measurement and act to keep the particle in that volume.

The algorithm can be implemented on a variety of physical setups. The instrument used in this work consisted of an inverted epifluorescent microscope (Axiovert 200M, Zeiss) augmented with confocal optics. A 488 nm diode laser (FiberTec II, Blue Sky Research) was expanded and collimated to overfill the back aperture of the objective lens (C-Apochromat 63x/1.20 W Corr M27, Carl Zeiss). The collimated beam, after being reflected off the dichroic filter (HW625/30m, Chroma), was focused onto the sample by the objective lens. A three-axis piezoelectric nanopositioner (Nano-PDQ, Mad City Labs) with a 50 × 50 × 25 μm range of motion held the sample above the objective lens. The generated fluorescence was imaged by the objective and passed through the dichroic filter and tube lens where it was focused onto a 75 μm-diameter pinhole (P75S, Thorlabs) and avalanche photodiode (SPCM-AQR-14, Perkin Elmer). For diagnostics, a pellicle beam splitter (CM1-BP133, Thorlabs) split approximately one-third of the fluorescence onto a CCD camera (EXi Aqua, QImaging); the CCD was used for finding regions of interest and was not used during tracking. The control algorithm was implemented on the field programmable gate array (FPGA, Spartan 6, Xilinx) of a real-time embedded controller (CompactRIO 9076, National Instruments) which recorded the digital photodiode pulses and calculated commands for the nanopositioner according to the tracking algorithm. A diagram of the instrument is shown in Fig 1.

 

Fig. 1 Schematic illustrating the confocal microscope used in this work. A 488 nm diode laser is expanded and collimated to overfill the back aperture of the objective lens. The beam reflects off the dichroic and passes through the objective lens which focuses it onto the fluorescent sample. A three-axis piezoelectric nanopositioner displaces the sample relative to the fixed beam. The resulting fluorescence is imaged by the objective, passed through the dichroic, and focused onto a 75 μm-diameter pinhole and avalanche photodiode. To assist with diagnostics, a pellicle beam splitter divided one-third of the fluorescence onto a CCD camera; the CCD was used for finding regions of interest and was not used during tracking. The control algorithm was implemented on the FPGA of a real-time embedded controller which sampled the photodiode pulses and calculated commands for the nanopositioner.

Download Full Size | PDF

3. Data analysis: particle trajectory and motion model estimation

Recall that the output of a particle tracking experiment using the tracking method of Sec. 2 is a finite collection of intensity measurements at their respective measured locations - the tracker does not provide an estimate of the particle position. Notationally, the output is a collection ${\mathbf{x}}_{\mathbf{s},1:N}\triangleq \{{\mathbf{x}}_{\mathbf{s},1},{\mathbf{x}}_{\mathbf{s},2},\dots ,{\mathbf{x}}_{\mathbf{s},N}\}$ of N positions (where x_s,k denotes the measured focal volume position at discrete time index k) and a collection $I_{1:N}\triangleq \{I_1,I_2,\dots ,I_N\}$ of N corresponding measurements (where I_k denotes the intensity measurement at discrete time index k). As noted above, this allows for estimation of the particle trajectory and experimental parameters to occur offline using a variety of algorithms.

The estimation procedure has two primary aims. The first is to provide an estimate of the N positions ${\mathbf{x}}_{\mathbf{p},1:N}\triangleq \{{\mathbf{x}}_{\mathbf{p},1},{\mathbf{x}}_{\mathbf{p},2},\dots ,{\mathbf{x}}_{\mathbf{p},N}\}$ given the measurements (where x_p,k denotes the position of the particle at time step k.). Specifically, the estimation procedure provides an approximation of the posterior probability density of the particle positions given the measurements: p(x_p,1:N |I_1:N, x_s,1:N). From this density, any metric may be calculated to infer the position of the particle, including the mean and variance. The second aim is to infer model-based parameters that describe the experimental process (e.g., the motion of the particle or the measurement noise). The estimation procedure described here provides an approximate maximum likelihood (ML) estimate of all unknown parameters θ, where the ML estimate ${\widehat{\theta }}_{\text{ML}}$ is defined by

We assume the particle location x_p,k moves stochastically in discrete time according to the probability density

Since the particle location is unknown and cannot be measured directly, intensity measurements I_k ∈ ℝ are taken according to a different probability distribution

Lastly, we assume the initial location of the particle is distributed according to

In this work, we assume the particle undergoes a combination of anisotropic diffusion and directed motion. Thus, the probability distribution for the x axis is given by

We further assume that the observation I_k is taken at the focal volume position x_s,k and is corrupted by both shot noise and uniform background noise; for simplicity, we assume the noise in the measured focal volume position is negligible. Thus, the observation density is

The initial condition x_p,1 is assumed to be normally distributed with unknown mean μ = [μ_x, μ_y, μ_z]^T and variance $\sum =\text{diag}\left({\sigma }_x^2,{\sigma }_y^2,{\sigma }_z^2\right)$ such that

Thus, there are six parameters to be estimated pertaining to the initial condition model: μ_x, μ_y, μ_z, ${\sigma }_x^2$, ${\sigma }_y^2$, ${\sigma }_z^2$.

These probability densities above are of known form but are parametrized by unknown parameters θ. To obtain estimates of these parameters, we would ideally calculate the maximum likelihood (ML) estimate

To use the EM algorithm, we consider the joint likelihood of the particle locations and the measurements, namely p_θ(x_p,1:N, I_1:N). Instead of maximizing over the likelihood as was done in Eq. (12), the EM algorithm considers instead the quantity

Here, $Q(\theta ,{\widehat{\theta }}_e)$ is a weighted averaging of the logarithm of the joint density with respect to the particle location given some estimate of θ, namely ${\widehat{\theta }}_e$. As was originally presented in [27]), the simplification of Eq. (13) with respect to Eqs. (6)–(8) yields

It was shown in [28] that for each choice of ${\widehat{\theta }}_{e+1}$ such that $Q({\widehat{\theta }}_{e+1},{\widehat{\theta }}_e)>Q({\widehat{\theta }}_e,{\widehat{\theta }}_e)$, the likelihood $p_{{\widehat{\theta }}_{e+1}}(I_{1:N})$ subsequently increases. Thus, one way to guarantee an increase is to choose

The repetition of these two steps, namely the calculation of Eq. (13) and the subsequent maximization Eq. (16), constitute the mechanics of the EM algorithm.

Calculation of Eq. (13) is difficult in most cases because determining an analytical representation of the posterior density $p_{{\widehat{\theta }}_e}({\mathbf{x}}_{\mathbf{p},1:N}|I_{1:N})$ may be intractable even if the joint density p_θ(x_p,1:N, I_1:N) is known. In the context of this work, the nonlinear and non-normal observation model Eq. (10) prevents us from expressing the posterior density as a tractable analytical expression. Consequently, we use the Monte Carlo approximation

Calculation of the approximate posterior densities Eq. (17) is broken into two steps. The first step, known as the “filtering step,” calculates the approximate posterior densities $p_{{\widehat{\theta }}_e}({\mathbf{x}}_{\mathbf{p},k}|I_{1:k})$ for k = 1, …, N; here, we employ the widely used Sampling Importance Resampling (SIR) algorithm [29]. The second step, known as the “smoothing step,” operates backward in time on the posterior densities produced in the filtering step to produce the densities in Eq. (17); here, we employ the Forward-Filtering Backward-Smoothing (FFBS) algorithm [30].

The SIR algorithm calculates the approximate posterior density via importance sampling. Starting at k = 1, M random samples are generated from the initial condition model Eq. (8),

Here, ${\tilde{\mathbf{x}}}_{k|j,e}^i\in {ℝ}^3$ denotes the ith randomly generated estimate of the particle position (i = 1, …, M) at EM iteration e at time step k and given measurements I_1:j; the “p” subscript, denoting particle position, is implied and is omitted to simplify notation. Additionally, the tilde denotes that the sample has not yet been resampled, a process which is discussed below. Instead, if k > 1 then the M random samples are generated according to the motion model Eq. (6),

The FFBS algorithm is a deterministic procedure (unlike SIR, which is stochastic) which calculates the posterior densities $p_{{\widehat{\theta }}_e}({\mathbf{x}}_{\mathbf{p},k}|I_{1:N})$ for k = 1, …, N. To do so, it operates backward in time on the estimates produced by SIR and produces a new series of “smoothed” weights; the algorithm does not produce a new series of locations and takes ${\mathbf{x}}_{k|k,e}^i={\mathbf{x}}_{k|N,e}^i$. Note that at k = N, the smoothed weights are equivalent to those produced by SIR. Proceeding backward in time, the smoothed weights are calculated according to

With the SMC approximation Eq. (17), here produced by the SIR and FFBS algorithms, the function $\widehat{Q}(\theta ,{\widehat{\theta }}_e)$ may be calculated and a new set of parameters may be generated from

Given these newly generated parameters, a “refined” posterior density can then be calculated (again, here via SIR and FFBS) and a “refined” set of parameters may be subsequently estimated. This iterative procedure constitutes the SMC-EM algorithm.

In the context of the models pertaining to confocal extremum seeking, namely Eq. (9) and Eq. (10), the implementation of the SIR and FFBS algorithms are straightforward: generating new samples in Eq. (20) and Eq. (21) involve sampling from a normal distribution, weighting the samples in Eq. (22) is done by evaluating a Poisson probability mass function, and calculating the smoothed weights Eq. (24) and Eq. (25) involves the evaluation of a normal probability density function. Determining the maximizing parameters in Eq. (26) is also straightforward but tedious. For the initial condition model, the maximizing mean and variance are given as

Note that Eq. (30) possesses no analytical solution and must be evaluated numerically; in this work, the tangent-hyperbolas method was employed [31].

4. Experimental demonstration

4.1. Characterization of the point spread function

Although the tracking algorithm requires no specific knowledge of the PSF aside from it possessing a global maximum, localizing the particle given the experimental data using the SMC-EM algorithm described in Sec. 3 relies on an accurate measurement of the PSF.

For the experiments described here, the PSF of the confocal microscope was measured by performing a three-dimensional raster scan of a fixed, fluorescent point emitter. Quantum dots (Qdot 625, Life Technologies) were adhered to a coverslip using Qmount mounting media (Qdot Qmount, Life Technologies) by following the manufacturer’s protocol. The slide of fixed quantum dots was attached to the sample holder and the laser power was set to 7 μW. A raster scan was performed with a resolution of 30 nm in x and y and 50 nm in z with a range of 2 μm in all three axes. Intensity measurements were acquired with a 1 ms accumulation period.

We assumed the PSF took the form of a three-dimensional, rotated Gaussian function. The rotation is needed to capture misalignment between the optical axis and the z–axis of the nanopositioning stage used to move the sample relative to the focal volume. Mathematically, the PSF model was

The following parameters were fit to the acquired data using MATLAB’s nlinfit method: F₀, x_c, y_c, z_c, σ_x, σ_y, σ_z, ψ_x, ψ_y, ψ_z. The relevant estimated parameters were ${\widehat{\sigma }}_x=216\text{nm}$, ${\widehat{\sigma }}_y=270\text{nm}$, ${\widehat{\sigma }}_z=533\text{nm}$, ${\widehat{\psi }}_x=11.3°$, ${\widehat{\psi }}_y=-52.2°$, and ${\widehat{\psi }}_z=131.6°$.

The resulting fit can be visualized in three dimensions (Fig. 2) and in two dimensions (Fig. 3). Specifically, Fig. 2 compares the data and the model side-by-side and defines three mutually-orthogonal planes (denoted by the colors green, orange, and magenta) that are formed by rotating the standard Cartesian unit vectors centered about the estimated particle center using the estimated rotations. These planes are shown in cross-section in Fig. 3.

 

Fig. 2 Three-dimensional PSF measurements (left) and the corresponding Gaussian model (right), calculated by a nonlinear least-squares fit. Three planes (magenta, orange, and green) are shown as cross sections in Fig. 3. These planes were determined by the ZYX rotation in Eq. (4.1) with the three rotations ( ${\widehat{\psi }}_x$, ${\widehat{\psi }}_y$, ${\widehat{\psi }}_z$) given by the least-squares fit. The intensity values in the measured PSF are normalized by the maximum measured intensity value, and the intensity values in the model PSF are normalized by the peak intensity value calculated by least-squares fit.

Download Full Size | PDF

 

Fig. 3 Two-dimensional point spread function measurements (top row) and their corresponding Gaussian models (bottom row) through the three planes depicted in Fig. 2.

Download Full Size | PDF

4.2. Tracking experiment

To demonstrate the applicability of this method, unconjugated quantum dots (QD) were tracked in an agarose hydrogel. To prepare a sample, 1.745 g of water and 0.062 g of powdered agarose (A9539, Sigma-Aldrich) were combined and heated until boiling. Subsequently, 1.924 g of glycerol (G31-1, Fisher Scientific) was added to the mixture. In the first two experiments described below, 200 nL of a solution (2 nM per 50 mM borate) of carboxyl QDs (Qdot 625, Life Technologies) was then added and the mixture was stirred for approximately one minute, poured onto a glass coverslip which was then adhered to a sample slide, and allowed to cool until it reached room temperature. In the third experiment described below, approximately 1 μL of the carboxyl QDs solution was placed on the cover slip followed by approximately 20 μL of the agarose-gylcerol micxture. This was manually stirred for approximately one minute, adhered to a second, larger cover slip, sealed, and then allowed to cool to room temperature. The two preparations produced similar results.

The hydrogel sample was fixed to the piezoelectric nanostage via a sample holder. The laser power was set to 25 μW as measured after the beam expander and prior to the collimating optics. The intensity was measured using a 1 ms accumulation period. This duration was selected to be fast relative to the (expected) motion while still providing significant signal in each integration period. It is, however, a user-selectable parameter limited primarily by the output emission rate and the capabilities of the avalanche photodiode; for example, faster tracking can be achieved by reducing the accumulation rate and, if needed, increasing excitation power to compensate. Prior to tracking, the background noise was measured to be approximately 4 counts per millisecond in a vacant section of the sample. A portion of the sample with a concentration of less than 0.002 QDs/μm³ was selected, as calculated by estimating the number of visible QDs on the CCD image and assuming a depth of field of 1 μm from the size of the PSF. Tracking was initiated when the accumulated intensity reached a value higher than 40 counts / ms and was subsequently disabled when the intensity dropped to zero. The tracking algorithm parameters were set to have a scanning radius of R = 50 nm, a convergence gain of K_p = 5×10⁻⁴, and angular frequencies of (7, 15, 22) Hz. Individual particles were selected through observation on the CCD camera. Note that given the lack of blinking in the observed intensities, tracked particles were likely an agglomeration of QDs rather than individual dots.

The results of one tracking run are shown in Figs. 4 and 5. The experimental output of length 50 s was processed by the SMC-EM algorithm using 50 Monte Carlo estimates. As described in Sec. 3, we assumed the particle moved via a combination of unconfined anisotropic diffusion and flow; thus, the motion parameters to be estimated included three diffusion coefficients and three velocities. The peak intensity of the PSF was assumed unknown and was estimated in addition to the six aforementioned motion parameters. We note that the length of the trajectory was not due to loss of tracking during the experiment but by the amount of memory required by the existing implementation of the SMC-EM algorithm. Our current version trades off computation speed for memory usage; much larger runs could be processed if a different tradeoff were made. A much longer trajectory is considered in the third experiment below, though it was analyzed by breaking it apart into smaller segments. After 1000 iterations of the SMC-EM algorithm, the estimated diffusion coefficients were (0.139, 0.082, 0.087) μm²/s and the estimated velocities were (−0.009, −0.001, 0.003) μm/s, in x, y, and z, respectively; the peak intensity was 1376 counts / ms. This peak intensity corresponds to a measurement directly on top of the particle; the actual measured intensity (top left plot in Fig. 4) is significantly lower due to the separation between the confocal volume and the particle. Note that the estimated diffusion was anisotropic and significantly larger in the x–direction; this is likely due to confinement of the tracked QD in a pore of the hydrogel with ashape allowing more motion along this axis. The estimated particle location, depicted in blue in Fig. 4, was determined by calculating the weighted sample mean of the estimated posterior density as a function of time; this is plotted relative to the measured focal volume position, depicted in black. Since the SMC-EM algorithm returns an estimate of the entire posterior distribution, the variance on the estimate can also be calculated; in this tracking run the average localization error was (16.1, 14.3, 13.9) nm. Additionally, Fig. 4 shows the measured intensity, in black, the estimated intensity, in blue, and the residual difference between the estimated and actual measurement. The estimated intensity was determined by evaluating the Gaussian model as a function of the difference of the measured focal volume position and the estimated particle location and the residual calculated as the difference between the measured and estimated intensity values. The mean of the residual was essentially zero at −1.24 counts/ms and displays the shot-noise character of the measured signal. A three-dimensional plot of the inferred particle position, subsampled by a factor of five for clarity in the image, is given in Fig. 5.

 

Fig. 4 The inferred three-dimensional position of a quantum dot (blue) diffusing in a hydrogel relative to the position of the focal volume (black) which followed the particle in real-time using the extremum seeking method described in this work. The particle position was inferred by employing the SMC-EM algorithm. The left three plots show the (top) x, (middle) y, and (bottom) z time-series while the right three plots show the (top) measured intensity, (middle) estimated intensity, and (bottom) residual intensity, given by the difference between the measured and estimated intensity values. The mean of the residual, shown as a solid blue line, was −1.24 counts / ms.

Download Full Size | PDF

 

Fig. 5 The inferred three-dimensional trajectory of a quantum dot diffusing in a hydrogel as parametrized by time and subsampled from the estimated trajectory by a factor of five. Three planar cross-sections are also shown. The quantum dot was tracked in a confocal microscope using the method presented in this work; the resulting particle position was inferred using the SMC-EM algorithm.

Download Full Size | PDF

The results of a second, shorter tracking run are shown in Figs. 6 and 7. The 3-D trajectory in Fig. 7 is subsampled by a factor of two for clarity of the figure. Once again, the SMC-EM algorithm with 50 Monte Carlo estimates was used to estimate the particle trajectory, motion model parameters, and peak intensity. After 1000 iterations of the algorithm, the estimated diffusion coefficients were (0.0082, 0.0081, 0.019) μm²/s and the estimated velocities were (0.0203, −0.0073, 0.11) μm/s in x, y, and z, respectively while the peak intensity was 108.9 counts / ms. This particle was significantly less bright than in the previous tracking run, yielding more variation between the measured and estimated intensities as seen in Fig. 6. The average position error throughout the run was (7.50, 7.00, 10.4) nm.

 

Fig. 6 The inferred three-dimensional position of a quantum dot (blue) diffusing in a hydrogel relative to the position of the focal volume (black) which followed the particle in real-time using the extremum seeking method described in this work. The particle position was inferred by employing the SMC-EM algorithm. The left three plots show the (top) x, (middle) y, and (bottom) z time-series while the right three plots show the (top) measured intensity, (middle) estimated intensity, and (bottom) residual intensity, given by the difference between the measured and estimated intensity values. The mean of the residual, shown as a solid blue line, was 2.97 counts / ms.

Download Full Size | PDF

 

Fig. 7 The inferred three-dimensional trajectory of a quantum dot diffusing in a hydrogel as parametrized by time and subsampled from the estimated trajectory by a factor of two. Three planar cross-sections are also shown. The quantum dot was tracked in a confocal microscope using the method presented in this work; the resulting particle position was inferred using the SMC-EM algorithm.

Download Full Size | PDF

As a final example, the results of a much longer run of 2000 s in duration are shown in Figs. 8 and 9. Since this run was too long for analysis using the existing SMC-EM implementation, it was divided into 100 segments, each of 20 s, which were then each individually analyzed with 80 Monte Carlo estimates. Under the assumption that the local environment in the hydrogel is uniform, this allowed for the calculation of the variance on the diffusion coefficient estimation. Approximately 700 iterations of the SMC-EM algorithm were performed with segments ranging from 680 to 720 iterations. The resulting combined time traces, measured intensity, estimated intensity, and residual, are shown in Fig. 8. Note that since each 20 s segment was analyzed independently, there are often large errors at the end points. For clarity of presentation in the figures, trajectory points at the trajectory ends are omitted with the understanding that quantitative results hold inside each 20 s segment but not necessarily across segments. The residual is also larger relative to the signal than in the runs shown in Figs. 4 and 6. This is likely due to a realignment of the optical path in the confocal microscope prior to this experiment, leading to a small change in the PSF of the instrument relative to the model used for estimation. Four 40 second portions of the entire trajectory, each subsampled by a factor of four, are shown in 3-D in Fig. 9. (Note that the entire trajectory consisted of two million data points and showing it in 3-D would require a subsampling high enough that the plot would no longer be representative of the original trajectory.)

 

Fig. 8 The inferred three-dimensional position of a quantum dot (blue) diffusing in a hydrogel relative to the position of the focal volume (black) which followed the particle in real-time using the extremum seeking method described in this work. The long tracking run was broken into 20 s segments, the particle position was inferred on each segment using the SMC-EM algorithm, and then the segments stiched together. The left three plots show the (top) x, (middle) y, and (bottom) z time-series while the right three plots show the (top) measured intensity, (middle) estimated intensity, and (bottom) residual intensity, given by the difference between the measured and estimated intensity values. The mean of the residual, shown as a solid blue line, was 2.97 counts / ms.

Download Full Size | PDF

 

Fig. 9 Four sections of the inferred three-dimensional trajectory of a quantum dot diffusing in a hydrogel as parametrized by time and subsampled from the estimated trajectory by a factor of four. Three planar cross-sections are also shown. The quantum dot was tracked in a confocal microscope using the method presented in this work; the resulting particle position was inferred using the SMC-EM algorithm.

Download Full Size | PDF

The estimated parameters across all segments are shown in Fig. 10 leading to the overall estimates of (0.0082 ± 0.0014, 0.0082 ± 0.0015, 0.0134 ± 0.0019) μm²/s for the diffusion coefficients and (0.006 ± 0.0342, −0.0027 ± 0.0308, 0.0005 ± 0.0401) μ/s for the velocities in x, y, and z respectively. The estimated intensity was 82.35 ± 17.56 cnts/ms. The position error over the entire run was (10.53, 10.12, 12.72) nm.

 

Fig. 10 SMC-EM estimates of (left) diffusion coefficients and (right) velocities for each of the 20 s segments of the run shown in Figs. 8. The overall motion parameters were (0.0082 ± 0.0014; 0.0082 ± 0.0015; 0.0134 ± 0.0019) μm²/s for the diffusion and (0.006 ± 0.0342, −0.0027 ± 0.0308, 0.0005 ± 0.0401) μ/s for the velocities in x, y, and z, respectively.

Download Full Size | PDF

The extremum seeking algorithm is a straightforward, simple-to-implement algorithm that requires little-to-no modification to existing instrumentation. It operates directly on the intensity signal, allowing for very high sample rates that are limited only by the capabilities of the avalanche photodiode. If desired, it would be straightforward to add a second channel in a second color for simultaneous wide-field detection to provide information about the particle’s environment.

Funding

National Science Foundation (NSF) (0845742,1352729,1509084).