{'fraction_non_alphanumeric': 0.04406355667207259, 'fraction_numerical': 0.01858152870968084, 'mean_word_length': 3.8647450110864745, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 36, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 42, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Stochastic Compartmental Model for Fast Axonal Transport

In this paper we develop a probabilistic micro-scale compartmental model and use it to study macro-scale properties of axonal transport, the process by which intracellular cargo is moved in the axons of neurons. By directly modeling the smallest scale interactions, we can use recent microscopic experimental observations to infer all the parameters of the model. Then, using techniques from probability theory, we compute asymptotic limits of the stochastic behavior of individual motor-cargo complexes, while also characterizing both equilibrium and non-equilibrium ensemble behavior. We use these results in order to investigate three important biological questions: (1) How homogeneous are axons at stochastic equilibrium? (2) How quickly can axons return to stochastic equilibrium after large local perturbations? (3) How is our understanding of delivery time to a depleted target region changed by taking the whole cell point-of-view?


Introduction
In all cells, one finds that proteins, membrane-bound organelles, and other structures (e.g. chromosomes) are transported from place to place at speeds much higher than diffusion. Though these transport processes are fundamental to cell function, many of the underlying mechanisms, organizational principles, and regulatory features remain unknown. Axonal transport is one of the best studied systems because the transport is basically onedimensional since axons are long and narrow. There are two speeds of axonal transport. Fast transport goes at speeds of roughly 0.2 to 0.5 meters/day 2 [20] [24], while slow transport goes at approximately 1 millimeter/day, the rate of axon growth and regeneration [5] [20]. The biology and principles of slow transport are not yet clear [5], but the basic mechanisms of fast axonal transport were discovered in the 1980s [2][3] [22] [34]. The model in this paper refers to fast axonal transport, which we will henceforth call axonal transport.
The axonal transport apparatus consists of vesicles which form reversible chemical bonds with motor proteins that bind reversibly to microtubules which run parallel to the long dimension of the axon [1]. When the vesiclemotor protein complex is assembled on the microtubule, the complex steps stochastically with step size approximately 8 nanometers for kinesin and dynein and 10 nanometers for myosin [6][10] [16] [31]. The vesicles enter from the cell body on microtubules and then detach and reattach to the transport mechanism at random times.
In this paper we propose a spatial Markov-chain compartmental model based on these dynamics. We will assume independence of the interactions, and exponential wait times between events. While we address the validity of these assumptions in the Discussion section, we consider this a useful "firstorder" approximation that permits study of the dynamics from both the perspective of individual vesicles as well as that of the full spatial system. Such a model unifies some earlier modeling efforts and can accommodate both qualitative and quantitative experimental data observed on multiple scales.
In much experimental work in the 1970s and 1980s, radio-labeled amino acids were put into the cell bodies continuously or for a few hours. The amino acids were incorporated into proteins that were packaged into vesicles and put on the transport system so that at later times radioactivity could be seen moving progressively down the axons. In the continuous infusion case, one would see a wave of radioactivity with a sharp but slowly spreading wavefront propagating at constant velocity down the axon. In the case of infusion for a few hours one would see at long times a slowly spreading pulse of radioactivity that looked normally distributed. It was to understand this behavior that Reed and Blum constructed PDE models for axonal transport [3][25] [26]. These models did not have traveling wave solutions, but the data certainly looked like approximate traveling waves. In [27] it was shown by a perturbation theory argument that, in the asymptotic limit where the unbinding and binding rates k 1 and k 2 get large, the solution approaches a slowly spreading traveling wave or a normal pulse. Recently, in a series of papers, Friedman and co-workers have introduced new PDE models and proved these results rigorously [12] [13] [14] [15].
Probabilistic models for axonal transport were introduced and used for simulations already in the 1980s [30] [32]. However, rigorous work began with Lawler [21] in 1995 and was continued by Brooks [4] who used a continuous time stochastic model to show that the distribution of an individual particle is a spreading Gaussian at large times. Brooks also proved tail estimates for the central limit theorem and used them to estimate the error from normal.

Summary of Results
In this paper we revise the existing probabilistic models in order to observe randomness in the system as a whole rather than exclusively from the particle's point of view. Our goal is to exploit experimental observations at multiple scales in order to make predictions about the perceived homogeneity of material along the length of the axon. To this end, we propose in Section 2 a continuous-time Markov chain queueing model for the axonal transport system and estimate the order of magnitude of the various parameters. The precise values of the parameters will vary for different molecular motors, for various types of cargo and for various animal species.
In Section 3 we take the individual vesicle point-of-view to recover previous results [4] [26] and show that this model does unify and extend the existing probabilistic and PDE models. Using standard results from renewal theory we calculate the mean velocity and the near-Gaussian wavefront spreading of the law of the vesicle's location. Since we assume independence of the particle interactions, the law of an individual is equivalent to the distribution of an ensemble of particles released at the same time. Therefore the PDE governing a spatial-continuum limit of the law of a single particle is the same as the PDE studied in [26] (see Section 3.3).
Subsequently in Section 4 we adopt the full spatial system perspective to quantify stochasticity along the length of the axon. We begin by calculating in Proposition 4.1 the stationary distribution of a flow-through system that has sustained input from the nucleus, while particles are removed upon reaching the distal end. The stationary distribution has a product Poisson structure which allows for seamless transition between spatial scales. Via the coefficient of variation we give a precise characterization in Section 4.2 of the experimenter's qualitative perception of homogeneity at the millimeter scale.
The final two sections deal with model predictions. In Section 5 we study the non-equilibrium dynamics of recovery. Due to the product structure of the law of the transient dynamics, behavior is determined by the 2Ndimensional ODE governing the means. From this we estimate the timescale of return to equilibrium as a function of the lengthscale of interest. In Section 6 we consider the effect of deposition of material in the cell membrane. We find that there should be an exponential loss of material along the length of the cell, and conclude that a new biological mechanism must be discovered to account for the homogeneity that is a hallmark of axonal transport experimental observations.

The model and its parameters
Let L be the length of the axon, divided evenly into N lateral sections each of length δ, equal to the step size of the motor protein. Within each section, we disregard any further spatial geometry and take the particles to be in one of two states: • an on-transport state with mean lateral velocity v, or • an off-transport state with lateral velocity 0.
We use a 2N -dimensional continuous time Markov chain to model the particle dynamics: • Q i (t) is the number of particles in the on-transport state in section i, • P i (t) is the number of particles in the off-transport state in section i.
Each Q i and P i has the natural numbers N = {0, 1, 2, . . .} as its state space. In the simplest model, we consider the following transitions and rates, • Lateral transport: (Q i , Q i+1 ) → (Q i − 1, Q i+1 + 1) at rate rQ i (t), where r = v/δ; • Switch from on-transport to off-transport: • Switch from off-transport to on-transport: (P i , Q i ) → (P i − 1, Q i + 1) at rate k 2 P i (t); • Production of new particles: (Q 1 ) → (Q 1 + 1) at rate q 0 r; • Removal of particles at distal end: (Q N ) → (Q N − 1) at rate rQ N (t).

Figure 1: Double Chain
The lateral transport rate, r, is inversely proportional to the length scale so that the mean number of particles per unit length is invariant with respect to rescaling δ. The rate of production, q 0 r, ensures that this mean number per unit scales with q 0 . A graph of the model is depicted in Figure 1.
In order to ensure the Markov property, we use exponential random variables for the waiting times between transition events. Specifically, we mean that after a given event we assign a new independent random variable to each of the 3N + 2 possible next events, exponentially distributed with the appropriate rate parameter. The system of values updates according to the transition associated with the minimum of these waiting times. Then we create a new set of exponential random variables and the process proceeds as before. 6 This is an exactly solvable model. The advantage of computing explicit formulas for quantities that can be observed in experiments is that the experimental data can then be used to determine the parameter values in the model. For the characterization of the approximate wavefront speed and spreading in Section 3 and the homogeneity calculations in Section 4 we need order-of-magnitude estimates for the parameters. Actual parameter values will certainly differ depending on the particular neural tissue and the particular particles being transported. However, we can get order-ofmagnitude estimates from existing data.
First we recall that fast transport has been observed to travel at speeds of 0.2 to 0.5 m per day. We can assume that the average velocity of particles while physically bound to microtubules is roughly 1 m per day, or 10 −6 m/s. We have already stated the assumption that the length scale of the individual steps satisfies δ ∼ 10 −8 m. This implies that the rate parameter should be r = v/δ = 10 2 s −1 .
We now turn our attention to the on-off rates rates k 1 and k 2 . These can be determined from experimentally observed run lengths on the transport system. Indeed, Dixit et al. [9] show that a typical run along microtubules for dinein and kinesin is on the order of 10 −6 m. We can compare this with the theoretical run length of the model to determine off-rate k 1 . Within the model, at each step on the transport mechanism the particle has a binary decision to jump laterally along the transport with probability r/(k 1 + r), or to jump off with probability k 1 /(r + k 1 ). The number of jumps along the transport system before jumping off is therefore geometrically distributed on the set {0, 1, . . .} with success probability r r+k 1 . It follows that average number of steps in the run is r k 1 and therefore the average run length is r k 1 × 10 −8 m. Setting this equal to the average experimental run length of 10 −6 from [9], we see that r k 1 ∼ 10 2 , implying that k 1 ∼ 1s −1 . As we will see in the computation of the stationary distribution in Section 4.1, the ratio of the expected number of particles on the track to those off the track is k 1 k 2 . Dixit [9] found that approximately 75% of the particles were motile so this ratio is approximately equal to 3. Since k 1 ∼ 1s −1 we see that k 2 ∼ 1 3 . It remains to estimate q 0 . We will see in Proposition 4.1 that the mean number of particles per unit length is (1 + k 1 k 2 )q 0 = 4q 0 . Of course, axons have a large variety of diameters and larger axons will have more vesicles per unit length so one expects a range of values for q 0 . However, examination of a large number of electron micrographs of axonal cross-sections (see for example [17], Fig. 3; [18]; [23]), which are typically 100 nm thick enables one to estimate the number of vesicles per 100 nm segment. This number is typically in the range of 10 to 100 which implies that there are 1 to 10 vesicles per "box" in our model. Therefore q 0 is in the range 0.25 to 2.5, for various axons.
We remark that we are ignoring some aspects of the physics and the biology of axonal transport. We are not including diffusion of the vesicles off the track. We are treating the microtubule track as though it were a single continuous entity from one end of the axon to the other, when it fact it consists of numerous, separated, microtubule fragments. And, we are ignoring retrograde transport and the details of the motor proteins. Nevertheless, this simple model will enable us to investigate the homogeneity questions that are the main goal of this paper.

Dynamics from the Particle Perspective
In this section we calculate properties of the stochastic dynamics by using standard theorems from queuing theory. In the δ → 0 limit, the law of the location of a single particle converges to the Green's function of a linear partial differential equation. This enables us to obtain, as a special case, the asymptotic behavior of the PDE models for axonal transport in various asymptotic limits.

The active transport mode
We first consider the simple case where the particle starts at X 0 = 0 and stays exclusively in active transport mode. Let X t ∈ {0, δ, 2δ, . . . , L} be the lateral position of a particle at time t and let n t be the number of jumps made by the particle as of time t. Observe, X t = δn t . Proposition 3.1. Let k 1 = k 2 = 0, and r = v/δ > 0, then X t ∼ P ois(vt). In particular, the mean velocity of the particle is given by 1 t E[X t ] = v. For any given t ≥ 0, in the limit as δ → 0 the position of the particle satisfies where B is a standard Brownian motion.
Proof. Since X t = δn t where n t t≥0 is a Poisson process with rate r = v/δ the set out results follow from the law of large numbers (LLN) and the central limit theorem (CLT) for a Poisson process. 8

The on/off dynamics
We now consider a particle which undergoes transitions from on-transport to off-transport state and back. Denote again by X t ∈ {0, δ, 2δ, . . . , L} the lateral position of a particle at time t and let n t be the number of lateral transition jumps made by the particle as of time t. Observe that the particle will spend only a fraction of its time in active transport and hence the lateral speed of the particle should be slower than before.
Proposition 3.2. Let k 1 , k 2 > 0, and let r = v/δ > 0, then the mean velocity of the particle Proof. Consider the time τ for a particle to make one step of lateral transport. Before doing so a particle performs m switches from on-transport to off-transport and back, where m is distributed as a Geometric( r k 1 +r ) variable. The amounts of time a particle takes for each switch τ Q→P i from on-transport to off-transport are iid variables with Exponential(k 1 +r) distribution. Likewise, the amounts of time a particle takes for each switch τ P →Q i from off-transport back to on-transport are iid variables with Exponential(k 2 ) distribution, and are independent of the other switching times. Then with E[τ ] = (k 1 + k 2 )/(k 2 r). Let n t be the number of times a particle makes a step of lateral transport until time t. By the Renewal Theorem , and the result follows from X t = δn t and δ k 2 r k 1 +k 2 = k 2 v k 1 +k 2 .
Proposition 3.3. Let k 1 , k 2 > 0, and r = v/δ > 0. In the limit as δ → 0 the position of a particle satisfies Proof. The number of lateral transport steps (n t ) t≥0 is a renewal chain with E[τ ] = (k 1 + k 2 )/k 2 r and Var[τ ] = ((k 1 + k 2 ) 2 + 2k 1 r)/r 2 k 2 2 . By the LLN and CLT for Renewal chains and now the result follows from X t = δn t , δ 1 The fact that an individual particle starting at x = 0 will have the distribution given by Proposition 3.3 at long times was also obtained by Brooks [4] who proved tail estimates on the CLT and found that the remainder is O( 1 √ t ). In our approach, one can obtain the same error estimate using the asymptotic analysis of renewal chains, see, for example, [7].

Connection to Partial Differential Equations Models.
In order to demonstrate the connection between our model and the PDEs seen in [27][26] [14], we make a formal calculation regarding convergence of the law of the location of a single particle in the full system.
Let δ N = 1 N and let X N (t) denote the lateral position of a particle in a system with N boxes. Define q N (x, t) := P{X N (t) = x, and the particle is on track} p N (x, t) := P{X N (t) = x, and the particle is off track.} We suppress N in the notation. By definition of the Poisson process, for small values of h we have which we may reorganize as Formally, we take the Taylor expansion of q in space: . Taking limits as h → 0 and δ → 0 yields the system of PDEs When k 1 = 0 = k 2 , the limiting PDE is simple linear transport: (∂ t + v∂ x )q(x, t) = 0. The initial condition q(x, 0) = δ 0 (x) corresponds to the density of a single particle at the origin at t = 0. The time evolution via simple linear transport is translation of the delta function, while the time evolution via the equations (1) and (2) will have a spreading profile. This is seen in Propositions 3.1 and 3.3 since the variance of the fluctuations goes to 0 in Proposition 3.1 but not in Proposition 3.3 as δ → 0.
In the experiments described in the introduction one sees "approximate" traveling waves of radioactivity in the axons in the sense that there is a slowly spreading wave front moving at constant velocity away from the cell body. The equations (3) and (4) are linear and do not have solutions that are bounded traveling waves. However, it was shown by a perturbation theory argument in [27][26] that as ε → 0 the solutions of are to leading order This asymptotic form is valid for small ε, that is for large k 1 and k 2 . However, if we set q(x, t) = q ε ( x ε , t ε ) and p(x, t) = p ε ( x ε , t ε ), then q and p satisfy (3) and (4), so the solutions of (3) and (4) behave like approximate traveling waves for large t and large x whether or not k 1 and k 2 are large. The results suggested by these perturbation theory arguments have been proven rigorously by Friedman and coworkers [12][13] [14] [15].
In the case that radiolabeled particles enter the axon only for a short time, 0 ≤ t ≤ T , the same perturbation analysis shows that for large times the solution of (3) and (4), where q(0, t) = 0 for t > T is to leading order a spreading Gaussian with mean and variance as given in Proposition 3.3. This corresponds to what is seen experimentally and arises asymptotically in the solutions of (3) and (4) by convolving the asymptotic solution corresponding to the initial condition q(x, 0) = δ 0 (x) with the indicator function χ [0,T ] (x).
4 Dynamics from the spatial system perspective 4.1 The spatial system in equilibrium We are now ready to characterize the steady state dynamics induced by continually adding particles from the nucleus and removing them when they reach the distal end of the cell.
Proposition 4.1. If the incoming rate of particles at Q 0 is q 0 r, then the system has the following stationary distribution Proof. Given the value of the input rate Q 0 = q * the generator of the process (Q 1 , P 1 ) in the first section is If we first use f (q, p) = Q 1 (t) then use f (q, p) = P 1 (t) and take expectations we get a system of ODE's governing the change in E[Q 1 ], E[P 1 ] over time indicating that in equilibrium E[ For purposes that will soon be clear we first assume that Q 0 has a constant Pois(λ * ) distribution over time rather than just being equal to λ * . Let π(q * , q, p) = π λ * (q * ) ⊗ π λ Q (q) ⊗ π λ P (p) where π λ are Pois(λ) distributions with rates λ * , λ Q = λ * , and λ P = λ * k 1 k 2 respectively. To show that π(q * , q, p) is a stationary distribution for the process (Q 1 , P 1 ) we need to check that A q * f (q, p)π(q * , q, p) = 0 for any choice of function f ∈ D(A q * ).
In the inner two sums, for a fixed value of q * , the factor multiplying f (q, p) for any q, p ∈ N 2 comes only from terms involving {q − 1, q, q + 1} and {p − 1, p, p + 1} and equals e −(λ * +λ Q +λ P ) times Summing over q * the factor multiplying f (q, p) becomes We've just shown that if the distribution of input particles Q 0 in equilibrium is Pois(λ * ) then the stationary distribution of Q 1 , P 1 is independent of the distribution of Q 0 and is π λ Q ⊗ π λ P where λ Q = E[Q 0 ], and λ P = E[Q 0 ]k 1 /k 2 . It is easy to see in the above calculation that if the input Q 0 was indeed constant q 0 then the stationary distribution of (P 1 , Q 1 ) would again be π λ Q ⊗ π λ P with rates λ Q = q 0 and λ P = q 0 k 1 /k 2 respectively. Now, since the stationary distribution of Q 1 is the input into (Q 2 , P 2 ) process, the more general calculation shows that in equilibrium the stationary distribution of (Q 2 , P 2 ) is independent of the stationary distribution for Q 1 and is π λ Q ⊗π λ P with rates λ Q = E[Q 1 ] = q 0 and λ P = E[Q 1 ]k 1 /k 2 = q 0 k 1 /k 2 again. It follows by induction that the stationary distributions for {(Q i , P i )} are independent and identically distributed as π λ Q ⊗ π λ P , λ Q = q 0 , λ P = q 0 k 1 /k 2 . This is also an example of a clustering process satisfying the detailed balance conditions with linear rates discussed in Sec. 8.2 of [19].

Homogeneity of the axons at equilibrium
Recall that δ = 10nm, roughly the step size of motor proteins, and that axons can be up to one meter in length. Thus we are interested in phenomena on all the length scales 10 ν δ, where ν = 1, 2, . . . 8. Let ∆ = 10 ν δ; we want to determine how similar different segments of the axon of size ∆ are. Let Q ∆ and P ∆ denote the numbers of on-track and off-track particles in a segment of length ∆.
In equilibrium, Q ∆ and P ∆ are both sums of 10 ν independent Poisson random variables with parameters λ Q = q 0 and λ P = k 1 k 2 q 0 , respectively. Therefore the distributions of Q ∆ and P ∆ are Poisson with parameters 10 ν λ Q and 10 ν λ P , respectively. The mean and the variance of the number of particles in the segment of length ∆ is 10 ν (λ Q + λ P ). To see how homogeneous different slices of length ∆ are, we consider the coefficient of variation, c ∆ , which is the standard deviation divided by the mean.
As indicated in Section 2, q 0 is in the range 0.25 to 2.5 in different axons. For illustrative purposes here, we will assume q 0 = 1. Since k 1 /k 2 = 3, we see that the scale-dependent coefficient of variation c ∆ = 1/(2 √ 10 8 ∆). Therefore at the ten nanometer scale the coefficient of variation is simply 1/2. At the micron scale c ∆ = 1/20 and at the millimeter scale c ∆ = 0.5 × 10 −5/2 . The cutoff between "high variance" and "low variance" distributions is usually considered to be when the coefficient of variation is near 1, so by this standard the axon is extremely homogeneous in its length at large scales.

Approaching Equilibrium
We have seen above that the axon is very homogeneous at stochastic equilibrium on a space scale down to micrometers. One of the beautiful properties of the system of transport with reversible binding is that if it is locally out of equilibrium, the dynamics will automatically take it back to equilibrium. This is of fundamental importance to the biological function of the system because it means that the axon will automatically "repair" itself without central control of the repair process. How good is this mechanism? If a segment of the axon is far away from equilibrium, how long does it take to get back close to equilibrium? We study this question first for a single location and then use the estimates derived to scale the results to segments of any length.
Proposition 5.1. Let (Q, P )(0) = (0, 0) and let the constants a > 0 and p ∈ (0, 1) satisfy the relationship a 2 p|λ ∞ | > 1 where λ ∞ is the equilibrium vector of (Q, P )(t). Then there exists a t * such that t ≥ t * , In fact, the choice Proof. We begin by noting that for any given β ∈ (0, 1), we may choose t * > 0 such that for all t ≥ t * , the vector of means λ(t) Applying Chebyshev's Inequality, and observing that the variance of a Poisson random variable is equal to its mean, we conclude that Since the initial condition for both P and Q are less than their respective equilibrium values, each are monotonically increasing functions and the above reduces to for all t > t * . In order to satisfy the requirement that the right hand side must be less than p, we solve for β and find It remains to study the convergence of the mean and the appropriate choice of t * . The dynamics of the mean vector λ(t) are given by the ODĖ where e 1 is the unit vector (1, 0) and The solution to (6) is where λ ∞ = q 0 rA −1 1 e 1 = q 0 (1, rk 1 k 2 ). This yields the estimate where α is the smaller of the eigenvalues of A 1 , namely Noting that α > 0, t * may be chosen so that e −αt * ≤ aβ, i.e.
We now suppose that the whole axon is in statistical equilibrium except for a segment of length δ × 10 ν m in which we will assume that there are no particles either on or off the track. We want to know how long will it take for this segment to get back close to equilibrium. Of course, Proposition 5.1 covered the case ν = 0. We are interested in ν = 1, . . . , 8. We imagine that the axon is broken up into 10 8−ν segments of length δ × 10 ν m. In this rescaled system, the unbinding and binding rates per particle, k 1 and k 2 remain the same, as well as the mean on-transport velocity v. In order to retain this mean velocity, the rate of lateral stepping must be decreased tõ r = r × 10 −ν .
The ODE for the mean vector of the rescaled system is given by We note that the last term in (7) contains an r rather than anr. This is because the input rate is unchanged while the exit rate is diminished. The resulting equilibrium value is therefore rescaled as well, Both components of this vector are of order 10 ν , as expected.
We will use the parameters discussed in Section 2: k 1 = 1, k 2 = 1 3 , v = 10 −6 m/s, r = 10 2 s −1 . This impliesr = 10 2−ν s −1 . As in our analysis in Proposition 5.1, the time to equilibrium,t * , is proportional toα −1 whereα satisfies:α We are most interested in the cases ν ≥ 3 (the segment has length ≥ 10 microns). Then k 2r is small compared to k 1 . Ignoring constants and restricting our attention to the leading order terms gives Thus,t * ∼ 10 ν−2 s. Thus, for a 10 micron segment (ν = 3) the time to equilibrium is about 10 seconds and for a 1 millimeter segment (ν = 5) the time to equilibrium is 1000 seconds or 15 minutes. Of course, the time to get close to equilibrium depends on the parameter p that represents what we mean by "close."

Transport with Deposition
One major feature of the biology we have ignored is that some particles in the off-transport state are actually deposed permanently in the membrane or are used for some other purpose. Such a phenomenon is easy to add to the model proposed above, but we must address a new qualitative feature of the results. Namely there is an exponential loss of material as we move toward the distal end. Figure 2: Fast Axonal Transport with Deposition In addition to deposition, we recognize that particles are only useful in the membrane for a random amount of time with a given half-life, after which they decompose and get transported back to the cell body. Neglecting this retrograde transport, we augment the previous model as follows: at each δlength lateral section let D i (t) be the number of particles deposed in section i, with N = {0, 1, 2, . . .} as the state space for D i . In addition to the previous dynamics we also define the following transition rates The constant τ in the decay rate is the average lifetime of a deposed particle.

Dynamics from the particle perspective
As before, we can use the renewal theorem to determine the average speed of particles that are not deposed, as well as the spreading of the wavefront. Proposition 6.1. Let k 1 , k 2 , k 3 > 0, τ < ∞ and let r = v/δ > 0, then the mean velocity of the particle that is not deposed into the membrane by time Moreover, in the limit as δ → 0 the position of a particle satisfies Remark 1. Notice the mean velocity of particles is higher than when there is no deposition A simple explanation for this is that particles that are still in lateral transport had returned to on-transport back from off-transport before managing to be deposed into the membrane, meaning that their time off-transport was shorter than it would have been if they did not have to beat the exponential clock for deposition.
Proof. When there is deposition then by time t each particle has either already deposed, or is still performing lateral transport. A particle that is deposed into the membrane has velocity zero from the time of its deposition on. A particle that can still be laterally transported must have avoided deposition until time t.
If we proceed along the lines of the proof without deposition, then the only difference now is that the amount of time a particle takes for each switch from off-transport back to on-transport are i.i.d. exponential variables with parameter k 2 +k 3 rather than k 2 . Hence the amount of time τ for a particle to make a step from one section to another has E[τ ] = (k 1 +k 2 +k 3 )/((k 2 +k 3 )r), and by the Renewal Theorem the number of steps a particle that can still be laterally transported so the result now follows by the Renewal LLN and CLT for n t and the fact that and 6.2 Dynamics from the spatial system perspective Once again we compute the stationary distribution of the flow-through system, this time with deposition. Proposition 6.2. Suppose the incoming rate of particles at Q 0 is q 0 r. Then the stationary distribution of the system is Proof. Since (Q i , P i ) do not depend on D i the proof for the stationary distribution of (Q i , P i ) follows along the same lines as in the case when there is no deposition. We first check the stationary distribution for the process (Q 1 , P 1 ) in the first section. For given value of Q 0 = q * the generator of the process (Q 1 , P 1 ) is given by The system of ODE's governing the change in Hence, in equilibrium A q * f (q, p)π(q * , q, p) = 0 can be verified exactly as before. Since now the mean of Q 1 differs from the mean of Q 0 by a multiplicative factor ρ < 1, and this is the mean of the input into Q 2 , it follows by induction that the sequence of rates for {Q i , P i } Poisson random variables decays geometrically in i by a factor ρ i . This immediately implies that in equilibrium the magnitude of the variances also decays geometrically along the line of progression of lateral transport.
To find the stationary distribution for D i note that if the number of particles at P i = p * the generator for D i is Corollary 6.3. When in equilibrium the mean and the variance of the number of particles along the axon decay exponentially, with a loss ratio and the same loss ratio holds for the off-transport and the deposed particles.
Remark 2. Notice that the rate of loss is depends on the ratio of reaction rates k 1 /k 2 , which we earlier estimated to be approximately 3, as well as Observing that for any small δ, the approximation (1 + δx The same ratio holds for the other two states as well

Homogeneity despite deposition
It remains to estimate the rates k 3 and τ from experimental data and then make a formal calculation to address the volume of material loss due to deposition and turnover as a function of length scale. To be precise, suppose we want to know the percentage loss over a length ∆ = 10 ν δ. Then from the Corollary 6.3 the fraction of material loss is approximately So, for example, at the meter scale, ν = 8, and so the loss ratio is exp(−k 1 k 3 /v(k 2 + k 3 )). We now consider the parameters of a specific example, sodium pumps being deposed in the cell membrane. We assume that the half life of a single pump is approximately a week, which means τ = 6 × 10 5 s. It remains to estimate k 3 . According to [36], there are roughly 1000 sodium pumps per micron 2 of membrane surface area. If we take the diameter of the neuron to be 5 microns, then the membrane surface area for one 10 nm slice is about 0.15 µm 2 . So there are roughly 150 sodium pumps per slice.
It is important to note that each vesicle carries roughly 100 sodium pumps and so the stationary distribution for the quantity D i is on average 1/100th of the number of sodium channels which are deposed into the membrane at the location i. Taking this into account and recalling our earlier parameter estimations, k 1 = 1, k 2 = 1/3, v = 10 −6 and q 0 = 1 we estimate the average number of sodium pumps in the membrane of the first 10 nm slice to satisfy 150 = 100 q 0 k 1 k 3 τ /(k 2 + k 3 ) which implies k 3 = 1 × 10 −6 .
From the above we see that the loss ratio at the meter scale is which is to say that over the length of an entire meter long axon, we expect the concentration of sodium pumps in the membrane near the distal end to be roughly a 1/25 of that found near the nucleus. It is striking to note that over a length 0.1 m, however, the loss ratio is e −.3 = .74. This implies that material loss due to deposition and turnover becomes an issue exactly over the range of lengths of typical axons: 0.1 m to 1 m. At the lower end of this range, the axon is very homogeneous. At the higher end, one predicts significant inhomogeniety.

Discussion
In this paper, we used a spatial Markov chain model to unify several previous approaches to modeling fast axonal transport. After inferring the orders of magnitude of the parameters of the model, we have used standard techniques from renewal theory and linear Markov chain theory to give a precise description of the most striking qualitative feature of axonal transport: homogeneity along the length of the axon. Furthermore we predict the timescale of recovery to this homogeneous state after removing the vesicles from a segment as a function of the length scale of the segment. In creating the model, we have become aware of two features of the biology which have not been adequately studied: 1) Deposition and decay of proteins; and 2) characterization of the process by which unbound vesicles reattach to the transport mechanism.
Regarding deposition and decay, we mentioned in Section 6 that the model unequivocally predicts a loss of material down the length of the axon. Furthermore, it seems that certain biological features that we have chosen to ignore -such as retrograde transport and a location-dependent unbinding rate [9] -would only serve to increase the non-homogeneity. In the case of sodium pumps, it is not immediately clear if there is a contradiction between deposition loss and homogeneity, but we do see a sensitive dependence on the length of the axon. Furthermore, it stands to reason that for other types of cargo the loss may significant. We believe that this raises an interesting biological issue, and mention a few possible mechanisms for overcoming the material loss: first, it is possible that proteins can be synthesized in axons [33]; second, there may be some signal that prevents material from deposing in a region that is already "full;" and third, it may be that for some materials "deposition" is a reversible process.
Regarding reattachment, we note that we have modeled all event wait times as exponential random variables, but this is certainly a simplification. As an example, the stepping process of kinesin is a well-studied though still not completely understood phenomenon. Much work has focused on assessing the dependence of the mean rate of translocation on both the load and the local concentration of ATP [35] [29]. Implicit in this analysis is the assumption of exponential wait times with state dependent rate parameters. However, when fitting to data and matching dispersion information the authors in [11] found it necessary to generalize the wait time distribution. This was followed by more detailed models for which it was shown that load carrying could in fact regularize the stepping times of kinesin motors [28] [8].
Generalizing the wait times does not significantly affect the particle perspective results of Sections 3.2 and 6.1. This is because renewal theory, which depends on independence of the wait times, is robust with respect to such changes. However, the spatial system results may be affected in that the full process is no longer Markov. While one may still expect a something like a stationary distribution for the flow-through system (where new particles enter from the nucleus and particles are removed from the system at the distal end), the analysis of the approach to equilibrium may change.
In the absence of direct observation of the phenomenon, we refrain from this more detailed analysis.
In light of this known need for generalized wait times in the stepping process, it seems likely that detailed observation of the rebinding process will call for new models as well. Recall that when a vesicle unbinds from a microtubule it is unclear whether it typically rebinds to the same microtubule or if it explores the region significantly via diffusion before finding a different microtubule to bind to. In the latter case, a more appropriate model for rebinding time would be to solve some kind of first passage time problem and use that distribution for the rebinding wait.  Abstract. In this paper we develop a probabilistic micro-scale compartmental model and use it to study macro-scale properties of axonal transport, the process by which intracellular cargo is moved in the axons of neurons. By directly modeling the smallest scale interactions, we can use recent microscopic experimental observations to infer all the parameters of the model. Then, using techniques from probability theory, we compute asymptotic limits of the stochastic behavior of individual motor-cargo complexes, while also characterizing both equilibrium and non-equilibrium ensemble behavior. We use these results in order to investigate three important biological questions: (1) How homogeneous are axons at stochastic equilibrium? (2) How quickly can axons return to stochastic equilibrium after large local perturbations? (3) How is our understanding of delivery time to a depleted target region changed by taking the whole cell point-of-view?
1. Introduction. In all cells, one finds that proteins, membrane-bound organelles, and other structures (e.g. chromosomes) are transported from place to place at speeds much higher than diffusion. Though these transport processes are fundamental to cell function, many of the underlying mechanisms, organizational principles, and regulatory features remain unknown. Axonal transport is one of the best studied systems because the transport is basically one-dimensional since axons are long and narrow. There are two speeds of axonal transport. Fast transport goes at speeds of roughly 0.2 to 0.5 meters/day [27] [33], while slow transport goes at approximately 1 millimeter/day, the rate of axon growth and regeneration [6] [27]. The biology and principles of slow transport are not yet clear [6], but the basic mechanisms of fast axonal transport were discovered in the 1980s [2][3][29] [43]. The model in this paper refers to fast axonal transport, which we will henceforth call axonal transport.
The axonal transport apparatus consists of vesicles which form reversible chemical bonds with motor proteins that bind reversibly to microtubules which run parallel to the long dimension of the axon [1]. When the vesicle-motor protein complex is assembled on the microtubule, the complex steps stochastically with step size approximately 8 nanometers for kinesin and dynein and 10 nanometers for myosin [7] [12][18] [40]. The vesicles enter from the cell body on microtubules and then detach and reattach to the transport mechanism at random times.
In this paper we propose a spatial Markov-chain compartmental model based on these dynamics. We will assume independence of the interactions, and exponential wait times between events. While we address the validity of these assumptions in the Discussion section, we consider this a useful "first-order" approximation that permits study of the dynamics from both the perspective of individual vesicles as well as that of the full spatial system. Such a model unifies all earlier deterministic and stochastic modeling efforts and can accommodate both qualitative and quantitative experimental data observed on multiple scales.
In much experimental work in the 1970s and 1980s, radio-labeled amino acids were put into the cell bodies continuously or for a few hours. The amino acids were incorporated into proteins that were packaged into vesicles and put on the transport system so that at later times radioactivity could be seen moving progressively down the axons. In the continuous infusion case, one would see a wave of radioactivity with a sharp but slowly spreading wavefront propagating at constant velocity down the axon. In the case of infusion for a few hours one would see at long times a slowly spreading pulse of radioactivity that looked normally distributed. It was to understand this behavior that Reed and Blum constructed PDE models for axonal transport [3][34] [35]. These models did not have traveling wave solutions, but the data certainly looked like approximate traveling waves. In [36] it was shown by a perturbation theory argument that, in the asymptotic limit where the unbinding and binding rates k 2 and k 1 get large, the solution approaches a slowly spreading traveling wave or a normal pulse. Recently, in a series of papers, Friedman and co-workers have introduced new PDE models and proved these results rigorously [14][15][16] [17].
Probabilistic models for axonal transport were introduced and used for simulations already in the 1980s [39] [41]. However, rigorous work began with Lawler [28] in 1995 and was continued by Brooks [5] who used a continuous time stochastic model to show that the distribution of an individual particle is a spreading Gaussian at large times. Brooks also proved tail estimates for the central limit theorem and used them to estimate the error from normal. Independently, Bressloff [4] developed a discrete stepping model and performed an analysis under the assumption that the rate of unbinding and binding to transport is fast relative to lateral velocity over the length scale of interest. The author derived a characterization of the spreading wavefront of a particle entering at the nucleus and traveling to the distal end. This model served as the basis for later investigations by Newby and Bressloff [31,32] wherein the authors characterize the axonal transport system as an intermittent search for hidden targets.
In this paper we revise the existing probabilistic models in order to study randomness in the system as a whole rather than exclusively from the point of view of an individual particle. Our goal is not only to recover and generalize previous results, but also to investigate three specific, biologically important, properties of the whole stochastic system.

