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.
It is well known that the cortical microtubule (hereafter abbreviated to MT) cytoskeleton in plant cells plays a decisive role in controlled cell expansion and oriented cell division, which together drive the plant morphogenesis [1]. In the absence of MT organizing centers like centrosomes [2], higher plants establish an ordered array of MTs at the cell cortex, the so-called cortical array (CA) [3]. Recent studies have revealed that cell shape may have influence on the orientation of the CA [4]. On the other hand, the orientation of the CA controls cell expansion and cell anisotropy, by guiding the deposition of cellulose synthase complexes along the MTs [5–8]. Through this coupling, the CA in turn can influence the cell shape, essentially setting up a morphogenetic feedback loop. This loop is possibly also amplified by a mechanical feedback mechanism discussed in [9]. Therefore, understanding both anisotropic cell expansion and oriented cell division, requires understanding the formation of the ordered CA from an initially disordered state just after cell division [10]. For non-growing cells, the influence of cell shape on the formation of the CA will be static in nature and thus determined by geometrical features alone. A significant variety in these features is observed between different cell types (see Fig 1). For growing cells, the evolving cell shape may generate a corresponding dynamic influence on the CA formation process. If growth is slow compared to the collective dynamical time scale of the MTs, a quasi-static approach which samples cell shapes at different time points may still be a reasonable approximation. MTs are highly dynamic and filamentous protein polymer aggregates, and form one of the principal components of the plant cytoskeleton [11]. MTs have structurally two distinct ends—a minus-end and a plus-end. The plus-end can dynamically switch from a growing state to a shrinking state or vice-versa. Switching of a MT plus-end from a growing state to a shrinking state is called catastrophe while the reverse switching of a shrinking state to a growing state is called rescue. This phenomenon of reversible switching of MT plus-ends between two states is called dynamic instability. On average, the minus-end of an unstabilized MT continually is in a shrinking state. Thus, the combination of overall growth at the plus-end and shrinkage at the minus-end seemingly moves a MT as a whole. This motion is called treadmilling and has been observed in both in vitro [12–14] and in vivo [15]. In contrast to animal cells, plant cells do not have a well defined MT organizing center. Instead MT activity is dispersed over the whole cell cortex, driven by the localized nucleation of new MTs by γ-tubulin complexes [16]. The cortical MTs are confined to a thin layer of cytoplasm just inside the plasma membrane of the plant cell and are attached to the cell envelope, ensuring that MTs do not translate or rotate as a whole [17,18]. In spite of their fixed attachment to the cell cortex, cortical MTs do show mobility which is due to treadmilling motion [17,19]. Two dimensional attachment to the cell cortex allows MTs to interact with each other via collisions, which occur when the polymerizing plus-end of a growing MT encounters a pre-existing MT. Depending on the value of the collision angle, three different possible events are observed [20]. For shallow angles (≲ 400), a growing MT bends toward the direction of the MT encountered and this kind of adaptive event is called zippering. For steeper angles (≳ 400), the encounter may lead to a so-called induced-catastrophe, where the initially growing MT switches to a shrinking state. Alternatively, the growing MT may slip over the one encountered, leading to a crossover event (see Fig 2). In vivo imaging of cortical MTs has revealed that they nucleate at the cortex, either from isolated nucleation complexes or from pre-existing MTs, and gradually develop from an initially sparse and disorganized state into a final ordered array over a time period of an hour [17,18,21–24] after the previous cytokinesis. MTs have a finite lifetime and ultimately disappear by shrinking to zero length. Plant cells within tissues typically have well defined relatively flat faces, bordered by cell-edges of significantly higher curvature [25]. In root epithelial cells microtubules have been observed to undergo catastrophes when they encounter these cell-edges with a probability that increases with increasing curvature of the cell-edge [26]. We will denote these as edge-catastrophe events. While existing experimental studies inform us about the molecular events and key parameters that are involved in CA formation [17,19,24,27,28], it represents a complex emergent phenomenon that is the macroscopic result of a large number of stochastic microscopic events involving large numbers of MTs. For a full mechanistic understanding of biological processes of this type, mathematical and computational modelling has by now been recognised as indispensable [29]. From the outset, starting with the seminal work of Dixit & Cyr [20], attempts have been made to model CA formation (for reviews on the different approaches involved please consult [30,31]). From this work a consensus hypothesis has emerged to explain the spontaneous order observed in the CA. This mechanism, also dubbed ‘survival of the aligned’ [32], can account for the statistically robust spontaneous alignement of MTs under well-characterized conditions, and appears to be consistent with observations. We therefore adopt this hypothesis as the starting point in the present work. As far as the influence of cell geometry on CA formation is concerned, modelling studies have so far been limited to two dimensional (2D) approaches where the effect of a closed shape in three dimension (3D) was mimicked by boundary conditions [33], or geometrically simple shapes where boundaries are readily described by mathematical expressions, such as cubes [34] and cylinders [35]. Also, some of the simulations reported relied on finite time step brownian dynamics, at most using a Gillespie-type event-sampling algorithm [36] for the spontaneous changes in MT dynamical state. This severely limits the time-efficiency of these simulations, and makes the acquisition of results for wide ranges of parameters with sufficient statistics a very difficult task. Tindemans et al. [35] therefore developed a fully event-driven algorithm (for a general introduction see [37], Chapter 6. 1), in which not only the spontaneous stochastic state changes but also the collisions between MTs are performed by directly sampling from a dynamic list of potential future collision events and their associated waiting times, which is constantly updated. This achieves a speed up of at least three orders of magnitude with respect to standard finite time step algorithms, and eliminates errors in the location of the collision points as well. Here we show how the algorithm of [35] can be implemented efficiently on arbitrary triangulated surfaces, such as those obtained by segmenting 3D reconstructed confocal microscopy images of cells. Strikingly, the additional computational overhead incurred by having to determine the correct passage from triangle to triangle in a manner consistent with the 3D path, is almost completely balanced by the speed up of the collision predicting algorithm due to the strong localization of the search area. In this way the new implementation is almost as fast as the original implementation, and allows the rapid exploration of many different parameter settings and/or geometries each sampled independently many times to obtain statistically reliable averaging over the stochastic ensemble. In addition, the MT dynamical parameters can be chosen independently within each triangular domain and passage probabilities can be associated with the shared edges between the triangular domains. This allows us to easily implement biologically relevant effects such as the observed edge-catastrophes at locations of strong curvature, or potential differences in local MT dynamics due to developmental distinctions between cell-faces, e. g. between faces created by recent cell divisions and faces that have expanded by growth.
In Fig 3, we illustrate the procedure of transforming experimentally obtained confocal image into simulation input. First, using the image processing software morphographX [38], we segment the experimentally obtained confocal images and extract the different cell shapes. For segmentation, confocal image stacks (tif) were Gaussian blurred using sigma value 0. 3 μm, subsequently we applied the ITK watershed auto-seeding with level threshold value in the range 500–1500. Then, we approximated the segmented cells by creating triangulated surface meshes using marching cubes 3D with cube size in the range 0. 5–1. 0 μm. Using the surface mesh processing software meshLab [39], we tag the different cell-edges and cell-faces, which are represented by appropriately coloured associated triangles. Finally, all this information is stored in a single mesh file, where each triangle has a cell-edge tag with the values of the associated cell-edge angles and cell-face tags (see S1 File, Sec. SI. 1). Confinement of the dynamics of cortical MTs to the surface (cortex) of plant cell effectively reduces the system to two spatial dimensions. In our modelling framework, we exploit this advantage by treating the cell shape as a 2D surface, embedded in the ambient 3D space. Because we are working with a triangulated approximation of surfaces, we developed an algorithm to establish a connectivity graph of the triangles [40] in order to appropriately propagate MT dynamics between adjacent triangles. Using a combination of rotations, implemented by quaternions [41], and z-axis translation, we first transform the 3D-triangulated surface from the 3D x − y − z space to the 2D x − y plane (see Fig 4). Each individual triangle has a unique quaternion operator which is fully determined by its normal vector in 3D, which operates on all its vertices simultaneously. Different triangles have different normals, and hence end up in different, typically disjoint, 2D positions after rotation, irrespective of whether they had shared edges or not, but this position is irrelevant to the dynamics within each triangle. We simulate MT nucleation, growth, shrinkage and collision events within each individual triangle (see Fig 5). In order to nucleate new MTs on the surface, we first randomly select a triangle from the entire set and generate uniformly distributed points within the selected triangle [42]. The initial MT growth direction is made isotropic by choosing a random angle of nucleation in the range [0,2π]. The persistence length of a MT is much longer than its average length, typically of the order of a few μm’s, [43,44] and therefore much larger than the typical dimensions of our triangles (well below 1μm), allowing us to model the individual MT within a triangle as an elongating straight line segment. The location of the collision point between two MTs is then readily determined by calculating the intersection point of their respective trajectories [45,46]. To deal with the various MT events, we implement an event-driven simulation technique [35]. The events are separated into two categories: stochastic events and deterministic events. The stochastic events associated with MTs are independent and result from Poissonian random processes specified by the nucleation rate rn, the spontaneous-catastrophe rate rc and the rescue rate rr. The deterministic events associated with the MTs include collisions, disappearance, intersection with triangle-edges, and simulation control events (i. e extracting simulation output at fixed time intervals). In between events, the plus-ends of MTs either grow or shrink with velocities v+ and v− respectively, while the minus-end retracts with the treadmilling speed vtm, so that all length changes between events are readily computed. We use triangle-edge to triangle-edge links (obtained from the connectivity graph) to propagate MT dynamics from one triangle to its relevant neighbour. For example, in Fig 6 ΔABC is connected with ΔADC through the shared edge AC, which are linked. The proper transition of a MT from ΔADC to ΔABC through the edge AC on the 3D surface requires, in the 2D plane, Using an affine transformation [47] on the coordinates of the growing MT plus-end tip at point P′ (x′, y′) of edge AC (ΔADC), we translate (see Fig 6 left panel) the tip to the edge AC (ΔABC). Next, we rotate (see Fig 6 right panel) the MT trajectory with rotation angle, θ R = arccos (A C → (Δ A B C). A C → (Δ A D C) | A C → | 2) (1) This rotation in the 2D plane implements the assumption that a MT grows along the straightest possible path on the 3D surface, also known as a geodesic. Although the mechanical energy of bending could conceivably also play a role in determining the path of the MT, potentially causing a more complex coupling to the geometry, we have, in the absence of sufficient experimental data discriminating these effects, chosen to opt for this arguably simplest rule. During transition of a MT from one triangle to another triangle, we implement the probability of local edge-catastrophe based on the local value of triangle-edge angle between the adjacent triangles (see S1 File, Sec. SI. 2). We also use a cell-face tag, which is assigned to the triangles during the input surface mesh creation (see S1 File, Sec. SI. 3), to implement local differences in stability of MTs e. g. by changing the local spontaneous-catastrophe rate. More generically, our simulation framework allows arbitrary domain specific parametrisation, with the smallest domain being a single triangle. This flexibility allows full freedom to incorporate experimental observations on local MT dynamics or test additional relevant hypotheses. Finally, whenever required we apply the inverse quaternion transformation on the set of 2D triangles and the MT segments present on them, to reconstruct the 3D surface with the MT configurations realised. Through this transformation, segments of a MT which appear to be disconnected in the set of 2D triangles, get (re) connected on the 3D surface. As illustrated in Fig 7, the segment M 1 M 2 ¯ is attached to ΔABC and the segment M 2 M 3 ¯ is attached to ΔADC in the 2D plane. After the inverse quaternion mapping, these two segments get connected in the 3D surface and represent a single MT with end points (M1, M3). To measure the degree of order and the orientation of the MT array, we define two order parameters. The first parameter is a scalar (Q (2) ), which measures the average degree of ordering of the MTs, and the second a vector Ω ^, which indicates the global orientation of the array. Both these parameters are derived from a single tensorial order parameter Q (2) (see S1 File, Sec. SI. 4), Q (2) being the absolute value of the smallest eigenvalue of this tensor and Ω ^ the corresponding normalized eigenvector, i. e. | Ω ^ | = 1. For completely random orientation of MTs Q (2) ≈ 0 and for well organised MTs that form an array, Q (2) ⪅ 1. The direction of Ω ^ is perpendicular to the average local orientation of the MTs (see Fig 8), i. e. perpendicular to the plane onto which the total projection of individual MTs is maximal on average. In Fig 9, we present an overview of the entire simulation approach. In order to obtain a suitable set of MT dynamical parameters for our simulations, we make use of the control parameter G, introduced and used in [48,49]. This parameter encapsulates the effect of all six single MT dynamics parameters (speeds v+, v− and vtm, dynamical instability rates rc, rr and nucleation rate rn) into a single magnitude that controls the frequency of MT interactions which are required to ensure a spontaneously ordered state. Its explicit form identifies it as the ratio between two length scales G = − (2 (v + − v t m) (v − + v t m) r n (v + + v −) ) 1 3 (r r v − + v t m − r c v + − v t m) − 1 = − l 0 l a v g. (2) Here l0 is the MT-MT interaction length, and lavg is the average length of MTs in the absence of any interaction effects. The statement that G is the single control parameter was only strictly derived for spatially homogeneous and open, boundary-free domains. To allow for possible effects of a closed cell geometry, we chose to independently vary the two factors involved: l0 by varying the MT nucleation rate rn and lavg through varying the spontaneous-catastrophe rate rc. Moreover, and unless explicitly mentioned, we implement homeostatic control of the MT growth speed v+ by implementing a finite tubulin pool (see S1 File, Sec. SI. 6 for details), as this speeds up the relaxation towards a steady state. The remaining simulation parameters are taken from [35]. We tested the parameter choice on a spherical surface with radius r ≈ 6 μm, inspired by the typical dimensions of a Arabidopsis thaliana embryonic cell. This is at the smaller end of the spectrum of plant cell sizes and therefore the most stringent test on finite size effects. We found that we could obtain robustly ordered arrays with Q (2) ≳ 0. 70 at a state point with l0 = 2. 05 μm and lavg = 417. 5 μm, corresponding to G ≈ −0. 005 (see Fig 10). The fact that G < 0 ensures bounded growth regime, where the length of MTs is always finite. All parameters used are summarised in S1 File, Sec. SI. 5 (Table 1). It is important to note that lavg is the theoretically expected value of average MT length assuming absence of any interaction effects. However, in simulation MTs will interact with each other leading to changes in effective values of MT dynamics parameters, i. e. modification in the values of over all catastrophe rates due to effect of induced-catastrophe. These dynamic changes in MT dynamics parameter will also modify the average MT length to a value l ¯ different from lavg. In the coming sections, average MT length will correspond to the simulated value of average MT length, i. e. l ¯. For further clarification on this point, variation in the value of l ¯ with respect to lavg is described in S1 File, Sec. SI. 6. Finally, we did not take any cell cycle related variations in MT dynamics into account, as this work focusses on understanding the relative contributions of the basic MT-MT interaction rule, shape anisotropy, edge-catastrophe and differential MT stabilization at developmentally different cell-faces on CA formation in interphase. We note, however, that such time-dependent effects could readily be implemented within our framework. The simulation code is written in C++, and compiled using GNU GCC 4. 8. 4 compiler. Development and incidental runs were done on a standard PC equipped with Linux. Production runs were performed on the Dutch National Grid infrastructure comprising 10000 distributed compute nodes. Individual runs of 10 hours of biological time took between 20 seconds to 120 seconds wall clock time, depending on the parameter settings.
We performed a set of simulations validating our computational framework. We first tested to what extent the triangulation of otherwise smooth surfaces influences the final results. To do so we considered the case of approximating a perfect sphere, varying both the number of triangles employed, as well as the triangulation method. The key desiderata in this case are (i) that in steady state, due to the spherical symmetry, there should be no bias in the orientation Ω ^ of the ordered array, and (ii) that the value of the scalar order parameter Q (2) does not depend on the nature of the triangulation. Next, in order to test the system as a whole against a known, and geometrically non-trivial case, we reproduced previously reported results on a cubic surface [50]. We are now in a position to use our framework to study the role of cell shape on MT array formation. We first do this in a geometrically simple setting, which allows us to systematically disentangle the role the various shape-related factors involved. Then we turn real cell shapes, provided by the leaf pavement cells of Arabidopsis thaliana, Hedera helix and Nicotiana benthamiana.
Arabidopsis thaliana and Nicotiana benthamiana plants were grown as described in [76]. Leaves from these plant species were harvested and incubated in Renaissance staining as described in [77], samples were then mounted in water to visualized their cell shape. 35S: : MAP65-GFP expression in Arabidopsis thaliana leaves was analysed using 3 days old seedlings and 4 weeks old plants respectively. Expression assays in Nicotiana benthamiana leaves were performed using 4 weeks old plants. Agrobacterium containing 35S: : TUB-mCHERRY were infiltrated as described in [78]. Arabidopsis thaliana leaves were stained with Propidium idodide (10μg/ml) and imaged using a C- Apochromat 40 X/1. 20 W Korr water immersion objective of a LSM 880 airyscan Zeiss confocal microscopes. GFP was excited using an argon laser 488nm, and fluorescence emission was detected from 500 to 540 nm. Renaissance staining for Nicotiana benthamiana leaves was detected at 415–440 nm excitation and a 405 nm beam splitter. For 35S: : TUB-mCHERRY localization in Nicotiana benthamiana, infiltrated leaves were mounted in water and mCHERRY was detected with 543 nm excitation and 488/543/633 beam splitter.
In contrast to animal cells, plant cells are encased in a rigid cell wall providing among others the necessary rigidity to contain the turgor pressure that allows the plant as a whole to raise itself against gravity. In order to maintain this mechanical integrity while growing requires plant cells to rigorously control their expansion. In this they are aided by a plant-unique intracellular structure called the cortical array, an organised set of long filamentous protein polymers called microtubules attached to the inside of the cell membrane. Experiments and modelling studies have revealed how collisions between these highly dynamic polymers promote the spontaneous alignment of microtubules observed in cells. However, it was hitherto unknown how the shape of cells influences this self-organisation process and helps in determining the overall orientation of the cortical array. The latter is key in deciding the axis of cell expansion and the orientation of cell division planes. Here we develop, validate and apply a simulation framework for studying the dynamics of interacting microtubules on arbitrary triangulated surfaces that can be extracted from three-dimensional reconstruction of real cells from confocal microscopy images. This allows us to start unravelling the delicate interplay between cell shape, effects of sharp cell-edges and developmental cues on the organisation of the cortical array.
plant anatomy
microtubules
pavement cells
cell cycle and cell division
plant cell biology
cell processes
brassica
condensed matter physics
plant science
model organisms
experimental organism systems
plants
flowering plants
cellular structures and organelles
cytoskeleton
research and analysis methods
arabidopsis thaliana
nicotiana
leaves
physics
eukaryota
plant cells
plant and algal models
cell biology
biology and life sciences
cellular types
physical sciences
nucleation
organisms
A computational framework for cortical microtubule dynamics in realistically shaped plant cells