Dynamics and length distributions of microtubules with a multistep catastrophe mechanism

Regarding the experimental observation that microtubule (MT) catastrophe can be described as a multistep process, we extend the Dogterom–Leibler model for dynamic instability in order to discuss the effect that such a multistep catastrophe mechanism has on the distribution of MT lengths in the two regimes of bounded and unbounded growth. We show that in the former case, the steady state length distribution is non-exponential and has a lighter tail if multiple steps are required to undergo a catastrophe. If rescue events are possible, we detect a maximum in the distribution, i.e. the MT has a most probable length greater than zero. In the regime of unbounded growth, the length distribution converges to a Gaussian distribution whose variance decreases with the number of catastrophe steps. We extend our work by applying the multistep catastrophe model to MTs that grow against an opposing force and to MTs that are confined between two rigid walls. We determine critical forces below which the MT is in the bounded regime, and show that the multistep characteristics of the length distribution are largely lost if the growth of an MT in the unbounded regime is restricted by a rigid wall. All results are verified by stochastic simulations.


Introduction
Assembly and disassembly of microtubules (MTs) is fundamental for many processes in the cell. The dynamics of MT polymerization is commonly described by the term dynamic instability [1]: a MT stochastically switches between a growing and a shrinking state in which it polymerizes or depolymerizes, respectively. The transition from growth to shrinkage is called catastrophe, the reverse process rescue. A first mathematical description of dynamic instability was provided by the Dogterom-Leibler model [2,3], which is based on four constant parameters: a (de)polymerization velocity for the growing and the shrinking state, and two transition rates for the occurrence of catastrophes or rescues. Depending on these parameters, the MT exhibits two dynamical regimes: a regime of bounded growth with zero mean growth velocity and a stationary exponentially decaying length distribution, and a regime of unbounded growth with a non-zero mean growth velocity.
Later on, Odde et al. [4] and, more recently, Stepanova et al. [5] and Gardner et al. [6] found that the durations of the growth intervals are not distributed exponentially as one would expect for a constant catastrophe rate. Instead, the measured distributions could be well described by a multistep process, which means that the MT ages and the catastrophe rate increases during growth. While it was concordantly reported from control groups of in vitro experiments that a MT has to pass approximately three steps to undergo a catastrophe [4,6,7], it was also shown in the same experiments that the number of steps depends on concentrations of kinesins [6] or MT-targeting agents [7].
The underlying mechanism of MT aging is still under debate and several microscopic models have been proposed. For instance, it was suggested that a catastrophe is triggered by a certain number of "sub-catastrophes" of single protofilaments [8,9,10]. A chemomechanical approach led to the conclusion that the MT tip becomes more tapered during growth, which promotes catastrophe [11]. Another model, which included Brownian dynamics of single tubulin molecules, revealed that MT aging might be a much more complex stochastic process relying on a fluctuating MT tip and an increasing number of curled protofilaments [12]. Chemomechanical stochastic growth models on the dimer level revealed that a catastrophe could be triggered by a "nucleus" of three neighboring protofilaments shrinking by more than 6 dimers, such that its GTP-cap is removed and its ends reach into the GDP-body of the MT [13,14].
In this paper, we do not concentrate on the microscopic details of MT aging but on the consequences that a multistep catastrophe mechanism has for the distribution of MT lengths. For that purpose, we extend the empirical Dogterom-Leibler model by subdividing the growing state into an arbitrary number n of sub-states a MT has to pass to undergo a catastrophe. In the bounded growth regime, the stationary form of the resulting master equations has to be solved numerically, except for the case that MTs can not be rescued. However, taking advantage of the results of Jemseena and Gopalakrishnan [15], who set up and analyzed master equations for dynamic instability with an age-dependent catastrophe rate, we are able to compute exact values for the mean MT length and higher moments and to provide an approximation that comprises the key characteristics of the length distribution for arbitrary numbers n of sub-states. While Jemseena and Gopalakrishnan made up heuristic functions to directly fit the age-dependency of catastrophe, our work is based on the model that catastrophe is a multistep process with equal transition rates for each step. Our main results are similar to those obtained by Jemseena and Gopalakrishnan: the duration of MT growth becomes less stochastic if more sub-steps are necessary to induce a catastrophe, which results in a more narrow length distribution with a lighter tail. In particular, the stationary distribution has a maximum if rescues are possible, i.e., the MT has a most probable length greater than zero, in contrast to the monotonically decreasing exponential distribution that follows from a single-step catastrophe. Going beyond the work of [15], we also examine the regime of unbounded growth, where the MT lengths approach a Gaussian distribution as in the case of a single-step catastrophe [2,3] but with a variance that decreases with the number of sub-steps that are necessary to trigger a catastrophe.
Furthermore, we apply our results for multistep MT dynamics to models of MTs under an opposing force and in a rigid confinement [16,17]. An opposing force suppresses MT growth and can counteract the effects of an increased number of required catastrophe steps n such that the critical force, below which the MT is in the bounded regime, becomes a monotonic function of n. A rigid confinement to maximum MT lengths L truncates the determined probability distributions to the accessible interval [0, L] accompanied by a finite probability to find the MT stalled at the confining boundary. Moreover, the confinement forces MTs from the unbounded regime to a stationary length distribution, which, however, does not exhibit multistep characteristics but is exponentially increasing.

