# Local temporal image correlation spectroscopy and Bayesian simulation technique for sparse estimation of diffusion in fluorescence imaging

Anca Caranfil, Yann Le Cunff, Charles Kervrann

► To cite this version:

Anca Caranfil, Yann Le Cunff, Charles Kervrann. Local temporal image correlation spectroscopy and Bayesian simulation technique for sparse estimation of diffusion in fluorescence imaging. 2022. [hal-03540375](https://hal HAL Id: hal-03540375)

https://inria.hal.science/hal-03540375

Preprint submitted on 23 Jan 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Local temporal image correlation spectroscopy and Bayesian simulation technique for sparse estimation of diffusion in fluorescence imaging

Anca Caranfil<fnref n="1" />,<fnref n="2" />, Yann Le Cunff<fnref n="2" />, and Charles Kervrann<fnref n="1" />

**ABSTRACT:** The dynamics and fusion of vesicles during the last steps of exocytosis are not well established yet in cell biology. An open issue is the characterization of the diffusion process at the plasma membrane. Total Internal Reflection Fluorescence Microscopy (TIRFM) has been successfully used to analyze the coordination of proteins involved in this mechanism. It enables to capture dynamics of proteins with high frame rate and reasonable SNR values. Nevertheless, we lack methodological approaches that can analyze and estimate diffusion in local small area at the scale of a single diffusing spot within cells. Accordingly, we propose a novel correlation-based method for local diffusion estimation. s starting point, we consider the diffusion model well established in biophysics. We derive an explicit parametric model which is further fitted to time-correlation signals computed from regions of interest (ROI) containing individual spots. This modeling and the proposed Bayesian estimation framework is well appropriate to represent isolated diffusion events and is robust to noise, ROI sizes and localizaiton of spots in ROIs. In contrast to previous fluorescence correlation methods focused on the characterization of the mobility of a population of molecules, no pre-treatment and little data are required in our approach. The performance of our method is shown on both synthetic and real TIRM images depicting Transferrin Receptor proteins.

## INTRODUCTION

There is a growing interest for the use of live cell imaging for cancer research and analysis of the dynamical behavior of biological structures. In cell biology, diffusion measurements are commonly used to compare several compartments within a cell or in different cells (i.e. neuronal synapses, Hela cells, etc...). At the scale of a single cell, intracellular dynamics may correspond to active transport of proteins but the main mode of transport is often free diffusion, i.e. proteins undergo Brownian motion. Several techniques have been developed to quantify the dynamics of fluorescently-tagged (i.e., Green Fluorescent protein (GFP)) proteins and to estimate diffusion coefficients from microscopy data. Single-particle tracking (SPT) approaches are especially recommended as the trajectory of a single molecule provide a quantitative description of motion in space and time. Given tracks obtained by nearest neighborhood algorithms or more sophisticated tracking methods [8], the mean-square displacement (MSD) of tracks is generally used to interpret and detect free diffusion, confined diffusion and directed flow [17]. The analysis is based on the slope of MSD curves or over time. The theoretical limits of SPT in particular for confined diffusion have been established in [5]. In practical experiments, this quantitative analysis is reliable for long tracks (i.e., tracks longer than 10 points) or for a set of short and homogeneous tracks.

Moreover, the applicability of SPT-based approaches is subject to several limitations, including labeling of single molecules with suitably high signal/noise marker particles and tracking errors with imperfect algorithms.

A second approach consists in analyzing the fluorescence recovery after its photobleaching in an specified area by a high intensity laser pulse. Fluorescence recovery after photobleaching (FRAP) analysis is performed directly by measuring the mean intensity in the photo-bleached region over the image sequence. The extracted fluorescence recovery curves are normalized to 1 for the pre-bleached intensity. Furthermore and all the recovery curves are fit to the theoretical model by using a nonlinear least square algorithm [18]. When these techniques are applied to estimate diffusion in large cell compartments (e.g. cytoplasm), confined diffusion is neglected for the range size (10µm to 1µm).

Fluorescence correlation spectroscopy (FCS) [9], [12] may be also considered as a valuable approach to provide quantitative information on single-molecule dynamics without individual object tracking. FCS generally analyzes temporal fluctuations of genetically encoded fluorescent molecules in a region of interest (ROI). It has been proved to provide estimates calculated from FCS curves, close to those obtained with SPT which needs individual molecule trajectories. FCS is a semi-local approach and assumes spatial and temporal stationarity of

<footnote n="1">¹ Inria Rennes, UMR144-CNRS, Institut Curie, PSL Research / Serpico Team-Project, 35042 Rennes Cedex France</footnote>

<footnote n="2">² IGDR UMR6290-CNRS, Faculty of Medicine, University of Rennes 1/ CeDRE Team, 35043 Rennes Cedex, France</footnote>
dynamics in the ROI.

Another popular family of methods for diffusion analysis is based on correlation measurements under the hypothesis of temporal stationarity of fluorescence signals [10], [11]. Temporal Image Correlation Spectroscopy (TICS) is derived from fluorescence correlation techniques [13]. This is a powerful method for diffusion determination, only based on post-processing of the image sequence and which does not require any molecule tracking. Spatio-Temporal Image Correlation Spectroscopy (STICS) is an extension of TICS and also widely used to recover physical parameters such as directional flow or diffusion parameters of moving molecules (see Appendix A in Supplementary Information). TICS and STICS are generally applied to small ROIs (or blocks of pixels of a regular grid) with a minimal size to analyze the sample particle displacements and to save time computing. More recently, 2D spatiotemporal cross correlation (pCF) of the fluorescence fluctuations between pairs of points was proposed in [14] to provide to accurately detect barriers of diffusion and to determine diffuse routes used by moving molecules. However, pCF needs several computational approximations to save memory and time and to provide a 2D anisotropic diffusion map from several dozens thousands of temporal frames acquired in TIRF microscopy with a frame rate 1ms/frame.

In this paper, we address the problem of estimating locally constant but spatially varying diffusion coefficient often representative of local change of the medium in time-lapse fluorescence microscopy. We propose a Bayesian framework (Approximate Bayesian Computing (ABC) [4]) for robustly estimating 2D local diffusion of fluorescently-tagged membrane proteins close to the plasma membrane. We show that the proposed method complements previous techniques by providing distinct advantages of local analysis over TICS counterparts on ROIs and non-overlapping pixel blocks. The main advantage of our Bayesian Temporal Correlation Spectroscopy (BayesTICS) method dedicated to the analysis of single diffusing spot are the following one: i/ formal analysis of diffusions from temporal correlation measurements; ii/ robust estimation over local ROIs. BayesTICS methods is robust to variable signal-to-noise ratios (SNR), the choice of the size of the analysis window, and the location the spot of interest in the ROI. Moreover, the method does not require accurate measurement of the point spread function (PSF) of the imaging system. Unlike previous methods based on non-linear model fitting [3], BayesTICS is based on simulation of time-correlation signals and provide the posteriori distribution of diffusion coefficient.

The remainder of the paper is organized as follows. In the next section, we describe the diffusion and image models with initial conditions, the closed-form solution and approximations, and the ABC-based algorithm to

estimate local diffusion coefficients from time-correlation signals computed over ROIs. In the experimental section, we demonstrate the robustness of BayesTICS to SNRs on simulated and real images depicting Transferin Receptor (TfR) dynamics at the plasma membrane during the late steps of exocytosis i.e. diffusion in the plasma membrane after vesicle fusion, as investigated in [9, 10, 11]. We evaluate the sensitivity of nuisance parameters (i.e., size of ROIs) and the algorithm parameters, and the accuracy of BayesTICS on simulated data. We believe that our approach further extends the capability of correlation spectroscopy methods and provides, to our knowledge, a new tool for the quantitative characterization of the spatiotemporal dynamics of molecules in living cells.

## METHOD

### Diffusion model and image generation

Diffusion of transmembrane proteins is mainly given by lateral diffusion in the membrane [1]. This phenomenon may be described by the Fick's second law [16]:

$$ \frac{\partial C(x, t)}{\partial t} = D \Delta C(x, t) \quad (1) $$

where $C(x, t)$ is the concentration of molecules at time $t$ and at location $x \in \Omega \subset \mathbb{R}^2$, $\Omega$ is the region of interest (or image domain), and $\Delta$ is the Laplacian operator. In our modeling framework, the diffusion coefficient $D$ in (1) is supposed to be constant over the support $\Omega$. In addition, we assume that all the proteins are concentrated at the center $x_0$ of $\Omega$, at initial time $t = t_0$, that is $C(x, t_0) = C_0\delta(x - x_0)$ where $\delta$ is the Kronecker symbol and $C_0$ is the initial concentration of molecules. With the aforementioned initial conditions, the closed form solution is well established and corresponds to:

$$ C(x, t) = \frac{C_0}{4\pi(t - t_0)D} \exp\left(-\frac{\|x - x_0\|_2^2}{4(t - t_0)D}\right), \quad \forall t > t_0. \quad (2) $$

Let $f = \{f(t)\}, t \in \{0, \dots, T - 1\}$ be a temporal series of fluorescence images acquired with TIRFM and depicting fluorescent spots. The intensity $f(t)$ is assumed is proportional (with factor $B$) to a convolution (denoted *) of the microscopic number density (or concentration) $C(t)$ and the instrumental point spread function $h$:

$$ f(t) = B(h * C(t)), \quad (3) $$

where $B = \rho\epsilon Q$, $\rho$ is the efficiency of the instrument to collect photons, $\epsilon$ is the molecular absorption coefficient, and $Q$ is the quantum yield of the fluorophore. Hence, the intensity $f(x, t)$ at time $t$ and location $x$ is then given by the convolution with the PSF function multiplied by a constant depending on the microscope and fluorophore. If we assume that the PSF is approximated with a 2D
Gaussian function with an isotropic bandwidth $\sigma_{PSF}$ in the
lateral direction, we get [3]:

$$
f(x, t) = \frac{A_0}{2D(t-t_0) + \sigma_{\text{PSF}}^2} \exp \left( - \frac{\|x-x_0\|_2^2}{4D(t-t_0) + 2\sigma_{\text{PSF}}^2} \right) \quad (4)
$$

with $A_0 = C_0 B / 2\pi$.

**Local Temporal Image Correlation Spectroscopy**

In this section, we describe our correlation-based method,
inspired from FCS, to compute local diffusion in fluo-
rescence images. First, we compute the temporal auto-
correlation of images to capture temporal fluctuations of
intensities.

**Uniform background** First, we are interested in locally estimating the diffusion of an isolated spot over a uniform background. In that case, the process associated with the image sequence is not stationary. Accordingly, the autocorrelation function is real time dependent and not only delay dependent; it is therefore impossible to compute the spatiotemporal average with only one realization of the underlying diffusison process. Therefore, we use the following expectation formula for the autocorrelation function:

$$
G_1(t, \tau) = \frac{\langle f(t) f(t + \tau) \rangle}{\langle f \rangle_t^2} \quad (5)
$$

where $\langle \rangle$ and $\rangle_t \rangle$ denote the spatial and spatio-temporal averages respectively, and $\tau$ is the temporal lag. By substituting the expression (4) to $f$ in (5), we get (see Appendix in Supplementary Information):

$$
G_1(t, \tau) \approx \frac{|\Omega|}{4\pi(D\tau + 2Dt + \sigma_{\text{PSF}}^2)}. \quad (6)
$$

Note that if consider the inverse of the autocorrelation
function, $G_1^{-1}(t, \tau)$ is a linear function of unknown pa-
rameters $D$ and $\sigma_{\text{PSF}}$:

$$
G_1^{-1}(t, \tau) \approx \frac{4\pi}{|\Omega|} D\tau + \frac{8\pi}{|\Omega|} Dt + \frac{4\pi}{|\Omega|} \sigma_{\text{PSF}}^2. \quad (7)
$$

**Non-uniform and cluttered background** In order to account for non-uniform background and the potential presence of neighboring spots diffusing in the region of interest $\Omega$, let us assume that the spatio-temporal average $\langle f \rangle_t$ is non-zero. Hence, we can use the following expectation formula for the intensity fluctuation autocorrelation function:

$$
G_2(t, \tau) = \frac{\langle \delta f(t) \delta f(t + \tau) \rangle}{\langle f \rangle_t^2} \quad (8)
$$

where $\delta f(t) = f(t) - \langle f \rangle_t$ is the fluctuation in fluores-
cence intensity around the temporal average value. T By
replacing $f$ in (8) with the expression given in (4), we
obtain:

$$
G_2(t, \tau) \approx G_1(t, \tau) + |\Omega|(K_1(t, \tau) + K_2(t)), \quad (9)
$$

with

$$
\begin{align*}
K_1(t, \tau) &= \frac{1}{4\pi D(T-t)} \ln \left( \frac{D\tau + 2Dt + \sigma_{\text{PSF}}^2}{D\tau + D(T+t) + \sigma_{\text{PSF}}^2} \right) \\
K_2(t) &= \frac{2DT + \sigma_{\text{PSF}}^2}{4\pi D^2(T-t)^2} \ln \left( \frac{2DT + \sigma_{\text{PSF}}^2}{D(T+t) + \sigma_{\text{PSF}}^2} \right) + \\
&\qquad \frac{D(T+t) + \sigma_{\text{PSF}}^2}{4\pi D^2(T-t)^2} \ln \left( \frac{2Dt + \sigma_{\text{PSF}}^2}{D(T+t) + \sigma_{\text{PSF}}^2} \right).
\end{align*}
$$

The model $G_2(t, \tau)$ can be considered as an extension of
the model $G_1(t, \tau)$.

In (6) and (9), the ROI size $|\Omega|$ is a multiplicative factor and can be interpreted here as a scaling parameter. We shall see that the ROI size does not influence the estimation of $D$ in our experiments provided there is only spot in the ROI. Unlike (6), the total time of observation $T$ appears in $K_1(t, \tau)$ and $K_2(t)$ and may influence the estimation of $D$. Nevertheless, if $T$ is large enough, it turns out

$$
\lim_{T \to \infty} G_2(t, \tau) = G_1(t, \tau).
$$

Finally, the extended model $G_2$ can be interpreted as an
“out of equilibrium” state of the observed system, that
allows us to take into account perturbations such as noisy
data, non-uniform background or other spots diffusing in
the image.

**Bayesian estimation of model parameters**

In this section, we present a Bayesian approach to estimate
$D$ and $\sigma_{\text{PSF}}$.

**Principles of Approximate Bayesian Computation.** As the estimation of *D* and *σ*<sub>PSF</sub> involved in (9) given image data is highly non-linear, we resort to another estimation approach when the likelihood (or data fidelity term) is not tractable or too complex. The so-called Approximate Bayesian Computation (ABC) [4] was especially designed to address this issue and relies on stochastic simulation that generates samples and selects those follow the posteriori distribution. This approach allows us to compute both the maximum a posteriori estimator and the a posteriori expectation from the selected samples.

Formally, let us assume the data *z* are discrete and are
generated from a model with parameters *θ* whose prior is
denoted by *p*(*θ*). The posterior distribution of interest is
defined as

$$
p(\theta|z) = \frac{p(z|\theta)p(\theta)}{p(z)}
$$

where $p(z) = \int p(z|\theta)p(\theta)d\theta$. The success of the rejection ABC method [4] depends on the fact that the underlying stochastic process is easy to simulate for a given set of parameters and the likelihood $p(z|\theta)$ is complicated as in (9). The ABC procedure can be summarized as below:

1) Generate $\theta$ from $p(\cdot)$;
2) Simulate $z'$ from the model with parameter $\theta$;

