An alternative framework for fluorescence correlation spectroscopy.

Experimental data with elongated confocal volume shapes

Here, we apply our method on experimental traces captured with an elongated confocal volume that we approximate by a cylinder. To do so, we apply our method on fluorescent beads (with average diameter of 45 nm) diffusing in water. We benchmark our estimated diffusion coefficients against the Stokes–Einstein prediction and results from FCS. In particular, Fig. 4 illustrates our method’s performance in the analysis of traces too short to be meaningfully analyzed by FCS. The FCS formulation we used here can be found in Supplementary Note 3 and additional results can be found in Supplementary Fig. 10 . Fig. 4 Traces of free fluorescent beads using an elongated confocal volume. a Experimental fluorescent intensity trace used in FCS. b Autocorrelation curve of the trace in a and best theoretical fit. c Portion of the trace in a to be used as the input to FCS and our method. d Autocorrelation curve of trace in c . e Posterior probability distribution over the diffusion coefficient estimated from the trace in c . The Stokes–Einstein prediction is denoted by a green line. Traces shown in a and c are acquired at 100 μs for a total period of 5 μs and 0.1 μs, respectively. The laser power used to generate the trace a is 100 μW (measured before the beam enters the objective). The estimate of diffusion coefficient resulting by autocorrelation fitting in a matches with the Stokes–Einstein prediction (i.e., 10.5 μm 2 s −1 ); whereas, in d is almost 10-fold higher (~ 10 μm 2 s −1 )

Experimental data with elliptical confocal volume shapes