Classical single-step model
The Dogterom-Leibler model [3] is an empirical coarse-grained MT model that describes dynamic instability as stochastic switches between phases of growth and shrinkage with constant velocities. Both catastrophe and rescue are described as single-step or Poisson process as sketched in figure 1(a). The model requires four parameters: growth and shrinkage velocities v ± as well as catastrophe and rescue rates ω c and ω r . The probability densities p ± (x, t) of the length x of a growing (+) or a shrinking (−) MT obey the following Fokker-Planck equations (FPEs): The mean velocity of the MT tip (time-averaged over many catastrophe and rescue cycles) is given by and can serve as an order parameter for the transition between bounded and unbounded regimes of MT growth. With a reflecting boundary at x = 0, i.e., with the MT undergoing a forced rescue as soon as it shrinks back to zero length, a negative parameter V < 0 leads to a zero mean velocity and bounded MT growth. Then, the overall probability density p(x, t) = p + (x, t) + p − (x, t) converges to a stationary exponential distribution with a finite mean length For V > 0, the parameter V is always identical to the mean growth velocity and MT growth is unbounded without a stationary length distribution. The probability density approaches a time-dependent Gaussian distribution for the length diffusivity.

Multistep model
In the following, we extend the Dogterom-Leibler model in a way that takes account of the experimental observation that MT catastrophe is a multistep process [4,5,6], see figure 1(b). For this purpose, we introduce the number of steps n a MT has to pass to undergo a catastrophe. As long as a MT has not passed all n steps, it continues growing with v + . Therefore, the growing state can be divided into n sub-states i = 1...n. In the experiments in [4,5,6], the conclusion that MT catastrophe is a multistep process was derived from the observation that growth durations τ + are gamma-distributed, with the gamma function Γ(n) = ∞ 0 t n−1 e −t dt [18]. A gamma distribution implies that each catastrophe step occurs with the same step rate ω and that backward steps are not allowed, i.e., the states i = 1...n are passed in a prescribed order as sketched in figure 1(b). Fits to the experimental data gave n ∼ 3 [4,6,7] but values of n depend on concentrations of kinesins [6] or MT-targeting agents [7]. Since rescue is still described as a single-step process, MT dynamics is now characterized by a set of five parameters n, v − , v + , ω r and ω.
Introducing n sub-states gives rise to a time-or age-dependent catastrophe rate via the general relation which holds for an arbitrary probability density p τ + (τ + ) of growth durations τ + , and where P τ + (τ ) = ∞ τ p τ + (t) dt is the survival probability of a growing MT: the growing shrinking Figure 1. Models of dynamic instability. The MT is either in a growing or in a shrinking state, in which its tip moves with constant velocities v + or v − . (a) In the classical Dogterom-Leibler model, catastrophe is a single step process with a constant rate ω c . (b) When catastrophe is modeled as a multistep process, the growing state is divided into n sub-states, and a growing MT has to pass n sub-steps, each with rate ω, before it undergoes a catastrophe. This n-step process can be summarized by means of an age-dependent catastrophe rate ω c (τ ).
catastrophe rate is only observed within the ensemble of surviving MTs. For a gamma distribution (6), we have P τ + (τ ) = Q(n, ωτ ) = Γ(n, ωτ ) Γ(n) (8) with the upper incomplete gamma function Γ(n, x) = ∞ x t n−1 e −t dt and its corresponding regularized form Q(n, x) [18]. This leads to a catastrophe rate ω c (τ ), which is monotonically increasing with time τ . For small times τ ω −1 , the final catastrophe is unlikely, because n − 1 prior sub-steps in a prescribed order are necessary that occur each with probability ωτ resulting in ω c (τ )/ω ≈ (ωτ ) n−1 /(n − 1)!. For large times τ ω −1 , n − 1 prior sub-steps have passed almost certainly with probability 1 − (n − 1)/ωτ , and the rate for the final catastrophe approaches ω, i.e., ω c (τ )/ω ≈ 1 − (n − 1)/ωτ . We note that this asymptotic form of the catastrophe rate with an algebraic approach to unity is different from the exponential approach assumed by Jemseena and Gopalakrishnan [15].
With the growth duration τ + , also the length gain x + = v + τ + during one growth interval is gamma distributed: On average, a MT grows for a duration of τ + = nω −1 , and its tip covers a distance of x + = v + nω −1 during that interval. Two generic situations will be of interest: (i) A comparison of models that yield the same (experimentally observed) average growth duration τ + but feature a different number of sub-states n. Then, we have to simultaneously use a n-dependent catastrophe step rate ω ∝ n, i.e., if more sub-states have to be passed to catastrophe we have to increase the step rate accordingly to get the same overall average growth duration or distance.
(ii) Models for experiments with microtubule associated regulating proteins, such as the catastrophe-promoting MCAK [6], where the catastrophe step rate ω remains fixed while the number of sub-steps n can be manipulated.
Together with the mean shrinking duration ω −1 r and distance v − ω −1 r , we deduce the mean tip velocity analogously to (2): Again, the sign of V n determines whether MT growth is bounded and a stationary state exists. MT growth is stabilized and may leave the bounded regime in scenario (ii) if only the number of sub-steps n to trigger a catastrophe is increased (n > n c = ωv − /ω r v + ) while the other parameters remain constant. If we compare models with the same average growth duration τ + but different number of sub-states n in scenario (i), on the other hand, we also have to use a catastrophe rate ω ∝ n, and V n becomes independent of n in (10). Then, also the sign of V n and, thus, the MT growth regime (bounded or unbounded) is unchanged by n. For a general mathematical description, we assign a probability density p i (x, t) to each sub-state i = 1...n of a growing MT. The total growing state density is given by p + (x, t) = i p i (x, t). The stochastic time evolution of the probability densities is described by a system of n + 1 FPEs: Due to the reflecting boundary, the probability current density has to vanish at x = 0. Furthermore, in any stationary state (∂ t p i (x, t) = 0), the current density is constant in space, as can be seen by summing up the FPEs (11): Together with j(x = 0) = 0, this implies that, in a steady state, the probability current density has to vanish everywhere. With the resulting relation we can eliminate p − (x) in the stationary FPEs and achieve and the abbreviations r = ω r /v − and c = ω/v + . Due to the reflecting boundary condition, a growing MT with length 0 must be in state 1, which provides the initial condition for the FPEs (15): Therewith, the solution of the FPE (15) can formally be expressed as

