Single trajectory characterization via machine learning

RF is an architecture based on decision trees. A decision tree is an efficient non-parametric method widely used for classification and regression problems [30]. The basic idea consists in producing recursive binary splits of the input space, so that the samples with the same label are grouped together. The criterion to produce the splits is based on a homogeneity measure (usually, the information entropy) of the target variable within each of the obtained groups. In regression problems, a commonly used criterion is to select the split that minimizes the mean squared error; this recursive process continues until some stopping rule is satisfied, e.g. a common one is to consider that a tree node can be split if it contains more than a given number of samples; therefore, the minimum number of samples required to split a tree node should be adjusted in order to control the size of the tree, thus preventing overfitting. Once a decision tree is obtained, the output for unseen samples is computed just passing them through the nodes of the tree, where a decision is made with respect to which direction to take. Finally, a terminal tree node is reached, where the output is obtained.

A RF is a tree-based ensemble method, which builds several decision tree models independently and then computes a final prediction by combining the outputs of the different individual trees [31]. In particular, the ensemble is produced with single trees built from samples drawn randomly with replacement (bootstrap) from the training set. An additional randomness is added when splitting a tree node because the split is chosen among a random subset of the input variables, selected in this case without replacement, instead of the greedy approach of considering all the input variables. Due to this randomization, the bias of the ensemble is slightly higher than that of a single tree, but the variance is decreased and the model is more robust to variations in the dataset.

RF is a very powerful, state-of-the-art technique for both regression and classification problems, usually outperforming not only single decision trees but also sophisticated models, as shown in a thorough comparison study [32].

The training dataset is built out of numerical simulations of trajectories from various kinds of theoretical models. As a natural choice, we included three of the best-known and used models that can give rise to anomalous diffusion: CTRW [17], FBM [18] and Lévy walks (LW) [33]. In addition, we included trajectories from the annealed transient time model (ATTM) [23]. In the ATTM, a diffuser performs a random walk but stochastically changes the diffusion coefficient at random times. Both the diffusion coefficient and the time at which the diffusivity changes are drawn from distributions with a power law behavior [23]. Its time-averaged MSD shows a linear scaling, but the model has a regime in which it displays weak ergodicity breaking. The ATTM has been shown to reproduce the features observed for the diffusion of a cell membrane receptor [12], one of the experimental datasets analyzed.

Our aim is to design a method that can be used to accurately characterize heterogeneous trajectories without having to calculate other parameters or using a priori knowledge. In order to be able to analyze data coming from any possible spatiotemporal scale, we designed a preprocessing procedure that properly rescale the data. We implemented the following procedure, chosen among other plausible normalization techniques as it gives rise to the best results in terms of accuracy:

• 1.  
  We use one of the models above to simulate the trajectory of a particle during t_max time steps. The result is a vector of positions, ${\boldsymbol{X}}=({{\boldsymbol{x}}}_{1},{{\boldsymbol{x}}}_{2},\,\ldots ,\,{{\boldsymbol{x}}}_{{t}_{\max }})$.
• 2.  
  This vector is transformed into a vector of distances traveled in an interval of time T_lag, i.e. ${\boldsymbol{W}}=({\rm{\Delta }}{{\boldsymbol{x}}}_{1},{\rm{\Delta }}{{\boldsymbol{x}}}_{2},\,\ldots ,\,{\rm{\Delta }}{{\boldsymbol{x}}}_{J-1})$,, where $J={t}_{\max }/{T}_{\mathrm{lag}}$. We define ${\rm{\Delta }}{{\boldsymbol{x}}}_{i}$ as

  $$\begin{eqnarray}&&{\rm{\Delta }}{{\boldsymbol{x}}}_{i}=\left|{{\boldsymbol{x}}}_{{{iT}}_{\mathrm{lag}}}-{{\boldsymbol{x}}}_{(i+1){T}_{\mathrm{lag}}}\right|.\end{eqnarray} \tag{ 2 }$$
