The self-organization of plant microtubules in three dimensions enable stable cortical localization and sensitivity to external cues

Many cell functions rely on the ability of microtubules to self-organize as complex networks. In plants, cortical microtubules are essential to determine cell shape as they guide the deposition of cellulose microfibrils, and thus control mechanical anisotropy in the cell wall. Here we analyze how, in turn, cell shape may influence microtubule behavior. Using a computational model of microtubules enclosed in a three-dimensional space, We show that the microtubule network has spontaneous configurations that could explain many experimental observations without resorting to specific regulation. In particular, we find that the preferred localization of microtubules at the cortex emerges from directional persistence of the microtubules, combined with their growth mode. We identified microtubule parameters that seem relatively insensitive to cell shape, such as length or number. In contrast, microtubule array anisotropy depends strongly on local curvature of the cell surface and global orientation follows robustly the longest axis of the cell. Lastly, we found that the network is capable of reorienting toward weak external directional cues. Altogether our simulations show that the microtubule network is a good transducer of weak external polarity, while at the same time, it easily reaches stable global configurations. Author summary Plants exhibit an astonishing diversity in architecture and shape. A key to such diversity is the ability of their cells to coordinate and grow to reach a broad spectrum of sizes and shapes. Cell growth in plants is guided by the microtubule cytoskeleton. Here, we seek to understand how microtubules self-organize close to the cell surface. We build upon previous two-dimensional models and we consider microtubules as lines growing in three dimensions, accounting for interactions between microtubules or between microtubules and the cell surface. We show that microtubule arrays are able to adapt to various cell shapes and to reorient in response to factors such as signals or environment. Altogether, our results help to understand how the microtubule cytoskeleton contributes to the diversity of plant shapes and to how these shapes adapt to environment.


