A nu-space for image correlation spectroscopy: characterization and application to measure protein transport in live cells

We present a new generalized theoretical framework for ICS methods [11] that guides simple derivations of the new modes of analysis. We assume as input an image time series collected via fluorescence microscopy. The imaged region contains N fluorescent point emitters following their own transport dynamics and photophysics. We assume that the time to acquire an image is negligible, so that all spatial coordinates are sampled simultaneously. We also neglect molecules with fast dynamics which will not correlate on the time scale of imaging, as they will only produce noise in our correlation function. We can represent the acquired image pixel intensities as the following sum over all fluorescent molecules:

$$\begin{equation} I(\vec{r},t)= \sum_{i=1}^N B_i(t) B_0 \delta\big(\vec{r}-\vec{r}_i(t)\big) \otimes_{\vec{r}} I_{\mathrm{p}}(\vec{r}), \end{equation} \tag{ 1 }$$

where B₀ represents the brightness of the probe (including all constant terms related to the optical microscopy), B_i(t) its photophysics (0 or 1 depending on the bleaching and/or blinking of molecule i at time t), $\vec {r}_i(t)$ the position of particle i at time t and $I_{\mathrm {p}}(\vec {r})$ is the point spread function (PSF) of our imaging system and the subscript on the convolution specifies which dimension we are convolving. Note the use of vectors so that the formalism developed may be used for two-dimensional (2D) or three-dimensional imaging. Since $\vec {r}_i(t)$ is a random variable, $I(\vec {r},t)$ will also be a random variable. We use the correlation function

$$\begin{equation} R(\vec{\xi},\tau)= \frac{1}{L_{\vec{r},t}} I(\vec{r},t) \otimes_{\vec{r},t} I(-\vec{r},-t), \end{equation} \tag{ 2 }$$

where $L_{\vec {r},t}$ total number of pixels/voxels (e.g. L_x × L_y × L_t) and $\vec {\xi }$ the spatial lag variable. Since the PSF is only a function of space, it is useful to consider separately spatial and time variables. By inserting equation (1) into (2) and taking the Fourier transform in space, we obtain

$$\begin{eqnarray} &&\fl R(\vec{\xi},\tau) = \frac{B_0^2}{L_{\vec{r}}} \sum_i \mathcal{F}_{\vec{r}}^{-1} \left\langle \tilde{I}_{\mathrm{p}} \mathrm{e}^{2\pi\mathrm{ i} \vec{k}\centerdot\vec{r}_i(t)} \tilde{I}_{\mathrm{p}}^\star \mathrm{e}^{-2\pi \mathrm{i} \vec{k}\centerdot\vec{r}_i(t+\tau)} B_i(t)B_i(t+\tau) \right\rangle_t \nonumber \\ &&+ \frac{B_0^2}{L_{\vec{r}}} \sum_{i\neq j} \mathcal{F}_{\vec{r}}^{-1} \left\langle \tilde{I}_{\mathrm{p}} \mathrm{e}^{2\pi \mathrm{i} \vec{k}\centerdot\vec{r}_i(t)} \tilde{I}_{\mathrm{p}}^\star \mathrm{e}^{-2\pi \mathrm{i} \vec{k}\centerdot\vec{r}_j(t+\tau)} B_i(t)B_j(t+\tau) \right\rangle_t , \end{eqnarray} \tag{ 3 }$$

where $\tilde {I}_{\mathrm {p}}$ is the Fourier transform of the PSF, $\mathcal {F}^{-1}$ the inverse Fourier transform operator and we have separated the self and the cross terms. We will now make two simplifying assumptions. We will suppose that the particles are independent and the photophysics is only time dependent and hence is independent of the position. We can therefore rewrite the correlation function as a product of averages:

$$\begin{equation} R(\vec{\xi},\tau) = \frac{B_0^2 N}{L_{\vec{r}}} \mathcal{F}_{\vec{r}}^{-1} \left\langle |\tilde{I}_{\mathrm{p}}|^2 \mathrm{e}^{2\pi\mathrm{ i} \vec{k}\centerdot(\vec{r}_i(t)-\vec{r}_i(t+\tau))} \right\rangle_t \Theta(t,\tau) +\left\langle I \right\rangle ^2, \end{equation} \tag{ 4 }$$