Next, we apply our method on experimental time traces derived from single-molecule fluorescence confocal microscopy. In our setup, we monitor Cy3 dyes, which diffuse freely in a mixture of water and glycerol. We benchmark our estimated diffusion coefficients against two values: those predicted by the Stokes–Einstein formula 39 , which is parametrized by physical quantities such as temperature and viscosity; and those estimated by FCS. To analyze the data using FCS reliably, we used the full (6 min) trace available. In benchmarking, we obtained and analyzed measurements with different: (i) numbers of Cy3 dyes inside the effective volume (tuned by varying Cy3 concentration); (ii) trace lengths; (ii) diffusion coefficients (tuned by adjusting the viscosity of the solution); and (iii) molecular brightness (tuned by adjusting the laser power). Just as before, the slower a molecule diffuses, the more time it spends in the vicinity of the confocal volume, so the more photons are collected, thereby leading to sharper posterior estimates for the diffusion coefficient; as seen on Fig. 5a . Fig. 5 Estimating diffusion coefficients of free Cy3. a Posterior probability distributions of diffusion coefficients of free Cy3 in different concentrations of glycerol/water mixture. The legend labels the posteriors according to FCS estimates of long time traces. For clarity, posteriors are normalized to maximum 1 and the horizontal axis is shown in logarithmic scale. Also, the 95% confidence intervals are shown by highlighted regions. Posteriors are obtained from the analyses of time traces acquired at 100 μs for total periods of 100 ms. Different diffusion coefficients are obtained by varying the amount of glycerol from 99 to 50% in the glycerol/water mixture. b Posterior probability distributions deduced from traces acquired at 100 μs with total trace lengths of 5 × 10 2 , 1 × 10 3 , 5 × 10 3 , 1 × 10 4 , and 5 × 10 4 time steps. For the sake of comparison, exact values and FCS estimates are also shown and, for clarity, the vertical axis is shown in logarithmic scale. Error bars in the FCS curve are produced by analyzing multiple windows in the initial trace. To estimate the diffusion coefficient within less than a factor of 2 of the true value, it is typical for FCS to require ≈ 50 × more data than our method. c – f Posterior probability distributions over the diffusion coefficients of traces generated by different laser powers (25 and 100 μW, respectively) with different concentrations of Cy3 (100 p m and 1 n m , respectively) in a glycerol/water mixture of 94% glycerol. For the sake of comparison, FCS estimates, shown with dashed lines, are obtained from traces each 6 min long In Fig. 5b , we illustrate the effect of different time trace length. In this case, we reach the same accuracy as FCS with 100 × times less data. Consistent with the synthetic data shown earlier, we obtain a broader posterior over diffusion coefficients when the number of dyes inside the effective volume is low and sharper posteriors for higher numbers of dyes. For example, in Fig. 5c–f , we illustrate the effect of different dye concentrations where a trace with stronger signal, anticipated when concentrations are higher, leads to better diffusion coefficient estimates (and thus sharper posteriors) on traces of equal length owing to the higher number of labeled molecules inside the confocal volume. The corresponding full joint posteriors of Fig. 5 are shown in Supplementary Figs 6 and 7 . In general, a posterior’s sharpness depends strongly on the number of molecules in a time trace, their respective locations, and thus their photon emission rates. As the molecular population near the center of the confocal volume may exhibit strong fluctuations, the width of the posterior may also fluctuate from trace to trace, especially when the individual traces are short. Thus, the individual posteriors become sharper only on average as we move to higher numbers of molecules inside the effective volume or molecular brightness. To test our method beyond free beads and dyes, we used labeled proteins, namely freely diffusing streptavidin labeled by Cy3. Similar to the previous cases, we tested a range of concentrations, diffusion coefficients, and laser powers. Figure 6 summarizes characteristic results and compares our analyses against the results of FCS, which is applied on longer time traces (6 min). As can be seen, even in this case our method provides acceptable estimates of the diffusion coefficient with 100 × times fewer data points than FCS. Fig. 6 Estimating diffusion coefficients of free streptavidin. a Posterior probability distributions of diffusion coefficients of free streptavidin labeled by Cy3 in different concentrations of glycerol/water mixture. The legend labels the posteriors according to FCS estimates of long time traces. For clarity, posteriors are normalized to maximum 1 and the horizontal axis is shown in logarithmic scale. Also, the 95% confidence intervals are shown by highlighted regions. Posteriors are obtained from the analyses of time traces acquired at 100 μs for total periods of 100 ms. Different diffusion coefficients are obtained by varying the amount of glycerol from 94 to 0% in the glycerol/water mixture. b – e Posterior probability distributions over the diffusion coefficients of traces generated by different laser powers (25 and 100 μW, respectively) with different concentrations of Cy3 (100 p m and 1 n m , respectively) in a glycerol/water mixture of 94% glycerol. For the sake of comparison, FCS estimates, shown by dashed lines, are obtained by traces each 6 min long

Methods Model overview

