Microtubule simulations in plant biology: A field coming to maturity

The plant cortical microtubule array is an important determinant of cell wall structure and, therefore, plant morphology and physiology. The array consists of dynamic microtubules inter-acting through frequent collisions. Since the discovery by Dixit and Cyr (2004) that the outcome of such collisions depends on the collision angle, computer simulations have been indispensable in studying array behaviour. Over the last decade, the available simulation tools have drastically improved: multiple high-quality simulation platforms exist with specific strengths and applications. Here, we review how these platforms differ on the critical aspects of microtubule nucleation, flexibility, and local orienting cues; and how such differences affect array behaviour. Building upon concepts and control parameters from theoretical models of collective microtubule behaviour, we conclude that all these factors matter in the debate about what is most important for orienting the array: local cues like mechanical stresses or global cues deriving from the cell geometry.


Introduction
Many mechanical and other material properties of plant parts and plant-based products derive from the structure of the cell walls that surround all plant cells.Deposition of the main load-bearing component, cellulose microfibrils, is primarily guided by cortical microtubules [11,30,56].Plant cortical microtubules were discovered in 1963 [23,44], about twenty years before the first reports of cortical microtubules in specific cell types of animals (starfish oocytes in 1984 [54] and, less convincingly, Tubifex eggs in 1981 [65]).The cortical microtubules themselves are thin, stiff and dynamic protein structures that grow and shrink by polymerization and depolymerization from their ends, respectively, with frequent switches between the growing and shrinking state called catastrophe and rescue.They are attached to the inside of the cell membrane and, therefore, interact through frequent collisions [16,22,75].Ever since the discovery by Dixit and Cyr in 2004 [20] that the outcome of such collisions depends on the relative collision angle e often simplified as bundling/zippering/entrainment for angles < 40 and crossover or induced catastrophe for larger angles [3,15,16,25] e theoretical and simulation studies have been indispensable in understanding several key processes from generic array alignment [3,15,18,25,33,70,71], to (re)orientation, e.g., of dark grown hypocotyls in response to blue light [47,52,59,60], and formation of additional structure like the banded patterns in protoxylem [38,39,63].The cortical localization of these microtubules has been a great advantage for modelling approaches, as the system can be described very realistically as a strictly 2D surface 1  [3,10,15,25,45,60,63,70,71], as well as the experimental quantification of microtubule dynamics [22,23,75].
The first debate, whether bundling or induced catastrophes are the most important for spontaneous alignment of the array [3,25,70], was initially settled in favour of induced catastrophes with help of a strong theoretical analysis [33,70] popularized as "survival of the aligned": spontaneous alignment occurs with sufficient interactions per microtubule lifetime, eor, alternatively worded, when the ratio G ¼ Àl 0 =l between interaction length scale (l 0 ) and intrinsic (average) microtubule length ðl Þ is sufficiently large2 e where the threshold for sufficient (G * ) depends on Fourier coefficients of the "interaction function" ehow collision outcomes depend on relative collision anglese specifically, on both the overall probability of induced catastrophes and how strongly these are focused on large angle collisions.Colloquially, the transition point depends on how strongly and specifically microtubules deviating from a majority direction are penalized (for more explanation on G and G * , see [17,70]).
The theory was later extended to explain the alignment promoting effect of katanin severing at crossovers [46,77] by understanding severing events as delayed and, therefore, weaker-acting induced catastrophes [18].Analytical theories are useful tools for developing a mechanistic understanding of the system, but are restricted to the analytically tractable.The "survival of the aligned"-theory has been derived either with retracting minus ends [15,71] (hybrid treadmilling [64]) or with bundling [33,70], but not (yet) with both at the same time.With both processes combined in simulations, bundling suddenly increases the regime of alignment beyond what is achieved by induced catastrophes alone [17,71].The control parameters (like G) that are derived with the analytical models, however, remain powerful tools in designing insightful simulation strategies and interpreting their results [16,17,63].Interestingly, the observation of [33,70] that bundling did not affect the alignment transition was recently also reported for "cytosolic" interacting microtubules in full 3D [29], notably, again with static minus ends [33,70].This historical consideration clearly shows that model choices, which often consider analytical and/or computational tractability, can have a major impact on conclusions.
There often is a trade-off between biological realism and computational cost.By educated guess, current simulation platforms differ at least hundred-fold in computational cost, meaning computational cost itself impacts which questions can be addressed [17].
We compare how recent key simulation models for plant microtubules differ on multiple relevant aspects (Figure 1).We have recently found that realistic nucleation and microtubule flexibility matter for array behaviour [38,39,58].The mode of membrane anchoring affects the behaviour of individual microtubules including their apparent flexibility as well as the strength of their interaction [6, 48,69].Finally, modelling studies differ greatly in the mechanisms they use to adjust microtubule orientation towards mechanical stresses in the cell wall [31], the actual mechanism of which remains a wideopen question.
CorticalSim works with an efficient event-driven algorithm, focusing on cortical microtubules on 2D surfaces and modelling them as sets of connected line segments.Its computational efficiency makes it ideal for running large numbers of simulations, providing more statistically robust conclusions and bringing new questions within reach [17].Cytosim operates on full 3D simulation domains, modelling microtubules as semiflexible objects using discrete time steps.This model is highly realistic and computationally demanding, originally developed for studying relatively small numbers of microtubules in animal cells.Finally, Tubulaton is a 3Dbased model with discrete time steps, treating microtubules as composed of fixed-length subunits.Its computational cost is in between CorticalSim and Cytosim.In terms of realism and biological detail, CorticalSim and Tubulaton are comparable overall, with different platforms being more realistic in different aspects.
All simulation models differ from the analytical theory in that they are based on discrete particles (microtubules) with explicitly stochastic behaviour, whereas the analytical models [29,33,70] describe continuous distributions of microtubule densities (as a function of length, orientation) using differential equations, without explicit stochasticity.Consequently, the analytical theory does not cover the spatial fluctuations in, e.g., microtubule density that arise from their interactions.Simulations confirm, however, that an exact prediction for G * can be obtained in the limit of weakly interacting microtubules [71].
The three models we review here are on the scale of many interacting microtubules within a single cell.Whereas detailed simulations of individual microtubules e.g., Refs.[1,6,41,69] study the mechanisms behind microtubule dynamics, whole-cell models explore their consequences, and are in turn too detailed for whole tissue processes.Nevertheless, they may be used to derive coarse-grained descriptions of array behaviour to improve tissue-level models like Refs.[31,35].

