A stochastic model for intracellular active transport

We develop a stochastic model for an intracellular active transport problem. Our aims are to calculate the probability that a molecular motor reaches a hidden target, to study what influences this probability and to calculate the time required for the molecular motor to hit the target (Mean First Passage Time). We study different biologically relevant scenarios, which include the possibility of multiple hidden targets (which breed competition) and the presence of obstacles. The purpose of including obstacles is to illustrate actual disruptions of the intracellular transport (which can result, for example, in several neurological disorders [11]). From a mathematical point of view, the intracellular active transport is modelled by two independent continuous-time, discrete space Markov chains: one for the dynamics of the molecular motor in the space intervals and one for the domain of target. The process is time homogeneous and independent of the position of the molecular motor. Raluca Purnichescu-Purtan Univeristy Politehnica of Bucharest e-mail: raluca.purtan@gmail.com Irina Badralexi Univeristy Politehnica of Bucharest e-mail: irina.badralexi@gmail.com Citation: Raluca Purnichescu-Purtan, Irina Badralexi, A stochastic model for intracellular active transport, in S. Markov (Editor), Conference Proceedings of Biomath 2018, Biomath Forum, Sofia, 2018, pp. 1-13, http://dx.doi.org/10.11145/texts.2018.10.277 Copyright: c © 2018 Purnichescu-Purtan et al. This article is distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

This study is related to intracellular transport, which is the directed movement of substances within a cell.Microtubules (microscopic hollow tubes) and microfilaments constitute a part of a cell's cytoskeleton.These represent the roads in intracellular transport.
An important element in vesicle transport is motor proteins.These proteins bind to vesicles and organelles and move them along the microtubule or microfilament network.The motor proteins which move vesicles along microtubules are called kinesins and dyneins, and those which move vesicles along microfilaments belong to the myosin family.
Issues with intracellular transport can have serious consequences.They have been linked to neurological disorders such as Alzheimers disease ( [11], [21]).
In the field of molecular biology, such models were designed to describe active transport of reactive chemicals in cells ( [15]), promoter protein searching for a specific target site on DNA ( [7], [9]), transport of mRNA ( [8]).
Depending on the context of the activity which is being modeled, one can consider a one dimensional case, which means that the movement only occurs from left to right or vice-versa ( [4]), or a higher dimensional case (two dimensional or three dimensional), where the movement is less restricted ( [5]).
The intermittent motions that occur at the microscopic level of reaction kinetics within biological cells are modelled in a one-dimensional framework.Previous models were based on the assumptions that the search begins from a random point within a specified interval and that the target is reached upon entering the target interval ([3], [6], [15]).Also, many of the studies regarding intermittent search problems in microbiology consider bidirectional transport and the condition that the target can be reached from a specified type of movement -ballistic (anterograde, with constant velocity) or during a diffusive phase ( [8], [16], [2], [13], [10], [18]).
Unlike the previous studies, in our research, the unidirectional active transport is modelled by two independent continuous-time, discrete space Markov chains: one for the dynamics of the molecular motor in the space intervals (outside the target domain) and one for the domain of target.The novelty of our study consists of considering two types of motions inside and outside the target domain: a brownian motion (from which the target can be reached) and a state of active motor-driven transport along microtubules or microfilaments.We also consider an imperfect detection of the target, the scenario of multiple targets (competitivity) and the presence of obstacles.A stochastic algorithm was also developed.
For each scenario we calculate the probability that a molecular motor reaches the target and calculate the time required for the molecular motor to hit the target (Mean First Passage Time -MFPT).
2 Intracellular transport problem − one target domain

General framework for one target domain
Consider a single motor-driven particle moving along a one-dimensional path of length L (unidirectional or anterograde movement).A "hidden target" is located at a fixed (known) location on the path.The term "hidden" means that the molecular motor may detect the target only when it enters into the target domain which lies fully within the interval I = [0, L], i.e., there is some positive constant ε > 0 for which (a − ε, a + ε) ⊂ I.The mathematical framework for our stochastic model is based on the following general assumptions: A1.For the molecular motor we consider two movement regimes, denoted by: • "u" -uniform one-dimensional movement with velocity v > 0; • "b" -one-dimensional Brownian motion (modeled as Continuous Time Random Walk).
A2.The particle motion is subjected to random decision moments which, independent from the time points, could change between the movement regimes; A3.The detection of the target (the searcher is "absorbed" by the target) can be done only from Brownian motion, in the interval (a − ε, a + ε) ⊂ I, and may take place with a given probability p (we consider an "imperfect detection" there is always a possibility that the target remains undetected); A4.The first movement regime (at time t 0 = 0) is the uniform one; A5.There is a finite time of observation [0, T ] (the particle is absorbed by the target, leaves the path though the right point or the observation time elapses).