3) Calculate a distance $\rho(z, z')$ between $z'$ and $z$;

4) Accept $\theta$ with probability $\frac{\pi_e(z, z')}{c}$ and return to 1.

Here $c$ is a constant chosen to guarantee that $\pi_e(z, z')$ defines a probability:

$$ \frac{\pi_e(z, z')}{c} = \begin{cases} 1 & \text{if } \rho(z, z') \le \epsilon, \\ 0 & \text{otherwise.} \end{cases} $$

The data $z$ and $z'$ are both real valued arrays of matching dimension and $\rho(z, z')$ evaluates the distance between the simulated and observed data. This algorithm gives the exact inference for the model which assumes a uniform measurement error in the region defined by the cut-off value $\epsilon$.

There are several advantages to this rejection method, among them the fact that they are usually easy to code, and they generate independent observations. This approach requires design of a suitable metric $\rho(\cdot)$ as well as a choice of $\epsilon$. As $\epsilon \to \infty$, accepted samples come from the prior, and $\epsilon \to 0$ accepted samples follow the target posterior distribution $p(\theta|z)$. The choice of $\epsilon$ reflects the balance between computability and accuracy, and for a given $\rho$ and $\epsilon$, accepted samples are independent and identically distributed from $p(\theta|\rho(z, z') \le \epsilon)$.

The next step is to calculate posterior expectation defined as $\hat{\theta}_{\text{MMSE}} = E(\theta|z) = \int p(\theta|D) \theta d\theta$, which is known as the minimum mean square error (MMSE) estimator. The simplest way to approach to compute the MMSE estimator is to draw samples $\{\theta_i\}_{i=1,\dots,N}$, from $p(\theta|z)$ using the previous algorithm and approximate $\hat{\theta}_{\text{MMSE}}$ using the sum $N^{-1} \sum \theta_i$. Instead of a uniform average, there is a direct extension to estimate a weighting average of samples $\theta_i$ by considering Epanechnikov kernels [15]. Meanwhile, the samples can be exploited to compute the posterior distribution for which the maximum mode equals the maximum a posteriori (MAP) estimator $\hat{\theta}_{\text{MAP}}$.

**ABC algorithm and implementation details.** In our context, the data $z$ is the autocorrelation computed from the experimental image sequence and $z'$ is the autocorrelation simulated from the model (6) and (9) with parameter $\theta = (D, \sigma_{\text{PSF}})^T$. We implemented the ABC rejection-based method with 0-1 cut-off and the distance $\rho(z, z') = \|z - z'\|_2^2$.

The prior distribution of both parameters $D$ and $\sigma_{\text{PSF}}$ is assumed to be uniform and depend on the application. Typically, $D \sim U[0.1, 2.0]$ and $\sigma_{\text{PSF}} \sim U[\sqrt{0.1}, \sqrt{2}]$. The number $N_{\text{samples}}$ is set to 100,000 samples. The rate of acceptation $\epsilon$ is set to 1%, that is the best 1000 samples are used to compute the posterior distribution, $\hat{\theta}_{\text{MAP}}$, and $\hat{\theta}_{\text{MMSE}}$.

The algorithm takes as input an autocorrelation signal $z$ computed from fluorescence intensities in the ROI as

follows:

$$ z(t, \tau) = \frac{1}{W \times H} \sum_{i=1}^{W} \sum_{j=1}^{H} \frac{(f_{i,j}(t + \tau) - \bar{f})(f_{i,j}(t) - \bar{f})}{\bar{f}_t^2} \quad (10) $$

where $(i, j)$ denotes a pixel in $\Omega$, $W$ and $H$ are the width and height of the ROI ($|\Omega| = W \times H$), $\bar{f} = |\Omega|^{-1} \sum_{i=1}^{W} \sum_{j=1}^{H} \sum_{t=0}^{T-1} f_{i,j}(t)$, and $\bar{f}_t = |\Omega|^{-1} \sum_{i=1}^{W} \sum_{j=1}^{H} f_{i,j}(t), t \in \{0, \dots, T-1\}$.

The algorithm is relatively fast as each sample $z'(t, \tau)$ from (6) or (9) takes about 0.01 second to generate. Moreover, all samples in a ROI are independent and can be generated in parallel. This means that the processing of multiple ROIs in an image sequence can be performed very efficiently.

## RESULTS

In this section, we evaluate the performance of BayesTICS on synthetic and real images depicting multiple diffusing spots at the plasma membrane observed in TIRF microscopy as illustrated in Fig. 1. We apply the BayesTICS algorithm to compute the MAP and MMSE estimators of diffusion coefficient by considering the two models (6) and (9), respectively. The ROI size depends on the spot size and here is set to 21 × 21 pixels windows from the analysis of real images. Meanwhile, the temporal series should be long enough to observe stationary fluctuations. In our experiments, $T > 100$ time points. As these choices may involve some arbitrariness, we assessed the sensitivity of $|\Omega|$ and $T$, as well as the position of the spot of interest in the ROI. In what follows, we show that changes in these parameters lead to essentially similar results and demonstrate the robustness of BayesTICS with respect to noise levels and cluttered background (i.e., presence of multiple spots in a ROI). We focus on the estimation of $D$ which is the parameter of interest.

### Evaluation on simulated data

First, in order to test the robustness of BayesTICS to noise, we simulated six image sequences depicting one diffusing spots in a ROI of 21×21 pixels, with a diffusion coefficient equal to 0.100 pixel/frame. The image intensities were then corrupted with Gaussian white noise with different standard deviations yielding image sequences with different signal-to-noise ratios. Figure 2 shows the 2D images at time $t = 0$ (left) extracted from the six image sequences and the autocorrelation plots (middle and right). We reported the results of BayesTICS, that is the posterior distributions, $\hat{D}_{\text{MAP}}$, and $\hat{D}_{\text{MMSE}}$ for the two autocorrelation models $G_1$ and $G_2$, respectively. The results suggest that the performance of estimators $\hat{D}_{\text{MAP}} \in [0.104, 0.126]$ pixel/frame and $\hat{D}_{\text{MMSE}} \in [0.130, 0.134]$ pixel/frame.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Microscopy image depicting diffusing pHluorin-tagged spots at the plasma membrane observed in TIRF microscopy. (courtesy of PICT facility, UMR144-CNRS Institut Curie)</figcaption></figure>

pixel/frame for $G_1$ and $\hat{D}_{MAP} \in [0.100, 0.101]$ and $\hat{D}_{MMSE} \in [0.136, 0.138]$ for $G_2$, are not hardly influenced by the amount of noise in images. We observe that the MMSE values are slightly higher than the true diffusion coefficient ($D = 0.100$ pixel/frame) whatever the selected model $G_1$ or $G_2$.

Furthermore, we evaluated the influence of the ROI size, the spot position in ROIs, and the amount of noise in ROIs. We analyzed six different ROIs shown in Fig. 3. The spot of interest (red circle) is not located at the center of the ROI. The size of the ROI vary from small regions to large square/rectangle windows which may contain background and additional spots. The signal-to-noise ratio changes from low values (upper-left image in Fig. 3) to high values (lower-right image in Fig. 3). In this experiment, we obtained $\hat{D}_{MAP} \in [0.190, 0.255]$ pixel/frame and $\hat{D}_{MMSE} \in [0.242, 0.282]$ pixel/frame for $G_1$ while $D_{true} = 0.250$ pixel/frame. The variance and inaccuracy of the MAP estimator is induced by the combined effects of the three sources of perturbation: noise, ROI size and spot positioning. The MMSE estimators are surprisingly very close for the six test sequences and ROIs, even in the case of ROIs with additional diffusing spots. In Supplementary Information/Supplementary Figs. 1-3, we reported additional results by changing the ROI size and the spot location, in the absence of noise.

## Experiments on real data: Transferrin Receptor molecules diffusion at different locations at the plasma membrane

We evaluated the potential of the BayeSTICS algorithm on real data TIRFM 2D image sequences depicting Transferrin Receptor proteins (TfR) tagged with pHluorin (pH-sensitive probe) in M10 cells, at the late steps of exocytosis, that is when the vesicles fuse with the plasma membrane. Each sequence has 600 images of size $256 \times 256$ (see Fig. 1), acquired at 10 frames/second.

Local diffusion events are first detected at the plasma membrane with an appropriate algorithm (e.g., [3]) (or manually selected by an user) and several temporal series of ROIs with the same size $|\Omega| = 21 \times 21$ are automatically extracted as illustrated in Fig. 4. The non-uniform fluorescent background in $\Omega$ is approximately estimated as the median of the intensity values, and further subtracted to reduce bias in the estimation of diffusion. The initial time point $t_0$ is automatically found in our experiments (TIRFM, pHluorin-TfR tagging) by detecting the frame in which the fluorescence intensity is maximum in the sequence composed of frames. The frames with index $t < t_0$ are then discarded. Finally, we assume that the images are corruped by white Gaussian noise. A variance stabilizing transform (Generalized Anscombe transform [7]) is potentially applied to produce a normally distributed noise as recommended in [2].

We selected five ROIs (S1, S11, S12, S3, S6) in the image shown in Fig. 1 which contain a diffusing spot and variable background. We applied BayeSTICS to compute the posterior distributions of $D$ and $\sigma_{PSF}$ and the MAP and MMSE estimators by considering the two models $G_1$ and $G_2$. In the present situation, the MAP and MMSE estimators ($\hat{D}_{MAP} \in [0.10, 0.32]$, $\hat{D}_{MMSE} \in [0.14, 0.33]$) computed with the $G_2$ model are more consistent with those published in the literature [3]. Additional similar results on several ROIs confirmed that $G_2$ is more appropriate to handle complex contents as observed in real images.

## CONCLUSION

In this paper, we have proposed an original method to analyse the vesicle fusion to plasma membrane at the end of the exocytosis process, observed in TIRF microscopy. We focus on Transferrin receptor (TfR) as biological model, which is a transmembrane protein inserted in cellular membranes. We proposed BayeSTICS which is a local autocorrelation approach for diffusion quantification, established in the Bayesian setting (Approximate Bayesian Computing). After demonstrating the efficiency and accuracy of the method on simulated sequences, we successfully applied it to real TIRFM images depicting TfR proteins.
The proposed method is fast and can locally analyze
diffusion unlike previous methods. In the future, we plan
to use BayesTICS in large-scale studies in 2D-TIRFM,
and to extend the method to 3D multi-angle TIRFM [6]
to decipher the late steps of the exocytosis plateform
Rab11/Rab11-FIP (Rab11 family of interacting proteins).

ACKNOWLEDGEMENTS

We thank J. Salamero, L. Leconte, and the PICT-Cell and Tissue
Imaging Facility (member of the National Infrastructure France
BioImaging / ANR-10-INBS-04) for providing real TIRM image
sequences.

[14] L. Malacrida, P.N. Hedde, S. Ranjit, F. Cardarelli, and E. Gratton. Visualization of barriers and obstacles to molecular diffusion in live cells by spatial pair-cross-correlation in two dimensions. *Biomed. Opt. Express*, 9(1):303–321, 2018.

[15] P. Marjoram, P. Molitor, V. Plagnol, and S. Tavaré. Markov chain Monte Carlo without likelihoods. *Proceedings of the National Academy of Sciences*, 100(26):15324–15328, 2003.

[16] J. Philibert. One and a half century of diffusion: Fick, einstein, before and beyond. *Diffusion Fundamentals*, 2(1):1–10, 2005.

[17] H. Qian, M.P. Sheetz, and E.L. Elson. Single particle tracking: Analysis of diffusion and flow in two-dimensional system. *Bio-physical Journal*, 60:910–921, 1991.

[18] D.M. Soumpasis. Theoretical analysis of fluorescence photo-bleaching recovery experiments. *Biophysical Journal*, 41(1):95–97, 1983.

REFERENCES

[1] P.F.F. Almeida and W.L.C. Vaz. Lateral diffusion in membranes. *Handbook of biological physics*, 1(C):305–357, 1995.

[2] A. Basset, J. Boulanger, S. Salamero, P. Bouthemy, and C. Kervrann. Adaptive spot detection with optimal scale selection in fluorescence microscopy images. *IEEE Transactions on Image Processing*, 24(11):4512–4527, November 2015.

[3] A. Basset, P. Bouthemy, J. Boulanger, F. Waharte, J. Salamero, and C. Kervrann. An extended model of vesicle fusion at the plasma membrane to estimate protein lateral diffusion from TIRF microscopy images. *BMC Bioinformatics*, 18:352, 2017.

[4] M. A. Beaumont, W. Zhang, and D. J. Balding. Approximate Bayesian computation in population genetics. *Genetics*, 162:2025–2035, 2002.

[5] T. Bickel. A note on confined diffusion. *Phys. A, Stat.Mech. App.*, 377(1):24–32, 2007.

[6] J. Boulanger, C. Gueudry, D. Münch, B. Cinquin, P. Paul-Gilloteaux, S. Bardin, C. Guérin, F. Senger, L. Blanchoin, and J.Salamero. Fast high-resolution 3d total internal reflection fluorescence microscopy by incidence angle scanning and azimuthal averaging. *Proc. Natl. Acad. Sci. USA*, 111(48):17164–17169, 2014.

[7] J. Boulanger, C. Kervrann, P. Bouthemy, P. Elbau, J. B. Sibarita, and J. Salamero. Patch-based nonlocal functional for denoising fluorescence microscopy image sequences. *IEEE Trans. Medical Imaging*, 29(2):442–454, Feb 2010.

[8] N. Chenouard, I. Smal, F. De Chaumont, M. Maska, I.F. Sbalzarini, Y. Gon, J. Cardinale, C. Carthel, S. Coraluppi, M. Winter, A.R. Cohen, W.J. Godinez, K. Rohr, Y. Kalaidzidis, L. Liang, J. Duncan, H. Shen, Y. Xu, K. Magnusson, J. Jalden, H.M. Blau, P. Paul-Gilloteaux, P. Roudot, C. Kervrann, F. Waharte, J.-Y. Tinevez, S.L. Shorte, J. Willemse, K. Celler, G.P. Van Wezel, H.-W. Dan, Y.-S. Tsai, C. Ortiz De Solorzano, J.-C. Olivo-Marin, and E. Meijering. Objective comparison of particle tracking methods. *Nature Methods*, 11(3):281–289, 2014.

[9] M.A. Digman and E. Gratton. Lessons in fluctuation correlation spectroscopy. *Ann. Rev. Physical Chemistry*, 62:645–668, 2011.

[10] B. Hebert, S. Costantino, and P.W. Wiseman. Spatiotemporal image correlation spectroscopy (stics) theory, verification, and application to protein velocity mapping in living cho cells. *Bio-physical J.*, 88(5):3601–3614, 2005.

[11] L. Ji and G. Danuser. Tracking quasi-stationary flow of weak fluorescent signals by adaptive multi-frame correlation. *J. Microscopy*, 220(3):150–167, 2005.

[12] S.A. Kim, H. Sanabria, E. Digman, M.A. ad Gratton, P. Schwille, W.R. Zipfel, and M.N. Waxham. Quantifying translational mobility in neurons: Comparison between current optical techniques. *The Journal of Neuroscience*, 30(49):16409–16416, 2010.

[13] D.L. Kolin and P.W. Wiseman. Advances in image correlation spectroscopy: measuring number densities, aggregation states, and dynamics of fluorescently labeled macromolecules in cells. *Cell Biochem. Biophys.*, 49(3):141–164, 2007.
Theoretical diffusion coefficient $D_{th} = 0.10$

<figure><img src="image_3.png" /><figcaption>Fig. 2. Robustness of the BayesTICS method to noise level with a fixed window size. Six different ROIs with increasing amount of noise (from top to bottom) were simulated. **a**, The ROIs at time $t = 0$ are extracted from six simulated sequences composed of $T = 300$ frames (256 × 256 pixels images) and depicting 2D diffusing spots with a theoretical diffusion coefficient $D_{true} = 0.100$ pixel/frame. The size of ROIs is set to 21 × 21 pixels and the spots are located at the center of each ROI. **b, c**, The autocorrelation versus time lag plot are shown for $G_1$ (**b**) and $G_2$ (**c**) models. Each plot shows the observed autocorrelation (black curve), and two autocorrelation samples generated from the $G_1$ and $G_2$ models with the $\hat{D}_{MAP}$ and $\hat{D}_{MMSE}$ parameters (red curves). The estimated posterior distribution are displayed for $D$ and $\sigma_{PSF}$ in each case.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Fig. 3. Robustness of the BayesTICS method to noise level with variable spot position and window size. Six different cases with varying noise level, spot position and window size were tested. They were extracted from six simulated sequences of 2D diffusing spots, with 256 × 256 pixels window size, total length of 300 frames, a theoretical diffusion coefficient of 0.250 pixels/frame. For each case, the z-projection of the maximum intensity of the simulated stack and the autocorrelation versus time lag plot are shown. The spot of interest is market in purple, and can be anywhere in the ROI. Other spots can be in the ROI, diffusing or not. The size of the ROI has the following values : 39 × 39, 41 × 31, 26 × 25, 21 × 20, 48 × 43, 31 × 35 pixels. Each plot shows the computed autocorrelation from the data (dark gray), generated autocorrelations for the BayesTICS method (light gray), and the minimum (green) and mean (red) estimated autocorrelations. The two estimates $\hat{D}_{MAP}$ and $\hat{D}_{MMSE}$ for the diffusion coefficient are given on the corresponding plots.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Fig. 4. Evaluation of BayesTICS on real TIRF image sequences depicting Transferin Receptor proteins (TfR) tagged with pHluorin (pH-sensitive probe) in M10 cells. a, Five 21*times*21 pixels ROIs were selected. b, The observed and simulated autocorrelations versus time lag plot are shown for *G*₁ and *G*₂ models. c, Each plot shows the observed autocorrelation (black curve), and two autocorrelation samples generated from the *G*₁ and *G*₂ models with the *D̂*MAP and *D̂*MMMSE parameters (red curves). The estimated posterior distribution are displayed for *D* (left) and *σ*PSF (right) in each case.</figcaption></figure>
# SUPPLEMENTARY INFORMATION

## A. Temporal Image Correlation Spectroscopy (TICS)

The most popular methods for diffusion analysis are based on correlation measurements under the hypothesis of temporal stationarity of fluorescence signals [10], [11]. The so-called Spatio-Temporal Image Correlation Spectroscopy (STICS), derived from fluorescence correlation techniques, is widely used in fluorescence imaging to recover physical parameters such as directional flow or diffusion parameters of moving molecules. This method does not require any particle or object tracking and integrates the variations of fluorescence over space and/or time via correlation measures to access to information at the molecular level, such as diffusion coefficients or dominant flow speed and direction [10]. The generalized spatial and temporal correlation expression is defined as

$$r(\mathbf{w}, \tau) = \frac{1}{N - \tau} \sum_{t=1}^{N-\tau} \frac{\langle \delta I(\mathbf{x}, t) \delta I(\mathbf{x} + \mathbf{w}, t + \tau) \rangle}{\langle I(\mathbf{x}, t) \rangle \langle I(\mathbf{x}, t + \tau) \rangle} \quad (11)$$

where $I: \Omega \times [1, N] \to \mathbb{R}$ is an image sequence of $N$ frames with $\Omega$ the image domain, $\mathbf{w} \in \mathbb{R}^2$ are the spatial lags, $\tau \in [1, N]$ is the temporal lag, $\langle \cdot \rangle$ is the spatial average over a patch and $\delta I(x, y) = I(x, y, t) - \langle I(x, y, t) \rangle$. We point out that $r(\mathbf{w}, \tau)$ is not a normalized correlation criterion but enables to recover the biophysical parameters associated to density, motion of molecules, and diffusion coefficient [10].

For transport estimation, the goal is to estimate the translation vector corresponding to the correlation peak maximum. The static or immobile molecule population is usually filtered by local averaging and $r(\mathbf{w}, \tau)$ is computed by Fast Fourier Transform. A 2D Gaussian function is considered to estimate accurately the correlation peak over time [10] using a Levendberg-Marquardt optimization scheme. In the experiments, the analysis is performed on image blocks and the size of the blocks determines the scale of moving objects retrieved (see Fig. ??). The spatial lag between blocks is chosen to achieve an acceptable trade-off between spatial accuracy and computational time.

In a diffusive motion scenario, the following heat diffusion equation is satisfied: $I_t(\mathbf{x}) = D \Delta I(\mathbf{x})$ where $I_t$ is the temporal derivative of $I$, $D$ is the isotropic diffusion coefficient and $\Delta$ denotes the Laplacian operator. In order to estimate the scalar value $D$, the diffusion decay $\tau_d$ is derived from the following equation obtained by combination of (11) and the heat equation (see [13] for details) :

$$r(0, 0, \tau) = r(0, 0, 0)(1 + \tau/\tau_d)^{-1} + r_{\infty}(\tau) \quad (12)$$

where $r_{\infty}(\tau)$ is the long-time offset and $D = \overline{w_0}/4\tau_d$ ($\overline{w_0}$ denotes the temporal average of $\omega_0$ (laser beam size)).

## B. Influence of the ROI size and spot position in the ROI

### Influence of the ROI size

We simulated a sequence of 256. × 256 pixels noise-free images depicting multiple diffusing spots, with a diffusion coefficient equal to 0.100 pixels/frame (see Supplementary Fig. 1A). The spots appear at random locations and at different times. Four typical images are also shown in Supplementary Fig. 1B at time $t = 0, 16, 20$ and $26$. In our experiment, only a few subsequences corresponding to temporal ROIs composed of 150 frames are used.

We computed the MAP and MMSE estimators of $D$ by increasing the ROI sizes as illustrated in Supplementary Fig. 1C. The spot of interest is assumed to be located at the center of the ROI. The time $t = t_0$ corresponds to the frame where this spot appear at the cell surface. The two secondary diffusing spots appear at $t_0 + 15$ and $t) + 20$, respectively. The largest ROI covers the spot of interest and the two other spots diffusing in the neighborhood (upper left image of Fig. 1C). In our experiments, the smallest ROI contains only the spot of interest (lower right of Fig. 1C).

We clearly notice that the MAP estimator is not influenced by the ROI size: $\hat{D}_{MAP} \in [0.103, 0.107]$ pixel/frame while the ground truth is $D_{true} = 0.100$ pixel/frame. For ROIs larger than 21 × 21 pixels, the values for the MMSE estimator are in the range [0.121, 0.131] pixels/frame. Although less accurate than the MAP estimator, the MMSE estimator is more robust to variable ROI sizes. When the ROI size is smaller than 21 × 21 pixels, the MMSE estimator is $\hat{D}_{MMSE} = 0.147$ pixel/frame, and is negatively affected when the window size is too small. This result is not a
surprise since the MMSE estimator takes into account the spatio-temporal mean of the sequence; if the window size is small, there is not enough information to compute a consistent averaging. The plots show both autocorrelations corresponding to the MAP and MMSE estimators.

For the MAP estimator, the accuracy of the estimation is directly influenced by the fitting to the experimental curve, whereas, in the case of the MMSE estimator, the accuracy is not directly related to the fitting of the curves. The MAP autocorrelation fits better the experimental autocorrelation, as it is computed in order to minimize the distance between the two curves. The MMSE autocorrelation fits less the experimental autocorrelation, as it tends to balance between the curve to fit and the global information of the sequence. This should allow the MMSE estimator to be more stable for the cases that are not ideal: noisy data, other spots diffusing near the spot of interest, non centred ROI etc. As expected, in almost ideal conditions, the best estimation for the MAP estimator corresponds to the small window size case, with only the important information for the fitting, which is, the spot of interest. In contrast, for the MMSE estimator, the best estimations correspond to the large window sizes with more information in the ROI. For the parameter tested here, the autocorrelation has different ranges of values. The larger the ROI is, the higher the maximum value of the autocorrelation is, as well as the difference between the maximum and minimum values of the autocorrelation (see Supplementary Fig. 1C). The shape of the autocorrelation is the same in the six tested cases, and it corresponds to pure diffusion.

## Influence of the spot position in the ROI

The spot position may potentially influence the quality of the estimation. From the previous simulated sequence, a spot of interest at six different locations in the ROI was selected (see Supplementary Fig. 2). The spot is randomly located in the ROI in most cases. The ROI can also contain the full spot or only a part of the spot of interest. It can be also located near another diffusing spot. In waht follows, the size of the ROI is fixed in that experiment.

We applied the BayesTICS algorithm and compute the MAP and MMSE estimators. The results show that the MAP estimator is not influenced by the spot position : $\hat{D}_{MAP} \in [0.101, 0.108]$ pixel/frame. The MMSE estimator behaves in the same way, as long as the spot is entirely included in the ROI: $\hat{D}_{MMSE} \in [0.130, 0.137]$ pixel/frame. As before, the MMSE estimator is less accurate, but the estimation is robust the spot position. Nevertheless, we can notice that, when the spot is partially observed, the MMSE estimator performs differently. In that situation, the diffusion coefficient is overestimated: $\hat{D}_{MMSE} = 0.139$ pixel/frame. The plots in Supplementary Fig. 2 show both autocorrelations corresponding to the MAP and MMSE estimators. Note that, as the size of the ROI is fixed, the maximum value and range for the autocorrelation are the same for the six tested cases. The best estimator is the MAP estimator corresponds to the best fitting case, which is consistent with the observation in the previous scenario. The best fitting is attained when the spot of interest is isolated in the ROI. Interestingly, there is no preference for a particular spot position. This latter remark also applies to the MMSE estimator; the presence of additional spots in the ROI is favored by the MMSE estimator.

## Combined influence of the ROI size and spot position

Finally, it is important to jointly evaluate the sensitivity of the MAP and MMSE estimators for varying ROI sizes and spot locations in the ROIs. To assess this question, we selected four cases from the previously simulated sequence and computed the MAP and MMSE estimators with the BayesTICS method. Supplementary Fig. 3 shows the four tested cases, along with their respective autocorrelation plots and diffusion coefficient estimates. Various ROI sizes, and shapes were considered. The spot of interest was selected in different locations, and several other spots, appearing after the spot of interest, were diffusing in the ROI. We found $DMAP \in [0.101, 0.106]$ pixel/frame $\hat{D}_{MMSE} \in [0.1220.130]$ pixel/frame. The two estimators lie in the range of values observed in the case of each parameter being analyzed apart. This demonstrates that estimation is robust to the cumulated effect of the ROI size and spot position, for both estimators, as long as the two parameters are not in sensitive zone seen previously. Note that the MAP estimator gives more accurate results, despite several spots diffusing in the ROI, provided the spot of interest is not overlapped by the neighboring spots.
<figure><img src="image_6.png" /><figcaption>**Supplementary Fig. 1. Robustness of BayesTICS to window size (with centered spot).** A, Simulated noise-free image sequence (256 × 256 pixels, 300 frames) of 2D diffusing spots with a theoretical diffusion coefficient of $D_{true} = 0.100$ pixel/frame. The displayed image is the maximum intensity projection of the 3D stack. The ROIs used for estimating the diffusion coefficient are shown in red, and the spot of interest in purple. B, Four images of the sequence are shown at time $t = 0, 16, 20$ and $26$. The ROI and the spot of interest are also highlighted red (square) and purple (circle), respectively. C, Six ROI sizes were tested: $73 \times 57, 54 \times 43, 31 \times 28, 21 \times 19, 12 \times 11$, and $8 \times 7$ pixels windows. The ROI is centered on the spot of interest. Each plot shows the computed autocorrelation from the data (dark gray), generated autocorrelations for the BayesTICS method (light gray), and the minimum (green) and mean (red) estimated autocorrelations. The two estimates $\hat{D}_{MAP}$ and $\hat{D}_{MMSE}$ for the diffusion coefficient are given for each window size on their respective plots.</figcaption></figure>
Theoretical diffusion coefficient $D_{th} = 0.100$

<figure><img src="image_7.png" /><figcaption>**Supplementary Fig. 2. Robustness of BayesTICS to spot position (with fixed window size).** Six different cases were tested. Each ROI was extracted from a simulated image sequence (256 × 256 pixels, 300 frames) of 2D diffusing spots with a theoretical diffusion coefficient of $D_{true} = 0.100$ pixel/frame. For each ROI, we display the maximum intensity projection of the 3D stack and the plots of the autocorrelation versus time lag. The diffusing spot of interest is labeled in purple (circle). The spot can be located anywhere in the ROI. The ROI can contain supplementary spots diffusing potentially at the same time or not. The size of the ROI is fixed to 37 × 37 pixels. Each plot shows the computed autocorrelation from the data (dark gray), generated autocorrelations for the BayesTICS method (light gray), and the minimum (green) and mean (red) estimated autocorrelations. The two estimates $\hat{D}_{MAP}$ and $\hat{D}_{MMSE}$ for the diffusion coefficient are given on the corresponding plots.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Supplementary Fig. 3. Robustness of BayesTICS to spot position with variable window size.** Four different cases were tested. Each ROI was extracted from a simulated image sequence (256 × 256 pixels, 300 frames) of 2D diffusing spots with a theoretical diffusion coefficient of $D_{\text{true}} = 0.100$ pixel/frame, or each ROI, we display the maximum intensity projection of the 3D stack and the plots of the autocorrelation versus time lag. The ROIs are shown in red (rectangle) and the spot of interest in purple (circle). The window size varies, from up to bottom and to the left, as following: 81 × 72, 49 × 51, 63 × 39, and 44 × 31 pixels. The spot of interest can be located anywhere in the ROI. The ROI can contain supplementary spots diffusing potentially at the same time or not. Each plot shows the computed autocorrelation from the data (dark gray), generated autocorrelations for the BayesTICS method (light gray), and the minimum (green) and mean (red) estimated autocorrelations. The two estimates $\hat{D}_{\text{MAP}}$ and $\hat{D}_{\text{MMSE}}$ for the diffusion coefficient are given on the corresponding plots.</figcaption></figure>