Connecting macroscopic dynamics with microscopic properties in active microtubule network contraction

The cellular cytoskeleton is an active material, driven out of equilibrium by molecular motor proteins. It is not understood how the collective behaviors of cytoskeletal networks emerge from the properties of the network's constituent motor proteins and filaments. Here we present experimental results on networks of stabilized microtubules in Xenopus oocyte extracts, which undergo spontaneous bulk contraction driven by the motor protein dynein, and investigate the effects of varying the initial microtubule density and length distribution. We find that networks contract to a similar final density, irrespective of the length of microtubules or their initial density, but that the contraction timescale varies with the average microtubule length. To gain insight into why this microscopic property influences the macroscopic network contraction time, we developed simulations where microtubules and motors are explicitly represented. The simulations qualitatively recapitulate the variation of contraction timescale with microtubule length, and allowed stress contributions from different sources to be estimated and decoupled.


Introduction
Active matter is a class of materials held out of equilibrium by the local conversion of energy from a reservoir into mechanical work at the scale of the system's components arXiv:1706.10235v1 [q-bio.SC] 30 Jun 2017 [1]. Active matter can exhibit emergent behaviors, such as collective motion and pattern formation, on lengthscales much larger than the size of the system's constituents.
Here we consider cytoskeletal networks, which are living, active systems [2] composed of polar polymeric filaments and held out of equilibrium by molecular motor proteins which convert chemical energy into mechanical work. Cytoskeletal networks are responsible for a number of cell biological processes, including cell division and chromosome segregation. The active behaviors of cytoskeletal networks have been experimentally investigated in both simplified purified mixtures [3,4,5,6] and in complex systems in cells and cell extracts [7,8,9,10]. However, it remains poorly understood how the dynamics and architecture of these networks are shaped by the properties of their constituent filaments and motors.
The past several years have seen several studies of the dynamics and mechanics of microtubule (MTs) / motor protein suspensions (see [11] for a review). A predominant thread of experiments have focused on dense, nematically aligned phases of MT bundles and immersed layers which show dynamics driven by extension of material along the direction of orientational order [12]. It is believed that this extensile motion is driven by the polarity sorting of anti-aligned MTs by multimeric kinesin motors. Theoretical modeling of the material stresses produced by polarity sorting indeed yields active stress tensors that are anisotropic, and arise from ensemble averages of terms of the form αpp, where p is an MT orientation vector, and α < 0 gives extensile flows along the p direction [13]. Different dynamics can follow from different MT interactions, and clustering of MT ends is a generic mechanism leading to contraction [14]. This has been previously explored in other purified systems containing kinesin motor proteins [15,16]. In this work, as in our earlier work of stabilized MTs in meiotically-arrested Xenopus oocyte extracts [7] and in several related studies [17,18,19], we consider a spontaneously formed MT network whose dynamics is driven by minus-end bundling by dynein motors. This leads to an apparently isotropic contraction, and we proposed a continuum model [7], which quantitatively reproduced the experimental results. Crucially, the model helped explain why networks contract to a preferred final density regardless of sample geometry and motor concentration. The continuum model assumed that the active material stress tensors were isotropic and of the form f (ρ)I, with f < 0 a combination of contractile motor stresses and repulsive steric stresses.
Here, we explore the same system as in our previous work [7] further by varying the initial MT density and the length distribution of MTs and measuring the resulting effects on the contraction process. Consistent with the predictions of our active fluid model, we find that the final network density varies little across a wide range of conditions. Intriguingly, we also find a correlation between the timescale of the contraction process and the average length of MTs in the system. To gain insight into how this microscopic parameter influences the macroscopic contraction timescale, we developed simulations of the contraction process. These simulations recapitulate many key aspects of the experimental results, including the variation of contraction timescale with MT length, and allow access to many aspects of the system difficult to measure experimentally, such as the spatial distribution and nature of the material stress tensors. Taken together, these results further support the underlying assumptions and key predictions of our active fluid model, and provide a framework for future studies.

Preparation of Xenopus extracts
Extracts were prepared from freshly laid Xenopus oocytes as described previously [20]. Fresh extracts were sequentially filtered through 2.0, 1.2, and 0.2 µm filters before being flash frozen in liquid nitrogen and stored at -80 • C until use.

