Emergence of Metachronal Waves in Active Microtubule Arrays
EEmergence of Metachronal Waves in Active Microtubule Arrays
Stephen E Martin, Matthew E Brunner, and Joshua M Deutsch ∗ Department of Physics, University of California, Santa Cruz, CA 95064, USA Voltaic Inc. 2150 Shattuck Ave, (Dated: June 20, 2018)The physical mechanism behind the spontaneous formation of metachronal waves in microtubulearrays in a low Reynolds number fluid has been of interest for the past several years, yet is still notwell understood. We present a model implementing the hydrodynamic coupling hypothesis fromfirst principles, and use this model to simulate kinesin-driven microtubule arrays and observe theiremergent behavior. The results of simulations are compared to known experimental observationsby Sanchez et al.[1, 2]. By varying parameters, we determine regimes in which the metachronalwave phenomenon emerges, and categorize other types of possible microtubule motion outside theseregimes.
Metachronal waves refer to the synchronization of thin,flexible appendages that result in large-scale wavelikeformations. These appear in biological systems at themacroscopic scale (e.g. the motion of millipede legs)and at the microscopic scale (e.g. cilia in air pathways).On the microscopic level, metachronal waves are essen-tial components of several critical biological processes,from motility in microorganisms to mucus clearance inhuman bronchial tubes [3, 4]. If cilia are unable to effec-tively move and synchronize, the results are often severe– especially if the disorder is genetic [3]. Research intophysical explanations for cilia beating [5], and of sponta-neous metachronal behavior in cilia is ongoing and stillnot well understood [6, 7], although many have suggestedthat this phenomenon can be explained from hydrody-namic coupling between cilia [8–11].Recently, in some remarkable experiments, Sanchezet al. demonstrated metachronal wave behavior in an in vitro system[1, 2]. Microtubules (MTs) aggregatedinto bundles of length 10 − µ m due to the additionof polyethylene glycol [12]. Many of these bundles at-tached at one end to a fixed boundary forming densearrays. When exposed to a solution containing clustersof kinesin and ATP, sustained metachronal wave behav-ior between MT bundles (similar to that displayed bycilia and flagella) was observed. MT bundles were con-strained to move between two glass slides. It is surpris-ing that a system with such few ingredients could developcomplex behavior that so closely resembles biological sys-tems, which are made up from a much more complicatedmachinery. Proteomic analysis indicate that eukaryoticcilia are composed of many hundreds of proteins [13].Some important details of this in vitro system are stillunclear, most notably whether the MTs in this experi-ment are unipolar or of mixed polarity. Opposite polarityMTs will move past each other, causing separation intounipolar bundles [14, 15]. We present analytical and nu-merical arguments for unipolarity in S-IV. The surprisingmechanism for the motion of unipolar bundles describedhere, has not previously been given [1, 2], and we believethat the agreement between our model and experiments F kin F stiff F base FIG. 1. Conceptual illustration of forces acting on a singlepolymer that are not due to hydrodynamic interactions. Theblue vectors indicate the buckling forces due to kinesin walk-ers (tangent to polymer), the red vectors show the directionand relative magnitude of stiffness forces (in the direction of d r /ds ), and the green arrows indicate a restorative forcekeeping the base of the polymer approximately perpendicularto the binding surface. provides further evidence to support our proposed expla-nation.The general mechanism proposed is quite similar to themodel used to describe and simulate cytoplasmic stream-ing in Drosophila oocytes, and the fact that it can beadapted as such is in many ways a testament to its pre-dictive power. A fair amount of attention has been paidin recent years to the understanding of how metachronalwaves form in such arrays [16–19]. However, such modelsoften rely on assumptions about individual MT (or cilia)beat patterns and/or on phenomenology. The model wepropose makes no such assumptions (beyond some minorsimplifications), relying on first-principles fluid mechan-ics calculations. This is important, as it is not clear whyone would want additionally to impose oscillatory be-havior on individual MTs given the lack of a well definedinternal structure. a r X i v : . [ q - b i o . S C ] J un We now present a model for the simulation of theSanchez et al. system. A similar method has beenused successfully to simulate cytoplasmic streaming in
Drosophila oocytes[20], and is based on theoretical workcompleted several decades ago regarding the calculationof Stokes flows created by a point force (stokeslet) nearno-slip boundaries[21, 22]. A conceptual explanation ofthis mechanism is given below, and further details re-garding theory and implementation are given in supple-mentary materials S-II and S-III.An illustration of how MT bundles are simulated isgiven in Fig. 1. Each MT bundle is modeled as a chainof monomers (i.e. polymer) which are held an approxi-mately fixed distance from one another by a spring force.The base of each polymer is anchored to a single point,and the polymer at the base is kept roughly perpendic-ular to the anchoring surface. Let the polymer be de-scribed by the curve r ( s ), where r (0) is the location ofthe polymer base, and s is the arc length. We give thepolymer a stiffness by implementing an energetic cost ofbending proportional to curvature squared, which impliesa local force at s proportional to d r /ds . Additionally,monomers feel a “buckling” force due to the drag fromthe walking kinesin F kin = − f k d r /ds , which is parallelto the polymer and toward the polymer base. f k will de-pend linearly on the speed of the kinesin and the solventviscosity. This force continually adds energy to the sys-tem (making it active), and has been shown to be a goodrepresentation of the average drag force due to kinesinwalking along the microtubule away from the polymerbase [20].This kind of model for a single chain was first em-ployed to understand glide assay dynamics in two di-mensions [23]. In three dimensions, periodic waves de-velop whose dynamics have been analyzed in detail [20],and related theoretical work has recently also been per-formed [24]. However, scaling can be used to get therelevant length and timescales [23]. The average radiusof curvature depends on the strength of the bucklingforce f k , and the elastic constant of a filament charac-terizing its stiffness k stiff . The radius of curvature overquite a wide range of parameters can be shown to be R = ( k stiff / ( βf k )) / , where β ≈ .
05. Likewise, theangular frequency is ω = f k / ( νR ), where ν is the hydro-dynamic drag coefficient per unit length. Although thereis a fairly large experimental uncertainty in parametersused to model a Drosophila oocyte, this model finds quitegood agreement with the experimental time and lengthscales. R was predicted to be 25 − µ m, close to the16 . ± . µ m observed. Likewise, the time scale was pre-dicted to be 203 − ± z , which is physically sensi-ble when considering the geometry of the Sanchez et al.experiments. In this experiment, MT bundles were ob-served between glass slides, with a height H , of approx-imately 10 µ m, creating a narrow channel for which fluidcan flow. For this reason, we adopt a two dimensionalgeometry. In addition, the no-slip boundaries of theplates have a large impact on the hydrodynamic forcesbetween monomers[21, 22], which we give explicitly in thesupplemental materials S-II. Other close-range contactforces were also used (repulsion from anchoring surface,monomer-monomer repulsion), and these are explainedin the supplementary materials S-III.We can now address at the qualitative level the mech-anism by which we propose the metachronal waves ob-served by Sanchez et al. form. As kinesin walk away fromthe polymer bases, the polymers will tend to buckle. If apolymer is isolated, this buckling will lead to corkscrewmotion or periodic waves[20]. When placed in an array,however, nearby polymers will exert hydrodynamic forceson one another that tend to synchronize their motion. Ifthese hydrodynamic forces are sufficiently strong, thiscan cause a transition from disordered motion to alignedMTs and correlated motion.Despite the fact that this model was developed to ex-plain and simulate cytoplasmic streaming, its mechanismcan be easily adapted for related biological phenomena.Indeed, when the conditions of the Sanchez et al. ex-periment are simulated in the same way, we observemetachronal waves. It is not clear if this is formallya transition or a more continuous crossover effect, butthe results found make strong predictions that should betestable experimentally. In the following, we present theresults of these simulations and discuss the required con-ditions for metachronal wave formation.Videos of select simulations are included in the Sup-plementary Materials. Fig. 2 shows some still framesof simulated arrays demonstrating metachronal wave be-havior in both the planar and circular geometries.We characterize the behavior of each system using the (a)(b) FIG. 2. Simulated metachronal wave formation for 128-polymer arrays in (a) circular and (b) planar geometries. Inboth cases, k oseen = 0 . k stiff = 10 . H = 1. correlation function for the chain ends x ( i, t ), C (∆ i, ∆ t ) = h ∆ x ( i + ∆ i, t + ∆ t )∆ x ( i, t ) i , (1)where ∆ x ( i, t ) = x ( i, t ) − h x ( i, t ) i . The average is performed over all chain indices i , andtime t , after a period of equilibration. Figs. 3, 4, and5 show correlation functions for planar and a circulargeometries (for the circular geometry, the polar angle θ is the position variable rather than x ). In the follow-ing, we will discuss these and examine how the systemresponds to changes in the strength of the interactiontensor, k oseen , k stiff , and height H . It should be notedthat changes in the viscosity or kinesin velocity and den-sity (that affect f k ), can be absorbed into a rescaling oftime, and of k stiff . k oseen , has a dramatic effect on the type of wave behav-ior seen, or whether it is observed at all. This strengthis a function of the hydrodynamic effects of kinesin walk-ing along microtubules, and will depend on their densityand speed, as explained in detail in Ref. [20]. Fig. 3shows the correlation results of three 128-polymer simu-lations in the same circular geometry shown in Fig. 2(a) for three different values of k oseen . There is an overallstrengthening of the metachronal behavior as k oseen isincreased from 0 . .
2. The sign of the slope reflectsthe initial conditions of the system. Long lived wavestravel predominantly in a single direction over long timesscales resulting in a slope of the crests of the correlationfunction that can either be positive or negative. Similarcrests are seen in the analysis of the real experimentaldata [1]. With this circular geometry, the correlationfunction must be periodic, which is why it rises againwhen i becomes large.The polymer stiffness k stiff also has an interesting ef-fect on metachronal wave formation. Fig. 4 shows thecorrelation functions for k stiff = 5.0, 10.0, and 20.0in a planar geometry. While Figs. 4(a-b) are quali-tatively similar, we do see an apparent decrease in themetachronal wavelength. Figs. 4(c-d) show that if thepolymer is made too stiff, no metachronal behavior isobserved at all. In general, planar geometry appears tocause more coherence in the motion of the different bun-dles, and the correlation function is dominated by motionat the longest lengths and time scales.The distance between plates, H , has a considerableeffect on the dynamics as well. Longer range, more co-herent motion is observed when H is larger, and shortrange, less coherent motion when H is small. See Fig. 5.This is to be expected due to the strong screening effectthat these boundary conditions impose. Smaller H re-duces the hydrodynamic coupling, causing a decrease incoherence.When comparing these results to those of Sanchez etal., we find that the basic features agree. The videos in-cluded in the supplemental materials qualitatively mimicthe experimental videos, and the experimental correla-tion analysis agrees quite well with the simulations. Moreimportantly, this agreement between theory and exper-iment was reached from first principles. We only use ahandful of forces in our simulations, and each force hasa physical justification for being used.There are potential shortcomings of this model thatmay result in some differences between experiment andtheory. The first is that the experiments observe bundlesof microtubules that taper away from their base. Thehydrodynamics are not expected to be uniform along thelength of a chain. In addition these bundles will, for shortenough times, behave like rigid material, but for longertimes, because they are connected through walking ki-nesin molecules, will behave more as individual mictro-tubules with a greatly reduced elastic constant. On thetime scales of the motion, we expect to be in the latterregime. However the details of the hydrodynamics andelasticity in these bundles is still not understood experi-mentally. In fact, as we mentioned earlier, the polarity ofindividual MT’s is not known experimentally, and argu-ments for their unipolarity are given in the supplemen-tary information, S-IV. But still, the basic mechanism of (a) (b) (c) (d) FIG. 3. Full correlation functions for circular geometry with H = 1 and k stiff = 10, with k oseen = 0.1, 0.2, and 0.3 (a-c,respectively). The correlation function at ∆ i = 0 for all of these values of k oseen are shown in (d). (a) (b) (c) (d) FIG. 4. Full correlation functions for planar geometry with H = 1 and k oseen = 0 .
1, with k stiff = 5.0, 10.0, and 20.0 (a-c,respectively). The correlation function at ∆ i = 0 for all of these values of k stiff are shown in (d). dynamic buckling due to kinesin drag, and metachronalwaves being generated by hydrodynamic coupling, is ro-bust over a wide parameter range, so we believe thatthese complications, aside from unipolarity, will not al-ter the basics of our explanation.At a more technical level, there are other things thatmay make a slight difference to the results here. The bun-dles are constrained to move only in the xy -plane, andwhile it is true that MT motion is nearly 2-dimensional,there is some room in the z − direction that MT bundlescan occupy. Additionally, this model does not accountfor the fluid boundary condition at the anchoring surface.This may introduce some errors if a monomer becomes close ( ∼ H ) to the anchoring plane. However, because ofthe screening effects of the plates, this should not alterthe behavior at distances large compared to the plate sep-aration. We have tested for this by adding image chargesto the planar case, and found that their effects on corre-lations are small, as expected.In conclusion, we have developed a model for the spon-taneous formation of wavelike behavior in active polymerarrays that only requires two ingredients: semi-flexiblechains tethered to a surface, and motors walking fromtheir bases to their tips. The hydrodynamics in theirconfined geometry gives rise to metachronal waves thatappear remarkably similar to what is observed experi- (a) (b) (c) (d) FIG. 5. Full correlation functions for circular geometry with k stiff = 10, k oseen = 0 .
2, and H = 1.0, 0.5, and 0.1 (a-c,respectively). The correlation function at ∆ i = 0 for all of these values of H are shown in (d). mentally [1, 2]. There is no need to posit additional mech-anisms that force individual bundles to oscillate. Thisall happens as a consequence of Newton’s laws and fluidmechanics, allowing us to gain a better understanding ofhow metachronal waves form with considerable predic-tive power. As such, we have examined new parameterspaces and have demonstrated boundaries between differ-ent types of metachronal behavior and regimes in whichno metachronal behavior exists. It would be of great in-terest to test these predictions experimentally. Given thesimplicity and robust nature of this mechanism, and theubiquity of microtubules and kinesin in cells, it gives onefurther impetus to look for other places in biology wherethis kind of behavior can be found.J.M.D. thanks Bill Saxton, Itamar Kolvin, Alex Tayar,and Zvonimir Dogic for useful discussions. S.E.M. waspartially supported by the ARCS Foundation. This workwas also supported by the Foundational Questions Insti-tute
Videos at https://sites.google.com/ucsc.edu/joshdeutsch/metachronal-videos?authuser=0
Supplementary video S1:
128 microtubule bundles (length 16) with kinesin walkers in a circular geometry in afluid chamber with k oseen = 0 .
1, chamber height H = 1 . k stiff = 10. Supplementary video S2:
128 microtubule bundles (length 16) with kinesin walkers in a circular geometry in afluid chamber with k oseen = 0 .
2, chamber height H = 0 . k stiff = 10. Supplementary video S3:
128 microtubule bundles (length 16) with kinesin walkers in a circular geometry in afluid chamber with k oseen = 0 .
2, chamber height H = 0 . k stiff = 10. Supplementary video S4:
128 microtubule bundles (length 16) with kinesin walkers in a circular geometry in afluid chamber with k oseen = 0 .
2, chamber height H = 0 . k stiff = 10. Supplementary video S5:
128 microtubule bundles (length 16) with kinesin walkers in a circular geometry in afluid chamber with k oseen = 0 .
2, chamber height H = 1 . k stiff = 10. Supplementary video S6:
128 microtubule bundles (length 16) with kinesin walkers in a circular geometry in afluid chamber with k oseen = 0 .
3, chamber height H = 1 . k stiff = 10. Supplementary video S7:
128 microtubule bundles (length 16) with kinesin walkers in a planar geometry in a fluidchamber with k oseen = 0 .
1, chamber height H = 1 . k stiff = 5. Supplementary video S8:
128 microtubule bundles (length 16) with kinesin walkers in a planar geometry in a fluidchamber with k oseen = 0 .
1, chamber height H = 1 . k stiff = 10. Supplementary video S9:
128 microtubule bundles (length 16) with kinesin walkers in a planar geometry in a fluidchamber with k oseen = 0 .
1, chamber height H = 1 . k stiff = 15. Supplementary video S10:
128 microtubule bundles (length 16) with kinesin walkers in a planar geometry in afluid chamber with k oseen = 0 .
1, chamber height H = 1 . k stiff = 20. Supplementary video S11:
Supplementary video S12:
Pillar of 9 microtubules of opposite polarities, green MT’s have minus ends on surface,and blue MT’s have plus ends on the surface. There are sliding boundary conditions on the surface. This shows asimulation in a regime with sufficiently weak attractive interactions, f a = 1, where there is a twisting motion insidethe pillar but then the minus microtubules suddenly slide off of the plus ones, finally lying close to parallel with theplane of attachment. Figure S-1. Fluid speeds as a function of distance ρ from the stokeslet F = ˆ ı . Solid curves are calculated using the fullinteraction tensor (S-2) and dashed lines are the far-field approximation (S-3). (a) u x ( ρ ) along the line y = 0; (b) u x ( ρ ) alongthe line x = 0; (c) u y ( ρ ) along the line y = x . II. THE QUASI-2D INTERACTION TENSOR
The interaction tensor used in simulations is that of a stokeslet enclosed by two infinite parallel plates, as derivedby Liron and Mochon[1]. In general, the interaction tensor G is defined as the relationship between the fluid flow u ( r )and the stokeslet F which causes this flow: u ( r ) = F · G ( r ) (S-1)We assume the system is embedded in a viscous fluid with viscosity µ . For computational efficiency, we assume allmonomers to be only in the xy -plane, with parallel plates at z = ± H/
2. This reduces a three-dimensional problem totwo dimensions, as (a) the stokeslet is located in the xy -plane, (b) the stokeslet’s direction has no z -component, and(c) we only concern ourselves with flows in the xy -plane (see Fig. S-2). For this arrangement, it can be shown fromLiron and Mochon’s general result that the interaction tensor a displacement r (and ρ ≡ | r | ) from a single stokeslet F at the origin reduces to G ( r ) = H πµρ (cid:26)(cid:20) (cid:16) ρH (cid:17) S − ρH I (cid:21) I + (cid:20) π (cid:16) ρH (cid:17) S + 12 ρH I − (cid:16) ρH (cid:17) I (cid:21) r ⊗ r ρ (cid:27) (S-2)where S ≡ ∞ X n =0 ( − n h(cid:0) ρH (cid:1) + n i / S ≡ π ρH ∞ X n =0 ( − n h(cid:0) ρH (cid:1) + n i / I ≡ Z ∞ ξJ (cid:16) ρH ξ (cid:17) tanh ξ sinh ξ − ξ dξI ≡ Z ∞ ξ h J (cid:16) ρH ξ (cid:17) − J (cid:16) ρH ξ (cid:17)i tanh ξ sinh ξ − ξ dξ Here, J n is the Bessel function of the first kind. Because S and S do not converge rapidly as defined above, we alsomake use of the Poisson sums S = ∞ X k =0 K h π (2 k + 1) ρH i S = ∞ X k =0 (2 k + 1) K h π (2 k + 1) ρH i H F z = 0 H2
Figure S-2. Illustration of the geometry for which interaction tensor is derived in II. While this is a three-dimensional system,we constrain polymers to the xy -plane. F a F b F a→b F b→a Figure S-3. Illustration of hydrodynamic forces between two example monomers in a planar polymer array. The green forcesare the sum of non-hydrodynamic forces on the monomer (and by extension the force the monomer exerts on the surroundingfluid). ~F a → b and ~F b → a are the hydrodynamic forces on monomer b due to ~F a and the hydrodynamic force on monomer a dueto ~F b , respectively. where K n is the modified Bessel function of the second kind.In the far field, it can be shown that (S-2) approaches G ( r ) ≈ − H πµρ (cid:18) I − r ⊗ r ρ (cid:19) (S-3)Fig. S-1 shows plots of u ( r ) at selected locations, and compares the exact value from (S-2) to the far-field approxi-mation from (S-3).We can now make some conceptual observations regarding this interaction tensor and how it compares to theboundary-free Oseen tensor G : G ( r ) = 18 πµr (cid:18) I + r ⊗ r r (cid:19) First, we immediately notice a 1 /r dependence (rather than 1 /ρ ). This means forces without boundaries tend to bemore long-range, and boundaries result in long-range screening. Second, G is always positive, whereas this is not truefor the interaction tensor used here. One key implication of this is that flows created by a stokeslet are often flowingopposite its direction (e.g. Fig. S-1b). Both of these qualities may enhance metachronal behavior in the confinedsystem. Screening means that interactions between nearby polymers are most important, creating a “domino effect”from one polymer to the next rather than having motion more influenced by long-range interactions. The creation ofopposing flows means (among other things) that if one polymer is moving toward the anchoring surface, it may exerta force on many of its neighboring polymers away from the anchoring surface. This encourages wavelike behaviorrather than uniformity of beating motion. III. SIMULATION METHODS
The algorithm we implement is built on work that was used to simulate the mechanism behind cytoplasmic streamingin
Drosophila oocytes [2], and many of the methods and equations below are explained in detail in these papers. Thissoftware simulates an array of active microtubules tethered to a plane that works as follows and is explained in furtherdetail below.1. After an array of polymers is initialized, forces on all monomers are summed (described below, also see Fig. 1of main paper) and monomer position and velocity are updated using time step dt .2. This motion initiates complex flow in the surrounding fluid. The fluid flow is not simulated directly, but theresulting hydrodynamic forces from this flow are calculated via an Oseen tensor with corrections by Blake[3].This is illustrated in Fig. S-3.3. Forces on each monomer are summed, and monomer position and velocity are updated accordingly.4. Once updated, steps 2-3 are repeated.In the present work there were these differences:1. N Microtubules are confined to the xy -plane, with polymer bases separated by a distance l tethered either to aflat plate at y = 0 or to a circular boundary. For all presented results, N = 128. The geometry of this is shownin Fig. S-2.2. At the tethering point, a potential was added in order to keep the base monomer approximately orthogonal tothe boundary.3. Rather than the Blake correction to the Oseen tensor, we use the simplified Liron/Mochon interaction ten-sor described in supplemental Section II. We also investigated varying the hydrodynamic coupling parameter k oseen ≡ / (8 πµ ).Now we describe how the above was accomplished in more detail. Each polymer is composed of n = 16 monomers.The i th monomer position r i is updated a using a fourth order Runge Kutte integration of the equation d r i dt = u ( r i ) − k kin ( r i − − r i +1 ) (S-4)where dt is the time step (set to 0.003), k kin (set to 0.2) controls the strength of the kinesin force tangent to thepolymer ( F kin in Fig. 1), and u ( r i ) is the fluid velocity due to the motion of all other monomers as given by EquationS-1 and S-2 (which imparts the forces F a → b in Fig. S-3): u ( r i ) = X j = i F j · G ( r i − r j ) (S-5)Here, F j is the total force on the fluid due to the j th monomer. Because there are are no inertial effects when Re (cid:28) F j = T j + C j + Q j , (S-6)where • T j = k spr [( | r j − | − ‘ ) ˆ r j − + ( | r j + | − ‘ ) ˆ r j + ]with r j ± ≡ r j ± − r j , is the spring force keeping monomer separation approximately constant. For our simula-tions, k spr = 100 and ‘ = 1. In these simulations the separation between polymer bases defines above, l is equalto 4 ‘ . • C j = k stiff (2 r i − r i +2 − r i − )is the stiffness force which resists polymer bending. k stiff is varied in our simulations, but typically 5 ≤ k stiff ≤ • Q j = P j + B j + W j + P k H jk is the sum of miscellaneous conditional forces: – P j = k pin ( r j − h ˆ ) if ( j mod n ) = 1is the force on the base monomer of each polymer chain keeping it pinned to the anchoring surface. Forour simulations, we set k pin = 100 and h = 1. – B j = k pin ( r j − r j − − ‘ ˆ ) if ( j mod n ) = 2is the force on the second monomer in each polymer chain, keeping the base of each polymer approximatelyorthogonal to the anchoring surface ( F base in Fig. 1). For our simulations, we set k base = 100. – W j = k wall (cid:20) − (cid:16) d wall y j (cid:17) (cid:21) ˆ if y j < d wall is the repulsive force exerted by the anchoring plane on any monomer that gets close to the wall. For oursimulations, we set d wall = 0 . k wall = 100. – H jk = k rep (cid:20) − (cid:16) d rep | r j − r k | (cid:17) (cid:21) ( r j − r k ) if | r j − r k | < d rep is the repulsive force between monomers that are very close to one another. For our simulations, we set d rep = 0 . k rep = 1. IV. ANALYSIS OF UNIPOLARITY
The work of Sanchez et al. [4, 5] consists of a mixture of biotin-labeled kinesin-1 motors bound together to formclusters using multimeric streptavidin and taxol stabilized microtubules in a polyethylene-glycol solution with ATP.These form bundles of microtubules, some of which are adsorbed to air-water or air-glass interfaces, that point outfrom the interface forming a lawn of microtubule bundles. These bundles are flexible and show bending similar towhat is seen in the simulations described here in both the time scales, length scales, and correlations between differentbundles.The question that is not answered in the experimental work is the directionality of the microtubules inside a bundle.The microtubules forced into bundles by the polyethylene glycol (PEG) could be of mixed polarity so that some havetheir minus ends at the interface while others have their plus ends there. We will refer to microtubules with differentorientations as having different “polarities”, minus-ends against the interface as “minus” and those with oppositepolarity as “plus”.The problems with having a mixed polarity bundle are two fold. The first is that for a wide range of experimentalparameters, we expect mixed polarity bundles to be unstable [6, 7]. The second problem is that it is not clear thatmixed polarity bundles can give rise to the motion seen experimentally. We will analyze both problems below.
A. Instability of mixed polarity bundles
The first problem is that adjacent microtubules with different polarities will be linked by kinesin clusters that willapply equal and opposite forces to them. This will cause the minus microtubules to be pushed toward the interface,and the plus ones away from it. The forces from the kinesin act in parallel on a microtubule over its length which isof order 10 µm . The forces that these cause can be competitive with depletion forces caused by the PEG as we willnow see. A full analysis of this is not possible without more information about the details of the system such as thedensity of kinesin clusters and chain lengths of the PEG. However we can do a calculation to show that even withvery modest assumptions concerning kinesin density, expulsion of plus microtubules will take place.Depletion forces exert an osmotic pressure on microtubules and filaments. Each polymer excludes a roughly sphericalregion of order its radius of gyration R g . Entropic forces favor the separation of microtubules into bundles becauseless volume is excluded by the PEG. We will estimate the force acting on a single microtubule protruding from abundle. PEG is depleted in a region of size R g around the microtubule. The increase in free energy per unit areacaused by this depletion is of order pR g where the osmotic pressure is p = k B T ρ , and ρ is the number of polymersper unit volume. The increase in free energy dF , in raising the microtubule by a height dz , is dF = (2 πR m dz ) px .Here R m is the microtubule radius. If we assume that the polymers are close-packed around the microtubule to getthe maximum effect, then ρ = 1 / (4 πR g / f = dF/dz = (3 / R m k B T /R g . R m ≈ nm and conservatively taking R g = 1 nm , which is quite small for PEG, f = 81 pN The stall force ofkinesin is approximately 5 pN [8]. So only 16.2 kinesins are needed to overcome the depletion forces and expel thismicrotubule from the bundle.The minimum separation of kinesin on a microtubule is 8 nm and there are 13 tracks around its circumference.Because kinesin has a strong affinity for microtubules we expect a high density of bound kinesin. Therefore 16kinesins contributing to the force over a distance of 10 µm is over three orders of magnitude less dense than themaximum density attainable. This suggests that for a wide range of parameters, the microtubule bundles will becomeunipolar with minus-ends against the interface. B. Model of mixed polarity bundles
The second problem is that it is not clear that a mixed polarity bundle can give rise to the motion seen inexperiment. Here we analyze this possibility by using simulation methods similar to what was used previously tounderstand molecular motor dynamics [9]We assume that the microtubules are inextensible and that opposite polarity microtubules apply forces in equaland opposite directions. We discuss the different forces separately.First there is an effective attractive interaction between microtubules independent of their polarities induced bythe presence of PEG polymers. We choose a short range force so the monomers separated by a distance r within arange σ s will feel an attractive force due to depletion forces as discussed above. To simplify the expressions we use anormalized unitless distance ∆ ≡ r /σ s . The force between any two monomers for ∆ < f attr = f a ∆ (1 − ∆ ) r (S-7)where f a is the strength of the attractive interaction. The reason for choosing this functional dependence on ∆ wasto produce a force that was close to constant for ∆ < .
6, and then drop smoothly to zero, so as to work well withthe Runge Kutte algorithm.Second, we introduce an even shorter range repulsion between monomers that diverges at a hard core radius σ h and goes to zero at σ s : f rep = f r (cid:18) r − σ h − σ s − σ h (cid:19) r (S-8)where f r is the strength of the repulsive interaction.Third, we introduce an equal and opposite forces between monomers on opposite polarity microtubules that arewithin a distance σ s . The direction of the force is as follows. We compute the tangents to both monomers as( r i +1 − r i − ) /
2. Then we choose the direction t , to be the average of these two tangents. The magnitude of thekinesin force is f kin = f k (1 − ∆ ) (S-9)where f k is similar the symbol used previously in the main text and denotes the magnitude of the kinesin force.These forces are added to the elastic forces, viscous drag, and tension that must be introduced to conserve linklength and the equation of motion is iterated using a method for updating chains with constant link length [10, 11].We also tried two separate kinds of boundary conditions. First, tethering the chains to fixed points on the surfacewhich we will call “fixed” boundary conditions. Second, confining the chain ends to a two dimensional plane butletting the ends move within that plane, which we will call “sliding” boundary conditions.We tried a wide range of parameters, of different elastic constants, attractive interactions, number of microtubules,and boundary conditions. What we found is now summarized.For two chain bundles of opposite polarity we did find a set of parameters which showed movement of the bundlewith: f r = 10 . σ s = 2, σ h = 1, f a = 3, f k = 0 . k stiff = 100, and chain length of 20, see supplemental move S11.For larger bundle sizes, e.g. 9 chains, we did not find anything similar to experiments. With fixed boundaryconditions, and started as a pillar of parallel microtubules with slightly randomized directions, the chains would settledown to a pillar shape that would not change with time for sufficiently small attractive interactions f a , but whenthis became greater than a certain value that depends on elastic constant and other parameters, it would suddenlycollapse into a ball because this is more highly favored energetically.When we chose sliding boundary conditions, and for sufficiently weak attractive interactions, f a = 1 there was aregime where there was twisting motion inside the pillar but then the minus microtubules would suddenly slide off ofthe plus ones, finally lying close to parallel with the plane of attachment, see supplemental movie S12. It thereforeappears that a two microtubule bundle moves because of a strong anisotropy in forces seen in cross sections. In largerbundles, the forces through the bundle are more homogeneous which acts to stabilize them.We conclude that by direct physical modeling of a mixed polarity bundle, it is not clear if there are any reasonableparameters which show motion similar to what is seen in the experiments of Sanchez et al [4, 5].Note that the elastic constant of a microtubule in a bundle will depend strongly on the rate at which it is bent. Forvery short times, the bonds between different microtubules caused by kinesin binding will be fixed in position givingthe bundle the elastic constant of a cylinder of radius R which is ∝ R . However the oscillations here take place onminute timescales. In that case the individual kinesin molecules have velocities of order 1 µm/s so they unbind andmove very far on this time scale. This allows neighboring microtubules to move relative to each other, to eliminatestress. Therefore on sufficiently long timescales, this reduces the elastic constant of a microtubule to that of one inisolation. [1] N. Liron and S. Mochon, Journal of Engineering Mathematics , 287 (1976).[2] C. E. Monteith, M. E. Brunner, I. Djagaeva, A. M. Bielecki, J. M. Deutsch, and W. M. Saxton, Biophys. J. , 2053(2016).[3] J. R. Blake, Mathematical Proceedings of the Cambridge Philosophical Society , 303 (1971).[4] T. Sanchez, D. Welch, D. Nicastro, and Z. Dogic, Science , 456 (2011).[5] T. Sanchez and Z. Dogic, in Methods in enzymology , Vol. 524 (Elsevier, 2013) pp. 205–224.[6] K. Kruse and F. J¨ulicher, Physical Review Letters , 1778 (2000).[7] T. B. Liverpool and M. C. Marchetti, Physical Review Letters , 138102 (2003).[8] E. Meyh¨ofer and J. Howard, Proceedings of the National Academy of Sciences , 574 (1995).[9] J. M. Deutsch and S. E. Martin, Macromolecules , 6703 (2015).[10] J. Deutsch, Science , 922 (1988).[11] J. Deutsch and T. Madden, The Journal of Chemical Physics90