Choosing and changing the movement regime
We consider a homogeneous finite Markov chain 1 with the state space S = {u, b, d}, where u and b represent the movement regimes of the molecular motor and d is the absorbing state (the target is reached, the searcher remains in state d with probability one); We shall denote the cardinal of the set S by card(S) ; For the time interval [0, T ], we consider a division 0 = t 0 < t 1 < ... < t n = T (the points t 1 , ...,t n−1 and the number n of the division points are random, see below the connection with the holding time).At t 1 , ...,t n−1 the particle probabilistically "decides" (independent from the time points) if it remains in the previous movement regime or if it changes to the other regime (except from state d).If the state d is not reached by t n = T , we shall consider that the detection of the target failed.
The n -step trajectory of the Markov chain is s 0 , s 1 , ..., s n−1 , where s k ∈ S is the state for the time interval The state holding time is denoted by T (k) = t k+1 − t k , for k = 1, n − 1 (a random exponential variable with parameter λ ).The value of the parameter is connected to the values of the diagonal elements of the rate transition matrix Q (known): The transition probability matrix is: and the connection between the transition probabilities, the parameter λ and the transition rates is:

The movement of the molecular motor
The molecular motor will move with s k regime on the interval (t k ,t k+1 ], and we have two possibilities: First, if s k = u, the space position (the net displacement) of the "searcher" is a continuous function of time and it is obtained by integrating the velocity function, as a function of time: Alternatively, if s k = b, the time interval (t k ,t k+1 ] is divided into a sequence of times τ 0 , τ 1 , ..., τ m (with τ 0 = t k and τ m = t k+1 ); then, for any fixed τ i , i = 1, m, the space position (the net displacement) of the "searcher" is given by: where y i are values of a standard normal (Gaussian) random variable, Y ∈ N(0, 1).