Preparation of microfluidic devices
Microfluidic device templates were designed using AutoCAD 360 and Silhouette Studio software. Device templates were cut from 125 µm tape (3M Scotchcal) using a Silhouette Cameo die cutter, and adhered to petri dishes to create a master. PDMS (Sylgard 184, Dow Corning) was mixed at the standard 10:1 ratio, poured onto masters, degassed under vacuum, and cured overnight at 60 • C. The cured devices were then removed from the masters and inlets and outlets were created using biopsy punches. Devices and coverslips were treated with air plasma using a corona device, bonded, and loaded with passivation mix composed of 5 mg/mL BSA and 2.5% (w/v) Pluronic F-127 before overnight incubation at 12 • C.

Bulk contraction assay
The bulk contraction assay was performed as previously described [7]. Briefly, Alexa-647 labeled tubulin was added to 20 µL Xenopus extract at a final concentration of ≈1 µM. Then, 0.5 µL of taxol suspended in DMSO was then added to the extract to the indicated final concentration. Extracts were then loaded into passivated microfluidic devices, sealed with vacuum grease, and imaged using spinning disk confocal microscopy (Nikon TE2000-E microscope, Yokugawa CSU-X1 spinning disk, Hamamatsu ImagEM camera, 2x objective, Metamorph acquisition software). The time t = 0 corresponds to when imaging begins, typically < 1 min after taxol addition. Images were analyzed using ImageJ and custom MATLAB software. The (t) curves were fit using time points when (t) > 0.1.

Initial and final density estimation
In order to estimate network densities, we first assume that Alexa-647 labeled tubulin uniformly incorporates into MTs, and the overall concentration of tubulin is taken to be constant and equal to the previously measured value of 18 µM [21]. From this we can take, where M i is the tubulin mass in pixel i, β is a constant conversion factor between concentration in micromolar and fluorescence intensity, is the depth of the sample, A is the area of the pixel, V i = A is the volume corresponding to pixel i, N is the total number of pixels, I i is the measured intensity of pixel i, A is the area of pixel i, and I N is the measured intensity averaged over all pixels in the channel. From this we can infer, The intensity at each pixel in the network contains a contribution N i from polymerized tubulin and a contribution B i from monomeric, unpolymerized tubulin. The signal from monomeric tubulin is assumed to be constant and homogeneous throughout the channel. Thus, at the given time point where the background and network intensities are measured, the average concentration of polymerized tubulin in the network is given by, where N network is the number of pixels in the network at the time point, and I network is the intensity averaged over all pixels in the network.
To estimate the network density, the frame at the time closest to t = T c + τ was selected, where τ is the characteristic contraction time and T c is the lag time between the beginning of imaging and the beginning of network contraction [7]. This frame was corrected for inhomogeneous illumination and the camera's dark current, and analyzed as above to find ρ (T c +τ ). As we assume that no MTs are created, destroyed, or added to the network during the contraction process, the total mass of MTs in the network must be conserved. Thus, where V 0 , V Tc+τ , V F are the volume of the network at the beginning, at time T c + τ , and at the final state. Then, where W, H, L are the width, height and length of the network at different times, denoted by the subscripts as with the volume. We assume that, as the network is pinned at the channel's inlet and outlet, there is no change in volume along the channel's length, and thus L 0 = L Tc+τ . Combining Eqns. 1 and 2 in Results, .
We further assume that the change in the network's height follows the same functional form as the change in width and thus, For the final density of the network, similar reasoning can be used to show,

Measurement of MT length distributions
Stabilized MTs were dissociated from motor proteins and fixed as previously described [22]. MTs were allowed to assemble at room temperature for 30 minutes (unless otherwise noted) before 5 µL of extract was diluted into 50 µL of MT Dissociation Buffer (250 mM NaCl, 10 mM K-HEPES, pH 7.7, 1 mM MgCl 2 , 1 mM EGTA, and 20 µM taxol). After a 2 minute incubation, 100 µL of Fix Buffer (0.1% glutaraldehyde in 60% glycerol, 40% BRB80) was added and incubated for 5 minutes. 1 mL of Dilution Buffer (60% glycerol, 40% BRB80) was added to dilute the sample, and 2 µL of the diluted sample was spread between a slide and a 22×22 mm 2 coverslip. After waiting 30 minutes to allow the MTs to adhere to the coverslip, MTs were imaged using spinning disk confocal microscopy (Nikon TE2000-E microscope, Yokugawa CSU-X1 spinning disk, Hamamatsu ImagEM camera, 60x objective, Metamorph acquisition software). Active contours were fit to individual MTs using the ImageJ plugin JFilament [23], and MT lengths were determined from the contours using custom MATLAB software. For each taxol concentration, distributions of MT lengths were fit to a log-normal distribution to find the location parameter, µ, and the scale parameter, σ. In each case, the mode microtubule length for the log-normal distribution is given by

