# Unravelling Heterogeneous Transport of Endosomes

Nickolay Korabel  Department of Mathematics, The University of Manchester, M13
9PL, UK<EMAIL_ADDRESS>Daniel Han  Department of
Mathematics, The University of Manchester, M13 9PL, UK  School of Biological
Sciences, The University of Manchester, M13 9PT, UK  Biological Physics,
Department of Physics and Astronomy, The University of Manchester, M13 9PL, UK
Alessandro Taloni  CNR - Consiglio Nazionale delle Ricerche, Istituto dei
Sistemi Complessi, via dei Taurini 19, 00185 Roma, Italy  Gianni Pagnini  BCAM
- Basque Center for Applied Mathematics, Mazarredo 14, E-48009 Bilbao, Basque
Country - Spain  Ikerbasque - Basque Foundation for Science, Plaza Euskadi 5,
E-48009 Bilbao, Basque Country - Spain  Sergei Fedotov  Department of
Mathematics, The University of Manchester, M13 9PL, UK  Viki Allan  School of
Biological Sciences, The University of Manchester, M13 9PT, UK  Thomas Andrew
Waigh  Biological Physics, Department of Physics and Astronomy, The University
of Manchester, M13 9PL, UK<EMAIL_ADDRESS>
###### Abstract

A major open problem in biophysics is to understand the highly heterogeneous
transport of many structures inside living cells, such as endosomes. We find
that mathematically it is described by spatio-temporal heterogeneous
fractional Brownian motion (hFBM) which is defined as FBM with a randomly
switching anomalous exponent and random generalized diffusion coefficient.
Using a comprehensive local analysis of a large ensemble of experimental
endosome trajectories ($>{10}^{5}$), we show that their motion is
characterized by power-law probability distributions of displacements and
displacement increments, exponential probability distributions of local
anomalous exponents and power-law probability distributions of local
generalized diffusion coefficients of endosomes which are crucial ingredients
of spatio-temporal hFBM. The increased sensitivity of deep learning neural
networks for FBM characterisation corroborates the development of this multi-
fractal analysis. Our findings are an important step in understanding endosome
transport. We also provide a powerful tool for studying other heterogeneous
cellular processes.

##  I Introduction

In eukaryotic cells, endosomes play a major role in sorting and transporting
proteins and lipids that are taken in from the cell surface and need to be
delivered to lysosomes for degradation. They also sort and redistribute
material back to the cell surface or to the Golgi apparatus. Viruses and gene
delivery vehicles exploit endosomal transport during infection pathways and
transfection respectively [1]. Therefore, to design new antiviral drugs and
develop efficient gene therapies, it is imperative to study the mechanisms of
endosomal transport. The movement of endosomes along microtubules is powered
by attachment to cytoplasmic dynein and various kinesin motor proteins [2, 3].
Individual endosomes travel long intracellular distances ($\sim 1-10$ $\mu$m)
in short bursts of directed motility, interspersed with periods of diffusive
[4, 5] and sub-diffusive [6] motion. It is likely that the heterogeneous
character of endosome movement is influenced by cargo sorting and membrane
fission [3, 4, 5, 7].

In previous works, the movement of endosomes inside living cells has been
shown to follow anomalous diffusion [1, 6, 8, 9, 10], which is defined by the
non-linear growth of the mean squared displacement (MSD),
$\left<r^{2}\right>\sim t^{\alpha}$ ($\alpha\neq 1$) [11]. The anomalous
exponent, $\alpha$, characterizes the nature of the anomalous diffusion as
slower (sub-diffusion, $\alpha<1$) or faster (super-diffusion, $\alpha>1$)
than Brownian motion (diffusive, $\alpha=1$). The sub-diffusion exhibited in
endosomal movement can be caused by different mechanisms including: the
cytoplasmic viscoelasticity produced by microtubules and organelles (e.g., the
endoplasmic reticulum network); caging and crowding effects; and temporal
binding. In contrast, the super-diffusive motions of endosomes are driven by
molecular motors. Various anomalous diffusion models have been proposed to
address different mechanisms of biological motility [11, 12, 13, 14, 15, 16,
17, 18, 19, 20, 21]. Fractional Brownian motion (FBM), continuous time random
walk (CTRW) and fractional Langevin equation (FLE) are among the most popular
models. However, many previous models only use a constant generalized
diffusion coefficient and a constant anomalous exponent to describe endosomal
movement. This conflicts with the biological evidence of spatio-temporal
heterogeneity in endosomal dynamics generated by the numerous intracellular
processes that give rise to sub-diffusion and super-diffusion [3, 4, 5, 6, 7].

Recently, it has become clear that anomalous diffusion with a constant
anomalous exponent and a constant generalized diffusion coefficient fail to
adequately describe intracellular transport. With improved microscopy imaging
and tracking methods, the intrinsic spatial and temporal heterogeneity within
individual trajectories of biological processes can be elucidated, for
example: diffusion of proteins and lipids in the cell membrane [22]; colloidal
beads diffusing along linear tubes [23, 24]; the motion of micron-sized beads
in amoeba [25]; labeled messenger RNA molecules in living E. coli and S.
cerevisiae cells [26]; individual quantum dots moving in the cytoplasm of
living mammalian cells [27]; Rab5 and SNX1 tagged endosomal transport [6]; and
lysosomal dynamics in living cells [28, 6]. Heterogeneous Brownian transport
processes have only had sparse coverage in the literature including the broad
distribution of diffusivities in the dynamics of the pathogen-recognition
receptor dendritic cell-specific intercellular adhesion molecule 3-grabbing
nonintegrin [29]; and the exponential probability distribution of
diffusivities and the Laplace probability distribution of displacements in the
motion of RNA molecules in E. coli and S. cerevisiae cells [23, 24, 26]. These
discoveries led to the development of a new class of mathematical models
involving heterogeneous diffusion described as “Brownian yet non-Gaussian
diffusion”. Among these models, the most popular are the so-called
superstatistical Brownian motion [30, 31] and FBM [33, 32], which include a
population of diffusion coefficients, position dependent diffusion
coefficients [38, 39], and models involving diffusing diffusivities [34, 35,
36, 37] that were recently extended to FBM [42]. Furthermore, extending
heterogeneity to the consideration of the anomalous diffusion, power-law
distributions of both diffusion coefficients and displacements were found in
the motion of H-NS proteins [40] and recently modelled through the joint
fluctuations of both the anomalous diffusion exponents and the diffusion
constants of the FBM [41]. A relatively sparse literature also exists on
position and density dependent anomalous exponents [43, 44, 45] and
distributed-order diffusion equations [46].

In this article, we characterized anomalous transport in experimental
trajectories of endosomes inside eukaryotic cells. We have used traditional
statistical methods and a recently developed neural network (NN) analysis [6]
based on a deep learning feedforward neural network. Both approaches reveal
heterogeneous behaviour on a single trajectory level and the NN method
provides a dramatic improvement in sensitivity. Using a large ensemble of
experimental trajectories, we find an exponential probability distribution of
anomalous exponents, a power-law probability distribution of local generalized
diffusion coefficients, and a power-law probability distribution of
displacements, $X$ and displacement increments, $\Delta X$, for the endosomes.
We show that endosomal movement is described by heterogeneous FBM (hFBM). In
particular, the superstatistical ensemble [30, 31, 37] of FBM trajectories
which accounts only for the spatial heterogeneity explains the power-law
probability distributions of displacements and displacement increments of
experimental endosomal trajectories. However, it does not explain the mean
squared displacements of ensemble of experimental endosomes. On the contrary,
the spatio-temporal heterogeneous FBM with random switching from persistent
($H>0.5$) to anti-persistent regime ($H<0.5$) (in the mathematical literature
known as multifractional FBM [47]) with random Hurst exponent $H$ and random
generalized diffusion coefficient explains all experimental results. Our
findings and methodologies shed new light on the underlying mechanism of
endosomal transport in eukaryotic cells and improve our understanding of
complex heterogeneous intracellular processes.

##  II Results

Figure 1: The probability distribution of endosomal displacements is a power-
law. a. PDFs of displacements $X=x(t)-x(0)$ extracted from experimental
trajectories at $t=0.8,1,2,4,8$ seconds. b. PDFs of displacements scaled by
the standard deviation $\sigma_{X}$ in a log-linear scale. c. Same as in b in
a log-log scale. PDFs of $X/\sigma_{X}$ are best fitted with a power law
$(X/\sigma_{X})^{-2}$ (the solid curve) contrary to Gaussian (the dashed
curve) and Laplace probability distributions (the dashed-dotted curve).

Figure 2:  Probability distribution of endosomal displacement increments is a
(truncated) power-law. a. PDFs of displacement increments $\Delta X=x(t+\Delta
t)-x(t)$ for time intervals $\Delta t=0.8,1,2,4,8$ seconds. b. Distribution of
increments scaled by the standard deviation $\sigma_{dX}$ in a log-linear
scale. c. Same as in b in a log-log scale. The solid curve represents the fit
with the truncated power law, $\exp(-0.15\Delta X/\sigma_{\Delta X})(\Delta
X/\sigma_{\Delta X})^{-2}$. The dashed curve is a fit with a Gaussian
probability distribution and the dashed-dotted curve is a fit with a Laplace
probability distribution.

We studied 103,361 experimental trajectories of early endosomes in a stable
MRC5 cell line expressing GFP-Rab5, obtained from tracking wide-field
fluorescence microscopy videos (see [6] for experimental details). The
endosomes were tracked using an automated tracking software (AITracker, based
on a convolutional neural network) [48]. The duration of all trajectories,
$T$, has a good fit to a power law distribution, $\phi\left(T\right)\sim
T^{-1.85}$ (see Appendix B, Fig. S1). This distribution 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 [49] leading to the emergence of the power-law probability
distribution for the trajectories’ duration.

###  II.1 Endosomal trajectories have power-law distributed displacements and
increments

Figure 3:  EMSD and E-TMSD of experimental trajectories reveal ergodic
behaviour. Longer trajectories gradually transition to sub-diffusive
behaviour. a. The ensemble-averaged MSD (EMSD) and the time-averaged averaged
over all trajectories (E-TMSD) as a function of time interval, t. b. EMSD as a
function of time interval of trajectories longer than $T$ seconds. The super-
diffusive behaviour at the intermediate time scale is conserved with the
anomalous exponent greater or approximately equal to $\alpha=1.26\pm 0.02$.
Trajectories longer than $8$ s (curves $T>8$ s and $T>15$ s) become sub-
diffusive at longer time scales with the anomalous exponent $\alpha\sim 0.9\pm
0.01$. c. Distributions of anomalous exponents $\alpha_{i}$ and d
distributions of generalized diffusion coefficients $D_{\alpha_{i}}$ obtained
by fitting the TMSDs of single trajectories at different times indicated in
the legends using the time window $W=0.4$ seconds. Notice that
$D_{\alpha_{i}}$ is dimensionless (see Methods).

