Connecting macroscopic dynamics with microscopic properties in active microtubule network contraction
Peter J. Foster, Wen Yan, Sebastian Fürthauer, Michael J. Shelley, Daniel J. Needleman
CConnecting macroscopic dynamics with microscopicproperties in active microtubule networkcontraction
Peter J. Foster , , Wen Yan , , Sebastian F¨urthauer , MichaelJ. Shelley , , Daniel J. Needleman , John A. Paulson School of Engineering and Applied Sciences, FAS Center forSystems Biology, Harvard University, Cambridge MA, United States Center for Computational Biology, Flatiron Institute, Simons Foundation, NewYork, NY, United States Courant Institute of Mathematical Science, New York University, New York, NY,United States Department of Molecular and Cellular Biology, Harvard University, Cambridge MA,United States Equal contributionE-mail: [email protected], [email protected]
June 2017
Abstract.
The cellular cytoskeleton is an active material, driven out of equilibriumby molecular motor proteins. It is not understood how the collective behaviorsof cytoskeletal networks emerge from the properties of the network’s constituentmotor proteins and filaments. Here we present experimental results on networksof stabilized microtubules in
Xenopus oocyte extracts, which undergo spontaneousbulk contraction driven by the motor protein dynein, and investigate the effects ofvarying the initial microtubule density and length distribution. We find that networkscontract to a similar final density, irrespective of the length of microtubules or theirinitial density, but that the contraction timescale varies with the average microtubulelength. To gain insight into why this microscopic property influences the macroscopicnetwork contraction time, we developed simulations where microtubules and motorsare explicitly represented. The simulations qualitatively recapitulate the variation ofcontraction timescale with microtubule length, and allowed stress contributions fromdifferent sources to be estimated and decoupled.
Keywords : active matter, microtubule, motor protein, dynein,
Xenopus extract
1. Introduction
Active matter is a class of materials held out of equilibrium by the local conversion ofenergy from a reservoir into mechanical work at the scale of the system’s components a r X i v : . [ q - b i o . S C ] J un ontraction in active microtubule networks α pp ,where p is an MT orientation vector, and α < p direction[13]. Different dynamics can follow from different MT interactions, and clustering ofMT ends is a generic mechanism leading to contraction [14]. This has been previouslyexplored in other purified systems containing kinesin motor proteins [15, 16]. In thiswork, as in our earlier work of stabilized MTs in meiotically-arrested Xenopus oocyteextracts [7] and in several related studies [17, 18, 19], we consider a spontaneously formedMT network whose dynamics is driven by minus-end bundling by dynein motors. Thisleads to an apparently isotropic contraction, and we proposed a continuum model [7],which quantitatively reproduced the experimental results. Crucially, the model helpedexplain why networks contract to a preferred final density regardless of sample geometryand motor concentration. The continuum model assumed that the active material stresstensors were isotropic and of the form f ( ρ ) I , with f < ontraction in active microtubule networks
2. Methods
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 beingflash frozen in liquid nitrogen and stored at -80 ◦ C until use.
Microfluidic device templates were designed using AutoCAD 360 and Silhouette Studiosoftware. Device templates were cut from 125 µ m tape (3M Scotchcal) using a SilhouetteCameo 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, degassedunder vacuum, and cured overnight at 60 ◦ C. The cured devices were then removed fromthe masters and inlets and outlets were created using biopsy punches. Devices andcoverslips were treated with air plasma using a corona device, bonded, and loaded withpassivation mix composed of 5 mg/mL BSA and 2.5% (w/v) Pluronic F-127 beforeovernight incubation at 12 ◦ C. 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 ≈ µ M. Then, 0.5 µ L of taxol suspended in DMSO was then added to the extract tothe indicated final concentration. Extracts were then loaded into passivated microfluidicdevices, sealed with vacuum grease, and imaged using spinning disk confocal microscopy(Nikon TE2000-E microscope, Yokugawa CSU-X1 spinning disk, Hamamatsu ImagEMcamera, 2x objective, Metamorph acquisition software). The time t = 0 correspondsto when imaging begins, typically < (cid:15) ( t ) curves were fit using time pointswhen (cid:15) ( t ) > . In order to estimate network densities, we first assume that Alexa-647 labeled tubulinuniformly incorporates into MTs, and the overall concentration of tubulin is taken tobe constant and equal to the previously measured value of 18 µ M [21]. From this we ontraction in active microtubule networks µM = (cid:80) Ni =1 M i (cid:80) Ni =1 V i = (cid:80) Ni =1 βI i ( (cid:96)A ) (cid:80) Ni =1 (cid:96)A = β (cid:104) I (cid:105) N , where M i is the tubulin mass in pixel i, β is a constant conversion factor betweenconcentration in micromolar and fluorescence intensity, (cid:96) is the depth of the sample, A is the area of the pixel, V i = (cid:96)A is the volume corresponding to pixel i, N is the totalnumber of pixels, I i is the measured intensity of pixel i, A is the area of pixel i, and (cid:104) I (cid:105) N is the measured intensity averaged over all pixels in the channel. From this we caninfer, β = 18 µM (cid:104) I (cid:105) N . The intensity at each pixel in the network contains a contribution N i from polymerizedtubulin and a contribution B i from monomeric, unpolymerized tubulin. The signalfrom monomeric tubulin is assumed to be constant and homogeneous throughout thechannel. Thus, at the given time point where the background and network intensitiesare measured, the average concentration of polymerized tubulin in the network is givenby, (cid:104) ρ (cid:105) ( t ) = βN network N network (cid:88) i =1 N i = βN network N network (cid:88) i =1 [ I i − B ] , = β [ (cid:104) I (cid:105) network − B ] = 18 µM (cid:104) I (cid:105) network − B (cid:104) I (cid:105) N , where N network is the number of pixels in the network at the time point, and (cid:104) I (cid:105) 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 + τ wasselected, where τ is the characteristic contraction time and T c is the lag time betweenthe beginning of imaging and the beginning of network contraction [7]. This frame wascorrected for inhomogeneous illumination and the camera’s dark current, and analyzedas above to find (cid:104) ρ (cid:105) ( T c + τ ). As we assume that no MTs are created, destroyed, or addedto the network during the contraction process, the total mass of MTs in the networkmust be conserved. Thus, M network = (cid:104) ρ (cid:105) V = (cid:104) ρ (cid:105) ( t = T c + τ ) V T c + τ = (cid:104) ρ (cid:105) F V F , where V , V T c + τ , V F are the volume of the network at the beginning, at time T c + τ , andat the final state. Then, (cid:104) ρ (cid:105) (0) = (cid:104) ρ (cid:105) ( T c + τ ) V T c + τ V = (cid:104) ρ (cid:105) ( T c + τ ) W T c + τ H T c + τ L T c + τ W H L , where W, H, L are the width, height and length of the network at different times, denotedby the subscripts as with the volume. We assume that, as the network is pinned at thechannel’s inlet and outlet, there is no change in volume along the channel’s length, and ontraction in active microtubule networks L = L T c + τ . Combining Eqns. 1 and 2 in Results, (cid:15) ( t = T c + τ ) = (cid:15) ∞ (1 − e − ) = W − W ( t = T c + τ ) W = 1 − W T c + τ W , which simplifies to, W T c + τ W = 1 − (cid:15) ∞ (1 − e − ) . We further assume that the change in the network’s height follows the same functionalform as the change in width and thus, (cid:104) ρ (cid:105) (0) = (cid:104) ρ (cid:105) ( T c + τ )[1 − (cid:15) ∞ (1 − e − )] . For the final density of the network, similar reasoning can be used to show, (cid:104) ρ (cid:105) F = 1(1 − (cid:15) ∞ ) (cid:104) ρ (cid:105) = [1 − (cid:15) ∞ (1 − e − )] (1 − (cid:15) ∞ ) (cid:104) ρ (cid:105) ( T c + τ ) . Stabilized MTs were dissociated from motor proteins and fixed as previously described[22]. MTs were allowed to assemble at room temperature for 30 minutes (unlessotherwise 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 , 1 mM EGTA, and 20 µ Mtaxol). 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 dilutedsample was spread between a slide and a 22 ×
22 mm coverslip. After waiting 30 minutesto allow the MTs to adhere to the coverslip, MTs were imaged using spinning diskconfocal microscopy (Nikon TE2000-E microscope, Yokugawa CSU-X1 spinning disk,Hamamatsu ImagEM camera, 60x objective, Metamorph acquisition software). Activecontours were fit to individual MTs using the ImageJ plugin JFilament [23], and MTlengths were determined from the contours using custom MATLAB software. For eachtaxol concentration, distributions of MT lengths were fit to a log-normal distributionto find the location parameter, µ , and the scale parameter, σ . In each case, the modemicrotubule length for the log-normal distribution is given by l mode = e µ − σ .
3. Results
To further investigate the dynamics of contracting MT networks in
Xenopus oocyteextracts, 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 inletand outlet using vacuum grease to prevent evaporation, and imaged at low magnification ontraction in active microtubule networks
125 µm 0.9 mm m m A B
Time (min) ε ( t ) , F r a c t i on C on t r a c t ed C C ha r a c t e r i s t i c T i m e , τ , ( m i n ) ε , F i na l F r a c t i on C on t r a c t ed ∞ D Taxol Concentraction (µM) E Taxol Concentraction (µM) F i na l N e t w o r k D en s i t y , ρ ( µ M ) I n i t i a l N e t w o r k D en s i t y , ρ ( µ M ) I Figure 1. Microtubule networks undergo spontaneous bulk contraction (A)Schematic showing characteristic dimensions of the microfluidic chamber. (B) Timecourse of MT network contraction when 10 µ M taxol is added (Scale bar: 500 µ m) (C)Curves showing the fraction contraction, (cid:15) (t) as a function of time for varying taxolconcentration. Curves are mean ± s.e.m. (D) Characteristic time, τ , and final fractioncontraction, (cid:15) ∞ , as a function of taxol concentration. (E) Initial network density, ρ I ,(blue line) and final network density, ρ , (red line) as a function of taxol concentration.The red dashed line denotes the average value of ρ =49.6 µ M . (Methods). Within minutes of taxol addition, the MT networks were found to undergoa spontaneous bulk contraction as shown in our previous work [7] (Figure 1B, Movie1). The size of the network along the width of the channel, W ( t ), was measured byidentifying the pixels with high intensity (belonging to the network) and recorded as afunction of time. Then the fraction contracted (cid:15) ( t ) was calculated, defined as (cid:15) ( t ) = W − W ( t ) W , (1)where W is the width of the channel, typically 0.9 mm. This process was repeated forvarying final concentrations of taxol, and the (cid:15) ( t ) curves for each taxol condition wereaveraged together to produce master curves for each condition (Figure 1C). The (cid:15) ( t )curves were well fit by a saturating exponential function, (cid:15) ( t ) = (cid:15) ∞ [1 − e − ( t − Tcτ ) ] , (2)where (cid:15) ∞ 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 contractionprocess. Equation 2 was fit to the (cid:15) ( t ) curves for each experiment to extract values forthe characteristic time, τ , and the final fraction contracted, (cid:15) ∞ . Values of τ and (cid:15) ∞ were averaged for each taxol condition (Figure 1D). The characteristic timescale, τ , was ontraction in active microtubule networks Microtubule Length (µm) C P r obab ili t y Taxol Concentraction (µM) B M ode M i c r o t ubu l e Leng t h ( µ m ) A Microtubule Length (µm) P r obab ili t y P D F Figure 2. Microtubule length distributions vary with taxol concentration (A) Histograms of MT lengths measured for differing taxol concentrations. Inset:Example distrubution for 2.5 µ M taxol with log-normal fit (red line). (B) The modeMT length decreases with increasing taxol concentration. (C) Example histograms for2.5 µ M taxol, where MT networks were dissociated and fixed at either 5 minutes or 45minutes after taxol addition. found to increase with approximate linearity for taxol concentrations > . µ M, whilethe final fraction contracted, (cid:15) ∞ , decreased slightly with increasing taxol concentration.We next investigated how changing taxol concentration influenced both the initialdensity of the MT network, ρ I , and the final density of the MT network, ρ . In systemsof purified tubulin, the concentration of tubulin polymerized into MTs has been shownto increase with increasing taxol concentration [24]. Fluorescence intensity was used as aproxy for tubulin concentration, and was calibrated using previously measured values forthe total tubulin concentration in Xenopus extracts (Methods). While the initial densityof the MT network, ρ I , monotonically increased with increasing taxol concentration, thefinal density of the MT network displayed no obvious trend (Figure 1E), and values ofthe final network density vary from the mean final density of ρ = 49 . µ M by less than25%. 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 initialMT density, ρ I , with increasing taxol concentration appears to saturate for large taxolconcentrations. 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 valueby only 14% (std/mean), far less than the ≈ × increase seen in the characteristiccontraction timescale, τ , over this same range. Thus, τ does not scale proportionally tothe initial density of MTs, ρ I , and the change in contraction timescale must come fromsome other effect of changing taxol concentration.Previously, it has been shown in clarified Xenopus extracts that the size ofmotor-organized taxol-stabilized MT assemblies, termed “pineapples”, decreased withincreasing taxol concentration, presumably due to a decrease in the average length ofMTs [22]. Thus, we next sought to measure the length distribution of MTs in our systemfor varying taxol concentrations. We used a previously described method to dissociateand fix MTs [22] (Methods). Briefly, MTs were allowed to assemble for 5 minutes after ontraction in active microtubule networks coverslip. MTs were allowed to adhere for 30 minutes before imaging. Anactive contour was fit to MTs in the images [23], and the length of the active contourwas used as a measure for the MT length. This process was repeated for each taxolconcentration.Example histograms of MT lengths for each taxol concentration are shown in Figure2A. Visually, the peak of the distribution shifts towards smaller values of MT lengthfor increasing taxol concentration. Empirically, we find that the distributions of MTlengths are well fit by a log-normal distribution (Figure 2A, Inset). As fitting log-normaldistributions to the measured histograms allows a more robust estimate of the modeMT length than using a purely empirical estimate, we fit the MT length distributionsfor each taxol concentration in order to find µ and σ , the two parameters of the log-normal distribution, and used these parameters to estimate the mode MT length foreach condition. The mode MT length was found to decrease with increasing taxolconcentration, varying by a factor of ≈
4. Model and Simulation
Our experimental results suggest that the change in contraction timescale observed forchanging taxol concentrations may be due to changes in the MT length distribution.While our earlier continuum active gel theory for MT network contractions [7] capturesthe global contractile behavior of the network accurately, understanding the dependenceof its parameters on microscale properties, such as MT length, is beyond the scope of themodel. To resolve microscale changes and study their effects on the network’s emergentproperties, and to test the dependence of contraction timescale on MT lengths, we turnto simulation. Our simulation tracks the behavior of a suspension of fixed length MTsactuated by dynein motors which actively crosslink them and drag them through thefluid. ontraction in active microtubule networks In our modelling, we represent MTs as rigid spherocylinders that interact sterically andthrough motor protein coupling. In principle the fluid in which MTs are suspendedcouples their dynamics globally. Since solving the full Stokes problem numerically isprohibitively expensive we here neglect hydrodynamic many-body couplings. Insteadwe approximate the effects of the fluid by a local drag, that is, each MT acts as if itmoved 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 arguethat they can be safely ignored for our purposes.Suspensions of passive Brownian spherocylinders are well understood in both theoryand simulation, and our numerical framework (in Appendix A) is based on well-knownwork [25, 26, 27]. What sets our simulations apart from the previous work is the presenceof 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 amoving crosslinking end with stochastic binding and unbinding behaviors. The twoends are assumed to be connected by a Hookean spring with a rest length l D (cid:39)
40 nm,which is the size of a force-free dynein. Collisions amongst dyneins or between dyneinsand MTs are ignored for simplicity. In our model, dyneins have one fixed end, whichstays rigidly attached to one MT throughout the simulation, and one motor end whichstochastically 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 anyforce to it. Given the motor’s small size, the orientation of such a motor is dominatedby Brownian motion, i.e. is random. If a free motor head is closer than a characteristiccapture 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 closeenough 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 locationwe choose the perpendicular projection the dynein center onto the target MT. Boundmotor ends walk towards the minus end of the MT they bind to, and exert a springforce and torque ( F motor = κ D ( l D − l D ) e D , T motor = r D × F motor ), where l D , κ D are thedynein’s current length and its spring constant, respectively, e D is the direction of theHookean spring force along the direction of dynein, and r D is a vector pointing fromthe center of mass of a MT to the point where the dynein binds and applies a Hookeanspring force. Consistent with earlier work we impose a force velocity relation relation[13] and set the velocity v = v M (cid:0) − min( F motor /F stall , (cid:1) iff F motor < F stall (3)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 τ u = τ u exp (cid:0) − F motor / ( F stall − F motor ) (cid:1) , (4) ontraction in active microtubule networks (cid:39) exp (cid:0) − F motor /F stall (cid:1) for F motor (cid:28) F stall , see [29], and diverges for F motor → F stall . Here τ u is the force-free unbinding time. Finally, dyneins that reachthe minus end do not immediately detach, but remain at the minus end, i.e. v = 0,keeping the same unbinding frequency. Table 1.
The parameters fixed in simulations.
Parameter Explanation Value Reference l D dynein free length 40 nm Estimation F stall dynein stall force 1 pN Estimation[29] κ D dynein spring constant 1 pN µ m − Estimation[30] r cap dynein capture radius 80 nm = 2 l D Assumption v M dynein max walking velocity 1 µ m / s Estimation[31] τ b dynein binding timescale 10 − s Assumed to be short τ u dynein unbinding timescale 10 s Estimation[17, 31] η fluid viscosity 0 .
02 pN µ m − s − Estimation[32] k B T energy scale at 300 K 4 . × − pN µ m 300 KFor this work we explored how key aspects of the microscopic model change theemergent behavior of the network. In particular, we varied the number of dyneins affixedto each MT, D/M , and the length distribution of MTs.The full model has a large number of parameters. However, many can be fixedfrom estimates in the literature. We summarize in Table 1 the parameters that are keptunchanged throughout the simulations presented here. We also fix the number of MTsin our simulations to be 10 µ M, which is the approximate initial concentration found inexperiments with 5 - 25 µ M of taxol. Following the experimental results, the MT lengthdistribution is taken to be lognormal with parameters ( µ, σ ), with σ fixed at 0 .
5. Notethat in simulation we use shorter MTs than in the experiments due to limitations incomputational resources, see Appendix D.In the simulations presented here, each MT carries the same number,
D/M , of fixeddynein molecules, attached at random locations on the MT that carries them. Note thatthis is different from the assumptions made in earlier work [7], which has dyneins all fixedto the minus ends of MTs. This difference is motivated by the numerical observationthat networks with all dyneins at minus ends will precipitate into asters, while networkswith dyneins randomly affixed to MTs are more percolated and globally contract. It’sunclear whether this difference is due to an artifact in the current simulations. Furtherexploring 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 × × µ m,with periodic boundary conditions along the long direction, mimicking the experimentalgeometry. The initial state is a random arrangement of MTs, in which all dynein motor ontraction in active microtubule networks Figure 3. Motor proteins contract the MT network.
The snapshots are takenat different time, for the simulation case of L mean = 2 µ m. In this simulation D/M = 2which means two dyneins are randomply placed on each MT at the beginning of thesimulation. The MTs are colored in blue, the minus ends are marked by green, and thedyneins are colored in red. The simulation box is set to be periodic in the x direction(the horizontal direction in the figure), and unbounded in the y and z direction. Thesimulation box size = 12 × × µ m, and 1365 MTs and 2730 dyneins are tracked. ends are unbound. Figure 3 shows a representative example of the numerically observed contractions. Afteran initial phase of (cid:39)
30 s, in which MTs locally rearrange and cluster, bulk contractiondriven by dyneins begins. It takes ∼ s to create a single, fully percolated network.The system has approached a steady state at (cid:39)
500 s, where the network has contractedto a cylindrical ribbon along the periodic direction of the simulation box.To compare the simulations to experiments, we calculate the fraction contracted (cid:15) ( t ) = R − R ( t ) R , (5)from our simulations. Here, R ( t ) is the radius of a cylinder aligned with the periodicdirection, and centered around the MT center of mass of the system portion whichcontains 50% of the total MT mass. By fitting (cid:15) ( t ) to the saturation exponential,e.g. Eqn 2, we now obtain a contraction timescale τ , and a final fraction contracted (cid:15) ∞ ,which can be directly compared to experiments. As was done for fitting the experimental ontraction in active microtubule networks (cid:15) ( t ) < . τ depends on the initial system size. Sinceour numerical system is far smaller than the experimental one, we do not expect τ to quantitatively match experiments, but only to reproduce trends with varyingparameters. Using the parameters from [7] for the continuum model, which was donewith 2.5 µ M taxol, and a system size of 8 µ m, we expect a contraction time of τ ≈ (cid:15) ∞ , is easier to directly compare between experimentsand simulation, since it does not depend on system size. Encouragingly, for the caseof D/M = 1, the simulated contraction produces a similar final fraction contracted of (cid:15) ∞ ≈ .
55, as in experiment, (Figure 1D and Figure 4). Note that here only taxolconcentrations of 5, 10, 25 µ M should be directly compared since the other conditionschange the initial density of the system to values which differ significantly from oursimulations.To further test our numerical model, we investigated how the contraction timescalevaries 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 (cid:15) ∞ (Figure 4A). Further increasing D/M does not significantly change either the timescale τ , or the final fraction (cid:15) ∞ . We also tested whether our model can reproduce thechanges of contraction times for different length distributions of MTs, correspondingto different taxol concentrations, as seen in experiment. Indeed, as the mean length ofMTs is increased, the contractions speed up (Figure 4B, Table 2), yet the final fractioncontracted doesn’t change significantly. This is again consistent with experimentalresults, (Figure 1). We conclude that our simulations are consistent with key featuresof 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 inFigure 3. L mean ( µ m) 1.618 2.007 2.297 2.511 τ (s) 68.6 39.4 33.4 28.8Having gained confidence in our simulation tools, we next inspect important aspectsof the physics of the MT network which are experimentally difficult to access. We firstask how the system stress changes during the contraction process. We extract the stressfrom our simulations using the expressions detailed in Appendix C. In particular, wecan distinguish between the stress tensor contributions from steric collisions, Σ col , andthe active stress tensor generated by motors, Σ motor .Figure 5 shows Σ col , Σ motor , and (cid:15) ( 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, weobserved that the average stress tensors Σ col and Σ motor are dominated by their diagonal ontraction in active microtubule networks Figure 4. The contraction process varies with the dynein number
D/M and MT length L mean . Each (cid:15) ( t ) data curve is shown with its saturating exponentialfit (2). The numbers of MTs are set at the same value as the 25 µ M taxol case inexperiments, for all simulations discussed in this paper. More detailed setting aboutthe parameters and MT length distribution can be found in appendices. components, and the three diagonal components are within ∼
20% of each other. Wereport in Figure 5 the average of trace of the stress tensor: Σ col = (cid:0) Σ colxx + Σ colyy + Σ colzz (cid:1) ,and Σ motor is calculated similarly. Second, we find that at low D/M , i.e. 1 or 2, thegrowth of (cid:15) ( t ) and of both stress contributions occur on the same timescale. This breaksdown at larger D/M for which the stresses grow and saturate much faster than (cid:15) ( 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, andcollision stress for the
D/M = 2 and 4 cases. At 60 s with
D/M = 4, the MTsare quickly dragged by dyneins into dense local clusters while the clusters have notcontracted much toward the center. At this stage, Σ col and Σ motor have significantlyincreased but (cid:15) ( t ) remains low. In contrast, for D/M = 2 the distibution at t = 60sis much less clustered. Consistently for D/M ≤ Σ col and Σ motor ismore synchronized with (cid:15) ( t ). This observation suggests a possible continuum model topredict the behavior of the entire network based on an ‘equation-of-state’ relating thelocal microscopic parameters and the binding-walking-unbinding cycles of dyneins to thelocal mechanical stress tensors Σ col and Σ motor , as the driving ‘force’ of the contraction,which will be the subject of future work. ontraction in active microtubule networks
5. 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 initialdensity of the networks varies by a factor of ≈ ≈ ≈ ≈ Figure 5. The correlation between (cid:15) ( t ) and stress varies with D/M . Thefraction contracted is shown in red. The motor stress Σ motor and the collision stressΣ col are colored in dark and light blue, respectively. The stresses are calculated withthe method described in Appendix C, where the volume V is the volume of the entiresimulation cell of 12 × × µ m. Here the scalar stresses are calculated as the averageof diagonal components of the stress tensor: Σ col = (cid:0) Σ colxx + Σ colyy + Σ colzz (cid:1) , and Σ motor is calculated similarly. ontraction in active microtubule networks Figure 6. The spatial distributions of MT volume fraction and stress varywith the dynein number
D/M . Snapshots are taken at t = 60 s and t = 500 s forboth cases. On the left: L mean = 2 µ m and D/M = 2. On the right: L mean = 2 µ mand D/M = 4. The spatial distribution is calculated by spatial binning with a regularmesh, and the mesh grid size = 0 . µ m. Here the scalar stresses are calculated asmentioned in Figure 5. We built a simulation tool to reveal more microscopic information of the networkcontraction process, and found that by placing dyneins randomly on each MT, wecould simulate a contraction process consistent with experiments. We found that thecharacteristic contraction timescale, τ , is decreased by increasing only the MT length,which is consistent with the experimental result that shorter MTs at higher taxolconcentration contract slower. We also investigated the effect of dynein number onthe contraction process. Increasing the number of dynein per MT, D/M , from 1 to2 sped up the contraction significantly, but further increasing the dynein number didnot affect the contraction. Instead, an increase of
D/M generated a local clusteringstructure 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 motorstresses, and their correlation with the contraction process (cid:15) ( t ). This indicates that wecould possibly build a microscopic ‘equation-of-state’ to improve our previous coarse-grained model [7] to directly relate the stress terms in the model with microscopic motorbehaviors, which will be the focus of future work. ontraction in active microtubule networks
6. Acknowledgments
This work was supported by the Kavli Institute for Bionano Science and Technologyat Harvard University and National Science Foundation grants PHY-1305254, PHY-0847188, and DMR-0820484 (DJN), and partially supported by National ScienceFoundation 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 modelingMTs as rigid spherocylinders and motor proteins as Hookean springs with crosslinkingends. At this micron-sized scale the motion of objects in fluid is overdamped, and MTsare tracked with the following simple explicit Euler time stepping, where X i is the centerof mass and q i is the unit orientation vector of MT i :∆ X i = M i · (cid:0) F motori + F coli (cid:1) ∆ t + ∆ X Bi , (A.1)∆ q i = 1 γ Ri (cid:0) T motori + T coli (cid:1) × q i ∆ t + ∆ q Bi . (A.2)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 principalall MTs are fully coupled through the Stokes equation, but in this work we ignorethis coupling. Here, the M i of each MT i is approximated as if each MT is movingindividually by itself in unbounded fluid. M i takes a 3 × M − i = γ (cid:107) i q i q i + γ ⊥ i ( I − q i q i ) , (A.3)where I is the 3 × γ (cid:107) i and γ ⊥ i are the translational drag coefficientsfor motions parallel and perpendicular to the axis of the MT i , depending on the fluidviscosity η , the MT length L i , and the diameter D i . Drags are solvable from the slender-body theory[33]: γ ⊥ i = 8 πηL i b i + 2 , γ (cid:107) i = 8 πηL i b i , γ ri = 2 πηL i b i + 2) , (A.4)where b i = − (1 + 2 ln( D i /L i )) is a geometric parameter. For MTs we have L i (cid:29) D i ,and b i ≈ L i /D i ), and so we have γ ⊥ i ≈ γ (cid:107) i . In simulations, we used the samediameter D for all MTs.At each timestep the Brownian displacement and rotation ∆ X Bi and ∆ q Bi aregenerated as Gaussian random vectors. The means are zero, while the variances aregiven by the local fluctuation-dissipation relation: (cid:104) ∆ X Bi ∆ X Bi (cid:105) = 2 k B T M i ∆ t, (A.5) (cid:104) ∆ q Bi ∆ q Bi (cid:105) = 2 k B T ( I − q i q i ) ∆ t/γ ri . (A.6) ontraction in active microtubule networks F coli and torque T coli 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 theiraxes 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 overlappingand crossing each other. The original WCA potential is stiff and poses a severe limit onthe maximum timestep ∆ t in simulations. We used a softened WCA potential to enablelarger timesteps: d ≥ βD : F col ( d ) = − (cid:18) k B TD (cid:19) (cid:20) α + d/D ) − α + d/D ) (cid:21) , (B.1) d < βD : F colL ( d ) = F col ( βD ) + ∂F col ∂d | d = βD ( d − βD ) . (B.2)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 colL ( 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 wefix α = 0 . β = 0 .
95, to preserve the repulsion between MTs, but enable thesimulations to complete within reasonable simulation time. Also due to the use of a softrepulsive force instead of an accurate rigid interaction, the effective diameter of MTs areestimated 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). Sinceboth contributions are pairwise between MTs, and satisfy Newton’s third law, we cancalculate the stress tensor pair-by-pair with the virial theorem. Also, because the MTsare 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 × σ ( i,j ) = − (cid:0) x ( i ) − x ( j ) (cid:1) f ( i,j ) , (C.1)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 thetwo collision points on each MT. For the dynein binding force, x ( i ) and x ( j ) are the twobinding locations (two ends) of one dynein on the two MTs. If more than one dyneinsbind 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 (cid:0) x ( i ) + x ( j ) (cid:1) /
2, which is the center of the two forcing location. For a certain volume V ontraction in active microtubule networks σ ij in this volume: Σ colV = − V (cid:88) ∈ V σ ( i,j ) ,col , Σ motorV = − V (cid:88) ∈ V σ ( i,j ) ,motor . (C.2) Appendix D. The length distribution of microtubules
The length distribution of MTs is taken to be given by a lognormal distribution withparameters µ and σ : P ( L ) = e − (log( L ) − µ )22 σ √ πLσ , L mean = e µ + σ / , L mode = e µ − σ . (D.1)In the simulations, to detect the binding and collision events a nearest neighborsearch procedure is performed during each timestep, and long MTs significantly slowdown its efficiency. Also, due to the periodic boundary condition used in simulations, aMT with length of the simulation box size may interact with its own periodic image, andgenerates some unphysical results. Therefore in simulations the lognormal distributionis clipped at L max , and the mean length of the MTs changes: L mean = e µ + σ Erfc (cid:16) µ − log( L max )+ σ √ σ (cid:17) Erfc (cid:16) µ − log( L max ) √ σ (cid:17) . (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 boxwith size 12 × × µ m. We also fixed σ = 0 . Table D1.
The settings of MT length in simulations. e µ ( µ m) 1.0 1.5 2.0 2.5 3.0 L mean ( µ m) 1.123 1.618 2.007 2.297 2.511 L mode ( µ m) 0.779 1.168 1.558 1.947 2.336 Appendix E. References [1] Marchetti MC, Joanny JF, Ramaswamy S, Liverpool TB, Prost J, Rao M, et al. Hydrodynamicsof soft active matter.
Rev Mod Phys . 2013;85: 1143–1189. doi:10.1103/RevModPhys.85.1143[2] Prost J, J¨ulicher F, Joanny JF. Active gel physics.
Nat Phys . 2015;11: 111–117.doi:10.1038/nphys3224[3] Sanchez T, Chen DTN, DeCamp SJ, Heymann M, Dogic Z. Spontaneous motion in hierarchicallyassembled active matter.
Nature . 2012;491: 431–434. doi:10.1038/nature11591[4] Linsmeier I, Banerjee S, Oakes PW, Jung W, Kim T, Murrell MP. Disordered actomyosin networksare sufficient to produce cooperative and telescopic contractility.
Nature Communications .2016;7: 1–9. doi:10.1038/ncomms12615 ontraction in active microtubule networks [5] Surrey T, N´ed´elec F, Leibler S, Karsenti E. Physical properties determining self-organization ofmotors and microtubules. Science. 2001;292: 1167–1171. doi:10.1126/science.1059758[6] N´ed´elec FJ, Surrey T, Maggs AC, Leibler S. Self-organization of microtubules and motors. Nature.1997;389: 305–308. doi:10.1038/38532[7] Foster PJ, F¨urthauer S, Shelley MJ, Needleman DJ. Active contraction of microtubule networks. eLife . 2015;4: e10837. doi:10.7554/eLife.10837[8] Brugu´es J, Needleman D. Physical basis of spindle self-organization. P Natl Acad Sci Usa . 2014;111:18496–18500. doi:10.1073/pnas.1409404111[9] Naganathan SR, F¨urthauer S, Nishikawa M, J¨ulicher F, Grill SW. Active torque generationby the actomyosin cell cortex drives left–right symmetry breaking. eLife . 2014;3: e04165.doi:10.7554/eLife.04165[10] Mayer M, Depken M, Bois JS, J¨ulicher F, Grill SW. Anisotropies in cortical tension reveal thephysical basis of polarizing cortical flows.
Nature . 2010;467: 617–621. doi:10.1038/nature09376[11] Shelley M J 2016
Annual Review of Fluid Mechanics Nature . 2012;491: 431-434.[13] Gao T, Blackwell R, Glaser M A, Betterton M D and Shelley M J 2015
Physical Review E ISSN 1539-3755, 1550-2376[14] Belmonte J, Leptin M, N´ed´elec F. A Theory That Predicts Behaviors Of Disordered CytoskeletalNetworks bioRxiv . 2017. doi:10.1101/138537[15] N´ed´elec F, Surrey T. Dynamics of microtubule aster formation by motor complexes. ComptesRendus de l’Acadmie des Sciences - Series IV - Physics-Astrophysics. 2001;2: 841–847.doi:10.1016/s1296-2147(01)01227-6[16] Hentrich C, Surrey T. Microtubule organization by the antagonistic mitotic motors kinesin-5 andkinesin-14. The Journal of Cell Biology. 2010;189: 465–480. doi:10.1083/jcb.200910125[17] Tan R, Foster PJ, Needleman DJ, McKenney RJ. Cooperative Accumulation Of Dynein-Dynactin At Microtubule Minus-Ends Drives Microtubule Network Reorganization. bioRxiv .2017. doi:10.1101/140392[18] White D, Vries G, Martin J, Dawes A. Microtubule patterning in the presence of moving motorproteins.
Journal of Theoretical Biology . 2015;382: 81-90.[19] Taunenbaum ME, Vale RD, McKenney R. Cytoplasmic dynein crosslinks and slides anti-parallelmicrotubules using its two motor domains. 2013;2: e00943.[20] Hannak E, Heald R. Investigating mitotic spindle assembly and function in vitro using Xenopuslaevis egg extracts. Nat Protoc. 2006;1: 2305–2314. doi:10.1038/nprot.2006.396[21] Parsons SF, Salmon ED. Microtubule assembly in clarified Xenopus egg extracts. Cell MotilCytoskeleton. 1997;36: 1–11. doi:10.1002/(SICI)1097-0169(1997)36:1 < > Journal of Chemical Physics
Nature
Reports on Progress in Physics ontraction in active microtubule networks Yu C C, Mogilner A and Gross S P 2011
Proceedings of the National Academy of Sciences
Proceedings of the National Academy of Sciences
Science
Biophysical Journal Journal of Integral Equations and Applications The Journal of ChemicalPhysics
Journal of Fluid Mechanics758