One way to achieve spatial resolution using fluorescence imaging—and track single molecules—is to use wide-field illumination and collect measurements over multiple sensors (camera pixels). Here we propose another way that uses confocal measurements and a single sensor. Traditionally, confocal microscopy has been used to achieve high temporal resolution at the expense of spatial resolution. This is because it utilizes very few, and commonly just one, sensors to collect data. Yet confocal data encode spatial information. Here we show that non-uniformities in the shape of the confocal excitation volume can be exploited to achieve spatial resolution. To achieve this, we formulate a specialized hidden Markov model and adapt a forward filtering-backward sampling Markov chain Monte Carlo scheme to efficiently handle molecular motion within a symmetric confocal volume characteristically used in fluorescence correlation spectroscopy. Our method can be used for single confocal volume applications or incorporated into larger computational schemes for specialized, multi-confocal volume, optical setups.

Fluorescent microscopy1,2 (FM) creates contrast between fluorescent molecules and their background. Unlike other microscopy techniques, such as electron microscopy3 or scanning probe microscopy,4 FM is less invasive and the laser intensity used is often limited by the sensitivity of the analysis tools used downstream to extract information from fluorescence data.1,5,6

In particular, wide-field optical setups, commonly encountered in FM, employ an array of regularly spaced light sensors (pixels). On account of these multiple sensors, the location of fluorescing molecules can be pinpointed with high spatial accuracy, and when used in conjunction with super-resolution approaches such as stochastic optical reconstruction microscopy (STORM)7 or photoactivated localization microscopy (PALM),8,9 the labeled molecules can be localized with accuracy approaching ≈20-30 nm, well below the conventional diffraction limit.5 In this context, linking of localized molecules across a sequence of successive frames enables single particle tracking (SPT).10–13 Nevertheless, although wide-field setups are capable of reconstructing dynamics with high spatial accuracy, they typically suffer from poor temporal accuracy.5 This is due to both hardware and software-related reasons. On the hardware side, the acquisition of an entire frame over a large number of sensors needs to be read by the same camera electronics.14,15 On the software side, a sufficient number of photons needs to be collected within individual frame acquisitions for algorithms to successfully localize molecules. As such, frame rates used in SPT are typically low, ≈10 Hz.15–19 Such frame rates are appropriate when imaging stationary or otherwise slowly diffusing molecules in vivo. Additionally, since large portions of the specimen are continuously illuminated, wide-field FM may also cause significant photo-damage to fluorescing molecules or be photo-toxic to the specimen under observation.20,21

Confocal optical setups used in FM provide an alternative. By using a small number of sensors, such as avalanche photodiodes (APDs), or most often just one sensor,2 confocal microscopy reaches data acquisition rates ≈10 GHz. At such high rates, fast biochemical processes are readily monitored with temporal accuracy superior to wide-field FM, however, with compromised spatial resolution. Additionally, because long exposures are unnecessary, since molecules are not spatially localized and the illumination is targeted to only a small portion of the specimen, photo-damage caused by confocal illumination is minimal.

Existing SPT analysis,10–12 however, cannot be used to extract spatial information from confocal data. Fundamentally, this is because of the reduced spatial information captured by only one or very few sensors. Yet, a glimpse at confocal data reveals that spatial information is in fact encoded even within the signal provided by just one pixel or avalanche photodiode (APD) (see Fig. 1).

FIG. 1.

Illustration of a confocal setup. (a) A sample trace illustrating photon intensities collected over an experiment’s time course. [(b) and (c)] Lateral (b) and axial (c) cross sections of a 3D confocal volume (green) with the loci of equal intensity highlighted (purple). A point emitter, such as a fluorescent molecule, located at any point on the purple surface outputs the same average number of photons consistent with the intensities in (a).

FIG. 1.

Illustration of a confocal setup. (a) A sample trace illustrating photon intensities collected over an experiment’s time course. [(b) and (c)] Lateral (b) and axial (c) cross sections of a 3D confocal volume (green) with the loci of equal intensity highlighted (purple). A point emitter, such as a fluorescent molecule, located at any point on the purple surface outputs the same average number of photons consistent with the intensities in (a).

Close modal

In particular, due to the non-uniform illumination used in confocal setups, the location of a molecule is encrypted in the instantaneous photon emission rates and ultimately reflected in the acquired traces. Nevertheless, to extract spatial information from such traces, on account of the symmetry inherent to typical confocal volumes, we must account for the fact that multiple locations of a molecule within a confocal volume (for example, those at similar distances from the focal point for the confocal volume shown) may give rise to the same signal (see Fig. 1). In addition, we must account for molecular motion across the illuminated volume and noise.

Here, we specifically tackle these challenges and develop the necessary computational framework to obtain spatial resolution from confocal measurements. Briefly, we develop and validate a Bayesian method that can be used on raw measurements obtained from just one sensor with the characteristics of an APD, as found in most common confocal setups, and provide diffusion coefficient estimates, similar to bulk fluorescence correlation spectroscopy (FCS) or fluorescence recovery after photobleaching (FRAP). Complementary to these methods that extract only dynamical parameters, such as diffusion coefficients, we also extract trajectories of the molecules. Additionally, our method reconstructs not only the most probable trajectory similar to wide-field SPT,5,13,22 but full distributions over molecular trajectories that are consistent with the raw measurements. To keep the computational cost low, we develop a formulation based on hidden Markov models and adapt tools from Computational Statistics such as non-linear filtering and Markov chain Monte Carlo (MCMC) sampling that we describe in Sec. II.

In this section, we describe the theoretical foundation of our method. That is, we first introduce the necessary formalism and subsequently present the computational tools for its evaluation. Throughout this study, we assume that successive measurements, in a confocal experiment, are reported at times tk, for k = 1, …, K, where t1 marks the time of the first measurement and tK marks the time of the last one. Although typically measurements are reported at equal time intervals tktk−1, the formalism we describe below can handle non-equal time intervals too with no further modifications.