where we defined the photophysics correlation function $\Theta (t,\tau )= \left \langle B(t)B(t+\tau ) \right \rangle$, assumed ergodicity to eliminate the sum (time average over a single particle equals ensemble average) and where the second term simplified because N is large and the particles are independent. Note that we kept the index i because of ergodicity (it can be any particle). If we now define the autocorrelation of the PSF, $C_{\mathrm {PP}} (\vec {r}) = \mathcal {F}_{\vec {r}}^{-1} (|\tilde {I}_{\mathrm {p}}|^2)$, and $n = N/L_{\vec {r}}$ as the density of particles in the region of interest, we can write the following:

$$\begin{equation} R(\vec{\xi},\tau) = n B_0^2 \Theta(t,\tau) C_{PP} (\vec{r}) \otimes_{\vec{r}} \left\langle \delta \left( \vec{r} - \left( \vec{r}_i\left(t\right) - \vec{r}_i\left(t+\tau\right) \right) \right)\right\rangle +\left\langle I \right\rangle ^2 \end{equation} \tag{ 5 }$$

We set $\vec {r}_i(0)=0$ to allow easier interpretation of the results but there is no loss of generality. Switching to ensemble averaging:

$$\begin{eqnarray} &&\left\langle \delta \left( \vec{r} + \vec{r}_i\left(\tau\right) \right)\right\rangle = \int \delta \left( \vec{r} + \vec{r}_i\left(\tau\right) \right) P_{\vec{r}_\tau|\vec{r}_0} \left(\vec{r}_i\left(\tau\right) |\vec{r}_i\left(0\right)=0 \right) \nonumber\\ &&P_{\vec{r}_0} \left( \vec{r}_i\left(0\right)= 0 \right) \mathrm{d}\vec{r}_i(0) \mathrm{d}\vec{r}_i(\tau)= P_{\vec{r}_\tau|\vec{r}_0} \left(-\vec{r} ~| \vec{r}_i\left(0\right)=0 \right), \end{eqnarray} \tag{ 6 }$$

where we used the fact that the integral over all possible states equals 1 and the properties of the Dirac delta function. In other words, this function represents the probability of finding a particle at some point in space if you know its original position. We can therefore write what we will call the STICS master equation:

$$\begin{equation} R(\vec{\xi},\tau) = n B_0^2 \Theta(t,\tau) C_{PP} (\vec{r}) \otimes_{\vec{r}} P_{\vec{r}_\tau|\vec{r}_0} \left(-\vec{r} ~| \vec{r}_i\left(0\right)=0 \right) +\left\langle I \right\rangle ^2. \end{equation} \tag{ 7 }$$

This equation explains the behaviour of the correlation function for any type of underlying molecular dynamics. We will now derive the equations for simple dynamics, from which it can easily be generalized.

2.1.1. Image correlation spectroscopy (ICS)

We start by determining the amplitude of the correlation peak. Since convolving the Dirac delta in this case merely shifts the function at the origin:

$$\begin{equation} R(\vec{0},0) = n B_0^2 \Theta(t,0) C_{\mathrm{PP}}(\vec{0}) + \left\langle I \right\rangle ^2, \end{equation} \tag{ 8 }$$

with B₀ and n being the biophysical parameters of interest. Noting that $\left \langle I \right \rangle = N B_0 \Theta (t,0) C_{PP}(\vec {0})/L_{\vec {r}}$, we define a normalized correlation function:

$$\begin{equation} r(\vec{\xi},\tau) = \frac{R(\vec{\xi},\tau) }{ \left\langle I \right\rangle ^2} -1, \end{equation} \tag{ 9 }$$

so that

$$\begin{equation} r(\vec{0},0) = \frac{1}{n \Theta(t,0) C_{\mathrm{PP}}(\vec{0})}. \end{equation} \tag{ 10 }$$