Results
To further investigate the dynamics of contracting MT networks in Xenopus oocyte extracts, we added taxol to the extracts to stabilize MTs, as previously described [7]. Extracts were loaded into rectangular microfluidic devices ( Figure 1A), sealed at the inlet and outlet using vacuum grease to prevent evaporation, and imaged at low magnification (Methods). Within minutes of taxol addition, the MT networks were found to undergo a spontaneous bulk contraction as shown in our previous work [7] ( Figure 1B, Movie 1). The size of the network along the width of the channel, W (t), was measured by identifying the pixels with high intensity (belonging to the network) and recorded as a function of time. Then the fraction contracted (t) was calculated, defined as where W 0 is the width of the channel, typically 0.9 mm. This process was repeated for varying final concentrations of taxol, and the (t) curves for each taxol condition were averaged together to produce master curves for each condition ( Figure 1C). The (t) curves were well fit by a saturating exponential function, where ∞ is the final fraction contracted, τ is the characteristic contraction time, and T c is a lag time between the beginning of imaging and the beginning of the the contraction process. Equation 2 was fit to the (t) curves for each experiment to extract values for the characteristic time, τ , and the final fraction contracted, ∞ . Values of τ and ∞ were averaged for each taxol condition ( Figure 1D). The characteristic timescale, τ , was found to increase with approximate linearity for taxol concentrations > 2.5µM, while the final fraction contracted, ∞ , decreased slightly with increasing taxol concentration. We next investigated how changing taxol concentration influenced both the initial density of the MT network, ρ I , and the final density of the MT network, ρ 0 . In systems of purified tubulin, the concentration of tubulin polymerized into MTs has been shown to increase with increasing taxol concentration [24]. Fluorescence intensity was used as a proxy for tubulin concentration, and was calibrated using previously measured values for the total tubulin concentration in Xenopus extracts (Methods). While the initial density of the MT network, ρ I , monotonically increased with increasing taxol concentration, the final density of the MT network displayed no obvious trend ( Figure 1E), and values of the final network density vary from the mean final density of ρ 0 = 49.6µM by less than 25%. This is consistent with previous results arguing that contracting MT networks in Xenopus extracts contract to a preferred final density [7]. Also, the increase of initial MT density, ρ I , with increasing taxol concentration appears to saturate for large taxol concentrations. Considering the three largest taxol concentrations investigated here (5 µM, 10 µM, and 25 µM), the initial MT densities was found to vary from the mean value by only 14% (std/mean), far less than the ≈ 2.5× increase seen in the characteristic contraction timescale, τ , over this same range. Thus, τ does not scale proportionally to the initial density of MTs, ρ I , and the change in contraction timescale must come from some other effect of changing taxol concentration.
Previously, it has been shown in clarified Xenopus extracts that the size of motor-organized taxol-stabilized MT assemblies, termed "pineapples", decreased with increasing taxol concentration, presumably due to a decrease in the average length of MTs [22]. Thus, we next sought to measure the length distribution of MTs in our system for varying taxol concentrations. We used a previously described method to dissociate and fix MTs [22] (Methods). Briefly, MTs were allowed to assemble for 5 minutes after taxol addition before being diluted into a dissociation buffer containing 250 mM NaCl, which dissociates motor proteins and other MAPs from the MTs. After a 2 minute incubation, this mix was diluted again into a fixation buffer containing glutaraldehyde. MTs were fixed for 5 minutes, further diluted, and thinly spread between a slide and a 22 mm 2 coverslip. MTs were allowed to adhere for 30 minutes before imaging. An active contour was fit to MTs in the images [23], and the length of the active contour was used as a measure for the MT length. This process was repeated for each taxol concentration.
Example histograms of MT lengths for each taxol concentration are shown in Figure  2A. Visually, the peak of the distribution shifts towards smaller values of MT length for increasing taxol concentration. Empirically, we find that the distributions of MT lengths are well fit by a log-normal distribution (Figure 2A, Inset). As fitting log-normal distributions to the measured histograms allows a more robust estimate of the mode MT length than using a purely empirical estimate, we fit the MT length distributions for each taxol concentration in order to find µ and σ, the two parameters of the lognormal distribution, and used these parameters to estimate the mode MT length for each condition. The mode MT length was found to decrease with increasing taxol concentration, varying by a factor of ≈ 1.7 across the taxol conditions measured ( Figure  2B).
One potential concern is that the MT length distribution may vary over the 45 minute timescale of the contraction experiments. Potentially, this could be due to a number of factors, including MT depolymerization, severing of MTs, or other causes. To address this issue, we repeated our length distribution measurement where MT dissociation and fixation began 45 minutes after taxol addition. We find close agreement between the MT length distributions measured either 5 minutes or 45 minutes after taxol is added, indicating that the length distribution is approximately constant over the timescale of the contraction experiments ( Figure 2C).