Assuming that the time levels tk are known, we start from individual photon detection counts at those times, w ¯ = ( w 1 , w 2 , , w K ) . These counts are taken to be the raw experimental output. Here, w k is the number of photons detected between tk−1 and tk. The source of these photons is either a single molecule in the vicinity of the illuminated volume or some unspecified fluorescent background. Fluctuations in w ¯ arise from photon shot noise (stochastic photon emission) and the motion of the fluorescent molecule within the inhomogeneously illuminated confocal volume. We denote the location of the single molecule by r ¯ = ( r 1 , r 2 , , r K ) , where rk is the location of the molecule relative to the center of the confocal volume at time tk.

For clarity, in the following presentation, we first consider a simpler 1D confocal model, where rk = xk. We generalize to 3D, where rk = (xk, yk, zk), later.

A simple graphical summary of our 1D model is provided in Fig. 2. Below we explain the photon emission and molecular motion separately.

FIG. 2.

Graphical representation of the formulation used in the analysis of time fluorescence traces in 1D. One molecule’s position xk evolves over the course of the experiment with time points tk indexed by k = 1, 2, …, K. Here, xk denotes the molecule’s location in 1D at time tk. At each time point, an observation wk is recorded which combines photon emissions between tk−1 and tk from the molecule and also the background. The diffusion coefficient D determines the evolution of the molecular locations xk which, in turn, influences the photon emission rates and ultimately the recorded photon intensity wk.

FIG. 2.

Graphical representation of the formulation used in the analysis of time fluorescence traces in 1D. One molecule’s position xk evolves over the course of the experiment with time points tk indexed by k = 1, 2, …, K. Here, xk denotes the molecule’s location in 1D at time tk. At each time point, an observation wk is recorded which combines photon emissions between tk−1 and tk from the molecule and also the background. The diffusion coefficient D determines the evolution of the molecular locations xk which, in turn, influences the photon emission rates and ultimately the recorded photon intensity wk.

Close modal

1. Photon emission

To analyze photon count time series, w ¯ , we first model the photon emission probability. This probability captures observation noise and is necessary in linking photon count measurements to molecular locations. Here, we assume that photon counts w k between tk−1 and tk follow Poisson (shot noise) statistics5,23

w k | x k P o i s s o n t k t k 1 μ k ,
(1)

where μk is the instantaneous photon emission rate. Technically, this is the emission rate for the detected photons (i.e., μk excludes photons that are lost in the optical path from the molecule to the detector such as those straying away from the objective or blocked by dichroics). In turn, the instantaneous photon emission rate

μ k = μ b a c k + μ m o l P S F ( x k )
(2)

combines a time-independent background rate μback and a time-dependent molecule photon emission rate μmolPSF(xk). Here, the rate μmol represents the maximum rate of photon emissions of a single molecule which is achieved when the molecule is at the center of the confocal volume, i.e., xk = 0. It is worth mentioning that the parameter μmol in this study, known as the brightness at the center of the confocal volume, is different from the molecular photon emission rate commonly described in the literature. The latter quantity is the average photon count rate, i.e., μmolPSF(x) averaged over x in the effective volume.24,25 Besides μmol, our molecule’s photon emission also depends upon the location of the molecule within the confocal volume through PSF(xk) which combines the effects of excitation and emission point spread functions (PSFs) into a single term.2,26 Namely, similar to existing FCS approaches, we use a single Gaussian27 as an approximation to the product of the two separate point spread functions.28–30 In 1D, our approximation is

P S F ( x ) = exp 2 x 2 ω x 2 ,
(3)

where ωx is the distance (waist) from the focal point24,31,32 when the PSF(x) drops to e−2.

2. Molecular motion

Since the molecule moves stochastically inside the specimen under observation, to reconstruct the dynamics captured in w ¯ , we consider a probabilistic motion model. We assume that the molecule follows isotropic diffusion with its position’s probability density, p(x, t), at time t satisfying Fick’s second law33–35 

p t = D 2 p x 2 ,
(4)

where D is the molecule’s diffusion coefficient. To solve this equation, we assume that at time tk−1, the molecule’s location is xk−1, which provides the initial condition, p(x, tk−1) = δ(xxk−1). The solution at any later time t then takes the form

p ( x , t ) = 1 4 π ( t t k 1 ) D exp ( x x k 1 ) 2 4 ( t t k 1 ) D ,
(5)

which equals to the probability density of a normal random variable with mean xk−1 and variance 2(ttk−1)D. Thus, our motion model describing the position of the molecule is

x k | x k 1 , D N o r m a l x k 1 , 2 ( t k t k 1 ) D .
(6)

The quantities which we are interested in, for example, the diffusion coefficient D and the molecule’s trajectory x ¯ = ( x 1 , x 2 , , x K ) are introduced as model variables in the preceding formulation. To estimate values for D and reconstruct trajectories for x ¯ , we follow the Bayesian paradigm.5,36,37 Within the Bayesian paradigm, the variables D and x1, in particular, are parameters of the model and, as such, require prior probability distributions, which we describe next.

To ensure positive values for the diffusion coefficient and also to achieve conjugacy with the motion model, which facilitates the computations below, we use an inverse-Gamma prior probability distribution on D

D I n v G a m m a α D , β D ,
(7)

with hyper-parameters αD and βD.

Due to the symmetries in our confocal point spread function [Eq. (3)], i.e., a molecule at a location +x emits on average the same number of photons as a molecule at location −x, the intensities in w ¯ cannot distinguish these two locations. Hence, we place a prior probability distribution on the initial location x1 that respects this symmetry. Accordingly, in our framework, at the onset of the measuring period, the molecule is equally likely to be located at any of the positions ±x1. To facilitate the computations (see below), we use a symmetric normal distribution

x 1 S y m N o r m a l μ x , σ x 2 .
(8)

Specifically, the probability density of a S y m N o r m a l μ , σ 2 distribution is

p ( x ) = 1 2 N o r m a l x ; μ x , σ x 2 + 1 2 N o r m a l x ; + μ x , σ x 2 .
(9)

In summary, our entire Bayesian model in 1D consists of the following relations:

D I n v G a m m a α D , β D ,
(10)
x 1 S y m N o r m a l μ x , σ x 2 ,
(11)
x k | x k 1 , D N o r m a l x k 1 , 2 ( t k t k 1 ) D ,    k = 2 , , K ,
(12)
w k | x k P o i s s o n ( t k t k 1 ) μ k ,    k = 1 , , K ,
(13)