The general stochastic algorithm for one target domain
For the stochastic algorithm, we consider the necessary input parameters: L, T , a, ε, δ τ, p, λ and the rate transition matrix Q = (q i j ) i, j=1, 2 .In what follows, we briefly outline the steps of the algorithm.
Step 0: We construct the transition probability matrix P = (p i j ) i, j=1,2 , according to (1), for a fixed value of λ .
Step 1: The particle enters the path and is moving in the space interval According to the general assumption A4., at t 0 = 0 we have s 0 = u and x(0) = 0 (the space position of the molecular motor).The holding time in states 0 = u is generated using the Integral Inversion Theorem2 ([1]): where µ is a value of the random variable U ∈ Uni f (0, 1).During this time, the particle is moving uniformly with the constant speed v along the path.After this time elapses, the new state is chosen (u with the probability p 11 or b with probability p 12 ; a Monte-Carlo discrete sequence is used).The position of the particle after its exits from the state u is calculated according to (2).
Step 2: While x(t k ) ≤ a − ε, with k > 0: • Randomly select the new state, using the corresponding row of the probability transition matrix P; • Generate the sojourn time in the current state,T (0 • If the new state is b, generate a value y from a standard normal (Gaussian) random variable, Y ∈ N(0, 1) ; Generate a sequence of equal step times τ 0 , τ 1 , ..., τ m (within the sojourn time T (k)), with τ 0 = t k and τ m = t k+1 ; during each time step δ τ = τ i+1 − τ i , i = 1, m, the particle is moving in one dimension, according to (3); • Update the new position of the particle and the time.Step 4: The particle is moving in the space interval If the state d is not reached in Step 3 and the time has not elapsed, the motion of the molecular motor is simulated as in Step 2.
• If the molecular motor re-enters the target interval from state b, go to Step 3; • The algorithm stops when the particle reaches the end of the path, L, or the time has elapsed; the fact that the target was not reached is recorded.
Remark.The algorithm can be easily adapted for the case of different holding times for the states u and b (in the sequence where the holding time is generated, two different values, λ 1 and λ 2 should be used).

General settings for the simulations
Depending on the domain of interest, an intermittent search strategy has its own natural mechanism.The use of the algorithm in biology, for example, requires some conditions to be met.Thus, the net displacement of the particle in the state u always exceeds the maximum net displacement of the particle during the Brownian motion state b: Also, the length of the target domain has to be proportional to the maximum net displacement of the particle during the Brownian motion state.For our simulations, we consider: To ensure the adaptability and flexibility of the algorithm, space and time are measured in arbitrary units; 3 Intracellular transport problem − two target domains (competitivity) We shall consider the same general framework as in the "one target domain" case, but with two "hidden targets", located at two fixed but unknown locations on the path, denoted by a and b.The target domains are (a − ε 1 , a + ε 1 ) ⊂ I and (b − ε 2 , b + ε 2 ) ⊂ I, for some positive constants ε 1 , ε 2 > 0 and the "competitivity" is modelled within the overlapping target interval I 4 .If the molecular motor is in any of the non-overlapping target domains I 2 or I 3 , it can be absorbed with two different given probabilities p 1 and p 2 (we consider an "imperfect detection").
Fig. 2 Schematic presentation of the model system with two targets domain.
The general assumptions A1. − A5. are valid and the state spaces and the transition probability matrices for the Markov chains are different for each interval.
For intervals I 1 and I 5 , the state space is S = {u, b} and the transition probability matrix is In this case we left the values for the transition probabilities unchanged from state u (the first row).We consider the same regime of absorption for both targets (the probability of absorption is p 1 = p 2 = p * 23 ).If there is some relevant biological reason for considering different absorption rates, p 1 = p 2 , we should consider two different transition probability matrices of type P 2 .
For interval I 4 , the state space is S = {u, b, d 1 , d 2 } and the transition probability matrix is In this case, the absorption probabilities are p 1 = p * * 23 and p 2 = p * * 24 .These probabilities can be considered equal or different, depending on the biological model.

Intracellular transport problem − one target domain in the presence of one obstacle
We shall consider the general framework of the "one target domain" case, including an "obstacle" particle, with a fixed but unknown location on the path, denoted c, with an obstacle interval (c − ε 2 , c + ε 2 ) ⊂ I, for some positive constant ε 2 > 0. We shall assume that the obstacle domain is located on the path before the target domain (a − ε 1 , a + ε 1 ) ⊂ I (or after the target domain, according to a biologically relevant scenario).We consider that the molecular motor can hit the obstacle only from Brownian motion and the result of "collision" may result in retrograde uniform movement (or other biologically relevant movement could be considered); The general assumptions A2. − A5. are valid.A1. should be modified as follows: A1 * .For the molecular motor we consider three movement regimes, denoted by: • "u 1 " -uniform one-dimensional movement with velocity v 1 > 0 (anterograde); • "b" -one-dimensional Brownian motion; • "u 2 " -uniform one-dimensional movement with velocity −v 2 < 0 (retrograde); For intervals I 1 , I 3 and I 5 , the state space is S = {u 1 , b} and the transition probability matrix is: For interval I 2 , the state space is S = {u 1 , b, u 2 } and the transition probability matrix is: In this case, we left the values for the transition probabilities from state u 1 unchanged and we consider that from the retrograde uniform movement, the molecular motor can reach only the states u 1 or b.
For interval I 4 , the state space is S = {u 1 , b, d} and the transition probability matrix is:

Numerical simulations
Numerical simulations were conducted only in the case of one target domain and that of two target domains.Simulations for the movement of the molecular motor in the presence of obtacles, complete with interpretations and discussions in relevant biological situations, will be addressed in future work.We are first interested in the variation of the hitting probability and the MFPT due to the different location of the target domain.We choose the first position of a (center of the target domain) at a = 1 and "move it" with 0.5 units until the position a = 9.For every situation, ε = 0.0025 and N = 10000 runs.
We notice that there are no significant variations of the hitting probability due to a different location of the target, as long as the other parameters are fixed.However, the MFPT increases (linearly) with the distance from 0 to the target domain.(see Table 1) Next, we focus on the variation of the hitting probability and the MFPT due to the different velocities of the uniform movement.
Because the position of the target domain has no influence on the hitting probability, we choose the target domain (a − ε, a + ε) = (2.995,3).Table 2 The hitting probability and the MFPT depending on the position of the velocity of the uniform movement.

Velocity for the uniform movement
As expected, both the hitting probability and the MFPT increase as the velocity of the uniform movement decreases (see Table 2).
For the intervals [0,6) and (6.3,8], the computed transition probability matrix is: P 1 = 0.2857 0.7143 0.1429 0.8571 . For the target domains (6,6.2) and (6.06,6.3), the computed transition probability matrices are: For the "competitive" domain (6.06,6.2), the computed transition probability matrix is: We focused on comparing the hitting probabilities and the MFPT for the two targets, taking into account our simulation settings.
We observe that, for equal absorbtion probabilities, the location of the target is important (target a -which is closer to the starting point -has a greater hitting probability than target b -which is further from the starting point).Likewise, our simulations show that the MFPT is related to the location of the target.(see Table 3) Table 3 The hitting probability and the MFPT in a two target domain.
Remark.In both scenarios, regardless of the other input data, if the ratio between the probability of entering in the state b and the probability of entering in the state u (at a time point) is less or equal to one (in the transition probability matrices), then the probability of hitting the target is close to 0 (the search is inefficient).

Conclusion and future work
We develop and analyze a stochastic model of directed intermittent search for a hidden target on a one-dimensional path within the framework of continuous time, discrete state Markov chain.In order to gain flexibility and adaptability to different relevant biological situations, we modelled the stochastic movement of the molecular motor in the target domain with a different Markov chain and we have also considered the possibility of different holding times for the Markov processes involved.
In the Markov chain approach, we chose to define the transition probabilities using the transition rates and the holding times.This can be of great interest from a modelling point of view: one can compare the dynamics of the process under two different hypothesis: first, considering that the holding times in the distinct states are different (exponential random variables with different parameters λ ) and second, using the same exponential distribution for all holding times (with a single parameter λ ).
The novelty of our work consists of the framework formulated for three realistic situations for a molecular motor to reach its target.This target lies in a target domain (once the particle reaches this domain, it is "absorbed" by the target with a known probability).
The first of the situations describes the movement of the moelcular motor in search of a single target.For the second one, we considered one molecular motor particle which can be absorbed by two different targets (which are seen in competition).The third one introduces the search of one target in the presence of one obstacle (in this case, only the theoretical approach is presented; a more detailed study will be conducted in future work).
For our framework, we considered two main types of movement: a uniform movement with constant velocity and a brownian motion.In the case of an obstacle, one more movement is introduced: a retrograde uniform movement with constant velocity.
We deveopled a highly adaptable and flexible algorithm for the study of all three situations of intermittent search introduced.The main output of the algorithm is the hitting probability and the MFPT.
The proposed algorithm allows for imperfect detection of the target in its domain, which can be easily converted in a unit probability -perfect detection -according to different biologically relevant assumptions.
In the case of competitive targets, multiple scenarios can be configured through the modification of the ratio between the probabilities of absortion, the imperfect detection or the distinct holding times.This flexibility can lead to determining the parameters which significally influence the hitting probability and the MFPT.
There are many interesting connections to be made between the input parameters of the model and the output, depending on the biological data available and this will be our purpose in a future work.We will also give attention to a combination of multiple obstacles and multiple targets, in a relevant biological context, based on the stochastic algorithm presented in this paper.

Fig. 1
Fig.1Schematic presentation of the model system with one target domain.

Step 3 :
The particle is moving in the target domain, interval I 2 = (a − ε, a + ε) ⊂ [0, L].The first state in this interval is corresponding to the last "new" state from Step 2. in order to include the new possible state d (the absorbing state), we construct the new transition probability matrix: The values of the new transition probabilities can be completely different from the values in Step 1 or can be different just in the second row (p 21 = p * 21 and p 22 = p * 22 + p * 23 , for example); The detection of the target may take place with a given probability p = p * 23 ; As long as a − ε < x(t k ) < a + ε, with k > 0: • If the state d is not reached, see Step 2 (with the matrix P * instead of P); • If the state d is reached, the algorithm stops and the time and position of the particle are recorded; • If the time elapsed, the algorithm stops and the fact that the target was not reached is recorded.

P 1 = p 11 p 12 p 21 p 22 .
For intervals I 2 and I 3 , the state spaces are S = {u, b, d 1 } or S = {u, b, d 2 } and the transition probability matrix is

Fig. 3
Fig. 3 Schematic presentation of the model system with one target domain and one obstacle.