Introduction
Despite their amazing diversity in shapes, biological organisms share some common structural components at the cellular level. Among those, one of the best conserved proteins across eukaryotes, tubulin, assembles into protofilaments, which in turn form 25 nm nanotubes known as microtubules, usually made of 13 protofilaments. The network of microtubule is highly labile and can reshape itself in a matter of minutes. In plants, microtubules form superstructures before (the preprophase band), during (the spindle) and after (the phragmoplast) cell division. Plant microtubules also form dense and organized arrays at the periphery of the cell during interphase [1] and these arrays are known as cortical microtubules (CMTs). The behavior of CMTs has been studied extensively at the molecular level [2]. One of their main functions is to guide the trajectory of the transmembrane cellulose synthase complex and thus to bias the orientation of cellulose microfibrils in the wall. This in turn impacts the mechanical anisotropy of the cell wall and controls growth direction [3], [4] , [5] , [6]. This function explains why most mutants impaired in microtubule-associated proteins exhibit strong phenotypic defects [7].
Whereas this provides a clear picture on how microtubules impact cell shape, in turn how cell shape impacts microtubule behavior has been much less explored. There is evidence that shape-derived mechanical stress can bias microtubule orientation towards the direction of maximal tension, both at the tissue and single cell scales [8] , [9]. The molecular mechanism behind this response remains unknown. In addition, cell geometry may also affect microtubule behavior, independently of the presence of mechanical stress in cell walls. This is what we explore here. Before doing so, we first review the main features of the microtubule network.
The molecular basis for microtubule dynamics is rather well established. Consistent with the absence of centrosome in land plants, microtubule nucleation is dispersed in plant cells, as it occurs at the cell cortex [10], along existing microtubules during branching events [11], [10], and at the nuclear envelope [12]. As they grow, microtubules form stiff and polar structures. They can alternate growth, pause and shrinking at the so-called plus end [13], whereas they mainly shrink or pause at the minus end [14]. The combination of an average shrinkage at the minus end and dynamic instability at the plus end leads to an overall 3 50 55 60 65 displacement of the microtubule, also called hybrid treadmilling. This is the main cause for microtubule encounters, and thus for reshaping the microtubule network [13], [14].
When one microtubule encounters another microtubule, different outcomes can be observed: if the encounter angle is shallow, zippering can occur, i.e. the growing microtubule bends and continues its polymerization along the encountered microtubule, which leads to the creation of microtubule bundles; if the encounter angle is steep, crossover can occur, i.e. a microtubule polymerizing without deviating its trajectory and crossing over the encountered microtubule; or alternatively catastrophe is triggered, i.e. a rapid plus end shrinkage after contact with the encountered microtubule.
Such selective pruning of microtubules may explain how microtubules can form parallel arrays from initially random orientations, and conversely change the net orientation of their arrays over time, through a phase of randomisation [15], [16]. The presence or absence of different microtubules associated proteins (MAPs) can modulate the stability of microtubules or their capacity to form bundles and to self-organize. For instance, the microtubule severing protein Katanin accounts for most of the pruning events at crossover sites [17].
The microtubule network is a typical example of a self-organizing system, where properties of individual elements and their interactions induce specific and sometimes counter-intuitive global properties. To predict how regulation at the level of each microtubule can give rise to specific global outcomes, one can resort to computational models. Modeling approaches have been developed, simplifying microtubule interactions by restricting them to the plasma membrane, i.e. a simpler 2D space [18], [19]. In those agent-based models, several microtubule properties were coded and interactions between CMTs, based on these properties, were simulated. The outcome is an emergent network, whose characteristics can be analyzed. For instance, increased microtubule severing was predicted to generate a larger number of free microtubules, more amenable to bundle into aligned arrays [20] and this was observed in experiments [21].
So far, most of the microtubule models have been implemented in a 2D space. A major outcome of such models was to demonstrate that global orientations of the network can spontaneously emerge from the interactions between microtubules. Many combinations of parameters and behaviors have been studied: 4 75 80 85 90 95 instability at the plus end [22] , [23], role of zippering [20], [24], [25], [26], nucleation modes [26], [27], and severing [28]. Beyond the differences, the fact that a global orientation emerges in many of these combinations suggests that reaching one global orientation is a robust feature of microtubule networks.
Conversely, the diverse combinations of microtubule properties provide different scenarios for the finetuning of the network structure and stability of this emerging behavior.
Some aspects of cell geometry were related to microtubule behavior in certain simulations. Simulations showed how different directional biases in nucleation can induce an ordering of the array toward directions that are correlated to cell geometry [22] and [29]. Further, branching nucleation rules can elicit handedness of the global direction of microtubule arrays, provided that the branching is biased toward one direction [30].
Other studies used a simulation space where borders, analog to cell edges, induce more or fewer catastrophe events or are more or less permissive toward microtubule growth [31], [22], [30] and [15]. Most studies concluded that a global orientation of microtubules can be correlated to cell face orientations.
The contribution of the third dimension to microtubule behavior has started to be investigated.
Computational models for animal systems have focused on 3D considerations but the nucleation hypotheses are too different from that in plants to be transposed directly [32], [33]. Accordingly, fully 3D models suited for plants are still lacking: almost all existing studies have only considered microtubules living on surfaces embedded in 3D [34], [30], [31], [35]. In [22], a 2D model was extended into a full 3D model but it did not include cell boundaries, which yielded microtubules distributed over the whole simulated domain, in contrast with the cortical localization of microtubules in planta.
In this paper, we explore the influence of 3D cell shape on the basic properties of a dynamic microtubule network. We do not assume that microtubules live on a 2D surface; rather, we simulate a closed volume where microtubules are more or less free to grow in all directions. Anchoring to the membrane is not imposed by the model and instead becomes a variable in the model. Using this framework, we investigate how microtubule interaction with the membrane influences microtubule dynamics. Our study also addresses the relative contributions of cell shape, microtubule interactions, and external directional cues in network

