Self-reinforcing directionality generates truncated Lévy walks without the power-law assumption

Introduction. Transport in biology is crucial on multiple scales to maintain life and deviations from normal movement are hallmarks of disease and ageing Murray (2007); *murray2001mathematical; Okubo and Levin (2013). From organisms to sub-cellular molecules, their motion is usually described by persistent random walks with finite velocities (run and tumble models) Othmer et al. (1988); Hillen (2002); Tailleur and Cates (2008); Thompson et al. (2011); Mendez et al. (2010). The preference for these velocity jump models over space jump random walks arise due to the physical constraints of organisms not instantaneously jumping in space and an inertial resistance to changes in direction. In recent years, Lévy walks Zaburdaev et al. (2015) attracted attention in modelling movement patterns of living things Reynolds (2018), from sub-cellular Chen et al. (2015); Song et al. (2018); Fedotov et al. (2018); *korabel2018non to organism Raichlen et al. (2014); Ariel et al. (2015); Huda et al. (2018) scales. Until now, Lévy walks have been mostly described by coupled continuous time random walks (CTRW) Zaburdaev et al. (2015), various fractional PDEs Compte and Metzler (1997); Sokolov and Metzler (2003); Meerschaert et al. (2002); Baeumer et al. (2005); Becker-Kern et al. (2004); Uchaikin and Sibatov (2011); Magdziarz et al. (2015) and integro-differential equations Fedotov (2016); Taylor-King et al. (2016). These approaches require power-law distributed running times with divergent first and second moments as an ab initio assumption. However, in many cases this assumption is difficult to justify, leading to ongoing discussions about the origin of power-law distributed runs, the Levy walk observed in nature Jansen et al. (2012); *wosniack2017evolutionary; Reynolds (2018) and Lévy foraging hypothesis Viswanathan et al. (1999); Tejedor et al. (2012).

Experiments exhibiting Lévy-like motion cannot be modeled with pure power-law jump or running time distributions due to finite limits in physical systems Reynolds (2012). So theoretically, power-law distributions are often truncated Mantegna and Stanley (1994) or exponentially tempered Cartea and del Castillo-Negrete (2007); *sabzikar2015tempered. Specifically for active cargo transport in cells, it was discovered that the motion is composed of multiple short runs that self-organize into longer, truncated power-law distributed, uni-directional flights (truncated Lévy walks) Chen et al. (2015). Explanations for this phenomenon have been attempted in terms of self-reinforcing directionality generated by co-operative motor protein transport Chen et al. (2015). Yet, the question remains, can a persistent random walk model generate superdiffusion without power-law distributed runs through the self-organization of exponentially distributed runs?

In this Letter, we propose a new, spatially and temporally inhomogeneous persistent random walk model that generates exponentially truncated Lévy walks through self-reinforced directionality. This exponential truncation is not directly introduced but rather arises naturally from the same microscopic mechanism that generates the power-law distribution of flights. Our model differs from the traditional Lévy walk model because superdiffusion is generated by self-reinforcing directionality rather than the assumption of power-law flight distribution with infinite second moment. Simulated densities of uni-directional flight lengths from this model show excellent agreement with truncated power-law distributions observed in experiments Chen et al. (2015). It is worth noting that there are a few examples where the power-law assumption has not been used as a starting point: superdiffusion of ultracold atoms Kessler and Barkai (2012); *barkai2014area and a random walk driven by an ergodic Markov process with switching Fedotov et al. (2007).

Self-reinforcing directionality. Consider a particle moving to the right and left in 1D with constant speed $\nu$ and exponentially distributed running time with rate $\lambda$. In the standard alternating case, the particle would change directions with probability $1$. To consider the instance when there is a probability that the particle continues in the same direction as the previous run, we introduce the transition probability matrix Cox and Miller (1965):

where $q_{+}$ is the conditional transition probability that the particle will continue in the positive direction given it was moving in the positive direction before. Similarly, $q_{-}$ corresponds to the particle moving in the negative direction. The standard alternating case corresponds to $q_{+}=q_{-}=0$.

