A model for the orientational ordering of the plant microtubule cortical array
aa r X i v : . [ c ond - m a t . s o f t ] M a y A model for the orientational ordering of the plant microtubule cortical array
Rhoda J. Hawkins, ∗ Simon H. Tindemans, † and Bela M. Mulder FOM Institute for Atomic and Molecular Physics,Science Park 113, 1098 XG, Amsterdam, The Netherlands
The plant microtubule cortical array is a striking feature of all growing plant cells. It consists of amore or less homogeneously distributed array of highly aligned microtubules connected to the innerside of the plasma membrane and oriented transversely to the cell growth axis. Here we formulate acontinuum model to describe the origin of orientational order in such confined arrays of dynamicalmicrotubules. The model is based on recent experimental observations that show that a growingcortical microtubule can interact through angle dependent collisions with pre-existing microtubulesthat can lead either to co-alignment of the growth, retraction through catastrophe induction orcrossing over the encountered microtubule. We identify a single control parameter, which is fullydetermined by the nucleation rate and intrinsic dynamics of individual microtubules. We solve themodel analytically in the stationary isotropic phase, discuss the limits of stability of this isotropicphase, and explicitly solve for the ordered stationary states in a simplified version of the model.
PACS numbers: 87.10.Ed, 87.16.ad, 87.16.Ka, 87.16.LnKeywords: cortical array, cytoskeleton, microtubules, alignment, model
I. INTRODUCTION
Most plant cells grow by uniaxial expansion. Estab-lishing and maintaining this characteristic anisotropicgrowth mode requires regulatory mechanisms that arerobust, and, in addition, sensitive to the cell geometry.A major role in this process is played by microtubules,highly dynamic filamentous protein aggregates that formone of the components of the cytoskeleton of all eukary-otic organisms (see Ref. [1], chapter 16). In growingplant cells microtubules are confined to a thin layer ofcytoplasm just inside the cell plasma membrane. Herethey form the so-called cortical array, an ordered struc-ture formed by highly aligned (bundles of) microtubulesoriented transversely to the growth direction [2]. Thisstructure is unique to plant cells and there is evidencethat it controls the direction of cell expansion by guid-ing the mobile transmembrane protein complexes thatdeposit long cellulose microfibrils in the plant cell wall[3]. These cellulose microfibrils are the main structuralelements of the cell wall, which, for mechanical reasons,are also transversely oriented to the cell axis in growingcells [4]. An in vivo image of the array and a schematicare shown in Figure 1.
In vivo imaging of microtubules la-belled with fluorescent proteins in plant cells by severalgroups has shown how the cortical array is establishedboth following cell division and after microtubule depoly-merizing drug (oryzalin) treatment [2, 3, 5, 6, 7, 8]. Inthese studies microtubules are seen to nucleate at the cor-tex and then develop from an initially disorganized stateinto the transverse ordered array over a time period on ∗ Current address: UMR 7600, Universit´e Pierre et MarieCurie/CNRS, 4 Place Jussieu, 75255 Paris Cedex 05 France; Elec-tronic address: [email protected]; The first two authorscontributed equally to this work. † The first two authors contributed equally to this work. a)b)
Figure 1: (a) Fluorescently labelled microtubules in a TobaccoBY-2 cell expressing GFP:TUA6. Image courtesy of JelmerLindeboom, Wageningen University. (b) Schematic represen-tation of a plant cell cortical array. the order of one hour. The nature of the self-organizationprocess by which the specific spatial and orientationalpatterning of this cytoskeletal structure is achieved is asyet only partially understood and forms the subject ofthis work.An important aspect of the problem is the nature of lo-calization of the microtubules to the cortical region. Flu-orescence recovery after photo-bleaching (FRAP) experi-ments by Shaw et al. [7] showed that the microtubules arefixed in space, so any apparent mobility of microtubulesis due to ‘treadmilling’, the process of simultaneous poly-merization at one end and depolymerization at the otherend. So, cortical microtubules do not translate or rotateas a whole. The same authors also did not detect de-tachment or (re)attachment of microtubules to the cellcortex, apart from some growing ends of single micro-tubules moving out of focus and found no evidence formotors working in the cortical array. These experimentsindicate that the microtubules in the cortical array arefixed to the inside of the cell membrane. Electron mi-croscopy has also shown cross-bridges between corticalmicrotubules and the membrane [9]. It is therefore widelyassumed that there are linker proteins that anchor themicrotubules through the plasma membrane to the rigidcell wall, although their molecular identity is under de-bate [10, 11, 12, 13, 14]. Since the cortical microtubulesare effectively confined to a 2D surface, they can interactthrough ‘collisions’ that occur when the polymerizing tipof a growing microtubule encounters a pre-existing mi-crotubule. The resulting dynamical interaction eventswere first characterized by Dixit and Cyr [15] in tobaccoBright Yellow-2 (BY-2) cells. They observed three dif-ferent possible outcomes: (i) zippering : a growing micro-tubule bending towards the direction of the microtubuleencountered, which occurs only when the angle of inci-dence is relatively small ( . ◦ ) (ii) induced catastrophe :an initially growing microtubule switching to a shrinkingstate and retracting after the collision, an effect predom-inant at larger angles of incidence and (iii) cross-over : agrowing microtubule ‘slipping over’ the one encounteredand continuing to grow in its original direction.There are clearly many coupled mechanisms at workin this complex biological system contributing to the as-sembly and maintenance of this microtubule cortical ar-ray structure. We are interested in understanding whatare the main contributing factors and how their interplayleads to the observed orientational ordering. With thisaim we develop a coarse-grained model, incorporating allthe effects discussed above. Our emphasis on the plant-specific biological mechanism of the ordering in the cor-tical array distinguishes our approach from earlier work.Over the years, various models for self-organizationof cytoskeletal filaments (and polar rods in general)have been proposed [16, 17, 18, 19], and the model byZumdieck et al. [20] was applied to the plant cortex.However, in each of these models the filaments are as-sumed to have rotational and, in most cases, translationaldegrees of freedom. This is inconsistent with the fact thatthe plant cortical microtubules are stably anchored. In-spired by the experimental results of Dixit and Cyr [15],Baulin et al. [21] were the first to report on a two di-mensional dynamical system of treadmilling and collidingmicrotubules. Their focus was on establishing the mini-mal interactions needed to generate dynamical alignment.Using stochastic simulations they showed that a pausingmechanism, whereby a growing microtubule stalls againstanother microtubule until the latter moves away, can in-deed lead to ordering. Stalling, however, is not oftenobserved in the cortical array. Moreover their modellacks dynamic instabilities, i.e. catastrophes, both spon-taneous and induced, and rescues, and employs a formof deterministic microtubule motion, which is arguablyunrealistic in view of the observed dynamics.The outline of the paper is as follows. In section II weformulate our course-grained model starting from a de-scription of the dynamics of individual microtubules. Wethen construct the continuity equations that couple thedensities of growing, shrinking and inactive microtubule segments due to the intrinsic and collisional dynamics.In the steady state we can reduce the initial set of equa-tions to four coupled non-linear integral equations. Wethen perform a dimensional analysis to identify the rele-vant control parameter of the system. In section III wepresent the results of our model. We first solve the modelanalytically in the isotropic stationary state. Using a bi-furcation analysis we then determine the critical valuesof the control parameter at which the system develops or-dered solutions. We interpret these results in terms of thephysical parameters of microtubule segment length andmesh size. Finally, we formulate a minimal model withrealistic interaction parameters that we can solve numeri-cally to obtain all stationary ordered solutions. We closeby giving arguments for the stability of these solutions.The paper concludes with a discussion section. An ap-pendix outlines further details of the numerical solutiontechnique employed. II. MODELA. Description of the microtubules and theirdynamics
As described in the introduction we confine the config-uration of the microtubules to a 2D plane. Since collision-induced zippering events can cause microtubules to bendalong the direction of preexisting ones, we divide each mi-crotubule into distinct segments with a fixed orientation.We treat these segments as straight rigid rods. This isjustifiable since the persistence length l p of microtubulesis long ( ∼ mm) compared to the average length of a mi-crotubule ( ∼ µ m) and, as mentioned above, adhesionto the plasma membrane further inhibits thermal motion.Microtubules are known to be dynamic in that they arecontinually growing or shrinking by (de)polymerization.We use the standard two-state dynamic instability modelof Dogterom and Leibler [22] which assumes that eachmicrotubule has a ‘plus’ end, located on the final segmentof each microtubule, that is either growing (labelled by +)with speed v + or shrinking (labelled by − ) with speed v − .This plus end can switch stochastically from growing toshrinking (a so-called ‘catastrophe’) with rate r c , or fromshrinking to growing (a so-called ‘rescue’) with rate r r ina process known as dynamic instability.We model the creation of new microtubules with aconstant, homogeneous, isotropic nucleation rate r n inthe plane of the 2D model. In vivo nucleation appearsto occur at the cortex and has been observed to occurin random orientations unattached to pre-existing micro-tubules [7]. Although microtubules have also been ob-served to nucleate from by γ -tubulin complexes bindingto pre-existing microtubules [2, 23, 24] we ignore thispossibility for simplicity’s sake. By the same token wedisregard the possibility of the shrinking of microtubulesat their less active ‘minus’ end, leading to motion throughthe ‘treadmilling’ mechanism [25]. The initial segment of r n r c r r P c ( θ) z ( θ) -P c ( θ) P z ( θ)θ zipperinginducedcatastrophecrossover Collision dynamicsNucleation Dynamic instability v+v-
Segment counting
21 2 13 4active segmentsinactive segments catastropherescue
Figure 2: Schematic representation of the model interaction.The microtubule of interest is drawn in black and other micro-tubules that it encounters are in grey. The active segments ofthe black microtubule have an arrow head indicating growthor shrinkage whereas inactive segments end in the junctionwith the following segment depicted by a dot. See also thedescription of the parameters in table I. each microtubule therefore remains attached to the nu-cleation point in our model.We call the final segment of a microtubule, which con-tains the growing or shrinking tip, active and all theremaining ones, which do not change their length, in-active (labelled by 0). A cartoon of an an individualmicrotubule according to these definitions is depicted infigure 2. When a microtubule collides with another micro-tubule and experiences a zippering event, its active seg-ment is converted into an inactive segment, and a new ac-tive segment is created alongside the encountered micro-tubule. The inverse can also occur: if the active segmentshrinks to zero length, a previously inactive segment inanother direction can be reactivated. An induced catas-trophe event simply causes the growing active segmentto become a shrinking one, as is the case for spontaneouscatastrophes. Finally, a crossover results in the growing active segment continuing to grow unperturbed.In figure 3 we present the relative probabilities for zip-pering, induced catastrophes and crossovers as a result ofcollisions between microtubules, based on the data pro-vided by Dixit and Cyr [15]. We assume that there areno microtubule polarity effects, as they were not reported.The probabilities P z ( θ − θ ′ ), P x ( θ − θ ′ ) and P c ( θ − θ ′ ) forzippering, cross-overs and induced catastrophes respec-tively are therefore even functions of the angle difference θ − θ ′ defined by their values on the interval [0 , π ]. Inthis article we will use only the following minimal set ofproperties, which are qualitatively supported by the data.Firstly, zippering becomes less likely for increasing angleof incidence, and is effectively zero at θ − θ ′ = π , whichis reasonable as the energy associated with bending themicrotubule increases with angle. Secondly, the proba-bility for induced catastrophes monotonically increaseswith increasing angle of incidence, reaching a maximumat θ − θ ′ = π , consistent with observations that indicatethat a microtubule which is hindered in its growth willundergo a catastrophe at a rate that depends inverselyon its growth speed [26]. æ æ æ æ æ æ æ æ æà à à à à à à à à Π (cid:144) Π (cid:144) DΘ p r ob a b ilit y P c P z P x Figure 3: Probabilities for zippering, cross-overs and catas-trophes as deduced from the observations of [15] (combineddata from MBD-DsRed and YFP-TUA6 labelling). Light greyshaded region: fraction of zippering events. Dark grey shadedregion: fraction of induced catastrophes. White region: frac-tion of crossovers. Every data point is located at the centreof the corresponding bin, and the shaded regions have beenextended to the boundaries using horizontal lines. The cor-responding lowest order Fourier coefficients of the interactionfunctions are: ˆ c = 0 .
59, ˆ c = − .
36, ˆ z = 0 .
24 (computedusing numerical integration of the product of | sin θ | and apiecewise linear interpolation of the data). B. Continuum model
Since there are many ( ≈ - 10 ) microtubules, eachof which can have multiple segments, in the cortical ar-ray of a typical interphase plant cell we treat the systemusing a coarse-grained description. In this approach, in- Table I: Overview of all parameters and variables in natural dimensionsParameters Description Dimensions v + growth speed [length] / [time] v − shrinkage speed [length] / [time] r c catastrophe rate 1 / [time] r r rescue rate 1 / [time] r n nucleation rate [length] − [time] − P c ( θ ) probability of induced catastrophe upon collision 1 P z ( θ ) probability of zippering upon collision 1Synthetic parameters g = r r v − − r c v + growth parameter 1 / [time] u = 1 + v + v − speed ratio 1 c ( θ ) = sin( | θ | ) P c ( θ ) ←→ { ˆ c n } effective catastrophic collision probability 1 z ( θ ) = sin( | θ | ) P z ( θ ) ←→ { ˆ z n } effective zippering probability 1Dependent variables k ( θ ) microtubule length density [length] − [radian] − l ( θ ) average microtubule segment length [length] ˘ m + i ( l, θ ) , m − i ( l, θ ) , m i ( l, θ ) ¯ density of growing/shrinking/inactive segmentswith length l and direction θ [length] − [radian] − stead of individual microtubules, we consider local den-sities of microtubule segments. This approximation isreasonable as long as the length scale of an individualmicrotubule segment is small compared to the linear di-mensions of the cell. From the outset we assume that thesystem is (and remains) spatially homogeneous, and wewill eventually restrict ourselves to the steady-state solu-tions. In order to deal with the memory effect caused bythe isotropic nucleation, followed by subsequent reorient-ing zippering events, we need to keep track of the segmentnumber i , which starts at 1 for the segments connectedto their nucleation site and increases by unity at eachzippering event. Our fundamental variables are there-fore the areal number densities m σi ( l, θ, t ) of segments instate σ ∈ { , − , + } with segment number i having length l and orientation θ (measured from an arbitrary axis) attime t. These densities obey a set of master equationsthat can symbolically be written as ∂ t m + i ( l i , θ i , t ) = Φ growth + Φ rescue − Φ spont. cat. − Φ induced cat. − Φ zipper (1a) ∂ t m − i ( l i , θ i , t ) = Φ shrinkage − Φ rescue + Φ spont. cat. + Φ induced cat. + Φ reactivation (1b) ∂ t m i ( l i , θ i , t ) = +Φ zipper − Φ reactivation (1c)The flux terms Φ event couple the equations for the grow-ing, shrinking and inactive segments and between differ-ent values of i . Equations (1) must be supplemented bya set of boundary conditions for the growing segments at l = 0. For the initial segment ( i = 1) this reflects theisotropic nucleation of new microtubules, given by v + m +1 ( l = 0 , θ, t ) = r n π , (2) where r n is nucleation rate. For the subsequent segments i >
1, this ‘nucleation’ of growing segments is the resultof the zippering of segments with index i −
1. Defining ϕ zipper ( θ i − → θ i , l i − , t ) as the flux of i -segments withangle θ i and length l i zippering into angle θ i +1 at time t (this will be made explicit in equation (13)), we obtainthe boundary condition v + m + i ≥ ( l i = 0 , θ i , t ) = Z d l i − Z d θ i − ϕ zipper ( θ i − → θ i , l i − , t ) . (3)Generally, this leads to a qualitatively different bound-ary condition for every value of i . The model thereforeconsists of an infinite set of coupled equations, three forevery value of i . However, in section II C we will showthat in the steady state, this can be reduced to a finite setby summing over all segment indices i . In the following,we derive explicit expressions for each of the flux termsΦ event .
1. Growth and shrinkage terms: Φ growth , Φ shrinkage Φ growth in Equation (1a) corresponds to the length in-crease of the growing segments. For segment growth inisolation, the length increase in a small time interval δt is given by v + δt, where v + is the growth velocity, and wehave m + ( l + v + δt, θ, t + δt ) = m + ( l, θ, t ). By expandingthe left hand term to first order in δt, we find ∂ t m + i ( l, θ, t ) = − v + ∂ l m + i ( l, θ, t ) ≡ Φ growth (4)A similar derivation yields that ∂ t m − i ( l, θ, t ) = v − ∂ l m − i ( l, θ, t ) ≡ Φ shrink (5)where v − is the shrinking velocity.
2. Dynamic instability terms: Φ rescue , Φ spont. cat. Φ rescue and Φ spont. cat. in equations (1a) and (1b) cor-respond to the fluxes due to the spontaneous rescuesand spontaneous catastrophe respectively and are simplygiven by Φ rescue = r r m − i ( l i , θ i , t ) (6)Φ spont. cat. = r c m + i ( l i , θ i , t ) (7)where r r is the spontaneous rescue rate and r c is thespontaneous catastrophe rate.So far, we have described the first three terms of equa-tions (1a) and (1b) (growth, shrinkage and dynamic in-stability terms). Together, these fully describe a systemof non-interacting microtubules, in which also the bound-ary condition (3) vanishes due to the absence of zippering.In this special case we recover the well-known equationsintroduced by Dogterom and Leibler [22] (for i = 1).
3. Interaction terms: Φ induced cat. , Φ zipper An interaction can occur when a growing active micro-tubule segment collides with another segment, irrespec-tive of the latter’s state and length. This prompts thedefinition of the total length density k ( θ, t ) of all micro-tubule segments in direction θ at time t , given by: k ( θ, t ) = X i Z d l i l i ( m + i ( l i , θ, t )+ m − i ( l i , θ, t )+ m i ( l i , θ, t ))(8)The density of collisions of a microtubule segment grow-ing in direction θ with other segments in direction θ ′ isdetermined by the geometrical projection | sin( θ − θ ′ ) | k ( θ ′ , t ) , (9)where the presence of the sin( θ − θ ′ ) factor ensures thecorrect geometrical weighting reflecting the fact that par-allel segments do not collide. When a collision occurs,one of the three possible events, induced catastrophe( c ) , zippering ( z ) or cross-over ( x ) occurs, with probabil-ities P c ( θ − θ ′ ) , P z ( θ − θ ′ ) and P x ( θ − θ ′ ) respectively.These probabilities can (and in-vivo do, see Figure 3) de-pend on the relative angle θ − θ ′ between the incomingsegment and the ‘scatterer’. For convenience sake weabsorb the geometrical factor | sin( θ − θ ′ ) | into the prob-abilities, by defining f ( θ − θ ′ ) = | sin( θ − θ ′ ) | P f ( θ − θ ′ )for all events f ∈ { c, z, x } . The incoming flux of growingmicrotubule segments with given segment number, lengthand orientation is given by v + m + i ( l, θ, t ) . With these def- initions we can write the interaction terms asΦ induced cat. = v + m + i ( l i , θ i , t ) Z d θ ′ c ( θ i − θ ′ ) k ( θ ′ , t )(10)Φ zipper = v + m + i ( l i , θ i , t ) Z d θ ′ z ( θ i − θ ′ ) k ( θ ′ , t )(11)The analogous term for crossovers is not used, becausethe occurrence of a crossover event has no effect on thegrowth of a microtubule.
4. Reactivation term: Φ reactivation The reactivation term Φ reactivation corresponds to theflux of active microtubule segments, with segment num-ber i + 1 , that, by shrinking to length zero, reactivatea previously inactive segment, effectively undoing a pastzippering event and so create a new, shrinking, activesegment with segment number i. The incoming flux ofof such segments coming from a given direction θ i +1 isgiven by v − m − i +1 ( l i +1 = 0 , θ i +1 , t ). The reactivation fluxis given byΦ reactivation = Z dθ i +1 v − m − i +1 ( l i +1 = 0 , θ i +1 , t ) p unzip ( θ i , l i | θ i +1 , t ) , (12)where the ‘unzippering’ distribution p unzip ( θ i , l i | θ i +1 , t )gives the probability that the shrinking microtubule reac-tivates an inactive segment with orientation θ i and length l i . This distribution will be determined below.A microtubule that has zippered will take a certainamount of time τ to undergo a catastrophe and returnto the zippering location, where τ is a stochastic vari-able. The unzipppering flux from direction θ i +1 at time t consists of microtubules that had zippered at a range oftimes t − τ and have now returned to the zippering loca-tion. This implicitly defines an originating time distribu-tion p origin ( t − τ | θ i +1 , t ) for the returning microtubules.Furthermore, because the evolution of a microtubule be-tween the zippering event and its return to the same lo-cation does not depend on the previous segments, thesegment that is re-activated by a microtubule returningto the zippering position after a time τ should be selectedproportional to the ‘forward’ zippering flux at time t − τ .The forward flux ϕ zipper ( θ i → θ i +1 , l i , t ) of microtubuleswith length l i and angle θ i zippering into angle θ i +1 isdefined in accordance with equation (11) as ϕ zipper ( θ i → θ i +1 , l i , t ) = v + m + i ( l i , θ i , t ) z ( θ i − θ i +1 ) k ( θ i +1 , t ) . (13)At each of the originating times t − τ , the distributionof microtubules that zipper into the direction θ i +1 withlength l i and orientation θ i is given by p zip ( θ i , l i | θ i +1 , t − τ ) = ϕ zipper ( θ i → θ i +1 , l i , t − τ ) R d l ′ d θ ′ ϕ zipper ( θ ′ → θ i +1 , l ′ , t − τ ) . (14)The probability distributions p origin ( t − τ | θ i +1 , t ) and p zip ( θ i , l i | θ i +1 , t − τ ) can be combined to determine theunzippering distribution p unzip ( θ i , l i | θ i +1 , t ) = Z t d τ p origin ( t − τ | θ i +1 , t ) × ϕ zipper ( θ i → θ i +1 , l i , t − τ ) R d l ′ d θ ′ ϕ zipper ( θ ′ → θ i +1 , l ′ , t − τ ) , (15)where we assume the system evolved from an initial con-dition at t = 0 in which no microtubules were present.Clearly all the complicated history dependence of the sys-tem is hidden in the originating time distribution. How-ever, in the steady state situation we consider below, thetime-dependence drops out and the details of this distri-bution become irrelevant. C. The steady state
We now consider the steady state of the system ofequations we have formulated. Setting the time deriva-tives to zero, the sum of equations (1a) to (1c) yieldsΦ growth + Φ shrinkage = 0, which together with equation(4) and (5) implies ∂ l i (cid:2) v + m + i ( l i , θ i ) − v − m − i ( l i , θ i ) (cid:3) =0. Because physically acceptable solutions should bebounded as l i → ∞ , we obtain the length flux balanceequation v + m + i ( l i , θ i ) = v − m − i ( l i , θ i ) , (16)showing that the growing and shrinking segments have,up to a constant amplitude, the same orientational andlength distribution. This allows us to eliminate m − i ( l i , θ i )from (1a) to obtain ∂ l i m + i ( l i , θ i ) = m + i ( l i , θ i ) × (cid:26) g − Z d θ ′ [ c ( | θ i − θ ′ | ) + z ( | θ i − θ ′ | )] k ( θ ′ ) (cid:27) , (17)where the growth parameter g = r r v − − r c v + . (18)characterizes the behaviour of the bare, non-interacting,system in which microtubules remain bounded in lengthfor g < g ≥ . As the brack-eted factor on the right hand side of (17) does not dependon the segment length nor on the segment number, weimmediately obtain that m + i ( l i , θ i ) has an exponentiallength distribution m + i ( l i , θ i ) = m + i ( θ i ) e − l i /l ( θ i ) (19) where the average segment length l ( θ i ) in the direction θ i is given by1 l ( θ ) = − g + Z d θ ′ ( c ( θ − θ ′ ) + z ( θ − θ ′ )) k ( θ ′ ) . (20)The nucleation boundary conditions (2) and (3) are nowtransformed into independent nucleation equations thatare expressed in terms of the amplitudes m + i ( θ i ) v + m +1 ( θ ) = r n π , (21) m + i ≥ ( θ ) = k ( θ ) Z d θ ′ z ( θ ′ − θ ) l ( θ ′ ) m + i − ( θ ′ ) . (22)We now note that under these conditions equation (1c)is already satisfied, as can be explicitly checked by con-sidering Φ reactivation (12) and using that in the steadystate ϕ zipper does not depend on time and the integralover p origin is by definition equal to 1. This, in combi-nation with the results (16), (19) and (22), yields theidentity with Φ zipper . We therefore need an independentargument to fix the densities of the inactive segments. Toobtain this we use the steady state rule that populationsize = nucleation rate × average lifetime . Consider a newly‘born’ growing segment, created either by a nucleation ora zippering event. Its average life time is by definitionthe average time until it shrinks back to zero length, i.e.the average return time. Clearly this time only dependson its orientation θ , the steady state microtubule lengthdensity k ( θ ′ ) and the dynamical instability parameters,but not on the segment number. We therefore denote itby τ ( θ ). The steady state density of inactive segmentswith length l i , orientation θ and segment number i is thengiven by m i ( l i , θ ) = Z d θ ′ ϕ zipper ( θ → θ ′ , l i ) τ ( θ ′ ) , (23)where ϕ zipper is defined by equation (13), as inactive seg-ments are created by a zippering event. Because the onlylength-dependent term on the right-hand side is m + ( l i , θ ),it follows that the length-dependence of the inactive seg-ment distributions is proportional to those of the activesegments, i.e. m i ( l i , θ ) = m i ( θ ) e − lil ( θ ) . (24)At the same time, the total integrated length density ofsegments, both active and inactive, with segment number i + 1 in the direction θ ′ is given by N total i +1 ( θ ′ ) = Z d l i +1 × (cid:2) m i +1 ( l i +1 , θ ′ ) + m + i +1 ( l i +1 , θ ′ ) + m − i +1 ( l i +1 , θ ′ ) (cid:3) = Z d θ ′′ Z d l i ϕ zipper ( θ ′′ → θ ′ , l i ) τ ( θ ′ ) , (25)where the last equality follows from the fact that everysegment with index i + 1 has been created by a zipper-ing event of a segment with index i . We solve this for τ ( θ ′ ) and insert the result in (23), which, after expand-ing ϕ zipper (13), produces the following expression for m i ( θ ): m i ( θ ) = m + i ( θ ) Z d θ ′ z ( θ − θ ′ ) l ( θ ′ ) × (cid:2) m i +1 ( θ ′ ) + m + i +1 ( θ ′ ) + m − i +1 ( θ ′ ) (cid:3)R d θ ′′ z ( θ ′′ − θ ′ ) l ( θ ′′ ) m + i ( θ ′′ ) . (26)We use the nucleation equation (22) to replace the in-tegral in the denominator of the integrand on the righthand side of this expression. In addition, we define thequantity Q i ( θ ) through m i ( θ ) = Q i ( θ ) (cid:2) m + i ( θ ) + m − i ( θ ) (cid:3) = (cid:18) v + v − (cid:19) Q i ( θ ) m + i ( θ ) ≡ uQ i ( θ ) m + i ( θ ) (27)and equation (26) leads to the following recursion relationfor Q i ( θ ) Q i ( θ ) = Z d θ ′ z ( θ − θ ′ ) k ( θ ′ ) l ( θ ′ ) (1 + Q i +1 ( θ ′ )) . (28)We now argue that the ratio Q i ( θ ) is in fact in dependentof the segment number. Using the fact that the growing,shrinking and inactive segments have an identical expo-nential profile, it follows from (27) that Q i ( θ ) is equal tothe ratio between inactive and active segments Q i ( θ ) = m i ( θ ) um + i ( θ ) = N i ( θ ) N + i ( θ ) + N − i ( θ ) . (29)After a new microtubule segment has been created itwill generally spend some time in an active state andsome time time in an inactive state. The expected life-time τ ( θ ) can also be separated into the expected ac-tive and inactive lifetimes for any newly created segment: τ ( θ ) = τ active ( θ ) + τ inactive ( θ ). These lifetimes are nec-essarily proportional to the total number of active andinactive segments, so that Q i ( θ ) = τ inactive ( θ ) /τ active ( θ ).As we have argued before, these lifetimes do not dependon the segment number, and, hence, neither does Q i ( θ ).An alternative route to the same conclusion follows fromexpanding out the forward recursion in (28) to show that Q i ( θ ) can for every i formally be written as the same in-finite series of multiple integrals involving z ( θ − θ ′ ), k ( θ )and l ( θ ). We therefore write the self-consistency relation-ship Q ( θ ) = Z d θ ′ z ( θ − θ ′ ) k ( θ ′ ) l ( θ ′ ) (1 + Q ( θ ′ )) . (30)The final closure of this set of equations is providedby the definition of the length density (8) applied to the steady state k ( θ ) = X i Z d l i l i (cid:2) m + i ( l i , θ ) + m − i ( l i , θ ) + m i ( l i , θ ) (cid:3) = ul ( θ ) (1 + Q ( θ )) X i m + i ( θ ) . (31) D. Dimensional analysis
In order to simplify our equations for further analysisand to identify the relevant control parameter we per-form a dimensional analysis. We therefore introduce acommon length scale and rescale all lengths with respectto this length scale. For example, our primary variables m + i ( θ ) have dimension [length] − [radian] − . Taking ourcue from (21) and (31) we adopt the length scale l = (cid:18) π v + u r n π (cid:19) , (32)where the additional factor of π − within the parenthesesis added to suppress explicit factors involving π in thefinal equations. This definition allows us to define thedimensionless variables L ( θ ) = l ( θ ) /l (33a) K ( θ ) = πk ( θ ) l (33b) M + i ( θ ) = πm + i ( θ ) l (33c) G = gl . (33d)In the absence of interactions, (20) shows that the av-erage length l of the microtubule is given by l = − /g .This implies G = − l /l , meaning that, for G < G can be interpreted as a measure for the non-interactingmicrotubule length.In addition, we adopt the dimensionless operator nota-tion F [ h ] ( θ ) = 1 π Z π d θ ′ f ( θ − θ ′ ) h ( θ ′ ) , (33e)where F ∈ { C, Z } . We are now in a position to ex-press the equations in terms of the dimensionless quan-tities. Applying (33c) to the nucleation equations (21)and (22) yields expressions for the growing segment den-sities M + i ( θ ). Furthermore, the M + i ( θ ) for the differentsegment labels can be absorbed into a single microtubuleplus end density (density of active segments), given by T ( θ ) = uL ( θ ) ∞ X i =1 M + i ( θ ) . (34)Performing all substitutions, the final set of dimension-less equations readsSegment length1 L ( θ ) = − G + C [ K ] ( θ ) + Z [ K ] ( θ ) (35a)Density K ( θ ) = L ( θ )(1 + Q ( θ )) T ( θ ) (35b)Inactive-active ratio Q ( θ ) = Z [ LK (1 + Q )] ( θ ) (35c)Plus end density T ( θ ) = L ( θ ) + L ( θ ) K ( θ ) Z [ T ] ( θ ) (35d)with G = (cid:20) v + v − r n ( v + + v − ) (cid:21) (cid:16) r r v − − r c v + (cid:17) (35e)Looking at the resulting equations, we see that thesegment length L is determined by the intrinsic growthdynamics ( G ) and the collisions leading to induced catas-trophes and zippering. The segment length density K is the product of the plus end density, the ratio of allsegments to active segments (1 + Q ) and the average seg-ment length. The ratio Q of inactive to active segmentsis modulated by the zippering operator, and the plus enddensity T consists of contributions from direct nucleationand zippered segments. We only consider parameter re-gions with physically realizable solutions that have havereal and positive values for L , K , Q and T .Finally, we note that the interaction operators definedby (33e) are convolutions of the operand with the interac-tion functions c ( θ ) and z ( θ ). Both interaction functionsare symmetric and π -periodic, and can therefore be writ-ten in terms of their Fourier coefficients as f ( θ ) = ˆ f ∞ X n =1 ˆ f n cos(2 nθ ) (36)ˆ f n = 1 π Z π d θ f ( θ )cos(2 nθ ) (37)Using the identity cos( θ − θ ′ ) = cos( θ )cos( θ ′ ) +sin( θ )sin( θ ′ ) we find that the functions cos(2 nθ ) andsin(2 nθ ) are eigenfunctions of the operators C and Z ,with the Fourier coefficients ˆ c n and ˆ z n , respectively, aseigenvalues: F [cos(2 nθ )] = ˆ f n cos(2 nθ ) (38)This convenient property will be exploited in later sec-tions. III. RESULTSA. Isotropic solution
In the isotropic phase all angular dependence dropsout. Because C [1] = ˆ c and Z [1] = ˆ z we are left with the set of equations1¯ L = − G + (ˆ c + ˆ z ) ¯ K (39a)¯ K = ¯ L (1 + ¯ Q ) ¯ T (39b)¯ Q = ˆ z ¯ L ¯ K (1 + ¯ Q ) (39c)¯ T = ¯ L + ˆ z ¯ L ¯ K ¯ T (39d)where the overbar denotes quantities evaluated in theisotropic phase. Solving for ¯ Q and ¯ T and inserting thisinto equation (39b) readily gives¯ K = ¯ L (1 − ˆ z ¯ L ¯ K ) , (40)which can be combined with equation (39a) to yield thefollowing relationship between G and the density¯ K (cid:0) ˆ c ¯ K − G (cid:1) = 1 . (41)We see that the isotropic density is an increasing functionof the microtubule dynamics parameter G and does notdepend on the amount of zippering. This can be under-stood by the fact that zippering only serves to reorientthe microtubules, which has no net effect in the isotropicstate. In the absence of induced catastrophes (ˆ c = 0),the density ¯ K diverges as G ↑
0, consistent with the re-sult by Dogterom and Leibler [22]. In the presence ofinduced catastrophes a stationary isotropic solution ex-ists for all values of G , although this solution need notactually be stable. B. Bifurcation analysis
We now search for a bifurcation point by consideringthe existence of steady-state solutions which are smallperturbations away from the isotropic solution. Thesesolutions are parametrized as follows L = ¯ L (1 + λ ) (42a) K = ¯ K (1 + κ ) (42b) Q = ¯ Q (1 + χ ) (42c) T = ¯ T (1 + τ ) (42d)Inserting these expressions into (35), subtracting theisotropic solutions and expanding to first order in theperturbations gives λ = − ¯ N ( C [ κ ] + Z [ κ ]) (43a) κ = λ + τ + ˆ z ¯ N χ (43b) χ = 1ˆ z Z (cid:2) λ + κ + ˆ z ¯ N χ (cid:3) (43c) τ = λ + ¯ N (ˆ z κ + Z [ τ ]) , (43d)where ¯ N = ¯ L ¯ K . Note that in these equations, ¯ N hasbecome the control parameter instead of G . Using (43b)and exploiting the linearity of Z , we expand Z [ κ ] = Z [ τ ] + Z (cid:2) λ + κ + ˆ z ¯ N χ (cid:3) − Z [ κ ] (44)= 1¯ N ( τ − λ ) − ˆ z κ + ˆ z χ − Z [ κ ] . (45)Solving this for Z [ κ ] and inserting the result into equa-tion (43a), combined with (43b), yields the relation(1 − ˆ z ¯ N ) κ = − N C [ κ ] (46)In the absence of induced catastrophes ( C [ κ ] = 0; onlyzippering), a bifurcation can only occur if ˆ z ¯ N = 1, whichonly happens for diverging density and G = 0 (as seenfrom equations (40-41)), therefore in this case there is nobifurcation. In the generic case where induced catastro-phes are present, (46) can be satisfied only if κ ( θ ) is aneigenfunction of C . We know that the family of functionscos(2 nθ ), n ≥
1, are eigenfunctions of C with eigenvaluesˆ c n , and therefore get a set of bifurcation values for ¯ N ,one for each eigenvalue: N ∗ n = ( − c n + ˆ z ) − . In addi-tion, we know that the isotropic solution must be stableas G → −∞ , because in this limit the microtubules havea vanishing length and do not interact. Therefore, therelevant bifurcation point is that for the lowest value of G , corresponding with the most negative eigenvalue of C (see also section III E). Assuming that the inducedcatastrophe probability increases monotonically with thecollision angle, ˆ c is always the most negative eigenvalue,so N ∗ = 1 − c + ˆ z . (47)We now derive the location of this bifurcation point interms of the control parameter G . Denoting ¯ N = ¯ L ¯ K ,equation (40) can be transformed to ¯ N (1 − ˆ z ¯ N ) = ¯ L ,into which we can substitute G ¯ L = (ˆ c + ˆ z ) ¯ N − G giving G ¯ N (1 − ˆ z ¯ N ) = (cid:2) (ˆ c + ˆ z ) ¯ N − (cid:3) (48)Combining this with the result (47) yields G ∗ = ( − c ) / (cid:18) ˆ c − c − (cid:19) . (49)The implication is that the location of the bifurcationpoint as a function of the control parameter G is deter-mined entirely by the eigenvalues of the induced catastro-phe function c ( θ ). Like the density in the isotropic phase,the location of the bifurcation point, this time perhapsmore surprisingly, does not depend on the presence oramount of zippering. C. Segment length and mesh size
An attractive interpretation of the microtubule lengthdensity K ( θ ) is that it represents the density of ‘obstacles’ that are pointing in the direction θ as seen by a micro-tubule growing in the perpendicular direction. From theobstacle density we can define a mesh size ξ ( θ ) - the aver-age distance between obstacles. Taking into account thegeometrical factor sin( θ ), we obtain ξ ( θ ) = (cid:20) πl Z π d θ ′ | sin( θ − θ ′ ) | K ( θ ′ ) (cid:21) − . (50)In the case of the isotropic solution, this simplifies to¯ ξ = πl / (4 ¯ K ). Using this equality we can derive anexpression for the average microtubule length ¯Λ in theisotropic phase, expressed in units of the mesh size. Thelength of each segment is given by ¯ L and the number ofsegments per microtubule is given by (1 + ¯ Q ), so using(40) we find ¯Λ = l ¯ L (1 + ¯ Q )¯ ξ = 4 ¯ K / π (51)Inserting this result into (41) provides the relationshipbetween ¯Λ and GG = (cid:18) π ¯Λ (cid:19) (cid:18) π ˆ c ¯Λ4 − (cid:19) (52)As was the case for the density, we see that the micro-tubule length as a function of mesh size does not dependon the amount of zippering. However, it should be notedthat the mesh size is defined through the average distancebetween single microtubules. In real systems, zipperingwould naturally lead to bundling, which in turn producesa system that has a larger mesh size between bundles (seealso the discussion).Combining equations (51) and (47), the expression for¯Λ at the bifurcation point becomes¯Λ ∗ = − π ˆ c (53)Assuming a monotonically increasing induced catastro-phe probability P c ( θ ), we know that the minimum valuefor ˆ c is reached when every collision at an angle largerthan 45 ◦ leads to a catastrophe. From (53), we see thatthis implies Λ ∗ ≥ / (2 √ l . In the absence of catastrophic collisions,we find ¯Λ | ˆ c =0 = 4 π ( − G − ) = 4 π (cid:18) ll (cid:19) , (54)where l = − /g is the average length of the microtubules. l is therefore a measure of the microtubule length thatis required to enable a significant number of interactions(¯Λ = 4 /π for l = l ). If the free microtubule length l is(much) shorter than l , the system is dominated by the(isotropic) nucleations, keeping the system in an isotropicstate. On the other hand, when l ≫ l , the interactionsdominate and, depending on the interaction functions,the system has the potential to align.0 c π P c θ c π P c θ c π P c θ ααα G α G α G α G α G α G α α K total α K total α K total acb S S S a b c
20 0.2 Figure 4: Bifurcation diagrams for the simplified interaction functions using three different induced catastrophe parameters.The figures on the left depict the probability P c ( θ ) to induce a catastrophe upon collision, along with the corresponding valuesof ˆ c . The centre and right columns depict the corresponding bifurcation diagrams as a function of G , expressed in terms ofthe total density K total and the 2D nematic order parameter S , respectively, where K total = R K ( θ )d θ . The isotropic solutionsare by definition disordered, so S = 0, and their density is computed from (41). The bifurcation point is determined using(49), with ˆ c = − /
2. For each diagram, ordered solutions have been computed for ˆ z = 0 (black), ˆ z = 1 (blue) and ˆ z = 10(red). The solutions have been computed using the method discussed in appendix A. Solid lines indicate stable solutions anddashed lines indicate unstable solutions (see also section III E). Note that the case of ˆ c = 1 in the absence of zippering is asingular case where the stability cannot be determined, because non-isotropic solutions only exist for G = 0. This has beenindicated by a dotted line. The S -diagrams include the asymptotic limit point at G = 0 with absolute ordering (at infinitedensity). The labels a , b and c indicate the parameter values of the results depicted in figure 5. The fact that the solutions for S in the case ˆ c = do not reach the asymptotic point ( G = 0 , S = 1) is a consequence of the slowdown in convergence ofthe path-following method as G ↑ D. Ordered solutions for simplified interactionfunctions
To find solutions beyond the immediate vicinity ofthe bifurcation point, we are hampered by the fact thatthese solutions are part of an infinite-dimensional solu-tion space. In appendix A it is shown that the solutionscan be constrained to a finite-dimensional space by re-stricting the interaction functions c ( θ ) and z ( θ ) to a finitenumber of Fourier modes.In this section, we will define a set of simplified interac-tion functions by restricting ourselves to Fourier modesup to and including cos(4 θ ). These modes provide uswith just enough freedom for the model to exhibit richbehaviour. Using the fact that c (0) = z (0) = z ( π/
2) = 0,we find that ˆ z = 0 and that both ˆ z and ˆ c are deter-mined by the remaining parameters. Furthermore, we introduce an overall factor of α in both equations, allow-ing us to set ˆ c = − /
2, so that c ( π/
2) = α . We thusobtain a system that is fully specified by the parametersˆ c , ˆ z and α . c ( θ ) = α (cid:20) ˆ c −
12 cos(2 θ ) + 12 (1 − ˆ c )cos(4 θ ) (cid:21) (55a) z ( θ ) = α (cid:20) ˆ z − cos(4 θ )) (cid:21) (55b)For α = 1, ˆ c and ˆ z are the actual Fourier coeffi-cients of the interaction functions. Demanding that P c ( θ ) = c ( θ ) / sin( θ ) is monotonically increasing on theinterval [0 , π/
2] leads to the constraint34 ≤ ˆ c ≤
98 (56)1and ˆ z is a positive real number. Of course, the totalprobability of zippering and catastrophe induction maynot exceed 1, placing an upper bound on α . In the ab-sence of zippering, we have α ≤ α has no quali-tative effect on the results. This can be understood byrealizing that the set of equations (35) is invariant underthe substitutions C → α C L → α − / L Z → α Z K → α − / KG → α / G T → α − / T where the first substitution reflects the presence of α inthe definitions (55). The existence of this scaling relationimplies that all functional dependencies between any ofthese parameters and variables remain unchanged whensubjected to the inverse scaling. Explicitly, the relevantparameters become C /α , Z /α and α − / G and the vari-ables α / L , α / K and α / T . With this in mind, wehave used the arbitrary choice α = 1 for our numericalcalculations, indicating the appropriate scaling on theaxes of figures 4 and 5.Equation (49) indicates that, for the simplified inter-action functions, the bifurcation point is located in therange − ≤ G ∗ ≤
18 (57)and from (41) and (53) we find that K ∗ = α − / andΛ ∗ = 4 / ( απ ). We have used the numerical proceduredescribed in appendix A to determine the ordered solu-tions of (35), starting from the bifurcation point. Thishas been done for nine different parameter values. Forthe values of ˆ c we used the extreme values 3 / / G ∗ = 0. For eachof these three cases, we have varied the zippering param-eter ˆ z , choosing values of 0, 1 and 10. Figure 4 showsthe results, depicting both the total density of the sys-tem and the degree of ordering as a function of G . Thedegree of ordering is measured by the order parameter S , defined as S [ K ( θ )] = | R π d θ e i θ K ( θ ) | R π d θ K ( θ ) , (58)which is the standard 2D nematic liquid crystal orderparameter (yielding 0 for a completely disordered systemand 1 for a fully oriented system). E. Stability of solutions
The bifurcation constraint (46) indicates that the spaceof bifurcating functions κ n is spanned by the functionscos(2 nθ ) and sin(2 nθ ) for a given value of n ≥
1. Thesesolutions are therefore symmetric with respect to an ar-bitrary axis that we choose to place at θ = 0. Even after θ θ θ θ θ θ z a)b)c) α T α T α T α K α K α K G α z α z α Figure 5: Three stable ordered solutions that correspond tothe points labelled (a), (b) and (c) in figure 4. The (unstable)isotropic solutions for the same parameter values are indicatedwith a dashed line. The parameter values for (a) and (b) differonly in the value of G , whereas the parameter values for (b)and (c) differ only in the value of ˆ z . All results have beencalculated using the method described in appendix A. the restriction to this symmetry axis there are still twosolution branches emanating from the bifurcation point,differing in the sign of the coefficient of the perturbation.These branches correspond to solutions peaked around θ = 0 and θ = π/ (2 n ), respectively, that are identical ex-cept for this rotation. The symmetry of these solutionsindicates that the bifurcation is of the pitchfork type. Infigure 5 we plot the solutions with a maximum at θ = 0.The presence of a pitchfork bifurcation implies a lossof stability of the originating branch [27]. In our sys-tem, we know that the isotropic solution must be stablein the limit G → −∞ . Therefore, the local stabilityof the isotropic solution is lost at the first bifurcationpoint (for the lowest value of G ), corresponding to theeigenfunction cos(2 θ ). Because this eigenfunction is or-thogonal to the eigenfunctions related to the subsequentbifurcation points (cos(2 nθ ) , n > G . This also means thatthe solution branches originating at further pitchfork bi-furcations will be unstable near the isotropic solution. Inthis paper we restrict ourselves to the analysis of the firstbifurcation point and the corresponding ordered solutionbranch. Because the solutions on this branch alreadyhave the lowest symmetry permitted by the interactionfunctions, there are no further bifurcation points alongthis branch.2Generically [27, chap. IV], the branches of the initialpitchfork bifurcation are stable for a supercritical bifurca-tion (branches bending towards higher values of G ) andunstable for a subcritical bifurcation (branches bendingtowards lower values of G ). In addition, turning pointsin the bifurcating branches generally correspond to an ex-change of stability [28, page 22]. This analysis allows usto assign stability indicators to the bifurcation diagramsin figure 4, even in the absence of a detailed study of thetime-dependent equations (1). IV. DISCUSSION AND CONCLUSIONS
Based on biological observations, we have constructeda model for the orientational alignment of cortical micro-tubules. The model has a number of striking features.First of all, for a given set of induced catastrophe andzippering probabilities ( P c ( θ ) and P z ( θ )), it allows usto identify a single dimensionless control parameter G ,which is fully determined by the nucleation rate and in-trinsic dynamics of individual microtubules. This resultby itself may turn out to be very useful in comparingdifferent in vivo systems or the same system under differ-ent conditions or in different developmental stages. Forincreasing values of G , the isotropic stationary solutionsto the model show an increase both in density and inabundance of interactions, as measured by the ratio ofmicrotubule length to the mesh-size. Secondly, the bifur-cation point, i.e. the critical value of G ∗ of the controlparameter at which the system develops ordered station-ary solutions from the isotropic state, is determined solelyby the probability of collisions between microtubules thatlead to an induced catastrophe.Indeed, from the numerical solutions of the minimalmodel introduced in Section III D it appears, perhapssurprisingly, that the co-alignment of microtubules dueto zippering events, if anything, diminishes the degreeof order. These results identify the “weeding out” ofmisaligned microtubules — by marking them for earlyremoval by the induced switch to the shrinking state —as the driving force for the ordering process. Finally,in spite of not being able to directly assess the stabilityof the solutions in the time-domain, we have providedarguments that stable ordered solutions are possible forthe regime G <
0, i.e. where the length of individualmicrotubules is intrinsically bounded.For
G >
0, individual microtubules have the tendencyto grow unbounded, unless they are kept in check bycatastrophic collisions. Although (locally) stable solu-tions may exist for values of G that are not too large, forevery G > G = 0, for which the microtubules are perfectlyaligned ( S = 1) and the system is infinitely dense. Theexistence of this point can be understood by the fact thatthe alignment also serves to decrease the number of colli- sions, and in the limit of a perfectly aligned system, the(relative) number of collisions vanishes.How realistic is the model presented? To answerthis question we need to consider several known factorsthat have not been included. First of all, microtubulestypically can de-attach from their nucleation sites andthen perform so called treadmilling motion, whereby theminus-end shrinks at a more or less steady pace, whichis small compared to both the growth- and the shrinkingspeed of the more active plus end. In the case that no zip-pering occurs at all it is relatively easy to show that theeffect of treadmilling simply leads to a renormalization ofthe parameter G and the interaction functions c ( θ ) and z ( θ ), but leaving the qualitative behaviour of the modelidentical to the one discussed here. When zippering doesoccur, one expects the treadmilling to enhance the degreeof ordering in an ordered state, as over time it “eats-up”the, by definition less ordered, initial segments of eachmicrotubule. This effect is also consistent with the ob-servation in figure 5c that in the case with zippering theactive tips are on average more strongly aligned than theaverage segment. In fact, given that the comparison be-tween figures 5b and 5c also shows that, all else beingequal, zippering sharpens the orientational distributionof the active tips as compared to the case with no zipper-ing, it is conceivable that the combination of zipperingand treadmilling could lead to more strongly ordered sys-tems for the same value of the control parameter.Next it is known that in vivo severing proteins, suchas katanin are active in, and crucial to, the formation ofthe cortical array [29]. Although in principle the effect ofsevering proteins could be included in the model, it wouldpresent formidable problems in the analysis as well asintroduce additional parameters into the model for whichprecise data is lacking.Another effect that has not been taken into accountexplicitly is microtubule bundling. Whenever a micro-tubule zippers alongside another segment, they form aparallel bundle [30]. However, the coarse-grained natureof our model precludes the formation of bundles and onlyallows for alignment of the segments. This means thata microtubule that is growing in a different direction en-counters each microtubule separately rather than as asingle bundle. It is to be expected that the catastropheand zippering rates stemming from N individual colli-sions will be higher than those from a single collisionwith a bundle of N microtubules. Hence, in realisticsystems the event rate is likely to be lower than that pre-dicted by the model, or, in other words, corresponds toa lower value of the scaling parameter α . Bundles arealso thought to be more than just adjacently aligned mi-crotubules, because they may be stabilized through as-sociation with bundling proteins that could potentiallydecrease the catastrophe rate of individual microtubuleswithin a bundle (see [31]). This is a non-trivial effect thatshould be considered separately and is likely to lead toan increased tendency to form an ordered structure.Finally, our model implicitly assumes that there is an3infinite supply of free tubulin dimers available for incor-poration into microtubules. Although there is no defi-nite experimental evidence for this, it is reasonable toassume that in vivo there is a limit to the size of thefree tubulin pool. Such a finite tubulin pool would havemarked consequences for the behaviour of the model,because the growth speed, and possibly also the nucle-ation rate, are dependent on the amount of free tubu-lin, or equivalently the total density of microtubules k tot .To a first approximation the growth speed is given by v + ( k tot ) = v + ( k tot = 0) (cid:16) − k tot k max (cid:17) where k max the max-imally attainable density when all tubulin is incorporatedinto microtubules. This allows for stable states to de-velop even when G ( k tot = 0) >
0, because under thispool-size constraint the length of individual microtubuleswill always remain bounded, and the system will settleinto a steady state with G ( k tot ) < G ( k tot = 0). Thisbehaviour could provide a biologically motivated mech-anism by which a solution with a particular density isselected.To see whether our model, in spite of its approximatenature, makes sense in the light of the available data wefirst use the collision event probabilities obtained by Dixitand Cyr [15] (see figure 3) to obtain an estimate for thebifurcation value of the control parameter of G ∗ = − . G > G ∗ . Given the available data on the micro-tubule instability parameters in this same system takenfrom Dhonukshe et al. [10] and Vos et al. [32] we wouldpredict using the definition (35e) that this requires thenucleation rate of new microtubules to be larger than 0.05min − µ m − (Dhonukse) and 0.01 min − µ m − (Vos) re-spectively. Both these estimates for a lower bound onthe nucleation rate are reasonable as they imply the nu-cleation of order 10 microtubules in the whole cortexover the course of the build-up towards full transverseorder, comparable to the number that is observed.Finally, we should point out that our model so faronly addresses the question of what causes cortical mi-crotubules to align with respect to each other. Giventhat in growing plant cells the cortical array is invariablyoriented transverse to the growth direction, the questionof what determines the direction of the alignment axiswith respect to the cell axes is as, if not more, importantfrom a biological perspective. We hope to address thisquestion, as well as the influence of some of the as yetneglected factors mentioned above, in future work. Acknowledgments
The authors thank Kostya Shundyak, Jan Vos andJelmer Lindeboom for helpful discussions. SHT is grate-ful to Jonathan Sherratt for his comments on the stabil-ity of solutions. RJH was supported by a grant withinthe EU Network of Excellence “Active Biomics” (Con- tract: NMP4-CT-2004-516989). SHT was supported bya grant from the NWO programme “Computational LifeScience”(Contract: CLS 635.100.003). This work is partof the research program of the “Stichting voor Funda-menteel Onderzoek der Materie (FOM)”, which is finan-cially supported by the “Nederlandse organisatie voorWetenschappelijk Onderzoek (NWO)”.
Appendix A: NUMERICAL EVALUATION OFTHE ORDERED SOLUTIONS
The solutions to the set of equations (35) lie in aninfinite-dimensional solution space. This creates signifi-cant hurdles for the numerical search for solutions. Inthis section, we will see that it is possible to restrictthe solutions to a finite-dimensional space by imposingconstraints on the interaction operators C and Z . Inaddition, we present a method to follow the branch of or-dered solution in this finite-dimensional space, startingfrom the bifurcation point (49).We start by reformulating the set of equations (35) byreplacing L ( θ ) and T ( θ ) through the definitions S ( θ ) = 1 L ( θ ) , U ( θ ) = 1 K ( θ ) (cid:18) T ( θ ) L ( θ ) − (cid:19) . (A1)Following these substitutions, the interaction operatorsare all applied at the outermost level of the equations,enabling us to make use of their properties in Fourierspace. Explicitly, we obtain S ( θ ) = − G + C [ K ] ( θ ) + Z [ K ] ( θ ) (A2) Q ( θ ) = Z [ K (1 + Q ) /S ] ( θ ) (A3) U ( θ ) = Z [(1 + KU ) /S ] ( θ ) (A4)and K ( θ ) = 1 + Q ( θ ) S ( θ ) − U ( θ )(1 + Q ( θ )) . (A5)Denoting the Fourier components of S ( θ ), Q ( θ ) and U ( θ ),by ˆ s n , ˆ q n and ˆ u n , respectively, the interacting micro-tubule equations reduce to a (potentially infinite) set ofscalar integral equations:ˆ s n = − δ n, G + ˆ c n + ˆ z n π Z π d θ cos(2 nθ ) K ( θ )(A6a)ˆ q n = ˆ z n π Z π d θ cos(2 nθ ) K ( θ ) (1 + Q ( θ )) S ( θ ) (A6b)ˆ u n = ˆ z n π Z π d θ cos(2 nθ )(1 + K ( θ ) U ( θ )) S ( θ ) . (A6c)From the structure of these equations, we immediatelysee that we can greatly reduce the dimensionality of theproblem by setting a number of Fourier coefficients ˆ z n c n to zero. In other words, by restricting our spaceof interaction functions c ( θ ) and z ( θ ), the problem canbe reduced to a finite number of scalar equations.Applied to the simplified interaction functions intro-duced in section III D, we know that the sets of station-ary solutions form lines in the 8-dimensional phase spacespanned by the variables { ˆ s , ˆ s , ˆ s , ˆ q , ˆ q , ˆ u , ˆ u } and theparameter G . At least two of such solution lines exist,one corresponding to the isotropic solution and the otherto the ordered solution, and these lines intersect at thebifurcation point.Within this 8-dimensional space, we have used a nu-merical path-following method similar to the one de-scribed in [33, 34] that follows the ordered solutionbranch by searching for a local minimum in the root meanerror of the constituent equations (A6). The search for or-dered solutions is initiated at the bifurcation point, with coordinates S ∗ = ˆ z − c ( − c ) / , Q ∗ = − ˆ z c , U ∗ = ˆ z ( − c ) / , (A7)so that in the case of our simplified interaction model { G, ˆ s , ˆ s , ˆ s , ˆ q , ˆ q , ˆ u , ˆ u } = { c − , z + 1) , , , z , , z , } . (A8)The initial instability affects only the cos(2 θ ) mode. Thismode only appears in the equation for ˆ s and the remain-ing parameters are affected only by higher order correc-tions. For this reason we choose the initial direction ofthe path to be the unit vector in the ˆ s -direction and thepath is traced from there. [1] Bruce Alberts, Alexander Johnson, Julian Lewis, MartinRaff, Keith Roberts, and Peter Walter. Molecular Biologyof the cell . Garland Science, 4th edition, 2002.[2] D.W. Ehrhardt and S.L. Shaw. Microtubule dynamicsand organization in the plant cortical array.
Annu. Rev.Plant Biol. , 57:859–75, 2006.[3] A. Paradez, A. Wright, and D.W. Ehrhardt. Microtubulecortical array organization and plant cell morphogenesis.
Curr. Opin. Plant Biol. , 9(6):571–578, 2006.[4] D.J. Cosgrove. Growth of the plant cell wall.
Nat. Rev.Mol. Cell Biol. , 6:850–861, 2005.[5] G.O. Wasteneys and R.E. Williamson. Reassembly ofmicrotubules in nitella tasmanica - quantitative-analysisof assembly and orientation.
Eur. J. Cell Biol. , 50(1):76–83, 1989.[6] F. Kumagai, A. Yoneda, T. Tomida, T. Sano, T. Nagata,and S. Hasezawa. Fate of nascent microtubules organizedat the M/G1 interface, as visualized by synchronized to-bacco BY-2 cells stably expressing GFP-tubulin: time-sequence observations of the reorganization of corticalmicrotubules in living plant cells.
Plant Cell Physiol. ,42(7):723–732, 2001.[7] S.L. Shaw, R. Kamyar, and D.W. Ehrhardt. Sustainedmicrotubule treadmilling in Arabidopsis cortical arrays.
Science , 300(5626):1715–1718, June 13 2003.[8] J.W. Vos, B. Sieberer, A.C.J. Timmers, and A.M.C.Emons. Microtubule dynamics during preprophase bandformation and the role of endoplasmic microtubules dur-ing root hair elongation.
Cell Biol. Int. , 27(3):295, 2003.[9] A.R. Hardham and B.E. Gunning. Structure of corticalmicrotubule arrays in plant cells.
J. Cell Biol. , 77(1):14–34, 1978.[10] P. Dhonukshe, A.M. Laxalt, J. Goedhart, T.W.J.Gadella, and T. Munnik. Phospholipase d activation cor-relates with microtubule reorganization in living plantcells.
Plant Cell , 15(11):2666–2679, 2003.[11] J.C. Gardiner, J.D. Harper, N.D. Weerakoon, D.A.Collings, S. Ritchie, S. Gilroy, R.J. Cyr, and J. Marc.A 90-kd phospholipase d from tobacco binds to mi-crotubules and the plasma membrane.
Plant Cell ,13(9):2143–2158, 2001. [12] T. Hashimoto and T. Kato. Cortical control of plantmicrotubules.
Curr. Opin. Plant Biol. , 9(1):5–11, 2006.[13] T. Hamada. Microtubule-associated proteins in higherplants.
J. Plant Res. , 120(1):79–98, 2007.[14] V. Kirik, U. Herrmann, C. Parupalli, J.C. Sedbrook,D.W. Ehrhardt, and M. H¨ulskamp. Clasp localizes in twodiscrete patterns on cortical microtubules and is requiredfor cell morphogenesis and cell division in arabidopsis.
J.Cell Sci. , 120(24):4416–4425, 2007.[15] R. Dixit and R. Cyr. Encounters between dynamic cor-tical microtubules promote ordering of the cortical arraythrough angle-dependent modifications of microtubulebehavior.
Plant Cell , 16(12):3274–3284, 2004.[16] E. Geigant, K. Ladizhansky, and A. Mogilner. An integro-differential model for orientational distributions of f-actinin cells.
SIAM J. Appl. Math , 59:787–809, 1998.[17] K. Kruse, J.F. Joanny, F. J¨ulicher, J. Prost, and K. Seki-moto. Generic theory of active polar gels: a paradigm forcytoskeletal dynamics.
Eur. Phys. J. E , 16:5–16, 2005.[18] I.S. Aranson and L.S. Tsimring. Theory of self-assemblyof microtubules and motors.
Phys. Rev. E , 74(3):031915,2006.[19] V. R¨uhle, F. Ziebert, R. Peter, and W. Zimmermann. In-stabilities in a two-dimensional polar-filament-motor sys-tem.
Eur. Phys. J. E , 27:243–251, 2008.[20] A. Zumdieck, M. Cosentino Lagomarsino, C. Tanase,K. Kruse, B.M. Mulder, M. Dogterom, and F. J¨ulicher.Continuum description of the cytoskeleton: ring forma-tion in the cell cortex.
Phys. Rev. Lett. , 95(25):258103,2005.[21] V.A. Baulin, C.M. Marques, and F. Thalmann. Collisioninduced spatial organization of microtubules.
Biophys.Chemist. , 128:231–244, 2007.[22] M. Dogterom and S. Leibler. Physical aspects of thegrowth and regulation of microtubule structures.
Phys.Rev. Lett. , 70(9):1347–1350, 1993.[23] T. Murata, S. Sonobe, T.I. Baskin, S. Hyodo,S. Hasezawa, T. Nagata, T. Horio, and M. Hasebe.Microtubule-dependent microtubule nucleation based onrecruitment of gamma-tubulin in higher plants.
Nat. CellBiol. , 7(10):961–968, 2005. [24] T. Murata and M. Hasebe. Microtubule-dependent micro-tubule nucleation in plant cells. J. Plant Res. , 120(1):73–78, 2007.[25] R.L. Margolis and L. Wilson. Microtubule treadmilling:what goes around comes around.
BioEssays , 20:830 – 836,1998.[26] M.E. Janson, M.E. de Dood, and M. Dogterom. Dynamicinstability of microtubules is regulated by force.
J. CellBiol. , 161(6):1029–1034, 2003.[27] M. Golubitsky and D.G. Schaeffer.
Singularities andGroups in Bifurcation Theory , volume 1. Springer-Verlag,1984.[28] G. Iooss and D.D. Joseph.
Elementary Stability and Bi-furcation Theory . Springer-Verlag, 1980.[29] A. Roll-Mecak and R.D. Vale. Making more microtubulesby severing: a common theme of noncentrosomal micro-tubule arrays.
J. Cell Biol. , 175:849–851, 2006.[30] D.A. Barton, M. Vantard, and R.L. Overall. Analysis ofcortical arrays from Tradescantia virginiana at high res-olution reveals discrete microtubule subpopulations anddemonstrates that confocal images of arrays can be mis-leading.
Plant Cell , 20:982–994, 2008.[31] J. Gaillard, E. Neumann, D. van Damme, V. Stoppin-Mellet, C. Ebel, E. Barbier, D. Geelen, and M. Van-tard. Two Microtubule-associated Proteins of Arabidop-sis MAP65s Promote Antiparallel Microtubule Bundling.
Mol. Biol. Cell , 19(10):4534–4544, 2008.[32] J.W. Vos, M. Dogterom, and A.M.C. Emons. Micro-tubules become more dynamic but not shorter duringpreprophase band formation: a possible ”search-and-capture” mechanism for microtubule translocation.
CellMotil. Cytoskeleton , 57(4):246–258, 2004.[33] E.L. Allgower and K. Georg.
Introduction to NumericalContinuation Methods , volume 45 of
Classics in AppliedMathematics . SIAM, 2003.[34] P. Deuflard, B. Fiedler, and P. Kunkel. Efficient numeri-cal pathfollowing beyond critical points.