• 3.  
  To normalize the data, we divide ${\boldsymbol{W}}$ by its standard deviation to get a new vector $\hat{{\boldsymbol{W}}}$.
• 4.  
  Then, we do a cumulative sum of $\hat{{\boldsymbol{W}}}$ to construct a normalized trajectory $\hat{{\boldsymbol{X}}}$.

Summarizing, the previous procedure generates a new trajectory which is constructed via the normalized displacements of the original trajectory. This makes that the magnitudes of the resulting trajectories are comparable, no matter what were their original values. While the RF could be trained using $\hat{{\boldsymbol{W}}}$, our results show that training with $\hat{{\boldsymbol{X}}}$ gives indeed much better results. The same preprocessing is applied to both the simulated and experimental trajectories used in sections 2 and 3.

In order to predict the diffusion model underlying a certain trajectory, we construct a RF whose input is the normalized trajectory $\hat{{\boldsymbol{X}}}$, and the output is a number between 0 and N − 1 corresponding to the different models, with N the total number of models used in the training. Figure 3(a) shows the accuracy of the RF. Each line corresponds to a training dataset made up of different models. In the absence of data preprocessing (point marked as 'Raw' in the x-axis), the RF shows large accuracy. The accuracy drops significantly as T_lag increases, likely as a consequence of the removal of microscopical properties of the model, such as short-time correlations, hence preventing the RF from learning important features of them. This might lead to the conclusion that the filtering introduced by the preprocessing steps only limits the time resolution. This is obviously true for simulated data, obtained at the same scale, for which preprocessing is unnecessary. However, when dealing with experimental data of unknown spatiotemporal scale, such a preprocessing is of fundamental importance to be able to apply the same architecture and training dataset, in spite of the little loss of performance.

Zoom In Zoom Out Reset image size

In addition, the accuracy heavily depends on similarities among the models to be discriminated. For example, the accuracy obtained with a RF trained only with trajectories reproducing conceptually different models such as FBM and CTRW (triangular markers in figure 3(a)) is higher than the one obtained when including in the training models with similar characteristics, such as CTRW and ATTM, independently of T_lag (red circles and yellow squares in figure 3(a)).

A first approximation toward the characterization of the anomalous exponent can be based on a regression problem, in which the output of the RF is the value of the anomalous exponent α. The nature of the regression algorithm makes that the output of the RF is the continuous value which better satisfies the constraints learn during training.

To characterize the performance of the method, we calculate the prediction error ε of a trajectory as the absolute distance between the predicted exponent and the ground truth value. The percentage of trajectories $\bar{N}(\varepsilon )$ with a given error ε is represented in the bar plots of figure 3(c) for three different cases and a subdiffusive dataset including trajectories obtained from FBM, CTRW and ATTM. The case (i) considers trajectories with ${t}_{\max }={10}^{3}$ without noise, while the cases (ii) and (iii) show results for shorter and noisy trajectories (see discussion below). For case (i) the calculated mean absolute error (MAE) of the prediction of the anomalous exponent gives a value of 0.11. Moreover, the histogram showed in figure 3(c)(i) shows that for ∼80% of the trajectories, the output exponent lies within 0.1 from the true value.

A remarkable feature of the method is the possibility to correctly characterize very short trajectories. In figures 3(b) and (d), we show the ability of the RF to characterize short trajectories. In figure 3(b), we plot the accuracy in model discrimination as a function of the length of the trajectories, t_max. In figure 3(d), a similar study is done, now tracking the MAE of the RF trained to predict α. Although we observe an expected decrease of performance for short trajectories, both plots show that the RF is able to characterize trajectories as short as only 10 points. Quantitatively, when comparing trajectories of 10 points with larger ones, of 1000 points, the model discrimination accuracy only decreases by a factor of 8.2%, while the MAE decreases by a factor of 18%. Panel (ii) in figure 3(e) shows the error distribution when predicting α for ${t}_{\max }=100$.