which have the structure of a Bayesian hidden Markov model38–41 and are illustrated graphically in Fig. 2. This model leads to a posterior probability distribution p ( D , x ¯ | w ¯ ) that captures all variables of interest. The nonlinearity in the point spread function [Eq. (3)] with respect to the variables in x ¯ , however, exclude analytic forms of p ( D , x ¯ | w ¯ ) . For this reason, below we develop a computational scheme exploiting Markov chainMonte Carlo36,42 (MCMC) that can be used to generate pseudo-random samples ( D ( i ) , x ¯ ( i ) ) from p ( D , x ¯ | w ¯ ) , which can be used to obtain any statistic of interest such as mean values and error bars.

Our MCMC algorithm exploits a blocked Gibbs sampling scheme.23,36,42 Accordingly, posterior samples ( D ( i ) , x ¯ ( i ) ) are generated by updating each one of the variables involved sequentially by sampling conditioned on the other variable and the measurements w ¯ . Algorithmically, the steps involved in the generation of ( D ( i ) , x ¯ ( i ) ) are as follows:

  1. Generate the diffusion coefficient D(i) by sampling from the distribution p ( D ( i ) | x ¯ ( i 1 ) , w ¯ ) .

  2. Generate the trajectory x ¯ ( i ) by sampling from the distribution p ( x ¯ ( i ) | D ( i ) , w ¯ ) .

We present the fine details of these steps below.

1. Sampling the diffusion coefficient

As mentioned above, we update the diffusion coefficient D by sampling from p ( D | x ¯ , w ¯ ) , which, due to the specific dependencies of the variables in our formulation, simplifies to p ( D | x ¯ ) . Using Bayes’ rule p ( D | x ¯ ) p ( x ¯ | D ) p ( D ) and because of conjugacy of our prior on Eq. (7), the latter reduces to

D | x ¯ I n v G a m m a α D , β D ,
(14)

where the derivation and explicit formulas for α D and β D are given in the supplementary material.

2. Sampling the trajectory

As we have already mentioned, to update the trajectory x ¯ , we need to sample from p ( x ¯ | D , w ¯ ) , which we achieve through backward sampling.39,43,44 Namely, since the density p ( x ¯ | D , w ¯ ) factorizes as

p ( x ¯ | D , w ¯ ) = p ( x 1 | x 2 , D , w ¯ ) p ( x 2 | x 3 , D , w ¯ ) p ( x K 1 | x K , D , w ¯ ) × p ( x K | D , w ¯ ) ,
(15)

we sample x ¯ sequentially starting from xK and marching backward towards x1. However, before we begin backward sampling, we need to compute the individual densities p ( x K | D , w ¯ ) and p ( x k | x k + 1 , D , w ¯ ) . We compute these recursively using forward filtering.38–40,43,44 Namely, we start our density computations at x1 and march forward towards xK. Below we describe in detail each one of the involved steps.

a. Forward filtering stage.

By Bayes’ rule, the individual densities in Eq. (15) factorize as

p ( x k | x k + 1 , D , w ¯ ) = C k p ( x k + 1 | x k , D ) p ( x k | D , w 1 : k ) ,
(16)

where we use w 1 : k as a shorthand for w 1 , …, w k and Ck is a factor that does not depend on xk. Therefore, for the computation of the densities required to sample xk in backward sampling, it is sufficient to compute only filter densities p(xk|D, w1:k) since we already know p(xk+1|xk, D) from the motion model and can recover Ck by normalization, i.e., x k p ( x k | x k + 1 , D , w ¯ ) d x k = 1 , although in practice we need not compute Ck at all (see the sampling stage below).

To be able to perform the necessary forward filter computations efficiently,45 we use an approximate observation model46 and we replace Eq. (1) with

W k | x k N o r m a l T k ( x k ) , 1 ,
(17)

where W k = 2 w k + 3 8 and T k ( x ) denotes the Anscombe transform,47 which in our context is given by

T k ( x ) = 2 f k ( x ) + 3 8 1 4 f k ( x ) ,
(18)

with f k ( x ) = t k t k 1 μ b a c k + μ m o l P S F x . This transform approximates Poisson random variables by univariate normal random variables47 which are necessary for the computations described next. Our approximation is highly accurate for large fk(x), while acceptable accuracy is maintained so long as fk(x) > 4.

Now, as both our kinetics and (transformed) observation probabilities are normally distributed, i.e., Eqs. (6) and (17), respectively, we can use the existing theory of Kalman filters.48–50 Nevertheless, because the mean, i.e., T k ( x k ) , of the observation distribution depends nonlinearly on the location xk, to apply Kalman filters we need to linearize the observation probability density.51–54 To do so, in Sec. III, we try two common approaches: (i) extended Kalman filters (EKFs),51–55 which apply a local approximation to linearize T k ( x k ) [Eq. (17)] around selected values; (ii) unscented Kalman filters (UKFs),56–58 which apply a global approximation to linearize in a statistical manner the entire density of N o r m a l T k ( x k ) , 1 in Eq. (17).

Despite the linearization involved, neither naive EKFs nor naive UKFs serve as appropriate filter approximations to p(xk|D, w 1 : k ) in our framework. This is because both assume that the filters themselves are normal densities. This assumption is particularly problematic for confocal observations with ±x symmetry which implies equal probabilities for the molecule to exist in either the negative or positive side of the illumination profile. In 1D, under our symmetric prior in Eq. (8) and isotropic dynamics in Eq. (6), this implies that our filter p(xk|D, w 1 : k ) has two, equally weighted and symmetrically placed across the origin, modes. We therefore modify EKFs and UKFs to account for bimodal densities appropriate for symmetric confocal setups.

More specifically, unlike the conventional filters51–58 that consider a uni-modal filter using p x k | D , w 1 : k N o r m a l x k ; m k , c k , for properly chosen mk and ck, we compute filter densities of the form

p x k | D , w 1 : k S y m N o r m a l x k ; m k , c k ,
(19)

