Tracking single fluorescent particles in three dimensions via extremum seeking

: The ability to track single ﬂuorescent particles in three-dimensions with sub-di ﬀ raction limit precision as well as sub-millisecond temporal resolution has enabled the understanding of many biophysical phenomena at the nanometer scale. While there are several techniques for achieving this, most require complicated experimental setups that are expensive to implement. These methods can o ﬀ er superb performance but their complexity may be over-whelming to the end-user whose aim is only to understand the feature being imaged. In this work, we describe a method for tracking a single ﬂuorescent particle using a standard confocal or multi-photon microscope conﬁguration. It relies only on the assumption that the relative position of the measurement point and the particle can be actuated and that the point spread function has a global maximum that coincides with the particle’s position. The method uses intensity feedback to calculate real-time position commands that “seek” the extremum of the point spread function as the particle moves through its environment. We demonstrate the method by tracking a di ﬀ using quantum dot in a hydrogel on a standard epiﬂuorescent confocal microscope.


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][6][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][9][10][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][16][17].Tracking of particles diffusing at rates on the order of 0.1 µm 2 /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 2 /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 2 /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 threedimensional 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 2 /s at moderate fluorescence intensities with relatively slow piezo-actuated state-scanning, and at rates on the order of 90 µm 2 /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].

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 where θ(t) and φ(t) are auxiliary variables.The intensity measurement taken at time t and at location x s (t) is denoted by I(t).The user is free to select the four strictly-positive parameters, namely R, ω 1 , ω 2 , and K p , so that the desired tracking behavior is obtained.R defines the desired steady-state radial distance between the particle and the focal volume for a fixed, non-moving particle.The rate of convergence to the sphere is controlled by K p ; smaller values result in a slower convergence while larger values yield faster convergence but may cause instability and loss of tracking; additionally, choosing K p too large will amplify measurement noise and can potentially cause the sensor to exhibit random walk-like behavior.Upon reaching the steady-state sphere, the focal volume oscillates with three angular rates: We note that one of the axes, here denoted as z, oscillates with only one frequency, namely ω 1 .The other two axes, here denoted x and y, oscillate with both ω 1 − ω 2 and ω 1 + ω 2 frequencies.To better illustrate the role of these parameters, it is useful to look at the steady-state distance between the focal volume and the stationary particle center, given by where ᾱ and β are constants that depend on the controller parameters, the initial conditions, and the shape of the field.
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 with auxiliary variables defined by Assuming a relatively fast update frequency relative to the largest of the three orbital frequencies (i.e., ω 1 + ω 2 ), 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 ω 1 and ω 2 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.

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 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.
does not provide an estimate of the particle position.Notationally, the output is a collection x s,1:N x s,1 , x s,2 , . . ., x 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 {I 1 , I 2 , . . ., 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 x p,1:N x p,1 , x p,2 , . . ., x 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 θML is defined by θML = arg min and where p θ I 1:N , x s,1:N is the likelihood of the measurements parametrized by the unknown parameters θ ∈ R n θ .This procedure was presented in [23] where it was described in the context of tracking fluorescent particles using widefield fluorescence measurements; code is available at an online repository [26].The approach is general and can be applied to a wide variety of measurement techniques; in the remainder of this section we describe its use in the current setting where confocal measurements are used.We highlight again that estimation is independent of the tracking algorithm in Sec. 2 and is a post-processing step to analyze the acquired data.We assume the particle location x p,k moves stochastically in discrete time according to the probability density with the superscript denoting "motion model."The motion model may take any form as long as it obeys the Markov property and is parametrized by an unknown parameter θ.Modeling diffusion, for example, could take Eq. ( 6) as a multivariate normal distribution with mean x p,k and covariance matrix with 2D∆t on its diagonal with the diffusion coefficient D a parameter to be estimated.
Since the particle location is unknown and cannot be measured directly, intensity measurements I k ∈ R are taken according to a different probability distribution with the superscript denoting "observation model."The observation model may take any form as long as the current observation depends only on the current state and is parametrized by an unknown parameter θ.For example, Eq. ( 7) could take the form of a Poisson distribution to model optical shot noise.Lastly, we assume the initial location of the particle is distributed according to with the superscript denoting "initial condition model."It is often convenient to assume the initial location of the particle is, for example, normally distributed with unknown mean and variance.
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 where N [x](µ, σ 2 ) denotes a normal distribution with mean µ and standard deviation σ 2 evaluated at x.The two unknown parameters for this axis are the velocity V x and the diffusion coefficient D x .Note that the y and z axes are of identical form, and, because of the assumed independence among axes, the complete probability distribution p(x p,k +1 |x p,k ) is just the product of the individual distributions.Thus, there are six parameters to be estimated that are related to the motion model: 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 where P[x](λ) denotes a Poisson distribution with mean λ evaluated at x, N bgd denotes the average number of background counts per time interval, F PSF denotes the PSF model and the parameter G determines the maximum intensity of the PSF and is unknown.As described in Sec.4.1, in this work we assume that F PSF takes the form of a three-dimensional rotated Gaussian with parameters experimentally determined.Thus, there is only one parameter to be estimated pertaining to the observation model: G.
The initial condition x p,1 is assumed to be normally distributed with unknown mean µ = µ x , µ y , µ z T and variance Thus, there are six parameters to be estimated pertaining to the initial condition model: µ x , µ y , µ z , σ 2 x , σ 2 y , σ 2 z .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 θ = arg max θ log p θ (I 1 , . . ., I N ) (12) using the collection of scalar measurements {I 1 , . . ., I N }.Unfortunately, the probability density in Eq. ( 12) is unknown, so an approximation must be made.Consequently, we use the expectation maximization (EM) algorithm to form a numerical approximation of the ML estimate of θ.
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(θ, θe ) is a weighted averaging of the logarithm of the joint density with respect to the particle location given some estimate of θ, namely θe .As was originally presented in [27]), the simplification of Eq. ( 13) with respect to Eqs. ( 6)-( 8) yields where It was shown in [28] that for each choice of θe+1 such that Q( θe+1 , θe ) > Q( θe , θe ), the likelihood p θ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 θe (x 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 which simplifies Eq. ( 14) and Eq. ( 15) to where Q2 and where δ represents the Dirac delta function.Here, the posterior density p θe (x p ,k |I 1:N ) and the pairwise joint posterior density p θe (x p ,k , x p ,k +1 |I 1:N ) are approximated as a finite sum of weighted delta functions Eq. ( 17), and the weights and locations of the delta functions are determined by SMC methods.Since SMC methods are widely discussed in the literature, we will only briefly detail the two specific SMC algorithms employed in this work.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 θe (x 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, xi k | j,e ∈ R 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), for i = 1, . . ., M; note that the samples are randomly generated with respect to the resampled randomly generated samples from the previous timestep.For each individual timestep k = 1, . . ., N , each random sample i is weighted according to the observation model Eq.(7).In other words, for i = 1, . . ., M. The weights are then scaled so that their sum over i is unity and the estimates resampled according to their relative weights.In this work, M new estimates j were generated, with corresponding weights w j k |k ,e set to 1/M.To summarize, the SIR algorithm performs three operations at each time step k: first, M samples are randomly generated via Eq.( 20) (if k = 1) or Eq. ( 21) (if k > 1); secondly, the samples are weighted by the observation model according to Eq. ( 22); finally, M new samples are generated according to Eq. ( 23).The output of the SIR algorithm is the approximate posterior density p θe (x p ,k |I 1:k ) for k = 1, . . ., N .
The FFBS algorithm is a deterministic procedure (unlike SIR, which is stochastic) which calculates the posterior densities p θe (x 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 for k = (N − 1), . . ., 1.We also require an approximation of the pairwise joint density Eq. ( 17); these were calculated in [27] to be With the SMC approximation Eq. ( 17), here produced by the SIR and FFBS algorithms, the function Q(θ, θ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 for the x coordinate; the y and z coordinates are of identical form.For the motion model, the maximizing speed is and the maximizing diffusion coefficient is for the x axis; again, the y and z coordinates are of identical form.For the observation model, the maximizing gain Ĝe+1 is given by the condition with the expected intensity at the focal volume location x s ,k defined by Note that Eq. ( 30) possesses no analytical solution and must be evaluated numerically; in this work, the tangent-hyperbolas method was employed [31].

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 where x c = x c , y c , z c T is the particle center, ψ = ψ x , ψ y , ψ z T are the rotation angles of the PSF, and Σ = diag σ 2 x , σ 2 y , σ 2 z represents a diagonal matrix of PSF widths.The peak intensity of Eq. ( 32) is assumed to be one; the actual intensity is included as a parameter in the observation density Eq. (10).The rotation matrix R (ψ) describes a ZYX rotation about the particle center, The following parameters were fit to the acquired data using MATLAB's nlinfit method: 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 mutuallyorthogonal 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.

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 3 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 −4 , 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 2 /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 a shape 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.
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 esti- mated diffusion coefficients were (0.0082, 0.0081, 0.019) µm 2 /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.
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 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 2 /s for the diffusion coefficients and (0.006 ± 0.0342, −0.0027 ± 0.0308, 0.0005 ± 0.0401) µm/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.
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

Fig. 1 .
Photon counts x i k |k ,e = x i k | N,e .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 w i k | N,e = w i k |k ,e ψ y , ψ z .The relevant estimated parameters were σx = 216 nm, σy = 270 nm, σz = 533 nm, ψx = 11.3 • , ψy = −52.2• , and ψz = 131.6 • .

Fig. 2 .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 ( ψx , ψy , ψ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.

Fig. 4 .
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.

Fig. 5 .
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.

Fig. 6 .
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.

Fig. 9 .
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.