$C_{\mathrm {PP}}(\vec {0})$ is simply a measure of the area/volume of the PSF (its variance, or spread). Moreover, Θ(t,0) represents the average fraction of emitting particles during the time window. This is why $1/r(\vec {0},0)$ is often referred to as the number of particles per beam area (BA) or focal volume. When applied to a single frame, this technique is the original variant called ICS [12, 13]. The PSF is often approximated by a Gaussian of the form

$$\begin{equation} I_{\mathrm{P}} (r)= \mathrm{e}^{-2r^2/w_r^2} \end{equation} \tag{ 11 }$$

so that its autocorrelation is simply a Gaussian two times wider:

$$\begin{equation} C_{\mathrm{PP}}(\vec{\xi})= {\pi}\!w_r^2\exp(-|\vec{\xi}|^2/w_r^2) \end{equation} \tag{ 12 }$$

and that $C_{\mathrm {PP}}(\vec {0})= {\pi}\!w_r^2$ and w_r is the e⁻² radius of the Gaussian PSF.

2.1.2. Flow

A uniform flow is a deterministic process. If $\vec {v}$ is the molecular flow velocity, then

$$\begin{equation} P_{\vec{r}_\tau|\vec{r}_t} \left(-\vec{r} ~| \vec{r}_i\left(t\right)=0 \right) = \delta \left( \vec{r} + \vec{v}\tau \right). \end{equation} \tag{ 13 }$$

Thus, equation (7) becomes

$$\begin{equation} R_{\mathrm{flow}}(\vec{\xi},\tau) = n B_0^2 \Theta(t,\tau) C_{\mathrm{PP}} (\vec{\xi}+\vec{v}\tau) +\left\langle I \right\rangle ^2 . \end{equation} \tag{ 14 }$$

That is, the correlation function is equal to the autocorrelation of the PSF shifted by $-\vec {v}\tau$ and multiplied by the photophysics function.

2.1.3. Diffusion

We shall now consider uniform diffusion in d dimensions. The probability equation (6) is equivalent to the Green's function of the heat/diffusion equation:

$$\begin{equation} P_{\vec{r}_\tau|\vec{r}_t} \left(-\vec{r} ~| \vec{r}_i\left(t\right)=0 \right) = \frac{1}{\left( 4 {\pi}\!D \tau \right)^{d/2}} ~ \mathrm{e}^{-|\vec{r}|^2/4D\tau}. \end{equation} \tag{ 15 }$$

Therefore, convolving with (15) has the effect of spreading the correlation function as τ increases. For example considering the Gaussian PSF with 2D diffusion, equation (7) becomes

$$\begin{equation} \fl R_{\mathrm{diff}}(\vec{\xi},\tau) = n B_0^2 C_{PP}(\vec{0}) \Theta(t,\tau) \frac{w_r^2}{(w_r^2+4D\tau)} \exp\left(\frac{-|\vec{\xi}|^2}{w_r^2+4D\tau}\right) +\left\langle I \right\rangle ^2 .\end{equation} \tag{ 16 }$$

2.1.4. Biased diffusion

One can think of a particle undergoing simultaneous diffusion and flow as a more realistic biological model of directed molecular transport. In that case, the probability distribution for biased diffusion is

$$\begin{equation} P_{\vec{r}_\tau|\vec{r}_t} \left(-\vec{r} ~| \vec{r}_i\left(t\right)=0 \right) = \frac{1}{\left( 4 {\pi}\!D \tau \right)^{d/2}} ~ \mathrm{e}^{-|\vec{r}+\vec{v}\tau |^2/4D\tau}, \end{equation} \tag{ 17 }$$

which is the same as equation (15) but with an average velocity $-\vec {v}$. Applying the same assumptions as in the last section:

$$\begin{equation} \fl R_{\mathrm{diff}+\mathrm{flow}}(\vec{\xi},\tau) = n B_0^2 C_{\mathrm{PP}}(\vec{0}) \Theta(t,\tau) \frac{w_r^2}{ (w_r^2+4D\tau)} \exp\left(\frac{-|\vec{\xi}+\vec{v}\tau |^2}{w_r^2+4D\tau}\right) +\left\langle I \right\rangle ^2 . \end{equation} \tag{ 18 }$$