In order to model self-reinforcing directionality using the matrix $\boldsymbol{Q}$, we introduce relative times, $t^{+}/t$ and $t^{-}/t$, that the particle moves in the positive and negative direction during time $t$. The key point is to define the probabilities in (1) as

Here we introduce the persistence probability, $w$, which parameterizes the extent that changes of direction are affected by relative times. If $w=1/2$ then the transition probabilities in both directions are the same: $q_{+}=q_{-}=1/2$. If $w>1/2$ then the matrix $\boldsymbol{Q}$ has repetition compulsion properties: the longer a particle spends moving in the one direction, the greater the probability to maintain directionality. Figure 1 illustrates self-organization of individual runs into long uni-directional flights.

Let us show that the conditional transition probabilities (2) can be expressed in terms of particle position $x=x_{0}+\nu\left[t^{+}-t^{-}\right]$. Since $t=t^{+}+t^{-}$, the relative times can be written as

Then substituting (3) into (2), we find that the probabilities $q_{+}$ and $q_{-}$ depend on $x$ and $t$,

Since our random walker is moving with finite speed $\nu$, the relation $(x-x_{0})/\nu t\leq 1$ must hold. Once again, the formula (4) shows self-reinforcing directionality since if $\alpha>0$ ($w>1/2$) then the transition probability $q_{+}(x,t)$ is an increasing function of particle position $x$ and the opposite for $q_{-}(x,t)$.

The equations for probability density functions (PDFs) of particles moving right ($+$) and left ($-$), $p_{\pm}(x,t)$, can be written in terms of transition matrix $\boldsymbol{Q}$ as

Rearranging and taking the limit $\Delta t\rightarrow 0$, we can write the equations for $p_{\pm}(x,t)$ as

Note these equations (6) can be rewritten in terms of space and time dependent switching rates

These switching rates will be used to show how the exponentially truncated power-law distribution of flights arise from exponentially distributed runs. Using standard methods Mendez et al. (2010); Schnitzer (1993); Thompson et al. (2011) and equation (4), we can obtain the partial differential equations (PDEs) for total density $p(x,t)=p_{+}+p_{-}$ and flux $J(x,t)=\nu p_{+}-\nu p_{-}$: $\partial p/\partial t=-\partial J/\partial x$ and $\partial J/\partial t=-\nu^{2}\partial p/\partial x-\lambda(J-\alpha xp/t)$ (See Supp. Mat.). The initial conditions are

where $u\in[0,1]$ is the probability that the intial velocity is $\nu$. Finally, we can find the hyperbolic PDE with a non-homogeneous in space and time advection term

The advection term of Eq. (9) is unconventional because it depends on the initial position $x_{0}$. Furthermore, if the initial conditions are symmetric, $u=1/2$ (see (8)), then the average drift is zero. Clearly, (9) is a modification of the classical telegraph or Cattaneo equation Kac (1974); Schnitzer (1993); Mendez et al. (2010), with a time and space dependent advection term. In what follows, we will show that this additional term generates superdiffusion. In fact, a generalized Cattaneo equation generating superdiffusion has been formulated using the Riemann-Liouville fractional derivative Compte and Metzler (1997). The advantage of Eq. (9) over fractional PDEs is that it is far simpler and does not require integral operators in time. To the authors’ knowledge, the hyperbolic PDE (9) is the first formulation of truncated Lévy walks without integral operators Fedotov (2016). In the diffusive limit, when $\lambda\rightarrow\infty$ and $\nu\rightarrow\infty$ such that $\nu^{2}/\lambda$ is a constant, (9) becomes the governing advection-diffusion equation for the continuous approximation of the elephant random walk Schütz and Trimper (2004); *paraan2006exact; *da2013non; *da2020non. Note that the system of equations (6) with transition rates (4) can also be mapped to the hyperbolic model for chemotaxis Hillen and Stevens (2000); *stevens1997aggregation with an unorthodox external stimulus $S(x,t)$ (see Supp. Mat.)

Moment analysis. Now, we show that the variance for the underlying random process, $x(t)$, exhibits superdiffusive behavior: $\text{Var}\{x(t)\}\propto t^{2\alpha}$ with $1/2<\alpha<1$. The moments of random walk position, $\mu_{n}(t)=\int_{-\infty}^{\infty}x^{n}p(x,t)dx$, can be found from (9) as

