Effects of aging in catastrophe on the steady state and dynamics of a microtubule population
EEffects of aging in catastrophe on the steady state and dynamics of a microtubulepopulation
V. Jemseena ∗ and Manoj Gopalakrishnan † Department of Physics, Indian Institute of Technology Madras, Chennai 600036, India (Dated: November 10, 2018)Several independent observations have suggested that catastrophe transition in microtubules isnot a first-order process, as is usually assumed. Recent in vitro observations by Gardner et al.[M. K. Gardner et al., Cell , 1092 (2011)] showed that microtubule catastrophe takes place viamultiple steps and the frequency increases with the age of the filament. Here, we investigate, vianumerical simulations and mathematical calculations, some of the consequences of age dependenceof catastrophe on the dynamics of microtubules as a function of the aging rate, for two differentmodels of aging: exponential growth, but saturating asymptotically and purely linear growth. Theboundary demarcating the steady state and non-steady state regimes in the dynamics is derivedanalytically in both cases. Numerical simulations, supported by analytical calculations in the lin-ear model, show that aging leads to non-exponential length distributions in steady state. Moreimportantly, oscillations ensue in microtubule length and velocity. The regularity of oscillations,as characterized by the negative dip in the autocorrelation function, is reduced by increasing thefrequency of rescue events. Our study shows that age dependence of catastrophe could function asan intrinsic mechanism to generate oscillatory dynamics in a microtubule population, distinct fromhitherto identified ones.
PACS numbers: 87.17.Aa, 87.16.Ln, 05.40.-a
I. INTRODUCTION
The dynamic instability in microtubules [1, 2], in par-ticular, the catastrophe transition by which a growing fil-ament abruptly starts shrinking, has been the subject ofextensive experimental and theoretical studies over sev-eral decades. The recently discovered phenomenon ofaging in microtubule catastrophe by Gardner et al. [3]brings to light some new aspects of the catastrophe tran-sition. Contradicting the long-standing view that (in theabsence of depolymerizing proteins) a filament is equallylikely to undergo catastrophe at any instant of time, itwas observed that the probability of a microtubule toundergo catastrophe depends on how long it has beengrowing; the older the microtubule is, the higher theprobability to undergo switching. Irrespective of tubu-lin concentration, the measured catastrophe frequencyexhibited a nonlinear dependence on the age of the mi-crotubule, increasing linearly during the early stage ofgrowth and then approaching a steady state ( Fig. 3.Eof Ref. [3]). Also, the presence of depolymerizing pro-teins was observed to affect the aging process; for ex-ample, Kip3p accelerated the rate of aging while in thepresence of MCAK, catastrophe was observed to becomea first order process. Until this observation was made, itwas generally believed that microtubule catastrophe is asingle step stochastic event with no memory associatedwith it and hence usually modeled as a first order pro-cess. Within this frame work, many theoretical modelshave been set up in an attempt to understand the intrin- ∗ [email protected] † [email protected] (corresponding author) sic dynamics of microtubules [4–7]. It has been shown byVerde et al. [4], the microtubule length distribution fol-lows a simple exponential decay in the bounded growthregime, when the switching rates are constants.Available experimental and simulation studies suggestthat the origin of the aging phenomenon lies in the multi-filament structure of a microtubule [3, 8]. A microtubuletypically consists of 13-14 protofilaments wrapped into acylindrical structure. During the early stages of growth,when all the protofilaments are of nearly equal length,the entire filament is structurally stable because of stronglateral bonding between protofilaments. As the filamentgrows longer, the tip of the microtubule becomes moreand more tapered, weakening the lateral bonding be-tween protofilaments and makes the filament more sus-ceptible to undergo catastrophe. In this view, the ”struc-tural defects” responsible for the switching of the fila-ment grows with time, until finally the filament under-goes catastrophe. An alternate view is that catastropherequires the loss of guanosine triphosphate (GTP) tip ina minimum number of protofilaments, and hence its ki-netics could be visualized as a multi-step process, eachstep requiring the completion of the previous one for itsmaterialization [3, 9, 10].The possibility that the kinetics of plus end in micro-tubule dynamics may be non-first order was suggestedfirst by Odde et al. nearly two decades ago [11], and morerecently by Stepanova et al. [12]. In Ref. [11], based ona probabilistic analysis of the experimentally observedgrowth time distributions, the authors showed that theplus and minus ends of a microtubule follow differentkinetics. For the minus end, the growth time distribu-tion could be described as exponential decay, character-istic of first order kinetics, whereas for the plus end, itturned out to be non-exponential; the experimental data a r X i v : . [ q - b i o . S C ] J u l in this case fitted well to gamma distribution. Since thegamma distribution characterizes a multi step process,it was inferred that a growing microtubule (at the plusend) has to go through a sequence of events in order thatcatastrophe occurs. A new parameter was introduced todescribe dynamic instability, apart from ν c (catastrophefrequency), ν r (rescue frequency) , v g (growth velocity)and v s (shrinking velocity), viz., the shape parameter r ofthe growth time (or equivalently, time until catastrophe,hence also called catastrophe time) distribution (gammadistribution). The value of the shape parameter that gavereasonably good fit with the experimental data is r ∼ θ ∼ . − . Later, Odde and Buettner[13] showed that oscillations in state (growth/shrinkage)arise as a natural consequence of non-first order kineticsof catastrophe and rescue processes.In the present paper, we develop a mathematical modelof dynamic instability, wherein the evolution of catastro-phe frequency with the age of the filament is explicitlytaken into account. Motivated by experimental obser-vations of Gardner et al.[3], we assume here that age-dependent catastrophe frequency is given by the phe-nomenological expression ν c ( τ ) = ν max c (cid:0) − e − λτ (cid:1) , (1)where ν max c is the asymptotic value, λ is referred to as theaging rate, with τ being the age , the time spent by thefilament in the growing state after the last rescue event.In this sense, we assume the age of a filament as beinggiven by an internal clock in the filament, the readingof which is reset to zero at every catastrophe, and whichstarts “ticking” as soon as rescue happens. For simplicity,and for lack of strong evidence to the contrary, we assumethat rescue events follow first-order kinetics characterizedby rate ν r . A fit of Eq.1 to the experimental data in [3]is shown in Fig. 1.The motivation for choosing the above expression isprimarily its mathematical simplicity, which makes itsuitable for further calculations. A similar, but not ex-actly identical, mathematical form for aging of catastro-phe also emerges naturally within a recently proposedmodel[10] where the transition is modeled as a first pas-sage process, such that a threshold number of protofil-aments are required to lose their individual GTP tipsby hydrolysis for it to happen. For details, we refer thereader to Supplementary Information [14]. In Ref. [3],by comparison, it is assumed that the catastrophe timedistribution (to be discussed in Sec. II) is represented bygamma distribution, similar to [11]. Here too, the catas-trophe transition is visualized as a first passage event,and requires the occurrence of a threshold number r ofunderlying events, but each event is assumed to occurwith the same rate θ . In this case, the precise nature ofthe underlying events is not specified in detail. Both the aging time τ (sec) ν c ( s ec - ) Experimental data βτ + ∆ Eq.1
FIG. 1. A fit of the experimental data given in figure 1. Din Ref. [3] , using the expression for age-dependent catastro-phe frequency given by Eq.1. A simpler linear fit given byEq.7 is also shown, see Sec. II for discussion. The best fitparameter values for the nonlinear curve given by Eq.1 are ν max c = 0 . − and λ = 0 . − , while for the linearcurve β = 0 . − and ∆ = 0 . − . exponential aging function in Eq.1 as well as the corre-sponding function derived from the gamma distributionfit the experimental data almost equally well; see Fig. 1in Supplemental Material [14].A brief summary of our results are as follows: as adirect consequence of Eq.1, the distribution of time in-tervals corresponding to the growth phase is found tobe non-monotonic (Section II). The Fokker-Planck equa-tions for microtubule dynamics in the presence of aging(set up in Sec. III) lead to a number of analytical re-sults, presented in Sec. IV. In particular, the stationarylength distribution is non-exponential, approaching half-Gaussian form in the limit λ →
0, while the standardexponential decay is recovered in the limit λ → ∞ . Ex-act mathematical expressions for length autocorrelationfunction are derived in Laplace space, but in general,explicit inversion is found difficult. Numerical simula-tions (Sec. IV) show that the autocorrelation functionpossesses a negative lobe, signifying oscillatory behav-ior, the half-period of which is measured as a functionof aging rate λ, ν max c and v g /v s , the ratio of growth andshrinkage velocities. Numerically computed velocity au-tocorrelation function also shows similar properties. Fornear-experimental parameter values, the half-period oflength oscillations is found to be in the range of ∼ − II. CATASTROPHE TIME DISTRIBUTION
In this paper, we introduce the age of the filament τ asa variable that describes the state of a growing filament(in addition to its length) irrespective of the origins ofthe aging effect. Although rescue events are not reportedin the experiments by Gardner et al. [3], here we takea step a forward: we assume that if rescue events arepresent, after every rescue, a filament starts to grow as afresh filament without any structural defects. Age is re-set to zero at each catastrophe event, which again startsto grow linearly with time upon the rescue of the fila-ment. Since our objective is to explore the consequencesof age-dependent catastrophe on dynamic instability, wedo not consider explicitly microscopic events such as ad-dition/removal of monomers and hydrolysis; these enterthe formalism implicitly through catastrophe and rescuerates[6, 10, 15, 16].In the present study, we use the two state continuumapproach model proposed by Verde et al. [4], where thedynamics is captured by a pair of coupled first order par-tial differential equations at a macroscopic level. Manymathematical models have been formulated in the past,based on this two state approach to explore dynamicinstability and its applications, eg., microtubule oscilla-tions [17, 18], dynamics under confinement[19, 20], searchfor targets[20–22] and the effects of regulators [23, 24]. p r ob a b ilit y d e n s it y ( m i n - ) λ =5min -1 λ =0.5min -1 λ =0.05min -1 β =0.05min -2 α ν c m a x T c FIG. 2. Probability distribution function for catastrophe time(growth time) given by Eq. 3, for different λ values with ν max c = 0 . − . The dotted curve gives the catastrophetime distribution for the linear model with β = 0 . − .In the inset the scaled average catastrophe time for the expo-nential model, given by Eq. 10 as a function of α is shown. Let us denote by T the time until catastrophe (referredto as the catastrophe time ), starting from a rescue eventat T = 0. Then, the probability distribution for T isgiven by f c ( T ) = ν c ( T ) exp (cid:32) − (cid:90) T ν c ( τ ) dτ (cid:33) . (2)After substituting Eq.1 in Eq.2, we find f c ( T ) = ν max c e α (cid:20) (1 − e − λT ) exp[ − ( ν max c T + αe − λT )] (cid:21) , (3)for which the following limiting cases are of special inter-est: f c ( T ) ∼ ν max c e α exp[ − ν max c T ] , λT (cid:29) f c ( T ) ∼ βT e − βT , λT (cid:28) , (5)where we have defined the parameters α = ν max c λ ; β = ν max c λ. (6)Note that the latter regime ( λT (cid:28) ν c ( τ ) = βτ, (7)which we shall call the “linear model”, to distinguish itfrom the “exponential model” given by Eq.1 (note thateven the linear model fits the experimental data reason-ably well, see Fig.1). Both models will be studied in thispaper. The linear model turns out to be more amenableto mathematical analysis; see, for instance, Sec. IV B forexplicit analytical expression for length autocorrelationfunction.The distributions in Eq.3 and Eq.5 are non-monotonicfunction of T (see Fig.2), with the maxima, correspond-ing to the most probable switching time, located at T m = 1 λ ln (cid:20) α (1 + 2 α − √ α ) (cid:21) (exponential model) ,T m = 1 √ β (linear model) . (8)In the limit α (cid:29)
1, the expression for the exponentialmodel approaches that of the linear model, as expected.The average time until catastrophe is calculated as T c = (cid:90) dT f c ( T ) T. (9)Using Eq.3 and Eq.9, we find, for the exponentialmodel, T c = α e α ν max c ∞ (cid:88) n =0 ( − n α n n ! (cid:20) α + n ) − α + ( n + 1)) (cid:21) , (10)with the limiting behavior T c ∼ /ν max c as λ → ∞ , corre-sponding to the constant catastrophe case. For the linearmodel in Eq.7, we find, similarly, T c = 14 (cid:114) πβ . (11)The non-monotonic behavior of f c ( T ), by virtue of the“tightening” of the direction reversal times (see Fig. 2),opens up the possibility of oscillatory dynamics, which weshall explore in detail in Sec. IV B later using autocor-relation functions. Numerical simulations indicate thatoscillations occur when the aging rate (in the exponentialmodel) and rescue frequency are sufficiently small. III. CONTINUUM EQUATIONS FOR THETWO-STATE STOCHASTIC MODEL
According to the two-state stochastic model, a givenmicrotubule exists either in the growing state or shrink-ing state during its life time. A third state called pausestate is also possible where the filament neither growsor shrinks [2]. In the present study, we ignore the exis-tence of the pause state, as well as the structural detailsof the filament, although such details may implicitly af-fect the catastrophe rate. A microtubule nucleates froma nucleation center at a rate ν (‘birth’) and it may dis-appear completely by shrinking to length zero (‘death’).In between birth and death, a microtubule switches be-tween growing and shrinking states; switch from growingstate to shrinking state (catastrophe) takes place at rate ν c ( τ ) (given by Eq.1, including the limiting linear form)and the reverse transition takes place at rate ν r (rescue)which is taken to be a constant. The growth velocityis denoted by v g , the shrinkage velocity by v s , and thelength by x .Let us now denote by G (cid:48) j ( x, τ, t | x , dxdτ the proba- bility to find a filament in the growing state with lengthand age in the range [ x, x + dx ] and [ τ, τ + dτ ] respectively,at time t . Similarly, G j ( x, t | x ) dx gives the probabilityto find a filament in the shrinking state with length inthe range [ x, x + dx ] at time t . In both cases, x is theinitial length and j = 1 ,
0, indicates the state of growthor shrinkage at t = 0. Note that the age variable τ is rel-evant only in the growing state, and is therefore presentonly in one set of Green’s functions, for which the primednotation is used.The variables x and τ evolve according to the followingequations: dxdt = v g ; dτdt = 1 (growing state) ,dxdt = − v s ; τ = 0 (shrinking state) . (12)The Green’s functions then satisfy the following dif-ferential equations, which suitably generalize the corre-sponding equations for constant catastrophe rate, withaging absent[4, 5, 7]: ∂G (cid:48) j ( x, τ, t | x , ∂t = − v g ∂G (cid:48) j ( x, τ, t | x , ∂x − ∂G (cid:48) j ( x, τ, t | x , ∂τ − ν c ( τ ) G (cid:48) j ( x, τ, t | x ,
0) + ν r G j ( x, t | x ) δ ( τ ) , (13) ∂G j ( x, t | x ) ∂t = v s ∂G j ( x, t | x ) ∂x + ∞ (cid:90) dτ ν c ( τ ) G (cid:48) j ( x, τ, t | x , − ν r G j ( x, t | x ) , x ≥ . (14)Note that resetting of age to zero in the shrinking stateis achieved through the introduction of delta functionin Eq.13.The Green’s functions satisfy the normalizationcondition: ∞ (cid:82) [ G j ( x, t | x ) + G j ( x, t | x )] dx = 1, where G j ( x, t | x ) ≡ ∞ (cid:90) G (cid:48) j ( x, t, τ | x , dτ. (15) IV. RESULTSA. Length and age distribution
Many microtubule-dependent functions in the cell de-pend on the spatial organization of microtubules, andhence, a quantity of primary interest is their length distri-bution. In the case of constant switching rates, Verde etal. [4] have shown that a microtubule population reachesa steady state with a stationary length distribution, when ν c v s > ν r v g , (16) in which case the lengths are exponentially distributed,with mean length (cid:104) x (cid:105) = ( ν c /v g − ν r /v s ) − . To under-stand how these results are modified in the presence ofage-dependence in catastrophes, we start by defining thestationary length distributions P (cid:48) ( x, τ ) ≡ lim t →∞ G (cid:48) j ( x, τ, t | x , ,P ( x ) ≡ lim t →∞ G j ( x, t | x ) , (17)which satisfy the following steady state equations, ob-tained by putting the LHS to zero in Eq.13 and Eq.14. v g ∂P (cid:48) ( x, τ ) ∂x = − ∂P (cid:48) ( x, τ ) ∂τ − ν c ( τ ) P (cid:48) ( x, τ ) , + ν r P ( x ) δ ( τ ) (18) v s ∂P ( x ) ∂x = − ∞ (cid:90) ν c ( τ ) P (cid:48) ( x, τ ) dτ + ν r P ( x ) . (19)To solve Eq.18 and Eq.19 we assume, for simplicity,that nucleation rate is very large such that nucleation ofa new filament takes place instantaneously once a micro-tubule disappears at that site; hence the net flux at theboundary is zero. The resulting flux balance conditionat the boundary x = 0 can be written as v s P ( x ) δ ( τ ) | x =0 = v g P (cid:48) ( x, τ ) | x =0 (20)and we define, J ≡ v s P ( x ) | x =0 . We now Laplace-transform x → p in Eq.18 and Eq.19 by defining˜ P (cid:48) ( p, τ ) = ∞ (cid:90) e − px P (cid:48) ( x, τ ) dx ; ˜ P ( p ) = ∞ (cid:90) e − px P ( x ) dx. (21)After applying the above transformations in Eq.18 andEq.19, followed by simplifications (see Appendix A fordetails), we arrive at the expression˜ P (cid:48) ( p, τ ) = (cid:20) J v s pv s p − ν r (1 − ζ ( p )) (cid:21) × exp (cid:18) − v g pτ − τ (cid:90) ν c ( τ (cid:48) ) dτ (cid:48) (cid:19) , (22)where the function ζ ( p ) is given by Eq.A8. After inte-grating out τ from Eq.22, we obtain˜ P ( p ) = (cid:20) J v s pv s p − ν r (1 − ζ ( p )) (cid:21) η ( p ) , (23) η ( p ) is given by Eq.A9. Similarly,˜ P ( p ) = J (1 − ζ ( p )) v s p − ν r (1 − ζ ( p )) . (24)Both Eq.23 and Eq.24, together with the normaliza-tion property of the probability distribution functions(PDFs) are used to find the expression for J for specificforms of ν c ( τ ). Exponential aging:
We use the form of ν c ( τ ) given inEq.1 to calculate ζ ( p ) and η ( p ) from the integrals givenby Eq.A8 and Eq.A9. The expression for J , fixed usingnormalization, i.e., (cid:80) j =0 , ˜ P j (0) = 1, turns out to be J = ν max c v s − ν r v g [1 + e α α − α γ ( α + 1 , α )]( v g + v s )[1 + e α α − α γ ( α + 1 , α )] . (25)Here γ ( α + 1 , α ) is the lower incomplete gamma func-tion, defined [25] γ ( a, z ) = z (cid:90) e − t t a − dt, (26) and the dimensionless constant α is defined in Eq.6. Inthe limit α → λ → ∞ ), the term in the square bracketapproaches unity. Steady state is guaranteed only if thenumerator of the expression in Eq.25 is positive, thusfor age-dependent catastrophe of the form in Eq.1, thecondition given by Verde et al. [4] is modified as ν max c v s > ν r v g [1 + e α α − α γ ( α + 1 , α )] . (27)In the limit λ → ∞ the above condition reduces to ν max c v s > ν r v g .Unfortunately, it is not possible to invert the Laplace-transforms in Eq.23 and Eq.24 explicitly, except for thespecial (but experimentally relevant) case ν r = 0. Aftersubstitution of the required ζ ( p ) and η ( p ) from Eq.B1and Eq.B2 respectively, we find that the length distribu-tion P ( x ) = (cid:80) j =0 , P j ( x ) has the form (for details seeAppendix B) P ( x ) = A exp (cid:18) − ν max c v g x − αe − λvg x (cid:19) , (28)with the normalization constant given by A = ν max c e α v g (cid:18) e α α − α γ ( α + 1 , α ) (cid:19) . (29)The state-specific distributions are given by P ( x ) = φP ( x ) and P ( x ) = (1 − φ ) P ( x ) where φ = v g / ( v g + v s ). The various limiting behaviours of the distributionin Eq.28 are as follows: P ( x ) ∝ exp (cid:18) − λν max c v g x (cid:19) , λ (cid:28) ν max c P ( x ) ∝ exp (cid:20) − (cid:18) ν max c v g x + e − ν max cvg x (cid:19)(cid:21) , λ = ν max c P ( x ) ∝ exp (cid:18) − ν max c v g x (cid:19) , λ (cid:29) ν max c . (30)The first two moments of the distribution in Eq.28 areof interest: these are (cid:104) x (cid:105) = A v g λ ∞ (cid:88) n =0 ( − n α n n ! (cid:20) α + n ) (cid:21) (31)and (cid:104) x (cid:105) = 2 A v g λ ∞ (cid:88) n =0 ( − n α n n ! (cid:20) α + n ) (cid:21) . (32) Linear aging : It is also interesting to study separatelythe implications of linear increase in catastrophe withage. In this case, the current at the boundary turns outto be J = v s − ν r v g (cid:113) π β ( v s + v g ) (cid:113) π β . (33)The inverse Laplace transform of Eq.23 and Eq.24,with η ( p ) and ζ ( p ) given by expressions Eq.B7 and Eq.B8respectively in Appendix B, leads to the length distribu-tion (in this case, half-Gaussian) P ( x ) = (cid:115) βπv g e − β v g x , (34)which is identical to the first limiting case in Eq.30.Some parallels with the related problem of length-dependent catastrophe are worth mentioning here. Inthe absence of rescue, length of a filament is given by v g τ , i.e, x ∝ τ , suggesting that in this case, the age-dependent and length-dependent catastrophes are effec-tively the same, although the underlying mechanisms arequite different. Experimental evidence suggests that theintrinsic dynamics of microtubules can be regulated bythe presence of microtubule associated proteins. A spe-cial class of proteins known as depolymerizing proteinsare capable of enhancing the catastrophe rate. A widelystudied example is the microtubule depolymerizer Kip3p,belonging to the kinesin 8 family, which is known to in-crease the catastrophe rate in a length-dependent man-ner [26]. Several authors have studied deploymerizer-induced length-dependent depolymerization and catas-trophe [24, 27–30]. Tischer et al. [24] have shown that inthe presence of length-dependent transition rates (bothrescue and catastrophe) the length distribution developsa peak, implying a tighter regulation of lengths. In fissionyeast experiments reported in Tischer et al.[31] whererescue was found to be absent but catastrophe increasesnearly linearly with length of the filament, it was pre-dicted that length distribution is half-Gaussian, identicalto Eq.34. Age distribution : The steady state age distribution inthe growing state is defined as ρ ( τ ) = ∞ (cid:90) P (cid:48) ( x, τ ) dx, (35)where P (cid:48) ( x, τ ) is given by inverse Laplace transform ofEq.22 with respect to p . The result (which is independentof ν r ), is ρ ( τ ) = λα α exp( − ν max c τ − αe − λτ ) γ ( α, α ) (exponential model)(36) with limiting behaviors ρ ( τ ) ∝ e − ν max c τ , λτ (cid:29) ρ ( τ ) ∝ e − βτ , λτ (cid:28) . (37)In the linear model, the age distribution is half-Gaussian at all times, consistent with Eq.34: ρ ( τ ) = (cid:114) βπ e − βτ (linear model) . (38) B. Autocorrelation functions
We will now investigate further on the effects of aging,and the consequent change in the catastrophe time dis-tribution, on correlation functions. In general, we expectthat the observed peak in the distribution of catastrophetime (Fig.2) would make the catastrophe transition moredeterministic and predictable, and hence may confer anoscillatory nature to microtubule dynamics. Microtubuleoscillations is a well-studied phenomenon experimentally[32–36], but here, the oscillations generally refer to peri-odic changes in free tubulin concentration due to poly-merization and depolymerization of filaments. The ex-ception is the work by Odde and Buettner[13]; here,the authors considered dynamic instability as a two-statejump process with gamma-distributed residence times ineach state. The autocorrelation function for the growthstate was shown to exhibit oscillations.In this study, we first looked at the autocorrelation inlength of a microtubule assembly. In terms of the Green’sfunctions G ij (Eq.13 and Eq.14) and probability distri-butions P j defined in Sec. III and Sec. IV.A respectively,the partial autocorrelation functions for length is definedas (cid:104) x (0) x ( t ) (cid:105) ij = (cid:90) dx (cid:90) dxx x G ij ( x, t | x ) P j ( x ) . (39)Next, by summing Eq.39 over all the initial and finalstates, we define C (cid:48) x ( t ) = (cid:88) i,j (cid:104) x (0) x ( t ) (cid:105) ij , (40)in terms of which the normalized autocorrelation is givenby C x ( t ) = C (cid:48) x ( t ) − (cid:104) x (cid:105) σ x , (41)where σ x = (cid:104) x (cid:105) − (cid:104) x (cid:105) is the variance in length. Thenormalization ensures that C x (0) = 1.The calculations of the autocorrelation function for theexponential model turned out to be cumbersome, andhence explicit results were derived only in the limit λ →∞ (constant catastrophe). For the linear model, explicitresults could be obtained in Laplace space. We presentour results for these cases now.
1. Case 1: Constant catastrophe, λ → ∞ When we put λ = ∞ strictly, the Green’s functions G j depends only on x , and not on τ . Then, Eq.13 andEq.14 reduce to ∂G j ( x, t | x ) ∂t = − v g ∂G j ( x, t | x ) ∂x − ν max c G j ( x, t | x )+ ν r G j ( x, t | x ) , (42) ∂G j ( x, t | x ) ∂t = v s ∂G j ( x, t | x ) ∂x + ν max c G j ( x, t | x ) − ν r G j ( x, t | x ) . (43)Eq.42 and Eq.43 are solved by taking Laplace transformboth with respect to time and space (see equations Eq.C1and Eq.C2 in Appendix C). Using Eq.41, and by putting ν r = 0 in Eq.C1 and Eq.C2, after a series of calculationswe arrive at the following simple expression for lengthautocorrelation function: C x ( t ) = ( v s e − ν max c t − v g e − vsvg ν max c t ) v s − v g ( ν r = 0) . (44)In particular, in the limit v s → v g the above expressionbecomes, C x ( t ) = e − ν c t (1 + ν max c t ) ( ν r = 0 , v g = v s ) . (45)
2. Case 2. Linear aging, ν c ( τ ) = βτ Eq.13 and Eq.14 are solved by taking transforms, with ν c ( τ ) given by Eq.7 (see Appendix C). By putting ν r =0 in Eq.C9 and Eq.C10, and using Eq.40, we find theLaplace transform of the unnormalized correlation as˜ C (cid:48) x ( s ) = v g sβ + v s e a s erfc( as ) s e b s erfc( bs ) (cid:20) v g (cid:114) πβ − v g sβv s e c s erfc( cs ) (cid:21) − v s v g s (cid:114) πβ + v g β (cid:114) π β e a s erfc( as ) e c s erfc( cs ) , (46)where the constants a, b, c are defined as a = 1 √ β ; c = v g v s √ β ; b = a + c. As consistency checks for the expression in Eq.46, wenote thatlim s →∞ s ˜ C (cid:48) x ( s ) = lim t → (cid:104) x (0) x ( t ) (cid:105) ≡ (cid:104) x (cid:105) = v g λν max c , (47)whilelim s → s ˜ C (cid:48) x ( s ) = lim t →∞ (cid:104) x (0) x ( t ) (cid:105) ≡ (cid:104) x (cid:105) = 2 v g πλν max c , (48)both of which agree with direct calculations, startingfrom the distribution in Eq.34.It turns out difficult to extract more information aboutthe time dependence of the function C (cid:48) x ( t ) from the trans-form in Eq.46, as the latter does not have a standardform which can be explicitly inverted. Therefore, to un-derstand it better, we turn to numerical simulations; theresults are presented in the next subsection.Lastly, the velocity autocorrelation ( numerical resultswill be discussed in the next subsection) is studied bydefining a state variable S ( t ), which takes value +1 or − t . Then the instan-taneous velocity of the filament is v ( t ) = v g (cid:18) S ( t )2 (cid:19) − v s (cid:18) − S ( t )2 (cid:19) (49)and the variable S ( t ) switches from +1 to − ν c ( τ ), and the reverse transition occursat rate ν r . Therefore, the velocity autocorrelation canbe easily expressed in terms of the state autocorrelationfunctions, defined as C (cid:48) S ( t ) = (cid:88) i,j (cid:104) S (0) S ( t ) (cid:105) ij , C S ( t ) = C (cid:48) S ( t ) − (cid:104) S (cid:105) σ S , (50)where σ S = 1 −(cid:104) S (cid:105) is the variance in S ( t ). In particular,when v s = v g , v ( t ) = v g S ( t ), and hence the normalizedvelocity autocorrelation C v ( t ) = C S ( t ) itself. C. Numerical simulations
Given the difficulties in inverting the Laplace trans-forms, it became necessary to carry out stochastic nu-merical simulations to complete our study. Fixed timestep (rather than Gillespie [37]) algorithm is used in sim-ulations as the distribution of catastrophe time is non-exponential when aging is present. In the simulations,each of the microtubule is allowed to evolve until a max-imum time, independent of the others. We choose a timestep dt such that 0 < Rdt <
1, where R can be ν c ( τ ), ν r and ν . The nucleation rate is chosen to be very largecompared to the other rates in the simulation in orderto achieve the boundary condition given by Eq.20. Amicrotubule in the growing (shrinking) state persists in µ m)00.010.020.030.040.05 p r ob a b ilit y λ =0.5min -1 λ =5min -1 µ m)00.0050.010.015 p r ob a b ilit y Analytical result ν r =0.1min -1 ν r =0 FIG. 3. Microtubule length distribution for the exponentialmodel. The parameter values are ν = 6min − , ν r = 0, ν max c =0 . − , v g = 1 µ m min − , v s = 1 . µ m min − . In theinset, a comparison of length distribution for ν r = 0 and ν r = 0 . − for λ = 0 . − is shown. The line is shownfor the analytical result given by Eq.28. µ m)00.050.10.15 p r ob a b ilit y β =0.5 min -2 β =5 min -2 µ m)00.0050.010.0150.02 p r ob a b ilit y Analytical result ν r =0.1min -1 ν r =0 FIG. 4. Microtubule length distribution for the linearmodel. The parameter values are ν = 6min − , ν r = 0, v g = 1 µ m min − , v s = 1 . µ m min − . In the inset, a compar-ison of length distribution for ν r = 0 and ν r = 0 . − for β = 0 . − is shown. The line is shown for the analyticalresult given by Eq.34. that state and elongates by v g dt (shrinks by v s dt ) untilit encounters catastrophe (rescue).We studied length distribution, average length as wellas length and velocity autocorrelation functions. Bothexponential and linear catastrophes are studied. Fromthe experiments reported in [3], for a tubulin concentra-tion of 12 µ M, we estimated that λ ≈ . − and ν max c ≈ . − , obtained by fitting the experimentaldata to the function in Eq.1, see Fig.1. Keeping thesevalues as a reference point in the parameter space, inthe simulations, we varied λ (keeping ν max c fixed at0 . − ) in the range 0 . − − . For simulationswith the linear aging model, we varied β in Eq.7 inthe range 0 . − − . Although rescue events werenot observed in the experiments reported in [3], weconsidered the effects of non-zero rescue separately. Thesimulation results are discussed in detail now. < x > ( µ m ) λ =0.05min -1 λ =0.5min -1 λ =5min -1 ∆ x ( µ m ) FIG. 5. Time evolution of average microtubule length forthe exponential model, with the inset showing time evolutionof standard deviation in length. The parameter values are ν = 6min − , ν r = 0, ν max c = 0 . − , v g = 1 µ m min − , v s = 1 . µ m min − . < x > ( µ m ) β =0.05min -2 β =0.5min -2 β =5min -2 ∆ x ( µ m ) FIG. 6. Time evolution of average microtubule length forthe linear model, with the inset showing time evolution ofstandard deviation in length. The parameter values are ν =6min − , ν r = 0, v g = 1 µ m min − , v s = 1 . µ m min − . (i) Length distribution and mean length
The length distribution for exponential (Fig.3) and lin-ear (Fig.4) aging models are shown, along with the cor-responding analytical results. The main plots show theresults with zero rescue; the effect of including small non-zero rescue frequency can be seen in the insets.The time evolution of the average length is shown inFig.5 (exponential) and Fig.6 (linear). As expected in-tuitively, the average length increases with the increasein the “memory” of catastrophe (smaller λ or β ). Itis observed from these plots that aging results in non-monotonic time evolution of average length for suffi-ciently small λ (or β , in the linear model), indicating thepossibility of oscillations. For the linear case, the peakin the curve remained noticeable even for large values of β (Fig.6).(ii) Length autocorrelation C x ( t ) λ =5min -1 λ =0.5min -1 λ =0.05min -1 C x ( t ) FIG. 7. Microtubule length autocorrelation function for theexponential model. The parameter values are ν = 6min − , ν r = 0, ν max c = 0 . − , v g = 1 µ m min − , v s =1 . µ m min − . In the inset, oscillatory behavior of the au-tocorrelation function corresponding to the parameter value λ = 0 . − is highlighted. C x ( t ) β =5min -2 β =0.5min -2 β =0.05min -2 C x ( t ) FIG. 8. Microtubule length autocorrelation function for thelinear model. The parameter values are ν = 6min − , ν r = 0, v g = 1 µ m min − , v s = 1 . µ m min − . In the inset, oscillatoryparts of the autocorrelation functions are shown enlarged. Simulation results show that, in the presence of ag-ing in catastrophe, the length autocorrelation has a pro-nounced (negative) minimum, followed by much weakerhigher order extrema, characteristic of damped oscilla-tions; see Fig.7 and Fig.8. For the exponential model(see Fig.7), the minimum loses depth and eventually dis-appears with increasing λ (eg., for ν max c = 0 . − , al-most no trace of oscillations are seen for λ > . − ).For the linear model, however, oscillations in the auto-correlation are seen to be present for a wide range of β values studied, spanning three orders of magnitude.When rescue is switched on, the autocorrelation func-tion loses its oscillatory character gradually in both mod-els. This is shown in Fig.9 for exponential (main plot)and linear models (inset) respectively. In the exponen-tial model, given the parameters v g = 1 µ m min − , v s =1 . µ m min − , ν max c = 0 . − and λ = 0 . − ,the autocorrelation correlation appears to follow multi- C x ( t ) C x ( t ) ν r =0 ν r =0.005min -1 ν r =0.01min -1 ν r =0.02min -1 ν r =0.05min -1 ν r =0.1min -1 FIG. 9. Effect of rescue frequency on length autocorrelationfunction. The main plot is shown for exponential model ( λ =0 . − , ν max c = 0 . − ) and inset for linear model( β =0 . − ). The parameter values are ν = 6min − , v g =1 µ m min − , v s = 1 . µ m min − . λ (min -1 )101001000 h a l f- p e r i od ( m i n ) Exponential β (min -2 ) 0.1110100 Linear g /v s )5101520 ExponentialLinear(a) (b)
FIG. 10. Half period of length autocorrelation function.( a ) corresponds to half-period as a function of λ , with in-set showing the same for the linear model as a function of β . The parameter values are ν = 6min − , v g = 1 µ m min − , v s = 1 . µ m min − and for the exponential model ν max c =0 . − . The dashed line is a fit of β − (Eq.51). ( b ) corre-sponds to half-period as a function of v g /v s . The parametervalues are ν = 6min − , ν max c = 0 . − , λ = 0 . − (exponential model) and β = 0 . − (linear model). exponential decay (without the negative lobe) when therescue frequency ν r is nearly 0 . − ; in the lin-ear model, this occurs at ν r = 0 . − , given β =0 . − .We also estimated the half-period of the oscillations inlength (about its mean value) from the autocorrelationfunction, as the position of the first minimum in C x ( t ),determined from simulation results. Fig.10 ( a ) shows thedependence of the half-period on the aging rates λ and β (inset): the logarithmic plots reveal that the half-period ∝ β − in the linear model, which is expected. Based onthe expression for autocorrelation function in Eq.46, wepropose the following scaling law for the linear model,when ν r = 0:0 L e ng t h ( µ m ) ν r =0.1min -1 ν r =0.01min -1 ν r =0.01min -1 ν r =0.1min -1 time(min) (c) (b)(d)(a) ν r =0.01min -1 ν r =0.1min -1 (e) (f) FIG. 11. Microtubule trajectories from stochastic simula-tion, with the dashed horizontal line denoting average lengthin the respective models studied. ( a ) and ( b ) are for the ex-ponential model with λ = 0 . − and ν max c = 0 . − .( c ) and ( d ) are for the linear model with β = 0 . − .( e ) and ( f ) for the λ → ∞ case with ν max c = 0 . − .All the trajectories are simulated by fixing v g = 1 µ m min − , v s = 1 . µ m min − . Only a and c correspond to oscillatorydynamics, as characterized by the negative lobe in the auto-correlation function. T o = 1 √ β f (cid:18) v g v s (cid:19) , (linear model) (51)where T o is the half-period of oscillations. It is clear that f ( x ) is a increasing function of x ; as v s → ∞ , the micro-tubule will crash to the origin as soon as it undergoescatastrophe, and the period is determined by averagetime spent in the growing state, specified by Eq.10. As v s is reduced, this time interval increases, which adds tothe oscillation period. This is confirmed in simulations.Fig.10( b ) shows the half-period plotted as a function of v g /v s in both exponential and linear models, which showsthat the half-period (for fixed λ and β ) is an increasingfunction of v g /v s .For the exponential model, scaling considerations sug-gest that T o = 1 √ β g (cid:18) λν max c , v g v s (cid:19) , (exponential model) (52)such that lim x → g ( x, y ) = f ( y ) for consistency withEq.51. Oscillations are present only when λ/ν max c is sufficiently small. A few sample trajectories of amicrotubule undergoing dynamic instability, under theage-dependent catastrophe are shown in Fig.11.(iii) Velocity autocorrelation
The normalized velocity autocorrelation function,computed from the simulations, also shows oscillatorybehavior (Fig.12) for sufficiently small values of λ (andall values of β studied in this paper). Also, increasing ν r lifts up the negative lobe, weakening the oscillations.Interestingly, a comparison with Fig.9 reveals that, for C S ( t ) time(min) -0.500.51 C S ( t ) ν r =0 ν r =0.01min -1 ν r =0.05min -1 FIG. 12. State autocorrelation function C S ( t ) against time t with a boundary at x = 0. The main plot is shown for thelinear model with β = 0 . − and inset for exponentialmodel with λ = 0 . − . The other parameter values are ν = 6min − , ν max c = 0 . − , v g = 1 µ m min − , v s =1 µ m min − . the same value of λ , the threshold rescue frequency fordisappearance of oscillations is higher for velocity auto-correlation when compared to length, implying that fora certain range of parameter values, direction reversalsare more regularly spaced in time when compared to ex-cursions of length above and below the mean.Our results for velocity/state correlation are similarto the results reported by Odde and Buettner [13]. Inthis paper, the authors modeled dynamic instability asa two state jump process, but switching from one stateto another state takes place in a sequence of steps suchthat residence times in both the states are characterizedby (identical) gamma-distributions. Qualitatively, theeffects of age-dependent (studied here) and multi-step(studied in Ref.[13]) catastrophes are similar; however,in our model, only catastrophe is assumed to be age-dependent while rescue is treated as a first order process. V. SUMMARY AND CONCLUSIONS
Motivated by recent experimental observations on “ag-ing” in microtubule catastrophe, in this paper, we haveformulated a mathematical model to investigate how thestatistical properties of a microtubule population getsmodified by it. Aging is characterized by a parameter λ , with units of inverse time, which we refer to as theaging rate in our model. Our results show that agingaffects statistical properties of a microtubule populationin important ways. The steady state length distribu-tion is no longer a simple exponential decay as in theconstant catastrophe case, while the condition for its ex-istence is altered. The length and velocity autocorre-lation functions develop a negative lobe for sufficientlysmall aging and rescue rates, which signifies oscillatorydynamics in the population. We characterized the os-cillations using analytical calculations and scaling argu-1ments, and checked their consistency with numerical sim-ulations. Velocity (or growth state) oscillations in micro-tubule dynamics were studied earlier by Odde and Buet-tner [13] using a model with identical gamma-distributedresidence times in each state (i.e, both catastrophe andrescue were assumed to undergo aging in a similar way).Our study is similar in spirit, but with the following dif-ferences: (i) two different forms of aging are used forcatastrophe, while rescue is treated as a first order pro-cess (in the absence of evidence to the contrary) (ii) westudy how the steady state length distribution is modifiedin the presence of aging (iii) the length autocorrelationfunction is studied in detail, and we have derived an ex-plicit mathematical expression for the same (in the linearaging model) in Laplace space.Microtubule length oscillations are likely to be rele-vant in the context of the well-known phenomenon ofchromosome oscillations, although the latter is certainlya more complex process involving many other proteinsand regulatory mechanisms. Nevertheless, we feel that itis interesting to make a broad comparison of our resultswith experimental observations on chromosome oscilla-tions. In eukaryotic cells, soon after the onset of mitosis,each chromosome pair exists in a “mono-oriented” statefor a while, i.e., pulled/pushed only from one side of themitotic spindle [38–40] by microtubules attached to itthrough macromolecular structures called kinetochores.A successful and efficient mitosis requires all the chro-mosomes to be “bi-oriented” (i.e., bound to microtubulebundles emanating from the opposite poles simultane-ously) at the spindle equator before segregation. It isknown from experimental observations that both mono- oriented and bi-oriented chromosomes undergo regularoscillatory motion.Skibbens et al. [38] have showed that mono-orientedchromosomes in Newt lung cells exhibit non-sinusoidaloscillatory motion with a time period of 200 − ∼ ∼ −
40s [41, 42]. With the availableinformation, it is unclear whether aging and relatedeffects in microtubule catastrophe plays a role in thesephenomena; however, given the fact that the observedperiod of oscillations of mono-oriented chromosomes[38, 40] is close to the what was found in our simulations(Fig.10), we feel that further investigations in thisdirection could be worthwhile.
ACKNOWLEDGMENTS
V.J thanks Erwin Frey for useful discussions during avisit to LMU, Munich. The authors acknowledge P.G.Senapathy Centre for Computing Resources, IIT Madrasfor computational support. The authors would also liketo thank Dr. M.K. Gardner for sharing the relevant ex-perimental data reported in [3]. [1] T. Mitchison and M. Kirschner, Nature , 232 (1984).[2] A. Desai and T. Mitchison, Annu. Rev. Cell Dev. Biol. , 83 (1997).[3] M. K. Gardner, M. Zanic, C. Gell, V. Bormuth and J.Howard, Cell , 1092 (2011).[4] F. Verde, M. Dogterom, E. Stelzer, E. Karsenti, and S.Leibler, J. Cell Biol. , 1097 (1992)[5] M. Dogterom and S. Leibler, Phys. Rev. Lett. , 1347(1993).[6] H. Flyvbjerg, T. E. Holy and S. Leibler, Phys. Rev. Lett. , 2372 (1994); Phys. Rev. E , 5538 (1996).[7] D. J. Bicout, Phys. Rev. E , 6656 (1997).[8] C. E. Coombes, A. Yamamoto, M. R. Kenzie, D. J. Odde,and M. K. Gardner, Current Biol. , 1342 (2013).[9] H. Bowne-Anderson, M. Zanic, M. Kauer and J. Howard,Bioessays , 452 (2013).[10] V. Jemseena and M. Gopalakrishnan, Phys. Rev. E ,032717 (2013).[11] D. J. Odde, L. Cassimeris and H. M. Buettner, Biophys.J. , 796 (1995).[12] T. Stepanova et al., Curr. Biol. , 1023 (2010).[13] D. J. Odde and H. M. Buettner, Biophys. J. , 1189(1998). [14] See Supplemental Material (available with the journal)for details of calculations.[15] T. Antal, P. L. Krapivsky, S. Redner, M. Mailman andB. Chakraborty, Phys. Rev. E , 041907 (2007).[16] R. Padinhateeri, A. B. Kolomeisky and D. Lacoste, Bio-phys. J. , 1274 (2012).[17] B. Houchmandzadeh and M.Vallade, Phys. Rev. E ,6320 (1996).[18] E. Jobs, D. E. Wolf and H. Flyvbjerg, Phys. Rev. Lett. , 519 (1997),[19] B.S. Govindan and W. B. Spillman, Phys. Rev. E ,032901 (2004).[20] B. M. Mulder, Phys. Rev. E , 011902 (2012).[21] T. E. Holy and S. Leibler, Proc. Natl. Acad. Sci. USA , 5682 (1994).[22] M. Gopalakrishnan and B. S. Govindan, Bull. Math. Biol. , 2483 (2011).[23] V. Yadav and S. Mukherji, Phys. Rev. E , 062902(2011).[24] C. Tischer, P. R. Wolde and M. Dogterom, Biophys. J. , 726 (2010).[25] M. Abramowitz and I. A. Stegun, Handbook of Mathe-matical Functions (Dover, New York,1972). [26] V. Varga et al., Nat. Cell. Biol. , 957 (2006); Cell 138,1174 (2009).[27] B.S. Govindan, M. Gopalakrishnan and D. Chowdhury,Europhys. Lett. , 40006 (2008).[28] G. A. Klein, K. Kruse, G. Cuniberti and F. J¨ulicher,Phys. Rev. Lett. , 108102 (2005).[29] A. Melbinger, L. Reese and E. Frey, Phys. Rev. Lett. , 258104 (2012).[30] H. S.Kuan and M. D. Betterton, Phys. Biol. , 036004(2013).[31] C. Tischer, D. Brunner and M. Dogterom, Mol. Syst.Biol. , 250 (2009).[32] M. F. Carlier, R. Melki, D. Pantaloni, T. L. Hill, and Y.Chen, Proc. Natl. Acad. Sci. USA , 5257 (1987).[33] R. Melki, M. F. Carlier and D. Pantaloni, The EMBOJournal, , 2653 (1988).[34] H. Obermann, E. M. Mandelkow, G. Lange, and E. Man-delkow, J. Biol. Chem. , 4382 (1990).[35] E. M. Mandelkow and E. Mandelkow, Cell Motil. Cy-toskeleton , 235 (1992).[36] A. Marx and E. E. Mandelkow, Eur. Biophys J. , 405(1994).[37] D. T. Gillespie, J. Phys. Chem. , 2340 (1977).[38] R. V. Skibbens, V. P. Skeen and E. D. Salmon, J. CellBiol. , 859 (1993).[39] J. R. Mcintosh, E. L. Grishchuk and R. R. West, Annu.Rev. Cell Dev. Biol. , 193 (2002).[40] T. M. Kapoor, M. A. Lampson, P. Hergert, L. Cameron,D. Cimini, E. D. Salmon, B. F. McEwen, and A. Khod-jakov, Science , 388 (2006).[41] ] A. C. Amaro, C. P. Samora, R. Holtackers, E. Wang,I. J. Kingston, M. Alonso, M. Lampson, A. D. McAinsh,and P. Meraldi, Nature Cell Biol. , 319 (2010).[42] N. Mchedlishvili, S. Wieser, R. Holtackers, J. Mouysset,M. Belwal, A. C. Amaro, and P. Meraldi, J. Cell. Sci. , 906 (2012) Appendix A: PDFs in transformed space
On taking Laplace transform of Eq.18 and Eq.19 withrespect to x , we get ∂ ˜ P (cid:48) ( p, τ ) ∂τ + [ v g p + ν c ( τ )] ˜ P (cid:48) ( p, τ ) = δ ( τ )[ ν r ˜ P ( p ) + J ] , (A1)˜ P ( p ) = J − Φ( p ) v s p − ν r , (A2)where Φ( p ) = ∞ (cid:82) ν c ( τ ) ˜ P (cid:48) ( p, τ ) dτ . Substituting for ˜ P ( p )from Eq.A2 in Eq.A1, ∂ ˜ P (cid:48) ( p, τ ) ∂τ + [ v g p + ν c ( τ )] ˜ P (cid:48) ( p, τ ) = δ ( τ ) (cid:20) J v s p − ν r Φ( p ) v s p − ν r (cid:21) . (A3)Eq.A3 is a first order non-homogeneous partial differ- ential equation in τ . We assume a solution of the form,˜ P (cid:48) ( p, τ ) = h ( p, τ ) exp (cid:18) − v g pτ − τ (cid:90) ν c ( τ (cid:48) ) dτ (cid:48) (cid:19) . (A4)After substituting Eq.A4 in Eq.A3 we get, h ( p, τ ) = J v s p − ν r Φ( p ) v s p − ν r . (A5)Therefore the total solution to Eq.A3 is given by˜ P (cid:48) ( p, τ ) = (cid:20) J v s p − ν r Φ( p ) v s p − ν r (cid:21) × exp (cid:18) − v g pτ − τ (cid:90) ν c ( τ (cid:48) ) dτ (cid:48) (cid:19) . (A6)Using the above equation we get,˜Φ( p ) = J v s pζ ( p ) v s p − ν r (1 − ζ ( p )) , (A7)where ζ ( p ) = ∞ (cid:90) ν c ( τ ) exp (cid:18) − v g pτ − τ (cid:90) ν c ( τ (cid:48) ) dτ (cid:48) (cid:19) dτ. (A8)The function η ( p ), which appears in Eq.23 of the maintext is defined as, η ( p ) = ∞ (cid:90) exp (cid:18) − v g pτ − τ (cid:90) ν c ( τ (cid:48) ) dτ (cid:48) (cid:19) dτ. (A9)As mentioned in the main text, all further analysis isdone by putting ν r = 0. In this case, Eq.23 and Eq.24 inmain text become, ˜ P ( p ) = J η ( p ) (A10)and ˜ P ( p ) = J − ζ ( p ) v s p . (A11) Appendix B: Length distribution a. Exponential aging
In this case, ζ ( p ) and η ( p ) are evaluated from the inte-grals given by Eq.A8 and Eq.A9 respectively, with ν c ( τ )substituted from Eq.1. ζ ( p ) = ν max c v g p + ν max c − v g pv g p + ν max c e α × α − vgp + ν max cλ γ (cid:18) v g p + ν max c λ + 1 , α (cid:19) , (B1)3 η ( p ) = 1 λ e α α − vgp + ν max cλ γ (cid:18) v g p + ν max c λ , α (cid:19) . (B2)Using the recurrence relation [25], γ ( a + 1 , b ) = aγ ( a, b ) − b a e − b , (B3)we get a relation between ζ ( p ) and η ( p ) as given below, ζ ( p ) = 1 − v g p η ( p ) . (B4)Therefore, Eq.A11 can be rewritten as˜ P ( p ) = v g v s J η ( p ) . (B5)Using the inverse Laplace transform, L − [ a − p γ ( p, a )] = e − ae − x , (B6)we find the Laplace inverse of Eq.A10 and Eq.B5 with η ( p ) substituted from Eq.B2. The resulting expressionfor the total PDF is given by Eq.28 in the main text. b. Linear aging The integrals given by Eq.A8 and Eq.A9 are evalu-ated using Eq.7. The expressions for ζ ( p ) and η ( p ) are,respectively, given as, ζ ( p ) = 1 − (cid:114) π β v g p e v gp β erfc (cid:20) v g p √ β (cid:21) , (B7) η ( p ) = (cid:114) π β e v gp β erfc (cid:20) v g p √ β (cid:21) . (B8)From Eq.B7 and Eq.B8, it can be seen that Eq.B4 holdstrue in the linear model as well. A useful inverse Laplacetransform to invert Eq.A10 and Eq.A11 in this case isgiven as [25], L − (cid:20) e a p erfc[ ap ] (cid:21) = e − x a a √ π (B9)and the final expression for the P ( x ) is given by Eq.34in the main text. Appendix C: Length autocorrelation c. Constant catastrophe: the limit λ → ∞ Conditional probabilities obtained from Eq.42 andEq.43 after taking Laplace transform with respect tospace and time, are given by G j ( p, s | x ) = e − px [( s − v s p ) δ j + ν r ] + v s ( s − v s p ) G j ( x = 0 , s | x )( s − v s p + ν r )( s + v g p + ν max c ) − ν max c ν r , (C1) G j ( p, s | x ) = e − px [( s + v g p ) δ j + ν max c ] − v s ( s + v g p ) G j ( x = 0 , s | x )( s − v s p + ν r )( s + v g p + ν max c ) − ν max c ν r . (C2)We use the convergence property of the conditional PDFsto find G j ( x = 0 , s | x ). Thus, at the boundary, expres-sions for the Green’s functions are given by G j ( x = 0 , s | x ) = ( s + v g β s ) δ j + ν max c v s ( s + v g β s ) e − β s x , (C3)where β s = − A ( s )2 + (cid:113) A ( s ) + B ( s ), with A ( s ) = s ( v s − v g ) + v s ν max c − ν r v g v g v s , (C4) B ( s ) = s ( s + ν max c + ν r ) v g v s . (C5) In the limit ν r = 0, Eq.C3 becomes, G j ( x = 0 , s | x ) = s ( v s + v g ) δ j + v s ν max c sv s ( v s + v g ) e − svs x . (C6)In the steady state, the PDFs take the simple form, P ( x ) = J v g e − ν max c xvg , (C7) P ( x ) = J v s e − ν max c xvg , (C8)where J = ν max c v s / ( v g + v s ) . d. Linear aging Solutions of Eq.13 and Eq.14 for the linear model intransformed space are respectively given as, G j ( p, s | x ) = e − px [( s − v s p ) δ j + ν r ] + v s ( s − v s p ) G j ( x = 0 , s )[ s − v s p + ν r (1 − ζ ( p, s )] η ( p, s ) , (C9) G j ( p, s | x ) = e − px [ δ j − ζ ( p, s )( δ j − − v s (1 − ζ ( p, s )) G j ( x = 0 , s )[ s − v s p + ν r (1 − ζ ( p, s )] , (C10)with ζ ( p, s ) = 1 − (cid:114) π β ( s + v g p ) e ( s + vgp )22 β erfc (cid:20) s + v g p √ β (cid:21) (C11)and η ( p, s ) = (cid:114) π β e ( s + vgp )22 β erfc (cid:20) s + v g p √ β (cid:21) . (C12) Using similar procedure mentioned above, we fix G j ( x =0 , s | x ) in the limit ν r = 0 as G j ( x = 0 , s | x ) = e − sx vs [ v s − (1 − δ j ) s ( vs + vg ) ψ ( s )] sv s ( v s + v g ) ψ ( s ) , (C13)where ψ ( s ) = (cid:114) π β exp (cid:18) ( v s + v g ) s v s β (cid:19) erfc (cid:20) ( v s + v g ) sv s √ β (cid:21) ..