A computational framework for cortical microtubule dynamics in realistically shaped plant cells

Plant morphogenesis is strongly dependent on the directional growth and the subsequent oriented division of individual cells. It has been shown that the plant cortical microtubule array plays a key role in controlling both these processes. This ordered structure emerges as the collective result of stochastic interactions between large numbers of dynamic microtubules. To elucidate this complex self-organization process a number of analytical and computational approaches to study the dynamics of cortical microtubules have been proposed. To date, however, these models have been restricted to two dimensional planes or geometrically simple surfaces in three dimensions, which strongly limits their applicability as plant cells display a wide variety of shapes. This limitation is even more acute, as both local as well as global geometrical features of cells are expected to influence the overall organization of the array. Here we describe a framework for efficiently simulating microtubule dynamics on triangulated approximations of arbitrary three dimensional surfaces. This allows the study of microtubule array organization on realistic cell surfaces obtained by segmentation of microscopic images. We validate the framework against expected or known results for the spherical and cubical geometry. We then use it to systematically study the individual contributions of global geometry, cell-edge induced catastrophes and cell-face induced stability to array organization in a cuboidal geometry. Finally, we apply our framework to analyze the highly non-trivial geometry of leaf pavement cells of Arabidopsis thaliana, Nicotiana benthamiana and Hedera helix. We show that our simulations can predict multiple features of the microtubule array structure in these cells, revealing, among others, strong constraints on the orientation of division planes.

In a triangulated surface, such faces (i.e. F 1 , F 2 ) as well as their adjacent edge (i.e. E), will be approximated by a large set of triangles. Therefore, E will be region composed of many triangles. As illustrated in Fig SI.1, we give a unique color code f 1 (blue) to Faces F 1 (blue colour)and F 2 (green colour), meet along the dotted line, (B) Using this dotted line as reference, an edge E (red colour) is detected, which is composed of multiple triangles belonging to either F 1 or F 2 . Right panel: Triangles with edge color (e c ; red) and one of the face colour (f 1 ;blue or f 2 ;green), are identified via the respective face colour only. the triangles that belong to F 1 . Similarly, we give a unique color code f 2 (green) to the triangles that belong F 2 . The triangles belonging to the edge (E) are assigned an additional color code e c (red). Now, by our definition all the triangles with color code e c will also have either a color code f 1 or f 2 . Based on these triangles with dual color code, we calculate the average (approximate) value of the edge angle as where,n (f1,ec) i is the normal to the i th triangle with dual color code (f 1 , e c ) andn (f2,ec) j is the normal to the j th triangle with dual color code (f 2 , e c ). N 1 and N 2 are total number of triangles with dual color code (f 1 , e c ) and (f 2 , e c ), respectively.
We tested the accuracy of this approximation in calculating the edge angles for regular shapes such as cylinder, cube or isosceles pyramid etc., where the edge angles are known in advance. The calculated value of the average edge angles, θ avg E deviates only up to ≈ 5% from the expected values for these shapes. Therefore, this approximation approach can be used as working formula for edge angle calculations in triangulated meshes. As the calculation of θ avg E is a a fairly laborious process, we aim to fully automate it in the future.