From the two-dimensional experimental trajectories
$\mathbf{r}\left(t\right)=\\{x\left(t\right),y\left(t\right)\\}$, we extracted
the displacements, $X=x\left(t\right)-x\left(0\right)$,
$Y=y\left(t\right)-y\left(0\right)$ and displacement increments, $\Delta
X=x\left(t+\Delta t\right)-x\left(t\right)$, $\Delta Y=y\left(t+\Delta
t\right)-y\left(t\right)$ and calculated their distributions. Endosomal
trajectories turned out to be isotropic since the Probability Density
Functions (PDFs) of the $X$ and $Y$ components of the displacements are
similar. Therefore, we show PDFs of the $X$ component only. The resulting
plots are shown in Fig. 1 and 2. We find that endosomes display strongly non-
Gaussian power-law distributions of displacements and displacements increments
along trajectories. The PDFs of scaled displacements $X/\sigma_{X}$ where
$\sigma_{X}$ is the standard deviation of X, collapse onto a single master
curve for $t>0.8$ s (Fig. 1b, c) and are well fitted with the power-law
$\left(X/\sigma_{X}\right)^{-2}$. The form of the PDFs of endosomal
displacements does not change in time within the time interval $0.8<t<8$ s.
This is different from other observations of heterogeneous systems where it
was found that the Laplace distribution of displacement changes its form to a
Gaussian distribution at longer times [23, 24]. For time $t>8$ s, the PDFs of
displacements become noisy due to a small number of remaining trajectories.
Distributions of displacement increments along trajectories $\Delta X$ exhibit
similar power-law form as shown in Fig. 2a. Scaled with the standard deviation
$\sigma_{\Delta X}$, the PDFs of displacement increments collapse into a
single curve, which has a truncated power law behaviour with the same exponent
as the PDFs of the displacements (Fig. 2b, c). Non-Gaussian distributions of
displacements and displacement increments are the characteristic feature of
heterogeneous diffusion processes (for a recent review see [37]). Commonly,
Laplace probability distribution of displacements in heterogeneous diffusion
is observed,
$\mathrm{PDF}\left(X\right)\sim\exp{\left(-\left|X\right|/2\sigma^{2}\right)}$
[37]. For endosomes, we find power-law probability distributions of
displacements and increments. To our knowledge, power-law probability
distributions were previously reported only for trajectories of individual
histone-like nucleoid-structuring proteins [40].

###  II.2 Experimental endosomal trajectories display ergodic but anomalous
mean-squared displacements

The ensemble-averaged MSD (EMSD) of experimental trajectories shows the
initial sub-diffusive-like behaviour (Fig. 3) which could be attributed to the
static measurement errors [50]. At time scale $0.2<t<2$ s, molecular motors
make the growth of EMSD super-diffusive with the anomalous exponent
$\alpha\sim 1.26\ \pm 0.02$ and the ensemble generalized diffusion coefficient
$D_{\alpha}\sim 0.082\pm 0.002$ $\mu$m${}^{2}/s^{\alpha}$. Similar behaviour
was observed for the transport of micron-sized beads along microtubules [51]
and endosomes filled with labeled nonviral DNA-containing polyplexes ranging
from $100$ to $200$ nm in diameter [10] and for single cell trajectories [52].
At longer time scales ($t\ >\ 2$ s), the EMSD become sub-diffusive (Fig. 3a)
with the anomalous exponent $\alpha\sim 0.5\pm 0.1$. The time averaged MSD
along single trajectories (averaged over the ensemble of all trajectories,
E-TMSD) is similar to EMSD behaviour. EMSD and E-TMSD are slightly different
both in the super-diffusion regime ($t<2$ s) and in the sub-diffusive regime
($t>2$ s), perhaps due to limited statistics. This ergodic behaviour
(equivalence of EMSD and E-TMSD) suggests that the endosomal movement is best
described by the FBM process. However, the sub-diffusion at long times does
not fit the standard FBM model with constant $\alpha$ and $D_{\alpha}$. As we
show below, this apparent sub-diffusion regime is spurious [49].

###  II.3 Time heterogeneous FBM explains spurious sub-diffusive behaviour of
endosomal ensemble

To investigate further the behaviour of MSDs, we calculated EMSDs by
considering ensembles of experimental trajectories with duration longer than a
certain time $T$ (Fig. 3b). The super-diffusive behaviour at the intermediate
time scale is conserved with approximately the same anomalous exponent,
$\alpha=\ 1.26\ \pm 0.02$. For longer trajectories, the generalized diffusion
coefficient becomes smaller, which again suggests that slower moving endosomes
remain in focus. Notice that the long-time behaviour of EMSD is different from
the apparent sub-diffusion observed in Fig. 3a ($\alpha\sim 0.5\ \pm 0.1$). In
fact, for trajectories longer than $T=8$ seconds, EMSDs remain super-diffusive
at intermediate time scales much longer and well into the apparent sub-
diffusive regime of the EMSD obtained from all trajectories. Moreover, the
longtime dynamics gradually transition to a sub-diffusive regime which is
different from the sub-diffusive exponent of the EMSD of all trajectories.
Trajectories of duration longer than $15$ seconds become sub-diffusive with
the exponent $\alpha\sim 0.9\ \pm 0.05$ (Fig. 3b). The behaviour of the
E-TMSDs obtained from trajectories of different duration $T$ is similar to
EMSD, therefore they are not shown. To check that the apparent sub-diffusive
regime at long time scales could originate from the super-diffusive movement
of endosomes, we simulated an ensemble of independent super-diffusive FBM
trajectories, each with the same Hurst exponent $H\ =\ 0.6$, but with a
different generalized diffusion coefficient. The generalized diffusion
coefficients were chosen to be correlated with the time length of trajectory -
for longer trajectories we chose smaller generalized diffusion coefficients.
This choice is reasonable since in the experimental setup fast endosomes with
bigger diffusivities quickly leave the focus area and therefore have shorter
duration. On the contrary, slowly moving endosomes with smaller diffusion
coefficients stay longer in focus and have longer duration. If it were the
case of a standard FBM motion, the EMSD would grow super-diffusively with the
anomalous exponents $\alpha=2H=1.2$. Hence, we simulated an ensemble of
heterogeneous FBM (hFBM) trajectories of different length $T$. The numerical
results are shown in Fig. 4, confirming that the sub-diffusion regime emerges
at long time scales even though all hFBM trajectories were purely super-
diffusive. This result shows that time heterogeneous FBM motion explains the
spurious sub-diffusion behaviour of MSDs at longer time scales. It also
implies that one cannot simply conclude the overall behaviour of the ensemble
of endosomes just by looking at the EMSD or the E-TMSD behaviour, since both
are masked by the apparent sub-diffusion regime at long time scales caused by
the heterogeneity in the trajectory lengths. Therefore, it is necessary to
look closer into the local dynamics of individual endosomal trajectories.

Figure 4: Simulations of ensemble of hFBM trajectories with coupling between
diffusion coefficients and duration of trajectories explains the apparent sub-
diffusive behaviour of experimental endosome trajectories at longer time
scales. Different curves correspond to EMSDs calculated for simulated hFBM
trajectories of duration longer than $T$ seconds. Individual hFBM trajectories
were simulated with constant Hurst exponent $H=0.6$ and constant generalized
diffusion coefficients. For standard FBM motion, this corresponds to the
super-diffusive growth of EMSD with the anomalous exponent $\alpha=1.2$ shown
by the dashed line. For hFBM, the duration of trajectories was random and the
generalized diffusion coefficient for each trajectory was chosen depending on
its duration. For longer trajectories, the generalized diffusion coefficient
was chosen to be smaller (see Methods for details). The apparent sub-diffusive
behaviour at longer time scales with the anomalous exponent $\alpha=0.5$ is
shown as the dashed-dotted line. The EMSD of experimental endosomal
trajectories is shown for comparison.

Figure 5: Experimental endosomal trajectories display time independent
exponential probability distributions of local anomalous exponents
$\alpha^{L}_{i}$ and power-law probability distributions of local generalized
diffusion coefficients $D_{\alpha^{L}_{i}}$. The local TMSD (L-TMSD) was
calculated in a time window $W=1.2$ s at different time $t$. Different curves
correspond to PDFs of $\bar{\alpha}$ and $D_{\alpha^{L}_{i}}$ calculated at
time $t$ given in the figure, see panel a and b respectively. Notice that
$D_{\alpha^{L}_{i}}$ is dimensionless (see Methods). The dashed lines
represent exponential fit in a and power law fit in b.

Figure 6: NN analysis reveals spatio-temporal heterogeneity of experimental
endosomal trajectories. NN confirms an exponential distribution of anomalous
exponent and a power-law distribution of generalized diffusion coefficients.
Displacements $\mathbf{r}(t)-\mathbf{r}(0)$ (an example for a sample
trajectory is shown in a) were analysed with NN and the local time dependent
Hurst exponent $H(t)$ was estimated. b. $H(t)$ for a sample trajectory
depicted in c. d. Probability density functions of local anomalous exponents
$\alpha_{i}^{NN}=2H$ estimated with the NN. e. PDFs of generalized diffusion
coefficients $D_{\alpha^{NN}_{i}}$. Different curves in d and e correspond to
PDFs estimated at $t=0.2,0.4,0.8,1,2,4,8,10,12,14$ s.

###  II.4 Distributions of time-dependent anomalous exponents and generalised
diffusion coefficients