Model and Simulation
Our experimental results suggest that the change in contraction timescale observed for changing taxol concentrations may be due to changes in the MT length distribution. While our earlier continuum active gel theory for MT network contractions [7] captures the global contractile behavior of the network accurately, understanding the dependence of its parameters on microscale properties, such as MT length, is beyond the scope of the model. To resolve microscale changes and study their effects on the network's emergent properties, and to test the dependence of contraction timescale on MT lengths, we turn to simulation. Our simulation tracks the behavior of a suspension of fixed length MTs actuated by dynein motors which actively crosslink them and drag them through the fluid.

Model Description
In our modelling, we represent MTs as rigid spherocylinders that interact sterically and through motor protein coupling. In principle the fluid in which MTs are suspended couples their dynamics globally. Since solving the full Stokes problem numerically is prohibitively expensive we here neglect hydrodynamic many-body couplings. Instead we approximate the effects of the fluid by a local drag, that is, each MT acts as if it moved through a quiescent fluid. Note that for a dense and highly percolated network, one expects long range hydrodynamic interactions to be screened and hence we argue that they can be safely ignored for our purposes.
Suspensions of passive Brownian spherocylinders are well understood in both theory and simulation, and our numerical framework (in Appendix A) is based on well-known work [25,26,27]. What sets our simulations apart from the previous work is the presence of motor proteins, which have been previously considered with other simulation tools [28]. We model dynein motors as having a non-moving, fixed crosslinking end and a moving crosslinking end with stochastic binding and unbinding behaviors. The two ends are assumed to be connected by a Hookean spring with a rest length l 0 D 40 nm, which is the size of a force-free dynein. Collisions amongst dyneins or between dyneins and MTs are ignored for simplicity. In our model, dyneins have one fixed end, which stays rigidly attached to one MT throughout the simulation, and one motor end which stochastically binds, unbinds, and walks along other MTs. When the motor end is free, the dynein is transported with the MT bound to its fixed end without applying any force to it. Given the motor's small size, the orientation of such a motor is dominated by Brownian motion, i.e. is random. If a free motor head is closer than a characteristic capture radius r cap to another MT, it binds with a probability P b = ∆t/τ b , where ∆t is time step size, and τ b is the characteristic binding time. Should a free end be close enough to several MTs to bind them, its total binding probability P b stays the same, and one of the candidate MTs is chosen randomly. Finally, for the binding location we choose the perpendicular projection the dynein center onto the target MT. Bound motor ends walk towards the minus end of the MT they bind to, and exert a spring force and torque (F motor = κ D (l D − l 0 D )e D , T motor = r D × F motor ), where l D , κ D are the dynein's current length and its spring constant, respectively, e D is the direction of the Hookean spring force along the direction of dynein, and r D is a vector pointing from the center of mass of a MT to the point where the dynein binds and applies a Hookean spring force. Consistent with earlier work we impose a force velocity relation relation [13] and set the velocity and v = 0, otherwise. Here, v M and F stall are the motors free velocity and stall forces, respectively. All bound motor ends can stochastically unbind, with a characteristic time which becomes exp −F motor /F stall for F motor F stall , see [29], and diverges for F motor → F stall . Here τ 0 u is the force-free unbinding time. Finally, dyneins that reach the minus end do not immediately detach, but remain at the minus end, i.e. v = 0, keeping the same unbinding frequency. For this work we explored how key aspects of the microscopic model change the emergent behavior of the network. In particular, we varied the number of dyneins affixed to each MT, D/M , and the length distribution of MTs.