The previous section was derived as the autocorrelation in a single channel. Let us briefly state the important differences when analysing the cross-correlation of distinct signals. Suppose we are imaging different kinds of fluorophores in two separate emission detection channels. We can rewrite equation (3), but summing only over interacting particles

$$\begin{eqnarray} &&\fl R_{ab}(\vec{\xi},\tau) = \frac{B_{0a} B_{0b}}{L_{\vec{r}}} \sum_i \mathcal{F}_{\vec{r}}^{-1} \left\langle \tilde{I}_{Pa} \mathrm{e}^{2\pi \mathrm{i} \vec{k}\centerdot\vec{r}_{ia}(t)} \tilde{I}_{Pb}^\star \mathrm{e}^{-2\pi \mathrm{i} \vec{k}\centerdot\vec{r}_{ib}(t+\tau)} B_{ia}(t)B_{ib}(t+\tau) \right\rangle_t +\left\langle I_a \right\rangle \left\langle I_b \right\rangle. \nonumber \\ \end{eqnarray} \tag{ 19 }$$

Since the photophysics is independent between channels, it does not contribute

$$\begin{eqnarray} &&\fl R_{ab}(\vec{\xi},\tau) = \frac{B_{0a} \left\langle B_{a}(t)\right\rangle B_{0b} \left\langle B_{b}(t)\right\rangle N_{ab}}{L_{\vec{r}}} C_{ab} (\vec{r}) \otimes_{\vec{r}} P_{\vec{r}_{\tau,b}|\vec{r}_{t,a}} \left(-\vec{r} ~| \vec{r}_a \left(t\right)=0 \right) + \left\langle I_a \right\rangle \left\langle I_b \right\rangle, \nonumber \\ \end{eqnarray} \tag{ 20 }$$

where we defined $C_{ab} (\vec {r})= \mathcal {F}^{-1} \{ \tilde {I}_{Pa} \tilde {I}_{Pb}^\star \}$ as the cross-correlation between the two PSFs. This function will not necessarily be centred at zero depending on the alignment of the two channels. Improper alignment will merely spatially shift the correlation function, thereby automatically correcting for the misalignment of the two channels. The conditional probability has to be interpreted as: where is particle b likely to be found at a later time, knowing the position of a particle a at some time. Finally, there is a need for a double normalization in order to extract the number of interacting molecules N_ab. Indeed, the normalized correlation peak will give

$$\begin{equation} r_{ab}(\vec{0},0) = \frac{ \left\langle B_{a}(t)\right\rangle \left\langle B_{b}(t)\right\rangle n_{ab}C_{ab}(0)}{n_a n_b \left\langle B_{a}^2(t)\right\rangle\left\langle B_{b}^2(t)\right\rangle C_{aa}(0) C_{bb}(0)}. \end{equation} \tag{ 21 }$$

Thus

$$\begin{equation} \frac{r_{ab}(\vec{0},0)}{r_{aa}(\vec{0},0)r_{bb}(\vec{0},0)} = \left\langle B_{a}(t)\right\rangle \left\langle B_{b}(t)\right\rangle n_{ab} C_{ab}(\vec{0}). \end{equation} \tag{ 22 }$$

Not all tagged proteins will have the same dynamics so they can be grouped in independent populations. For example, a significant fraction of the proteins may be immobilized within the cell, while another population diffuses and another is being actively transported (i.e. flowing). Assuming no interaction between populations, we can do the same derivation as in section 2.1 but summing over individual subgroups to obtain

$$\begin{equation} R(\vec{\xi},\tau)= \sum_k N_k R_k(\vec{\xi},\tau) + \left\langle I \right\rangle ^2, \end{equation} \tag{ 23 }$$

where we assumed for simplicity that they have the same brightness (different brightness would simply add a weighting factor). Defining the fraction of molecules in a population f_k = N_k/N:

$$\begin{equation} r(\vec{\xi},\tau)= \sum_k f_k r_k (\vec{\xi},\tau), \end{equation} \tag{ 24 }$$

where

$$\begin{equation} r_k(\vec{\xi},\tau) = \frac{\Theta(t,\tau)}{n C_{\mathrm{PP}}(0)^2\Theta(t,0)^2} C_{\mathrm{PP}} (\vec{r}) \otimes_{\vec{r}} P_{\vec{r}_{\tau,k}|\vec{r}_{t,k}} \left(-\vec{r} ~| \vec{r}_k\left(t\right)=0 \right). \end{equation} \tag{ 25 }$$

The total correlation function is therefore a sum of the individual functions weighted by their respective fraction. This is why a dominant population representing a vast majority of the molecules can overshadow the dynamics of the minority population(s). This effect has motivated the developments described in this paper. Note that it is possible to extract f_k by normalizing by the total peak amplitude

$$\begin{equation} \frac{r_k(0,0)}{r(0,0)} = f_k. \end{equation} \tag{ 26 }$$

Since populations with different dynamics will be segregated into different regions of the frequency space, we now consider ν-space ICS (nICS) as a way to separate the immobile population. nICS can be understood as the spatial correlation of the temporal power spectrum. Taking the time Fourier transform of the mean corrected equation (7) gives the nICS master equation:

$$\begin{equation} R(\vec{\xi},\nu) = n B_0^2 (C_{\mathrm{PP}} (\vec{r}) \otimes_{\vec{r}} \mathcal{F}_t\lbrace P_{\vec{r}_\tau|\vec{r}_t} \left(-\vec{r} ~| \vec{r}_i\left(t\right)=0 \right) \rbrace ) \otimes \tilde{\Theta}(\nu). \end{equation} \tag{ 27 }$$

The correlation function is now convolved by the reciprocal of the photophysics function. This suggests that nICS will be more sensitive to the photophysics. However, in many cases photobleaching can be 'detrended' by adding a random number whose mean is equal to the required number to flatten the image time series [14, 15]. Furthermore, on the slower imaging time scale sampling the molecular transport, fluorescent protein blinking will appear only as a lower effective quantum yield.

2.4.1. Flow

As the time Fourier transform of the conditional probability is not trivial, we only show the equation for a flow with a Gaussian PSF:

$$\begin{eqnarray} &&\fl R_{\mathrm{flow}}(\vec{\xi},\nu) = n B_0^2 C_{\mathrm{PP}}(\vec{0}) \frac{\sqrt{\pi}w_r}{|v|} \exp \left(\!- \left(\!\frac{\pi \nu}{v/w_r} \!\right)^2 \right) \exp \left(\!- \left(\!\frac{v_x \eta - v_y \xi}{vw_r} \!\right)^2 \!\right) \exp \left(\! -2 \pi \mathrm{i} \nu \frac{\vec{v} \cdot \vec{\xi}}{v^2} \!\right) .\nonumber \\ \end{eqnarray} \tag{ 28 }$$

By taking the limit of v → 0, we obtain a Dirac delta for the plane ν = 0, which indicates that the immobile population will be concentrated in a single plane. While working with finite data sets, it is important to use proper methodology for power spectrum estimation [16]. That is, the computed power spectrum will be a convolution of equation (27) and the time Fourier transform of the temporal window function.

2.4.2. Biased diffusion

For biased diffusion with a Gaussian PSF, the correlation function is

$$\begin{eqnarray} &&\fl R_{\mathrm{diff}+\mathrm{flow}}(\vec{\xi},\nu) = n B_0^2 C_{\mathrm{PP}}(\vec{0}) \left( \int_{-\infty}^{\infty}\mathrm{d}\tau \mathrm{e}^{-2\pi \mathrm{i} \nu \tau} \frac{w_r^2}{ (w_r^2+4D|\tau |)} \exp\left(\frac{-|\vec{\xi}+\vec{v}\tau |^2}{w_r^2+4D|\tau |}\right) \right) \otimes \tilde{\Theta}(\nu). \nonumber \\ \end{eqnarray} \tag{ 29 }$$