Taking the Laplace transform and solving the resulting differential equations for the first and second moments using the initial conditions, $\mu_{1}(0)=\mu_{2}(0)=0$, $d\mu_{1}(0)/dt=\nu(2u-1)$, $d\mu_{2}(0)/dt=0$ we can obtain $\hat{\mu}_{1}(s)=\nu(2u-1)s^{-1-\alpha}(s+\lambda)^{\alpha-1}$ and $\hat{\mu}_{2}(s)=\frac{2\nu^{2}}{\lambda(2\alpha-1)}\left[s^{-1-2\alpha}(s+\lambda)^{2\alpha-1}-s^{-2}\right]$ (for details see Supp. Mat.). The inverse Laplace transform gives $\mu_{1}(t)=\nu(2u-1)t\prescript{}{1}{F}_{1}(1-\alpha,2,-\lambda t)$ and $\mu_{2}(t)=2\nu^{2}t\left[\prescript{}{1}{F}_{1}(1-2\alpha,2,-\lambda t)-1\right]/(\lambda(2\alpha-1))$, where $\prescript{}{1}{F}_{1}(a,b,z)$ is the Kummer confluent hypergeometric function. In the long time limit ($s\rightarrow 0$), we obtain

and

Clearly, the random walk exhibits superdiffusive behavior for $1/2<\alpha<1$. The variance, $\text{Var}\left\{x(t)\right\}=\mu_{2}(t)-\mu_{1}(t)^{2}$, is

where $A=\frac{2}{(2\alpha-1)\Gamma(2\alpha+1)}-\frac{(2u-1)^{2}}{\Gamma^{2}(\alpha+1)}$. Results of Monte Carlo simulations are in perfect agreement with the analytical solution of the second moment (see Supp. Mat.).

Truncated Lévy walk. In what follows, we show analytically how self-reinforcing directionality generates a truncated Lévy walk with exponentially tempered power-law distributed flights. Let us find the conditional distribution of flights given the particle changes direction at the position $x_{*}$ and time $t_{*}$ (illustrated in Fig. 1 where the particle changes velocity from $-\nu$ to $\nu$). First, we calculate the conditional survival function, $\Psi_{\pm}(\tau|x_{*},t_{*})$, for flights moving in the positive ($+$) and negative ($-$) direction. This function gives us the probability that the random duration of a flight, $T$, is greater than $\tau$. The rates, $\lambda_{\pm}$ (7) define the conditional switching rate along the particle trajectory starting at position $x_{*}$ and time $t_{*}$ as $\lambda_{\pm}(\tau|x_{*},t_{*})=\lambda(1-q_{\pm}(x_{*}\pm\nu\tau,t_{*}+\tau))$. If we rearrange using $\alpha=2w-1$, then

where $\gamma_{\pm}=(w-\frac{1}{2})\lambda\left(t_{*}\mp\frac{x_{*}}{\nu}\right)$. Using this, we can find the conditional survival function by $\Psi_{\pm}(\tau|x_{*},t_{*})=\exp\left(-\int_{0}^{\tau}\lambda_{\pm}(s|x_{*},t_{*})ds\right)$. Then, we can see that the constant term in (14) produces an exponential tempering factor and the term inversely proportional to the running time, $\tau$, generates a power-law with exponent $\gamma_{\pm}$. Explicitly, we can write

The conditional distribution of flight lengths $l=\nu\tau$ is given by: $F_{\pm}(l|x_{*},t_{*})=1-\Psi_{\pm}(l/\nu|x_{*},t_{*})$. Each uni-directional flight is drawn from a truncated power-law distribution dependent on the position $x_{*}$ and time $t_{*}$, where the particle changes direction. Clearly, this demonstrates the difference between our new model and the traditional Lévy walk model where the power-law distribution of flights is an a priori requirement and not dependent on $x_{*}$ and $t_{*}$. Furthermore, our new model is able to generate superdiffusion (12) without an infinite second moment of flight times (15) through self-reinforcing directionality, unlike traditional Lévy walk models.