Here we describe the formulation and mathematical foundation of our model. Our overarching goal is to start from an experimental time series of photon counts, documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$bar w = (w_1,w_2,...,w_K)$$end{document} w ¯ = ( w 1 , w 2 , . . . , w K ) where w k denotes the photon intensity assessed at time t k (which includes both background photons as well as photons derived from the labeled molecules of interest), and derive estimates of kinetic quantities such as molecular locations with respect to the center of the confocal volume as well as diffusion coefficients. To derive estimates for the desired quantities, we need to compute intermediate quantities which include: (i) molecular brightness; (ii) background photon emission rate; and, most importantly, (iii) the unknown population of moving molecules and their relative locations with respect to the center of the confocal volume. Below we explain each one of these in detail. Computational details and a working implementation of the entire method are available in the Supplementary Notes 3 and 4 . For convenience, we summarize our notation, abbreviations and mathematical definitions in Supplementary Tables 2 – 4 . Model description The starting point of our analysis is the raw data, namely the photon counts. As our current focus is on deducing dynamical information on timescales exceeding ≈ 1 μs, we ignore triplet state and photon anti-bunching effects that occur on vastly different timescales 16 , 55 , 56 . At the timescale of interest, individual photon detections, assuming saturation is not reached, happen stochastically and independently from each other. Accordingly, the total number of photon counts w k between successive assessments follows Poisson 15 , 27 (shot noise) statistics 1 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$begin{array}{*{20}{l}} {w_ksim {mathrm{Poisson}}left( {(t_k - t_{k - 1})left( {mu _{back} + mathop {sum}limits_n {mu _k^n} } right)} right)} hfill end{array}$$end{document} w k ~ Poisson ( t k - t k - 1 ) μ b a c k + ∑ n μ k n where μ back is a background photon emission rate and documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$mathop {sum}limits_n {mu _k^n}$$end{document} ∑ n μ k n gathers the photon emission rates documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$mu _k^n$$end{document} μ k n from individual fluorescent molecules that we index with n = 1, 2, …. The number of molecules involved in the summation above is to be determined. This is the key reason we invoke Bayesian non-parametrics in the model inference section (see below). As we only collect a small fraction of the total photons emitted by the fluorescent molecules, as we describe above, in our framework documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$mu _k^n$$end{document} μ k n coincides with the emission rate of detected photons, as opposed to the true photon emission rate, which might be larger. Each rate documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$mu _k^n$$end{document} μ k n depends on the position documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$(x_k^n,y_k^n,z_k^n)$$end{document} ( x k n , y k n , z k n ) of the corresponding molecule relative to the center of the confocal volume as well as other features such as laser intensity, laser wavelength, quantum yield, and camera pinhole size 57 . Similar to other studies 40 , 58 , 59 , we combine all these effects into a characteristic point spread function (PSF) that combines excitation and emission PSFs 2 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$mu _k^n = mu _{{mathrm{mol}}}{mathrm{PSF}}(x_k^n,y_k^n,z_k^n).$$end{document} μ k n = μ mol PSF ( x k n , y k n , z k n ) . The parameter μ mol represents the molecular brightness and, as we discuss in the Supplementary Note 3 , it is related to the maximum photon emission rate of a single molecule that is located at the center of the confocal volume. Specific choices of PSF models, such as Gaussian or Gaussian-Lorentzian, are also detailed in the Supplementary Note 3 and a comparison of the different PSF models is shown in Supplementary Fig. 5 . Finally, we associate individual molecular locations across time by adopting a motion model. Here we assume that molecules are purely diffusive and arrive at 3 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$begin{array}{l}x_k^nsim {mathrm{Normal}}left( {x_{k - 1}^n,2(t_k - t_{k - 1})D} right)\ y_k^nsim {mathrm{Normal}}left( {y_{k - 1}^n,2(t_k - t_{k - 1})D} right)\ z_k^nsim {mathrm{Normal}}left( {z_{k - 1}^n,2(t_k - t_{k - 1})D} right)end{array}$$end{document} x k n ~ Normal x k - 1 n , 2 ( t k - t k - 1 ) D y k n ~ Normal y k - 1 n , 2 ( t k - t k - 1 ) D z k n ~ Normal z k - 1 n , 2 ( t k - t k - 1 ) D where D denotes the diffusion coefficient, which we assume is the same for all molecules. As we explain in the Supplementary Note 4 , these probabilities result directly from the diffusion equation. A graphical summary of the entire formulation is shown on Fig. 8 . Fig. 8 Graphical representation of the formulation. A population of model molecules, labeled by n = 1, 2,…, evolves over the course of the experiment that is marked by k = 1, 2, …, K . Here, documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$x_k^n$$end{document} x k n , documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$y_k^n$$end{document} y k n , documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$z_k^n$$end{document} z k n denote the location, in Cartesian space, of molecule n at time t k ; μ mol denotes the brightness of an individual molecule; and μ back denotes the background photon emission rate. During the experiment, only a single observation w k , combining photon emissions between t k −1 and t k from every molecule and background is recorded at every time step. The diffusion coefficient D determines the evolution of the molecular locations which, in turn, influence the photon emission rates and ultimately the recorded photon intensity w k . Auxiliary variables b n , or “loads”, and corresponding prior weights q n , are introduced in order to estimate the unknown population size. The dashed arrows apply for the 3D Gaussian and 2D-Gaussian-Lorentzian PSFs; while, in the case of the 2D-Gaussian-Cylindrical, there is no dependency of the measurements w k on the documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$z_k^n$$end{document} z k n coordinates of the molecules (see the Supplementary Note 3 for the definitions of the PSFs) In addition, in the Supplementary Note 5 and Supplementary Fig. 16 , we illustrate how this motion model can be generalized to capture more than one diffusion coefficients.