Note that this equation includes the case for pure diffusion by taking $\vec {v}=0$. In practice, data sets can be fitted numerically to the time Fourier transform of equation (18) with the appropriate imaging parameters. Since we are evaluating the power spectrum from the correlation function with finite discrete sampling, we have to select a window function with a real positive Fourier transform to avoid spectral leakage of the negative components [16].

Since different populations are concentrated in different parts of the frequency spectrum, another approach to separating them is to filter them out before computing the correlation function. This approach has the advantage that a dominant population like the immobile one will be removed before computation of the correlation function, therefore minimizing leakage of this component. However, it will affect and, to some extent, distort all other populations. We now propose the use of a first-order Butterworth [17] (infinite impulse response (IIR)) filter that does not depend on the length of the time series analysed, introduces few artefacts in the analysis and can be tuned to filter specific dynamics on different time scales.

For simplicity, we show the effects of the filter on the zero spatial lag correlation function for a single flow with a Gaussian PSF. The full derivation is shown in appendix A. The correlation function of the output of the filter, ϕ_F (τ), is

$$\begin{equation} \fl \phi_{F}(\tau) = \mathrm{e}^{-T^2} - \frac{1}{2} \mathrm{e}^{R (R-2\,T)} \sqrt{\pi } R (1-\mathrm{Erf}[R-T]+\mathrm{e}^{4 R T} (1-\mathrm{Erf}[R+T])) ,\end{equation} \tag{ 30 }$$

where we introduced two dimensionless variables, a ratio R = πw_rν_c/v and a normalized time T = tv/w_r. As w_r/v is the time scale of fluctuation, R is simply a ratio between the cutoff frequency of the filter and the fluctuation frequency. This function is plotted for different values of R and T in figure 2.

Zoom In Zoom Out Reset image size

The amplitude of the correlation function decays rapidly as R increases and the effects on the shape of the function are quite minimal (figure 2). The correlation function can still be approximated by a Gaussian. Moreover, the error on the fitted width is not very important in STICS as we are interested in finding the peak centroid at each time lag.

The same development still holds for a diffusing population, but the calculations are more complicated due to the shape of equation (15). For the sake of brevity, we will simply state the result for the amplitude of the filtered 2D diffusion correlation function, ϕ_F,d(0).

$$\begin{equation} \phi_{{F},d}(0) = 1 - S \mathrm{e}^{S} {E}_1(S), \end{equation} \tag{ 31 }$$

where we used a definition of the exponential integral, ${E}_1(x)=\int _{x}^{\infty } \mathrm {e}^{-t}/t~ \mathrm {d}t$ and defined another dimensionless ratio, S = 2πν_cw²_r/4D. A flowing population and a diffusing one sharing a similar time scale (4Dw²_r = v/w_r) will be filtered similarly (not shown).

A simple alternative to filtering is to use the fact that the only time asymmetric component in the correlation function is the flow. Effectively, a diffusing and/or immobile population behaves the same way backwards in time, as opposed to the flow that will move in the opposite direction. Therefore, one could define an asymmetric correlation function:

$$\begin{equation} R_{\mathrm{asym}}(\vec{\xi}, \tau) = R(\vec{\xi}, \tau) - R(\vec{\xi}, -\tau). \end{equation} \tag{ 32 }$$

This is equivalent to removing the correlation of the time reversed intensity. It is a simple way to isolate the flow without disturbing it. However, the correlation function is still computed with all the other populations, which can give a noisy asymmetric correlation if the other populations are dominant. Furthermore, photobleaching and other non-stationary effects will be enhanced and should be corrected by 'detrending'.

In order to test the techniques developed, computer simulations were run using the custom Matlab (The MathWorks, Natick, MA) routine previously described [18]. These simulations model key aspects of fluorescence microscopy image time series acquisition. A list of typical simulation parameters is found in table 1.