In addition, exponential truncation (tempering) is usually introduced by simply multiplying the power-law jump or waiting time densities by an exponential factor with an additional parameter, leading to tempered fractional calculus Cartea and del Castillo-Negrete (2007); *sabzikar2015tempered. The advantage of our model is that both exponential tempering and power-law flight distribution are generated through a single microscopic mechanism involving self-reinforced directionality and the subsequent analysis of (9) is more convenient than tempered fractional calculus. An interesting feature of our model is that despite the exponential truncation of flights, the variance (13) still exhibits superdiffusive behavior for $1/2<\alpha<1$ ($3/4<w<1$). Numerical simulations confirm exponentially truncated power-law conditional survival function (15) for flight times (see Supp. Mat.).

Finally, we can obtain numerically the PDF, $p(l)$, of flight length, $l$, without the condition that flights begin at $x_{*}$ and $t_{*}$. In the case of strong directionality as $w\rightarrow 1$, $p(l)$ should approach a pure power-law since the truncation length $\nu/\lambda(1-w)\rightarrow\infty$. Figure 2 confirms this and shows that, when $w=0.999$, a power-law density is recovered for more than two orders of magnitude. Next, we provide evidence that the truncated power-law PDF shows excellent correspondence to published data on uni-directional endosome flights.

Experimental data and numerical simulations. Truncated Lévy walk behavior is observed experimentally in intracellular transport Chen et al. (2015), which arises from the self-organization of exponentially distributed runs, $x_{i}$, into uni-directional flights, $x_{f}$ Chen et al. (2015). Here, we switch from $l$ to $x_{f}$ to avoid confusion between theoretical and experimental flights. Until now, there had been no governing PDE like (9) and an underlying persistent random walk model to describe this phenomenon. Experimentally, the authors Chen et al. (2015) report power-law tails in the flight-length density, $p(x_{f})\propto x_{f}^{-2}$. Figure 3 demonstrates that numerical simulations of our model are able to generate the power-law tails for flight length density and emulate the experimental data on the whole $x_{f}$ scale using reasonable parameters. Furthermore, the parameters of our new model, such as persistence probability $w$ or $\alpha$, rate $\lambda$ and speed $\nu$ can be easily found by comparing the exact analytical formula of second moment with experimental mean squared displacements.

Bimodal densities and transition from diffusive to superdiffusive regime. Surprisingly, our truncated Lévy walk model in the long time limit leads to bimodal densities in the superdiffusive regime ($1/2<\alpha<1$). This phenomenon does not exist for classical superdiffusive Lévy walks. Bimodal densities (Lamperti distributions) only appear in the ballistic regime for Lévy walks with a divergent first moment for running times Uchaikin and Sibatov (2009); Klafter and Sokolov (2011); Uchaikin and Sibatov (2011); Zaburdaev et al. (2015); Froemberg et al. (2015). In the superdiffusive case, the density for Lévy walks is Gaussian in the central part with power-law tails Klafter and Sokolov (2011). This is completely different to densities for (9) (see Fig. 4). Density peaks for $p(x,t)$ in (9) occur at $|x|<\nu t$. Figure 4 shows the emergence of a bimodal density for $t=10^{2}$ and $w>w_{c}$, where $w_{c}=3/4$. The density $p(x,t)$ exhibits two distinct long time behaviors: it is Gaussian for $\alpha<1/2$ ($w<w_{c}$) and bimodal for $\alpha>1/2$ ($w>w_{c}$). The form of the PDF also shows a non-equilibrium phase transition (see Supp. Mat.). Note that similar bimodal densities are observed in velocity random walks with interacting particles Lutscher et al. (2002); *fetecau2010investigation; *carrillo2014non; Fedotov and Korabel (2017).