Setup and Parameters
The full model has a large number of parameters. However, many can be fixed from estimates in the literature. We summarize in Table 1 the parameters that are kept unchanged throughout the simulations presented here. We also fix the number of MTs in our simulations to be 10 µM, which is the approximate initial concentration found in experiments with 5 -25 µM of taxol. Following the experimental results, the MT length distribution is taken to be lognormal with parameters (µ, σ), with σ fixed at 0.5. Note that in simulation we use shorter MTs than in the experiments due to limitations in computational resources, see Appendix D.
In the simulations presented here, each MT carries the same number, D/M , of fixed dynein molecules, attached at random locations on the MT that carries them. Note that this is different from the assumptions made in earlier work [7], which has dyneins all fixed to the minus ends of MTs. This difference is motivated by the numerical observation that networks with all dyneins at minus ends will precipitate into asters, while networks with dyneins randomly affixed to MTs are more percolated and globally contract. It's unclear whether this difference is due to an artifact in the current simulations. Further exploring the effects of the localization of motors is a future research direction.
All of our simulations are performed in a simulation box of size 12 × 8 × 8 µm, with periodic boundary conditions along the long direction, mimicking the experimental geometry. The initial state is a random arrangement of MTs, in which all dynein motor  Figure 3 shows a representative example of the numerically observed contractions. After an initial phase of 30 s, in which MTs locally rearrange and cluster, bulk contraction driven by dyneins begins. It takes ∼ 100s to create a single, fully percolated network. The system has approached a steady state at 500 s, where the network has contracted to a cylindrical ribbon along the periodic direction of the simulation box.

Simulation Results
To compare the simulations to experiments, we calculate the fraction contracted from our simulations. Here, R(t) is the radius of a cylinder aligned with the periodic direction, and centered around the MT center of mass of the system portion which contains 50% of the total MT mass. By fitting (t) to the saturation exponential, e.g. Eqn 2, we now obtain a contraction timescale τ , and a final fraction contracted ∞ , which can be directly compared to experiments. As was done for fitting the experimental results, we ignore the very early stages of contraction ( (t) < 0.1), since they mostly mirror the initial rearrangement of the network and not the global contraction dynamics. Earlier work [7] has shown that τ depends on the initial system size. Since our numerical system is far smaller than the experimental one, we do not expect τ to quantitatively match experiments, but only to reproduce trends with varying parameters. Using the parameters from [7] for the continuum model, which was done with 2.5 µM taxol, and a system size of 8 µm, we expect a contraction time of τ ≈ 90 s, which is close to the time scales on which our numerical system contracts.
The final fraction contracted, ∞ , is easier to directly compare between experiments and simulation, since it does not depend on system size. Encouragingly, for the case of D/M = 1, the simulated contraction produces a similar final fraction contracted of ∞ ≈ 0.55, as in experiment, ( Figure 1D and Figure 4). Note that here only taxol concentrations of 5, 10, 25 µM should be directly compared since the other conditions change the initial density of the system to values which differ significantly from our simulations.
To further test our numerical model, we investigated how the contraction timescale varies with the number of motors in the system. As expected from [7], increasing D/M to 2 drastically changes the contraction timescale, with only minor changes in ∞ ( Figure 4A). Further increasing D/M does not significantly change either the timescale τ , or the final fraction ∞ . We also tested whether our model can reproduce the changes of contraction times for different length distributions of MTs, corresponding to different taxol concentrations, as seen in experiment. Indeed, as the mean length of MTs is increased, the contractions speed up ( Figure 4B, Table 2), yet the final fraction contracted doesn't change significantly. This is again consistent with experimental results, (Figure 1). We conclude that our simulations are consistent with key features of the experiments. Table 2. The fitted characteristic contraction time τ for different lengths of MTs at D/M = 2. All simulations start from the same network density as the one shown in  Having gained confidence in our simulation tools, we next inspect important aspects of the physics of the MT network which are experimentally difficult to access. We first ask how the system stress changes during the contraction process. We extract the stress from our simulations using the expressions detailed in Appendix C. In particular, we can distinguish between the stress tensor contributions from steric collisions, Σ col , and the active stress tensor generated by motors, Σ motor . Figure 5 shows Σ col , Σ motor , and (t), for L mean = 2 µm and different numbers of D/M , as a function of time. From this data, we make four key observations. First, we observed that the average stress tensors Σ col and Σ motor are dominated by their diagonal components, and the three diagonal components are within ∼ 20% of each other. We report in Figure 5 the average of trace of the stress tensor: Σ col = 1 3 Σ col xx + Σ col yy + Σ col zz , and Σ motor is calculated similarly. Second, we find that at low D/M , i.e. 1 or 2, the growth of (t) and of both stress contributions occur on the same timescale. This breaks down at larger D/M for which the stresses grow and saturate much faster than (t). Third, increasing D/M leads to increasing Σ motor , while Σ col rapidly saturates. Fourth, as Σ col saturates, the overall contraction timescale ceases to decrease (see also Figure 4). These observations point to the local structure formation changing significantly for larger D/M , which we further explore. Figure 6 shows the spatial distribution of MT volume fraction, motor stress, and collision stress for the D/M = 2 and 4 cases. At 60 s with D/M = 4, the MTs are quickly dragged by dyneins into dense local clusters while the clusters have not contracted much toward the center. At this stage, Σ col and Σ motor have significantly increased but (t) remains low. In contrast, for D/M = 2 the distibution at t = 60s is much less clustered. Consistently for D/M ≤ 2 the growth of Σ col and Σ motor is more synchronized with (t). This observation suggests a possible continuum model to predict the behavior of the entire network based on an 'equation-of-state' relating the local microscopic parameters and the binding-walking-unbinding cycles of dyneins to the local mechanical stress tensors Σ col and Σ motor , as the driving 'force' of the contraction, which will be the subject of future work.