Using TMSD curves, we calculate the time-dependent anomalous exponent
$\alpha_{i}\left(t\right)$ and the time-dependent generalized diffusion
coefficient $D_{\alpha_{i}}\left(t\right)$ of single experimental endosome
trajectories. We did this by calculating TMSDs of single trajectories first,
and then fitting TMSDs in a time window $\left(t-W/2,t+W/2\right)$ with the
power-law function. From the fit we estimated $\alpha_{i}\left(t\right)$ and
$D_{\alpha_{i}}\left(t\right)$ (see Methods). The distribution of
$\alpha_{i}\left(t\right)$ (Fig. 3c) is exponential for $t<0.4$ s and
(truncated) Gaussian-like at longer time scales. We note that negative
exponents $\alpha_{i}<0$ are fitting artefacts which have no physical meaning.
These artefacts arose from the finite duration of experimental trajectories
and small window size $W$. For larger window size, the fitting artefacts
gradually diminish (see Appendix B, Fig. S2). The peak of the Gaussian fit
(dashed-dotted curve in (Fig. 3c) corresponds to $\alpha_{i}\sim 0.9$ which
agrees with the EMSD behaviour at this time scale (Fig. 3b). PDFs of
$D_{\alpha_{i}}$ (Fig. 3d) are independent of time and are best fitted with a
power law (shown in Fig. 3d as the dashed line). In contrast to the constant
ensemble anomalous exponent $\alpha$ and the constant generalized diffusion
coefficient $D_{\alpha}$, anomalous exponents $\alpha_{i}\left(t\right)$ and
generalized diffusion coefficients $D_{\alpha_{i}}\left(t\right)$ characterize
endosomal motion in time, but also in space. Indeed,
$\alpha_{i}\left(t\right)$ and $D_{\alpha_{i}}\left(t\right)$ contain some
amount of spatial averaging since they were calculated within a finite time
window of the TMSD curves, obtained collecting displacements at the time scale
$t$ from any point in the endosomal trajectory (see the Methods for the TMSD
definition). On top of this, we already noticed that the distributions of the
time-dependent anomalous exponents $\alpha_{i}\left(t\right)$ have fitting
artefacts. Therefore, in the next subsection we proceed to examine the local-
time characteristics of individual trajectories by estimating the anomalous
exponent and generalized diffusion coefficient straight from endosomal
trajectories.

###  II.5 Endosomes display exponential distributions of local anomalous
exponents and power law distributions of generalised diffusion coefficients

To examine the local characteristics of individual trajectories, we calculated
the local TMSDs (L-TMSD) from each experimental trajectory at the time $t$.
This is defined as the TMSD of a single trajectory’s portion contained within
the time window $\left(t-W/2,t+W/2\right)$ [25]. L-TMSDs were fitted with a
power law function (see Methods) and the local anomalous exponents
$\alpha_{i}^{L}\left(t\right)$ and the local generalized diffusion
coefficients $D_{\alpha_{i}^{L}}\left(t\right)$ were estimated. In general,
$\alpha_{i}^{L}$ and $D_{\alpha_{i}^{L}}$ are different from $\alpha_{i}$ and
$D_{\alpha_{i}}$. On the time scale $0.2\ <\ t\ <\ 8$ s, the PDFs of
$\alpha_{i}^{L}$ and $D_{\alpha_{i}^{L}}$ do not depend on time (Fig. 5a, b).
They also do not depend on window size (see Appendix B, Fig. S3). Most
importantly, $\alpha_{i}^{L}$ and $D_{\alpha_{i}^{L}}$ of experimental
endosomal trajectories are best fitted with exponential and power law
functions respectively (Fig. 5a, b). Averaged over all trajectories, the
values of $\alpha_{i}^{L}$ and $D_{\alpha_{i}^{L}}$ at different times are
markedly different from the values of $\alpha$ and $D_{\alpha}$ estimated from
the EMSDs (see Appendix B, Fig. S4). PDFs of local anomalous exponents
$\alpha_{i}^{L}$ (Fig. 5a) show that despite super-diffusive motion being
present at all times, the average local anomalous exponent $\alpha_{i}^{L}$ is
sub-diffusive. The values of the local anomalous exponents and generalized
diffusion coefficients exhibit non-linear positive correlations (see Appendix
B, Fig. S5) ($\alpha_{i}^{L}=a\log{D_{\alpha_{i}^{L}}}+b$ with $a=0.42$,
$b=1.22$ and Pearson’s correlation coefficient $r=0.8$).

Figure 7: Power-law probability distributions of displacements and increments
of experimental endosomal trajectories are the same as for the heterogeneous
ensemble of FBM trajectories. Probability distributions of displacements
$X=x(t)-x(0)$ scaled by the standard deviation $\sigma_{X}$ a and increments
$\Delta X=x(t+\Delta t)-x(t)$ scaled by the standard deviation $\sigma_{\Delta
X}$ b for the hFBM trajectories with the generalized diffusion coefficients
distributed according to a power law. The PDFs of displacements and increments
for experimental trajectories overlap with PDFs for hFBM. The dashed and
dashed-dotted curves are the fits with the power law a, truncated power law b
and Laplace PDFs.

To assess the local characteristics of a single endosome motion using a
different method, we applied the neural network (NN) [6]. Using NN we
estimated the local Hurst exponent $H$ of experimental trajectories which is
related to the local anomalous exponent, $H=\alpha/2$. The advantage of the NN
is that it does not make use the fit of TMSD to estimate $H$ and it is a more
sensitive estimator [6]. Since endosomes MSD appears to be ergodic (see Fig.
3a), we trained the NN on FBM trajectories and estimated the local Hurst
exponents $H\left(t\right)$ for each endosomal trajectory (see [6] for details
of the implementation and training of the NN). An example of a heterogeneous
endosomal trajectory and the local Hurst exponent estimated using the NN is
shown in Fig. 6a-c. The figure shows that the endosome movement switches
between slow anti-persistent ($H<0.5$) and fast persistent movements $H>0.5$.
The PDFs of $\alpha_{i}^{NN}$ estimated using the NN (Fig. 6d) are
qualitatively similar to the distributions of the local anomalous exponents
(Fig. 5a) although the exponential decay of PDFs of $\alpha_{i}^{NN}$ is
faster. Moreover, the PDFs of $\alpha_{i}^{NN}$ display maxima around
$\alpha_{i}^{NN}\sim 0.7$. We stress the fact that, although the power law
fitting of L-TMSD may produce non-physical values of $\alpha_{i}^{L}$
($\alpha_{i}^{L}<0$ or $\alpha_{i}^{L}>2$), the estimate of $\alpha_{i}^{NN}$
is free of this drawback. The NN methods could directly infer only the value
of the anomalous exponent. For the estimate of the generalized diffusion
coefficient $D_{\alpha_{i}^{NN}}$ we proceeded as follows. Handling the same
endosome trajectories portions used in the L-TMSD analysis, we fit the
resulting MSD with the law $D_{\alpha_{i}^{NN}}\;t^{\alpha_{i}^{NN}}$. As a
consequence, $\alpha_{i}^{NN}$ and $D_{\alpha_{i}^{NN}}$ are non-linearly
positively correlated as $\log{\alpha_{i}^{NN}}=c\log{D_{\alpha_{i}^{NN}}}+d$
with $c=0.03$, $d=0.12$ and the Pearson’s correlation coefficient $r=0.54$
(see Supplementary Note, Fig. S4). Importantly, by comparing the distributions
of $D_{\alpha_{i}^{L}}$ and $D_{\alpha_{i}^{NN}}$ (Figs. 5b, 6e), we can
conclude that the power-law decay appears robust to the specific methodology
used. At the same time, stressing the difference among the PDFs of the local
anomalous exponents in Fig. 5a and Fig. 6d and the inherent fitting problems,
we can assert that the NN method is more reliable.

Figure 8: Spatio-temportal heterogeneous FBM model (ST-hFBM) explains MSDs,
power-law probability distributions of displacements and increments of
experimental endosomal trajectories. a EMSD and E-TMSD of ST-hFBM trajectories
compared with the EMSD of experimental trajectories. Different curves
correspond to different duration of trajectories $T$ used to calculate EMSD. b
Probability distributions of displacements $X=x(t)-x(0)$ scaled by the
standard deviation $\sigma_{X}$ and c increments $\Delta X=x(t+\Delta t)-x(t)$
scaled by the standard deviation $\sigma_{\Delta X}$ for the ST-hFBM
trajectories compared with the corresponding PDFs of experimental
trajectories.

###  II.6 Endosome movement is a spatio-temporal heterogeneous fractional
Brownian motion

The ergodic behaviour of the anomalous MSDs of experimental endosomal
trajectories (Fig. 3a) and the stationary (time independent) property of their
scaled by the standard deviation displacements and increments distributions
(Figs. 1b,c and 2b,c) agree with the FBM model. On the other hand, the power
law form of distributions of displacements and increments and the PDFs of
local anomalous exponents and generalized diffusion coefficients are not
captured by the standard FBM. Below we show that heterogeneous FBM (hfBm)
explains all the experimental data for endosomal transport in eukaryotic
cells.

####  II.6.1 Spatially heterogeneous fractional Brownian motion explains PDFs
of displacements and increments, but fails to capture MSDs of the endosome
ensemble

A simple form of heterogeneous FBM behaviour which is only spatially (but not
temporally) heterogeneous could be realised via an ensemble of independent FBM
trajectories with Hurst exponent $H$ and generalized diffusion coefficient
$D_{H}$, different for each trajectory and constant in time [41]. $H$ was
drawn from an exponential probability distribution, while $D_{H}$ was
distributed in a power law fashion in analogy to Fig. 5 and Fig. 6d,e (see
Methods for details). This simplified picture corresponds to the spatial
inhomogeneity of FBM motion, while no time heterogeneity was implemented.
Figure 7 demonstrates an excellent agreement between distributions of
displacements and increments of the spatially heterogeneous ensemble of FBM
trajectories and PDFs obtained from experimental trajectories. Analytical
calculations (see Appendix A and Appendix B, Fig. S6) confirm the existence of
the scaling behaviour (collapse of PDFs on a master curve) of the PDFs of
displacements of hFBM. The velocity autocorrelation functions of the ensemble
of hFBM trajectories also show similar qualitative behaviour compared to the
experimental trajectories (see Appendix B, Fig. S6). However, the EMSD of the
hFBM ensemble does not agree with the EMSD of experimental trajectories (see
Appendix A, Eq. (8) for analytical derivation). This discrepancy arises
because of the simplified picture of the spatially heterogeneous FBM approach
which completely disregards the trajectories’ time-heterogeneity. Therefore,
we ought to consider a spatio-temporal heterogeneous FBM.

####  II.6.2 A spatio-temporal heterogenous FBM best describes experimental
endosome motion

Opposite to the spatially heterogeneous FBM, the Hurst exponent of a single
spatio-temporal heterogeneous FBM (ST-hFBM) trajectory is not constant but
randomly switches in time between anti-persistent $H=0.25$ value and
persistent $H=0.6$ value (see Methods). Such a dichotomous model corresponds
to random switching of endosomal movement between active and passive states of
motion as observed in experiment (see Fig. 6b). Results of simulations of ST-
hFBM qualitatively agree with the experimental data (Fig. 8). The EMSDs,
distributions of displacements and increments of simulated ST-hFBM
trajectories fit well the EMSDs and corresponding distributions of
experimental trajectories. This supports our assumption that the endosomal
motion is best described by a spatio-temporal heterogeneous FBM. The
quantitative deviations suggest that there could be more than two states in
endosome motion. The multiple active states most probably will correspond to
endosomal movement powered with different number/type of molecular motors [58,
59, 60].

In summary, our results show that including spatio-temporal heterogeneity into
the FBM captures essential statistical characteristics of endosomal motion
(Figs. 1, 2, 3). The analysis of local dynamic properties reveals a robust
character of endosomal dynamics and a broad spectrum of local anomalous
exponents and generalized diffusion coefficients. The distributions of local
anomalous exponents and generalized diffusion coefficients do not depend on
time and best described by the exponential and power-law PDFs respectively
(Figs. 5, 6). The dependence of the generalized diffusion coefficient on
trajectory duration in the hFBM also explains the apparent sub-diffusive
behaviour seen at long-time scales (Fig. 4). Moreover, a simple spatially
heterogeneous FBM accurately reproduces the power-law distributions of
displacements and increments of endosomes (Fig. 7). By incorporating both the
spatial and temporal heterogeneity found in experiments (Fig. 6), we show that
the ST-hFBM describes experimental endosomal trajectories. It reproduces the
MSDs (including the ergodic behaviour in terms of the equivalence of EMSD and
E-TMSD), displacement and increment probability distributions of endosomes
(Fig. 8).

##  III Methods

###  III.1 Statistical analysis of experimental trajectories

####  III.1.1 EMSD, $\alpha$ and $D_{\alpha}$

The ensemble-averaged mean squared displacement (EMSD) of 2D experimental
trajectories is defined in a standard way:

$\mbox{EMSD}(t)=\frac{\left<\mathbf{r}\right>^{2}(t)}{l^{2}},$ (1)

where $l$ is the length scale which we choose $l=1$ $\mu$m and

$\left<\mathbf{r}\right>^{2}(t)=\left<(x_{i}(t)-x_{i}(0))^{2}+(y_{i}(t)-y_{i}(0))^{2}\right>,$
(2)

The angled brackets denotes averaging over an ensemble of trajectories,
$\left<A\right>=\sum_{i=1}^{N}A_{i}/N$, where $N$ is the number of
trajectories in the ensemble.

By fitting the EMSD to the power law function, we extracted the anomalous
exponent $\alpha$ and the generalized diffusion coefficient $D_{\alpha}$

$\mbox{EMSD}(t)=4D_{\alpha}\left(\frac{t}{\tau}\right)^{\alpha}.$ (3)

The time scale $\tau=1$ sec and the length scale $l=1$ $\mu$m are introduced
in order to make the generalized diffusion coefficient $D_{\alpha}$
dimensionless [62]. Notice that $\alpha$ and $D_{\alpha}$ are constants which
characterize averaged transport properties of ensemble of endosomal
trajectories.

####  III.1.2 TMSD, E-TMSD, $\alpha_{i}(t)$ and $D_{\alpha_{i}}(t)$

The subscript $i$ indicates that $\alpha_{i}(t)$ and $D_{\alpha_{i}}(t)$ are
extracted for an individual trajectory with index $i$ as described below.
$\alpha_{i}(t)$ and $D_{\alpha_{i}}(t)$ are also calculated at time $t$ (so,
they depend on time) in contrast to $\alpha$ and $D_{\alpha}$ which are
constants calculated for ensemble of all trajectories (see the definition
above). First, we calculate the time-averaged mean squared displacement (TMSD)
of a individual trajectory $\\{x_{i},y_{i}\\}$ of a duration $T$:

$\mbox{TMSD}_{i}(t)=\frac{\overline{\delta^{2}(t)}}{l^{2}},$ (4)

where where $l$ is the length scale which we choose $l=1$ $\mu$m and

$\overline{\delta^{2}(t)}=\frac{\int_{0}^{T-t}\Delta_{i}(t,t^{\prime})dt^{\prime}}{T-t},$
(5)

where

$\Delta_{i}(t,t^{\prime})=(x_{i}(t^{\prime}+t)-x_{i}(t^{\prime}))^{2}+(y_{i}(t^{\prime}+t)-y_{i}(t^{\prime}))^{2}.$
(6)

TMSDs of individual trajectories can be averaged further over the ensemble of
trajectories to get the ensemble-time-averaged MSD (E-TMSD):

$\mbox{E-TMSD}(t)=\left<\mbox{TMSD}_{i}(t)\right>.$ (7)

By fitting the TMSDs to the power law function in a time window
$(t-W/2,t+W/2)$, we extract the anomalous exponent $\alpha_{i}(t)$ and the
generalized diffusion coefficient $D_{\alpha_{i}}(t)$ at time $t$:

$\mbox{TMSD}_{i}(t)=4D_{\alpha_{i}}\left(\frac{t}{\tau}\right)^{\alpha_{i}},\;\;\;(t-W/2,t+W/2).$
(8)

####  III.1.3 Local TMSD (L-TMSD), $\alpha^{L}_{i}(t)$ and
$D_{\alpha^{L}_{i}}(t)$

We calculate the local TMSD (L-TMSD) e.g. TMSD is calculated in a time window
$(t-W/2,t+W/2)$. The window is then shifted along a trajectory and L-TMSD is
calculated in each window. The subscript $i$ indicates that
$\alpha^{L}_{i}(t)$ and $D_{\alpha^{L}_{i}}(t)$ are extracted for an
individual trajectory with index $i$. By fitting the first $10$ L-TMSDs points
to power law functions, we extracted the time dependent local anomalous
exponents $\bar{\alpha}(t)$ and the time dependent local generalized diffusion
coefficients $D_{\bar{\alpha}}(t)$:

$\mbox{L-TMSD}_{i}(t)=4D_{\alpha^{L}_{i}}\left(\frac{t}{\tau}\right)^{\alpha^{L}_{i}}.$
(9)

Here we introduced the time scale $\tau=1$ sec and the length scale $l=1$
$\mu$m in order to make the local generalized diffusion coefficients
dimensionless [62]. Notice the difference between the TMSD fit in Eq. 8 and
the L-TMSD fit in Eq. 9. In Eq. 8 we calculate the TMSD first and then fit it
in the time window $(t-W/2,t+W/2)$. In Eq. 9 we first calculate the L-TMSD
(TMSD in the time window $(t-W/2,t+W/2)$) and then fit the first $10$ points
of L-TMSD.

The time dependence of the EMSD determine the character of the movement of
ensemble of endosomes as normal diffusion for $\alpha=1$, sub-diffusion for
$\alpha<1$ and super-diffusion for $\alpha>1$. Similarly the time dependence
of the L-TMSD determine the character of local movement of a single endosome
as normal diffusion for $\alpha^{L}_{i}=1$, sub-diffusion for
$\alpha^{L}_{i}<1$ and super-diffusion for $\alpha^{L}_{i}>1$.

###  III.2  Simulation of spatially heterogeneous FBM trajectories (Fig. 4)

FBM trajectories were simulated in MATLAB (The MathWorks, Nantick, MA) ffgn
function which uses the Fast Fourier Transform method to generate FBM. Each
trajectory is given a constant Hurst exponent, $H$, and generalized diffusion
coefficient, $D$ [41]. In simulations shown in Fig. 4, the Hurst exponents
were constant ($H=0.6$) for all trajectories. The generalized diffusion
coefficients were assigned from the power law distribution $f(D)=\gamma
D_{min}^{\gamma}D^{-1-\gamma}$ with $\gamma=0.5$ as it was observed for
endosomal trajectories (see Fig. 5b). The PDF of $D$ was truncated at
$D_{min}=0.001$ and $D_{max}=7$. To reproduce the sub-diffusive behaviour of
MSDs at longer time scales, the duration of the hFBM trajectories $T$ was
coupled with their generalized diffusion coefficients. Specifically, we chose
$T=5D^{-0.6}$ s which was truncated by $T_{min}=0.04$ s and $T_{max}=70$ s.
So, longer trajectories have smaller generalized diffusion coefficient $D$.

###  III.3  Simulation of spatially heterogeneous FBM trajectories (Fig. 7)

FBM trajectories were simulated with the same method as detailed above. In
contrast with simulations shown in Fig. 4, in Fig. 7 the Hurst exponents were
constant in time but random for each trajectory. They were assigned from the
exponential distribution $p(H)=3.72\exp(-3.72H)$ which corresponds to the PDF
of anomalous exponents $\mbox{PDF}(\alpha)=1.86\exp(-1.86\alpha)$ observed for
experimental endosomal trajectories (see Fig. 5a). Hurst exponents are
bounded, $H\in[0,1]$. The generalized diffusion coefficients were assigned
from the power law distribution $f(D)=\gamma D_{min}^{\gamma}D^{-1-\gamma}$
with $\gamma=0.5$ as it was observed for endosomal trajectories (see Fig. 5b).
The PDF of $D$ was truncated at $D_{min}=0.1$ and $D_{max}=100$.

###  III.4 Simulation of spatio-temporal heterogeneous FBM trajectories (Fig.
8)