where S y m N o r m a l m k , c k denotes the symmetric normal distribution with the density given in Eq. (9) (see also Fig. 3). Our approximation, as we have explained earlier, arises due to the linearization in the observation distribution. As EKF and UKF compute mk and ck differently, we describe the procedures, adapted to Eq. (19), separately in the supplementary material.

FIG. 3.

A symmetric filter probability density arises due to the symmetry inherent to the confocal volume and isotropy of free diffusion. Here, the shaded region indicates the illumination profile. Vertical dashed lines indicate the mean values of the filter’s modes in the positive and negative sides of the illuminated region. Horizontal dashed lines indicate the variances of the two modes.

FIG. 3.

A symmetric filter probability density arises due to the symmetry inherent to the confocal volume and isotropy of free diffusion. Here, the shaded region indicates the illumination profile. Vertical dashed lines indicate the mean values of the filter’s modes in the positive and negative sides of the illuminated region. Horizontal dashed lines indicate the variances of the two modes.

Close modal
b. Backward sampling stage.

After computing the filter densities p(xk|D, w 1 : k ) in the forward filtering stage, either through EKF or UKF, we are able to sample the locations of the molecule using backward sampling as in Eq. (15). Specifically, according to Eq. (16), the sampling densities on the individual locations become S y m N o r m a l x K ; m K , c K for the last time level and S y m N o r m a l x k ; m k , c k × N o r m a l x k + 1 ; x k , 2 D ( t k + 1 t k ) for the preceding ones,38 where mk and ck are the filter parameters computed in the forward filtering stage. The resulting distributions are

x K 1 2 N o r m a l m K , c K + 1 2 N o r m a l + m K , c K , x k 1 1 + exp + 2 m k x k + 1 c k + 2 D t k + 1 t k N o r m a l x k + 1 c k 2 D t k + 1 t k m k c k + 2 D t k + 1 t k , 2 D c k t k + 1 t k c k + 2 D t k + 1 t k + 1 1 + exp 2 m k x k + 1 c k + 2 D t k + 1 t k N o r m a l x k + 1 c k + 2 D t k + 1 t k m k c k + 2 D t k + 1 t k , 2 D c k t k + 1 t k c k + 2 D t k + 1 t k , k = 1 , , K 1 ,
(20)

which are directly sampled in our Gibbs sampling scheme.

The full case of a 3D confocal volume, illustrated graphically in Fig. 4, is a straightforward generalization of the 1D case we described earlier. More precisely, the full Bayesian model, in the form of a hidden Markov model, is

D I n v G a m m a α D , β D ,
(21)
x 1 S y m N o r m a l μ x y , σ x y 2 ,
(22)
y 1 S y m N o r m a l μ x y , σ x y 2 ,
(23)
z 1 S y m N o r m a l μ z , σ z 2 ,
(24)
x k | x k 1 , D N o r m a l x k 1 , 2 ( t k t k 1 ) D ,    k = 2 , , K ,
(25)
y k | y k 1 , D N o r m a l y k 1 , 2 ( t k t k 1 ) D ,    k = 2 , , K ,
(26)
z k | z k 1 , D N o r m a l z k 1 , 2 ( t k t k 1 ) D ,    k = 2 , , K ,
(27)
w k | x k P o i s s o n ( t k t k 1 ) μ k ,    k = 1 , , K ,
(28)

with an instantaneous rate μk = μback + μmolPSF(xk, yk, zk). In principle, our treatment is general enough to treat any PSF; though, for simplicity, we restrict our discussion here to a 3D Gaussian of the form

P S F ( x , y , z ) = exp 2 x 2 ω x y 2 + y 2 ω x y 2 + z 2 ω z 2 ,
(29)

where ωxy and ωz are the minor and major semi-axis of the 3D Gaussian.24,31,32

FIG. 4.

A graphical representation of the formulation used in the analysis of a time fluorescence trace in 3D. One molecule evolves over the course of the experiment t1, t2, …, tK. Here, xk, yk, zk denote the location of the molecule in Cartesian coordinates at time tk. More details are provided in the caption of Fig. 2.

FIG. 4.

A graphical representation of the formulation used in the analysis of a time fluorescence trace in 3D. One molecule evolves over the course of the experiment t1, t2, …, tK. Here, xk, yk, zk denote the location of the molecule in Cartesian coordinates at time tk. More details are provided in the caption of Fig. 2.

Close modal

For the 3D model above, the corresponding posterior is p ( D , x ¯ , y ¯ , z ¯ | w ¯ ) . Due to isotropy of the kinetics, we are able to sample each of x ¯ , y ¯ , and z ¯ separately using forward filters and backward sampling similar to the 1D case. Accordingly, we generate posterior samples ( D ( i ) , x ¯ ( i ) , y ¯ ( i ) , z ¯ ( i ) ) with the following Gibbs steps:

  1. Generate the diffusion coefficient D(i) by sampling from the distribution p ( D ( i ) | x ¯ ( i 1 ) , y ¯ ( i 1 ) , z ¯ ( i 1 ) , w ¯ ) .

  2. Generate the trajectory x ¯ ( i ) by sampling from the distribution p ( x ¯ ( i ) | D ( i ) , y ¯ ( i 1 ) , z ¯ ( i 1 ) , w ¯ ) .

  3. Generate the trajectory y ¯ ( i ) by sampling from the distribution p ( y ¯ ( i ) | D ( i ) , x ¯ ( i ) , z ¯ ( i 1 ) , w ¯ ) .

  4. Generate the trajectory z ¯ ( i ) by sampling from the distribution p ( z ¯ ( i ) | D ( i ) , x ¯ ( i ) , y ¯ ( i ) , w ¯ ) .

To benchmark our model, we apply it to synthetic datasets w ¯ generated through common pseudorandom simulations33,59–62 that mimic data generated in confocal microscopy. Then, by applying our sampling scheme, we estimate the diffusion coefficient D and reconstruct the trajectory r ¯ using the model posterior p ( D , r ¯ | w ¯ ) . We sample and compare the posterior in three cases: (i) exact, where we do not use any of the approximations in the filter computation, (ii) EKF, and (iii) UKF. Exact sampling is described in thesupplementary material. It is worth noting that exact sampling is inefficient for real applications due to its exceedingly long running time which is the reason we only used it as a benchmark. For all cases, we compare the results with the ground truth diffusion coefficient D and trajectory r ¯ used to generate the synthetic data.