Discussion
Here we examined how the dynamics of bulk MT network contraction in Xenopus extracts vary with taxol concentration. We find that using different taxol concentrations, the networks contract to approximately the same final density, even though the initial density of the networks varies by a factor of ≈ 3.6. While the final density of these networks is constant, we find that the timescale of the contraction process increases with increasing taxol concentration. In the high taxol regime, where the initial network density varies by only ≈ 14%, we find that the timescale increases by a factor of ≈ 2.5, arguing that the timescale is not simply varying proportionally to the initial network density. The length distributions of the MTs were measured for each taxol condition, and it was found that the average MT length decreases with increasing taxol concentration by a factor of ≈ 1.7 across the tested conditions. We postulate that it is this change in MT length that drives the changes in contraction time, and turn to simulation to test this idea.  We built a simulation tool to reveal more microscopic information of the network contraction process, and found that by placing dyneins randomly on each MT, we could simulate a contraction process consistent with experiments. We found that the characteristic contraction timescale, τ , is decreased by increasing only the MT length, which is consistent with the experimental result that shorter MTs at higher taxol concentration contract slower. We also investigated the effect of dynein number on the contraction process. Increasing the number of dynein per MT, D/M , from 1 to 2 sped up the contraction significantly, but further increasing the dynein number did not affect the contraction. Instead, an increase of D/M generated a local clustering structure at the early stage of the contraction, and then the local clusters contract. This process is clearly revealed by the spatial-temporal variations of collision and motor stresses, and their correlation with the contraction process (t). This indicates that we could possibly build a microscopic 'equation-of-state' to improve our previous coarsegrained model [7] to directly relate the stress terms in the model with microscopic motor behaviors, which will be the focus of future work.

Acknowledgments
This work was supported by the Kavli Institute for Bionano Science and Technology at Harvard University and National Science Foundation grants PHY-1305254, PHY-0847188, and DMR-0820484 (DJN), and partially supported by National Science Foundation Grants DMR-0820341 (NYU MRSEC), DMS-1463962, and DMS-1620331 (MJS). We thank Bryan Kaye for useful discussions.