In contrast to the spatially heterogeneous FBM trajectories, the anomalous
exponent, $H$, and the generalized diffusion coefficient, $D$, of a single
spatio-temporal heterogeneous FBM trajectory are not constant in time. $H$
randomly switches between $H=0.2$ and $H=0.6$ with switching rates $4/$s and
$1/$s. At $t=0$, we choose $H=0.2$ or $H=0.6$ with equal probability. The
diffusivities $D$ of the persistent regime $H=0.6$ were assigned from a power
law distribution $f(D)=\gamma D_{min}^{\gamma}D^{-1-\gamma}$ with $\gamma=0.5$
as it was observed for experimental endosomal trajectories (see Fig. 3). The
PDF of $D$ was truncated at $D_{min}=0.0004$ and $D_{max}=20$. The
diffusivities of the anti-persistent regime ($H=0.2$) were chosen as $0.2D$.
In order to reproduce the apparent subdiffusive behaviour of the experimental
MSDs at longer times (see the discussion in the Results and Fig. 4), we couple
the duration of trajectories with their diffusivities. We chose $T=13D^{-0.6}$
s which was truncated by $T_{min}=0.04$ s and $T_{max}=40$ s.

##  IV Discussion

It has been known for a long time that complex endosomal movement in
eukaryotic cells consists of random switching from fast directed movement
powered by molecular motors interspersed with periods of slow movement
characterized by anomalous diffusion, rather than Brownian motion [4]. From
our results, we find that standard anomalous diffusion with a single constant
anomalous exponent and a single generalized diffusion coefficient for the
ensemble cannot adequately describe the transport of endosomes. Instead, a
time-space heterogeneous FBM with an anomalous exponent that switches between
super-diffusive motion $\alpha>1$ and sub-diffusive motility $\alpha<1$ can do
so. We conclude this because the local dynamics of endosomes display a whole
spectrum of values (Figs. 5, 6). The anomalous exponents were found to follow
an exponential probability distribution while the generalized diffusion
coefficients are described by a power-law probability distribution. This
result is robust since it does not depend on the method used to obtain these
distributions. Both standard statistical analysis and the newly developed
neural network [6] trained to detect FBM corroborate these results.

Another characteristic of the heterogeneous diffusion is the non-Gaussian
distributions of displacements and displacement increments. Often, Laplace-
like exponentially decaying densities were observed [37], which were
accompanied by apparent Brownian diffusion and exponential probability
distributions of diffusion coefficients. In our data, we observed anomalous
diffusion and power-law probability distributions of displacements and
generalized diffusion coefficients respectively for experimental trajectories
of endosomes. Theoretically, the distribution of displacements of spatially
heterogeneous ensemble of FBM trajectories [36] with the distribution of
diffusion coefficient of a power law form $p\left(D_{\alpha}\right)\sim
D_{\alpha}^{-1-\gamma}$, have power law tails of the form
$P\left(X\right)\simeq\left|X\right|^{-1-2\gamma}$ (see Supplementary Note and
Ref. [36]). For experimental endosomal trajectories we find the power law
distribution of local generalized diffusion coefficients and estimated
$\gamma\sim 0.5\ \pm 0.1$ (Fig. 5b). This gives $1+2\ \gamma\sim 2$ for the
exponent of the power-law tail of power-law distributions which agrees well
with the exponents found in the experimental PDFs (Figs. 1, 2).

What are the factors promoting the appearance of the power-law distributions
of the generalized diffusion coefficients? There are 3 possible reasons. (i)
Power-law distribution of diffusion coefficients may indicate that endosomal
networks exhibit scale free properties [63], many small endosomes, and few
large endosomes. Indeed, it was found in single-molecule experiments within
the cell that fluctuations of molecular size generate a heterogeneous
diffusion process with non-Gaussian distributions [26, 40, 53]. Hence
differences in their diameters generate distinct diffusive properties. (ii)
Power-law distribution of local diffusivities of endosomes could be related to
a power-law probability distribution of the empirical velocities of motor
proteins which drive endosomes. This power law distribution of velocities
originates from random attachments and detachments of motors to microtubules
[55]. (iii) The power law distribution of $D_{\alpha_{i}^{L}}$ could be the
non-specific interactions of endosomes with the endoplasmic reticulum and
other organelles. Indeed, non-specific interactions were recently shown to
generate heterogeneous diffusion of nanosized objects in mammalian cells [49].
It is possible that each or the interplay of all scenarios contribute to the
observed power-law probability distributions of local generalized diffusion
coefficients of endosomes. However, further experiments are required to verify
this.

Can other anomalous diffusion models describe endosomal transport? There are
three popular anomalous diffusion models for intracellular transport: the
continuous time random walk, the generalized Langevin equation and Fractional
Brownian motion. The equivalence of the ensemble and time averaged MSDs
(averaged over all trajectories) shows that the endosomal transport is ergodic
which rules out the continuous time random walk model as it is non-ergodic.
The remaining FBM and generalized Langevin equation models are both ergodic
but differ by their physical assumptions. The generalized Langevin equation
model assumes local equilibrium while the FBM model describes out of
equilibrium motion. Since molecular motors drive dynamics of endosomes out of
equilibrium, we can also rule out the generalized Langevin dynamics. This
suggests that the spatio-temporal heterogeneous FBM is the most appropriate
anomalous diffusion model to describe endosomal dynamics. The ensemble of
spatio-temporal hFBM trajectories display ergodicity e.g., the equivalence of
the ensemble and ensemble-time averaged MSDs (TMSDs averaged over all
trajectories).

In conclusion, for the first time we have shown that the anomalous
intracellular transport of endosomes is described the spatio-temporal
heterogeneous ensemble of FBM motions. We find that endosomal displacements
and increments both have strongly non-Gaussian power laws distributions.
Analysing local endosomal dynamics, we find that it is described by a spectrum
of exponentially distributed anomalous exponents and the power-law distributed
generalized diffusion coefficients. Such heterogeneity of endosomal transport
has implications for sorting and delivering molecules for different
biochemical reactions in different locations inside living cells.
Heterogeneous dynamics of endosomes in space and time has a huge impact on
these diffusion-limited reactions: it broadens the distribution of the first-
passage time to a reaction event and increases the likelihood of both short
and long trajectories hitting targets [4, 56, 57]. Both effects haves a big
impact on biochemical reaction rates, for example signalling of Rab proteins
on endosomes. The dynamic properties of other intracellular organelles, such
as mitochondria and lipid droplets, share many similarities with endosomal
transport [58, 59, 60, 61] and we anticipate their transport properties will
also be heterogeneous.

###### Acknowledgements.