1. 1D model

From Figs. 5 and 6, we immediately note that the posteriors over positions are sharper, when the molecule is closer to the center of the illuminated region in all cases, i.e., whether exact, EKF, and UKF. This occurs because the molecules emit more photons near the center of the confocal volume which, in turn, leads to cleaner data, thereby sharpening our posterior estimates.

FIG. 5.

Posterior estimates over the diffusion coefficient in 1D. [(a1)–(a3)] The posterior probability distribution over the diffusion coefficient using the exact sampler, EKF and UKF, respectively, applied to the trajectory shown in Fig. 6. The ground truth of the diffusion coefficient is shown by the dashed lines. To generate the trajectory, we use a single diffusive molecule with a diffusion coefficient of 10 μm2/s passing through an illumination/excitation profile (PSF) of size ωx set at 0.3 μm. The molecular μmol and background μback photon emission rates are 105 and 103 photons/s, respectively.

FIG. 5.

Posterior estimates over the diffusion coefficient in 1D. [(a1)–(a3)] The posterior probability distribution over the diffusion coefficient using the exact sampler, EKF and UKF, respectively, applied to the trajectory shown in Fig. 6. The ground truth of the diffusion coefficient is shown by the dashed lines. To generate the trajectory, we use a single diffusive molecule with a diffusion coefficient of 10 μm2/s passing through an illumination/excitation profile (PSF) of size ωx set at 0.3 μm. The molecular μmol and background μback photon emission rates are 105 and 103 photons/s, respectively.

Close modal
FIG. 6.

Posterior estimates of the molecule trajectory in 1D. (a) A sample trace illustrating the number of photons collected over time. [(b1)–(b3)] Trajectory estimates of the exact sampler, EKF, and UKF, respectively, applied to the trace shown in (a). [(c1)–(c3)] Difference between the median of the trajectories in (b1)–(b3) and the ground truth trajectory. The colored background represents the gradient of the illumination profile with the maximum at the center. See the caption of Fig. 5 for details on the parameters used to generate the trace in (a).

FIG. 6.

Posterior estimates of the molecule trajectory in 1D. (a) A sample trace illustrating the number of photons collected over time. [(b1)–(b3)] Trajectory estimates of the exact sampler, EKF, and UKF, respectively, applied to the trace shown in (a). [(c1)–(c3)] Difference between the median of the trajectories in (b1)–(b3) and the ground truth trajectory. The colored background represents the gradient of the illumination profile with the maximum at the center. See the caption of Fig. 5 for details on the parameters used to generate the trace in (a).

Close modal

As shown in Figs. 5 and 6, the posteriors of the diffusion coefficient learned by EKF and UKF are close to the exact posterior learned by the exact random walk sampler. Also, the estimation of the molecule’s trajectory when the molecule is close to the center of the confocal volume matches, almost exactly, the ground truth.

We note that, in Fig. 6, while the molecule finds itself either on the positive or negative side of the 1D illumination profile, the posterior has a symmetric probability density around the origin. Thus, our framework resolves where the molecule is with respect to the confocal center but cannot resolve whether it is on the left- or right-hand side of the origin.

In Fig. 7, we compare the results obtained using our symmetric (bi-modal) EKFs and UKFs described in Sec. II with those obtained by conventional (uni-modal) EKF and UKF filters without any further modifications. A glance at the trajectory estimates suggest that, indeed, uni-modal filters as approximations to the true bi-modal filters are inadequate for confocal applications with symmetric confocal volumes around the origin. The inadequacy is particularly pronounced in the case of UKF.

FIG. 7.

Comparison of symmetrized (bi-modal) and conventional (uni-modal) EKFs and UKFs. (a) A sample trace illustrating the number of photons collected over time. [(b1) and (b2)] Trajectory estimates using uni-modal and bi-modal EKFs, respectively, applied to the trace shown in (a). [(c1) and (c2)] Trajectory estimates using uni-modal and bi-modal UKFs, respectively, applied to the trace shown in (a). See the caption of Fig. 5 for details on the parameters used.

FIG. 7.

Comparison of symmetrized (bi-modal) and conventional (uni-modal) EKFs and UKFs. (a) A sample trace illustrating the number of photons collected over time. [(b1) and (b2)] Trajectory estimates using uni-modal and bi-modal EKFs, respectively, applied to the trace shown in (a). [(c1) and (c2)] Trajectory estimates using uni-modal and bi-modal UKFs, respectively, applied to the trace shown in (a). See the caption of Fig. 5 for details on the parameters used.

Close modal

2. 3D model

As the confocal volume exhibits symmetry in all 3 dimensions, for a 3D Gaussian confocal volume (centered by convention at the origin), posteriors are now invariant under a transformation that leaves the normalized distance

d k = x k ω x y 2 + y k ω x y 2 + z k ω z 2
(30)

unchanged. In particular, this distance is symmetric under xk ↦ −xk, yk ↦ −yk, and zk ↦ −zk.

Next, we benchmark our method in 3D. Again, we compare the results of EKF and UKF in estimating the trajectory and diffusion coefficient of the molecule. Figure 8 demonstrates a sharper posterior as molecules move toward the center of the confocal volume consistent with the higher photon emission rates observed in this region.

FIG. 8.

Posterior results in 3D. (a) A sample trace illustrating the number of photons collected over time. [(b1) and (b2)] Trajectory estimates using EKFs and UKFs, respectively, applied to the trace shown in (a). [(c1) and (c2)] The posterior probability distribution over the diffusion coefficient using EKFs and UKFs, respectively, applied to the trace in (a). The ground truth of the diffusion coefficients is shown by the dashed lines. To generate the trace in (a), we use a single diffusive molecule with a diffusion coefficient of 10 μm2/s passing through a confocal volume of sizes ωxy and ωz set at 0.3 μm and 0.7 μm, respectively. The molecular μmol and background μback photon emission rates are 105 and 103 photons/s, respectively.