Unless otherwise noted, these parameters were used for all simulations. Those that were kept constant for all simulations presented in this paper are indicated in bold. Note that we chose the parameters to closely match what we usually have in real cellular imaging data (i.e. [5]). For each data point shown, 32 simulations were run and the results are indicated either as mean ± standard deviation (STD) or as boxplots. The focus was on accuracy rather than precision since we analysed small regions of interest (ROI) and typically many of these can be fit inside a single cell. Dashed lines were included in the graphs to indicate the acceptable region where the relative error is less than 20%, which is less than the typical broad biological distribution encountered with small finite populations of cells.

For the IIR filtering, a highpass Butterworth filter was converted numerically through bilinear transformation. The data were then filtered pixel by pixel. The correlation function was computed and we rejected long time lags where the maximum of the function was of similar amplitude to random fluctuations. A Gaussian was fit for each time lag, and the slope of a linear fit of their centroid positions was used to determine the velocity. To extract the flow density, we used a linear regression on the fitted volume (height times BA) to obtain the zero lag peak height (as the correlation at τ = 0 will contain other components). This amplitude was normalized by the total amplitude obtained with ICS, to obtain a flow fraction. For biased diffusion, we did a nonlinear fit of the peak amplitudes with equation (18).

For the asymmetric analysis, we subtracted the space–time correlation function of the time-reversed data from the normal one. This created a negative peak in the direction of the flow. We then fitted two Gaussians of equal but opposite positions and amplitudes for each time lag. We rejected the first lags where the two peaks were too close to be distinguished and used the same criterion as in IIR filtering to reject the long lags. We then proceeded similarly to the IIR filter to obtain the flow velocity and density. For the diffusion coefficient, we used the slope on the linear fit of the peak width according to equation (18).

For the nICS analysis, we use a Hamming window half the time length of the signal. We then averaged the nICS correlation function obtained for three half-overlapping time windows. A nonlinear regression of the whole correlation function with equation (28) was used to extract the velocity. This equation was convolved in frequencies with the Fourier transform of the autocorrelation of the window. The peak amplitude was normalized to the total ICS peak amplitude to extract the flow fraction. For biased diffusion analysis, we fitted the time-Fourier transform of (18) (windowed by the autocorrelation of the Hamming window) to extract the different fit parameters.

Simulations including photobleaching were 'detrended' as specified in the following references [14, 15]. Briefly, the mean intensity of the images in the time series was fit to a mono-exponential decay. Random numbers were added to each pixel so that the average intensity would stay constant over time.

CHO.B2 cells (from Rudi Juliano and co-workers [19]) were cultured in DMEM medium (Invitrogen). Cells were transfected with Lipofectamine and Lipofectamine 2000 (Invitrogen) according to the manufacturer's instructions and α6-GFP and mCherry-paxillin were expressed as previously described [5, 20]. αL-GFP and β2-GFP were made by in-frame replacement in αL-YFP and β2-YFP, kind gifts of T. Springer.

Image time series were acquired on an Olympus IX70 inverted microscope (1.45 NA oil 60 ×  PlanApo objective) with a Ludl modular automation controller (Ludl Electronics Products) controlled by Metamorph (Molecular Devices). The 488 nm laser line of an argon ion laser and the 543 nm laser line of helium–neon laser (Melles Griot) were used to excite GFP and mCherry, respectively, and a dual emission filter (z488/543) was used to separate emissions. Images were acquired via a charge-coupled device camera (Retiga Exi, Qimaging).

The laminin (LM) and ICAM ligands were adsorbed to coverslips pre-coated with 1 mg ml⁻¹ poly-l-lysine, which has been previously reported to improve adsorbtion [21, 22].

To verify the theory of the frequency filter (equations (30) and (31)), we simulated simple flows and diffusion while changing the characteristic dimensionless ratio. We changed R/S by varying one parameter and keeping the other constant (ν_c = 2.5 mHz, v = 2 μm min⁻¹ / D = 0.01 μm² s⁻¹ and w_r = 0.4μm) and computed the remaining fraction of molecules by ICS. The results, shown in figure 3, agree well with the theoretical function, except for a general higher bias. This bias is probably due to the fact that we are still fitting a Gaussian to the correlation function. Therefore, the two negative sidelobes (as seen in figure 2) create an apparent greater 'Gaussian amplitude'.