The authors thank: NK and SF acknowledge financial support from EPSRC Grant
No. EP/V008641/1. DH acknowledges financial support from the Wellcome Trust
Grant No. 215189/Z/19/Z. GP is supported by the Basque Government through the
BERC 2018-2021 programs and by the Spanish Ministry of Economy and
Competitiveness MINECO through the BCAM Severo Ochoa excellence accreditation
SEV-2017-0718.

Figure S1: Probability density function of the duration of endosomal
trajectories, $T$. The dashed line represents the power law fit, $PDF(T)\sim
T^{-1.85}$.  Figure S2: Probability density function of anomalous exponents,
$\alpha_{i}$, extracted from experimental endosomal trajectories does not
depend on window size $W$ in the region $\alpha_{i}\in(0,2)$. Different curves
correspond to different window size. PDFs were obtained at $t=0.1$ s for
different window size $W$ indicated in the legend. Negative values of
$\alpha_{i}$ and $\alpha_{i}>2$ which are artefacts of fitting of TMSDs
gradually diminish with increased window size $W$.

Figure S3: Probability distributions of local anomalous exponents
$\alpha^{L}_{i}$ and local generalized diffusion coefficient
$D_{\alpha^{L}_{i}}$ estimated from local TMSDs (LTMSDs) are robust with
respect to changing window size $W$. Figures (A), (B) show PDFs for window
size $W=0.6$ s and (C), (D) for window size $W=0.2$ s.

Figure S4: The average local anomalous exponent $\left<\alpha^{L}_{i}\right>$
(A) stays constant for $t<8$ s while the average local generalized diffusion
coefficient $\left<D_{\alpha^{L}_{i}}\right>$ (B) decrease with time. For
$t>8$ s both $\left<\alpha^{L}_{i}\right>$ and
$\left<D_{\alpha^{L}_{i}}\right>$ decrease with time. Different curves
correspond to window size $0.2$ s, $0.6$ s and $1.2$ s. The error bars in both
panels represent standard deviation.  Figure S5: Anomalous exponents $\alpha$
and generalized diffusion coefficients $D$ are positively correlated. (Top
row) Scatter plots of $\alpha$ against $D$ calculated from each trajectory
using two different methods: power-law fitting to MSDs; and neural network
estimation. (Bottom row) Kernel density estimation heat maps for the scatter
plots in the top row.

Figure S6: Normalized ensemble averaged velocity auto-correlation function
(EVACF) of experimental trajectories (A) and heterogeneous ensemble of
noiseless FBM trajectories with exponential distribution of Hurst exponents
and power law distribution of generalized diffusion coefficients (B). Values
of $\delta$ are given in the legend. A common behaviour is observed: negative
peaks at $t=\delta$ gradually disappearing with time. The FBM trajectories
were generated without measurement noise.

Figure S7: Distributions of displacements (A) and distributions of
displacements scaled by their standard deviation (B) for the superstatistical
ensemble of FBM trajectories with the exponentially distributed anomalous
exponents $\alpha$ and power-law distributed generalized diffusion
coefficients $D$. PDFs were obtained by numerical integration of Eq. (5) at
different times indicated in the legends. The distributions of displacements
scaled by their standard deviation in (B) show collapse in the central part.

## Appendix A

###  .1 Derivation of the power-law distribution of displacements and
increments for heterogeneous ensemble of FBM

For simplicity we consider the 1D case which is easily generalized to 2D. For
the FBM with a constant generalized diffusion coefficient $D$ and anomalous
exponent $\alpha$, the distribution of displacements is Gaussian [20]:

$G(x,t|D,\alpha)=\frac{1}{\sqrt{4\pi
Dt^{\alpha}}}\exp\left(-\frac{x^{2}}{4Dt^{\alpha}}\right).$ (1)

In a superstatistical ensemble, each trajectory is assumed to explore a local
spatial patch which is characterized by its own constant diffusivity $D$ and
anomalous exponent $\alpha$. The values of $D$ are drawn from the PDF $f(D)$
and values of $\alpha$ from the PDF $g(\alpha)$. The resulting distribution of
displacements of the heterogeneous ensemble of FBM trajectories is given by
the convolution:

$P(x,t)=\int_{0}^{2}\int_{0}^{\infty}G(x,t|D,\alpha)f(D)g(\alpha)dDd\alpha.$
(2)

We consider the PDFs of $\alpha$ and $D$ which were suggested by the
experimental data. Specifically, we choose the power law PDF of $D$:

$p(D)=\frac{\gamma D^{-1-\gamma}}{D_{-}^{-\gamma}-D_{+}^{-\gamma}},$ (3)

and the exponential PDF of $\alpha$:

$h(\alpha)=\frac{e^{-\alpha/\left<\alpha\right>}}{\left<\alpha\right>(1-e^{-2/\left<\alpha\right>})}.$
(4)

From experimental trajectories we estimated $\gamma\sim 0.5$ and
$\left<\alpha\right>\sim 0.53$ (Fig. 5). The distribution of displacements Eq.
2 then reads:

$P(x,t)=\frac{\gamma
4^{\gamma}}{\sqrt{\pi}(D_{-}^{-\gamma}-D_{+}^{-\gamma})(1-e^{-2/\left<\alpha\right>})}\int_{0}^{2}e^{-\alpha\left(\frac{1}{\left<\alpha\right>}-\frac{\gamma\ln{t}}{2}\right)}\left[\Gamma\left(\gamma+\frac{1}{2},\frac{x^{2}}{4D_{+}t^{\gamma}}\right)-\Gamma\left(\gamma+\frac{1}{2},\frac{x^{2}}{4D_{-}t^{\gamma}}\right)\right]d\alpha.$
(5)

Fig. S7 (A) illustrates PDFs of $x$ at different time which were obtained by
numerical integration of Eq. 5. The PDFs of the scaling variable $X/\sigma$ in
Fig. S7 (B) show the data collapse in the central part.

Considering the Fourier-Laplace transforms and expanding integrals in Eq. (2)
at long time and space regime, the power law behaviour of the distribution of
displacements is [42]:

$P(x,t)\simeq|x|^{-1-2\gamma}.$ (6)

Since for the distribution of increments for standard FBM is also Gaussian
[20], a similar derivation holds for the power law distribution of increments
for the superstatistical ensemble of FBM trajectories.

###  .2 Derivation of the ensemble averaged mean squared displacement (EMSD)
for superstatistical ensemble of FBM

The EMSD is defined as:

$\left<x^{2}(t)\right>=\int_{-\infty}^{\infty}x^{2}P(x,t)dx.$ (7)

Using the expression for $P(x,t)$ and PDFs $f(D)$, $h(\alpha)$, we obtain the
following expression:

$\left<x^{2}(t)\right>=\frac{2\left<D\right>}{1-\left<\alpha\right>\ln{t}}\frac{1-e^{-2/\left<\alpha\right>}t^{2}}{1-e^{-2/\left<\alpha\right>}},$
(8)

where

$\left<D\right>=\frac{\gamma}{1-\gamma}\frac{(D_{+}^{1-\gamma}-D_{-}^{1-\gamma})}{(D_{-}^{-\gamma}-D_{+}^{-\gamma})}.$
(9)

###  .3  The ensemble and the time averaged velocity auto-correlation
functions for superstatistical ensemble of FBM

The velocity auto-correlation function (ACF) averaged over an ensemble of
trajectories (EACF) is calculated as:

$C_{v}^{\delta}(t)=\left<\vec{v}(t)\vec{v}(0)\right>.$ (10)

where $\vec{v}=\frac{\vec{r}(t+\delta)-\vec{r}(t)}{\delta}$. The time averaged
ACF along a single trajectory (TACF) is defined as:

$\overline{C_{v}^{\delta}(t)}=\frac{\int_{0}^{T-t-\delta}\vec{v}(t^{\prime}+t)\vec{v}(t^{\prime})dt^{\prime}}{T-t-\delta}.$
(11)

## Appendix B: Supplementary figures

Figure S1: the probability density function of the duration of experimental
endosomal trajectories, $T$, is shown. The power law, $PDF(T)\sim T^{-1.85}$,
fits well the data.

Figure S2: The PDFs of anomalous exponents $\alpha_{i}$ extracted from
experimental endosomal trajectories are robust. The Figure shows PDFs obtained
at $t=0.1$ s for different window size $W$ indicated in the legend. In the
region $\alpha_{i}\in(0,2)$, PDFs stay the same. The negative values of
$\alpha_{i}$, which are artefacts of fitting of TMSDs, gradually diminish with
increased window size $W$.

Figure S3: The PDFs of local anomalous exponents $\bar{\alpha}$ and local
generalized diffusion coefficient $D_{\bar{\alpha}}$ estimated from local
TMSDs (LTMSDs) of experimental endosomal trajectories are robust. PDFs of
$\bar{\alpha}$ and $D_{\bar{\alpha}}$ obtained for different window size $W$
are the same.

Figure S4: The average over all trajectories local anomalous exponents
$\left<\bar{\alpha}\right>$ (panel (A)) stay constant for $t<8$ s while the
average local generalized diffusion coefficients
$\left<D_{\bar{\alpha}}\right>$ (panel (B)) decrease with time. This suggests
that longer trajectories have smaller generalized diffusion coefficients, see
Fig. 4 and Fig. 8a in the main text. For $t>8$ s both
$\left<\bar{\alpha}\right>$ and $\left<D_{\bar{\alpha}}\right>$ decrease with
time which could be due to limited statistics.

Figure S5: Anomalous exponents $\alpha$ and generalized diffusion coefficients
$D$ are positively correlated. (Top row) Scatter plots of $\alpha$ against $D$
calculated from each trajectory using two different methods: power-law fitting
to MSDs; and neural network estimation. (Bottom row) Kernel density estimation
heat maps for the scatter plots in the top row.

Figure S6: Normalized ensemble averaged velocity auto-correlation function
(EVACF) of experimental trajectories (A) and heterogeneous ensemble of FBM
trajectories with exponential distribution of Hurst exponents and power law
distribution of generalized diffusion coefficients (B). Values of $\delta$ are
given in the legend. A common behaviour is observed: negative peaks at
$t=\delta$ gradually disappearing with time. The FBM trajectories were
generated without measurement noise.

Figure S7: Distributions of displacements (A) and distributions of
displacements scaled by their standard deviation (B) for the superstatistical
ensemble of FBM trajectories with the exponentially distributed anomaous
exponents $\alpha$ and power-law distributed generalized diffusion
coefficients $D$. PDFs were obtained by numerical integration of Eq. (5) at
different times indicated in the legends. The distributions of displacements
scaled by their standard deviation in (B) show collapse in the central part.