FIG. 8.

Posterior results in 3D. (a) A sample trace illustrating the number of photons collected over time. [(b1) and (b2)] Trajectory estimates using EKFs and UKFs, respectively, applied to the trace shown in (a). [(c1) and (c2)] The posterior probability distribution over the diffusion coefficient using EKFs and UKFs, respectively, applied to the trace in (a). The ground truth of the diffusion coefficients is shown by the dashed lines. To generate the trace in (a), we use a single diffusive molecule with a diffusion coefficient of 10 μm2/s passing through a confocal volume of sizes ωxy and ωz set at 0.3 μm and 0.7 μm, respectively. The molecular μmol and background μback photon emission rates are 105 and 103 photons/s, respectively.

Close modal

To benchmark our model across different time scales, we apply it on molecules diffusing with a range of diffusion coefficients (from slow to fast). Here, an intermediate diffusion coefficient is 1 μm2/s, where the molecule emits appreciably (say, an average of at least 4 photons per time point) for an average of 1000 time points. From Fig. 9, we see that when molecules diffuse rapidly, that is, when D is on the order of 10 μm2/s, the (symmetrized) EKF provides a poor estimate of the diffusion coefficient due to the local linearization employed, while the UKF provides a reasonable result due to its global approximation.

FIG. 9.

Posterior results in 3D at different diffusion coefficients. [(a1)–(a3)] Sample traces illustrating the number of photons collected over time by a molecule with diffusion coefficients of 0.1, 1, and 10 μm2/s, respectively. [(b1) and (b2)] The posterior probability distribution over the diffusion coefficient using EKFs applied to the traces in (a). [(c1) and (c2)] The posterior probability distribution over the diffusion coefficient using UKFs applied to the traces in (a). See the caption of Fig. 8 for details on the parameters used.

FIG. 9.

Posterior results in 3D at different diffusion coefficients. [(a1)–(a3)] Sample traces illustrating the number of photons collected over time by a molecule with diffusion coefficients of 0.1, 1, and 10 μm2/s, respectively. [(b1) and (b2)] The posterior probability distribution over the diffusion coefficient using EKFs applied to the traces in (a). [(c1) and (c2)] The posterior probability distribution over the diffusion coefficient using UKFs applied to the traces in (a). See the caption of Fig. 8 for details on the parameters used.

Close modal

We now extend our 3D model to allow for different illumination profiles.2 In particular, we consider a cylindrical PSF shape (i.e., a major semi-axis far exceeding in size the minor semi-axis for the 3D Gaussian) in addition to the ellipsoidal one (a 3D Gaussian with major and minor semi-axes of the same order) previously considered. That is, for a cylindrical volume, we consider the case z/ωz ≪ 1, and thus, the PSF reduces to

P S F ( x , y , z ) = exp 2 x 2 ω x y 2 + y 2 ω x y 2 .
(31)

Since, for a cylindrical shaped volume, the observation probability is the same for any location of the molecule along the z-axis, now it is only possible to measure the radial distance of the molecule with respect to the z-axis. As such, for a molecule at (xk, yk, zk), we can estimate trajectories in terms of the following distance:

R k = x k 2 + y k 2 ,
(32)

where, Rk is the radial distance of the molecule respect to the z-axis.

Figure 10 shows that, even in this new geometry, we can capture the normalized distance of the molecule with respect to the z-axis (Rk) as well as its diffusion coefficient.

FIG. 10.

Posterior results in 3D under a cylindrical PSF. (a) A sample trace illustrating the number of photons collected over time. [(b1) and (b2)] Trajectory estimates using EKFs and UKFs, respectively, applied to the trace shown in (a). [(c1) and (c2)] The posterior probability distribution over the diffusion coefficient using EKFs and UKFs, respectively, applied to the traces in (a). The ground truth of the diffusion coefficient is shown by dashed lines. More details are provided in the caption of Fig. 8.

FIG. 10.

Posterior results in 3D under a cylindrical PSF. (a) A sample trace illustrating the number of photons collected over time. [(b1) and (b2)] Trajectory estimates using EKFs and UKFs, respectively, applied to the trace shown in (a). [(c1) and (c2)] The posterior probability distribution over the diffusion coefficient using EKFs and UKFs, respectively, applied to the traces in (a). The ground truth of the diffusion coefficient is shown by dashed lines. More details are provided in the caption of Fig. 8.

Close modal

Although both filters provide reasonably good estimates of the trajectory, the UKF again provides better estimates of the diffusion coefficient than the EKF.

In this study, we introduce a new framework to decode spatial information from the most common illumination profiles used across confocal microscopy.2 As the confocal volumes exhibit symmetry, i.e., multiple molecular positions generate identical signals, this creates modeling and computational challenges. To address these, as well as to address the inherent non-linearities in the observation process introduced by the PSFs dependence on the molecular location, we develop new approaches for non-linear filtering. In particular, we developed multi-modal generalizations of both extended and unscented Kalman filters.51–58 With these filters, we estimate not only molecular trajectories and diffusion coefficients but also errors around all of our estimates as, within the Bayesian paradigm, we obtain full posteriors over all quantities of typical interest in experiments employing confocal microscopy.

Our novel framework enables us, for the first time, to achieve spatial resolution from a conventional single-focus confocal setup. Although several experimental methods derive spatial information from confocal microscopy,18,63,64 all such methods rely on complex optical setups with multiple confocal volumes. It is now possible to generalize the approach presented here to treat multi-confocal volumes and even generalize it to the case of complex photophysics65 or multiple molecules within the illuminated region. The latter is a critical generalization as, if more than one molecule is present in the confocal volume, using our method unmodified would result in an estimate of only an average trajectory (as opposed to separate trajectories for each molecule). In a companion study,66 we show that the treatment presented herein, under the appropriate extensions, indeed carries over many molecules diffusing simultaneously through the confocal volume.

