Local Analysis of Heterogeneous Intracellular Transport: Slow and Fast Moving Endosomes

Trajectories of endosomes inside living eukaryotic cells are highly heterogeneous in space and time and diffuse anomalously due to a combination of viscoelasticity, caging, aggregation and active transport. Some of the trajectories display switching between persistent and anti-persistent motion, while others jiggle around in one position for the whole measurement time. By splitting the ensemble of endosome trajectories into slow moving subdiffusive and fast moving superdiffusive endosomes, we analyzed them separately. The mean squared displacements and velocity auto-correlation functions confirm the effectiveness of the splitting methods. Applying the local analysis, we show that both ensembles are characterized by a spectrum of local anomalous exponents and local generalized diffusion coefficients. Slow and fast endosomes have exponential distributions of local anomalous exponents and power law distributions of generalized diffusion coefficients. This suggests that heterogeneous fractional Brownian motion is an appropriate model for both fast and slow moving endosomes. This article is part of a Special Issue entitled: “Recent Advances In Single-Particle Tracking: Experiment and Analysis” edited by Janusz Szwabiński and Aleksander Weron.


Introduction
Intracellular transport of organelles, such as endosomes, has been described by anomalous diffusion caused by different mechanisms [1,2]. Various models have been proposed to describe it, such as fractional Brownian motion (FBM), continuous time random walks and fractional Langevin equations [3]. However, which of these models is the best is a current topic of much debate.
To decipher which mechanism is at work and determine the appropriate mathematical model to describe it, a large ensemble of trajectories is necessary. Modern experimental techniques facilitate the tracking of large ensembles of intracellular objects for considerable amounts of time. Therefore, the extraction of meaningful statistical information from trajectories is becoming an important issue. The traditional statistical analysis of trajectories includes quantification of ensemble evolution in time and space using the ensemble-averaged mean squared displacements (EMSD), time-averaged MSD (TMSD), probability density functions of displacements and correlation functions. As the accessible measurement time in experiments increases with better live-cell microscopy techniques, the accurate analysis of single trajectories has become possible [4]. New methods of trajectory analysis were developed, such as local time-averaged MSD [5], first passage probability analysis [6][7][8] and time-averaged diffusion coefficients [9].
Recently, the intracellular transport of endosomes in eukaryotic cells was shown to be described by spatiotemporal heterogeneous fractional Brownian motion (hFBM) with nonconstant Hurst exponents [40]. By analyzing the local motion of endosomes, we found that it is characterized by power-law probability distributions of displacements and displacement increments, exponential probability distributions of local anomalous exponents and powerlaw probability distributions of local generalized diffusion coefficients. In this paper, we split the ensemble of endosomes into slow and fast moving vesicles, which is the main difference between this study and that of [40]. This splitting allows us to study sub-ensembles separately in addition to studying the ATP driven active transport of endosomes. In particular, there is the central question: What is the appropriate mathematical model to describe the subdiffusive transport of slow moving endosomes? By analyzing locally the slow and fast endosomal trajectories, we find that both are characterized by exponential distributions of anomalous exponents and power-law distributions of generalized diffusion coefficients. This suggests that hFBM is an appropriate model for both slow and fast endosomes.
Endosome trajectories are composed of segments of active and passive motion, and therefore they could be further decomposed into directed runs and random motion. We segmented endosomal trajectories in this way in [10]. In this study, we separated endosome trajectories into superdiffusive trajectories and subdiffusive trajectories for their whole duration. Subdiffusive trajectories do not contain segments of directed movement and cannot be segmented further into active and passive motion. In contrast, fast superdiffusive trajectories can be further segmented. We leave the segmentation of fast trajectories into directed runs and random motion for future work.