Force dependent growth velocities and catastrophe step rates
Here, we introduce force dependencies of the growth velocity and the catastrophe rate as similarly done in [17]. The effective velocity v + of a growing MT results from the attachment of tubulin dimers to the MT tip with rate ω on and their detachment with rate ω off : where d = 8 nm/13 ≈ 0.6 nm is the effective dimer size. An opposing force F modifies the attachment rate by a Boltzmann factor [19], yielding the force-dependent growth velocity v where we introduced the dimensionless force f = F/F 0 with F 0 = k B T /d ≈ 7 pN and the thermal energy k B T = 4.1 pN nm. As a consequence, MT growth is stalled (v + (f ) = 0) at the dimensionless stall force In the classical picture of dynamic instability, a growing MT is stabilized by a GTP-cap that is formed by GTP-tubulin dimers continuously added to the MT tip [1]. The GTP-tubulin dimers tend to be hydrolyzed after incorporation into the MT lattice so that the GTP cap may vanish and the MT undergoes a catastrophe. Since this is less likely if the addition of new GTP-tubulin dimers is fast, the catastrophe rate is expected to be negatively correlated with the growth velocity. Different microscopic models of dynamic instability were proposed that may be employed to deduce the velocitydependence of the catastrophe rate [20,21,8,9]. Here, we use the phenomenological approach from [17] that is based on the experimental observation that the mean growth duration τ + is a linear function of growth velocity, with a = 13.8 s 2 nm −1 and b = 20 s [22]. For different multistep catastrophe mechanisms that yield the same experimentally observed average growth duration τ + = nω −1 but feature a different number of substates n (scenario (i)), this implies which conforms with (22) independently of the choice of n. Therefore, the step rate ω(v + ) inherits a force dependence from the growth velocity via (20). Because the experimental observations suggest sub-steps that occur with the same rate ω(v + ), also their force-dependence should be identical. Moreover, we assume that the number of catastrophe steps n is a structural feature of the MT that depends neither on the growth velocity nor on the applied force. Then, we find the force-dependent step rate with v + (f ) from (20).

Confinement between rigid walls
Above, we assumed unrestricted MT growth. In the following, we introduce a rigid wall that confines the MT length to a box between the reflecting boundary at x = 0 and the position of the wall x = L. Similar models for confined MTs with a single-step catastrophe were discussed in [16,17]. If a MT tip reaches the wall, its growth is stalled (v + = 0). A MT that reaches the wall in growing state i has to pass all n − i subsequent states to undergo a catastrophe and to leave the wall. Due to the velocity dependence of the catastrophe step rate (23), reaching the wall induces catastrophes with an increased step rate ω L ≡ ω(v + = 0).
We introduce probabilities Q i to find the MT stalled at the wall while it is in growing state i. Since the MT leaves the wall instantaneously after a catastrophe, we can assume that Q − = 0. The time evolution of Q i (t) is given by While the first terms describe the catastrophe steps of a MT that is already stalled, the expressions v + p i (L, t) provide the probability currents onto the wall of a MT in state i. The total probability to find a MT tip at the wall is given by yields the total probability current onto the wall. In the stationary state (∂ t Q i (t) = 0), (25) is solved by which adds up to the total probability In the interior of the confining box (x < L), MT dynamics is still described by the FPEs (11). Therefore, the stationary probability density of the length of a confined MT is given by p conf (x) = p + (x) + p − (x) + Qδ(L − x) and has to satisfy the normalization condition Since Q is a linear combination of p i (L), both Q and p ± (x) are proportional to the constant of integration p 0 introduced in the initial condition (17). Therefore, normalization can be achieved by adjusting p 0 , just like in the unconfined case.