Model inference

The quantities that we want to estimate, for example the diffusion coefficient D , molecular locations through time documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$(x_k^n,y_k^n,z_k^n)$$end{document} ( x k n , y k n , z k n ) , molecular brightness μ mol and background photon emission rate μ back , or the molecular population, are introduced as model variables in the preceding formulation. To estimate values for these variables, we follow the Bayesian paradigm 15 , 28 , 38 , 59 . Variables such as D , μ mol , and μ back are parameters of the model and, as such, require priors. Choices for these priors are straightforward and, for interpretational and computational convenience, we adopt the distributions described in the Supplementary Note 4 . In addition, we must place priors on the initial molecular locations, documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$(x_1^n,y_1^n,z_1^n)$$end{document} ( x 1 n , y 1 n , z 1 n ) , i.e., the locations of the molecules at the onset of the measurement period. Specifying a prior on initial molecular locations also entails specifying a prior on the molecular population. In particular, to allow the dimensionality or, alternatively, the complexity of our model to fluctuate based on the number of molecules that contribute to the fluorescent trace, we abandon traditional Bayesian parametric priors and turn to the non-parametric formulation described below. Before we proceed any further, we recast Eq. ( 2 ) as 4 documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$mu _k^n = b^nmu _{{mathrm{mol}}}PSF(x_k^n,y_k^n,z_k^n).$$end{document} μ k n = b n μ mol P S F ( x k n , y k n , z k n ) . The newly introduced variables b n , one for each model molecule, may take only values 1 or 0. In particular, the possibility that b n = 0, coinciding with the case where molecules do not contribute to the observation, allows us to introduce an arbitrarily large number of molecules, technically an infinite number. With the introduction of b n , we can estimate the number of molecules that contribute photons (termed “active” to distinguish them from those that do not contribute termed “inactive”) simultaneously with the rest of the parameters simply by treating each b n as a separate parameter and estimating its value (of 1 for active molecules and 0 for inactive ones). To estimate b n , we place a prior b n ~ Bernoulli( q n ) and subsequently a hyperprior on q n in order to learn precisely how many model molecules are active. For the latter, we choose q n ~ Beta( A q , B q ) with hyperparameters A q and B q . Both steps can be combined by invoking the newly developed Beta-Bernoulli process 36 , 60 which is described in more detail in the Supplementary Note 4 . Once the choices for the priors above are made, we form a joint posterior probability distribution documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$$p(D,mu _{{mathrm{mol}}},mu _{{mathrm{back}}},{ x_k^n,y_k^n,z_k^n,b^n,q^n} _k^n|bar w)$$end{document} p ( D , μ mol , μ back , { x k n , y k n , z k n , b n , q n } k n ∣ w ¯ ) encompassing all unknown variables which we may wish to determine. The nonlinearities in the PSF with respect to variables documentclass[12pt]{minimal} usepackage{amsmath} usepackage{wasysym} usepackage{amsfonts} usepackage{amssymb} usepackage{amsbsy} usepackage{mathrsfs} usepackage{upgreek} setlength{oddsidemargin}{-69pt} begin{document}$${ x_k^n,y_k^n,z_k^n} _k^n$$end{document} { x k n , y k n , z k n } k n and the non-parametric prior on { b n , q n } n exclude analytic forms for our posterior. For this reason, we develop a computational scheme exploiting Markov chain Monte Carlo 38 , 61 that can be used to generate pseudo-random samples from this posterior. The main bottleneck of a naive implementation of our method, as compared with correlative methods, is its higher computational cost. As we explain in the Supplementary Note 4 , to have computations run on an average desktop computer, we adopt mathematical approximations (e.g., photon binning, Anscombe transform 62 and filter updates 63 , 64 ) that are tested on the synthetic data presented. Specifically, the time trace preparation is described in the Supplementary Note 3 and Supplementary Fig. 14 . A working implementation of the framework described in this study is provided in the source code and the graphical user interface (GUI) is shown on Supplementary Fig. 15 .