Summary of Results.
In Section 2, we create a continuous-time Markov chain queueing model for the axonal transport system. We show how to use experimental data to determine (or estimate) all the parameters of the model.
In Section 3, we take the individual vesicle point-of-view. We prove the asymptotic forms in [36] with rigorous error estimates. We show that in the limit as the compartment size becomes small our model becomes the probabilistic model of [5]. We also show that in the limit as the length of the axon and time become large (with the scale of axon length on the order of the squared scale of time) our model becomes the PDE model of [35]. Since we assume that particles are independent, the time evolution of the law of an individual will reflect the behavior of an ensemble of particles released at the same time.
In Section 4 we adopt the full spatial system perspective to quantify stochasticity along the length of the axon. We begin by calculating in Proposition 4.1 the stationary distribution of a flow-through system that has sustained input from the nucleus, while particles are removed upon reaching the distal end. The stationary distribution has a product Poisson structure which allows for seamless transition between spatial scales. With this mathematical model, we are able to make precise statements about three biologically important system properties. In the stationary distribution, the number of vesicles in each slice of the axon is independent and identically distributed (Proposition 4.1), however this in and of itself is not sufficient to account for the sense that samples taken for different parts of the cell "look the same." We compute in Section 4.2 the coefficient of variation for the number of vesicles in sections of different length and show that the coefficient of variation is low for all but the smallest length scales. In Section 4.3, we study the intermittent search problem posed by Newby and Bressloff [31] from the system point of view. Efficient transport to locations that need material must balance the speed of transport of material from the nucleus to the distal end of the cell with the rates of dissociation from the transport apparatus. We calculate the expected hitting time for a hidden target by all vesicles in the system. In so doing we encounter the counterintuitive result that while increasing the velocity of the motors while on transport increases the chance of any particular vesicle missing the target, the expected hitting time by the system actually decreases.
This hitting time approach is natural for needed material that is sparsely distributed throughout the axon, but when the needed cargo in question is more common, the time to replenishment is better addressed through the ODE approach that we develop in Section 4.4. Due to the product structure of the law of the transient dynamics, this non-equilibrium behavior is determined by the 2N -dimensional ODE governing the means. From this we estimate the timescale of return to equilibrium as a function of the length scale of interest.