Bounded growth
We start our investigation of the multistep MT model with the force-free case. In general, the solution p + (x) = exp(M x) p 0 of the FPEs (15) can only be evaluated numerically, e.g., by numerical diagonalization of the coefficient matrix M as described in Appendix A. However, for the experimentally relevant case [6] that a MT can not be rescued (ω r = r = 0) except for the forced rescue at the boundary x = 0, the solution can be expressed analytically: where Q(n, x) = Γ(n, x)/Γ(n) is the regularized form of the upper incomplete gamma function (see (8)).
In order to approach a solution of the general case (ω r > 0), we make use of the results of Jemseena and Gopalakrishnan [15], who calculated the Laplace transform of the steady state length distribution for the case of an age-dependent catastrophe rate ω c (τ ), where the age τ is the time that has passed since the last rescue event.
Given an arbitrary probability density p τ + (τ + ) of growth durations τ + , the associated age-dependent catastrophe rate is given by (7). Combining that with the results of Jemseena and Gopalakrishnan [15], we find with the Laplace transformp x + (s) of the probability density of growth distances p x + (x). The derivation is presented in Appendix B. To achieve a result for the n-step catastrophe process, we substitute the Laplace transform of p x + (x) from (9) into (32): If r = 0, the inverse Laplace transform yields the probability density from (30). For the general case with rescues (r > 0), there is no analytical result for the inverse Laplace transformp(s) → p(x). Nonetheless, we are able to compute exact results for the mean MT length x and the variance Var(x) = x 2 − x 2 by interpreting the Laplace transform as moment-generating function: The method can be easily extended to higher moments and cumulants, which are listed in Appendix C up to the fourth degree. Moments x m of the length distribution are dominated by large x and, thus, by the analytical properties of the Laplace transformp(s) for small s around s = 0. On the other hand, we can approximatep(s) for large s as Then, the inverse Laplace transform is possible and provides an approximation of p(x) for short MT lengths: where P (n, x) = 1 − Q(n, x) is the regularized lower incomplete gamma function [18].
In the following, we compare our analytical results with stochastic simulations that solve the equation of motion of the MT for fixed time steps ∆t and include the random occurrence of a catastrophe step or a rescue after each time step. Here, we consider situation (i) and want to compare models with the same average growth duration τ + and, thus, the same parameter V n in (10) but different number of sub-states n. Then, we have to use a catastrophe step rate proportional to n, with a constant ω 0 , and the MT remains in the bounded regime (V n < 0) independently of n. We use the following values: which are typical for bounded MT growth with negative V n and correspond to tubulin concentrations around c tub ∼ 10 µM ( [17,14] and references therein). Figure 2 shows the results in absence of rescue events (r = 0) except for the forced rescues at x = 0. The analytical predictions from (29) and (30) perfectly match with the results from the simulations. The overall probability densities are monotonically decreasing functions and converge towards a step function c 0 Θ(c −1 0 − x) for large n. This uniform distribution for an infinite step catastrophe process can be made plausible by considering the growth distances x + : Since the standard deviation of the gamma distribution (9) is given by ∆x + = √ nc −1 , the relative error of growth distances ∆x + / x + = 1/ √ n vanishes for large n. Moreover, as we assumed that c = nc 0 , also the absolute deviation decreases as 1/ √ n whereas the mean growth distance stays constant. Consequently, the more steps a MT has to pass to undergo a catastrophe, the more deterministic and predictable the length gain becomes. In the infinite step limit, the MT tip always covers the same distance during one growth interval. Then, in the absence of rescue events, a MT grows from x = 0 to exactly x = c −1 0 where it undergoes a catastrophe, and shrinks back to zero length where it is rescued again, finally resulting in a uniform distribution of MT lengths. Dynamically, repeated growth and shrinking by a sharp distance x = c −1 0 results in deterministic oscillations with a sharp period Similar effects have been discussed in [15]. Increasing the number of sub-states n at fixed mean growth duration thus sharpens the distribution of growth lengths and gives rise to quasi-deterministic MT length oscillations.
If rescues are possible (r > 0), the probability density functions are not monotonic anymore but increase exponentially for short MT lengths up to a maximum, see   (29) and (30) (red lines) match with the distributions measured in simulations (bars). (c) For an infinitestep process, the overall probability density converges to a step function.
figure 3(a). The exponential increase and the maximum are well described by the approximation (38). After the maximum, however, the approximation deviates from the real distribution. In that region, the probability densities measured in stochastic simulations are only fitted well by the numerical solution according to Appendix A, which is also the case for the single state densities p i (x) depicted in figure 3(b). If the number of steps n increases, the maximum becomes sharper and moves towards longer MT lengths up to x = c −1 0 . As we show in Appendix D, in the infinite step limit, the probability density approaches a piecewise defined function that initially grows exponentially as (c − r) exp(rx) until it has a step discontinuity at x = c −1 0 . Moreover, there are non-analyticities of higher order at each multiple of c −1 0 . Like in the absence of rescues, this behavior can be explained with the determinism of MT growth distances in the infinite step limit: Since now the rescue rate is greater than zero, a MT can be rescued before shrinking to zero length and is able to grow beneath the single growth distance c −1 0 . Nevertheless, a MT that grows from zero length after a forced rescue and undergoes a catastrophe at x = c −1 0 still shrinks back to zero with the probability exp(−r/c 0 ) = 74 %. These 74 % alone would result in a step function again, and only in 74 % of the growth cycles that start from x = 0, the MT reaches lengths x > c −1 0 , finally leading to the step discontinuity.
As it can be seen in figure 4, (35) and (36) correctly describe the mean length and its variance as measured in stochastic simulations. If the number of catastrophe steps and the step rate are increased proportionally, the mean length decreases by up to one half of the single-step value.
The length distributions in Figs. 2(c) and 3(c) exhibit a characteristic change of shape as a function of the sub-state number n. Measurements of complete length distributions will, thus, allow to draw conclusions on this parameter related to the catastrophe mechanism. Also the reduction of the mean length x in Fig. 4 to one half of the single-step value for increasing n will allow to determine the sub-state number n   (35) and (36).
if the average growth duration τ + , rescue rate and growth and shrinkage velocities v ± are known.