## References

  * [1] G. Seisenberger, M.U. Ried, T. Endreß, H. Büning, M. Hallek, C. Bräuchle, Real-time singlemolecule imaging of the infection pathway of an adeno-associated virus, Science 294, 1929-1932 (2001).
  * [2] E. Nielsen, F. Severin, J.M. Backer, A.A. Hyman, M. Zerial, Rab5 regulates motility of early endosomes on microtubules, Nat. Cell. Biol. 1, 376-82 (1999).
  * [3] O.J. Driskell, A. Mironov, V.J. Allan, P.G. Woodman, Dynein is required for receptor sorting and the morphogenesis of early endosomes, Nat. Cell. Biol. 9, 113-20 (2000).
  * [4] N. Flores-Rodriguez, S.S. Rogers, D.A. Kenwright, T.A. Waigh, P.G. Woodman, V.J. Allan, Roles of dynein and dynactin in early endosome dynamics revealed using automated tracking and global analysis, PLoS ONE 6, e24479 (2011).
  * [5] A.L. Zajac, Y.E. Goldman, E.L. Holzbaur, E.M. Ostap, Local cytoskeletal and organelle interactions impact molecular-motor-driven early endosomal trafficking , Curr. Biol. 23, 1173-80 (2013).
  * [6] D. Han, N. Korabel, R. Chen, M. Johnston, A. Gavrilova, V.J. Allan, S. Fedotov, T.A. Waigh, Deciphering anomalous heterogeneous intracellular transport with neural networks, eLife 9, e52224 (2020).
  * [7] J. Rink, E. Ghigo, Y. Kalaidzidis, M. Zerial, Rab conversion as a mechanism of progression from early to late endosomes, Cell 122, 735-49 (2005).
  * [8] K. Nørregaard, R. Metzler, C.M. Ritter, K. Berg-Sørensen and L.B. Oddershede, Manipulation and motion of organelles and single molecules in living cells, Chem. Rev. 117, 4342-75 (2017).
  * [9] K. Chen, B. Wang, and S. Granick, Memoriless self-reinforcing directionality in endosomal active transport within living cells, Nat. Mater. 14, 589 (2015).
  * [10] R.P. Kulkarni, K. Castelino, A. Majumdar, S.E. Fraser, Intracellular transport dynamics of endosomes containing DNA polyplexes along the microtubule network, Biophys. J. 90, L42 (2006).
  * [11] R. Metzler and J. Klafter, The random walk’s guide to anomalous diffusion: A fractional dynamics approach, Phys. Rep. 339, 1-77 (2000).
  * [12] R. Klages, G. Radons and I.M. Sokolov (eds.), Anomalous Transport: Foundations and Applications (Wiley VCH - Verlag, Weinheim, 2004).
  * [13] F. Höfling and T. Franosch, Anomalous transport in the crowded world of biological cells, Rep. Prog. Phys. 76, 046602 (2013).
  * [14] Y. Meroz and I.M. Sokolov, A toolbox for determining subdiffusive mechanisms, Phys. Rep., 573, 1-29 (2015).
  * [15] E. Barkai, Y. Garini and R. Metzler, Strange kinetics of single molecules in living cells, Physics Today 65, 8, 29 (2012).
  * [16] T.A. Waigh, The Physics of Living Processes: A Mesoscopic Approach (John Wiley & Sons, 2014).
  * [17] J. Klafter, I.M. Sokolov, First Steps in Random Walks: From Tools to Applications. Oxford University Press (Oxford University Press, 2011).
  * [18] I. M. Sokolov, Models of anomalous diffusion in crowded environments, Soft Matter 8, 9043-9052 (2012).
  * [19] P.C. Bressloff and J.M. Newby, Stochastic models of intracellular transport, Rev. Mod. Phys. 85, 135-196 (2013).
  * [20] R. Metzler, J.-H. Jeon, A.G. Cherstvy and E. Barkai, Anomalous diffusion models and their properties: non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking, Phys. Chem. Chem. Phys., 16 24128 (2014).
  * [21] V. Briane, M. Vimond, C. Kervrann, An overview of diffusion models for intracellular dynamics analysis, Brief. Bioinform. 21, 1136-1150 (2020).
  * [22] M. J. Saxton, Single-particle tracking: the distribution of diffusion coefficients, Biophys. J., 72, 1744-1753 (1997).
  * [23] B. Wang, S.M. Anthony, S.C. Bae and S. Granick, Anomalous yet Brownian, PNAS 106, 15160 (2009).
  * [24] B. Wang, J. Kuo, S.C. Bae, S. Granick, When Brownian diffusion is not Gaussian, Nat. Mater. 11, 481-485 (2012).
  * [25] D. Arcizet, B. Meier, E. Sackmann, J.O. Radler and D. Heinrich, Temporal Analysis of Active and Passive Transport in Living Cells, Phys. Rev. Lett. 101, 248103 (2008).
  * [26] T.J. Lampo, S. Stylianidou, M.P. Backlund, P. A. Wiggins and A.J. Spakowitz, Cytoplasmic RNA-Protein Particles Exhibit Non-Gaussian Subdiffusive Behavior, Biophys. J. 112, 532-542 (2017).
  * [27] A. Sabri, X. Xu, D. Krapf and M. Weiss, Elucidating the Origin of Heterogeneous Anomalous Diffusion in the cytoplasm of mammalian cells, Phys. Rev. Lett 125, 058101 (2020).
  * [28] Q. Ba, G. Raghavan, K. Kiselyov and G. Yang, Whole-cell scale dynamic organization of lysosomes revealed by spatial statistical analysis, Cell Reports, 23, 3591-3606 (2018).
  * [29] C. Manzo, J.A. Torreno-Pina, P. Massignan, G.J. Lapeyre, Jr., M. Lewenstein and M.F. Garcia-Parajo, Weak Ergodicity Breaking of Receptor Motion in Living Cells Stemming from Random Diffusivity, Phys. Rev. X, 5, 011021 (2015).
  * [30] C. Beck and E. D. B. Cohen, Superstatistics, Physica A, 322, 267 (2003).
  * [31] C. Beck, Superstatistical Brownian motion, Prog. Theor. Phys. Suppl. 162 29 (2006).
  * [32] A. Mackała and M. Magdziarz, Statistical analysis of superstatistical fractional Brownian motion and applications, Phys. Rev. E, 99, 012143 (2019).
  * [33] D. Molina-García, T.M. Pham, P. Paradisi, C. Manzo and G. Pagnini, Fractional kinetics emerging from ergodicity breaking in random media, Phys. Rev. E 94, 052147 (2016).
  * [34] M.V. Chubynsky and G.W. Slater, Diffusing diffusivity: a model for anomalous yet Brownian diffusion, Phys. Rev. Lett. 113, 098302 (2014).
  * [35] V. Sposini, A.V. Chechkin, F. Seno, G. Pagnini and R. Metzler, Random diffusivity from stochastic equations: comparison of two models for Brownian yet non-Gaussian diffusion, New J. Phys. 20, 043044 (2018).
  * [36] A.V. Chechkin, F. Seno, R. Metzler and I.M. Sokolov, Brownian yet Non-Gaussian Diffusion: From Superstatistics to Subordination of Diffusing Diffusivities, Phys. Rev. X 7, 021002 (2017).
  * [37] R. Metzler, Superstatistics and Non-Gaussian, Eur. Phys.; J. Special Topics 229, 711-728 (2020).
  * [38] A.G. Cherstvy and R. Metzler, Nonergodicity, fluctuations, and criticality in heterogeneous diffusion processes, Phys. Rev. E, 90, 012134 (2014).
  * [39] A. J. Spakowitz, Transient Anomalous Diffusion in a Heterogeneous Environment, Front. Phys. 7, 119 (2019).
  * [40] A.A. Sadoon and Y. Wang, Anomalous, non-Gaussian, viscoelastic, and age-dependent dynamics of histonelike nucleoid-structuring proteins in live Escherichia coli, Phys. Rev. E 98, 042411 (2018).
  * [41] Y. Itto and C. Beck, Superstatistical modelling of protein diffusion dynamics in bacteria, J. R. Soc. Interface 18, 20200927 (2021).
  * [42] W. Wang, A.G. Cherstvy, A.V. Chechkin, S. Thapa, F. Seno, X. Liu and R. Metzler, Fractional Brownian motion with random diffusivity: emerging residual nonergodicity below the correlation time, J. Phys. A: Math. Theor. 53, 474001 (2020).
  * [43] N. Korabel and E. Barkai, Paradoxes of subdiffusive infiltration in disordered systems, Phys. Rev. Let. 104, 170603 (2010).
  * [44] S. Fedotov and N. Korabel, Self-organized anomalous aggregation of particles performing nonlinear and non-Markovian random walks, Phys. Rev. E 92, 062127 (2015).
  * [45] S. Fedotov and D. Han, Asymptotic behavior of the solution of the space dependent variable order fractional diffusion equation: ultraslow anomalous aggregation, Phys. Rev. Lett. 123, 050602 (2019).
  * [46] T. Sandev, A.V. Chechkin, N. Korabel, H. Kantz, I.M. Sokolov and R. Metzler, Distributed-order diffusion equations and multifractality: Models and solutions, Phys. Rev. E 92, 042117 (2015).
  * [47] A. Ayache and M. S. Taqqu, Multifractional processes with random exponent, Publicacions Matemàtiques 49, 459-486 (2005).
  * [48] J.M. Newby, A.M. Schaefer, P.T. Lee, M.G. Forest and S. K. Lai, Convolutional neural networks automate detection for tracking of submicron-scale particles in 2D and 3D, PNAS 115, 9026-9031 (2018).
  * [49] F. Etoc, et al. Non-specific interactions govern cytosolic diffusion of nanosized objects in mammalian cells, Nat. Mater. 17, 740-746 (2018).
  * [50] S.C. Weber, M.A. Thompson, W.E. Moerner, A.J. Spakowitz and J.A. Theriot, Analytical Tools To Distinguish the Effects of Localization Error, Confinement, and Medium Elasticity on the Velocity Autocorrelation Function, Biophys. J. 102, 2443-2450 (2012).
  * [51] A. Caspi, R. Granek and M. Elbaum, Enhanced diffusion in active intracellular transport, Phys. Rev. Lett. 85, 5655 (2000).
  * [52] S. Thapa, N. Lukat, C. Selhuber-Unkel, A.G. Cherstvy and R. Metzler, Transient superdiffusion of polydisperse vacuoles in highly motile amoeboid cells, J. Chem. Phys. 150, 144901 (2019).
  * [53] W. He, H. Song, Y. Su, L. Geng, B.J. Ackerson, H.B.Peng and P. Tong, Dynamic heterogeneity and non-Gaussian statistics for acetylcholine receptors on live cell membrane, Nat. Commun. 7, 11701 (2016).
  * [54] C. Collinet, M. Stoter, C.R. Bradshaw, N. Samusik, J.C. Rink, D. Kenski, B. Habermann, F. Buchholz, R. Henschel, M.S. Mueller, W.E. Nagel, E. Fava, Y. Kalaidzidis and M. Zerial, Systems survey of endocytosis by multiparametric image analysis, Nature 464, 243-249 (2010).
  * [55] J. Hughes, S. Shastry, W.O. Hancock and J. Fricks, Estimating Velocity for Processive Motor Proteins with Random Detachment, J Agric Biol Environ Stat 18, 204-217 (2013).
  * [56] Y. Lanoiselée, N. Moutal and D.S. Grebenkov, Diffusion-limited reactions in dynamic heterogeneous media, Nature Commun. 9, 4398 (2018).
  * [57] A. Godec and R. Metzler, First passage time distribution in heterogeneity controlled kinetics: going beyond the mean first passage time, Sci. Rep. 6, 20349 (2016).
  * [58] N. Korabel, T.A. Waigh, S. Fedotov and V.J. Allan, Non-Markovian intracellular transport with sub-diffusion and run-length dependent detachment rate, PLoS ONE 13, e0207436 (2018).
  * [59] S. Fedotov, N. Korabel, T.A. Waigh , D. Han and V.J. Allan, Memory effects and Lévy walk dynamics in intracellular transport of cargoes, Phys. Rev. E 98, 042136 (2018).
  * [60] D.A. Kenwright, A.W. Harrison, T.A. Waigh, P.G. Woodman and V.J. Allan, First-passage-probability analysis of active transport in live cells, Phys Rev E 86, 031910 (2012).
  * [61] S. Detmer and D. Chan, Functions and dysfunctions of mitochondrial dynamics, Nat Rev Mol Cell Biol 8, 870-879 (2007).
  * [62] M. Götz, K.F. Hodeck, P. Witzel, A. Nandi, B. Lindner and D. Heinrich, Probing cytoskeleton dynamics by intracellular particle transport analysis, Eur. Phys. J. Special Topics 224, 1169-1183 (2015).
  * [63] L. Foret, J.E. Dawson, R. Villaseñor, C. Collinet, A. Deutsch, L. Brusch, M. Zerial, Y. Kalaidzidis, F. Jülicher, A general theoretical framework to infer endosomal network dynamics from quantitative image analysis, Curr Biol. 22, 1381-90 (2012).