Nucleation
Microtubules in the interphase cortical array are predominantly nucleated by g-tubulin ring complexes (g-TuRC) [50], with the exception of nearly empty arrays [45].The g-TuRC independent fraction decreases rapidly with increasing density, because the critical tubulin concentration for nucleation is much higher than for elongation [76].Combined with a finite pool of available tubulin, biases in, e.g., absolute nucleation orientation can have a big impact on array orientation [60], and an overall change in the net nucleation rate can change the length distribution of microtubules, e.g., with a direct effect on spindle length [49].On microtubules, g-TuRC nucleates more efficiently [51], with a distinct angular distribution relative to the parent [12,50]   Overview of major microtubule simulation platforms used for plant cells in recent years.References are only given if later publications differ in the relevant aspect from the key reference given at the top.+: implemented; -: not implemented.Grey shade in the left column means the aspect is comparable with the analytical model of Hawkins et al.Together, these findings demonstrate that microtubulebased nucleation has substantial effects on array behaviour.However, microtubule-based nucleation provides more realistic results than isotropic nucleation only when the inhomogeneity problem is avoided.

Microtubule membrane anchoring and array organization
Microtubules are rather stiff polymers, with in vitro persistence lengths (l p ) up to millimetres [34,62], which may contribute to their cortical localization [48].Freely rotating microtubules confined in small containers tend to minimize their bending energy via orientations that minimize their total curvature [14].Cortical microtubules, however, are anchored to the inside of the cell membrane, preventing translational movement.In the mathematical limit of immediate and continuous anchoring, microtubules would only minimize their curvature locally and, consequently, follow geodesic paths.E.g., the locally straight paths on a cylinder: perfectly transverse, longitudinal or a constant-pitched helix [6,71].Reality is intermediate, with anchoring by discrete proteins, giving growing microtubules a short "free" plus-end that could locally sense and minimize its curvature before anchoring.Depending on anchor spacing, this may lead to a substantial force away from geodesics, e.g., implying gradual reorientation from transverse to longitudinal on an infinitely long cylinder [6,69].Interestingly, stochastic attachment of the anchors can result in both smaller and larger deflection than deterministic anchoring at regular intervals with the same anchor density, for large and small average anchoring distance relative to curvature, respectively [69].Indeed, when MT anchoring is reduced by chemical or mutational targeting of Phospholipase Dd (PLDd), the frequency of longitudinal arrays increases the expense of transverse arrays [5,19,28].Note that these papers explain the effect via assumed increased exposure to cytoplasmic streaming rather than as an effect of bending energy minimization over a larger free length.An additional effect of PLDd disruption is an increase of disordered arrays [5,19,28] in line with the effect of reduced interaction intensity predicted by the "survival of the aligned"-theory [71], and confirmed by 3D simulations [48] comparing "weak" and "strong" anchoring to the membrane.The observation of hyperaligned arrays in the clasp-1 mutant in combination with reduced membrane anchoring [4] seems, at first glance, at odds with the theory.It may be explained, however, by the wider range of collision angles resulting in zippering [2,4], which seem to specifically reduce induced catastrophes for intermediate collision angles, leading to a larger aligned regime (more negative G * ) [17,18].
Most simulation models of interacting microtubules implicitly assume infinitely fast and dense anchoring.We expect that the effect of realistic anchoring in these models will be minor in most cases, but may be important in understanding the phenotypes of a few specific mutants, or in parameter regimes where the array is extremely sensitive to any kind of orienting bias.