Unbounded growth
In the regime of unbounded growth, there is no stationary solution which is why we have to analyze the time dependent FPEs (11). As described in detail in Appendix E, by use of the Fourier transform and the implicit function theorem, we are able to calculate a dispersion relation that is valid for long MTs. As V n > 0 in the unbounded regime, the limit of long MTs is reached after long times. We conclude from the approximated dispersion relation that the probability density of MT lengths approaches a Gaussian distribution as in (4), but with the mean velocity V n given by (10) and the diffusion constant D n = n(n + 1) 2 One can easily see that D n , and thus the standard deviation of MT lengths vanish for large n if the other parameters remain constant. Together with (10), this means that an increase of the number of catastrophe steps favors growth over shrinkage and thereby can make the MT leave the bounded regime, where a further increase of n favors growth even more until catastrophes are almost completely suppressed. In figure 5, the results are compared to stochastic simulations. Here again, we consider situation (i) that the catastrophe step rate is proportional to n so that the mean growth duration and the mean velocity are constant. Then, the diffusion constant D n does not vanish for large n but still decreases by up to one half of the single-step value. We use the same parameters as in (40) but with a ten times higher rescue rate in order to induce an unbounded state with V n > 0. At the beginning of the simulations both the initial and, since the MTs are still short, the boundary condition cause the length distribution to significantly deviate from the Gaussian approximation. After long times, however, the MTs gained so much length that the probability of reaching the boundary is negligible. Then, the approximation provides an accurate description of the simulation results.

Microtubule growth against a force
An opposing force f affects the growth velocity and the step rate as described in section 2.3. We start our investigation with a MT that grows against a constant force and distinguish the two scenarios introduced above for the step rate: (i) We consider an experimentally given mean growth duration τ + as in (22), and consider models with different numbers n of sub-steps that conform with this experimental data. Then we use ω(f ) = nω 0 (f ) as in (24). (ii) Experimentally, also the number of sub-steps n might be changed without affecting the catastrophe step rate ω(f ), e.g. by altering the MCAK concentration [6]. Then, we have a given ω(f ), for which we use ω(f ) = ω 0 (f ) with ω 0 (f ) from (24) to obtain agreement with the experimental data of [22] for n = 1.
For both scenarios, the opposing force influences the MT length distribution via a force dependent sub-step catastrophe rate ω(f ) and via a force dependent growth velocity v + (f ) as explained in section 2.3. Therefore, application of an opposing force does not give rise to novel shapes of MT length distributions but simply shifts parameters, in particular, the mean velocity parameter V n = V n (f ). This gives rise to other important effects, namely the existence of a critical force f c for the boundary between bounded and unbounded regime [17].
While we vary n, f and ω r in the following analysis, the other parameters are fixed throughout this section to v − = 200 nm s −1 , ω off = 6 s −1 [23], ω on = 39.33 s −1 , where ω on was chosen such that the growth velocity and the step rate used in the previous sections (see (40)) are reproduced for f = 0.
Because of the force dependence of growth velocity and step rate, also the mean velocity V n in (10) becomes force-dependent. In general V n decreases with increasing force f such that a critical force f c exists above which MT growth transitions from the unbounded to the bounded regime because it is suppressed by force. This critical force is determined from the condition V n = 0 and is always smaller than the stall force f stall = ln(ω on /ω off ) = 1.88, which is set by the condition v + = 0 so that the MT is not able to leave the boundary at x = 0 and, therefore, must be in the bounded regime.
For the first scenario (i) with a given τ + and ω(f ) = nω 0 (f ) as in (24), the velocity V n is n-independent. Likewise, the critical force f c is n-independent with Dynamics and length distributions of microtubules with a multistep catastrophe mechanism16 Therefore, measurements of the critical force will not allow to deduce information about the number of sub-steps n.
For the second scenario (ii), where the number of sub-steps n is changed without affecting the catastrophe step rate ω(f ) = ω 0 (f ) with ω 0 (f ) from (24), on the other hand, the velocity V n depends on n (see (10)). Therefore, we also find an n-dependent critical force above which the MT is in the bounded regime. For both scenarios, the critical force is shown in figure 6(a) with different rescue rates ω r . In the second scenario, which applies to experiments employing MCAK [6], a value of n could be determined from experiments measuring the critical force f c . If the parameter values lie in the unbounded regime in the absence of force, the critical force f c can be determined by increasing the force f such that MT growth becomes limited to a finite mean length x < ∞ for f = f c at the transition to the bounded regime.
Substituting v + and ω in (35) with their force-dependent versions, we determine force dependent mean lengths x , which are shown as function of f and n in figures 6(b) and (c), respectively. For the first scenario, we find a very weak n-dependence as in the absence of force (see (35) and figure 4). In the second scenario, where the number of sub-steps n is changed without affecting the catastrophe step rate, we find a much more pronounced n-dependence as also evidenced in figures 6(b) and (c). Again, we arrive at the conclusion that in the second scenario, a value for the sub-step number n could be determined from experiments measuring the mean length x (f ).
The asymptotic behavior of f c and x is summarized in table 1. In the second scenario, the critical force becomes a monotonic function of n, that grows logarithmically for small n and converges to f stall for n → ∞, see figure 6(a) and table 1. Moreover, the mean length grows as a function of n and diverges at n c = ω(f )v − /ω r v + (f ) indicating a switch to the unbounded regime, see figure 6(c). Analogously, we observe diverging mean lengths in figure 6(b) for both scenarios when the applied force approaches the critical force f c .