Microtubules become cortical in a 3D space because of their directional persistence and growth mode
We modelled microtubules as a set of line segments that nucleate, grow, shrink, and interact between them and with the cell surface represented as a triangular mesh (Methods). Nucleation occurs randomly at the surface. Growth occurs from the plus end with a small directional noise that is related to the persistence length of microtubules. Shrinkage occurs at the minus end. A microtubule that encounters a previously existing microtubule either changes direction to that of the previous microtubule (if the encounter angle is steep) or undergoes catastrophe (complete shrinkage). We considered two types of interaction with the cell surface: strong anchoring, whereby a microtubule that reaches the surface continues growing tangentially to the surface, and weak anchoring, whereby microtubules are prevented from leaving the cell interior. Typical   [6] and CSI1 was also proposed to stabilize the microtubule network [37].
Yet, the influence of CMT-CESA interactions on the microtubule network is still poorly understood. Thus, we took advantage of the 3D nature of our model to study the impact of the anchoring rule to the membrane In case of the strong anchoring, microtubule growth is biased so as to stay near the membrane, as if anchoring proteins were highly concentrated. As expected, strong interactions led to a surface-localized cortical zone with microtubules trajectories embedded in the plane parallel to the mesh (Fig 2A). We then tested the microtubule response when anchoring to the membrane is weak. In that case, microtubules can crawl in every direction, but sometimes as they encounter the membrane, their direction may be transiently tangent to the membrane. Even if weak anchoring allows microtubules to travel through all the volume of the cell, we find that such weak interaction with the membrane is enough to elicit the existence of cortical microtubules. In those conditions, the predicted distance between microtubules and membrane would be from 50 to 250 nm (Fig 2A). Therefore, the three-dimensional nature of our model helps us demonstrate that strong anchoring is not required for the presence of large populations of cortical microtubules in plant cells: the directional persistence, together microtubule growth mode, can cause such subcellular localization.
As the microtubules tend to stay at close to the membrane, they also bundle, at proportions varying from 30% to 75% ( Fig 2B). The ability for the microtubule network to generate a spontaneous bundled structure is consistent with 2D models. Strikingly, this effect is also present in the case of weak anchoring.

Strong anchoring decreases microtubule number and length, and increases microtubule bundling
Independently of the encounter rule, weak anchoring strength increases the total number and the size of 7 160 165 170 175 microtubules (by about 20%, in length and in number) when compared to strong anchoring. Weak anchoring strength also generated less bundling (reduction of about 20%). This effect could be due to the fact that a weak anchoring to the membrane allows microtubules to crawl inside the cell, thus diminishing the encounter probability. Consequently, microtubules weakly bound to the membrane have more space to grow and are less subject to blockage or bundling by other microtubules (Fig 2B). Lastly, we found that anchoring strength significantly affected the anisotropy of the microtubule arrays. High anchoring increased the anisotropy of the network by about 10% (Fig 4C).
The number and length of microtubules, as well as the proportion of microtubules in bundles, are relatively independent of cell shape Next, we used our model to determine the consequences of changing the cell shape on global properties of the network. We simulated the network in three main sharp shapes represented on Fig 1A. The results from the simulations indicate that these shapes only have a marginal effect on the number of microtubules (less than 10% difference), on the length of microtubules (5 to 15% difference) and on the proportion of bundles (less than 15% difference). Overall, elongated cells have more microtubules, more bundles, and longer microtubules, while cubes show the lower values (Fig 3).

Microtubule array anisotropy is not influenced by the cell global shape but by its curvature
We also analyzed the effect of cell shape on the microtubule array anisotropy, averaged over the cortex (see Methods). As microtubule array anisotropy is skewed towards low values (found inside the cell) we used a non-parametric test based on ranks for comparisons.
The type of cell shape did not appear to influence the global anisotropy of the microtubule network. This is interesting as it suggests that different cell types with various shapes do not require differential regulation of  However, local curvature significantly increased the anisotropy of the microtubule network, by 10% (p<0.001) (Fig 4B). Low convex curvatures that characterize flat cell faces generate less anisotropy in the arrays than higher curvatures. This effect was reduced in the case of elongated cells.