The model and its parameters.
Let L be the length of the axon, divided evenly into N = L/δ lateral sections each of length δ, equal to the step size of the motor protein. Within each section, we disregard any further spatial geometry and take the particles to be in one of two states: · an on-transport state that steps laterally at a rate r = v/δ per section, or · an off-transport state that does not step laterally. We use a 2N -dimensional continuous time Markov chain {(Q i (t), P i (t)), i = 1, . . . , N } to model the particle dynamics, where · Q i (t) is the number of particles at time t in on-transport state in section i, · P i (t) is the number of particles at time t in off-transport state in section i. · (Lateral transport) ; · (Switch from on-transport to off-transport) (Q i , P i ) → (Q i − 1, P i + 1) at rate k 2 Q i (t); · (Switch from off-transport to on-transport) . The lateral transport rate, r = v/δ, is inversely proportional to the length scale so that the mean number of particles per unit length is invariant with respect to rescaling δ. We will assume that the rate of production q 0 = δρ 0 , for some constant ρ 0 > 0, in order for the mean number of particles in each compartment to scale with the size δ of each compartment. This will imply that the mean number of particles per unit length scales as ρ 0 . A graph of the model is depicted in Figure 1.
In order to insure the Markov property, we use exponential random variables for the waiting times between transition events. Specifically, we mean that after a given On-and-off Transport Chain event we assign a new independent random variable to each of the 3N +2 possible next events, exponentially distributed with the appropriate rate parameter. The system of values updates according to the transition associated with the minimum of these waiting times. Then we create a new set of exponential random variables and the process proceeds as before.
The advantage of computing explicit formulas for quantities that can be observed in experiments is that the experimental data can then be used to determine the parameter values in the model. For the characterization of the approximate wavefront speed and spreading in Section 3 and the homogeneity calculations in Section 4 we need order-of-magnitude estimates for the parameters. Actual parameter values will certainly differ depending on the particular neural tissue and the particular particles being transported. However, we can get order-of-magnitude estimates from existing data.
First we recall that fast transport has been observed to travel at speeds of 0.2 to 0.5 m per day. We can assume that the average velocity of particles while physically bound to microtubules is roughly 1 m per day, or v = 10 −6 m/s. We assume that the compartment size scales as the length scale of the individual steps of the motor protein, so δ ∼ 10 −8 m. This implies that the rate parameter should be r = v/δ = 100s −1 .
We now turn our attention to the on-off rates rates k 2 and k 1 . These can be determined from experimentally observed run lengths on the transport system. Indeed, Dixit et al. [10] show that a typical run along microtubules for dinein and kinesin is on the order of 10 −6 m. We can compare this with the theoretical run length of the model to determine off-rate k 2 . Within the model, at each step on the transport mechanism the particle has a binary decision to jump laterally along the transport with probability r/(k 2 + r), or to jump off with probability k 2 /(r + k 2 ). The number of jumps along the transport system before jumping off is therefore geometrically distributed on the set {0, 1, . . .} with success probability r r+k2 . It follows that average number of steps in the run is r k2 and therefore the average run length is r k2 × 10 −8 m. Setting this equal to the average experimental run length of 10 −6 from [10], we see that r k2 ∼ 100, implying that k 2 ∼ 1s −1 . As we will see in the computation of the stationary distribution in Section 4.1, the ratio of the expected number of particles on the track to those off the track is k2 k1 . Dixit [10] found that approximately 75% of the particles were motile so this ratio is approximately equal to 3. Since k 2 ∼ 1s −1 we see that k 1 ∼ 1 3 . It remains to estimate q 0 . We will see in Proposition 4.1 that the mean number of particles per compartment is (1 + k2 k1 )q 0 = 4q 0 . Of course, axons have a large variety of diameters and larger axons will have more vesicles per unit length so one expects a range of values for q 0 . However, examination of a large number of electron micrographs of axonal cross-sections (see for example [20], Fig. 3; [21]; [30]), which are typically 100 nm thick enables one to estimate the number of vesicles per 100 nm segment. This number is typically in the range of 10 to 100 which implies that there are 1 to 10 vesicles per compartment in our model. Therefore q 0 is in the range 0.25 to 2.5, for various axons.
We remark that we are ignoring some aspects of the physics and the biology of axonal transport. We are not including diffusion of the vesicles off the track. We are treating the microtubule track as though it were a single continuous entity from one end of the axon to the other, when it fact it consists of numerous, separated, microtubule fragments. And, we are ignoring retrograde transport and the details of the motor proteins. Nevertheless, this simple model will enable us to investigate the homogeneity questions that are the main goal of this paper.
3. Dynamics from the Particle Perspective. In this section we calculate properties of the stochastic dynamics by using stochastic convergence theorems and stochastic averaging theorems from probability theory. We first see that, in the δ → 0 limit, the law of the location of a single particle corresponds to that of a particle with a piecewise linear Markov motion. We then show this law can be approximated by the Green's function of a linear partial differential equation. This enables us to obtain, as a special case, the asymptotic behavior of the PDE models for axonal transport in an either rate limiting or perturbed setting.
3.1. The active transport mode. We first consider the simple case where the particle starts at X δ 0 = 0 and stays exclusively in active transport mode. Let X δ t ∈ {0, δ, 2δ, . . . , L} be the lateral position of a particle at time t and let n δ t be the number of jumps made by the particle as of time t. Observe, X δ = δn δ . Lemma 3.1. Let k 2 = k 1 = 0, and r = v/δ > 0, then the position of the particle satisfies, for any time t < ∞ sup s≤t X δ s − vs −→ δ→0 0 a.s. and, for B a standard Brownian motion in distribution on the Skorokhod space of cadlag (right continuous left limited) functions.
Proof. Since n δ is a Poisson process with rate r = v/δ, defining N t := n δ δt we get a Poisson process N with rate v, and we have X δ · = δN ·/δ . Our results then follow directly from the functional law of large numbers (FLLN) and the functional central limit theorem (FCLT) for the Poisson process N .
3.2. The on/off dynamics. We now consider a particle which undergoes transitions from on-transport to off-transport state and back. Denote again by X δ t ∈ {0, δ, 2δ, . . . , L} the lateral position of a particle at time t and let n δ t be the number of lateral transition jumps made by the particle as of time t. Observe that the particle will spend only a fraction of its time in active transport and hence the lateral speed of the particle should be slower than before.
A non-compartmental stochastic model for axonal transport, introduced by Brooks in [5], is as follows. A particle can be in one of two states: · an on-transport state with deterministic lateral velocity v, or · an off-transport state with lateral velocity 0. We use a 2-dimensional Markov process to model the particle dynamics, where · X t be the lateral position of this particle at time t, · ξ t be the indicator for whether it is on (1) or off (0) transport at time t. Definition 3.2. (Stochastic non-compartmental model) Let (X t , ξ t ) t≥0 be a piecewise-linear Markov process with values in (R + , {0, 1}) started at (X 0 , ξ 0 ) = (0, 1) with the following dynamics: · (Switch from on-transport to off-transport) (X t , 1) → (X t , 0) at rate k 1 · (Switch from off-transport to on-transport) (X t , 0) → (X t , 1) rate k 2 · (Lateral travel) X t = t 0 vξ s ds The path of (X t ) t≥0 consists of alternating sequence of Exponential(k 2 ) stretches of time where the lateral position increases linearly with speed v, and Exponential(k 1 ) stretches of time where it remains constant. Proposition 3.3. Let k 2 , k 1 > 0, and r = v/δ > 0, then the position of the particle converges in distribution on the Skorokhod space of cadlag functions.
Proof. If we let ξ δ be the indicator of whether the particle in the compartmental model is on (ξ δ = 1) or off (ξ δ = 0) transport, then (X δ t , ξ δ t ) t≥0 is a strong Markov process. We will see that ξ δ is a continuous-time Markov chain on {0, 1} (independent of δ), and conditionally on ξ δ the transition law of X is easily expressed. Likewise, for the non-compartmental model above, (X t , ξ t ) t≥0 is a strong Markov process, with ξ the same continuous-time Markov chain on {0, 1} as ξ δ , and conditionally on ξ the change in X is easily given in terms of its linear speed and ξ.
We will start by showing that for any t > 0, (X δ t , ξ δ t ) converges to (X t , ξ t ) in distribution as δ → 0. We then show that the finite dimensional distributions of (X δ , ξ δ ) converge to those of (X, ξ). A tightness argument finally implies (Lemma 16.2 and Theorem 16.3 [22]) that (X δ , ξ δ ) converges to (X, ξ) in distribution on the Skorokhod space of cadlag processes. Suppose that initially the particle is on transport at x 0 , so (X δ 0 , ξ δ 0 ) = (x 0 , 1). The first time τ 1 = inf{t > 0 : ξ δ t = 0} at which the particle steps off transport has Exponential(k 2 ) distribution, irrespective of δ. The first subsequent increment in time σ 1 = inf{t > 0 : ξ δ τ1+t = 0} after which the particle steps back on transport has Exponential(k 1 ) distribution, irrespective of δ as well. This is repeated, and ξ δ is a simple continuous-time Markov chain on {0, 1} with transition rates k 2 and k 1 , from 1 → 0 and 0 → 1, respectively.
Until time τ 1 the particle behaves as if it were in active transport (k 2 = k 1 = 0) and conditionally on the value of τ 1 , for any 0 ≤ s ≤ τ 1 the lateral change in position over time s, X δ s − X δ 0 , is δn δ s where n δ s is a Poisson(rs) random variable. Moreover, by results of Proposition 3.1, conditionally on the value of τ 1 , we have To get the unconditioned law of X δ τ1 − X δ 0 , we observe that the number of boxes the particle traverses before it steps off n δ τ1 has a Geometric(k 2 /(k 2 + r)) distribution (note that a Poisson rate r process sampled at an Exponential(k 2 ) time independent of the process has this distribution). Since k 2 /δ(k 2 + r) → k 2 /v, as δ → 0, X δ τ1 −X δ 0 = δn δ τ1 converges in distribution to Exponential(k 2 /v) variable. Consider now the particle in the non-compartmental model. It is immediate from the definition of the model that ξ has the same law of a continuous-time Markov chain on {0, 1} with transition rates k 2 and k 1 , from 1 → 0 and 0 → 1, as ξ δ . Since we prove our convergence in law results by first conditioning on the values of ξ δ and ξ, we can without loss of generality henceforth assume ξ δ = ξ a.s., and drop its superscript. Suppose that initially the particle in the non-compartmental model is on transport at x 0 , so (X 0 , ξ 0 ) = (x 0 , 1). Conditionally on τ 1 , for all 0 ≤ s ≤ τ 1 , X s − X 0 = vs and X τ1 − X 0 = vτ 1 , hence unconditionally, X τ1 − X 0 is an Exponential(k 2 /v) random variable. We now have both sup 0≤s≤τ1 (X δ s − X δ 0 ) − (X s − X 0 ) → 0 a.s. and X δ τ1 − X δ 0 ⇒ X τ1 − X 0 . Between times τ 1 and σ 1 the particle in both models stays in place, so conditionally on values of τ 1 , σ 1 , and of X δ τ1 , X τ1 , sup τ1≤s≤τ1+σ1 (X δ s − X δ τ1 ) − (X s − X τ1 ) ≡ 0, and X δ τ1+σ1 −X δ 0 ⇒ X τ1+σ1 −X 0 . At time τ 1 +σ 1 , the same process starts over from initial values (X δ τ1 , 1) and (X τ1 , 1) in the compartmental and non-compartmental model, respectively.
Let τ 1 , σ 1 , τ 2 , . . . , be the sequence of time increments between the consecutive times when the particle in both models gets off and gets back on transport, let σ 0 = 0 and for i ≥ 1 Then (τ i ) i≥1 and (σ i ) i≥1 are independent sequences of i.i.d. Exponential(k 2 ) and Exponential(k 1 ) variables, respectively. For any t > 0, let η t be the number of times the particle in either model gets back on transport until time t, and η ′ t the number of times it gets off, Note that, η ′ t = η t iff ξ t = 1, and η ′ t = η t + 1 iff ξ t = 0. Letτ t be the last time before time t that the particle changed whether it was on or off transport, that is, τ t = sup{0 ≤ s ≤ t : ξ s − = ξ s }. Then, we have If η ′ t = η t , then from time τ t to t the particle is in active transport, and the same convergence argument as before implies that conditionally on the values of τ t and X δ τt , Also X δ τt − X δ 0 = δn δ τt where n δ τt is the number of boxes the particle traverses by time t. Conditionally on the value of η t , n δ τt is a sum of η t i.i.d. Geometric(k 2 /(k 2 + r)) random variables, hence X δ τt − X δ 0 converges in distribution to a sum of η t i.i.d. Exponential(k 2 /v) random variables. In the non-compartmental model, if η ′ t = η t , then conditionally on the values of η t and τ t , for τ t ≤ s ≤ t, X s − X τt = v(s − τ t ) and conditionally only on the value of η t , X δ τt − X 0 is a sum of η t i.i.d. Exponential(k 2 /v) variables. Hence, conditionally on η t and τ t , sup τt≤s≤t (X δ s − X δ τt ) − (X s − X τt ) → 0 a.s., and conditionally only on η t , X δ τt − X δ 0 ⇒ X τt − X 0 . If η ′ t = η t + 1, then from time τ t to t the particle is in both models stays in place, so conditionally on η t and τ t , sup τt≤s≤t (X δ s − X δ τt ) − (X s − X τt ) = 0 a.s.. Also, conditionally only on the value of η ′ t , n δ τt is a sum of η ′ t i.i.d. Geometric(k 2 /(k 2 + r)) variables, hence X δ τt −X δ 0 converges in distribution to a sum of η ′ t i.i.d. Exponential(k 2 /v) variables. In the non-compartmental model, if η ′ t = η t + 1, conditionally only on the value of η ′ t , X δ τt is a sum of η ′ t i.i.d. Exponential(k 2 /v) variables. Hence, conditionally only on η ′ t , X δ τt − X δ 0 ⇒ X τt − X 0 . Now, integrating over the possible values of η t ,η ′ t and τ t , we get that for any t ≥ 0, (X δ t , ξ δ t ) ⇒ (X t , ξ t ). Convergence of finite dimensional distributions follows from an iterative use of the Markov property of (X δ , ξ δ ) and (X, ξ), and the fact that the increments of both (X δ , ξ δ ) and (X, ξ) are stationary.
The fact that an individual particle will have the distribution given by Proposition 3.3 as the size of the boxes decreases means that our model is a microscopic version of the stochastic model used by Brooks [5], and that the hydrodynamic limit of our model as δ → 0 is equal to the macroscopic stochastic model from [5]. An approximation of the particle's position X t is obtained in [5] to be X t ≈ µt + √ tσZ as t → ∞, where µ = k 1 v/(k 2 + k 1 ), σ = 2k 2 k 1 v 2 /(k 2 + k 1 ) 3 and Z is a standard Normal variable. That approximation is valid only for large fixed values of t, while we next extend this result to give an approximation for the whole time trajectory of the particle's path. This is accomplished by the following functional central limit theorem for the position of the particle undergoing stochastic transport.
Proposition 3.4. Let X be the position of a particle following the piecewise linear Markov process from Proposition 3.3 started at X 0 = 0 on transport, then and, if B denotes a standard Brownian motion in distribution on the space of continuous functions.
Once we rescale the position for the particle by 1/ √ n and time by 1/n, ξ still changes at a much faster rate than the position of the particle X. The generator of the rescaled centered process (X n , ξ n ) defined as t , ξ n t := ξ nt is the closure of the operator ). We will use the stochastic averaging theorem (Theorem 2.1 of [26]) to show that the paths of the centered rescaled process converge in distribution to paths of a Brownian motion with a diffusion coefficient equal to Let h(ξ) be the function Then f n → f as n → ∞ and whereĀ is defined on D(Ā) = C 2 0 (R) bȳ Define a sequence of processes Then our earlier calculation implies that for any f ∈ D(Ā) is a sequence of martingales. Since f ∈ C 2 0 (R), ξ n t ∈ {0, 1} it is clear that In order to apply Theorem 2.1 of [26] on stochastic averaging it is only left to show the process X n satisfies the compact containment condition, that is for any t > 0 and ∆ > 0 there exists a compact set K ⊂ R such that inf n P{X n s ∈ K ∀s ≤ t} ≥ 1 − ∆ This follows from the fact that X n t + h(ξ n t )/ √ n is a sequence of martingales (let and choosing C (given on t and ∆) so that the right hand side of the inequality with n = 1 is less than ∆, shows that with K = [−C − h max , C − h min ] the compact containment condition holds for (X n ) n≥1 . Now Theorem 2.1 of [26] implies that X n ⇒ W in distribution on the Skorkhod space of continuous functions, where W is a process with generatorĀ, and consequently has the same distribution as 2k 2 k 1 v 2 /(k 2 + k 1 ) 3 B where B is a standard Brownian motion.

Connection to Partial Differential Equations Models.
In order to demonstrate the connection between our model and the PDEs seen in [36][35] [16], consider the process (X, ξ) of the particle following the piecewise linear Markov process, and for any x ≥ 0, t ≥ 0 let q(x, t) = P{X t ∈ dx, ξ t = 1}/dx, p(x, t) = P{X t ∈ dx, ξ t = 0}/dx denote the probability densities of the particle's location x on and off transport, respectively, over time. Kolmogorov forward equations for (X, ξ) imply that q and p satisfy the system of PDEs (3.2) When k 2 = 0 = k 1 , the limiting PDE is simple linear transport: (∂ t +v∂ x )q(x, t) = 0. The initial condition q(x, 0) = δ 0 (x) corresponds to the density of a single particle at the origin at t = 0. The time evolution via simple linear transport is translation of the delta function, while the time evolution via the equations (3.1) and (3.2) will have a spreading profile. This is clear from the macroscopic limits of (X δ , ξ δ ) as δ → 0. When k 2 = k 1 = 0 the particle never switches off from traveling on transport at speed v and is deterministic, as seen in Proposition 3.1. When k 2 , k 1 > 0 the particle follows a truely stochastic process (X, ξ) with a non-zero variance, as seen in Proposition 3.3 and Proposition 3.4.
In the experiments described in the introduction one sees "approximate" traveling waves of radioactivity in the axons in the sense that there is a slowly spreading wave front moving at constant velocity away from the cell body. Equations (3.1) and (3.2) are linear and do not have solutions that are bounded traveling waves. It was shown by a perturbation theory argument in [36][35] that as ε → 0 the solutions of subject to q ε (0, t) = q 0 , are to leading order where H satisfies the heat equation This asymptotic form is valid for small ε, that is for large k 2 and k 1 . However, if we set q(x, t) = q ε ( x ε , t ε ) and p(x, t) = p ε ( x ε , t ε ), then q and p satisfy (3.3) and (3.4), so the solutions of (3.1) and (3.2) behave like approximate traveling waves for large t and large x whether or not k 2 and k 1 are large. These results have been proven rigorously by Friedman and coworkers [14][15] [16] [17].
4. Dynamics from the spatial system perspective.
4.1. The spatial system in equilibrium. We are now ready to characterize the steady state dynamics induced by continually adding particles from the nucleus and removing them when they reach the distal end of the cell.
Proposition 4.1. Let ({(Q i (t), P i (t)), i = 1, . . . , N }) t≥0 be the number of particles in the axonal transport system with compartments of size δ, on and off transport, respectively. Suppose the rate of production of particles from the source is rq 0 = vρ 0 . Then this Markov chain has the product-form stationary distribution Proof. Since the production rate is rq 0 , the generator of the process (Q 1 , P 1 ) is If we use f (q, p) = Q 1 (t) and f (q, p) = P 1 (t), and take expectations, we get a system of ODEs governing the change in indicating that in equilibrium in the first section the mean numbers of on-transport particles and off-transport particles are E[Q 1 ] = q 0 , and E[P 1 ] = E[Q 1 ] k2 k1 = q 0 k 2 /k 1 , respectively. Let π q0 (q, p) = π λQ (q)⊗ π λP (p) be a product of two independent Poisson distributions with rates λ Q = q 0 , and λ P = q 0 k 2 /k 1 respectively. To show that π q0 (q, p) is a stationary distribution for the process (Q 1 , P 1 ) we need to check that In the two sums the factor multiplying f (q, p) for any (q, p) ∈ N × N comes only from terms involving {q − 1, q, q + 1} and {p − 1, p, p + 1} and equals e −(λQ+λP ) times Thus, in equilibrium the input rate for (Q 2 , P 2 ), which is rQ 1 , has a Poisson distribution with mean rq 0 , and is independent of P 1 . Let π q0 (q 1 ) ⊗ π λQ (q 2 ) ⊗ π λP (p 2 ) be a product of three Poisson distributions with rates q 0 , λ Q = q 0 , and λ P = q 0 k 2 /k 1 respectively. To show that this is a stationary distribution for the process (Q 1 , Q 2 , P 2 ) we need to check that for any choice of function f ∈ D(A q1 ). For each fixed value q 1 , according to our previous calculation the inner two sums give 0, so the whole sum is 0.
Thus, in equilibrium, (Q 2 , P 2 ) have the distribution π λQ (q) ⊗ π λP (p) with λ Q = q 0 , λ P = q 0 k 2 /k 1 as well, and are independent from (Q 1 , P 1 ). It follows by induction that the stationary distributions for {(Q i , P i )} are independent and identically distributed as π λQ (q) ⊗ π λP (p), with λ Q = q 0 , λ P = q 0 k 2 /k 1 . This is also an example of a clustering process satisfying the detailed balance conditions with linear rates discussed in Sec. 8.2 of [25].
We point out that the mean number of particles both on and off transport is (1 + k2 k1 )q 0 = (1 + k2 k1 )ρ 0 δ, scaling with the size of a compartment. To obtain the mean number of particles per unit length we add particles in ≈ 1/δ compartments, and the mean number of particles per unit length is (1 + k2 k1 )ρ 0 independent of the choice of compartment size.
One immediate consequence is the analogous result for the number of particles in the stochastic non-compartmental model at any location along the axon. Namely, suppose the particles move according to the piecewise-linear Markov process (X, ξ) with a Poisson rate ρ 0 influx of new particles at location 0. Then, at any location 0 < x < L along the axon, the numbers of particles (Q (x,x+dx) , P (x,x+dx) ) on and off transport, respectively, have the stationary distribution Pois(ρ 0 dx) ⊗ Pois( k2 k1 ρ 0 dx) where for any x 1 , . . . x k ∈ (0, L), {Q xi } 1≤i≤k and {P xi } 1≤i≤k are mutually independent. We note that this result would not have been obvious to notice without going through the compartmental model first, yet its consequences for prediction and analysis of the long term stochastic behavior of the system are quite powerful.

4.2.
Homogeneity of the axons at equilibrium. Recall that δ = 10nm, roughly the step size of motor proteins, and that axons can be up to one meter in length. Thus we are interested in phenomena on all the length scales 10 ν δ, where ν = 0, 1, 2, . . . 8. Let ∆ = 10 ν δ; we want to determine how similar different segments of the axon of size ∆ are. Let Q ∆ and P ∆ denote the numbers of on-track and off-track particles in a segment of length ∆.
In equilibrium, Q ∆ and P ∆ are both sums of 10 ν independent Poisson random variables with parameters λ Q = q 0 and λ P = k2 k1 q 0 , respectively. Therefore the distributions of Q ∆ and P ∆ are Poisson with parameters 10 ν λ Q and 10 ν λ P , respectively. The mean and the variance of the number of particles in the segment of length ∆ is 10 ν (λ Q + λ P ). To see how homogeneous different slices of length ∆ are, we consider the coefficient of variation, c ∆ , which is the standard deviation divided by the mean.
As indicated in Section 2, q 0 is in the range 0.25 to 2.5 in different axons. For illustrative purposes here, we will assume q 0 = 1. Since k 2 /k 1 = 3, we see that the scaledependent coefficient of variation c ∆ = 1/(2 √ 10 8 ∆). Therefore at the ten nanometer scale the coefficient of variation is simply 1/2. At the micron scale c ∆ = 1/20 and at the millimeter scale c ∆ = 0.5 × 10 −5/2 . The cutoff between "high variance" and "low variance" distributions is usually considered to be when the coefficient of variation is near 1, so by this standard the axon is extremely homogeneous in its length at large scales.

4.3.
Balance between efficient transport and targeted delivery. The preceding characterization of the transport apparatus enables us to address questions concerning whole cell function. One core issue is that intracellular transport must simultaneously accommodate two functional demands: some material, such as the enzymes used to construct neurotransmitters, must be transported from the soma to the axon terminal in a timely manner; whereas other cargo, such as sodium channels, need to be delivered to unspecified locations as needs arise throughout the length of the cell. The tradeoff between these two goals is clear. If a typical vesicle spends the vast majority of its time in transport mode, the mean velocity will be close to the mean on-transport velocity, but any needs that arise in the central part of the cell will be neglected. On the other hand, if a typical vesicle spends too much time off transport, presumably available for use if needed locally, then it will take substantially longer to traverse the entire cell.
Recently Bressloff and Newby [31,32] modeled particles that are created near the nucleus that then undergo intermittent search (being in search mode while off transport and not in search mode when on transport) for a target hidden somewhere along the axon. Their model is a non-compartmental individual particle model with the additional feature that vesicles can move backwards as well as forwards. They compute the probability that the particle is successful and conditioning on success, the mean first hitting time. With our system-wide model we can accommodate the observation the if a given vesicle misses the target, another vesicle with similar cargo will pass by before too long. We will assess this hitting time under two assumptions about the density of relevant material. Our standing assumption q 0 ∼ 1 is appropriate for types of cargo that are found densely throughout the cell. In this setting, the wait time is essentially just the time it takes for one of several nearby vesicles to unbind from transport in the target region. A more interesting case is a setting where the needed cargo is sparsely distributed, say q 0 ∼ 10 −3 . In this setting, if the first cargo to reach the target region fails to unbind, there will be significant time before the next arrival. As we will see, we can still assess the trade-off intrinsic between risking a target miss and diminishing the time of the next arrival.
To make the discussion concrete we define a target region R n = {i * +1, . . . , i * +n} where i * ∈ 1, . . . , N − n. At time zero, we take the system {Q i (0), P i (0)} to be drawn from the stationary distribution described by Proposition 4.1 conditioned on the event that for all i ∈ R n , Q i (0) = P i (0) = 0. We introduce the hitting time H n := inf{t > 0 : i∈Rn P i (t) > 0}, which marks the first time a particle is off transport while in R n . Since computing the mean of H n is analytically intractable, we introduce another hitting time H ′ n , stochastically dominating H n , that nevertheless reflects the essential tradeoff between maximizing mean velocity and making detachment from transport likely in the target region.
Let I * := {i ∈ {1, . . . , i * } | Q i (0) + P i (0) > 0} be the set of all non-empty sections of the cell at time 0. Among the particle in these sections some will be "successful" in that they will detach from transport in the target region, while others will be unsuccessful. Labeling each particle in these terms, we decompose I 0 into locations with successful particles I s * and unsuccessful particles I u * . We are interested in the time H ′ n at which the rightmost successful particle, that is, a particle starting from position i m = max{I s * }, detaches while in R n . If S = 0 then we define i m = max{∅} := 0. For the position of this particle we use the notation (X δ t ) t≥0 from Section 3, where X δ 0 = δi m , X δ t ∈ {δi m , . . . , δN }, together with the indicator (ξ δ t ) t≥0 , ξ δ t ∈ {0, 1} of whether the particle is on or off transport, respectively. Let Note that H ′ n = inf{t > 0 : X δ t ∈ δR n }, since particles cannot skip sections, and always enter a section on transport. If we were to analogously define a sequence of times {H i n } i∈I s * where for each i ∈ I s * , H i n = inf{t > 0 : X δ t ∈ δR n } with X δ 0 = δi and ξ δ 0 ∈ {0, 1}, then the exact first hitting time of the target region will satisfy H n = min{H i n : i ∈ I s * }. Hence, clearly H n ≤ H im n ≡ H ′ n . Since the exact distribution of H i n is complicated, E[H n ] is analytically intractable, and instead we focus on finding a simple expression for E[H ′ n ]. Our point is that, in the sparse material limit (q 0 small), the rightmost particle becomes increasingly likely to be the first successful particle to detach in the target region, and H n approaches H ′ n .
The computation of E[H ′ n ] requires computing the time it takes a particle to travel a certain axonal distance, given by the following.
Lemma 4.2. Let the initial position of the particle be (X δ 0 , ξ δ 0 ) = (0, 1) and let L * ∈ {1, . . . , L} be given where L is the total length of the axon. Then, the time Proof. Since (X δ t , ξ δ t ) t≥0 is a Markov process with generator it follows that M 1 t := X δ t − rδ t 0 ξ δ s ds and M 2 t := ξ δ t + t 0 k 2 ξ δ s ds − t 0 k 1 (1 − ξ δ s )ds are both martingales. Since both M 1 and M 2 have bounded increments and E[T L * ] < ∞, the optional stopping theorem implies that and our claim follows. We also note that the same computation holds for the mean time a particle in the stochastic non-compartmental model (X t , ξ t ) t≥0 takes to reach a distance L * , as the two martingales used in the proof depend only on v and not on δ.
We next compute the hitting time H ′ n of a particle started at location i m ∈ I * .
where p n = 1 − ( r k2+r ) n and λ n = k2+k1 k1 q 0 i * p n . Proof. The proof of Proposition 4.1 shows that conditioning on the values of P i (0), i ∈ R n does not affect the law of {(Q i (t), P i (t)), i ≤ i * }, hence for any t ≥ 0, they form two mutually independent sequences of Pois(q 0 ) and Pois( k2q0 k1 ) random variables.
Let S denote the number of successful particles between sites 0 and i * at time zero. The total of number of particles at time 0 at these sites that are either on or off transport is distributed as a Pois( k2+k1 k1 q 0 i * ) variable. We next compute the probability p n that any particle once it reached the target region is "successful" in detaching there. This probability can be written p n := n i=1 p(i) where p(i) is the probability it first detaches at location i * +i. Since p(i) = ( r k2+r ) i−1 k2 k2+r is the chance a particle gets off transport in the i-th compartment, p n = n i=1 p(i) = 1 − ( r k2+r ) n . The probability that a particle is successful does not depend on which location it was at time 0, and whether it was on or off transport at that time. Hence, S is is distributed as a Pois( k2+k1 k1 q 0 i * p n ) variable, and conditioned on the value S, the set of locations of the particles at time 0 is distributed as a set of S draws from the uniform distribution on {1, . . . , i * }. Since i m is the maximum of S samples from a Uniform{1, . . . , We now decompose H ′ n = H ′ e,n + H ′ o,n , where we let H ′ e,n is the time it takes a successful particle initially at location i m to enter the region R n , and H o,n is the time it takes any successful particle after it enters the region to get off transport. Thus, If a successful particle starts at location X δ 0 = i m ∈ {0, . . . , i * } on transport ξ δ 0 = 1, then the time it takes to enter the region R n is by Lemma 4.2 If it starts at location X δ 0 = i m ∈ {1, . . . , i * } off transport ξ δ 0 = 1 then it takes an additional Exponential time with mean 1/k 1 for it to get back on transport at the same location, so E H ′ e,n |(X δ 0 , ξ δ 0 ) = (i m , 0) = E H ′ e,n |(X δ 0 , ξ δ 0 ) = (i m , 1) + 1 k1 , and since we assume new particles always enter the system on transport To calculate E H ′ o,n note that if a particle first gets off transport in the i-th compartment then the time of its travel until this point is a sum of i independent and identically distributed exponential random variables with parameter k 2 + r. The probability a successful particle first gets off transport in the i-th compartment is p(i) pn . Hence, the time a successful particle takes to get off transport once it enters the region R has the mean and our claim follows.
To complete our analysis, we wish to characterize this result in terms of length along the axon and independent of the stepping size δ. To this end, we fix a length ℓ and for a given compartment size ℓ, we let n = ⌈ℓ/δ⌉, and for the start of the region we let i * = ⌈ℓ * /δ⌉. As such, as δ → 0 the limiting region becomes R ℓ = (ℓ * , ℓ * + ℓ). Then, under the assumption that q 0 /δ → ρ 0 we have It remains to interpret this result with respect to the parameter choices we have made. First, we must set a value for the size ℓ of the target region. For this purpose we note that the typical size of a Node of Ranvier -a gap in the myelin sheath of a myelinated axon where sodium channels are concentrated in the cell membrane -is approximately one micron. 1 As seen in the preceding proof there are three contributions to E[H ′ ], and we next analyze them in terms of their dominance for the overall value. We begin with the contribution from the last term, which is the time it takes for a successful particle to detach once it has reached the target region. By our choice of ℓ, the recurring ratio ℓ/v is one. Along with the earlier assumption that k 2 = 1, the entire term simplifies to (e − 2)/(e − 1) = 0.4 s. The multiplicative factor 1 − e −λ ℓ preceding the first term of (4.1) results from a boundary effect: when ℓ * is very close to zero, it is very unlikely there are any particles already in the system between the soma and the target region. When the target region is in the middle of the axon, this contribution from particles in that section of the axon at time 0 is significant.
When q 0 ∼ 1 as assumed earlier, then ρ 0 = q 0 /δ ∼ 10 8 . Then if ℓ * > 10 −8 , which corresponds to the target region being just one motor step down the length of the axon, we have ρ 0 ℓ * ∼ 1 and λ ℓ > 2.3 so 1 − e −λ ℓ > 0.9. Looking at the first term inside the parenthesis, we note that vρ 0 ∼ 10 2 , while p ℓ = 1 − exp(−ℓk 2 /v) ≈ 0.6 is the probability that any given vesicle will be successful in detaching from transport in the target region, so 1/(vρ 0 p ℓ ) ∼ 10 −2 . Meanwhile, the term k 2 /(k 1 (k 2 + k 1 )), which is the expected time to bind to transport if a successful particle happens to be off transport at time 0, is 1/4 for the assumed values of k 2 and k 1 . Therefore under the q 0 ∼ 1 assumption, both this and the final term contribute significantly to the hitting time.
However, in the sparse material regime, say q 0 ∼ 10 −3 , we have then ρ 0 ∼ 10 5 , and the factor 1 − e −λ ℓ > 0.9 when ℓ * > 10 −5 . That is, if the target region is at least 1/100, or at least 1000 segments, down the length of a 10 −3 m axon. However, now the rate limiting factor is the wait time for the first successful vesicle to arrive in the target region, which is captured by the term 1/vρ 0 p. The product vρ 0 ≈ 0.1 measures the average rate at which new particles should arrive and together with the probability of success p ℓ ≈ 0.6, we now have 1/vρ 0 p ℓ ≈ 16s. So, in the sparse material regime the first summand in E[H ′ ] giving the mean time for arrival of the particles to the target region dominates.
It is interesting to consider what happens to the arrival rate term under perturbations of the various parameters. In particular, we draw the reader's attention to changes in v. When viewing intermittent search in terms of a single particle, the probability of finding the target is strictly decreasing in v. Higher velocity seems to be the enemy of finding the target. Indeed, from a system point of view, this implies that the system will require more trials before a successful particle arrives in the target region. What the equation (4.1) gives us is the ability assess how much more quickly the trials will happen. In fact, the function v(1 − exp(−ℓk 2 /v)) is increasing in v. Therefore the entire expected wait time E[H ′ ] is actually decreasing in v. That is to say, while the particles are less likely to succeed, they will be arriving more rapidly enough to counterbalance the lost time. We believe that this kind of quantitative analysis will prove fruitful in future study when coupled with more details of the biology of deposition of materials in the cell membrane.

Approaching Equilibrium.
We have seen above that the axon is very homogeneous at stochastic equilibrium on a space scale down to micrometers. One of the beautiful properties of transport with reversible binding is that if it is locally out of equilibrium, the on-off dynamics can return the system to equilibrium on a much faster time scale than waiting for new material to arrive from the nucleus. Furthermore, as is discussed in [45], one of the key goals of biophysics investigation is the discovery of behaviors for which biochemical regulation is not necessary. This is of fundamental importance to the biological function of the system because it means that the axon will automatically "repair" itself without central control of the repair process.
How good is this mechanism? If a segment of the axon is far away from equilibrium, how long does it take to get back close to equilibrium? To investigate this question, we imagine that the axon is at stochastic equilibrium except for some segment R, where the total number of particles on and off transport, Q i (0)+ P i (0) = 0, ∀i ∈ R, is zero at time 0. Let a and ε be given small numbers and suppose that λ ∞ is the mean vector for (Q i , P i ) at stochasttic equilibrium. We want to compute (an upper bound for) the time t * so that that is, the probability that (Q i , P i ), ∀i ∈ R is significantly different from λ ∞ is very small. In the applications below we will choose a = 0.1 and ε = 0.05, and we will see that a 10 micron segment can recover in about 10 seconds, while a 1 millimeter segment will take about 1000 seconds or 15 minutes to recover. We study this question first for a single location R = {i}, and then use the estimates derived to scale the results to segments of any length.
It remains to study the convergence of the mean and the appropriate choice of t * . The dynamics of the mean vector λ(t) are given by the ODE where e 1 is the unit vector (1, 0) and A 1 = k 2 + r −k 1 −k 2 k 1 .
The solution to (4.2) is λ(t) = λ ∞ + e −A1t (λ(0) − λ ∞ ) where λ ∞ = q 0 rA −1 1 e 1 = q 0 (1, k2 k1 ). This yields the estimate |λ(t) − λ ∞ | ≤ e −A1t λ ∞ ≤ e −αt |λ ∞ | where α is the smaller of the eigenvalues of A 1 . Noting that α > 0, t * may be chosen so that e −αt * ≤ aβ, i.e. t * = α −1 ln(1/(aβ)) is an upper bound for the time to be close to equilibrium with high probability. In order to calculate return to equilibrium at various scales, we now suppose that the whole axon is in statistical equilibrium except for a segment R of length ∆ = δ10 ν in which we will assume that there are no particles either on or off the track. Proposition 4.4 covered the case ν = 0. We are interested in ν = 1, . . . , 8. We imagine that the axon is broken up into 10 8−ν segments of length ∆. In this rescaled system, the unbinding and binding rates per particle, k 2 and k 1 remain the same, as well as the mean on-transport velocity v. In order to retain this mean velocity, the rate of lateral stepping must be decreased tor = r10 −ν .
The ODE for the mean vector of the rescaled system is given by d dtλ (t) = −Ã 1λ (t) + q 0 re 1 (4.3) We note that the last term in (4.3) contains an r rather than anr. This is because the input rate is unchanged while the exit rate is diminished. The resulting equilibrium value is therefore rescaled as well,λ ∞ = q 0 rÃ −1 1 e 1 = q 0 r r 1 k2/k1 = q 0 10 ν 1 k2/k1 . Both components of this vector are of order 10 ν , as expected.
For the given parameter values (with |λ ∞ | > 5 × 10 3 in particular), the constant of proportionality ln ε|λ ∞ |/(a ε|λ ∞ | − 1) is contained in the interval (2.3, 2.5) and does not have a impact on how the relaxation time scales with ν. In terms of analyzingα, we note that k 1r is small compared to k 2 . To leading order,α ∼ (k 1r )/(k 2 + k 1 +r) ∼ 10 2−ν s −1 . It follows thatt * ∼ 10 ν−2 s. Thus, for a 10 micron segment (ν = 3) the time to recover is about 10 seconds and for a 1 millimeter segment (ν = 5) the time to recover is 1000 seconds or 15 minutes. The time to recover depends, of course, on the parameter ε that represents what we mean by "close." We also note that one can compute various measures of time to recover using the PDE models discussed in Section 3.
5. Discussion. In this paper, we created a spatial Markov chain model for studying various aspects of fast axonal transport. Previous models that use PDEs treat the velocity of transport as constant when particles are attached to the fast transport system. Since it is known that transport along the microtubules is itself stochastic, it is important to have a fully stochastic model. Our model allows us to unify and extend previous work. In Section 3.2 we show that from the particle perspective as the compartment size δ tends to zero, our model converges in distribution on the Skorokhod space of càdlàg functions to the piecewise-deterministic model analyzed by Brooks [5]. Namely, we show that the paths of particles in our model converge to those of particles in a stochastic non-compartmental model. The argument proceeds by an explicit computation of the finite-dimensional distributions and a tightness argument. In Proposition 3.4 we give a rigorous probabilistic proof of why the paths of particles follow "approximate traveling waves" described by other authors [36,5,14]. This proof is based on stochastic averaging arguments which show that a functional central limit theorem holds on the space of continuous functions for the paths of particles as the compartment size decreases. The diffusion of particles around their mean position can consequently be approximately described jointly for all time by a Brownian motion with the appropriate diffusion coefficient.
In Section 2, we show how to use existing experimental data to indentify (ranges for) all the parameters of our model. In light of this, we can use the model to investigate several important biological questions. These are based on describing the spatial distribution of multiple particles in our model. In Section 4.1 we derive the stationary distribution for the number of particles in different compartments on and off transport along the axon. This gives an explicit description of the stochasticity of the system that is present even after a long time. In Section 4.2 we derive estimates for how homogeneous the axon is on different spatial scales. In Section 4.3 we study a question introduced by Bressloff, by providing a stochastic quantity which describes the balance the system needs to achieve between rapid transport that brings new material quickly and efficient local search that improves time of delivery to a target. Finally, in Section 4.4, we use the model to calculate the length of time that it would take for axonal segments of different lengths to recover to near stochastic equilibrium after they have been depleted of vesicles.
In our stochastic compartmental model all event wait times are assumed to be exponential random variables, but this is certainly a simplification. As an example, the stepping process of kinesin is a well-studied though still not completely understood phenomenon. Much work has focused on assessing the dependence of the mean rate of translocation on both the load and the local concentration of ATP [44] [38]. Implicit in this analysis is the assumption of exponential wait times with state dependent rate parameters. However, when fitting to data and matching dispersion information the authors in [13] found it necessary to generalize the wait time distribution. This was followed by more detailed models for which it was shown that load carrying could in fact regularize the stepping times of kinesin motors [37] [9]. Generalizing waiting times would significantly affect our results. Since the particle position process is no longer Markov, we no longer have the direct connection to the previous results stated in Section 3.1.-3.3., nor can we use the stationary distribution employed in Section 4.1 and used for addressing the biologicals questions in Sections 4.2.-4.4. In light of the known need for generalized wait times in the stepping process, it seems likely that detailed observation of the rebinding process will call for new mathematical models as well. Recall that when a vesicle unbinds from a microtubule it is unclear whether it typically rebinds to the same microtubule or if it explores the region significantly via diffusion before finding a different microtubule to bind to. In the latter case, a more appropriate model for rebinding time would be to solve some kind of first passage time problem and use that distribution for the rebinding wait.
An important aspect of the biology of axonal transport is not included in the model presented here, namely the local deposition and eventual degradation of transported materials. For example, sodium channels and sodium pumps are synthesized at the soma, transported down the axon and deposited in the axonal membrane, ether uniformly as in an unmyelinated axon or at the nodes of Ranvier in a myelinated axon. Channels and pumps are proteins with half-lives on the order of days to weeks. The present model can be extended to include a deposition compartment at each location, and, clearly, the processes of deposition and subsequent degradation will cause the mean number of particles both on and off transport to be monotone decreasing as one moves down the axon. How inhomogeneous this makes the axon will depend on the details of deposition and degradation rates. Our preliminary calculations indicate that long axons, such as the meter-long axons in human sciatic nerve, would be quite inhomogeneous. This is an important biological issue because it is controversial whether the machinery for protein synthesis (i.e. ribosomes) exist in axons [42]. We have also not included retrograde transport or the fact that some axons may have location-dependent unbinding rates [10]. All of these issues will be the subject of future work.