Microtubule flexibility
From experimental images, cortical microtubules appear less straight than expected with millimetre l p .Correspondingly, recent simulation studies of interacting cortical microtubules differ in whether microtubules are modelled as straight line segments in between interactions [9,10,58,60,63], or as (simplified) semiflexible objects (finite l p ) [21, 29, 38, 48, 49, see also Table 1].The immediate effect of introducing this flexibility is that the resulting arrays look more "natural".How else does flexibility affect array behaviour?
Obviously, the impact will depend on the value of l p , but few relevant measurements exist.Mirabet et al. [48] and Durand-Smet et al. , because the minimal number of interactions for alignment not only has to occur within the intrinsic microtubule length, but also within the (in this case) much shorter l p .As the same G value can be obtained for both relatively long, weakly interacting and short, strongly interacting microtubules [16,17], we expect the largest impact of flexibility for weakly interacting microtubules, i.e., when crossover probabilities exceed those of bundling and induced catastrophes.Actual l p in plant cells is likely at least 60e100 mm, but hard to quantify, as in vitro gliding assays at high and low density show that microtubule interactions severely reduce the apparent l p [62].
In protoxylem simulations, flexibility had a mixed impact on patterning of the array: on the one hand, the array could more easily adopt a pattern that suboptimally "fits" the microtubule array.On the other hand, microtubule density may be lost by microtubules "wandering" into the gaps, with band density possibly dropping below the threshold for alignment with isotropic nucleation [38].In metaxylem, local tuning of l p via MAP70-5, which reduced l p 2e3 fold in vitro [62], allows for more circular gaps compared to map70-5, map70-1 knockouts [57,62].Functionally, the l preducing effect of MAP70-5 allows for a sharp delineation of the gaps by ring-shaped bundles and reduces the range over which the gap impacts the overall array orientation, particularly if MAP70-5 is enriched at the gap boundaries.From a physical perspective, MAP70-5 reduces the energy penalty for inserting a circular gap in an aligned array.
Flexibility could also affect the overall orientation of the array in certain parameter regimes.The shape of the cell tends to favour a limited set of array orientations, of which one could be selected via penalties for crossing edges (e.g., "a clasp mutant") or cell faces with different microtubule dynamics (e.g., hormone-regulated), etc [10,17,58].Microtubule flexibility reduces how far information about cell geometry can propagate through the array, restricting the maximum cell size for such mechanisms.We expect flexibility to enhance the general trend that "local" cues become relatively more important with weaker "global" cues, like when increasing cell size [58].An interesting open question is whether semiflexible microtubules would affect the number of favourable orientations on realistic cell shapes, which tend to have fewer geodesic paths than the idealized cell shapes used in many theoretical studies.
Overall, microtubule flexibility could affect cortical array behaviour in many ways, but whether potential effects are relevant strongly depends on the ill measurable value of l p .