The average microtubule orientation is strongly influenced by cell shape
In order to determine how the shape of the cell influences the global orientation of the network, we measured First, we observed that in case of square shaped faces, most of the microtubules align along the cell face diagonal, i.e. the longest path (Fig 5). This occurred whatever the anchoring strength and the shape of faces at the side. Second, we observed a strong correlation between the long axis of the cell and microtubule orientation, and this correlation was higher for the elongated cell. Indeed, the microtubule distribution was centered around an angle of 0. This effect is higher in the case of strong anchoring, whereas in the case of weak anchoring, secondary peaks show that the diagonals are also overrepresented. These simulation results indicate that the microtubule network is able to read the longest axis of the cell and orient toward that axis, by default. Interestingly, cortical microtubules become longitudinal in hypocotyl cells when growth stops [5], [38], [44], ], suggesting that they may adopt their configuration "by default" in that situation. Arrow showing the reference for angle measurement for the elongated shape (top) and for the square or the cubic shape (bottom). The shapes are not at the same scale.

Small external directional bias is sufficient to affect the orientation of the whole network
Next, we investigated the robustness of microtubule network behavior. Microtubule arrays entirely reorganize during cell division [39]. Light and hormones can also completely reorient the microtubule network within minutes [40], [15], [16], suggesting that the constraints on microtubules cannot be too strong to allow such rapid reorganization in vivo. We tested whether our model provides such adaptability, using the case of external, directional cues (Fig 6). We investigated the effect of an external cue, assuming that, as microtubules grow in the cell, their direction is biased toward the cue, with a specific weight. We used an ellipsoidal cell in which the membrane exhibits a directional circumferential cue perpendicular to the main axis of the cell. In the absence of such cue, microtubules are on average parallel to the main axis of the cell. Strikingly, our simulations indicate that even a very low bias (~0.1%) could disrupt the main orientation of the network. When the weight reaches 1% or more, microtubules massively reorient toward the direction of the cue.
These results indicate that, despite an apparent robust organization, the microtubule network remains 10 240 245 250 255 260 extremely sensitive to directional cues. As such, it is capable of reading slight directional cues and generating a strong polarity toward that direction. This ability to read directional cues is probably linked to the self-enhancement of microtubule orientation through their interactions: As more microtubules orient toward a direction, they prevent growth of microtubules in the perpendicular direction.
Our simulations indicate that the network should exhibit two behaviors concerning its polarity. When no external directional cue is present, the network orients toward the main axis of the cell, and generates a polarity that is a direct reading of the global shape. When a directional cue is present, the network reorients so as to emphasize the direction of the cue.

The default state of the microtubule network
The impact of external cues on the microtubule network is well characterized. However, because real cells and tissues are never really devoid of external cues, the behavior of the microtubule network by default in a plant cell remains an open question. Our work provides some clues to address this question, taking into account the 3D shape of the cell and with the most minimal set of parameters for microtubule behavior. We found that microtubule directional persistence largely determines the subcellular localization and orientation of the microtubule network in various cell shapes. We also identified parameters that seemed relatively insensitive to cell shape, such as microtubule length or number. Last, despite the constraints of microtubule stiffness and growth mode on the organization of the microtubule network, we found that such parameters allow high sensitivity to directional cues, even when such cues go against the default orientation. Altogether, this provides a conceptual framework to dissect the exact contribution of microtubule regulators to the microtubule network organization, in relation to 3D cell shape.

Cortical localization is an emerging property of microtubule stiffness and growth mode in convex cells
Based on TEM images where microtubules are often seen very close to the plasma membrane, it is assumed that anchoring of microtubules to the plasma membrane is relatively strong. However how this anchoring 11 270 275 280 285 290 would be mediated remains unknown. Many biochemical studies have been performed to extract proteins that would link the plasma membrane to cortical microtubules, and so far the only published candidate is a phospholipase [41], for which no follow-up results have been obtained to our knowledge. Other links have been put forward, such as CLASP [31] or CELLULOSE SYNTHASE INTERACTIVE PROTEIN 1 (CSI1) [36], [42], but they might represent rather indirect regulators of the link between microtubules and plasma membrane. Does this mean that microtubule could be cortical without any anchoring module? Our results suggest that in most of the measured cell shapes, microtubules do not need to be strongly connected to the membrane to remain cortical. This prediction implies that modulating anchoring would affect self-organizing properties tangentially to the cell surface, rather than modulating the density of microtubules inside the cell.
This would allow molecular regulators to modify the microtubule organization without directly affecting the rate of cellulose deposition.

