Statistical analysis of particle trajectories in living cells

Recent advances in molecular biology and fluorescence microscopy imaging have made possible the inference of the dynamics of molecules in living cells. Such inference allows us to understand and determine the organization and function of the cell. The trajectories of particles (e.g., biomolecules) in living cells, computed with the help of object tracking methods, can be modeled with diffusion processes. Three types of diffusion are considered: (i) free diffusion, (ii) subdiffusion, and (iii) superdiffusion. The mean-square displacement (MSD) is generally used to discriminate the three types of particle dynamics. We propose here a nonparametric three-decision test as an alternative to the MSD method. The rejection of the null hypothesis, i.e., free diffusion, is accompanied by claims of the direction of the alternative (subdiffusion or superdiffusion). We study the asymptotic behavior of the test statistic under the null hypothesis and under parametric alternatives which are currently considered in the biophysics literature. In addition, we adapt the multiple-testing procedure of Benjamini and Hochberg to fit with the three-decision-test setting, in order to apply the test procedure to a collection of independent trajectories. The performance of our procedure is much better than the MSD method as confirmed by Monte Carlo experiments. The method is demonstrated on real data sets corresponding to protein dynamics observed in fluorescence microscopy.

Figure 1

Representative trajectories from (a) simulated data and (b) a Rab11a protein sequence in a single cell (courtesy of UMR No. 144, CNRS, Institut Curie, PICT IBiSA). For the simulated data in (a), trajectory (1) is Brownian, trajectory (2) is from Brownian motion with drift, trajectory (3) is from fractional Brownian motion (with the parameter $\mathfrak{h}>1/2$), trajectory (4) is from an Ornstein-Uhlenbeck process, and trajectory (5) is from fractional Brownian motion ($\mathfrak{h}<1/2$). The parameters of the processes are given in Table 4.

Figure 2

Classification rule for motion modes from MSD. The solid line is the theoretical MSD of the standard Brownian motion. The dash-dotted lines are the bounds defined by [19]: $t\to t^{\beta }$, $\beta =0.9$, and 1.1. If $\widehat{\beta }$, the estimation of $\beta$, is such that $0.9<\widehat{\beta }<1.1$ it is classified as Brownian motion. The dashed lines are the pointwise high-probability interval of 95% associated with the empirical MSD curve for a standard Brownian motion trajectory of length $n=30$. The bounds of the interval are the 2.5% and 97.5% empirical quantiles of (2) and are computed by Monte Carlo simulation from $10001$ Brownian trajectories of size $n=30$.

Figure 3

Monte Carlo estimate of the power of the test at level $\alpha =0.05$ according to the trajectory length $n$ and the parameter associated with the following parametric alternatives: (a) Brownian motion with drift [with the parameter $v=(v_1,v_2)$ such that $v_1=v_2]$, (b) the Ornstein-Uhlenbeck process (parameter $\lambda$), and (c) fractional Brownian motion (parameter $\mathfrak{h}$). We use $10001$ Monte Carlo replications to compute each point of the power curves.

Figure 4

Box plots of the estimated exponent $\widehat{\beta }$. From left to right, we plot the box plot of $\widehat{\beta }$ corresponding to Brownian trajectories ($\mathfrak{h}=1/2$), subdiffusive fractional Brownian trajectories ($\mathfrak{h}=0.13$), and superdiffusive fractional Brownian trajectories ($\mathfrak{h}=0.85$). The trajectory size is $n=30$. Here $\widehat{\beta }$ is obtained by linear regression from the log(MSD) curve. We use the first 25 points of the MSD curve to estimate $\beta$. The dashed lines correspond to the true $\beta$ for each situation. We recall that in the case of fractional Brownian motion of the parameter $\mathfrak{h}$ we have $\beta =2\mathfrak{h}$.

Figure 5

Monte Carlo estimate of the average power against the proportion of true null hypotheses $m_0/m$ in the collection of hypotheses tested: (a) $m=100$ hypotheses and (b) $m=200$ hypotheses.