Importantly, one has also to take into account that the experimental trajectories have a limited localization precision, that results into Gaussian noise. Therefore, it is necessary to test the robustness to noise of the RF. To this end, we trained the RF with trajectories simulated as described before and then we try to predict the anomalous exponent of trajectories belonging to the same dataset, but whose positions ${\boldsymbol{X}}$ were perturbed with noise to obtain the dataset ${{\boldsymbol{X}}}^{(n)}$

$$\begin{eqnarray}&&{{\boldsymbol{x}}}_{i}^{(n)}={{\boldsymbol{x}}}_{i}+{\mu }_{i}(\mu ,{\sigma }_{n}),\end{eqnarray} \tag{ 3 }$$

where μ_i(μ, σ_n) is a random number retrieved from a Normal distribution with mean μ = 0 and variance σ_n. The results obtained for training with FBM, CTRW and ATTMs are presented in figure 3(d). The RF shows a great robustness against noise. For σ_n < 1, the MAE appears almost unaffected. When increasing σ_n, we see that the MAE increases, as expected, but even for large σ_n the MAE is still reasonable.

Following the same structure of the previous section, we start by discriminating the diffusion model that can be associated to datasets (i)–(iii). The results are reported in table 1, showing a high rate of correct classification for the dataset (i). For the experimental data in datasets (ii) and (iii), we do not dispose of ground truth values, thus we compare our results with those of previous analysis, performed with alternative methods. For the trajectories of Experiment 1, we found that the algorithm largely assign them to the FBM, in strong agreement with previously reported results based on the concept of variation [14]. The data of Experiment 2are mainly assigned to the ATTM model. This model was shown to reproduce features observed in these data, such as subdiffusion and weak ergodicity breaking [12]. Moreover, a little fraction of trajectories are classified as CTRW. As previously mentioned, CTRW and ATTM share similar features (such as time subordination), increasing the difficulty in discriminating between them. This appears to be the main source of error in the results.

To obtain further insights on the study of the diffusion, we used the RF to extract the anomalous exponents. For the first dataset (i), based on simulations, we generated trajectories having a broad range of subdiffusive trajectories, namely α ∈ [0.2, 1]. Then, we used the trained RF to predict the value of the anomalous exponent and evaluated the error as the absolute value of the difference between the actual and predicted α. The results are reported in the histogram of figure 4(i) and display a distribution similar to the one obtained for the training/testing data. Thus, we run the same procedure on the experimental data. For the two datasets, in figures 4(ii)–(iii) we report the values obtained for the anomalous exponent α. The histogram of the α obtained for the trajectories of Experiment 1 (dark blue) shows mainly subdiffusive values, peaked in the range 0.6– 0.8. This is in good agreement with the original paper [16], where α was estimated by means of two different approaches as 0.7 and 0.77. However, the method also classifies a percentage of the trajectories as having α = 1. Importantly, the performance of the method can be further improved by taking advantage of the results of the model discrimination discussed above and shown in table 1. In fact, when the latter classification indicates that most of the trajectories follow a specific diffusive model, one can train the algorithm with a dataset composed only of trajectories simulated with that model. This kind of training produces exponent values in the same range, but largely reduce the fraction of those associated to α = 1, as shown in figure 4(ii) (light blue).

Zoom In Zoom Out Reset image size

Last, in figure 4(iii) we plot the distribution of exponents obtained for the Experiment 2. The subdiffusive values show a large number of occurrences in the 0.8–0.9 range, compatible with the exponent 0.84 calculated in previous studies [12]. Noteworthy, due to the nonergodic nature of the data, in the original paper α could only be calculated from the ensemble-averaged MSD, whereas the RF is able to determine this exponent from single trajectories.

We summarize here the different parameters used to train the RF models whose results are presented in this work. Common to all are the number of ensembled trees (100), the train/test dataset ratio (0.8/0.2) and the size of training set (1, 2 · 10⁵ trajectories). For the details of each figure, see table A1.