For $-1<\alpha<1/2$ ($0<w<3/4$), the variance (13) corresponds to the diffusive regime with the effective diffusion coefficient $D=\nu^{2}/\lambda(1-2\alpha)$. For alternating velocity random walks, the conventional diffusion coefficient would be $D_{0}=\nu^{2}/\lambda$. As $\alpha\rightarrow 1/2$, the effective diffusion coefficient tends to infinity and this indicates the transition from the diffusive to superdiffusive regime. For the long time limit of (9) in the diffusive regime when $\partial^{2}p/\partial t^{2}$ becomes negligible, the solution is Gaussian: $p(x,t)=(4\pi Dt)^{-1/2}\exp\left(-(x-\mu_{\alpha}t^{\alpha})^{2}/4Dt\right),$ where $\mu_{\alpha}=\nu(2u-1)\lambda^{\alpha-1}/\Gamma(\alpha+1)$. Figure 4 shows the excellent agreement between numerical simulations and the Gaussian solution for $w=0.6$.

Self-reinforced directionality in endosome movement. Now, the question remains: what is the underlying biological mechanism explaining self-reinforced directionality in intracellular transport? To answer this question, we suggest a simple illustration of the possible origin for self-reinforced directionality. Let us consider endosomal motility, which is governed by the adaptor complexes on its membranes, the most notable being Rab GTPases Stenmark (2009). These adaptors facilitate attachment to dynein and kinesin motors Urnavicius et al. (2018) and, therefore, dictate the positioning and motility of endosomes in the cell Stenmark (2009); Gindhart and Weber (2009). To simplify this vastly complex process, consider that the endosome contains a number, $n_{-}$, of adaptor proteins that attach to kinesin leading to transport towards the cell periphery. Similarly, $n_{+}$ is the number of adaptor proteins that attach to dynein and facilitate transport towards the cell nucleus. Then, when an endosome happens to attach to a microtubule, the simplest assumption about probabilities $q_{+}$ and $q_{-}$ in (1) would be $q_{+}=n_{+}/(n_{+}+n_{-})$ and $q_{-}=n_{-}/(n_{+}+n_{-})$. However, due to the complexity of endosomal transport we can introduce a weight, $w\in[0,1]$, such that $q_{\pm}=wn_{\pm}/(n_{+}+n_{-})+(1-w)n_{\mp}/(n_{+}+n_{-}).$

From the very beginning of endocytosis until degradation, endosomes undergo a maturation process, including the association of proteins, such as Rab5 and PI(3)K Huotari and Helenius (2011). Recent work has show that effectors of Rab5 display distinct spatial densities Villaseñor et al. (2016) suggesting that $n_{-}$ and $n_{+}$ are functions of the time spent running towards, $t^{-}$, or away, $t^{+}$, from the cell center. So, we assume that $n_{\pm}=n_{\pm}^{0}+at^{\pm}$ with $a$ being some constant rate. The more an endosome moves in towards the nucleus, the more $n_{+}$ increases and vice versa. This reinforcement rule is similar to that of discrete reinforced random walks Davis (1990); Stevens and Othmer (1997). Neglecting $n_{\pm}^{0}$, this formulation is exactly what leads to the repetition compulsion property in (2), since then $n_{-}/(n_{+}+n_{-})=t^{-}/t$ and $n_{+}/(n_{+}+n_{-})=t^{+}/t$.

Summary. In this Letter, we developed a persistent random walk model with finite velocity that generates superdiffusion without the ab initio assumption of power-law distributed run times with infinite second moment. Our model shows that truncated power-law distributed flights can be generated from exponentially distributed runs through self-reinforcing directionality. A governing hyperbolic PDE (9) for particle probability density was derived along with exact solutions for the first and second moments. The theory is able to explain the experimentally observed self-organization of exponentially distributed runs into unidirectional flights leading to exponentially truncated Lévy walks Chen et al. (2015). We showed excellent agreement between the density of flight lengths from numerical simulations and in vivo cargo transport experiments. In the superdiffusive regime, numerical simulations of particle densities show bimodal densities (aggregation), which is a new phenomenon not seen in the classical linear hyperbolic or Lévy walk models. We believe that our methodology can be used to model migrating cancer cells Mierke et al. (2011); *mierke2013integrin; Huda et al. (2018), T-cell motility Harris et al. (2012), human foraging Raichlen et al. (2014), front propagation phenomena Méndez et al. (2004); *fedotov2015persistent, first passage time problems Redner (2001); Condamin et al. (2007); *angelani2014firstand viruses mobility inside cells Lagache and Holcman (2008); *greber2006superhighway; *dodding2011coupling.