A high curvature increases anisotropy
In our simulations, we observed differences in anisotropy originating from the sharpness of face-face contacts and curvature of the cell. Cells that are more curved exhibit increased array anisotropy compared to flatter ones. This result is interesting, as epidermal cells possess faces with very different curvatures. In the L1 layer of the shoot apical meristem for example, the outermost wall has a stronger curvature than all the walls separating the cell from its neighbors. Pavement cells also display various local curvature values.
Based on our results, higher anisotropy would be produced in L1 cells without specific regulation. Such results are difficult to account for in models with 2D surfaces embedded in 3D [31].
Similarly, changes in the curvature of the epidermal wall can occur through changes in internal pressure, which in return influence the anisotropy of the network [43]. An increase in the microtubule organization could be a result of an increased pressure, that would increase the curvature of the epidermal cell. Further experimental work is required to investigate the role of curvature on microtubule behavior and its relation to mechanical stress in the epidermis. Our simulations show that the cell aspect ratio has an impact on the global orientation of the network. The predicted default behavior of microtubules is their alignment parallel to the long axis of a cell, due to the directional persistence of microtubules associated with their bending stiffness. This is in agreement with previous models with microtubules living on 2D surfaces embedded in 3D [34], [35], where the default orientation is longitudinal for long cylinders.
This default state was observed with microtubules polymerizing in vitro inside elongated chambers [34]. In slowly growing cells of the hypocotyl, microtubules are oriented along the long axis of the cell, whereas microtubules are circumferential in rapidly elongating cells [38], [44]. Our model suggests that directional cues are needed to avoid this default orientation is growing cells.
At the boundary between the shoot apical meristem and the primordia, cell division leads to cell shapes that are elongated along the axis of the boundary; our model predicts that microtubules will be oriented along the same direction by default, amplifying their response to mechanical stress [8].

A weak directional bias is sufficient to change the orientation of the microtubule network
In this study, we show that the microtubule network is oriented by default along the longest axis of the cell.
However, microtubules in plants often show supracellular orientation, independently of cell shape, a behavior that has been ascribed to tissue-level signals, notably mechanical stress [8]. Moreover, It has been demonstrated that inside the cell, microtubules orientation relates to polarity markers such as proteins from the PIN FORMED and RHO OF PLANTS families [45], [46]. Simulations have assessed how localized membrane heterogeneity could result in a biased orientation of the microtubule network [31], [30]. In this study, we show that a weak directional cue influencing microtubule growth rapidly modifies the orientation of the network towards the direction of the cue. As such, the network behaves as a perception tool translating  [16].

Perspectives
The shape of the cell has little influence on mean length, number of microtubules or bundles proportions. In addition, anisotropy of the network is not highly correlated to changes in global cell shape. This prediction of a robust network suggests that plant cells do not need specific regulations to compensate for their great variations in shapes. Accordingly, the microtubule network appears as a good polarity system, with a default orientation and a high sensitivity to directional cues. It was recently shown that, for global polarity to emerge in a tissue, the only requirement is the existence of internal cellular polarity [ 47]. In this work we show that the microtubule network is adapted to such requirement.
Overall, microtubules and associated proteins form a complex self-organizing system that is difficult to comprehend without resorting to models. The results obtained here demonstrate that our three-dimensional model provides a framework to test hypotheses on the regulation of the microtubule cytoskeleton in plant cells.
The model given here is a beginning to a more complete analysis. We have not yet incorporated microtubule severing [40], [28], microtubule branching [11], and the possible effects of connections between cortical microtubules and cellulose fibrils outside of the cell as mediated by the cellulose synthase complex [48].
Severing in particular has been shown to be key to microtubule reorientation [28] following mechanical signals [49]. We also have not included limiting levels of tubulin [35], which could affect overall microtubule number. Altogether, we expect our model to help progress in understanding how microtubule self-organization integrates three-dimensional cell shape and directional cues and how microtubuleassociated proteins modulate this integration. 14