SI.2 Implementation of edge-catastrophe
Consider a MT entering an edge region E with average edge angle θ avg E and attempting to pass from face F 1 to face F 2 . We can virtually construct its potential future trajectory M, using the triangle-to-triangle propagation rules. To mark the actual transition point from the one face F 1 to the other face F 2 , we consider a curve E which connects the location of maximal curvature within the edge region between F 1 and F 2 , which are provided the meshing software [39]. At the intersection point P between M and E we construct the unit tangent vectorsm P andê P to these curves (see Fig SI.2). This allows us to define the incidence angle The incidence angle determines whether the degree of curvature the MT will experience passing from the one face to the other; maximal for θ i ≈ π 2 or minimal θ i ≈ π 2 . We then define the effective bending angle (θ b ) the MT experiences in passing from F 1 to face F 2 as i.e. the geometrical expression one obtains if the the smooth parts of faces F 1 and F 2 were joined directly at a sharp, cusp-like edge. This form guarantees that the maximum achievable value of the bending angle is θ max b = θ avg E , which occurs at θ i = π 2 , and the minimum value is θ min b = 0, which occurs at θ i = 0. This bending angle θ b is then used as input for the actual edge-catastrophe calculation. Total edge-catastrophe Geometrically speaking, for a completely flat edge θ b = 0, here MT will make a transition through the edge without any edge-catastrophe, i.e. probability of edge-catastrophe P E = 0. An edge with maximum bending (θ b = π) will completely prohibit passage of a MT (P E = 1). Based on these two extreme case scenarios, we formulate the following relation of geometry based probability of edge-catastrophe However, experimental observation [26] suggests that at an edge of steep curvature with a value of θ b ≈ π 2 , MTs are already prohibited from crossing that edge, i.e receive full PLOS 2/10 edge-catastrophe, i.e we are allowed to set P E = 1 at θ b = π 2 . Accordingly we define a working form of the probability of edge-catastrophe (P W ork where 0 E cat 2, is an edge-catastrophe multiplier, specific to a certain edge and its value depends on the biomechanical properties of that edge. The experimental observation that drives the foundation of Eq SI.6 is well fitted by E cat = 2.

Local distribution of edge-catastrophe
For a pair of triangles (T 2 ) at their shared edge (k), we calculate the local edge angle θ are respective normals to T  through their shared edge k andt (k) is a unit vector along this edge.
Considering the local incidence angle θ and T (k) 2 , we define the effective local bending angle where,m (k) is a unit vector along the growth direction of a MT passing from T respectively, and both pointing either inwards or outwards to the triangulated surface.
Using Eq SI.2 and Eq SI.8, we distribute the purely geometry based total edge-catastrophe P Geom E (see Eq SI.5) among each pair of triangles (T 2 ) by a weight factor at their Taking into account of the edge-catastrophe multiplier E cat (see Eq SI.6), we define the effective local probability of edge-catastrophe (SI.10)

SI.3 Implementation of MT stabilization
We illustrate the implementation technique of MT stabilization in Fig SI.4. All the triangles belonging to a face F 1 are assigned with a rate of spontaneous catastrophe r F1 c and the triangles belonging to an adjacent face F 2 are assigned with a rate of spontaneous Fig SI.4. Implementing differential MT stabilization. When a MT passes from a triangle T i of a face F 1 to a triangle T j of another face F 2 , we update its spontaneous In a triangle, a lower value of spontaneous catastrophe (r c ) results in a higher average length of the associated MT segments, leading to a longer lifetime of the segments, in comparison to another triangle with higher value of r c . This enhancement in lifetime can be interpreted as endowing enhanced stability of the MT segments by the corresponding triangle. This way, we implement the idea of difference in stability of MTs at different faces of the triangulated surface, via face specific different rate of spontaneous catastrophes.

SI.4 Order parameter tensor
Consider a (piecewise) smooth oriented 2-manifold M embedded in R 3 , which we call the surface. Except for possibly a set of areal measure zero, we can endow a tangent space at each of its points σ with an orthonormal reference frame defined by two unit vectorsê 1 (σ) andê 2 (σ), using the prescription that the outward normal to the surface is given byn (σ) =ê 1 (σ) ×ê 2 (σ).
We follow the convention that lower in italic indices take on the values {1, 2} and we adopt Einstein summation convention. Components of the tensor q (m) (σ) upto order two are explicitly given by where . . . σ denotes orientational average, i.e. if distribution of possible orientations at σ is coded into a normalized orientational distribution ψ σ (α), then the orientational average is On the local tangent space, the tensor components are symmetric and traceless for m ≥ 2.
We can average the local order parameter tensors over the full surface to extract the global order parameter tensor where ρ(σ) is the local areal density of the orientable constituents present on the elementary area dA(σ). For m ≥ 2, symmetric and traceless properties of q (m) remain conserved by the surface averaging process, i.e Q (m) is also symmetric and traceless.

Extension to triangulated surface
Consider a triangulated three dimensional surface, where a bunch of rod like particles (line segments) are distributed over each triangle of the surface. For a triangle E, at its PLOS 5/10 center we assign a tangent plane with an orthonormal reference frame defined by the unit vectorsê 1 ,ê 2 andn =ê 1 ×ê 2 . Lets assume on the plane of this triangle, a particle of length l i (E) is oriented with angle θ i , with respect to the direction ofê 1 . In the local tangent frame, corresponding orientation vector (Eq SI.11) of this particle iŝ (SI.14) In matrix form If q (2) (E) is the local order parameter tensor associated with the triangle E, then the associated tensor components (Eq SI.12) are where l max is the total number of particles present in the triangle E. We avoid explicit influence of geometry on MT order and array orientation, and also keep the resulting tensor trace less by setting q This local order parameter tensor q (2) (E) is defined at the local tangent frame attached to the triangle E. To express this tensor with respect to the global coordinate frame, first we express the basis vectors (ê 1 ,ê 2 ,n =ê 1 ×ê 2 ) of the local tangent frame with respect to the global frame (x,ŷ,ẑ) aŝ 1ẑ , In matrix form

PLOS
6/10 Using the rotational transformation, we express the local order parameter tensor q (2) (E) with respect to the global coordinate frame as For all the triangles of the three dimensional polyhedral surface, averaging over the local order parameter tensors q r (2) , we will get the global order parameter tensor (Eq SI.13) where L j is length of all the particles present on the triangle E j and N is the total number of triangle used to approximate the polyhedral surface. In simulation, we will refer the global order parameter tensor Q (2) as the order parameter tensor.

SI.5 Simulation parameters
treadmilling speed 0.01µm sec −1 r r rescue rate 0.007 sec −1 r c spontaneous-catastrophe rate variable: 0.003-0.02 sec −1 r n nucleation rate variable: 0.001-0.01 µm −2 sec −1 ρ tub finite tubulin pool density (length equiv.) 10 µm −1 θ z angle of zippering 40 0 p ind−cat probability of induced-catastrophe 0.5 p cross probability of crossover 0.5 Table 1. Overview of the MT dynamics parameters and variables with their default values (if applicable). The dynamic instability parameters were used from [19] and the value of v tm has been approximated from [17]. The nucleation rate (r n ) has been chosen primarily from [33], however also varied to maintain sufficient number of MTs, such that the system always reach to an ordered state. The angle of zippering (θ z ) and probabilities for zippering (p z ), crossovers (p cross ) and induced-catastrophes (p ind−cat ) were deduced from [20] , combining data from MBD-DsRed and YFP-TUA6 labeling. In the biological measurements, no distinction was made between spontaneous-catastrophe (r c ) and induced-catastrophe, resulting in a likely high value by an unknown factor. Based on this observed uncertainty and the linear dependency of G on r c , we used r c as the primary mean to vary G.

SI.6 Finite tubulin pool effect
In most of the simulations presented in this article, we used a finite tubulin pool, which means the speed of a growing MT plus-end at any time t, is dependent on the collective total length L(t) of all MTs in the system where L max = ρ tub A is the length equivalent of the total amount of tubulin on the surface of area A and finite tubulin density ρ tub , and v + is the simulation input MT plus-end growth speed under the assumption of an infinite tubulin pool (see Table 1).

7/10
Generically, following Eq SI.17, the speed of the growing plus-end v + tub (t) will be less than v + . Such a decrease in growth speed of MT plus-end will lead to a corresponding decrease in simulated average MT length (l). As a consequence, the actual value of the control parameter (G) will tend towards a more negative value than the initially set value. First we considered an infinite tubulin pool (ρ tub = ∞ µm −1 ) and systematically varied the value of G to achieve a steady state value ofl with sufficient degree of MT order Q (2) . As illustrated in Fig SI.5 (A), we identified a regime, −0.06 G ≤ −0.12, which meet both the requirements. Here, G was calculated by using v + , which by default assumes infinite tubulin pool. Next, we considered a system with finite tubulin pool (ρ tub = 10 µm −1 ) and initiated simulation of MT dynamics with G ≈ −0.005, which was again calculated by using v + which assumes infinite tubulin pool. Over time the system 'eats up' the available tubulin pool, resulting in a modified value of MT plus-end growth speed from v + to v + tub (t) (see Eq SI.17). Calculation of G by using steady state value of v + tub (t = T max ) resulted in a modified value from G ≈ −0.005 to G ef f ≈ −0.05, where both corresponds to the set of same input MT parameter values and is the working domain in Fig 10. T max sets the upper limit of simulation time which assured steady state MT order and array formation. The interesting point to note here is the value of control parameter G ef f ≈ −0.05 with finite tubulin pool to be equivalent to that with infinite tubulin pool, i.e. −0.06 G ≤ −0.12, suggesting that effective value of working domain remains invariant under the assumption of finite tubulin pool. In Fig 10, the values of l avg which are in order of 10 2 µm, may at first sight seem unrealistically large. They are, however, the bare (theoretical) values in the absence of both pool limitations effects, as well as additional induced-catastrophes by the MT interactions. Here, they are used solely for the purposes of establishing a proper parametrization of the simulations in which all the relevant effects are included. The actual simulated average MT length (l) are in the order of a few µm's, as illustrated in Fig SI.5 (B).

SI.7 Triangulation effect on MT-array orientation distribution
In terms of the polar angles (θ, φ) the infinitesimal element of surface on the sphere is given by dS(θ, φ) = sin θdθdφ. Introducing the cartesian coordinates we see that the infinitesimal element of surface area in terms of these coordinates becomes the Euclidean area measure dS(x, y) = dxdy. We can therefore homogeneously map the surface of the sphere to the rectangular planar domain [−1, 1] × [0, 2π] ∈ R 2 . Using