Figure 6

Box plots of the $p$ value $p_{30}$ [Eq. (23)] under $H_1$ and $H_2$. We simulate a set of trajectories ${\mathcal{X}}_m$ with $m=100$ and $m_0=20$ according to the simulation scheme described in Sec. 5. We plot the box plot of the $p$ values $p_{30}^{(i:m)}$ corresponding to each true alternative hypothesis $H_1$ and $H_2$. The green (orange) line is the threshold $h=p^{\left(k^{*}\right)}$ obtained by the first step of Procedure 1 (first step of adaptive Procedure 1). The null hypothesis is rejected if the $p$ value is lower than $h$. The black line is the level $\alpha =5\%$.

Figure 7

Map of the classification of the trajectories of the Rab11a sequence with (a) the standard multiple-test Procedure 1, (b) the adaptive version of Procedure 1, (c) MSD, and (d) a single-test procedure. The color code is blue for Brownian motion, red for superdiffusion, green for subdiffusion, and cyan for immobile particles (for the MSD method only).

Figure 8

Descriptive statistics of the 166 trajectories of the RAb11a sequence. (a) Histogram of the trajectories length. (b) The MSD curves of the first 50 trajectories. We plot only the first 30 points of the curves. (c) Histogram of the exponent $\beta$ obtained from a linear regression of the log(MSD). The green and red vertical lines correspond to the thresholds $\beta =0.9$ and $\beta =1.1$, respectively, proposed by [19] for classifying the trajectories. (d) Histogram of the diffusion coefficient (in $\mu {\text{m}}^2{\text{s}}^{-1}$) estimated with the estimator (19) which corresponds to the first point of the MSD curve.

Figure 9

Histogram of the test statistic $T_n$ obtained from the Rab11a sequence. The green (red) vertical line represents the quantile of order 2.5% (97.5%) of the asymptotic distribution of $T_n$.

Figure 10

Box plots of the proportions of Brownian motion, subdiffusion, and superdiffusion computed from 12 Rab11a sequences obtained with the single-test procedure (blue), Procedure 1 (cyan), the adaptive Procedure 1 (violet), and the MSD method (orange). Here Br stands for Brownian motion, Sb for subdiffusion, and Sp for superdiffusion.

Figure 11

Monte Carlo estimate of the power of the test at level $\alpha =0.05$ for different percentages of missing points (0, 20%, and 40%) along trajectories of size $n=30$. We give the power of the test for different parametric alternative hypotheses: (a) Brownian motion with drift [with the parameter $v=(v_1,v_2)$ such that $v_1=v_2]$, (b) the Ornstein-Uhlenbeck process (parameter $\lambda$), and (c) fractional Brownian motion (parameter $\mathfrak{h}$). We use $10001$ Monte Carlo replications to compute each point of the power curves. We note that we limit the study to a proportion of missing points less than 50%. Beyond 50%, we can consider that the time resolution between two observations $\mathrm{\Delta }$ has changed. Consequently, it is no longer a problem of missing points and we can restrict our study to the case where there is less than 50% missing points.

Figure 12

Monte Carlo estimate of the power of the test at level $\alpha =0.05$ for trajectories of size $n=30$ according to the signal-to-noise ratio and the parameter associated with the following parametric alternatives: (a) Brownian motion with drift [with the parameter $v=(v_1,v_2)$ such that $v_1=v_2]$, (b) the Ornstein-Uhlenbeck process (parameter $\lambda$), and (c) fractional Brownian motion (parameter $\mathfrak{h}$). We use $10001$ Monte Carlo replications to compute each point of the power curves.

Figure 13

Monte Carlo estimate of the power of the test at level $\alpha =0.05$ for a CTRW with different parameters settings: (a) $\mathrm{\Delta }/\tau =1$ and estimated power curves for different values of $n$ and (b) $n=30$ and estimated power curves for different values of $\mathrm{\Delta }/\tau$. We use $10001$ Monte Carlo replications to compute each point of the power curves.