Microtubule dynamics under a rigid confinement
So far, we assumed that the MT can grow freely to arbitrary lengths. However, the more realistic situation both in vitro and in vivo is that the MT is confined to a finite length. In the following, we investigate how multistep MT dynamics is affected by a rigid wall that restricts the MT length to x ≤ L as described in section 2.4. In contrast to the reflecting boundary at x = 0, the MT tip stays at the wall after reaching the length L until it has completed the remaining catastrophe steps with an increased step rate ω L = ω(v + = 0). In the stationary state, the probability Q to find the MT stalled at the wall (x = L) is given by (27) while the probability density in the interior (0 < x < L) is still defined by the FPEs (15). Therefore, the probability density p conf (x) of a confined MT follows from truncating the unconfined probability density at x = L and weighting it with 1−Q in order to fulfill the normalization condition (28), see figure 7(a). As for an unconfined MT, the FPEs have to be solved numerically in the general case with rescues (ω r , r > 0), see Appendix A. However, in presence of a confinement to 0 < x < L, an analytical solution for moments is no longer possible via the Laplace transform because the Laplace transform is defined over the whole half-space 0 < x.
Again, we distinguish the two scenarios (i) of a step rate that is proportional to n (ω = nω 0 ) and (ii) of a constant step rate (ω = ω 0 ), applying the force-free parameters from (40). Then, according to (23), the step rate at the wall is ω L = nω 0 (v + = 0) = n/b or ω L = ω 0 (v + = 0) = 1/b, respectively. In the proportional Table 1. Asymptotic behavior of the critical force f c and the mean MT length x in various limits, cf. figure 6.

Limit
Asymptotic behavior where the MT stays in the bounded regime for any n, the length distribution only retains the characteristics from figure 3 with a maximum and a sharp decrease thereafter if the position L of the wall lies beyond the mean growth length x + = c −1 0 = 5.9 µm, see figure 7(a). If L falls below the mean growth and the mean shrinkage length x − = r −1 = 20 µm, the MT length will approach a roughly uniform distribution as the MT is likely to both grow and shrink from boundary to boundary without catastrophes or rescues in the interior.
For a constant catastrophe step rate (ii), the MT switches from the bounded to the unbounded regime at n = n c = 3.4. The growth of a MT in the unbounded regime is now confined by the rigid wall so that the MT lengths approach a stationary distribution as in the bounded regime. Since V n > 0 in the formerly unbounded regime, the MT tends to the wall at x = L and greater MT lengths are more likely, see figure 7(a). It is important to note that the multistep characteristics of the length distribution, in particular the maximum, are lost if V n > 0. Instead we find exponentially increasing probability densities because the distance from the right boundary that can be reached during a shrinkage-growth-cycle that starts from the wall is dominated by the single-step rescue. In order to restore a distribution with a maximum, rescue had to be a multistep process.
The behavior of the length distributions also becomes manifest in the stall probabilities Q and the mean lengths x depicted in figures 7(b-e). In the proportional step rate scenario (i), Q decreases with n ( figure 7(b)) since the lighter tailed distributions make it less likely for the MT to reach the wall. With a constant step rate (ii), on the other hand, Q approaches 1 for two reasons: firstly, leaving the bounded For ω = nω 0 (top and center), the decisive factor for the shape of the distribution is whether the right boundary lies beyond the maximum, which coincides with the mean growth length x + = c −1 0 = 5.9 µm for n → ∞. For ω = ω 0 (bottom), the MT switches from the bounded to the unbounded regime between n = 3 and n = 4 (n c = 3.4). Accordingly, the MT lengths are shifted from being cumulated near the left boundary towards the right one. The total weight in the interior of the box decreases from n = 10 to n = 100 and further as the MT tends to be stalled at the wall for a longer time, see (b). (b,c) Probability Q to find the MT at the wall and mean MT length x as function of n for L = 1, 10, 100 µm. For ω = nω 0 , Q decreases with n due to the lighter tails of the length distributions. Since the MT is always in the bounded regime, the wall at x = L becomes irrelevant for L → ∞ and the mean length converges to (35), which is depicted in figure 4. For ω = ω 0 , a larger n shifts the MT length towards the right boundary and increases the residence time at the wall. Consequently, Q approaches 1 and x approaches L for n → ∞. (d,e) Probability Q to find the MT at the wall and mean MT length x as function of L for n = 2, 3, 4, 10, 100. If the MT is in the bounded regime (ω = nω 0 or ω = ω 0 and n = 2, 3), the wall can be neglected for large L so that Q vanishes and x approaches the mean length of a free MT from (35). In the unbounded regime (ω = ω 0 and n = 4, 10, 100) the MT tends to the wall and is unlikely to shrink to zero length if L is large. Therefore, Q converges towards a finite value and x grows linearly with L. regime by increasing n drives the MT towards the wall; secondly, the residence time at the wall is longer if more steps are required to leave it. As a consequence of Q → 1, the mean MT length approaches L if n is increased with a constant step rate, see figure 7(c).

