Non-Markovian Persistent Random Walk Model for Intracellular Transport

: Transport of vesicles and organelles inside cells consists of constant-speed bidirectional movement along cytoskeletal ﬁlaments interspersed by periods of idling. This transport shows many features of anomalous diffusion. In this paper, we develop a non-Markovian persistent random walk model for intracellular transport that incorporates the removal rate of organelles. The model consists of two active states with different speeds and one resting state. The organelle transitions between states with switching rates that depend on the residence time the organelle spends in each state. The mesoscopic master equations that describe the average densities of intracellular transport in each of the three states are the main results of the paper. We also derive ordinary differential equations for the dynamics for the ﬁrst and second moments of the organelles’ position along the cell. Furthermore, we analyse models with power-law distributed random times, which reveal the prevalence of the Mittag-Lefﬂer resting state and its contribution to subdiffusive and superdiffusive behaviour. Finally, we demonstrate a non-Markovian non-additivity effect when the switching rates and transport characteristics depend on the rate of organelles removal. The analytical calculations are in good agreement with numerical Monte Carlo simulations. Our results shed light on the dynamics of intracellular transport and emphasise the effects of rest times on the persistence of random walks in complex biological systems.


Introduction
Transport processes within eukaryotic cells bear striking similarities to the intricate networks found in large cities.Cytoskeletal networks, analogous to roads and highways, crisscross cells, and facilitate the transport of various cargoes.Intracellular transport is a fundamental biological process involving the movement of molecules, proteins, and organelles from one location to another within cells.It relies on a combination of passive diffusion and ATP-driven movement along cytoskeletal filaments.This process is crucial to maintaining cellular structure, function, and survival, and is involved in essential cellular activities such as growth, division, and signaling.Dysregulation of intracellular transport has been implicated in numerous diseases, including cancer.Therefore, understanding the mechanisms and regulation of intracellular transport is of great importance for both basic cell biology research and clinical applications.In neurons, long-range intracellular transport is even more critical, as its malfunctioning has been linked to adult-onset neurodegenerative diseases, including Alzheimer's, Parkinson's, and Huntington's disease [1].
Vesicles and organelles have been observed to undergo bidirectional motion along cytoskeletal filaments, periodically pausing for finite time intervals.These cargoes utilise molecular motors of opposite polarity in rapid succession [2,3], facilitating their transportation over long distances [4].The two major families of motors involved in this process are kinesins and dyneins [5,6].Kinesins transport vesicles along microtubules toward the plus ends, facilitating material transport from the cell interior toward the cortex (anterograde movement).Dyneins move material toward the microtubule minus ends, moving from the cell periphery to the cell interior (retrograde movement).Cytoplasmic dynein predominantly facilitates the minus-end directed long-range transport along microtubules (MTs) [7], and its activity is regulated by multifunctional adaptors such as dynactin [8].
Various modeling approaches have been reviewed to understand intracellular transport [2,9,10], including Markovian models such as reaction-diffusion models [11,12], Brownian ratchet models [13], random walk models [14,15], intermittent search processes [16,17], and exclusion models [18][19][20].While most of these approaches treat microtubules implicitly, some models explicitly consider them, allowing the investigation of the influence of cytoskeleton topology on intracellular transport [21,22].A Markovian model of particle diffusion with spatially varying diffusivity shows non-ergodic behaviour [23] frequently observed experimentally in intracellular transport, which was also modelled using non-Markovian fractional Brownian motion [24].Several multi-state Markovian models exhibit anomalous diffusive regimes before transitioning into normal diffusion at longer times [22,25,26], or even subdiffusion due to topological trapping at longer times [21] or due to non-Markovian run-length-dependent detachment rate of cargo from a microtubule [27].The effect of advection on diffusive particles was studied using the mobile-immobile model where the particle switch between mobile and immobile states with finite rates [28] and with power-law and mixed trapping time distributions [29].The non-Markovian mobile-immobile model without advection was used to model the diffusion of excitons in layered perovskites and transition metal dichalcogenides [30].
Given the complexity of intracellular space, cytoskeleton topology, and various interactions within cells, it is not surprising that intracellular transport deviates from standard Brownian diffusion and displays anomalous diffusion [31][32][33][34][35][36][37].In cellular and molecular biology, anomalous diffusion is also highly heterogeneous [38].Experimental studies have observed different non-Brownian regimes of diffusion at various time scales [24,27,[39][40][41][42][43][44][45][46][47][48][49][50][51][52].The intracellular environment is highly crowded with macromolecules, subcellular compartments, and confinement domains, leading to subdiffusion at small time scales (less than tens of milliseconds) [53][54][55][56] which transitions to superdiffusion associated with the activity of molecular motors.Whether this superdiffusion persists at longer time scales or if it represents an intermediate regime before transitioning into further normal or subdiffusion regimes is still a subject of debate.Also, the long-time subdiffusion could originate from different mechanisms.At longer times, the influence of cytoskeleton topology [21] and the heterogeneity of the transport process come into play, leading to spurious subdiffusion [24].In spite of this, the anomalous nature of intracellular transport suggests that non-Markovian models could be viable.
The non-Markovian nature of intracellular transport has been modelled using random walks [27], fractional Brownian motion [24,51,52] and Levy walks [57] and combinations of FBM and CTRW [46] and CTRW and Levy walk [49].The later stochastic model consists of altering phases of active motion with constant velocity interchanged with periods of passive dynamics.Depending on the distribution of times in the two phases of motion, the model shows transitions from normal diffusive to superdiffusive and to ballistic behavior [58][59][60][61].Recently, Han et al. [62] (see also [63]) proposed a three-state self-reinforcing random walk model involving Mittag-Leffler distributed rest: Here, p + and p − are PDFs of particles moving in the positive and negative direction.The running times are exponentially distributed with the rate λ.The PDF p 0 (x, t) describes the rest states drawn from the Mittag-Leffler distribution with parameter β and D 1−β t denotes the Riemann-Liouville fractional derivative given by: Three densities satisfy the following normalisation condition: Parameters r + and r − describe the probabilities of transitions from the rest.In the self-reinforcing random walk model [62], these probabilities depend on time and space.This three-state self-reinforcing random walk model is characterised by superdiffusive behaviour at intermediate times and by subdiffusion at a longer time scale.
Clearly, in reality, the speeds of cargoes in anterograde and retrograde directions are not equal, since they are powered by different motors.Therefore, in this paper, we analyse intracellular transport using a one-dimensional three-state persistent random walk model with different speeds of cargoes in anterograde and retrograde directions v + and v − that generalises the system (1)- (3).We suggest a model in which cargo transitions between moving states and resting states with zero velocity are non-Markovian.We consider the case when the rate of transitions depends on the running time.In order to take into account the random duration of cargo trajectories observed in experiments [24,52], we introduce a removal rate that terminates trajectories at random times.This effect was not considered in a recent paper [64].The constant removal rate corresponds to the exponentially distributed duration of trajectories.The structure of the paper is as follows.In Section 2, we set up the Three-State Model with a residence time variable and derive the equations for the structured densities in two active states and one passive state.In Section 3, the mesoscopic master equations for densities in all three states are introduced.Specific examples of the Markovian Three-State Model and the Three-State Model with Mittag-Leffler distributed resting times are considered, and equations are derived for the first and second moments.The analytical results of Section 3 are compared in Section 4, in which we perform Monte Carlo simulations.Finally, we discuss our results in Section 5 and give a summary and conclusions in Section 6.

Three-State Transport Model with a Residence Time Variable
We consider a model in which cargo randomly switches between three possible states: Two states in which the cargoes move with a constant speed ν + > 0 (state A) or ν − > 0 (state R) and a resting state (state 0) in which the cargoes do not move.In this section, we aim to develop a stochastic intracellular transport model to describe the non-Markovian dynamics of intracellular transport.
We introduce the densities that are functions of running and residence time τ [65][66][67][68].Consider n A (x, t, τ) and n R (x, t, τ) the densities of cargoes at point x at time t whose residence time in an active state (state A or state R, respectively) lies in the interval (τ, τ + dτ).The density of cargoes in a resting state is n 0 (x, t, τ).The balance equations for structured densities can be written as [66]: where λ A (τ), λ R (τ) and λ 0 (τ) are the switching rates that depend on τ and θ A , θ R and θ 0 are the constant removal rates.We assume that the residence time of all cargoes at t = 0 equals zero.In this case, the initial conditions are: where p 0 j (x) represents the initial densities (0 ≤ τ ≤ t).Three escape (switching) rates, λ A (τ), λ R (τ), λ 0 (τ) appearing in ( 5)-( 7) can be defined as [69]: where ψ j (τ) are PDFs of the residence times in corresponding states and Ψ j (τ) = ∞ τ ψ j (u)du are the survival functions.Our approach is very flexible and can be used for various forms of the residence time probability distribution.This flexibility is important, since recently, the non-exponential PDF of the residence time was found in experiments [27].
The boundary conditions at zero running time (τ = 0) can be written as follows: These equations describe the switching process to the states A, R, and resting state 0. If the cargo moves in the positive direction, it can switch with rate λ A (τ) to the opposite direction with the probability α R , to the resting state with the probability α 0 , or continue in the same direction with the probability α A .The cargo moving in the negative direction can switch with rate λ B (τ) to the opposite direction with the probability β A , to the resting state with the probability β 0 , or continue in the same direction with the probability β R .Finally, for the resting cargo, it can switch with rate λ 0 (τ) to cargo moving in the positive direction with the probability γ A , to cargo moving in the negative direction with the probability γ R , or remain at rest again with the probability γ 0 .Obviously, The product λ j (τ)n j (x, t, τ) gives the escape rate corresponding to a particular residence time τ.If we denote the escape rates from the states A, R and 0 by i A (x, t), i R (x, t) and i 0 (x, t), correspondingly.They can be obtained by integrating λ j n j over variable τ from 0 to t [67]: It follows from Equations ( 10)-( 12) and ( 14) that We solve Equations ( 5)-( 7) using the method of characteristics: For active states (0 For resting state (0 One can see that all solutions ( 18)-( 20) contain an exponential factor e − τ 0 λ j (s)ds , which is recognised as the survival function Ψ j (τ): Using equations ( 10)-( 12), the solutions in ( 18) can be rewritten using the survival function from ( 21) and the escape rates (switching terms) as follows: Notice that the residence time PDF, ψ j (τ), is related to the switching rate, λ j (τ), as follows [69]:

Mesoscopic Master Equations
Now, we introduce the mesoscopic densities for all three states: p A (x, t), p R (x, t) and p 0 (x, t).The functions p A and p R are the densities for cargoes moving in a positive/negative direction and p 0 is the density of cargoes with a resting state.These densities can be obtained by integrating the structured densities n j (x, t, τ) over residence time variable τ: By differentiating (26) with respect to time t and using ( 5)- (7), we obtain the mesoscopic master equations for p A , p R and p 0 : Using ( 10)-( 29), we rewrite the system of equations using the escape rates i A , i R , and i 0 : The escape rates i A , i R , and i 0 can be expressed in terms of p A , p R and p 0 , respectively, as follows (see Appendix A for the details of the derivation): where H i (t) is the memory kernel defined by its Laplace transform [70]: Here, ψj (s) and Ψj (s) are the Laplace transforms of the residence time PDF ψ j (τ) and the survival function Ψ j (τ), respectively.
It is clear from ( 33)-( 35) that the escape rates i A (x, t), i R (x, t) and i 0 (x, t) depend on the rates of organelles removal.This is the non-Markovian non-additivity effect that leads to the dependence of the transport characteristics on the kinetics of organelles removal [70][71][72].Note that the equations above can be applied for the analysis of the phenomenon of the migration-proliferation dichotomy in gliomas [73][74][75], anomalous transport in spiny dendrites [70].
Next, we will examine various distributions for residence time PDFs, ψ k (τ), including power-law distributions.We will derive the equations for the moments of the random walk position where p(x, t) = p A (x, t) + p R (x, t) + p 0 (x, t).

Markovian Three-State Model
In this subsection, we consider the Markovian case when the switching rates for all three states, λ j (τ), are constant.The residence time's PDFs are exponential: with the Laplace transform ψj (s) = λ j λ j +s .The escape rates i A (x, t), i R (x, t) and i 0 (x, t) do not involve memory effects and can be written as follows: The balance equations for the mean densities p A (x, t), p R (x, t) and p 0 (x, t) are: A similar Markovian model describing macroscopic intracellular transport of vesicles and organelles has been introduced in the classical paper [11].Note that without the resting state, the system of ( 40), ( 41) can be reduced to the generalised telegraph equation [66].

Three-State Model with Mittag-Leffler Distributed Resting Times
Now, let us consider another example involving anomalous behavior.For simplicity, we assume that the removal rate θ 0 for the rest state is zero.We consider the case when the cargoes move along microtubules with running times that follow an exponential distribution with rate λ in the active states and pause for a Mittag-Leffler (ML) distributed residence time in the resting state.In this case, the PDF for the residence time exhibits a power law such that the mean residence time does not exist.The survival function can be written in terms of the one-parameter ML function, E µ (.), as: where τ 0 is a time scale.The Laplace transforms of the survival function Ψ 0 (τ) and residence time PDF ψ 0 (τ) are: Then, from (35), we obtain the rate i 0 as: where D 1−µ t is the Riemann-Liouville fractional derivative defined in (3).The system of integro-differential Equations ( 30)-(32) become: The system of fractional PDEs ( 46)-( 48) can be also rewritten in terms of Caputo fractional derivatives.In [62], the simplified version of ( 46)-( 48) has been studied, in which the authors derived an equation with a telegraph operator; see Equation (8) in [62].

Moments Equations
In this subsection, we calculate the first two moments defined in (37).We put θ A = θ R = θ 0 = 0 and introduce the following functions: The first moment, m 1 (t), and second moment, m 2 (t) can be determined as: where m iA (t), m iR (t), m i0 (t) can be found from the system of equations: From these equations, one can find the equations for the first and second moments Let us find the first moment m 1 .For simplicity, we consider the case when The equations for m 0A , m 0R and m 00 take the form: We take the Laplace transform of ( 56)-( 58), set p A (x, 0) = δ(x), p R (x, 0) = 0 and p 0 (x, 0) = 0 and obtain, in the long-time limit (s → 0), expressions for m0A From the normalisation condition (4), we find: and in the limit s → 0, we obtain: Finally, from ( 59) and ( 61) we obtain the first moment in the long time limit It follows from ( 62) that the anomalous rest state becomes dominant leading to subdiffusive motion for µ < 1/2.Now, let us find the second moment for the case when the first moment is zero; that is, The main aim is to show again that the rest state is dominant in the long-time limit that leads to subdiffusion.The equation for the second moment is where m 1A , m 1R obey the equations We take the Laplace transform and obtain in the limit s → 0 Performing straightforward calculations, we obtain Taking inverse Laplace transform, we obtain the second moment in the long time limit Clearly, this formula gives us the subdiffusive motion.
In the next sections, we performed the Monte Carlo simulations.

Monte Carlo Simulations
Monte Carlo simulations of the random walk follow standard procedure: (1) Set initial conditions x 0 = 0 and t 0 = 0.The initial state was randomly selected corresponding to ) For the Markovian Three-State Model, generate an exponentially distributed random time τ 0 = −1/λ ln(1 − p) where p is a uniformly distributed random number in [0, 1).
In the Non-Markovian Three-State Model, resting times were drawn from the ML distribution, and running states were generated using exponentially distributed ran-dom times.For the resting state, generate the Mittag-Leffler random number τ using the Matlab mlrnd function (Guido Germano (2023).Mittag-Leffler random number generator (https://www.mathworks.com/matlabcentral/fileexchange/19392-mittagleffler-random-number-generator(accessed on 1 September 2023)), MATLAB R2020b Central File Exchange).( 3) Update position and time to x 1 = x 0 + v 0 τ 0 , t 1 = t 0 + τ 0 , respectively.Update state by randomly selecting new velocity 2) and (3) until the predefined simulation time T max is reached.Thus, a single trajectory x(t) is generated.(5) Repeat steps ( 1)-( 4) N = 10 4 times to generate an ensemble of N trajectories.Note that trajectories have random durations due to step ( 2). ( 6) The ensemble of trajectories is then analysed by calculating the distribution of positions at a given T max .We estimated these distributions using histograms.To quantify the anomalous diffusion in the Three-State Model, we calculated the moments of these distributions as a function of time.In particular, the first and the second moments were calculated as m Here, (i) denotes the index of a trajectory.
To achieve the simulation results for the Markovian Three-State Model, we compute the variance m 2 (t) − m 2 1 (t) and distributions of walkers with parameters λ = 1, v + = 1 and v − = 1.The results are shown in Figure 1.We chose equal probabilities of switching from resting to running state γ A = γ R = γ 0 = 1/3.As expected, the PDFs of the walker's positions follow the Gaussian distribution.The corresponding second moments show ballistic regimes at short times and asymptotic linear growth.For the Non-Markovian Three-State model, we also find good agreement with the theoretical predictions.First, we calculated the first moment m 1 (t) of random walker position as a function of time with zero removal rates We considered two values of anomalous exponent µ = 0.3 and µ = 0.7.As shown in Figure 2, the numerics agree with the analytical Equation ( 62) without fitting.
The ML-distributed resting states make the distributions of the walker's position non-Gaussian.An example of the distributions of random walkers' position with the resting states generated from the ML distribution with the scaling parameter µ = 0.7 is shown in Figure 3.At long times, the MSD is growing sub linearly with time, m 2 (t) ∼ t β with β = µ.At the same time, it was previously found that the two-state model with one running state and the ML resting state has the long-time behaviour of the variance m 2 (t) − m 2 1 (t) ∼ t β with β = 2µ.So, for µ < 0.5 the variance grows subdiffusively with time, and for µ > 0.5 the growth is superdiffusive.
We note that the two-state model is equivalent to the three-state model with the speed of one running state equal to zero.The following question arises: what would be the long-time behaviour of the three-state model with ML-distributed resting times with unequal retrograde and anterograde speeds?We calculated the variance of the position for the three-state model with ML-distributed resting times and showed that the model has subdiffusion or superdiffusion at long times depending on the parameters.At long times, the variance depends on cargo speeds and is well fitted by the power-law function: For v + = v − , we recover the subdiffusive behaviour β = µ < 1 predicted by the Equation (70).The dependence of β/µ − 1 on the speed v + (with v − kept constant) for the three-state model with ML-distributed resting state is shown in Figure 4a.The asymptotic behaviour of different models given by the dependence of exponents of MSDs or variances at long times, β, is shown in Figure 4a.In contrast to the three-state model with equal retrograde and anterograde cargo speeds v + = v − , the three-state model with unequal retrograde and anterograde cargo speeds, v + = v − can have subdiffusive behaviour and superdiffusive behaviour.

Discussion
In cells, vesicles are transported by dyneins from the cell periphery to the cell interior (retrograde transport) and kinesins transport them from the cell interior toward the cortex (anterograde transport).Dyneins and kinesins move at different velocities.Thus, in contrast to previous works [62,64], here, we considered different cargo velocities in retrograde and anterograde directions.As a result, the first moment of the cargo position m 1 (t) is non-zero but grows linearly with time for the Markovian model or as a power law in the case of the non-Markovian model.In experiments, this could be used to distinguish between Markovian and non-Markovian dynamics.
The variance of cargo position m 2 (t) − m 2 1 (t) bears additional information about the cargo dynamics.For the non-Markovian model with unequal velocities, the variance can grow subdiffusively and superdiffusively depending on the anomalous exponent µ of the distribution of the duration of resting states.Combining predictions for the first moment, variance, and distribution of cargo's displacements, it would be possible to verify the non-Markovian three-state model using the ensemble of experimental trajectories.

Summary and Conclusions
We have introduced a non-Markovian persistent random walk model for intracellular transport that incorporates the removal rate and finite speed of propagation.This model includes random transitions between two active states (anterograde and retrograde movements) and one resting state, with the probability of switching depending on the amount of time the organelle spends in each state.We have derived new mesoscopic integro-differential equations for densities of organelles inside the cell by using the structural densities obeying Markov differential equations of the first order.We have found a non-Markovian non-additivity effect when the switching rates and transport characteristics depend on the rate of organelles removal.Ordinary differential equations describing the dynamics for the first and second moments of the organelles' position along the cell have been derived.A model of stochastic transport of organelles which takes into account experimentally observed Mittag-Leffler resting time distribution has been analysed in detail.In particular, it determines the subdiffusive and superdiffusive behaviour of the first moment.Notice that the fractional Riemann-Liouville derivative naturally arises as the result of the Mittag-Leffler distribution of residence times in the resting state.For a general form of resting time distribution, one should not expect the escape rate to be expressed with any kind of fractional derivative.Rather, it will be given by an integral form defined by Equation (34) A numerical Monte Carlo algorithm has been developed to calculate moments and organelle mesoscopic densities in cells, taking into account non-Markovian effects caused by the dependence of transition rates on residence times.
In this paper, we propose the mathematical model of intracellular transport.This model is based on experimental observations (see, for example [2,6]) that intracellular transport demonstrates so-called stop-and-go behaviour and is characterised by long time sub-diffusive behaviour.In this work, we modelled the stop-and-go behaviour as switching from active to resting states and showed that the ML distributed time spent in resting states leads to sub-diffusion at long times.So, experimentally, it would be possible to test the predictions of our model, e.g., whether the resting times are described by the ML distribution.In fact, in our previous work [27,57], we found that the distribution of resting times indeed might be better fitted by a power law (notice that the long-time asymptotic of the ML distribution is t −1−µ ).However, more experimental data are needed to make this comparison more reliable.
In this paper, we used the assumption of constant cargo speeds.In experiments, intracellular organelles might be described based on the distribution of speed.Our constant speed assumption could be viewed as utilising the first moment of this distribution.In the future, it will be interesting to extend our work by taking into account the distribution of cargo speeds [51] and the heterogeneity of intracellular transport [24].

Figure 1 .
Figure 1.Markovian three-state model.(a) The distribution of random walkers is calculated at T max = 10 with zero removal ratesQ A = Q R = Q 0 = 0 (circles) and Q A = Q R = Q 0 = 0.1 (squares).Other parameters are given in the text.The dashed line represents the Gaussian distribution with Figure 1.Markovian three-state model.(a) The distribution of random walkers is calculated at T max = 10 with zero removal ratesQ A = Q R = Q 0 = 0 (circles) and Q A = Q R = Q 0 = 0.1 (squares).Other parameters are given in the text.The dashed line represents the Gaussian distribution with σ 2 = 2DT max .(b) The corresponding second moment m 2 (t) as a function of time (blue and red curves) grows ballistically for short times and switches to linear growth with the effective diffusion coefficient D.

Figure 3 .
Figure 3. Non -Markovian three-state model.(a) The distribution of random walkers is calculated at T max = 10 with zero removal ratesQ A = Q R = Q 0 = 0 (circles) and Q A = Q R = Q 0 = 0.1 (squares), τ 0 = 1, λ = 1 and v + = 1, v − = 0.1.The dashed line shows the Gaussian distribution for comparison.(b) The corresponding mean squared displacements (blue and red curves) grow ballistically for short times and switch to sub-linear growth m 2 (t) ∼ t β with β = µ and the exponent

Figure 4 .
Figure 4. Three-state model with ML-distributed resting times can have subdiffusive or superdiffusive asymptotic behaviour.(a) The dependence β/µ − 1 as a function of the speed v + (with v − kept constant) for the three-state model with ML-distributed resting state.The dashed line corresponds to Equation (71) and symbols show the results of numerical simulations.(b) The diagram of the asymptotic behaviour of different models is given by the dependence of exponents of MSDs or variances at long times, β, as functions of the exponent of ML distribution of resting states, µ.The line with the biggest slope corresponds to the two-state model (that is, the three-state model with either v + = 0 or v − = 0) with ML-distributed resting states.It has subdiffusive behaviour β < 1 for µ < 1/2 and superdiffusive behaviour β > 1 for µ > 1/2.The lowest line corresponds to the three-state model with equal retrograde and anterograde cargo speeds v + = v − .The line with an intermediate slope corresponds to the three-state model with unequal retrograde and anterograde cargo speeds, v + = 1 and v − = 0.25.