Data acquisition

Synthetic data: We obtain the synthetic data presented in the Results section by standard pseudo-random computer simulations 65 – 67 that mimic the common single-molecule fluorescence confocal setup. We provide details and complete parameter choices in the Supplementary Note 3 , Supplementary Tables 5 and 6 . Experimental data: For the experimental data acquired with elongated confocal volumes, a stock solution of Cy3B (mono-reactive NHS ester, GE Healthcare) was prepared by dissolving a small amount of solid in 1 mL of doubly-distilled water, and its concentration was determined from the absorbance of the solution using the extinction coefficient provided by the vendors. A 10 n m solution was then prepared by appropriate dilution of the stock and measured on a silicone perfusion chamber mounted on a glass coverslip. Fluorescent beads were purchased from ThermoFisher (Catalog number: F8792. Lot number: 1604237). The average diameter was 0.046 μm as indicated in the certificate of analysis provided by the vendors. Suspensions for FCS measurements were prepared by adding 3 μL of stock solution (9.4 × 1014 particles/mL) to 1 mL of water and sonicating the mixture for 20 min. Measurements were carried out using a home-built instrument. A 532 nm continuous-wave laser (Compass 215M-10, Coherent, Santa Clara, CA) was attenuated to 100 μW and focused onto an PlanApo 100 ×, 1.4 NA, oil-immersion, objective (Olympus, Center Valley, PA). Emitted fluorescence was collected using the same objective and then passed through a 50 μm pinhole to reject the out-of-focus light. The signal was detected using a silicon avalanche photodiode (SPCM-AQR-14; Perkin-Elmer, Fremont, CA). A bandpass filter (Omega 3RD560-620) in front of the detector was employed to further reduce the background signal and an ALV correlator card (ALV 5000/EPP, ALV-GmbH, Langen, Germany) was used to correlate the detected fluorescence signal. Data for our analysis were acquired with 100 μs resolution using a PCI-6602 acquisition card (National Instruments, Austin, TX). Measurements have been tested for saturation separately and shown on Supplementary Fig. 11 . For the experimental data acquired with elliptical confocal volumes, Cy3 dye and Cy3-labeled streptavidin solutions were prepared by suspending Cy3 or streptavidin in glycerol/buffer (pH 7.5, 10 m m Tris-HCl,100 m m NaCl and 10 m m KCl, 2.5 m m CaCl 2 ) at different v/v, to a final concentration of either 100 p m or 1 n m . The solutions were added onto a glass-bottomed fluid-cell, mounted on a custom designed single-molecule fluorescence confocal microscope 68 , 69 and a 532 nm laser beam was focused to a diffraction-limited spot on the glass coverslip of the fluid-cell using a × 60, 1.42 NA, oil-immersion objective (Olympus). The laser power was measured before the objective and the beam was reflected by a dichroic and focused by the objective on to the sample. The dichroic reflected 95% of the intensity on to the objective. Emitted fluorescence was collected by the same objective and focused onto the detection face of a Single Photon Avalanche Diode (SPAD, Micro Photon Devices) that has a maximum count rate of 11.8 Mc/s. A bandpass filter was placed in front of the detector to transmit only the fluorescence from Cy3 and to block the back-scattered excitation light. TTL pulses, triggered by the arrival of individual photons on the SPAD, were timestamped and recorded at 80 MHz by a field programmable gated array (FPGA, NI Instruments) using custom LabVIEW software 69 and initially binned at 100 μs.

Supplementary information Supplementary Information