Figures 7(d) and (e)
show Q and the mean MT length, respectively, as a function of L. Q(L) decreases monotonically starting from Q(0) = 1, where the wall coincides with the reflecting boundary at x = 0. The behavior for larger L depends on whether an unconfined MT with the same parameters would be in the bounded (V n < 0) or in the unbounded regime (V n > 0): In the bounded regime, it is unlikely that the MT reaches the wall if L is large, and, hence, Q(L) vanishes while the mean length converges to the value of an unconfined MT as given by (35). For V n > 0, on the other hand, the MT stays always in the vicinity of the wall so that its mean length becomes a linear function of L, and Q converges to a finite value because the probability density behaves as if the left boundary did not exist. We draw the general conclusion that the left (x = 0) or the right boundary (x = L) can be neglected for V n > 0 or V n < 0, respectively, if L is sufficiently large.
The results in figure 7 suggest that, in scenario (i), it might be possible to determine the number of catastrophe sub-steps n from the length distribution if the MT is confined to sufficiently large compartments L 10 µm so that the results for an unconfined MT can be applied. Then, also the probabilities Q to find the MT at the wall display a characteristic n-dependence, however, its absolute values are to small for an accurate measurement. In scenario (ii), which applies to experiments employing MCAK [6], a value of n could be determined more easily from length distributions and Q-measurement but also from experiments measuring the mean length x .

Discussion
Based on experimental results that characterize MT growth periods as gamma distributed and conclude that catastrophe is a multistep process [4,6], we extended the empirical Dogterom-Leibler model [2,3] in order to analyze the consequences a multistep catastrophe mechanism has for the distribution of MT lengths. The multistep process has two main effects on the growth durations of a MT, which also underlie the consequential changes in the length distributions: Firstly, if the number of catastrophe steps is increased while keeping the rate of a single step constant, the growth durations become longer and the MT may leave the bounded regime. Secondly, the growth periods and hence the length gain during one growth interval are less stochastic if more steps are necessary to trigger a catastrophe.
In the case of bounded growth, the stationary length distribution has a steep descent in the vicinity of the mean growth distance n/c. In absence of rescues, the steep descent follows an area where the probability density is only slowly decreasing and becomes nearly constant for large n. If rescues are allowed, the distribution is exponentially increasing and has a maximum before it decreases sharply. In both cases, the length distributions are lighter tailed than the exponential distribution resulting from a singlestep catastrophe, i.e., a multistep catastrophe reduces the number of MTs that are longer than the mean growth length. As a consequence, the mean MT length decreases by up to one half the single-step value if the number of catastrophe steps and the step rate are increased proportionally (see figure 4).
If the average growth duration τ + is known, we show that measurements of the shape of the bounded length distributions (Figs. 2(c) and 3(c)) or the mean MT length will give access to the number of sub-states n that best describe the experimental data (scenario (i)).
We also discussed the situation where MT dynamics is altered by MT regulators that do not only affect the velocities or the transition rates but the number of catastrophe sub-steps [6] (scenario (ii)). We conclude that by such a regulation, a MT acquires more potent ways to adapt to special situations inside the cell: while altering the classical four parameters only adjusts the range of MT lengths, which stay exponentially distributed in the single-step case (as long as they do not leave the bounded regime), variation of the additional parameter n changes the shape of the length distribution. As similarly discussed by Gardner et al. [6], this could be beneficial during mitosis. For instance, during prometaphase, the steep descent in the length distribution can appropriately limit the area that is explored by MTs in order to fasten search-andcapture of chromosomes [24]. In metaphase, accumulation of MT lengths around the maximum of the distribution may support the precise positioning of chromosomes in the metaphase plate and the maintenance of spindle length.
Stationary length distributions that have a maximum for short MT lengths before they apparently decrease exponentially and are similar to the ones in figure 3 were measured in several experimental studies [25,26,27,28,29,30]. Though some of these studies have already been cited as evidence for a multistep catastrophe mechanism [6,31], there are different reasons why this interpretation is dubious. In contrast to our model with a fixed minus and a dynamic plus end, the in vitro studies in [25,29] examined free MTs that could polymerize simultaneously at both ends. Therefore, the shape of the length distribution can be rationalized by a convolution of the respective exponential distributions at the plus and the minus end, which are both obeying single-step dynamics [29]. In [26], the deviation from an exponential distribution for short MT lengths is attributed to the image resolution being to low to detect very short MTs. The results in [27,28,30] seem to be more in line with our model, yet these publications do not provide a quantitative evaluation of the measured length distributions, which makes a valid conclusion difficult. Besides, [27,28] are in vivo experiments so that additional effects from microtubule associated proteins or spatial restrictions are likely.
In the regime of unbounded growth, the MT lengths approach a Gaussian distribution as in the single-step case but with a reduced variance. In vivo, the stabilization of MT growth due to an increase of the number of catastrophe steps might help interphase MTs, which have been shown to be in the unbounded regime [2,32], to reach the cell boundary. On the other hand, at the transition from interphase to mitosis, MT lengths are significantly reduced in order to prepare the mitotic spindle assembly [33,28]. The restructuring of the MT array may be supported by a reduction of the number of catastrophe steps, which destabilizes the MTs and shifts them to the bounded regime. This hypothesis is supported by the observation of Gardner et al. [6] that MCAK, which plays a key role for the control of MT dynamics during mitosis [34], promotes catastrophes by reducing the required steps from n = 3 to n = 1 and simultaneously keeping the step rate ω constant.
We also added a force-dependence of growth velocity and step rate to the model and analyzed how force affects multistep MT dynamics in the scenario of a step rate that is proportional to the number of required steps (ω(f ) = nω 0 (f ), scenario (i)) as well as for a step rate that does not depend on n (ω(f ) = ω 0 (f ), scenario (ii)). If the opposing force exceeds the critical force f c , the mean velocity V n changes its sign and the MT switches from the unbounded to the bounded regime.
With a step rate that is proportional to n, the mean velocity and, thus, the critical force are constant. This scenario is useful for determining the number of sub-states n if the mean growth duration τ + = nω 0 is known. Because the critical force is strictly n-independent, measurements of the critical force will not allow to deduce information about the number of sub-steps n.
When the step rate is independent of n, the force f and an increase of n counteract each other because the force suppresses growth and enhances catastrophe steps, thereby decreasing the mean growth duration and length, while a larger n results in an increased mean growth duration τ + = nω 0 . Thus, the critical force is a monotonic function of n and converges to the stall force, at which the growth velocity vanishes. Manipulating the step number without affecting the step rate is possible by a variation of the MCAK concentration [6]. In this situation measurements of the critical force f c give information on the value of n.
If the MT is confined to finite lengths L, for instance by the cell cortex, the MT will either tend to shrink back to zero length if it is in the bounded regime, or it will grow repeatedly against the confining boundary if MT growth would be unbounded without the confinement. In the latter case, the stationary length distribution has no maximum but increases exponentially up to the boundary. Since this is due to the single-step kinetics of the rescue, a maximum could be restored by a multistep rescue, which, however, has not been observed so far to our knowledge. As aforementioned, MTs are in the unbounded regime during interphase [2,32], where the microtubule organizing center is positioned by direct pushing and/or by dynein-mediated pulling interactions between the plus ends and the actin cortex [35,36,37,38]. Since it would be disadvantageous for these interactions if the MT distribution had a maximum and the most probable MT length lay in front of the cell cortex, a single-step rescue and an exponential length distribution are indeed favorable in this situation. Moreover, we argued above that interphase MTs might be in the unbounded regime as a result of an increased number of required catastrophe steps n. This mechanism would further support the interaction with the cell cortex since the probability Q to find the MT plus end at the confining boundary correlates positively with n if the step rate is constant.