Experimental Trajectories
We studied a large ensemble of two dimensional experimental trajectories, r(t) = {x(t), y(t)}, of early endosomes in a stable MRC5 cell line expressing GFP-Rab5. Trajectories were obtained from tracking wide-field fluorescence microscopy videos (see [10] for experimental details). We studied 103,361 experimental trajectories of early endosomes, the same data acquired in [10]. Three live-cell microscopy videos of MRC5 cells stably expressing GFP-Rab5 could be found in the Supplementary Material (https: //zenodo.org/record/5106450#.YPBsEuhKhPY, accessed on 23 July 2021). An example of experimental trajectories is shown in Figure A1. The endosomes were tracked using an automated tracking software (AITracker, based on a convolutional neural network) [41]. Currently, it is not yet feasible to determine the diameter of endosomes in these experiments, because they are diffraction limited. Thus, it was possible to track the centers of endosomes with sub-pixel accuracy, but not the sizes of the smaller endosomes (less than 200 nm). The duration of all trajectories, T, has a good fit to a power law distribution, T −1.85 [40], which is a manifestation of the heterogeneity of the trajectories. Slow moving endosomes stay longer within the observation volume and therefore have longer trajectories than fast moving endosomes, leading to the emergence of the power-law probability distribution for the trajectories' duration.

Splitting of Ensemble into Slow and Fast Moving Endosomes
We split ensemble of trajectories into slow and fast moving endosomes using the distance traveled by endosomes: Trajectories which possess active motion have periods of rapid increase or decrease of R ( Figure 1A). Fast trajectories which have active motion are defined as max{R(t)} > and slow trajectories which exhibit only passive motion are defined by max{R(t)} < . Here, max{R(t)} denotes the maximum values of R(t) attained in the time interval (0, t) and is the threshold. We choose the threshold = 0.25 µm. In the Appendix, we show that changing the threshold to = 0.2 µm in the splitting does not qualitatively change the results. Therefore, we define fast moving endosomes as those that, in the time interval (0, t), experienced at least one period of active motion and the maximum distance travelled from the origin exceeds the threshold of = 0.25 µm. Otherwise, an endosome is defined as slow moving. Small variations of the threshold value do not affect the EMSDs of slow and fast moving endosomes, which suggests that the splitting method is robust ( Figure A3). Changing the splitting threshold from max{R(t)} = 0.25 µm to max{R(t)} = 0.2 µm, the increase of the number of slow trajectories was 12%. Therefore, in addition to the method of splitting trajectories which uses the minimum travelled distance, we also tested a second method, which makes use of the time-dependent Hurst exponent H(t) neural network (NN) estimate at the single trajectory level [10]. The procedure is as follows: (1) estimate the time-dependent anomalous exponent α NN using the NN; (2) if the anomalous exponent α NN is superdiffusive α NN (t) > 1 for more than 4 consecutive time points, the endosome is considered as fast moving. Otherwise, the endosome is labeled as slow moving (see Figure A2). The correct implementation of the NN procedure requires a minimum time window [10] that is larger than the duration of some of the endosomal trajectories. Hence, short trajectories were discarded in this analysis. The similarity of the distributions of generalized diffusion coefficients ( Figures A3B and A4B) suggests that the chosen threshold max{R(t)} = 0.25 µm was reasonable. Alternative methods of binary classification could be performed using the first passage probability analysis [7] or implementing the normalized radius of gyration of each trajectory [42].

Ensemble and Time Averaged Mean Squared Displacements
From the two dimensional experimental trajectories r(t) = {x(t), y(t)}, we calculated the ensemble-averaged mean squared displacement (EMSD) as where l is the length scale which we choose l = 1 µm, where the angular brackets denotes averaging over an ensemble of trajectories, A = ∑ N i=1 A i /N and N is the number of trajectories in the ensemble.
By fitting the EMSD to power law functions, the anomalous exponent α and the generalized diffusion coefficient D α can be extracted using where α and D α are constants which characterize averaged transport properties of ensemble of endosomal trajectories. The time scale τ = 1 s and the length scale l = 1 µm are introduced in order to make the generalized diffusion coefficient D α dimensionless. The time-averaged mean squared displacement (TMSD) of an individual trajectory {x i , y i } of a duration T is calculated as: where l is the length scale, for which we chose l = 1 µm, and TMSDs of individual trajectories are averaged further over the ensemble of trajectories to get the ensemble-time-averaged MSD (E-TMSD): where the angular brackets denotes averaging over an ensemble of trajectories as before.

Local Analysis of Endosomal Trajectories
The time-local statistical analysis was implemented as follows. We considered only the portion of a single endosomal trajectory within a window of size W and centered around the time t, i.e., (t − W/2, t + W/2). We calculated the TMSD within this chunk of trajectory only: this is the reason for the acronym L-TMSD, i.e., the local TMSD. The experimental detection of the endosomal motion is achieved with the frame rate 1/∆t s −1 , hence t = i∆t (here, i = 0, 1, 2, . . . is the time index) and W = N∆t, with N > 10. The first 10 points of the L-TMSD were fitted with the power-law function where t = 10∆t. α L (t) and D α L (t) are the local anomalous exponent and generalized diffusion coefficient, respectively. We iterate this procedure by shifting the time window of a single ∆t (i → i + 1) until the end of the experimental endosomal trace, thus obtaining α L (t) and D α L (t) along the entire trajectory. Notice that α L and D L are not constants in time and they vary, being local properties of each endosomal trajectory.

The Time and Ensemble-Time Averaged Velocity Auto-Correlation Functions
The time averaged auto-correlation function (TVACF) along a single trajectory is defined as: where . TVACFs of individual trajectories are averaged further over the ensemble of trajectories to get the ensemble-time averaged VACF (E-TVACF): where the angular brackets denotes averaging over an ensemble of trajectories. The velocity autocorrelation function was suggested as a tool to distinguish between subdiffusion models [43].

Results
We split the ensemble of endosomes into slow and fast moving vesicles using the two methods described above (see Methods, Figures 1A and A2). For both slow and fast endosomes, the EMSDs and E-TMSDs show similar behavior, which suggests ergodicity (see Methods and Figure 1B). MSDs of slow endosomes are not increasing in time, which confirms that these trajectories have no active periods of motion. Surprisingly, we found that both EMSDs and E-TMSDs of slow endosomes are decreasing functions of time, which to our knowledge has never been observed before. We explain this behavior in terms of the coupling between the average diffusivities of slow trajectories and their duration (see Figure 4 and the discussion below). Conversely, MSDs of fast endosomes are increasing functions of time in the intermediate time scale (0.2, 2) s. The anomalous exponent extracted from EMSD or E-TMSD of fast endosomes is α 1, smaller than the anomalous exponent obtained by considering all trajectories without distinction into fast or slow, i.e., α 1.26. Notice that two subdiffusive regimes characterize the MSD time behavior for fast and all trajectories. The first, at small time scales (t ≤ 10 −1 s), can be attributed to the measurement errors [44][45][46]. The second, at longer time scales (t > 10 s), was shown to be spurious and originate from the coupling of the trajectories' duration and their diffusivities [40,47]. We suggest that, due to this coupling, the anomalous exponents deduced from the powerlaw fit of EMSD and E-TMSD, do not capture the essential characteristics of the endosome superdiffusive motility, nor shed light on its fundamental aspects. Therefore, to reveal the effect of the duration of trajectories on the statistical analysis, we consider only trajectories longer than a certain threshold T [40]. Figure 2A,B shows the EMSDs and E-TMSDs of slow and fast endosomes, considering only experimental trajectories with duration longer than T seconds (2 or 8 s). Unlike the slow moving endosomes, the MSDs of fast vesicles ( Figure 2B) present similar qualitative behaviors by choosing T = 2 s, T = 8 s or no T at all (all the fast molecules considered as in Figure 1B, black curve). However, in the intermediate regime, the superdiffusive behavior becomes more and more apparent, ∝ t 1.26 , and stable. In [40] we found that this process is described by the space-time heterogeneous FBM with the Hurst exponent H that randomly switches between persistent H > 0.5 and anti-persistent regimes H < 0.5, together with the coupling between the diffusivity and duration of trajectories which account for spurious subdiffusion at longer time scales. Moreover, the EMSD curves obtained for T = 2 s and T = 8 s deviates considerably from the corresponding E-TMSD curves. The MSDs of slow endosomes (Figure 2A) display very different, but ergodic, behavior. For 0.01 < t < 2 s, the MSDs of all slow endosomes decreases in time. On the other side, the MSDs of the sub-ensembles of slow endosomes with T = 2 s and T = 8 s reveal subdiffusive trends with α∼0.5. As in the case of fast moving endosomes, we argue that this behavior is due to the coupling between the diffusivity and duration of trajectories. Therefore, we attempt to confirm this hypothesis, by performing simulations of an ensemble of heterogeneous FBM trajectories with constant Hurst exponent H = 0.25 (see Figure 4).
The velocity auto-correlation functions (VACF) also confirm the effectiveness of this simple threshold splitting ( Figure 3A,B). Indeed, slow and fast endosomes have very different VACFs. Ensemble-time averaged VACFs (E-TVACFs) of fast endosomes ( Figure 3B) are positive as expected for superdiffusive motion. In contrast, E-TVACFs of slow endosomes have negative dips at t = τ and approach zero from negative values ( Figure 3A). Such behavior is characteristic of FBM and the generalized Langevin equation but cannot be reproduced by the CTRW model [3]. To verify that heterogeneous FBM describes slow moving endosomes, we simulated an ensemble of hFBM trajectories. Individual hFBM trajectories were simulated with constant Hurst exponent H = 0.25. For standard FBM, this would correspond to subdiffusive MSDs, r 2 (t) ∼t 2H ∼ t 0.5 . The duration of hFBM trajectories was drawn from the powerlaw distribution φ(T)∼T −1.85 , in accordance with the experimental evidence [40]. The generalized diffusion coefficients were chosen inversely proportional to the duration of trajectories, i.e., D∼T −0.6 . As shown in Figure 4, the EMSDs of hFBM trajectories agree well with the experimental data. We next implemented the local analysis [40] to better characterize the slow and fast endosomal dynamics. We calculated the local TMSDs (L-TMSD) for each experimental trajectory at various times t ( Methods). From the fit of L-TMSD to Equation (8), we extracted the local anomalous exponents α L (t) and the local generalized diffusion coefficients D α L (t) for slow and fast endosomes separately. The local anomalous exponents α L (t) and the local generalized diffusion coefficients D α L (t) appear to be positively correlated both for slow and fast endosomes (see Figure A6). The origin of these correlations is not known and will be investigated in future publications. In [40], we found that PDFs of local anomalous exponents and local generalized diffusion coefficients do not depend on the window size or the time t (stationary) and are best fitted with exponential and power law functions, respectively.
The PDFs of α L and D α L for slow and fast endosomes are shown in Figure 5A,B. In both cases, the PDFs of α L follow an exponential distribution, while those of D α L are best fitted with a power-law. However, the parameters characterizing the distribution shapes are very different. Furthermore the parameters for the fast endosomes' PDFs coincide with those found by considering all experimental trajectories [40]. This is in agreement with a heterogeneous FBM model of endosomal transport [40], which describes the endosome motion as FBM with non-constant Hurst exponents. Finally, we calculated propagators of experimental trajectories for slow and fast endosomes ( Figure 6). Using the power-law forms of distributions of local generalized diffusion coefficients of slow p S (D L )∼(D L ) −1−γ S and fast p F (D L ) ∼ (D L ) −1−γ F endosomes with γ S 1.7 and γ F 0.5 ( Figure 5), we fit the propagators with the propagators of hFBM, PDF(ξ)∼|ξ| −1−2γ with γ = γ S for slow endosomes and γ = γ F for fast endosomes (see Supplementary Note and [33]). For slow endosomes ( Figure 6A), we also compare the experimental PDFs with the analytical propagator for obstructed diffusion in two dimensions, ξ −0.108 exp(−|ξ| 1.65 ) [48].

Discussion
In this paper, we extend our investigation of the heterogeneous intracellular transport of endosomes based on the local analysis of experimental trajectories [40]. Individual endosomes move for long distances in a heterogeneous way with short bursts of directed motility, interspersed with periods of subdiffusive motion [49,50]. The heterogeneous character of this motion is also manifested as some endosomes are less motile than others. Some endosomes look as if they are jiggling in one position for the whole period of observation. Therefore, we split the ensemble of trajectories into slow and fast moving endosomes. The distinct time behavior of mean squared displacements and velocity autocorrelation functions confirm the effectiveness of these methods. The splitting allowed us to study passive subdiffusive and active superdiffusive transport of endosomes separately.
Comparing the behavior of fast endosomes (MSDs, VACFs and propagators) to the behavior of the entire ensemble, we find that they are most consistent with FBM models [40]. Therefore, we conclude that fast endosomes follow heterogeneous FBM [40]. The ergodicity (Figure 2A) and the VACF ( Figure 3A) suggest that slow endosomes are also described by the hFBM or heterogeneous generalized fractional Langevin equation motion. For slow endosomes, crowding and obstruction effects could also lead to subdiffusive behavior [2,4]. It is known that obstructed diffusion has many similarities with FBM such as stationarity of the increments and the equivalence of the time and ensemble MSDs [48,51]. The propagators provide a clear way to distinguish obstructed diffusion from FBM. Therefore, we calculated propagators of experimental slow endosomes and compared them with analytical prediction for the propagator of obstructed diffusion and prediction of heterogeneous fBM. The results shown in Figure 6 indicate that slow endosomes follow hFBM at longer time scales, while on smaller scales obstructed diffusion likely contributes to their subdiffusive behavior as well. Crowding effects remain as a possible source of anomalous diffusion of slow endosomes. Recently, in numerical simulations, lipids in crowded conditions of the membrane were shown to be multifractal and anomalous. The dynamics was no longer described by the mechanism consistent with the fractional Langevin equation or by any single known mechanism. Instead, the motion was found to be non-Gaussian and heterogeneous, yet maintains its ergodic properties [52], which is similar to what we observed for experimental trajectories of slow endosomes.
Both slow and fast endosomal trajectories are found to be highly heterogeneous in space and time. The spatial heterogeneity in the form of coupling between endosome diffusivity and duration of endosome trajectory explains the behavior of the MSDs. Longer trajectories have smaller generalized diffusion coefficients since in experiments slowly moving endosomes with smaller diffusion coefficients stay longer in the field of view, having longer durations. For slow and fast endosomes, we can conclude that EMSD and E-TMSD are not adequate to describe the large heterogeneity exhibited in space and time. Therefore, we applied a time local analysis of individual trajectories.
From the local analysis, we found that slow and fast endosomal trajectories are both characterized by exponentially distributed anomalous exponents and power-law distributed generalized diffusion coefficients. However, the parameters of these distributions are different. Although the factors that cause the power-law distributed generalized diffusion coefficients for slow and fast endosomes could be different, some common factors can exist. One of them could be the scale free properties of endosomal networks [53]. Hence, the differences in endosome diameters could generate distinct diffusive properties intrinsic to each endosome. Heterogeneous diffusion generated by the fluctuations of molecular size was found in single-molecule experiments within the cell [14,18,42]. Another common factor promoting power-law distributions of generalized diffusion coefficients could be non-specific interactions with the endoplasmic reticulum or other organelles and large intracellular structures. Recently, non-specific interactions were shown to generate heterogeneous diffusion of nanosized objects in mammalian cells [47].
Our analysis of endosomal transport would be valuable for both fundamental cell biology and nanomedicine applications such as drug and gene delivery. In these applications, nanoparticles are often used as cargo-carrying vesicles, which in turn utilize the endosomal network for their intracellular transport. For example, gold nanoparticles were shown to cluster inside endosomes and move via sub-and superdiffusion [54]. Our results would also be useful for the nanoparticle enhanced radiation therapy of cancer [55][56][57] where clusters of nanoparticles inside endosomes are used for dose enhancement.
In the future, we expect microscopy techniques will improve in tandem with tracking algorithms, providing datasets with larger ranges of time scales and improved res-olution. Thus, further subclassification of ensembles of endosomal tracks (beyond the binary fast and slow separation) will become possible towards the ultimate goal of single molecule specificity. Increasing the dynamic range (to submillisecond time scales) will allow the stepping motion of the motor proteins (kinesin and dynein) attached to microtubules to be connected with the spectra of α and D α for the fast moving endosomes at a fundamental level.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.   Figure A2: Two splitting methods used to separate endosome trajectories into slow and fast moving endosomes. The first method uses the maximum distance traveled R(t). The second method uses the time-dependent anomalous exponent H(t) estimated with the neural network. An example of the two trajectories is shown, which were classified as slow and fast by both methods. See the main text for details of the methods.    Figure A4: Comparison of distributions of anomalous exponents α NN and generalized diffusion coefficients D NN and local anomalous exponents α L and D L of slow moving endosomes. Anomalous exponents α NN were estimated using a neural network with window size 0.26 s. The generalized diffusion coefficients D NN were estimated by fitting the local TMSD of the trajectory with the power law D NN t α NN . The distribution of α NN has a maximum of 0.6 and decays faster than the distribution of α L . This may be because many short trajectories are missing in the NN analysis, since the NN could analyze trajectories with durations longer than its window size [10]. The distributions of generalized diffusion coefficients (right panel), on the other hand, are similar to each other.