* [11] Martin C. Cooper, Arnaud Lequen, and Frédéric Maris.  Isomorphisms between STRIPS problems and sub-problems.  In Christine Solnon, editor, CP 2022, volume 235 of LIPIcs, pages 13:1–13:16. Schloss Dagstuhl - Leibniz-Zentrum für Informatik, 2022.
  * [12] Bence Cserna, William J. Doyle, Jordan S. Ramsdell, and Wheeler Ruml.  Avoiding dead ends in real-time heuristic search.  In Sheila A. McIlraith and Kilian Q. Weinberger, editors, AAAI, pages 1306–1313, 2018.
  * [13] Adnan Darwiche and Pierre Marquis.  A knowledge compilation map.  Journal of Artificial Intelligence Research, 17:229–264, 2002.
  * [14] Stefan Edelkamp.  Planning with pattern databases.  In ECP, 01 2001.
  * [15] Salomé Eriksson, Gabriele Röger, and Malte Helmert.  Unsolvability certificates for classical planning.  In Laura Barbulescu, Jeremy Frank, Mausam, and Stephen F. Smith, editors, ICAPS, pages 88–97, 2017.
  * [16] Salomé Eriksson, Gabriele Röger, and Malte Helmert.  A proof system for unsolvable planning tasks.  In Mathijs de Weerdt, Sven Koenig, Gabriele Röger, and Matthijs T. J. Spaan, editors, ICAPS, pages 65–73, 2018.
  * [17] Richard Fikes, Peter E. Hart, and Nils J. Nilsson.  Learning and executing generalized robot plans.  Artificial Intelligence, 3(1-3):251–288, 1972.
  * [18] Richard Fikes and Nils J Nilsson.  STRIPS: A new approach to the application of theorem proving to problem solving.  Artificial Intelligence, 2(3/4):189–208, 1971.
  * [19] Maria Fox and Derek Long.  The detection and exploitation of symmetry in planning problems.  In Thomas Dean, editor, IJCAI, pages 956–961, 1999.
  * [20] M. R. Garey and David S. Johnson.  Computers and Intractability: A Guide to the Theory of NP-Completeness.  W. H. Freeman, 1979.
  * [21] Hector Geffner and Blai Bonet.  A Concise Introduction to Models and Methods for Automated Planning.  Morgan & Claypool Publishers, 2013.
  * [22] Moritz Göbelbecker, Thomas Keller, Patrick Eyerich, Michael Brenner, and Bernhard Nebel.  Coming up with good excuses: What to do when no plan can be found.  In Ronen I Brafman, Hector Geffner, Jörg Hoffmann, and Henry A Kautz, editors, ICAPS, pages 81–88, 2010.
  * [23] Malte Helmert.  Complexity results for standard benchmark domains in planning.  Artificial Intelligence, 143(2):219–262, 2003.
  * [24] Malte Helmert.  The fast downward planning system.  Journal of Artificial Intelligence Research, 26:191–246, 2006.
  * [25] Jörg Hoffmann, Peter Kissmann, and Alvaro Torralba.  “distance”? who cares? tailoring merge-and-shrink heuristics to detect unsolvability.  In Torsten Schaub, Gerhard Friedrich, and Barry O’Sullivan, editors, ECAI 2014, volume 263 of Frontiers in Artificial Intelligence and Applications, pages 441–446. IOS Press, 2014.
  * [26] Rostislav Horčík and Daniel Fišer.  Endomorphisms of classical planning tasks.  In AAAI, pages 11835–11843, 2021.
  * [27] Henry A Kautz and Bart Selman.  Planning as satisfiability.  In Bernd Neumann, editor, ECAI, pages 359–363, 1992.
  * [28] Richard E. Korf.  Macro-operators: A weak method for learning.  Artificial Intelligence, 26(1):35–77, 1985.
  * [29] Richard Earl Korf.  Learning to solve problems by searching for macro-operators.  PhD thesis, Carnegie Mellon University, USA, 1983.  AAI8425820.
  * [30] Songtuan Lin and Pascal Bercher.  Change the world - how hard can that be? On the computational complexity of fixing planning models.  In Zhi-Hua Zhou, editor, IJCAI, pages 4152–4159, 8 2021.
  * [31] Nir Lipovetzky, Christian J. Muise, and Hector Geffner.  Traps, invariants, and dead-ends.  In Amanda Jane Coles, Andrew Coles, Stefan Edelkamp, Daniele Magazzeni, and Scott Sanner, editors, ICAPS, pages 211–215, 2016.
  * [32] Mao Luo, Chu-Min Li, Fan Xiao, Felip Manyà, and Zhipeng Lü.  An effective learnt clause minimization approach for CDCL SAT solvers.  In Carles Sierra, editor, IJCAI 2017, pages 703–711, 2017.
  * [33] Alan K Mackworth.  Consistency in networks of relations.  Artificial intelligence, 8(1):99–118, 1977.
  * [34] Steven Minton.  Selectively generalizing plans for problem-solving.  In Aravind K. Joshi, editor, IJCAI, pages 596–599, 1985.
  * [35] Bharath Muppasani, Vishal Pallagani, Biplav Srivastava, and Forest Agostinelli.  On solving the rubik’s cube with domain-independent planners using standard representations.  arXiv, 2023.
  * [36] Nir Pochter, Aviv Zohar, and Jeffrey S. Rosenschein.  Exploiting problem symmetries in state-based planners.  In Wolfram Burgard and Dan Roth, editors, AAAI, 2011.
  * [37] Gabriele Röger, Silvan Sievers, and Michael Katz.  Symmetry-based task reduction for relaxed reachability analysis.  In ICAPS, pages 208–217, 2018.
  * [38] Tomas Rokicki and Morley Davidson.  God’s number is 26 in the quarter-turn metric.  https://web.archive.org/web/20220906131819/https://www.cube20.org/qtm/, 2022\.  Accessed: 2022-09-6.
  * [39] Francesca Rossi, Peter Van Beek, and Toby Walsh.  Handbook of Constraint Programming.  Elsevier, 2006.
  * [40] Karem A. Sakallah.  Symmetry and satisfiability.  In Armin Biere, Marijn Heule, Hans van Maaren, and Toby Walsh, editors, Handbook of Satisfiability, volume 336 of Frontiers in Artificial Intelligence and Applications, pages 509–570. IOS Press, second edition edition, 2021.
  * [41] Jendrik Seipp and Malte Helmert.  Counterexample-guided cartesian abstraction refinement for classical planning.  Journal of Artificial Intelligence Research, 62:535–577, 2018.
  * [42] Alexander Shleyfman and Peter Jonsson.  Computational complexity of computing symmetries in finite-domain planning.  Journal of Artificial Intelligence Research, 70:1183–1221, may 2021.
  * [43] Alexander Shleyfman, Michael Katz, Malte Helmert, Silvan Sievers, and Martin Wehrle.  Heuristics and symmetries in classical planning.  In Blai Bonet and Sven Koenig, editors, AAAI, pages 3371–3377, 2015.
  * [44] Silvan Sievers and Malte Helmert.  Merge-and-shrink: A compositional theory of transformations of factored transition systems.  Journal of Artificial Intelligence Research, 71:781–883, 2021.
  * [45] Silvan Sievers, Gabriele Röger, Martin Wehrle, and Michael Katz.  Structural symmetries of the lifted representation of classical planning tasks.  ICAPS Workshop on Heuristics and Search for Domain-independent Planning, 1:67–74, 2017.
  * [46] Sarath Sreedharan, Siddharth Srivastava, David E. Smith, and Subbarao Kambhampati.  Why can’t you do that HAL? explaining unsolvability of planning tasks.  In Sarit Kraus, editor, IJCAI, pages 1422–1430, 2019.
  * [47] Simon Ståhlberg, Guillem Francès, and Jendrik Seipp.  Learning generalized unsolvability heuristics for classical planning.  In Zhi-Hua Zhou, editor, IJCAI, pages 4175–4181, 2021.
  * [48] Marcel Steinmetz and Jörg Hoffmann.  State space search nogood learning: Online refinement of critical-path dead-end detectors in planning.  Artificial Intelligence, 245:1–37, 2017.
  * [49] Alvaro Torralba.  Sympa: Symbolic perimeter abstractions for proving unsolvability.  UIPC 2016 planner abstracts, pages 8–11, 2016.
  * [50] Julian R. Ullmann.  An algorithm for subgraph isomorphism.  Journal of the ACM, 23(1):31–42, 1976.
  * [51] Viktor N. Zemlyachenko, Nickolay M. Korneenko, and Regina I. Tyshkevich.  Graph isomorphism problem.  Journal of Soviet Mathematics, 29(4):1426–1481, 1985.

##  Appendix A Proofs of the relations between embeddings and projections

Lemma 3. _Suppose that $(\upsilon{},\nu{})$ is an embedding of $P^{\prime}$ in
$P$, and let $\Theta^{P^{\prime}}$ be the state-space of $P^{\prime}$. Let
$\alpha$ be the projection of $P$ over the set of fluents
$\upsilon\left(F^{\prime}\right)$ over which the embedding ranges, and let
$\Theta^{\alpha}$ be the abstract state space. Then_

  * •

_There exists a bijection_ $b$ between the states of $\Theta^{\alpha}$ and the
states of $\Theta^{P^{\prime}}$;

  * •

_For any two different states_ $s^{\alpha}_{1}\not=s^{\alpha}_{2}$ of
$\Theta^{\alpha}$, if there exists a transition from $s^{\alpha}_{1}$ to
$s^{\alpha}_{2}$, then there exists a transition from $b(s^{\alpha}_{1})$ to
$b(s^{\alpha}_{2})$ in $\Theta^{P^{\prime}}$

###### Proof.

With $P^{\prime}=\langle F^{\prime},I^{\prime},O^{\prime},G^{\prime}\rangle$,
let us denote $\Theta^{P^{\prime}}=\langle
2^{F^{\prime}},O^{\prime},T^{\prime},I^{\prime},{S^{G^{\prime}}}\rangle$ the
state-space of $P^{\prime}$, and $\Theta^{\alpha}=\langle
S^{\alpha},O,T^{\alpha},\alpha(I),\alpha(S^{G})\rangle$ the LTS associated to
the projection $\alpha$ over $\upsilon\left(F^{\prime}\right)$. We show that
any transition of $T^{\alpha}$ between two _different_ states $s_{1},s_{2}\in
S^{\alpha}\subseteq 2^{F}$ has an equivalent transition in $P^{\prime}$.

For a start, notice that $\alpha\circ\upsilon{}$ is a bijection from the set
$2^{F^{\prime}}$ of states of $\Theta^{P^{\prime}}$, to the set
$S^{\alpha}=2^{\upsilon\left(F^{\prime}\right)}$ of states of
$\Theta^{\alpha}$. Indeed, since for any $s\in 2^{F}$,
$\alpha(s)=s\cap\upsilon\left(F^{\prime}\right)$, $\alpha$ is surjective when
restricted to the domain $2^{\upsilon\left(F^{\prime}\right)}$, and thus a
bijection since $\alpha:2^{F}\rightarrow 2^{\upsilon\left(F^{\prime}\right)}$
(by definition of a projection, and using a cardinality argument). Likewise,
$\upsilon{}$ defines a bijection from $2^{F^{\prime}}$ to
$2^{\upsilon\left(F^{\prime}\right)}$ since $\upsilon{}$ is injective by
hypothesis (and thus bijective from $F^{\prime}$ to
$\upsilon\left(F^{\prime}\right)$). Let us denote
$b=\left(\alpha\circ\upsilon{}\right)^{-1}$ the bijection from $S^{\alpha}$ to
$2^{F^{\prime}}$.

