Length Regulation Drives Self-Organization in Filament-Motor Mixtures

Cytoskeletal networks form complex intracellular structures. Here we investigate a minimal model for filament-motor mixtures in which motors act as depolymerases and thereby regulate filament length. Combining agent-based simulations and hydrodynamic equations, we show that resource-limited length regulation drives the formation of filament clusters despite the absence of mechanical interactions between filaments. Even though the orientation of individual remains fixed, collective filament orientation emerges in the clusters, aligned orthogonal to their interfaces.

Cytoskeletal networks form complex intracellular structures. Here we investigate a minimal model for filament-motor mixtures in which motors act as depolymerases and thereby regulate filament length. Combining agent-based simulations and hydrodynamic equations, we show that resourcelimited length regulation drives the formation of filament clusters despite the absence of mechanical interactions between filaments. Even though the orientation of individual remains fixed, collective filament orientation emerges in the clusters, aligned orthogonal to their interfaces.
The microtubule cytoskeleton plays an important role in numerous cellular functions such as intracellular transport and cell division [1,2]. These complex processes are based on active remodelling of the cytoskeletal structure [3], which is mediated by the interaction of microtubules with a variety of microtubule associated proteins (MAPs) [4][5][6]. In addition to generating forces between microtubules [7], MAPs play an important role in regulating the length of individual microtubules by affecting the rates of their polymerisation kinetics from tubulin subunits [8][9][10][11]. How forces affect the large-scale selforganization of microtubules has been studied in detail both theoretically and experimentally [12][13][14][15][16][17][18]. In contrast, the role of length regulation has only been investigated in the context of individual filaments [9,[19][20][21][22][23][24][25][26], or of a globally accessible pool of constituents (tubulin and MAPs) [27][28][29][30]. However, recently the focus of interest is shifting to their role in many filament systems, as there is increasing experimental evidence that this regulatory function, in combination with the local availability of MAPs and tubulin, plays an essential role in the self-organization, scaling and maintenance of microtubule structures [31][32][33][34][35][36][37]. It remains an important open question how the interplay and spatial redistribution of these resources through cytosolic diffusion and transport along microtubules affects the self-organization of the microtubule cytoskeleton [38][39][40].
Here, we approach this question by studying the collective motor-filament dynamics in planar geometry with limited resources of tubulin units and molecular motors. These cytosolic resources are spatially redistributed by diffusion while filament-bound motors additionally move uni-directionally towards the filament plus-end where they act as depolymerases (Fig. 1). We show that the interplay of motor-catalyzed depolymerization and local resource availability leads to self-organization of the filament assembly into aster-like patterns. Those patterns show colocalisation of microtubule plus ends and polarity sorting at the interfaces of emerging filament clusters.
Model. -We propose an agent-based model that builds on current experimental and theoretical studies addressing the resource-limited length regulation of a sin- gle microtubule by the kinesin-8 homologue Kip3 from Saccharomyces cerevisiae. [26]. Specifically, we study filament dynamics in planar geometry containing a finite number of tubulin units (N T ), molecular motors (N M ), and filaments (N F ); see Fig. 1(b). Each individual filament i ∈ 1, . . . , N F is represented by a directed rigid rod with fixed minus-end position b i and fixed orientation θ i ∈ [0, 2π), which are drawn randomly from uniform distributions. The lengths l i (t) of the individual filaments are dynamic variables that change by polymerization kinetics at the plus end. When filaments shrink to zero length, they are assumed to regrow form the same minus-end position and with the same orientation; filament shrinkage to zero length, though, rarely occurs. In the cytosol, both motors and tubulin units diffuse freely with diffusion constants D M and D T , respectively. Cytosolic motors can bind with rate k on to any point that is within the binding radius r M along a filament; for details see Supplemental Material Sec. SII. Filament-bound motors move towards the filament plus-end at speed v m , where they catalyse filament depolymerization at rate δ [see Fig. 1(b)]. Upon depolymerization, the filament length is reduced by one tubulin unit (of length a) and both the plus-end-bound motor and the associated tubu-lin unit are released into the cytosol. Cytosolic tubulin within a distance r T of a filament plus-end, binds to it at the rate γ, increasing filament length by a [ Fig. 1(d)].
Single-filament dynamics. -Consider a cytosolic volume V 0 , containing a single filament and a finite number of tubulin units ρ T V 0 and motor proteins ρ M V 0 . For now, we assume for simplicity that the cytosolic concentrations c M and c T are spatially uniform; this assumption is relaxed when we discuss a spatially extended system with many filaments. The length change of the filament is determined by the antagonism between polymerization and depolymerization kinetics ∂ t l(t) = v g − v s with the growth and shrinkage velocity given by v g = a c T γ and v s = a m + (t) δ, respectively, where m + (t) denotes the density of motors bound to the plus end [22,25].
For biologically relevant parameter ranges, the motor dynamics are fast compared to filament growth and shrinkage [10,26]. This separation of time scales implies that for a given filament length, the motor density can be assumed to be in a quasi-steady state, where the total attachment flux of motors onto the filament, j on = k oncM l, and the off-flux due to depolymerization events at the plus end, j off =ṽ s /a, are in balance; quasi-steady states are indicated by a tilde. Thus, the depolymerization velocityṽ s = a k on lc M is determined by the cytosolic densityc M , which in turn is related to the filament-bound motor numberM via mass conservation ρ M V 0 =c M V 0 +M . In steady state, the filament-bound motor density exhibits an antenna profilẽ m(s) = koncM vm s [41], which is inferred from the transport equation ∂ t m(s, t) = −v m ∂ s m(s, t) + k on c M (t) [42,43], implyingM = koncM 2vm l 2 . Taken together, we find a relation for the shrinkage velocity in terms of the filament length and the total motor concentratioñ v s (l, ρ M ) = a k on lc M = a k on l ρ where we have defined the characteristic length scale l c := 2v m V 0 /k on . For filament lengths l < l c , the shrinkage velocity increases with the filament length, as would be expected with unlimited motor resources and has been observed experimentally [10,19]. Increasing the filament length beyond l c leads to a depletion of the cytosolic motor pool and thereby a decreased on-flux k oncM l. According to the flux balance condition, this reduces the off-fluxṽ s /a and thus the shrinkage velocity, so thatṽ s ∼ 1/l for l l c [see Fig. 2(a)]. The growth velocity v g can be written in terms of filament length l and total tubulin density ρ T using tubulin mass conservation (ρ T V 0 = c T V 0 + l/a) as v g (l, ρ T ) = γ (ρ T a − l/V 0 ). Together with the balance of filament growth and shrinkage, v g (l, ρ T ) =ṽ s (l, ρ M ), this determines the steady state length l * (ρ T , ρ M ) [ Fig. 2(a)] [44].
self-organization in a spatially extended system. -How does the length regulation of individual filaments play out in a spatially extended system where resources are shared by cytosolic diffusion between many filaments? In the limiting case where the cytosolic concentration is slowly varying on the scale of the (typical) filament length, the filaments can be treated as point-like objects carrying a tubulin mass proportional to their length l(x, t). The single filament dynamics can then immediately be generalized to a local length regulation dynamics with the local shrinkage speed given in terms of the local quasi-steady state approximation for the cytosolic motor density,ṽ s ( Eq. 1). The dynamics of the cytosolic tubulin concentration is governed by a reaction-diffusion equation where the local polymerization kinetics induces sinks and sources of cytosolic tubulin; here V 0 = V /N F denotes the cytosolic volume associated with a single filament. The total motor density is redistributed by cytosolic diffusion where we again used the local quasi-steady state approximation for the cytosolic motor densityc M (x, t). Taken together, Eqs. (2)-(4) form a closed set governing the system's dynamics in the long-wavelength limit. The stability of a spatially uniform state (l * , c * T ,ρ M ) against spatial perturbations can be studied using a linear stability analysis (see Supplemental Material Sec. SIV for details). Figure 3 shows a typical dispersion relation σ(q) for the eigenvalue with the largest real part and the ensuing stability diagram as a function ofρ M andρ T . For ρ T >ρ crit T (ρ M ) there is a band of unstable Fourier modes q ∈ (0, q max ) extending to long wavelengths (q → 0). It is instructive to first consider the particular limit of well-mixed cytosolic tubulin. Then, the marginal mode q max reduces to q 2 max = −ρ M ∂ lṽs | l * /(D McM ). This implies that there is an instability against spatial perturbations (lateral instability) only if ∂ lṽs | l * < 0. Moreover, the band of unstable modes narrows with increasing D M , showing that cytosolic motor diffusion attenuates the lateral instability. Relaxing the assumption of well-mixed cytosolic tubulin, i.e., explicitly accounting for tubulin diffusion, yields the critical ratio of diffusion constants D crit T /D M ≈ γ/(ak on V 0ρM ) in the limit l * l c . For physiological parameters, we find that the there is a lateral instability if the average number of motors per filament satisfiesρ crit M V 0 > 0.57 D M /D T . This condition is well met for biologically relevant motor concentrations as D M /D T ∼ 1/6 (see Supplemental Material Sec. SI).
The feedback mechanism underlying the lateral instability can be explained in terms of a mass-redistribution instability [45][46][47]. To simplify the argument, we set D M = 0 for the moment so that the total motor density remains invariant under the dynamics and therefore spatially uniform ρ M =ρ M , cf. Eq. (4). Consider now a small perturbation δl(x) added to the homogeneous state l * , while keeping the cytosolic tubulin concentration c T (x) = c * T initially constant [ Fig. 2(b,c)]. Since then v g = a γ c T initially remains uniform, the effect of δl(x) on the net growth velocity v = v g −ṽ s depends on the slope of the shrinkage velocity at l * . For ∂ lṽs | l * < 0, filaments grow (shrink) when they are long (short). This leads to an decrease (increase) of the cytosolic tubulin concentration [arrows in Fig. 2(b,c)] creating gradients in the cytosolic tubulin concentration that drive diffusive transport of tubulin mass towards regions of increased filament length. Since this tubulin mass redistribution leads to an increase of v g in regions where δl > 0, it promotes further filament growth there, i.e., the initial spatial perturbation δl(x) is amplified [ Fig. 2(b)]. In contrast, if the regulatory kinetics is such that the shrinkage velocity increases with filament length (∂ lṽs | l * > 0), the effect is opposite. Cytosolic tubulin diffusion then redistributes the tubulin mass to regions with shorter filaments, counteracting the original disruption. Taken together, one finds the condition ∂ lṽs | l * < 0 for a spatial instability that is driven by free tubulin diffusion, in accordance with the result of the linear stability analysis.
The above reasoning also explains why cytosolic diffusion of motor proteins mitigates the lateral instability. Regions with short filaments contain fewer binding sites for motors and thus the cytosolic motor concentration is high there. The opposite holds for regions with long filaments. This creates gradients, and thereby diffusive fluxes, of motors towards regions of long filaments. The resulting diffusive influx of motors increases the rate of filament depolymerization there and thus counteracts the instability driven by tubulin diffusion.
Agent-based simulations. -To study the spatiotemporal dynamics above the critical tubulin concentrationρ crit T (ρ M ) we perform agent-based simulations. While Eqs. (2)-(4) capture well the initial dynamics at the long-wavelength instability, they fail to give the correct dynamics once gradients begin to emerge at the small length scales (see Supplemental Material Sec. SIV). What this continuum theory lacks are effects due the spatial extent of the filaments which includes motor binding along the length of filaments as well as motion of each filament plus-end due to polymerization kinetics. Figure 4 shows a time sequence obtained from the simulations (see also Movie S1). First, regions with short (depletion zones) and long (clusters) filaments are formed, which corresponds to the initial dynamics described by the mass redistribution instability (Fig. 4, t = 30 min). Moreover, filament plus-ends start to accumulate at the interface between these zones. As the dynamics progresses, the depletion zones grow in size and the interfaces sharpen (Fig. 4, t = 60 min). At this time point, the filament-length distributions match on a qualitative level with experimental measurements [26] (see Supplemental Material Sec. SIII). Subsequently, the high density regions segregate into individual large scale filament clusters, which are characterized by sharp boundaries and strong co-localization of filament plusends at their periphery (Fig. 4, t = 180 min). This colocalization is caused by the movement of the filaments' plus-end due to polymerization dynamics that is directed to zones where the net growth rate changes sign, namely cluster interfaces. In the long run, the large filament clusters grow at the expense of the smaller ones, until eventually only a single cluster remains, which then develops into an aster-like structure (Fig. 4, t = 700 min).
Inside the clusters, the filaments exhibit net polar order that is aligned along tubulin-density gradients, i.e., orthogonal to the cluster boundaries. This is because the plus ends localized there belong predominantly to filaments whose minus end lies within the cluster's interior, implying an orientation orthogonal to the boundary on average ( Fig. 5 and Movie S2).
To quantify this effect, we monitor the density gradient ∇ρ T , the local net polarity p, and the angle θ enclosed between these vectors. Figure 5(d) shows the time evolution of the histogram P(θ) of the angle θ weighted by the product of the magnitudes of ∇ρ T and p to highlight the alignment of filaments near the cluster boundaries. The initially uniform distribution P(θ) evolves quickly into a peaked distribution centered around zero-indicating the onset of polar order-and subsequently sharpens slowly; see also snapshots in insets of Fig. 5(a). The onset of this polar order occurs simultaneously with the massredistribution instability, as can be seen from the comparison of the spatial averages p · ∇ρ T and |∇ρ T | 2 , which are coarse-grained measures of filament orientation and density gradients, respectively; Fig. 5(e). The polar order leads to advective flow of filamentbound motors out of clusters, which is balanced against diffusive influx caused by gradients in cytosolic motor . Fast binding of motors inside clusters together with advective motor transport leads to the depletion of motors in the cluster interior and the formation of sharp gradients in cytosolic motor concentration. Those gradients help to maintain the filament plus-end localization at the interfaces: Plus ends that protrude beyond the interface are subjected to an increased on-flux of motors, causing the filaments to shrink back. Conversely, plus-ends within the cluster are subjected to a reduced motor on-flux, causing them to grow towards the interface. Finally, the sharp cytosolic gradient leads to a shrinkage velocityṽ s that is independent of filament length because motor attachment occurs only in a narrow band at the interface. This is a collective effect and in contrast to the length regulation of a single filament, which is strongly length-dependent [cf. Eq. (1)]. The size regulation of clusters and a quantitative analysis of the final, aster-like, stationary state will be presented in a forthcoming publication [48].
Discussion. -Commonly, the spatial self-organization of a filament arrangement is attributed to motor proteins that reorient and move filaments by mechanical forces, such as dynein or kinesin-5 [7,18,49,50]. Here, we have shown that microtubule length regulation (through kinesin-8) in combination with resource limitation can lead to aster-like spatial patterns. The underlying instability is driven by diffusive redistribution of cytosolic tubulin mass. Such a mass-redistribution instability is a potential candidate to explain the emergent selforganization observed in cell extracts [51]. Notably, this self-organization is heralded by spatial patterns that emerge in the tubulin density, which have comparable morphology and wavelength (∼ 100 µm) as those we observe in our simulations (cf. Fig. 4). We also expect that our theory for resource-limited filament length regulation can be used to investigate heterogeneous growth dynamics in systems where spatial heterogeneities in filament length and/or density are imposed, e.g. by experimental design [39] or by upstream gradients [52,53]. In the system studied here, inhomogeneous filament polymerization occurs spontaneously, leading to inhomogeneous filament densities (lengths).
From a broader perspective, the conceptual model investigated here is in itself an interesting active matter system exhibiting self-organized patterns, polarity sort-ing, and coarsening. Combining (active) length regulation with mechanical interactions of filaments will be an important avenue for future research.
We would like to thank Henrik Weyer, Isabella Graf and Philipp Geiger for careful reading of the manuscript. We acknowledge financial support by the Deutsche Foschungsgemeinschaft through the Excellence Cluster ORI-GINS under Germany's Excellence Strategy (EXC-2094-390783311). 6

SI. PARAMETER ESTIMATION
Since experimental data are available for many of the parameters in our model [9,10,54,55], we can assign experimentally motivated numerical values to them. In the following we will show how to convert the experimentally measured rates into appropriate rates for our simulation. The experimental and the resulting model parameters are summarized in table I.
A typical microtubule consists of thirteen protofilaments. We simplify our analysis by considering only a single protofilament per microtubule, i.e., we neglect correlations between neighbouring protofilaments of a single microtubule. This is achieved by converting rates and resources to rates and resources per protofilament. Given a heightL z of our simulation environment, the total number of resources per protofilament is given by N M = ρ M L x L yLz /13 and N T =ρ T L x L yLz /13 for motor proteins and tubulin units respectively. We define L z :=L z /13 as the effective z-extent per protofilament, which corresponds to the effective height of our simulation box. For a given choice of the spatial extents of the simulation box the values for N T and N M will be determined from the experimentally given values for the densitiesρ T andρ M .
The motor attachment rate k on can be determined from the experimentally measured rate constant for motor attachment, which is given in units per concentration, per time, and per length. Varga et al. [10] report a value of k exp on = 0.4 nM −1 s −1 µm −1 = 0.66 µm 3 s −1 µm −1 , which was measured over a range of motor concentrations; note that we do not cancel one unit of µm to emphasize that k on is a rate per volume concentration and per length of the microtubule. The motor on-flux onto a filament of length l at a given motor concentration c M is than given by k exp on c M l. The motor concentration in experiments is measured in nM which is converted into a number densities as 1 nM ≈ 0.6 particles/µm 3 . Since Varga et al. used a TIRF setup, where motors can bind to and walk on approximately 5 protofilaments [56], the attachment rate per protofilament is given by k on = k exp on /5 ≈ 0.13 µm 3 s −1 µm −1 . For our agent-based simulations of the filament-motor mixture, we need to further convert this rate to a per-capita rate for a single motor protein. To this end, we first convert the experimental rate constant per volume concentration into a rate per area concentration by k 2D on = k exp on /L z . This is converted to a per-capita rate by specifying a reaction radius r M within which the motors can attach to the filament. Thus, the experimental value for the attachment rate constant per protofilament k exp on can be converted to a per-capita rate used in the simulations by k sim on = k exp on /(L z πr 2 M ). Given the experimental value for the attachment rate, this scaling relation gives us some freedom in choosing the effective height L z of the simulation box and the reaction radius r M . The height of the simulation box has to be chosen such that all concentrations (tubulin, motor and filament concentration) can be assumed to be well mixed in z−direction, i.e., L z λ c , where λ c is the wavelength of the initial instability (see Sec. SIV B). Moreover L z has to be chosen large enough such that the number of filaments and particles is high enough to prevent stochastic effects due to number fluctuations. As long as these constrains are fulfilled we are free to choose L z in a way convenient for our simulations. In particular, this means we choose L z small to keep 7 Experiment Theory Simulation Ref.

Motor parameters
Tubulin and filament parameters  . To convert experimental rates into rates required for our simulation we used rM = rT = 0.04 µm and Lz = 0.6 µm.
If not explicitly stated otherwise we used parameters as specified in the column Theory.
the particle numbers sufficiently low for the numerical simulations to be feasible. The reaction radius r M must be smaller than the average distance covered by a motor in the time interval ∆t by diffusion, where ∆t is the time increment of the simulation (see Sec. SII). Otherwise particles could cover an unphysically large distance in the time interval ∆t by successive attachment/ detachment events. The spontaneous polymerization rate γ of microtubules can be obtained similarly to the attachment rate. The polymerization velocity of microtubules has been measured to be v exp = 0.19 µm min −1 µM −1 [9], which corresponds to a polymerization rate constant γ exp = v exp /a = 0.38 µM −1 s −1 = 6.3 · 10 −4 µm 3 s −1 . Analogous to the case of the attachment rate constant, the polymerization rate constant depends on the volume concentration of tubulin units. We convert this into a per-capita polymerization rate by γ sim = γ exp /(L z πr 2 T ), where r T denotes the reaction radius of a tubulin unit.
For the simulation results shown in the main text, we choose r T = r M = 0.04 µm and L z = 0.6 µm. At average concentrationsρ T = 2.75 µM andρ M = 50 nM this results in a total number of N T ≈ 2.2 · 10 7 tubulin units and N M ≈ 4 · 10 5 motor proteins; as well as per-capita binding rates γ sim ≈ 0.22 s −1 and k sim on ≈ 40 s −1 µm −1 (see Table I). The diffusion constant of cytosolic tubulin was measured as D T = 6 µm 2 s −1 [55]. To the best of our knowledge, the cytosolic diffusion constant D M of Kip3 is not known. Since the molecular mass of Kinesin-8 is about 2-3 times larger than that of tubulin [57], we expect the diffusion constant of Kip3 to be of the order D M ∼ 1 µm 2 s −1 . As we discuss in Sec. SIV B below, the precise value of the motor diffusion is not relevant for the qualitative results of our model. Linear stability analysis predicts that slower motor diffusivities entail a shorter wavelength of the fastest growing mode (see Fig. S11). To keep the computational cost of our agent-based simulations reasonable, we choose D M = 0.5 µm 2 s −1 , which allows us to simulate a smaller system and thus keep the number of particles in the simulation volume manageable.

A. Single-filament simulation
To simulate the motor-mediated length regulation of a single filament, we follow a two-pronged approach: First, we employ the Gillespie algorithm [58] to simulate the full stochastic dynamics of motor proteins using a lattice gas model (TASEP), which has been shown to be a good model system for studying the motion of motor proteins on microtubules. [42,[59][60][61]. Since this exact simulation approach suffers from performance problems when studying many filaments coupled via a diffusive reservoir, we also implement an approximate simulation scheme in which filament-bound motors move deterministically and steric interactions between motor proteins are neglected. To guarantee its validity, the approximate simulation scheme is compared with the results of the exact Gillespie method in the relevant parameter range.
We develop our approximate simulation schema based on theoretic results of the TASEP-LK model [42,43]: In the limiting case of low densities, the dynamics of the mean motor proteins on the filament is described by the mean-field equation While this result follows strictly from full stochastic dynamics, it can also be obtained from deterministic dynamics where during the time interval ∆t the position s i of the motor proteins is updated as i.e. all motor proteins move ballistically with the same speed v m . Moreover, for sufficiently high depolymerization rates and low densities, TASEP-LK predicts that the depolymerization ratẽ v s is given by [26]ṽ Neglecting effects from particle exclusion, this reduces toṽ s = v m m(l) = a m + δ , i.e., instead of explicitly modeling each depolymerization event, one can simply move the motor forward and if s i > l depolymerize the filament. Both approximations have the key advantage that they can be performed in parallel in the case of a spatially extended system with many motor proteins. To ensure that several motors do not depolymerize the filament within one iteration step, the time step must be chosen sufficiently small. Throughout all of our simulations we choose ∆t = 0.01 s. With a motor speed of v m = 0.06 µms −1 , this means that 'double' depolymerization events will only occur if the distance between two consecutive motors is less then 6 · 10 −3 µm, which is less then the length a of a tubulin unit. Figure S6 shows a comparison between the filament length at steady state obtained by a full stochastic simulation and our approximate simulation scheme.
B. Spatially extended system

Filament dynamics
A filament in the spatially extended system is characterized by its minus end position b i , orientation θ i , and length l i . The positions and orientations of the minus ends are initially drawn from a uniform distribution and remain unchanged throughout the simulation. Free tubulin units are inizialized at random positions within the simulation box. The stochastic dynamics of the free tubulin units is implemented as follows: First, we check whether filament plus ends are within the distance r T of the tubulin unit. For each plus end (0 . . . k) within range, a reaction time τ k is drawn from an exponential distribution; τ k ∼ γ sim exp(γ sim t). We choose the smallest reaction time τ k that fulfills τ k < ∆t, where ∆t denotes the time increment of the simulation, and perform the respective growth event. If there is no filament plus end in range or if all reaction times are τ k > ∆t, we let the tubulin unit perform a free diffusive motion implemented by a Brownian dynamics algorithm [62]. This is the position x n i = (x n i , y n i ) of the tubulin unit is updated as Here r n are independent and identically distributed random variables with zero mean and A is an amplitude chosen such that the fluctuation dissipation theorem is satisfied [63]. Importantly, however, it is not necessary for the random numbers to be Gaussian-distributed [62]. Here, we choose uniform random numbers in the interval (−0.5, 0.5) as this has several implementationspecific advantages. Since the variance of uniformly distributed random variables is given by (r n ) 2 = 1/12 this results in A = √ 24D T ∆t. The pseudocode for the stochastic dynamics of cytosolic tubulin units is given in Algorithm 2.

Motor dynamics
The motors can be either filament-bound or free (cytosolic). We perform simulations of the stochastic dynamics of cytosolic motors analogous to that of cytocolic tubulin. First, we determine all potential binding partners, i.e., all filaments (0 . . . k) that intersect with a radius r M around the motor position x i = (x i , y i ) (see Fig. S7). Next, the chord length ∆l is calculated; for an illustration see Fig. S7(b). The per-capita attachment rate for a single motor protein attaching to the filament is then given by k sim on ∆l. Similar as for the free tubulin dynamics, the reaction times τ k ∼ k sim on ∆l exp(−k sim on ∆lt) are drawn from an exponential distribution and the reaction with the smallest reaction time satisfying τ k < ∆t is executed. If a reaction occurs, the motor starts at a random position within the chord length ∆l. If no reaction occurs, the motor protein performs free diffusion, which is implemented in the same way as for free tubulin units. The dynamics of filament-bound motors is implemented as discussed above for single filaments.

SIII. COMPARISON TO EXPERIMENTAL DATA
We choose the setup of our agent-based model comparable to the experimental system in Rank et al. [26]. The experiments in [26] were performed in a setup where GMP-CPP stabilized microtubules were pre-grown at a concentration ofρ T = 2 µM in a confined system of 100 µl. By varying the amount of time the microtubules were pre-grown the authors controlled In experiments it is hard to control the precise level of protein concentrations. For example, the concentration of "active" Kip3 could not be determined in [26] since the inactivation of Kip3 during the purification, snap-freezing and thawing process could not be quantified. In addition the authors expect Kip3 to form clusters, which would cause an additional amount of Kip3 that does not participate in the length regulation dynamics in particular for high motor concentrations. In addition to the uncertainties in the protein concentrations, there are measurement inaccuracies because, for example, short microtubules could not be distinguished from tubulin clusters in experiaments. When comparing experiments and agent based simulations we tried to account for this lack of resolution by only taking filament length l i > 0.5 µm into account. In particular as the protein concentrations that are actively involved in the length regulation process are not known in experiments a quantitative comparison between the experimental data and our agent based simulation is not possible. It is however possible to compare the results on a qualitative level.
In our agent-based simulations we observe a slow initial dynamics consistent with the theo- retical result of the linear stability analysis that predicts an initial growth rate of a perturbation which is on the order ∼ O(10 −4 s −1 ). Once the initial perturbation has grown nonlinear effects start playing a role and the dynamics speeds up significantly. Our agent based simulations show that the system has not yet reached its steady state after one hour of incubation with Kip3, it is rather in the initial phase of pattern formation (cf. Fig. 4 in the main text). The measured length distribution after 60 min of incubation with Kip3 therefor strongly depend on the initial condition (initial length distribution). In [26] the authors found that for narrow initial length distribution the filament length distributions did not change significantly from the initial length distribution after 60 min of incubation with Kip3 (cf. Fig. S20 in [26]). Those observations are consistent with our theory and observations from the agent based simulation as the initial dynamics is slow (in particular for low motor concentrations) as stated above. For broader initial length distributions the authors observed different behaviour at different concentrations of Kip3 (see upper panel in Fig. S8). At low motor concentrations no significant change in the length distribution was observed. At intermediate concentrations the measured length distribution was bi-modal and at high motor concentrations the microtubule length distribution was broad (∼ uniform). Those observations are consistent with measurements from our agent-based simulation (see lower panel in Fig. S8). Note the volume per microtubule V 0 was estimated to be 1.66 µm 3 in [26]. Here however we used V 0 = 0.5 µm 3 resulting in less resources per filament and therefor shorter microtubules (this was done for the shake of computational feasibility). We used slightly higher tubulin concentrations (ρ T ∼ 2.5 µM) in the agent based simulation to obtain similar length distributions. Moreover we find a comparable evolution of the measured microtubule length distribution; see Fig. S9.
Our agent-based simulations shows that the microtubule length distribution for narrow initial distributions and low motor concentrations also becomes broad but at significantly later times (t 60 min). Moreover our simulations suggest that the bi-modal distribution observed in experiments for intermediate motor concentrations is just a transient phenomenon. At later times we observe (in the agent-based simulation) that the length distributions for both intermediate and low motor concentrations becomes more reminiscent of what is observed for high motor concentrations (as long as the concentrations are chosen in a range where patterns form see Fig. S10(d)).
We are recognizing that the experimental data available are not sufficient to fully validate our model. A first experimental approach to validate our model would be to repeat the experiments described in [26] and to measure the length distributions over a longer period of time to verify the predictions by our agent based simulations, namely that the length distributions at all motor concentrations were pattern formation is observed become similar. To gain further progress in understanding spatiotemporal dynamics, experiments that spatially and temporally resolve tubulin density would be desirable.

SIV. POINT-LIKE FILAMENT APPROXIMATION
In this section, we provide a detailed analysis of the dynamics in the point-like filament approximation. For the reader's convenience, we repeat the equations that govern the dynamics in this approximation where l 2 c = 2v m V 0 /k on . These dynamics conserve the total numbers of tubulin units and motors The corresponding average densitiesρ T andρ M will be the main control parameters of interest in the following. We will also discuss the role of the diffusion constants D T and D M . We will first consider the homogeneous steady states, which are fixed points of the singlefilament dynamics. In particular, our analysis will show that these fixed points can be read off from a graphical construction in the l-c T phase plane and that there is a regime of bistability. Next, we will perform a linear stability analysis of the homogeneous steady states against spatial perturbations and discuss various limiting regimes and the role of the diffusion constants D T and D M . Finally, we will present numerical simulations of the equations that govern the dynamics in the point-like filament approximation. These simulations capture the patterns that initially emerge from the instability predicted from linear analysis. However, sharp gradients, on scales shorter than the filament lengths, rapidly emerge such that the underlying assumption of pointlike filaments is violated.

A. Homogeneous steady states
The homogeneous steady states are given by the fixed points of the polymerization kinetics determined by the local balance of polymerization and depolymerization and conservation of the total density of tubulin In the l-c T plane, the solutions to these equations are given by the intersections between the nullcline c T =ṽ s (l,ρ M )/(aγ) and the reactive phase space c T =ρ T − l/(aV 0 ); see Fig. S10(a,b). This graphical construction allows us to gain insight into the behaviour of the homogeneous steady states as a function of parameters. As long as the slope of the nullcline is larger than −1/(aV 0 ) for all l, there is only a single fixed point [ Fig. S10(a)]. If there is a section where the nullcline slope is more negative than −1/(aV 0 ), i.e., ∂ lṽs < −γ/V 0 , the reactive phase space can intersect the nullcline three times [ Fig. S10(b), case (iii)], giving three homogeneous steady states. The bistable region is delimited by saddle-node (also called fold or limit point) bifurcations where the reactive phase space is tangential to the nullcline, i.e. where the nullcline slope is ∂ lṽs = −γ/V 0 .
To determine the stability against these steady states against spatially homogeneous perturbations, we linearize Eq. (10a) under the mass conservation constraint c T + l/(aV 0 ) =ρ T which yields ∂ t δl = σ poly δl, with σ poly = γ/V 0 + ∂ l . Thus, the steady state that lies on the nullcline section with ∂ lṽs < −γ/V 0 is unstable against spatially homogeneous perturbations (empty disk) while the other two are stable (filled disks).
The locations of these limit point bifurcations can be estimated by simple approximations, which highlight the role of the various parameters for the location of the bistable regime in theρ M -ρ T diagram. Near the apex of the nullcline at l = l c [marked by a red, dashed line in Fig. S10(a,b)], its curvature is large such the slope ∂ lṽs reaches −γ/V 0 near the apex. Substituting l * = l c into the total density of tubulin and using aγc * T =ṽ s (l * ) yields This shows that there is a linear relation between the motor densityρ M and tubulin density at the apex,ρ apex T [red line Fig. S10(c)]. Note that this provides a good approximation of the upper edge of the bistable regime.
To estimate lower edge of the bistable regime, we approximate consider the limit l * l c . To lowest order, this yields l * ≈ aV 0ρT and ∂ lṽs ≈ ak on l c /l 2 . These approximations, substituted in to the criterion which is plotted as a teal line in Fig. S10(c). Together, the estimates Eq. (14) and Eq. (15) reveal how the boundaries of the bistable regime depend on the system parameters. Increasing the polymerization rate γ moves the bistable regime to lower tubulin densities and reduces its size while increasing the motor velocity v m move the bistable regime to higher tubulin densities. The attachment rate k on and volume per filament V 0 only affect the upper boundary of the bistable regime, moving it to higher tubulin densities when increased.

B. Linear stability analysis
Let us now analyse the stability of the homogeneous steady states against spatial perturbations. For a more compact notation, we combine the field variables into the vector u = (l, c T , ρ M ). The homogeneous steady states are then u * = (l * , c * T ,ρ M ). Writing the spatial perturbations in terms of Fourier modes where q denotes the wave vector, and linearizing Eqs. (10) for small δu q yields an eigenvalue problem for the growth rates σ(q) with the Jacobian (or stability matrix) In the following, we omit the subscript u * , with the understanding that all expressions are evaluated in the homogeneous steady state. For spatially uniform perturbations, q = 0, the last row of the Jacobian vanishes and the eigenvalues of J(0) are given by The two zero eigenvalues correspond to perturbations that would change the average total densitiesρ T andρ M , respectively. In a closed system with conservation of mass, however, these disturbances are not possible. Thus, the remaining eigenvalue σ 3 determines the stability to homogeneous perturbations that respect conservation of mass. One finds an instability for ∂ lṽs < − γ V0 , which confirms the geometric criterion (in the in the l-c T phase space) discussed above.
In general, the stability of a spatially uniform steady state against spatial perturbations with a wave number q is determined by the dispersion relation σ(q); a typical dispersion relation in the laterally unstable parameter regime is shown in Fig. 3(a) in the main text.
While the eigenvalues of the full 3×3 Jacobian Eq. (18) can in principle be found analytically, the resulting expressions are not insightful. We therefore determine the eigenvalues numerically. Notably, we find that the instability is always of long-wavelength type, i.e., the band of unstable modes extends to q → 0. Next, we will further analyze this long-wavelength limit. Further below, we will also consider the limit of well-mixed tubulin. In both these limits the resulting effective Jacobians are 2×2 matrices which allows us to find closed expressions for the instability criteria.

Long-wavelength limit
Since we numerically find that instability onset is always at long wavelengths, i.e. in the limit q → 0, we can determine the instability criterion by considering the limit of small wavenumbers q. In this long-wavelength limit, we can make the local quasi-steady state approximation where the local equilibrium (l * , c * T ) depends on the local total densities ρ T , ρ M . The dynamics for the tubulin density is obtained by adding the equations for l and c T , Eqs. (10)(a,b). Using the local quasi-steady state approximation then yields The respective long-wavelength (LW) Jacobian is given by and its eigenvalues read (as for any 2 × 2 matrix) Therefore, there is an instability (positive eigenvalue σ i ) either when the trace (σ 1 + σ 2 ) is positive or when the determinant (σ 1 · σ 2 ) is negative. Moreover, when the determinant is positive, det J LW > 0, the eigenvalues' real parts cross zero as a pair of complex conjugates when tr J LW crosses zero. This indicates an oscillatory instability (Hopf bifurcation) where the oscillation frequency is determined by the imaginary part of the eigenvalues. (Further away from the onset of instability, imaginary part vanishes for the fastest growing mode, so the instability loses it's oscillatory character; see Fig. 3(a) in the main text.) J LW (q) has a positive trace if We will discuss this instability criterion below. Since, det J LW > 0, the instability will be oscillatory near onset, as noted above.
A positive eigenvalue could also result from a negative determinant. However, numerically, we find that the determinant seems to be always positive, although we could not show this analytically. There is also a physical reasoning why the determinant should always be positive. First, observe that the diffusion constants can be factored out. The remaining term, which then determines the sign of the determinant is independent of the diffusion constants. Therefore, if this term were negative, it would imply an instability independent of the diffusion constants. This is at odds with physical intuition for a diffusion driven lateral instability, in particular, because sufficiently fast motor diffusion will always act to suppress lateral instability.
Let us now analyze the stability criterion Eq. (23). First consider the case D M = 0. Then the instability condition simply reads ∂ ρT c * T < 0, which can be understood as follows. When D M = 0, the total motor density remains constant under the dynamics and therefore spatially uniform (by choice of initial condition). The only dynamic variable in the long-wavelength limit is then the tubulin density, governed by where we used the chain rule for the second equality. This is a diffusion equation with effective diffusion constant D T ∂ ρT c * T . Thus, if ∂ ρT c * T is negative, there is an instability corresponding to effective "anti-diffusion." This is the basic mechanism underlying mass-redistribution instability in the long-wavelength limit [47]. This shows that the basic mechanism of the instability is the same as in protein-based pattern forming systems like the ones discussed in [45].
In the main text, we discuss the stability condition based on the nullcline slope ∂ lṽs . To relate the derivative ∂ ρT c * T to the nullcline slope ∂ lṽs , we apply the derivative ∂ ρT to the fixed point equations Eq. (13) which yields Thus, along sections of the nullcline where ∂ lṽs > −γ/V 0 (implying stability against homogeneous perturbations), the lateral instability condition ∂ ρT c * T < 0 is analogous to the slope condition ∂ lṽs < 0 discussed in the main text. In theρ M -ρ T plane, the lineρ apex T (ρ M ), marks the nullcline apex where ∂ lṽs = 0 [red line in Fig. S10(c)]. Above this line, the nullcline slope is negative, and therefore, the homogeneous steady state is laterally unstable. In other words, a sufficiently high tubulin density compared to the motor density is required for the instability to occur.
In the bistable regime, we have to distinguish between the two stable steady states. Notably, the nullcline slope is always negative for the large-l fixed point (see Fig. S10(c). Therefore, this fixed point always laterally unstable. In contrast the low-l fixed point is unstable only in a very narrow regime due to the high nullcline curvature near its apex. In the stability diagram in Fig. 3(b) in the main text, we show the stability of the large-l fixed point.
The above conditions for lateral instability are necessary and sufficient in the case D M = 0. Let us now turn to the case D M > 0. As we heuristically argued in the main text, motor diffusion generally counteracts lateral instability. Indeed, the stability threshold tr J LW = 0 obtained from the above linear stability analysis gives Substituting Eq. (25) for ∂ ρT c * T and using ∂ ρMcM =c M /ρ M [cf. Eq. (10d)], we can write this condition as determining the threshold value for the ratio of the diffusion constants D T /D M . For diffusivity ratios below this threshold, i.e. for too fast motor diffusion, the instability is suppressed [see red lines in Fig. S11]. In Fig. S10(d), contour lines show the instability threshold in theρ M -ρ T plane for several diffusivity ratios.
Notably, for sufficiently largeρ T , the threshold becomes independent ofρ T . The critical motor densityρ crit M in this regime is indicated by the he dashed green lines in Fig. S10(c,d). To obtain an estimate for this critical motor density D T /D M , we approximate Eq. (27) in the limit l * l c and solve forρ M :ρ This shows that the critical motor density is proportional to the diffusivity ratio D M /D T . A higher motor diffusivity requires a higher motor density for the instability to occur. At first glance, this may seem somewhat counterintuitive. The reason for this effect is that increasing the motor density increases the magnitude of the nullcline slope sinceṽ s ∝ ρ M . This, in turn, increases the growth rate of the tubulin-mass-redistribution instability [cf. Eq. (24)], which allows it to overcome the stabilizing effect of motor diffusion. Physically, a steeper nullcline slope means that gradients in the filament length lead to steeper gradients in cytosolic tubulin concentration.
Substituting the values for a, γ, and k on from Table I   Above, we have analyzed the long-wavelength limit, where the polymerization kinetics can be assumed to be in a local quasi-steady state. We now turn to the dynamics at short wavelengths, where cytosolic diffusion of tubulin is faster than the polymerization kinetics. Relaxation to the local steady state length happens at the rate σ poly = γ/V 0 + ∂ lṽs , which we derived in the stability analysis for homogeneous perturbations above. The rate of diffusive transport for modes with wavenumber q is given by D T q 2 . Thus, for wavenumbers q |σ poly |/D T , the cytosolic tubulin density can be assumed well-mixed. (A detailed discussion of this "reactionlimited regime" and the complementary "diffusion-limited regime" at large wavelengths is given in Ref. [47] in the context of mass-conserving two-component reaction-diffusion equations.) For spatial non-uniform perturbations (q = 0), a well-mixed (WM) cytosolic tubulin implies δc T = 0. Thus, the reduced Jacobian is obtained by removing the central row and column from the full Jacobian J, Eq. (18), and reads As above the eigenvalues of this 2×2 matrix can be obtained from its determinant and trace where we usedṽ s ∝c M ∝ ρ M [cf. Eq. (10d,e)] such that ∂ ρMcM =c M /ρ M and ∂ ρMṽs =ṽ s /ρ M . Thus, if −∂ lṽs > 0, there is a band of unstable modes q ∈ (0, q max ) with This shows that cytosolic motor diffusion suppresses a lateral instability on short length scales. In particular, the band of unstable modes vanishes in the limit D M → ∞. Therefore, an approximation in which both the cytosolic tubulin and the motors are assumed to be well-mixed will not reproduce the instability. Equation (32), derived under the assumption of well-mixed cytosolic tubulin, only approximates edge of unstable modes of the full Jacobian Eq. (18). This approximation is valid if q 2 max D T |σ poly |. Substituting the expressions and rearranging the terms yields the condition Comparing to Eq. (27) shows that the approximation is valid deep in the unstable regime, far 21 from the instability threshold.

C. Finite-element simulations
To complement linear stability analysis, we performed numerical simulations of the point-like filament dynamics, governed by Eqs. (10). Specifically, a finite element method implemented in the COMSOL Mulitphysics software was used (see point-like-filament-PDEs.mph). To ensure numerical stability, we add a small diffusion term (diffusion constant 1 · 10 −2 µm 2 s −1 ) to the dynamics of l(x, t). This is necessary because sharp interfaces emerge rapidly around the long-filament clusters that from the initial instability (Fig. S12). Because of these sharp gradients, the assumption of point-like filaments underlying these simulations is no longer valid. Still, it is instructive to discuss the dynamics in this regime as it shares some of the features with the agent-based simulations with spatially extended filaments.
The interfaces of clusters that emerge from the initial instability propagate such that the regions of long filaments become smaller. This propagation is driven by the diffusive influx of motors from the regions where filaments are short, and therefore, most motors are in the cytosol. The motors diffusive into the long-filament regions, where they rapidly attach near the interface. This drives depolymerization of filaments near the interface, causing the long-filament regions to shrink. The released tubulin units then diffuse in the cytosol and drive further growth of filaments in the long-filament regions. In fact, the propagation of interfaces is already indicated by the dispersion relation. As we discussed below Eq. (22), a positive determinant of J WM implies that the eigenvalues are a pair of complex conjugates near the onset of instability. In the dispersion relation, this means that σ(q) has a non-zero imaginary part near the zero crossing of its real part at q = q max [see Fig. 3(a) in the main text]. In a previous study on mass-conserving reaction diffusion systems, we found that the properties of interfaces can be inferred from the right edge of the dispersion relation [47]. Specially, non-zero imaginary part at q max indicates that the mode that defines the interface is propagating.
In the point-like filament approximation, interfaces will continue to propagate until the longfilament clusters have completely disappeared. Subsequently, new clusters will emerge from lateral instability. And these clusters will again collapse due to interface propagation, driven by diffusive flux of motors into the clusters. In contrast, in agent-based simulations with spatially extended filaments, advective transport of motors along filaments counteracts diffusive influx of motors into the clusters. This is because the net orientation of the filaments is aligned with the density gradients and leads out of the clusters. As a result of this advective motor transport, interface motion arrests eventually, thus producing the final, aster-like steady state structure (see Movie 1 and Fig. 4 in the main text). if r < k on l∆t (#cytosolic motors) then 10: s k =random real ∈ (0, l) if r < γ∆t (#cytosolic tubulin) then 15: l += a x i += √ 24D T ∆t · (random real ∈ (−0.5, 0.5)) 14: y i += √ 24D T ∆t · (random real ∈ (−0.5, 0.5)) Algorithm 3 Free motor dynamics if remembered filament ! = empty then 12: attach motor to remembered filament at random position within ∆l 13: set motor state to filament-bound 14: else 15: x i += √ 24D M ∆t · (random real ∈ (−0.5, 0.5)) 16: y i += √ 24D M ∆t · (random real ∈ (−0.5, 0.5))