Our current formulation can also be adapted for different physical scenarios, for different experimental setups, or for different biophysical/biochemical phenomena. As an example, in cases that the diffusive molecule is also influenced by drift, our motion model can not only be modified to incorporate drift but also to learn the drift velocity. Drift velocity, in turn, could resolve the symmetry problem (i.e., our inability to distinguish in which portion of the shell of a Gaussian PSF a molecule is located). Furthermore, we can now envision treating more complex PSFs32 or more complex motion models as well, in particular, cases where the diffusion coefficient itself is changing over time.67 Also, using the Stokes-Einstein equation,33 we could extract the temperatures through the diffusion coefficients.

As a comment on possible broader implications of this work, we note that, traditionally, in fluorescence microscopy, signals are produced using multi-pixel cameras that receive photons at different rates. Accordingly, images are only subsequently formed by interpreting differences among the signals recorded in each pixel. Therefore, naive interpretation of the acquired data suggest that just one pixel in fluorescence microscopy should not produce images or should not provide any spatial resolution at all since it lacks nearby pixels to compare against. Nevertheless, this is true only under the assumption of uniform sample illumination.

Yet confocal data, often acquired to investigate processes with higher temporal resolution, encode spatial information even if just one pixel, for example, a sensor-like APD, is used. This is fundamental because the sample is under non-uniform illumination and so spatial information is already encoded on the photon emission rate itself.

Put differently, signals in fluorescence microscopy may be generated by signal differences either present at the observation level (i.e., using multi-pixel cameras with different signals arriving in each pixel) or at the data-generation level (i.e., using non-uniform illumination). Here, we take a first step in the latter direction suggesting a strategy for achieving some degree of spatial resolution—while fully retaining the temporal resolution—of confocal methods from just one pixel.

See supplementary material for a detailed derivation of EKF and UKF formulas, a description of the exact sampler, and a summary of our notational convention. The source code of the model is also available.

S.P. acknowledges support from NSF CAREER Grant No. MCB-1719537.

1.
J. C.
Stockert
and
A.
Blazquez-Castro
,
Fluorescence Microscopy in Life Sciences
(
Bentham Science Publishers
,
2017
).
2.
J.
Art
, in
Handbook of Biological Confocal Microscopy
(
Springer
,
2006
) pp.
251
264
.
3.
M. A.
Hayat
 et al,
Principles and Techniques of Electron Microscopy. Biological Applications
(
Edward Arnold
,
1981
).
4.
B.
Voigtländer
,
Scanning Probe Microscopy
(
Springer
,
2016
).
5.
A.
Lee
,
K.
Tsekouras
,
C.
Calderon
,
C.
Bustamante
, and
S.
Pressé
,
Chem. Rev.
117
,
7276
(
2017
).
6.
F. W.
Rost
,
Fluorescence Microscopy
(
Cambridge University Press
,
1992
), Vol. 2.
7.
M. J.
Rust
,
M.
Bates
, and
X.
Zhuang
,
Nat. Methods
3
,
793
(
2006
).
8.
E.
Betzig
,
Science
313
,
1642
(
2006
).
9.
S. T.
Hess
,
T. P.
Girirajan
, and
M. D.
Mason
,
Biophys. J.
91
,
4258
(
2006
).
10.
K.
Jaqaman
,
D.
Loerke
,
M.
Mettlen
,
H.
Kuwata
,
S.
Grinstein
,
S. L.
Schmid
, and
G.
Danuser
,
Nat. Methods
5
,
695
(
2008
).
11.
M. J.
Saxton
and
K.
Jacobson
,
Annu. Rev. Biophys. Biomol. Struct.
26
,
373
(
1997
).
12.
H.
Shen
,
L. J.
Tauzin
,
R.
Baiyasi
,
W.
Wang
,
N.
Moringo
,
B.
Shuang
, and
C. F.
Landes
,
Chem. Rev.
117
,
7331
(
2017
).
13.
I.
Sgouralis
,
A.
Nebenführ
, and
V.
Maroulas
,
SIAM J. Imaging Sci.
10
,
871
(
2017
).
14.
A. J.
Theuwissen
,
Solid-State Imaging with Charge-Coupled Devices
(
Springer Science and Business Media
,
2006
), Vol. 1.
15.
F.
Huang
,
T. M.
Hartwich
,
F. E.
Rivera-Molina
,
Y.
Lin
,
W. C.
Duim
,
J. J.
Long
,
P. D.
Uchil
,
J. R.
Myers
,
M. A.
Baird
,
W.
Mothes
 et al,
Nat. Methods
10
,
653
(
2013
).
16.
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
,
Nano Lett.
10
,
4732
(
2010
).
17.
K.
Welsher
and
H.
Yang
,
Nat. Nanotechnol.
9
,
198
(
2014
).
18.
E. P.
Perillo
,
Y.-L.
Liu
,
K.
Huynh
,
C.
Liu
,
C.-K.
Chou
,
M.-C.
Hung
,
H.-C.
Yeh
, and
A. K.
Dunn
,
Nat. Commun.
6
,
7874
(
2015
).
19.
M. F.
Juette
and
J.
Bewersdorf
,
Nano Lett.
10
,
4657
(
2010
).
20.
J.
Icha
,
M.
Weber
,
J. C.
Waters
, and
C.
Norden
,
BioEssays
39
,
1700003
(
2017
).
21.
V.
Magidson
and
A.
Khodjakov
,
Methods in Cell Biology
(
Elsevier
,
2013
), Vol. 114, pp.
545
560
.
22.
N.
Chenouard
,
I.
Smal
,
F.
De Chaumont
,
M.
Maška
,
I. F.
Sbalzarini
,
Y.
Gong
,
J.
Cardinale
,
C.
Carthel
,
S.
Coraluppi
,
M.
Winter
 et al,