Let us now show that every _non-reflexive_ transition $\langle
s^{\alpha}_{1},o,s^{\alpha}_{2}\rangle\in T^{\alpha}$ of $\Theta^{\alpha}$ has
an equivalent transition $\langle
b(s^{\alpha}_{1}),\nu\left(o\right),b(s^{\alpha}_{2})\rangle\in\Theta^{P^{\prime}}$.
Let $\langle s^{\alpha}_{1},o,s^{\alpha}_{2}\rangle\in T^{\alpha}$ be a
transition of the abstract state space. By definition of the transitions
$T^{\alpha}$ of an abstraction, there exists a transition $\langle
s_{1},o,s_{2}\rangle\in T$ in the concrete state space, such that
$\alpha(s_{1})=s^{\alpha}_{1}$ and $\alpha(s_{2})=s^{\alpha}_{2}$.

First, suppose that $o$ does not satisfy condition (18), i.e. that
$(\textsf{eff}^{+}(o)\cup\textsf{eff}^{-}(o))\cap\upsilon\left(F^{\prime}\right)=\emptyset$.
Then since $s_{2}=s_{1}[o]$, and
$\alpha(s)=s\cap\upsilon\left(F^{\prime}\right)$, we have that

$\displaystyle\alpha(s_{2})$ $\displaystyle=\alpha(s_{1}[o])$
$\displaystyle=\left((s_{1}\setminus\textsf{eff}^{-}(o))\cup\textsf{eff}^{+}(o)\right)\cap\upsilon\left(F^{\prime}\right)$
$\displaystyle=\left(\left(s_{1}\cap\upsilon\left(F^{\prime}\right)\right)\setminus\left(\textsf{eff}^{-}(o)\cap\upsilon\left(F^{\prime}\right)\right)\cup\left(\textsf{eff}^{+}(o)\cap\upsilon\left(F^{\prime}\right)\right)\right)$
$\displaystyle=s_{1}\cap\upsilon\left(F^{\prime}\right)$ (since (18) is not
satisfied) $\displaystyle=\alpha(s_{1})$

As a consequence, the transition $\langle
s^{\alpha}_{1},o,s^{\alpha}_{2}\rangle$ is reflexive.

Now suppose that condition (18) is satisfied by $o$, i.e. that
$(\textsf{eff}^{+}(o)\cup\textsf{eff}^{-}(o))\cap\upsilon\left(F^{\prime}\right)\not=\emptyset$.
Then let us show that $\langle
s^{\prime}_{1},o^{\prime},s^{\prime}_{2}\rangle$ is a transition in
$\Theta^{P^{\prime}}$, where $s^{\prime}_{1}=\upsilon{}^{-1}(s_{1})$,
$s^{\prime}_{2}=\upsilon{}^{-1}(s_{2})$, and $o^{\prime}=\nu\left(o\right)$.
For this, we have to show that $\textsf{pre}(o^{\prime})\subseteq
s^{\prime}_{1}$ and that $s^{\prime}_{2}=s^{\prime}_{1}[o^{\prime}]$.

We start by showing that $o^{\prime}$ is applicable in $s^{\prime}_{1}$. Since
$\textsf{pre}(o)\subseteq s_{1}$ and
$\upsilon\left(\textsf{pre}(\nu\left(o\right))\right)\subseteq\textsf{pre}(o)\cap\upsilon\left(F^{\prime}\right)$
(since $(\upsilon{},\nu{})$ is an embedding), we have

$\displaystyle\upsilon\left(\textsf{pre}(\nu\left(o\right))\right)\subseteq
s_{1}\cap\upsilon\left(F^{\prime}\right)$

So, by the application of $\upsilon{}^{-1}$ on both sides,

$\displaystyle\textsf{pre}(\nu\left(o\right))$
$\displaystyle\subseteq\upsilon{}^{-1}(s_{1}\cap\upsilon\left(F^{\prime}\right))$
$\displaystyle\subseteq s^{\prime}_{1}$

We now show that $s^{\prime}_{2}=s^{\prime}_{1}[o^{\prime}]$, and we follow a
similar schema as above to prove this.

$\displaystyle s_{2}$ $\displaystyle=s_{1}[o]$
$\displaystyle=(s_{1}\setminus\textsf{eff}^{-}(o))\cup\textsf{eff}^{+}(o)$
$\displaystyle\upsilon{}^{-1}(s_{2})$
$\displaystyle=\upsilon{}^{-1}\left((s_{1}\setminus\textsf{eff}^{-}(o))\cup\textsf{eff}^{+}(o)\right)$
$\displaystyle=\upsilon{}^{-1}\left(\left((s_{1}\setminus\textsf{eff}^{-}(o))\cup\textsf{eff}^{+}(o)\right)\cap\upsilon\left(F^{\prime}\right)\right)$
$\displaystyle=\upsilon{}^{-1}\left(\left(\left(\left(s_{1}\cap\upsilon\left(F^{\prime}\right)\right)\setminus\left(\textsf{eff}^{-}(o)\cap\upsilon\left(F^{\prime}\right)\right)\right)\cup\left(\textsf{eff}^{+}(o)\cap\upsilon\left(F^{\prime}\right)\right)\right)\right)$
$\displaystyle=\upsilon{}^{-1}\left(\left(\left(s_{1}\cap\upsilon\left(F^{\prime}\right)\right)\setminus\upsilon\left(\textsf{eff}^{-}(\nu\left(o\right))\right)\right)\cup\upsilon\left(\textsf{eff}^{+}(\nu\left(o\right))\right)\right)$
(using (19) and (20))
$\displaystyle=\upsilon{}^{-1}\left(\left(\left(s_{1}\cap\upsilon\left(F^{\prime}\right)\right)\right)\setminus\textsf{eff}^{-}(\nu\left(o\right))\right)\cup\textsf{eff}^{+}(\nu\left(o\right))$
$\displaystyle\text{(since }\upsilon{}^{-1}\circ\upsilon{}=Id\text{)}$
$\displaystyle=\left(s^{\prime}_{1}\setminus\textsf{eff}^{-}(\nu\left(o\right))\right)\cup\textsf{eff}^{+}(\nu\left(o\right))$
$\displaystyle s^{\prime}_{2}$ $\displaystyle=s^{\prime}_{1}[o^{\prime}]$

As a consequence, we have $\langle
s^{\prime}_{1},\nu\left(o\right),s^{\prime}_{2}\rangle\in T^{\prime}$, which
is such that
$\langle\alpha(\upsilon\left(s^{\prime}_{1}\right)),o,\alpha(\upsilon\left(s^{\prime}_{2}\right))\rangle=\langle
s^{\alpha}_{1},o,s^{\alpha}_{2}\rangle\in T^{\alpha}$.

In conclusion, we have showed that every non-reflexive transition in
$\Theta^{\alpha}$ has an equivalent transition in $\Theta^{P^{\prime}}$, with
respect to the bijection $b=\left(\alpha\circ\upsilon{}\right)^{-1}$ from the
states of $\Theta^{\alpha}$ to the states of $\Theta^{P^{\prime}}$. ∎

Proposition 6 Suppose that $(\upsilon{},\nu{})$ is an embedding of
$P^{\prime}$ in $P$, and let $\Theta^{P^{\prime}}$ be the state-space of
$P^{\prime}$. Let $\alpha$ be the projection of $P$ over the set of fluents
$\upsilon\left(F^{\prime}\right)$ over which the embedding ranges, and let
$\Theta^{\alpha}$ be the abstract state space. Then if there exists a path
from the initial state to a goal state in $\Theta^{\alpha}$, there also exists
a path from the initial state to a goal state in $\Theta^{P^{\prime}}$.

###### Proof.

Suppose that there exists a path from the initial state to a goal state in
$\Theta^{\alpha}$. Let us use the notations of Lemma 3. Lemma 3 and its proof
give a bijection $b$ between states of $\Theta^{\alpha}$ and
$\Theta^{P^{\prime}}$ such that if a path between two states
$s^{\alpha}_{1},s^{\alpha}_{2}$ of $\Theta^{\alpha}$ exists, then there exists
a path from $b(s^{\alpha}_{1})$ to $b(s^{\alpha}_{2})$ in
$\Theta^{P^{\prime}}$. This is straightforward to show by induction.

Note that the initial state $s^{\alpha}_{I}$ of $\Theta^{\alpha}$ is not
necessarily mapped by $b$ to the initial state $I^{\prime}$ of
$\Theta^{P^{\prime}}$ (this is due to condition (23), which does not enforce
equality).

Two points are to be shown. First, we show that if
$s^{\alpha}_{G}\in\alpha(S^{G})$, then $b(s^{\alpha}_{G})\in S^{G^{\prime}}$
(i.e. that $b$ maps goals states of $\Theta^{\alpha}$ to goal states of
$\Theta^{P^{\prime}}$). Then, we show that there exists a state $s^{\alpha}\in
S^{\alpha}$ such that $b(s^{\alpha})=I^{\prime}$ and there exists a path in
$\Theta^{\alpha}$ from $s^{\alpha}$ to some goal state $s^{\alpha}_{G}$.

Let $s^{\alpha}_{G}\in\alpha(S^{G})$. By definition of $\alpha$ and $S^{G}$,
we have that $G\cap\upsilon\left(F^{\prime}\right)\subseteq s^{\alpha}_{G}$.
In addition, since $\upsilon{}$ is an embedding, we have (by (22)) that
$\upsilon\left(G^{\prime}\right)\subseteq
G\cap\upsilon\left(F^{\prime}\right)$. Combining the two equations, we obtain
$\upsilon\left(G^{\prime}\right)\subseteq s^{\alpha}_{G}$. Applying $\alpha$
on both sides (which is possible since $\alpha:2^{F}\rightarrow
2^{\upsilon\left(F^{\prime}\right)}$ and
$2^{\upsilon\left(F^{\prime}\right)}\subseteq 2^{F}$), we have
$\alpha\circ\upsilon\left(G^{\prime}\right)\subseteq s^{\alpha}_{G}$ (since
the projection $\alpha$ is idempotent). By applying
$b=\left(\alpha\circ\upsilon{}\right)^{-1}$ on both sides, we have that
$G^{\prime}\subseteq b(s^{\alpha}_{G})$, which shows that $s^{\alpha}_{G}$ is
mapped to a goal state in $\Theta^{P^{\prime}}$.

Let us now show that there exists a state $s^{\alpha}\in S^{\alpha}$ such that
$b(s^{\alpha})=I^{\prime}$ and there exists a path in $\Theta^{\alpha}$ from
$s^{\alpha}$ to some goal state $s^{\alpha}_{G}$. First remark that, if there
exists a path from $\alpha(I)$ to some $s^{\alpha}_{G}\in\alpha(S^{G})$ in
$\Theta^{\alpha}$, then there exists a path from any $s^{\alpha}\in
S^{\alpha}$ such that $\alpha(I)\subseteq s^{\alpha}$ to some
$s^{\prime\alpha}_{G}\in\alpha(S^{G})$. Let us consider
$s^{\alpha}_{I^{\prime}}=\alpha(\upsilon\left(I^{\prime}\right))\in
S^{\alpha}$. We have:

$\displaystyle s^{\alpha}_{I^{\prime}}$
$\displaystyle=\alpha(\upsilon\left(I^{\prime}\right))$
$\displaystyle=\upsilon\left(I^{\prime}\right)$
$\displaystyle\text{(Definition of }\alpha\text{)}$ $\displaystyle\supseteq
I\cap\upsilon\left(F^{\prime}\right)$ (By (23)) $\displaystyle=\alpha(I)$

So there exists a path from $s^{\alpha}_{I^{\prime}}$ to a goal state
$s^{\prime\alpha}_{G}$ in $\Theta^{\alpha}$. Using the results that stems from
the generalization of Lemma 3, as presented in the beginning of this proof, we
have that there exists a path from $I^{\prime}=b(s^{\alpha}_{I^{\prime}})$ to
some goal state $b(s^{\prime\alpha}_{G})$ in $\Theta^{P^{\prime}}$, which
concludes the proof.

∎