Zoom In Zoom Out Reset image size

4.2.1. Pure flow

We extracted the fraction of molecules flowing (f), the recovered velocity (v) and angle (θ) as a function of increasing numbers of immobile particles and noise (width factor of the analogue detection noise (WF)). As correlation techniques are generally not affected by white noise (except in the zero lag), we added correlated noise in the form of a fast diffusion on the time scale of analysis (D = 1 μm² s⁻¹). In this analysis, we identified the fraction of mobile particles and split the rest equally between immobile and diffusing particles. Unless otherwise noted, default simulation parameters were used (see table 1).

We simulated an immobile fraction of up to 90%, a WF of up to 20 (where the SNR is approximately one) and a mobile fraction down to 10%. In most of these difficult situations, the techniques did not show a significant decrease in accuracy or precision (figure B.1). This suggests that the new techniques are not greatly affected by these parameters. Note that the accuracy of all techniques is not very good for recovering the fraction of molecules flowing as it relies on ICS and this technique works poorly in these low sampling conditions (small ROIs [18]).

Zoom In Zoom Out Reset image size

4.2.2. Biased diffusion

One can think of biased diffusion (simultaneous diffusion and flow) as a more realistic model of directed molecular transport in a cell. Effectively, due to the stochastic nature at the molecular level in the low Reynolds number sub-cellular environment, biological transport is probably never perfect flow. We define here a transport bias coefficient, B = vw_r/4D, so that it is close to one when flow and diffusion are on the same time scale.

We compared the new techniques as a function of the transport bias coefficient (figure 4). For the fitted density, the IIR analysis is accurate over most of the range. For nICS, it becomes heavily biased upwards as the fit fails more frequently, converging to a null velocity and a very high density. The asymmetric analysis gets slowly biased towards zero. As for the velocities, all of the techniques start diverging at a ratio of around 0.6.

Zoom In Zoom Out Reset image size

We were also able to fit the diffusion coefficient of the flow, which provides more information about the biological process. For the IIR and asymmetric analysis, the fit is very approximate and only the order of magnitude estimate of the diffusion coefficient could be obtained for a B between 0.2 and 5. The two techniques have an opposite and similar bias. Effectively, we used different methods in the IIR and asymmetric analysis to extract the diffusion coefficient. nICS performs remarkably well as it is able to measure an accurate and precise diffusion coefficient down to a bias coefficient of 0.2, even if it fails to fit the velocity and density. At high bias coefficient, the relative error becomes large since the diffusion coefficient is very low.

We repeated all the simulations for a biased diffusion; setting D to 0.005 μm² s⁻¹ to obtain a B close to one (figure B.2). The techniques are generally less precise than for pure flow, but display acceptable accuracy over a wide range of parameters.

Zoom In Zoom Out Reset image size

It is important to note that the IIR filter shows a higher bias in velocity. For a biased diffusion, the correlation function will spread as well as move for increasing time lags. This creates an apparent spread in velocities in the direction of the flow. Therefore, the filter will remove more of the low 'apparent velocity' than the fast ones (as it is a highpass). This distorts the correlation function, pushing the centroid of the peak in the direction of the flow and creating a higher measured velocity. In this paper, we carefully chose the cutoff frequency so that the effect would be minimal, but over filtering can lead to velocity up to two to three times higher than the real one (not shown). This effect is not present in the other techniques and highlights the importance of carefully selecting a cutoff frequency for filtering and verifying with other approaches.

We also verified that we could 'detrend' the time series that were affected by photobleaching. The results, shown in figure B.3, demonstrate that the correction was successful at all of the levels tested. The precision of the techniques gets lower as the fraction of bleached molecules increases, however, the accuracy is still acceptable even after 90% of molecules have bleached during the time series, which would be an unacceptable level for most live-cell imaging experiments. It is not very surprising that the techniques are not sensitive to photobleaching, as the correction involves adding uncorrelated random numbers (i.e. noise), and white noise does not significantly affect correlation techniques at non-zero lags.

Zoom In Zoom Out Reset image size