Methods
The dynamical microtubule model was implemented in C++. Simulations were performed on Intel/AMD desktop computers running Debian and Ubuntu operating systems.

The microtubule network
Microtubules were coded as 3D multi-segment vectors of constant length. A ring of tubulin of the height of a dimer is represented as a unit vector in the simulation. Microtubule growth in the model occurs by adding one vector element at the plus end of the microtubule, at a position in the direction pointed to by the last vector. In plants, microtubules are considered static, their growth and shrinkage are the result of treadmilling processes. To code for microtubule directional persistence (which relates physically to bending stiffness), the direction at which a new vector is added to an existing microtubule changes by a small random amount.
Microtubule shrinkage occurs at the minus end by removing the first vector from the list.

Time and space correspondences
Typically, a cell has a width of several micrometers, and we take the unit of length as 8nm, the height of a

Nucleation and minus end rules
The microtubules are nucleated close to the cell surface at a rate of 1 to 20 microtubules/s. Once nucleated, the microtubules do not immediately shrink. At each of the time steps that follow nucleation, a microtubule has a probability n s to begin shrinkage.

Plus-end growth
Without any interactions with the membrane or with neighbouring microtubules, the new direction of growth is that of the previous vector (tubulin), modified by a random direction r d , typically deviating less than 10° 15 375 380 385 390 395 from the previous direction. This corresponds to a persistence length of a few micrometers. In the presence of directional cue, the new vector is computed as the weighted sum of its unbiased direction and the direction of the cue vector. When a microtubule reaches the plasma membrane, it follows one of two behaviors depending on the interaction strength with the membrane. a) strong anchoring: after interaction, the microtubule grows along to the surface with further new vectors added tangentially to the surface, or b) weak anchoring: the microtubule grows as in the case of free space, except that it cannot go out of the cell: if the new vector is outside the cell then the vector is projected on the plane tangential to the surface of the cell; due to the noise on vector orientation, the microtubule may leave cell surface.
When a plus end is closer than 25 nm to another microtubule, two situations can occur: a) if the angle between the microtubule and the encountered one is smaller than a threshold α, the microtubule begins zippering, and its direction becomes that of the previously existing microtubule; b) if the angle is larger than α, the microtubule begins to shrink at the plus end (catastrophe event) or stalls according to the encounter rule. α was set to 40° [20]

Cell shape and interaction with the membrane
The cell contour is described with a triangular mesh. Each vertex possesses the information of the vector normal to the surface, which is used during the simulation to calculate a local approximation of the tangential plane. It is possible to add information at the vertex level that can be read during the simulation and serve as extrinsic input. Inputs can be scalars or tensors. The distance to the membrane is calculated using the nearest point on the surface. At this point, the membrane is approximated by the plane perpendicular to the normal of the mesh. The distance between a tubulin ring (unit vector) and the membrane is calculated as the shortest 16 405 410 415 420 distance between the endpoint of the vector and this plane. A collection of standard cell shapes was generated for our simulations (Fig. 1A).

The simulation process
The simulation progress is made through discrete timesteps. At each timestep new vectors are added and vectors are removed from the simulation space. Collision tests are realized so as to implement the different growth or shrinkage situations. In order to increase the simulation speed, space is divided into subelements and the vectors are identified according to the subelement to which they belong, which diminishes the number of particles involved in the collision test (locality-sensitive hashing [51].

Visualization
To visualize microtubule density and orientations an image is created using a matrix of resolution rx,ry,rz. rz is larger relative to rx and ry (typically 10 to 1 ratio), which mimics the anisotropic quality of confocal microscope resolution. The simulation space is then screened. When a vector is located inside a cube of the matrix, the value of this cube is incremented by one, and the immediate neighbouring cubes are incremented by a lower number (typically 0.3). At the end of the process, a stack is formed where microtubules appear as blurred intensity signals. One can either visualize each subimage from the stack by moving along the z axis, or create a projection that sums the matrix along the z axis.

Quantification of anisotropy of microtubule arrays
The space is subdivided into cubes of arbitrary size, typically segmenting the structure into circa 27 cubes.