Appendix A. Algorithm for MT network simulation
The simulation program tracks the motion of MTs driven by motor proteins by modeling MTs as rigid spherocylinders and motor proteins as Hookean springs with crosslinking ends. At this micron-sized scale the motion of objects in fluid is overdamped, and MTs are tracked with the following simple explicit Euler time stepping, where X i is the center of mass and q i is the unit orientation vector of MT i: The forces F motor and F col stem from motor proteins and collisions between MTs, respectively, with T motor and T col being the corresponding torques. Finally ∆X B and ∆q B are the Brownian contributions to the motion of MTs. The mobility matrix M i describes the effect of hydrodynamic drag. In principal all MTs are fully coupled through the Stokes equation, but in this work we ignore this coupling. Here, the M i of each MT i is approximated as if each MT is moving individually by itself in unbounded fluid. M i takes a 3 × 3 symmetric matrix form because we modeled each MT as an axisymmetric spherocylinder: where I is the 3 × 3 identity matrix. γ i and γ ⊥ i are the translational drag coefficients for motions parallel and perpendicular to the axis of the MT i, depending on the fluid viscosity η, the MT length L i , and the diameter D i . Drags are solvable from the slenderbody theory [33]: where b i = − (1 + 2 ln(D i /L i )) is a geometric parameter. For MTs we have L i D i , and b i ≈ 2 ln(L i /D i ), and so we have γ ⊥ i ≈ 2γ i . In simulations, we used the same diameter D for all MTs.
At each timestep the Brownian displacement and rotation ∆X B i and ∆q B i are generated as Gaussian random vectors. The means are zero, while the variances are given by the local fluctuation-dissipation relation: We follow the method of Tao et al. [34] to generate these two random vectors. For Eq. A.1, the motor force and torque are described in the main text. The collision force F col i and torque T col i are calculated in the simulation as detailed in Appendix B.

Appendix B. Collision between microtubules
When two MTs i and j are close to each other, the shortest distance d ij between their axes is calculated. If d ij is smaller than the microtubule diameter D then a collision force F col is calculated with a WCA-type repulsive force to prevent them from overlapping and crossing each other. The original WCA potential is stiff and poses a severe limit on the maximum timestep ∆t in simulations. We used a softened WCA potential to enable larger timesteps: Where F col (d) is a shifted WCA potential with a shift parameter α. For α = 0, F col (d) becomes the usual WCA repulsive force. F col L (d) is a linearized continuation of F col (d) for d < βD, to prevent the blow-up of collision force at small d. In simulations we fix α = 0.1 and β = 0.95, to preserve the repulsion between MTs, but enable the simulations to complete within reasonable simulation time. Also due to the use of a soft repulsive force instead of an accurate rigid interaction, the effective diameter of MTs are estimated to be about ∼ 80% of the specified diameter D [34]. Therefore in simulations D is set to 30 nm to reproduce the true diameter 24 nm of MTs.

Appendix C. The calculation of stress and its spatial distribution
The stress can be calculated by the binding force and collision force in Eq. (A.1). Since both contributions are pairwise between MTs, and satisfy Newton's third law, we can calculate the stress tensor pair-by-pair with the virial theorem. Also, because the MTs are thin-and-long slender cylinders, the calculation can be further simplified [35].
For two MTs i and j, their contribution to the total system stress is a 3 × 3 tensor formed by the outer product: where r (ij) = x (i) − x (j) , is the vector pointing from the forcing location x (j) on MT j to x (i) on MT i. f (i,j) is the force from j to i. For the collision force, x (i) and x (j) are the two collision points on each MT. For the dynein binding force, x (i) and x (j) are the two binding locations (two ends) of one dynein on the two MTs. If more than one dyneins bind i and j, the contribution from each dynein is calculated and added independtly. Each σ (i,j) contributes to the total system stress tensor at the spatial location x (i) + x (j) /2, which is the center of the two forcing location. For a certain volume V in space, the average stress tensor in this volume is a simple arithmetic average of all σ ij in this volume: In the simulations, to detect the binding and collision events a nearest neighbor search procedure is performed during each timestep, and long MTs significantly slow down its efficiency. Also, due to the periodic boundary condition used in simulations, a MT with length of the simulation box size may interact with its own periodic image, and generates some unphysical results. Therefore in simulations the lognormal distribution is clipped at L max , and the mean length of the MTs changes: Erfc µ−log(Lmax) √ 2σ .

(D.2)
While the mode length does not change with L max . In simulations we used L max = 4 µm, which is half of the width of a simulation box with size 12 × 8 × 8µm. We also fixed σ = 0.5 to match the experiment data, and the actual simulation settings are detailed in Table D1.