Conclusion
In conclusion, a catastrophe mechanism that requires multiple steps has significant effects on the length distribution of MTs. Modifying the number of required catastrophe steps, e.g., by regulation via MCAK [6], allows to adapt not only the scale but also the shape of the length distribution to be beneficial for the present physiological situation. The multistep characteristics of MT length are retained under an opposing force; the critical force to suppress unbounded growth becomes sensitive to modification of the number of catastrophe steps. Confined MTs can continue to show characteristic maxima in the MT length distribution only in the bounded regime and if the compartment is sufficiently long; otherwise MT, length distributions tend to simple exponentially decreasing or increasing distributions in the bounded and unbounded case, respectively.

Appendix A. Numerical solution of the stationary Fokker-Planck equations
We solve the stationary FPE (15) by numerical determination of the eigenvalues λ j and -vectors v j of the matrix M . The solution can then be written as with S = ( v 1 , ..., v n ). The coefficients a j follow from the initial condition (17) at x = 0: with a = (a 1 , ..., a n ) T . Finally, the constant p 0 is determined by the normalization condition (note that Re λ j < 0 for any eigenvalue): If a MT is confined to x ≤ L by a rigid wall as described in section 2.4, the probability density in the interior of the box is still given by (A.1) and (A.3), but is accompanied by the probability Q to find the MT tip stalled at the wall, which is a linear combination of p i (L) in the stationary state, see (27). Therefore, for a confined MT, the constant p 0 follows from the normalization condition (28), which has to be fulfilled by the total probability density p conf (x) = p(x) + Qδ(L − x): Appendix B. Solution for an age-dependent catastrophe rate Jemseena and Gopalakrishnan [15] found that in case of an age-dependent catastrophe rate ω c (τ ) and a reflecting boundary at x = 0, the Laplace transformed overall probability densityp(s) of MT lengths is given bỹ Substituting the catastrophe rate ω c (τ ) = −∂ τ ln P τ + (τ ) that follows from an arbitrary distribution of growth durations (see (7)), we find Ω(τ ) = − ln P τ + (τ ), (B.5) The last line is valid for small s. Next, we substitute x = v + τ to show that ζ(s) is the Laplace transform of the probability density of growth distances: For the purpose of normalization, we use the approximated form of ζ(s) from (B.6): 1 =p(0) = J 0 x + s s − r x + s s=0 = J 0 x + 1 − r x + , (B.8) Substituting (B.7) and (B.9) into (B.1) finally results inp(s) as given in (32).