Nat. Methods
11
,
281
(
2014
).
23.
U.
Von Toussaint
,
Rev. Mod. Phys.
83
,
943
(
2011
).
24.
R.
Rigler
and
E. S.
Elson
,
Fluorescence Correlation Spectroscopy: Theory and Applications
(
Springer Science and Business Media
,
2012
), Vol. 65.
25.
T. A.
Laurence
, “
Photon-counting single-molecule spectroscopy for studying conformational dynamics and macromolecular interactions
,” Ph.D. thesis,
University of California
,
Berkeley
,
2002
.
26.
H.
Qian
and
E. L.
Elson
,
Appl. Opt.
30
,
1185
(
1991
).
27.
B.
Zhang
,
J.
Zerubia
, and
J.-C.
Olivo-Marin
,
Appl. Opt.
46
,
1819
(
2007
).
28.
E.
Wolf
, in
Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences
(
The Royal Society
,
1959
), Vol. 253, pp.
349
357
.
29.
B.
Richards
and
E.
Wolf
, in
Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences
(
The Royal Society
,
1959
), Vol. 253, pp.
358
379
.
30.
R.
Rigler
,
Ü.
Mets
,
J.
Widengren
, and
P.
Kask
,
Eur. Biophys. J.
22
,
169
(
1993
).
31.
M.
Born
and
E.
Wolf
,
Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light
(
Elsevier
,
2013
).
32.
S. F.
Gibson
and
F.
Lanni
,
J. Opt. Soc. Am. A
9
,
154
(
1992
).
33.
H. C.
Berg
,
Random Walks in Biology
(
Princeton University Press
,
1993
).
34.
T.
Hida
,
Brownian Motion
(
Springer
,
1980
), pp.
44
113
.
35.
A.
Einstein
,
Ann. Phys.
322
,
549
(
1905
).
36.
A.
Gelman
,
J. B.
Carlin
,
H. S.
Stern
,
D. B.
Dunson
,
A.
Vehtari
, and
D. B.
Rubin
,
Bayesian Data Analysis
(
CRC Press
,
Boca Raton, FL
,
2014
), Vol. 2.
37.
M.
Tavakoli
,
J. N.
Taylor
,
C.-B.
Li
,
T.
Komatsuzaki
, and
S.
Pressé
, “
Single molecule data analysis: An introduction
,” in
Advances in Chemical Physics, Volume 162
, edited by
S. A.
Rice
and
A. R.
Dinner
(
John Wiley & Sons
,
2017
), pp. 205-305.
38.
C. M.
Bishop
,
Pattern Recognition and Machine Learning
(
Springer
,
2006
).
39.
I.
Sgouralis
and
S.
Pressé
,
Biophys. J.
112
,
2021
(
2017
).
40.
L.
Rabiner
and
B.
Juang
,
IEEE ASSP Mag.
3
,
4
(
1986
).
41.
I.
Sgouralis
and
S.
Pressé
,
Biophys. J.
112
,
2117
(
2017
).
42.
C.
Robert
and
G.
Casella
,
Introducing Monte Carlo Methods with R
(
Springer Science and Business Media
,
2009
).
43.
O.
Cappé
,
E.
Moulines
, and
T.
Rydén
, in
Proceedings of EUSFLAT Conference
(
Springer
,
2009
), pp.
14
16
.
44.
T.
Rydén
 et al,
Bayesian Anal.
3
,
659
(
2008
).
45.
L. R.
Rabiner
,
Proc. IEEE
77
,
257
(
1989
).
46.
S. L.
Scott
,
J. Am. Stat. Assoc.
97
,
337
(
2002
).
47.
F. J.
Anscombe
,
Biometrika
35
,
246
(
1948
).
48.
R. E.
Kalman
,
J. Basic Eng.
82
,
35
(
1960
).
49.
R. E.
Kalman
and
R. S.
Bucy
,
J. Basic Eng.
83
,
95
(
1961
).
50.
H. W.
Sorenson
,
Advances in Control Systems
(
Elsevier
,
1966
), Vol. 3, pp.
219
292
.
51.
R.
Frühwirth
,
Nucl. Instrum. Methods Phys. Res., Sect. A
262
,
444
(
1987
).
52.
L.
Ljung
,
IEEE Trans. Autom. Control
24
,
36
(
1979
).
53.
T.
Song
and
J.
Speyer
,
IEEE Trans. Autom. Control
30
,
940
(
1985
).
54.
M.
Hoshiya
and
E.
Saito
,
J. Eng. Mech.
110
,
1757
(
1984
).
55.
M. Y.
Byron
,
K. V.
Shenoy
, and
M.
Sahani
, technical report, Department of Electrical Engineering, Stanford University, 2004.
56.
E. A.
Wan
and
R.
Van Der Merwe
, in
Adaptive Systems for Signal Processing, Communications, and Control Symposium 2000. AS-SPCC. The IEEE 2000
(
IEEE
,
2000
), pp.
153
158
.
57.
H. M.
Menegaz
,
J. Y.
Ishihara
,
G. A.
Borges
, and
A. N.
Vargas
,
IEEE Trans. Autom. Control
60
,
2583
(
2015
).
58.
J.
Stoer
and
R.
Bulirsch
,
Introduction to Numerical Analysis
(
Springer Science and Business Media
,
2013
), Vol. 12.
59.
O. C.
Ibe
,
Elements of Random Walk and Diffusion Processes
(
John Wiley & Sons
,
2013
).
60.
J.
Haile
,
I.
Johnston
,
A. J.
Mallinckrodt
,
S.
McKay
 et al,
Comput. Phys.
7
,
625
(
1993
).
61.
D. J.
Higham
,
SIAM Rev.
43
,
525
(
2001
).
62.
R.
Erban
and
S. J.
Chapman
,
Phys. Biol.
6
,
046001
(
2009
).
63.
J. A.
Germann
and
L. M.
Davis
,
Opt. Express
22
,
5641
(
2014
).
64.
L.
Lanzanò
and
E.
Gratton
,
Methods Appl. Fluoresc.
2
,
024010
(
2014
).
65.
I.
Sgouralis
,
S.
Madaan
,
F.
Djutanta
,
R.
Kha
,
R. F.
Hariadi
, and
S.
Pressé
,
J. Phys. Chem. B
123
,
675
(
2018
).
66.
S.
Jazani
,
I.
Sgouralis
,
O. M.
Shafraz
,
S.
Sivasankar
, and
S.
Pressé
,
bioRxiv
2018
,
426114
.
67.
S.
Thapa
,
M. A.
Lomholt
,
J.
Krog
,
A. G.
Cherstvy
, and
R.
Metzler
,
Phys. Chem. Chem. Phys.
20
,
29018
(
2018
).

Supplementary Material