Stress-sensing and orientation-dependent microtubule dynamics
Ever since the first discovery that cortical microtubules can align with local stress patterns in the cell wall [13,31,61], the question has been how this is possible.
Multiple biochemical pathways have been suggested for sensing mechanical stress, but as most proteins involved also have a major impact on microtubule dynamics, mutant phenotypes could also be explained as response problems [32,73].It has also been suggested that the microtubules themselves are sensitive to tensile stress [32], based on in vitro observations that microtubules under tensile stress polymerize faster [27,40,72].Theoretical work shows that this could indeed select for a specific direction [60].Interestingly, Ref. [73] notes that a few proteins e NEK6, SPIRAL2, CLASP e that affect the stability of microtubule ends reduce the sensitivity to local mechanical stress patterns, suggesting that plants actively reduce the risk of idiosyncratic morphological responses to local variations in wall mechanical stress [24,36,68,73].This is an interesting topic, worthy of further exploration using both single cell-and tissue-level simulations.
Interestingly, animal and in vitro systems show that, besides full severing, katanin can remove small numbers of tubulin subunits.These incomplete cuts subsequently heal with GTP-tubulin, forming temporary "rescue sites" [43,74].Detailed energetic calculations on single microtubules and in vitro experiments suggest that without microtubule-associated proteins, such GTP islands are the predominant source of rescues [1].If this "nibbling" of katanin is dependent on tensile stress on the microtubule, e.g., via an increased density of minor lattice defects, this could increase rescue frequencies for microtubules aligned with wall stresses.This idea currently is quite hypothetical.The alignment of microtubules with stress depends on katanin [24,61]; however, complete severing increases in the process [61] and, moreover, katanin severing critically amplifies microtubules with the new, opposing orientation in bluelight induced reorientation of the hypocotyl array [46,47,52,60].
Intuitively, any orientation-dependent bias on microtubule dynamics could orient the array, as long as the effect is strong enough and can overrule orienting effects like those from the cell geometry [10,58].Many biases have been tried: Mirabet et al. [48] use microtubule tips preferentially adjusting their growth direction towards a particular orientation (corresponding to maximal cortical tension on a similarly shaped pressurized cell).Schneider et al. [63] bias array orientation (for computational efficiency) by a transient bias on the orientation of the unbound fraction of microtubule nucleations.Analytical work shows that a difference in growth propensity can be selected for a desired orientation, whether modulated via growth speed, catastrophe rate or intrinsic nucleation rate [60].In a conceptually similar way, multiple studies have investigated patterning of the array via local modulation of the spontaneous catastrophe rate based on a hypothesized [38,39,63] or simulated [66] pattern of ROPs (and effectors).Note that with biases on dynamic instability, one has to carefully check whether microtubules remain in the bounded growth regime for all possible orientations to avoid artefacts [16,17].

Outlook
We have shown how theoretical concepts can increase our understanding of plant microtubule systems and help design informative and efficient simulation strategies (see also [17]).A big open question in the field is how the cortical array integrates different orienting factors (mechanical stress, geometry, specific proteins) under different conditions.A related exciting question is whether proteins like abovementioned NEK6, SPIRAL2 and CLASP tune the sensitivity of the array for specific cues, or affect the responsiveness of local array orientation in a more generic way.These questions are inherently quantitative, requiring both the increasingly quantitative nature of experiments [13,24,37,47,52], and matching array simulations e i.e., with sufficiently realistic mechanical behaviour (flexibility), katanin severing and microtubule-based nucleation e as these factors can all have a big impact on array sensitivity to local and global cues.
Future research may also require multi-level or multicomponent modelling strategies, like interaction with the actin cytoskeleton (for, e.g., division plane orientation [42] or the shaping of metaxylem pit boundaries [67]) or local modulation of microtubule dynamics in response to a dynamic cue (e.g., supracellular mechanical stresses in the cell wall or intracellular dynamics of ROPs and their effectors).Interaction with actin is conceptually most easily performed in Cytosim (already done in, e.g., Ref. [7]), whereas coupling with PDE3 models for mechanics or ROPs is computationally much easier in 2D models like CorticalSim or Tubulaton when restricted to a surface, and is already implemented in Ref. [66].

Declaration of competing interest
The authors declare no competing interests.

References
Papers of particular interest, published within the period of review, have been highlighted as:

Figure 1
Figure 1 [33], Tindemans et al.[70].Light blue shade in the right column means the aspect is (approximately) the same as in the analytical model of Hawkins et al.[33]  and Tindemans et al.[70].produced a much stronger sensitivity to cylindrical geometry than isotropic and global density-dependent nucleation[58].