A buckling instability and its influence on microtubule orientation in plant cells

TL;DR

This paper investigates how buckling instability influences microtubule orientation in plant cells, revealing that mechanical forces can cause microtubules to deviate from geodesic paths, affecting cell growth patterns.

Contribution

The study introduces an analytical model showing how elastic microtubules on cylindrical surfaces deflect due to curvature, impacting array formation and plant cell morphology.

Findings

Microtubules tend to deflect away from geodesics at high curvature.

Array formation is influenced by the balance between curvature and anchor density.

The buckling model generalizes classical Euler instability to biological structures.

Abstract

In growing plant cells, parallel ordering of microtubules (MTs) along the inner surface of the cell membrane influences the direction of cell expansion and thereby plant morphology. For correct expansion of organs that primarily grow by elongating, such as roots and stems, MTs must bend in the high-curvature direction along the cylindrically shaped cell membrane in order to form the required circumferential arrays. Computational studies, which have recapitulated the self-organization of these arrays, ignored MT mechanics and assumed MTs follow geodesics of the cell surface. Here, we show, through analysis of a derived Euler-Lagrange equation, that an elastic MT constrained to a cylindrical surface will deflect away from geodesics and toward low curvature directions to minimize bending energy. This occurs when the curvature of the cell surface is relatively high for a given anchor density. In the limit of infinite anchor density, MTs always follow geodesics. We compare our analytical predictions to measured curvatures and anchor densities and find that the regime in which cells are forming these cortical arrays straddles the region of parameter space in which arrays must form under the antagonistic influence of this mechanically induced deflection. Although this introduces a potential obstacle to forming circumferentially orientated arrays that needs to be accounted for in the models, it also raises the question of whether plants use this mechanical phenomenon to regulate array orientation. The model also constitutes an elegant generalization of the classical Euler-bucking instability along with an intrinsic unfolding of the associated pitchfork bifurcation.

Equations

$$\mathcal{K}[\phi ]={\int }_0^L((\cos \phi {)}^4+{\dot{\phi }}^2)ds,$$

$$\frac{1}{2}\ddot{\phi }+(\cos \phi {)}^3\sin \phi =0,\phi (0)={\phi }_0,\dot{\phi }(L)=0.$$

$$H(\phi ,p)=\frac{1}{4}{p}^2-(\cos \phi {)}^4.$$

$$\begin{array}{rl}\dot{\phi }& =\frac{p}{2},\\ \dot{p}& =-4(\cos \phi {)}^3\sin (\phi ).\end{array}$$

$$L={\int }_0^{{\phi }_f}\frac{1}{\sqrt{(\cos \phi {)}^4-(\cos {\phi }_f{)}^4}}d\phi$$

$$L_n=\frac{(1+2n)\pi }{2\sqrt{2}}.$$

$${\Delta }_L({\phi }_0,n)=\sum _{j=1}^n({\phi }_j(\ell )-{\phi }_j(0)).$$

$$|{\Delta }_L({\phi }_0,n)|\le \frac{L^2}{n}.$$

$$|{\Delta }_L({\phi }_0,n)|\le {\Delta }_L({\phi }_0,1),$$

$$\underset{n\to \infty }{\lim }{\Delta }_L({\phi }_0,n)=0.$$

$$C=\left\{\left[\begin{array}{c}\cos \theta \\ y\\ \sin \theta \end{array}\right]:\theta \in [0,2\pi ],y\in \mathbb{R}\right\},$$

$$\mathbf{x}(s)=\left[\begin{array}{c}\cos (\theta (s))\\ y(s)\\ \sin (\theta (s))\end{array}\right],$$

$$\mathbf{v}(s)=\dot{\mathbf{x}}(s)=\left[\begin{array}{c}-\sin (\theta (s))\\ 0\\ \cos (\theta (s))\end{array}\right]\dot{\theta }(s)+\left[\begin{array}{c}0\\ 1\\ 0\end{array}\right]\dot{y}(s)$$

$$\dot{\theta }(s)=\cos (\phi (s)),\dot{y}=\sin (\phi (s)).$$

$${\kappa }^2(s)=\dot{\theta }(s{)}^4+\ddot{\theta }(s{)}^2+\ddot{y}(s{)}^2=(\cos (\phi (s)){)}^4+\dot{\phi }(s{)}^2.$$

$$F_L({\phi }_0)=\{\phi \in C^1([0,L];\mathbb{R}):\phi (0)={\phi }_0\},$$

$$\mathcal{K}[\phi ]={\int }_0^L((\cos \phi {)}^4+{\dot{\phi }}^2)ds\phi \in F_L({\phi }_0).$$

$$\frac{d}{dr}\mathcal{K}[\phi$$

$$\displaystyle=2\dot{\varphi}(L)\xi(L)-2\int_{0}^{L}\big{(}2(\cos\varphi)^{3}\sin\varphi+\ddot{\varphi}\big{)}\xi\;\;ds,$$

$$\frac{1}{2}\ddot{\phi }+(\cos \phi {)}^3\sin \phi =0.$$

$$E=H(\phi (L),p(L))=-(\cos {\phi }_f{)}^4$$

$$p^2=4((\cos \phi {)}^4-(\cos {\phi }_f{)}^4),$$

$$L={\int }_0^{{\phi }_f}\frac{1}{\sqrt{(\cos \phi {)}^4-(\cos {\phi }_f{)}^4}}d\phi$$

$$\underset{{\phi }_f\to \frac{\pi }{2}}{\lim }L({\phi }_f)=\infty .$$

$$L({\phi }_f)\approx \frac{\pi }{2\sqrt{2}}+\frac{5\pi }{16\sqrt{2}}{\phi }_f^2.$$

$$\min L=\frac{\pi }{2\sqrt{2}}=:L_0$$

$$L_n=\frac{(1+2n)\pi }{2\sqrt{2}},n\in \mathbb{N},$$

$$\bar{\phi }(r)=\phi (Lr).$$

$$\bar{\cal K}_{L}[\bar{\varphi}]=\int_{0}^{1}\Big{(}L(\cos\bar{\varphi})^{4}+\frac{1}{L}(\bar{\varphi}^{\prime})^{2}\Big{)}\>dr.$$

$${\bar{\phi }}^{\prime \prime }+2L^2\bar{\phi }=0,\bar{\phi }(0)=0,{\bar{\phi }}^{\prime }(L)=0$$

$${\bar{\mathcal{K}}}_L[\bar{\phi }]>{\bar{\mathcal{K}}}_L[0]=L$$

Full text

{AMS}

92C80, 92C37, 92C10, 74L15

1 Introduction

Our understanding of the molecular basis of organismal development and morphology has been accelerating over the last couple decades under the influence of new techniques in molecular biology coupled with increasingly powerful theoretical and computational methods. One area that has benefited recently from these advances is plant development. Plant growth is heavily dependent on the manner in which cells deposit cellulose to form the cell wall because cell expansion occurs most easily in the direction perpendicular to the orientation of cellulose fibres [16]. Cellulose deposition is influenced, in turn, by ordered arrays of microtubules (MTs) that adhere to the inner surface of the cell membrane by anchoring to the cell wall [11]. Remarkably, unlike animal MT organization that is largely driven by centrosomes, for example, as seen in the mitotic spindle, these plant cortical arrays self-organize by a combination of cell geometry, distributed nucleation, and MT-MT interactions.

When a growing MT encounters an existing MT, the interaction can lead to one of three outcomes that together have been proposed to explain ordering: induced catastrophe, entrainment/zippering, or cross-over and subsequent severing [10]. Early computational and mathematical analysis demonstrated that models of cortical MTs on a two dimensional surface could self-organize into highly ordered arrays driven only by these MT-MT interactions [7, 3, 23, 12, 17, 21, 9]. Simulating such arrays on cylinders with inaccessible end-caps, the arrays align circumferentially as seen in growing plant cells [3, 12]. At the level of the whole cell, extensions of these models have been implemented in more complex geometries such as polyhedral domains [4] and more recently in finely resolved triangulations of actual cells [8]. Another recent model examined the dynamic relocalization of MTs from cytosol to cortex in abstract 3D geometries to assay the influence of cell shape on MT distribution and orientation, in the absence of MT-MT interactions [19]. These whole-cell models have helped clarify how array organization is globally coordinated across geometrically complicated domains and in some cases how MT-associated proteins can regulate that coordination [4, 8].

In all of these models, MT growth was assumed to follow geodesics of the cell surface and the influence of cell wall curvature deflecting MTs through bending in the tangent plane was ignored (in one stochastic model, MTs ought to follow geodesics on average [19]). Although some authors have given brief justification for this simplifying assumption [3, 23, 22], the issue has never been properly addressed in the context of this type of model.

The cortical-array organization of MTs had been studied in earlier work using in vitro experiments in fabricated cylindrical geometries containing free-growing MTs and, to understand these experiments, using thin-rod elasticity theory applied to theoretical rods in an analogous geometry [18]. Although the theory was accurate in its prediction of the MT configurations seen in the experiments, neither the experiments nor the theory were consistent with the configurations observed in analogous cylindrically shaped growing plant cells. As later pointed out, cortical MTs in plants are not free to move laterally on the cell membrane because they are anchored along their lengths [22]. In the latter article, the authors justified their own and others’ assumption that elasticity of MTs did not play a role in determining MT growth direction by invoking the stiffness of MTs and the frequency of anchoring. Essentially, they claimed that for unanchored MT tip lengths that are small compared to the local radii of curvature of the cell wall, MT deflection away from geodesics, induced by minimizing bending energy, would be negligible. While this claim is true, a rigorous treatment of this claim has never been given. Nor is it clear that plant cells are in the alleged scaling regime. The relevant radii of curvature of cell walls have been measured to be in the range of $7\ \mu m\pm 2.5\mu m$ (mean $\pm$ standard deviation) around longitudinal edges that MT arrays are seen to circumvent [4]. Distributions of unanchored MT tip lengths have also been measured and average around $3\ \mu m$ with tails extending up to $8\ \mu m$ [5, 2]. Although, on average, the reported unanchored lengths are smaller than the reported curvatures, the distributions are overlapping which makes predictions about the untested influence of mechanics on array orientation difficult to assess.

In this paper, we address this question directly: To what extent does minimization of bending energy allow the curvature of the cell wall to deflect the unanchored tips of cortical MTs away from geodesic paths?

2 Model equation derivation

We consider a single microtubule (MT) anchored at the inner surface of a cell membrane. In a plant cell, the cell wall is much more rigid than a MT so we assume the cell-wall shape is a fixed cylinder and the MT shape is determined by a process of bending-energy minimization. Anchoring proteins along the length of the MT keep it tight to the membrane but, as growth occurs at the plus end, a short length of the MT at that end, which we call the tip, is free to fluctuate thermally. While free, this tip can explore the surface of the membrane, and to some extent the cytosol, before getting pinned down by an anchor somewhere along the free length. We examine the problem of minimizing the bending energy of a single MT tip and consider implications for the organization of cortical MT arrays.

On a flat surface, the path of MT growth should be straight. When a surface has curvature, the MT bends. Previous work on MTs growing in such conditions assumed that MTs follow geodesics of the surface, as described above. However, when the surface has differing principal curvatures, we expect that the equilibrium configuration is one that deviates from a geodesic so as to turn the MT tip closer to the principal curvature direction with lower curvature (see Fig. 1 for illustration). More precisely, the shape of the MT between its plus end and the nearest anchor point, at the base of the tip, should be such that it minimizes its total bending energy given by the scalar curvature

[TABLE]

under appropriate boundary conditions. Here, the shape of the MT tip is described by the function $\varphi(s)$ which gives the angle formed between the circumferential direction of the cylinder and the tangent vector to the MT at a distance $s$ from the anchored end of the MT tip.

In Appendix A, we derive the functional (1) and the corresponding Euler-Lagrange equation:

[TABLE]

An abridged derivation of this equation has been published previously with different boundary conditions to describe free MTs in constrained geometries [18].

Hamilton’s equations

Analysing the problem in the phase plane is informative so we recast the Euler-Lagrange equation (2) as a Hamiltonian system. Using standard methods, we derive the Hamiltonian

[TABLE]

and the first order Hamiltonian system

[TABLE]

We illustrate the Hamiltonian vector field and some of its integral curves in Fig. 2. Solutions of the boundary value problem (BVP) $\varphi$ correspond to integral curves starting along the vertical line $\varphi=\varphi_{0}$ and ending on the horizontal line $p=0$ when $s=L$. This being a Hamiltonian system, solutions stay on a level curve of $H$.

3 Numerical methods

All numerical calculations were carried out using MATLAB®. To solve the EL equations (3) with boundary conditions $\varphi(0)=\varphi_{0}$ and $p(L)=0$ as required for Fig. 3, we integrated a family of initial conditions with $\varphi(0)=\varphi_{0}$ and $-2(\cos\varphi_{0})^{2}\leq p(0)\leq 2(\cos\varphi_{0})^{2}$ and selected the solution having $p(L)$ closest to zero. For Fig. 6, we integrated forward in $s$ starting with two left boundary conditions $(\varphi(0),p(0))=(\varphi_{0},0)$ and $(\varphi(0),p(0))=(\varphi_{0},2)$ using Matlab®’s fourth order Runge-Kutta scheme ode45 and used a bisection method on $p(0)$ to approximate the solution satisfying $p(L;p(0))=0$. To generate the bifurcation diagrams, we integrated equation (3) over a two-parameter family of $p(0)$ and $L$ values and approximated the diagram with a contour plot of $p(L)=0$.

4 Results

MTs growing parallel to the cylinder axis

The steady state solution $(\varphi(s),p(s))=(\pi/2,0)$ corresponds to a MT growing parallel to the axis of the cylinder and satisfies the right boundary condition $p(L)=0$ for all $L$. Since the curvature functional (1) is non-negative in general and vanishes in this steady state, this solution is a global minimizer of ${\cal K}$ and therefore a stable solution. This is, of course, clear from the physical point of view since any bending away from the axial direction increases both the cylindrical bending (the $\cos$ term in $\mathcal{K}[\varphi]$) and the bending of the MT in the tangent plane (the $\dot{\varphi}$ term). The symmetric solution for $\varphi_{0}=-\pi/2$ has similar properties but points in the opposite direction.

As suggested by the phase portrait, there is also a non-constant solution of the Euler-Lagrange equation associated with $\varphi_{0}=\pi/2$: it is the heteroclinic connection between $(\pi/2,0)$ and $(-\pi/2,0)$. However, for any finite $L$ it is not a solution of the minimization problem since it reaches $p=0$ only for $L\to\infty$: This heteroclinic corresponds to an infinitely long MT that comes from infinitely far down the cylinder, undergoes a half turn and heads back down the cylinder. The reverse heteroclinic is the mirror image solution on the cylinder.

A buckling instability for circumferentially growing MTs

Solutions having $\varphi_{0}=0$ have a richer structure. They correspond to a MT initially oriented circumferential to the cylinder, namely in the direction of maximal curvature on the cell. The constant $\varphi(s)=0,p(s)=0$ is a solution for any value of $L$. On the cylinder, this solution is an arc of a circle around the circumference of the cylinder.

Although intuitively a circumferential solution might seem to maximize the curvature, the situation is more subtle. In order for a MT tip, clamped in the circumferential direction at one end, to reduce its total bending energy by exploiting the flat axial direction along the cylinder, it must bend away from its initial direction. This induces an additional bending in the tangent plane on top of to the bending normal to the cell surface which is necessary to stay on that surface. Hence, it is natural to expect a transition in stability as a function of the length $L$, with the MT preferring the circumferential direction for small $L$, but bending away from it and towards the axial direction whenever $L$ is large enough to ensure a sufficiently long interval of MT heading in the flat direction.

In Appendix B, we show that for $L<L_{0}$ where $L_{0}=\pi/2\sqrt{2}\approx 1.1$, there is only the circumferential solution $\varphi(s)=0,\ p(s)=0$ to the EL equation with $\varphi_{0}=0$. Furthermore, we find that at $L=L_{0}$, there is a supercritical pitchfork bifurcation through which a symmetric pair of stable solutions appear and the zero solution goes unstable. Numerical approximations of the three solutions to the EL equation for a value of $L$ just above the bifurcation are shown in the top left and right panels of Fig. 3. The former shows the solutions in the phase plane and the latter shows them as rods on the cylinder.

We also derive an implicit expression for the angular direction of these solutions at the end of the tip, $\varphi_{f}(L)$:

[TABLE]

which can be inverted numerically to get the shape of the emerging branches. As $L$ increases, $\varphi(L)$ approaches $\pi/2$ meaning that this solution bends closer to the axial direction.

In addition to the pitchfork bifurcation at $L_{0}$, for each $n=1,2,3...$, there is a pitchfork bifurcation at

[TABLE]

The stability of the zero solution can be understood by analyzing the energy functional in a neighbourhood around it. We find that the zero solution is a minimum of the energy functional for $L<L_{0}$ and hence is stable there. At $L=L_{0}$, a single mode goes unstable and for all $L>L_{0}$ the zero solution remains unstable. This analysis is explained in Appendix B. To confirm that the branching solution that appears at $L_{0}$ is stable, we evaluate the energy functional numerically along the line defined by $\varphi_{u}(s)=A\sin(\sqrt{2}s)$ for several values of $L$ just above $L_{0}$. The mode $\varphi_{u}$ is the one that goes unstable at $L_{0}$ and the line spanned by it should be a close approximation of the heteroclinic connecting the zero solution to the branching solution for $L$ above but close to $L_{0}$. We find a pair of minima for some value of $\pm A$ that gradually move away from the zero solution with increasing $L$. This energy profile is shown in the bottom left panel of Fig. 3. Note the maximum at $A=0$ (unstable) and the two minima near $A=\pm 0.6$ (stable).

This result corresponds precisely to the intuitive explanation of the bifurcation given above. We now see that the tradeoff of tangential bending close to the anchor point in exchange for less bending farther away reaches a balance at $L_{0}$, with it being worthwhile above $L_{0}$ but not below.

The bifurcation diagram with stability denoted is outlined in Fig. 4.

Deflection of MTs growing close to circumferentially

A complete bifurcation analysis as carried out for the $\varphi_{0}=0$ case is not as simple for the case when $\varphi_{0}\neq 0$. However, we can still approximately substantiate the main message from the previous section: if MTs anchor frequently to the cell wall and the lengths of all segments are short then they grow nearly along geodesics. With $\varphi_{0}\neq 0$, symmetry is broken and the pitchfork unfolds into a stable steady state and a saddle-node pair. Fig. 5 shows the numerically calculated bifurcation diagram for the case of $\varphi_{0}=0.008$.

The biologically relevant feature is that for $L$ below $L_{0}$, $p(0)$, and hence $\varphi(L)$, is still close to zero for much of that interval. However, even for this small value of $\varphi_{0}$, below but close to $L_{0}$ there is already a significant deflection.

It has been commonly assumed in previous modelling work that higher frequency anchoring reduces the total deflection away from circumferential. The bifurcation diagram in Fig. 5 shows that the deflection approaches zero as $L$ approaches zero, a necessary but insufficient condition to justify this assumption. We must further ask whether these almost-zero deflections may add up to something significant over many anchored segments.

We compare the total deflection when the MT attaches only at $L$ with the total deflection when it attaches sequentially $n$ times to the surface with equal spacing $\ell=L/n$. We enforce a no-kink condition at each attachment point. For $1\leq j\leq n$, let $\varphi_{j}$ be the solution of the Euler-Lagrange equation on $[0,\ell]$ with boundary conditions $\dot{\varphi}_{j}(\ell)=0$ and $\varphi_{j}(0)=\varphi_{j-1}(\ell)$. Here we set $\varphi_{1}(0)=\varphi_{0}$. The total angular change is then given by the telescopic sum

[TABLE]

In this notation, $\Delta_{L}(\varphi_{0},1)$ corresponds to the scenario discussed above of a tip growing to length $L$ and being anchored only at the end without intermediate anchoring while $\Delta_{L}(\varphi_{0},n)$ corresponds to the scenario in which there are $n$ equally spaced anchors. In Appendix D, we prove that

[TABLE]

Let us explain how Eq. (5) supports the claim that high anchor frequency leads to near-geodesic MT growth for $\varphi_{0}\neq 0$. From (5) it is clear that for any $L$

[TABLE]

for $n$ large enough. In fact, we calculated $\Delta_{L}(\varphi_{0},n)$ numerically for $\varphi_{0}=0.1$ and several values of $n$ (see Fig. 6)

and found that $\Delta_{L}(\varphi_{0},n)$ appears to be a decreasing sequence in $n$ for any fixed $L$. Thus, more anchors means less deflection. Furthermore, Eq. (5) implies

[TABLE]

From this, we reach the same conclusion as in the case $\varphi_{0}=0$, although not as sharply. If $n$ is large enough, then the total angular change approaches zero for any fixed $L$. This conclusion brings us halfway to justifying the assumption that cortical MTs follow geodesics. The one issue to resolve, which we address in next section, is how long are the MT tip lengths in real cells.

Curvature and anchoring data suggest significant MT deflection

Recall that the length $L$ in our model is the nondimensional tip length with the background scale being the radius of the cylinder. To extract this from data, we need measurements of tip lengths and radii of curvature from actual cells. Tip lengths have been measured previously [5] and analyzed using a growth and anchoring model [2]. Typical radii of curvature have also been measured [4]. We used the mean ($7\mu$m), SEM ($0.2\mu$m), and sample size ($n=154$) of radii of curvature for longitudinal edges from root epidermal division zone cells from Fig. 3c of [4] to reconstruct an approximation of the original distribution in the data. Then, combining that with the actual distribution of tip lengths in wildtype cotyledon epidermal cells shown in Fig. 4 of [5] ($4.2\pm 1.6\mu$m), we generated a synthetic data set of nondimensional tip lengths by randomly sampling both distributions and calculating their ratio. Obviously, no MT whose tip length was measured in [5] ever crossed the longitudinal edges from which the curvatures were measured in [4] but these distributions represent reasonable scales for the required quantities.

The resulting distribution is shown as a histogram in Fig. 7.

As explained in the caption, although many MTs are predicted to stay close to geodesics, roughly one third of the “synthetic” MT tips are predicted by the model to undergo substantial deflection. This raises serious concerns for the assumption commonly made that cortical MTs follow geodesics, especially given that the deflection is away from the circumferential direction.

5 Discussion

The issue of how MT mechanics influences the organization of cortical arrays in plant cells has remained largely unexplored until now. Here, we found that MTs with sufficiently long unanchored lengths at their plus ends deflect from high to low curvature directions on the cortical surface of the cell. Given the important role of circumferential orientation of cortical MT arrays for cellulose deposition and hence plant growth and morphology, the density of anchors is thus a critical parameter that plant cells must control.

Specifically, we found that the free tip of a MT with its anchored end fixed exactly in the circumferential direction deflects significantly towards the axial direction when the tip length is larger than a critical length of $L_{0}\times R$ where $L_{0}=\pi/2\sqrt{2}\approx 1.1$ and $R$ is the radius of curvature in the circumferential direction. This deflection also occurs for MTs whose anchored ends are fixed close to the circumferential direction, even below the critical length. On the other hand, in the limit of short tip lengths, we proved that deflection goes to zero fast enough that, when added up over many anchored segments, the total deflection over a longer length of MT is still negligible.

Although we have shown that surface curvature can influence individual MT orientation, we have not addressed how this tendency plays out at the level of array organization and orientation tuning. Our analysis of previously published data on tip lengths and curvatures of cell walls indicates that a large fraction of MTs in the simulations of the ordering and orienting of cortical arrays should actually be subject to a deflecting force that is disruptive to circumferential organization. The effect of this force on array orientation needs to be tested in the full context of cortical MT-MT interactions and MT dynamic instability on complex cell geometries. Complicating the story of MT array orientation tuning, other factors have also been found to influence it. In particular, stress in the cell wall has been found to influence the orientation of MTs [14], possibly through the tension-dependent modulation of MT dynamic instability parameters which has been observed experimentally [13, 24] and conjectured to be the link between stress and MT orientation selection [15]. All these factors need to be tested in a comprehensive model of array formation and orientation to determine when, where and how each of these factors matter.

Interestingly, our analysis is independent of the flexural rigidity or Young’s modulus of the MT. In fact, a dimensional version of the energy functional given by Eq. (1) would have the flexural rigidity as a coefficient but it would fall out of the Euler-Lagrange equation due to homogeneity of the equation. The flexural rigidity would matter in two circumstances: (1) if we did not assume an infinite stiffness for the cell wall, the ratio of the elastic moduli of the cell wall and the MT would determine their relative deformations; (2) the probability that MT polymerization slows or stalls on a curved surface is dependent on curvature because a MT must deform to allow for the addition of subunits in the manner of an elastic ratchet [20]. Clearly, the ratio of elastic moduli is finite but large enough that cell wall deformation can be ignored. On the other hand, slowing or stalling of MT growth due to curvature is likely to occur. In the present work, we are focused on those MTs that continue to grow around corners and their speed is not relevant in our model, however, in a full-context model alluded to in the previous paragraph, this effect could be important.

Also relating to the flexural rigidity, one might ask how the energy scale of bending compares to the thermal energy scale. For a typical MT tip length ($4.2\ \mu$m) and radius of curvature ($7\ \mu$m), the bending energy is roughly $1000k_{B}T$ so we do not expect thermal fluctuations to make a large impact on the results described here.

Although polymerization speed is not a parameter in our model, it does provide for an interesting anchor-density regulatory mechanism. At slower polymerization speeds, anchor density is predicted to be higher even if the background rate of anchor formation is the same [2].

We have made several simplifying assumptions that may be worth revisiting for the problem we considered here and for those that we cannot consider given these simplifications. For example, in our model, new anchors attach at the very end of a MT tip. In reality, anchors could attach anywhere along the tip length, the impact of which on our quantitative predictions is not entirely clear. Cell geometry is another detail that requires further exploration. Although a cylinder is a reasonable first approximation of the shape of an elongating root cell, these cells tend to have polygonal rather than circular cross sections. Other plant cells, for example leaf pavement cells, have even more complicated geometries and microtubule orientation still appears to play a role in that complexity[1].

Our mathematical analysis generalizes the classical analysis of the buckling of a thin elastic rod under compression. The connection between these two is made more clear by considering an alternate physical scenario that leads to the same mathematical formalism discussed here. Suppose a thin elastic rod is clamped at one end to a flexible sheet and the sheet is slowly curled into a cylinder. If the clamped end is oriented perpendicular to the cylinder’s axis, the rod buckles through a pitchfork bifurcation when the product of the curvature of the sheet with the length of the rod is $L_{0}$, analogous to the classic rod buckling. In addition, this system provides an elegant unfolding of that pitchfork bifurcation when the orientation of the rod is not perpendicular to the cylinder axis.

In light of our findings, a few questions become a priority. Do simulations of cortical arrays still orient correctly under the deflecting influence of cell wall curvature? What is the actual distribution of MT tip lengths as MTs circumvent highly curved regions of the cell surface, especially near sharp edges where CLASP may be modulating deflection by promoting anchoring or where slowing of polymerization might be increasing anchor density? When MT arrays rearrange from a circumferential orientation during growth to more disordered or longitudinally oriented during non-growth phases, could the transition be driven by a regulated decrease in anchor density or increase in polymerization speed? How does cell geometry influence MT array organization in cells with more complicated geometry, for example leaf pavement cells, and is MT mechanics enough to explain the relationship?

Appendix A Derivation of the Euler-Lagrange equation

MTs on a cylindrical cell surface

We describe a plant cortical MT as a curve embedded within the cell surface, in the case described here, a cylinder. We make a few simplifying physical assumptions that allow us to treat each segment of MT between neighbouring anchor points independently. (1) We assume that anchors always attach the MT to the cell surface at the location of the plus end of a growing MT at the moment it anchors. (2) Upon attachment of an anchor at the plus end, the position and orientation of that point along the MT are henceforth fixed. This is equivalent to saying that anchors act as infinitely stiff positional and torsional springs. (3) As the MT grows beyond a new anchor, the shape of the new tip is determined by minimizing the total bending energy along the entire length of the MT tip under the assumption that it is on the cell surface at all points along the tip. In the calculations that follow, in light of assumption (2), we are always considering the shape of the MT tip and not the rest of the MT given that behind the foremost anchor the MT shape is fixed.

We consider curves constrained to lie on the (infinite) cylinder

[TABLE]

where units are chosen so that the cylinder has a (nondimensional) unit radius. All lengths in the model can thus be interpreted dimensionally by multiplying the nondimensional length by the actual radius of the cell in question. For a parametrized curve in $C$ given by

[TABLE]

the tangent vector ${\bf v}$ is given by

[TABLE]

so that $\|{\bf v}(s)\|^{2}=\dot{\theta}(s)^{2}+\dot{y}(s)^{2}$. Thus ${\bf x}$ is parametrized by arclength if $\dot{\theta}(s)^{2}+\dot{y}(s)^{2}=1$. Without loss of generality, we set $\theta(0)=y(0)=0$. The tangent vector and hence the curve is uniquely determined by a single function $\varphi(s)$ with the relations

[TABLE]

Note that these relations imply the identity $\dot{\varphi}(s)^{2}=\ddot{\theta}(s)^{2}+\ddot{y}(s)^{2}$. With this parametrization, the local curvature $\kappa^{2}(s)=\|\dot{\bf v}(s)\|^{2}$ is given by

[TABLE]

In order to keep the notation as simple as possible, we shall from now on drop the explicit $s$-dependence of all functions.

The Euler-Lagrange equation

The principle determining the MT’s shape is the minimization of the total curvature over paths of fixed total length $L$, the length of the MT tip. The initial orientation of the MT, $\varphi(0)=\varphi_{0}$ (equivalently described by $\dot{\theta}(0)$ and $\dot{y}(0)$), is determined by the orientation locked in by the anchor at the base of the MT tip. Let

[TABLE]

We are looking for the minimizers of the functional

[TABLE]

(as given in the main text as Eq. (1). To find the critical points of ${\cal K}$, we derive the associated Euler-Lagrange equation. We consider variations $\varphi(s)+r\xi(s)$ where $\xi\in F_{L}(0)$. The first variation in the direction of $\xi$ is

[TABLE]

where we have used $\xi(0)=0$ to eliminate one boundary term. If $\varphi$ is to be a critical point of $\cal K$, then this vanishes for all $\xi$, and in particular for all $\xi\in C^{\infty}_{c}((0,L))$ for which the first term vanishes. By the fundamental lemma of the calculus of variations, this implies that $\varphi$ is a solution of

[TABLE]

Once we know this, the first variation in the direction of $\xi$ reduces to $2\dot{\varphi}(L)\xi(L)$. This must vanish for all $\xi$ and hence $\dot{\varphi}(L)=0$.

In conclusion, $\varphi$ is a critical point for $\cal K$ if it solves (8) with boundary conditions $\varphi(0)=\varphi_{0}$ and $\dot{\varphi}(L)=0$.

Appendix B Bifurcation analysis and buckling instability

For $\varphi_{0}=0$, the constant $\varphi(s)=0,p(s)=0$ is a solution in $F_{L}(0)$ for any value of $L$. We seek to build a bifurcation diagram around this solution. First, we look for solutions having $\varphi_{0}=0$ but such that $p(0)>0$. If $\varphi_{f}=\varphi(L)$ is the value of the solution at the end point, the energy is given by

[TABLE]

where we used that $p(L)=0$ if a solution exists. In turn, this yields

[TABLE]

for all $s$. We now constrain ourselves to consider solutions that make a quarter rotation in the phase plane (starting on the line $\varphi=0$ and wrapping to the line $p=0$). For these, $L=\min\{s>0:p(s)=0\}$. Furthermore, $p(s)>0$ for all $s<L$ so that (9) and Hamilton’s equation $2\dot{\varphi}=p$ yield

[TABLE]

This defines a function $L(\varphi_{f})$ which gives the length of a MT tip required to ensure that the solution to equation (3) with boundary condition $\varphi(0)=0$ terminates at $\varphi(L)=\varphi_{f}$ when $p(L)=0$. We prove the following three facts about $L$ as a function of $\varphi_{f}$ in Appendix C:

$L$ is finite for $0\leq\varphi_{f}<\frac{\pi}{2}$. 2. 2.

$L$ is monotone increasing on $[0,\pi/2)$, with

[TABLE] 3. 3.

For $\varphi_{f}$ small,

[TABLE]

In particular,

[TABLE]

proving the following bifurcation picture. If $L\leq L_{0}$, the only solution with initial condition $\varphi_{0}=0$ is the steady state $\varphi=0$. If however $L>L_{0}$, there exists another branch of solutions having $p(0)>0$ and as $L$ increases, the direction at the end of the MT tip continuously deflects further from the circumferential direction to the axial direction.

Inverting $L(\varphi_{f})$ gives us a characterization of this second solution described by its final value $\varphi_{f}(L)$. Fig. 8 shows a numerical plot of $\varphi_{f}(L)$ and the quadratic approximation given above.

A third solution is characterized by the other inverse branch $-\varphi_{f}(L)$ and this corresponds to the symmetric solution in the phase plane starting at $p(0)<0$. Combining the zero solution with these two spatially varying solutions in a single diagram, we see that the system seems to undergo a pitchfork bifurcation at $L=\pi/(2\sqrt{2})$.

In fact, let us show that there is a sequence of solutions bifurcating from the zero solution in the direction $\sin(\sqrt{2}s)$ for

[TABLE]

For this, it is best to rescale the variable by setting $s=rL$ and define

[TABLE]

The curvature functional becomes

[TABLE]

While the Lagrangian now depends explicitly on $L$, the function space is $L$-independent, being $F=F_{1}(0)$ for all $L$. Now, the Euler-Lagrange equation can be cast as $G(L,\bar{\varphi})=0$, where $G$ acts on the Banach space $\mathbb{R}\times F$. By the inverse function theorem, $G(L,\cdot)$ is invertible in a neighbourhood of the point $\varphi=0$ as long as $G^{\prime}(L,0)$ is invertible. Hence, it suffices to check that the kernel of $G^{\prime}(L,0)$ is trivial. This is nothing else than showing that the linearized Euler-Lagrange equation

[TABLE]

has no solution in $F$. The kernel of the operator $\partial^{2}+2L^{2}$ is the intersection of $F$ with the two-dimensional space of solutions of the differential equation, namely $\mathrm{span}\{\sin(\sqrt{2}s),\cos(\sqrt{2}s)\}$. It is non-trivial if and only if $\sqrt{2}L=\pi/2+n\pi$ for any $n\in\mathbb{Z}$, in which case it is one-dimensional since the cosine never lies in $F$. Standard bifurcation arguments, see e.g. [6], then yield (i) neighbourhoods of $(L_{n},0)$ in which all solutions are given by two branches, namely the trivial one and a continuous branch bifurcating in the direction of $\sin(\sqrt{2}s)$ and intersecting the trivial one only at $(L_{n},0)$, and (ii) that there are no other bifurcation points along the zero solution.

Let us finally discuss the stability of the zero solution. Note that $\bar{\cal K}_{L}[0]=L$. We claim that the zero solution is stable for $L<L_{0}$, by which we mean that

[TABLE]

for all $L<L_{0}$ and all $\bar{\varphi}\neq 0$ in a neighbourhood of [math] in $F$. Let us first observe that Eq. (13) implies that for any fixed function $\bar{\varphi}\neq 0$, we have that $\bar{\cal K}_{L}[\bar{\varphi}]>L$ for $L$ small enough. Conversely, since $\int_{0}^{1}(\cos\bar{\varphi})^{4}dr<1$ for any $\bar{\varphi}\neq 0$ the reverse inequality $\bar{\cal K}_{L}[\bar{\varphi}]<L$ holds for $L$ large enough. In fact, since

[TABLE]

the function $L\mapsto\bar{\cal K}_{L}[\bar{\varphi}]-L$ is strictly decreasing. Hence, for each $\bar{\varphi}$ there is a unique $L(\bar{\varphi})$ for which $\bar{\cal K}_{L(\bar{\varphi})}[\bar{\varphi}]-L(\bar{\varphi})=0$. By the implicit function theorem, $L(\bar{\varphi})$ is a continuous function.

But then, the claim follows from the already established fact that there is no bifurcation point for $L<L_{0}$. Indeed, suppose that there is a path $[-1,0)\ni\sigma\mapsto\bar{\varphi}_{\sigma}$ such that $\bar{\varphi}_{\sigma}\to 0$ as $\sigma\to 0^{-}$ and such that $L(\bar{\varphi}_{\sigma})\to L_{*}<L_{0}$ as $\sigma\to 0^{-}$. Then by construction $\bar{\cal K}_{L(\bar{\varphi}_{\sigma})}(\bar{\varphi}_{\sigma})=L(\bar{\varphi}_{\sigma})\to L_{*}=\bar{\cal K}_{L_{*}}(0)$ and hence $\bar{\cal K}_{L_{*}}$ has a flat direction: there is a bifurcation at $(L_{*},0)$, which is a contradiction.

Appendix C Properties of $L(\varphi_{f})$

In order to simplify notation, we set $\varphi_{f}=x$. Setting $\varphi=xt$, we obtain

[TABLE]

First of all, if $0<x<\frac{\pi}{2}$, the denominator is ${\cal O}((t-1)^{1/2})$ and hence integrable at $t=1^{-}$, so that $L$ is finite. If $x=\frac{\pi}{2}$ however, the integrand is of order $(t-1)^{-2}$ and hence not integrable. In order to show monotonicity, we compute

[TABLE]

where $G(x)=(\cos x)^{4}+2x(\cos x)^{3}\sin x$. To show that $L^{\prime}$ is positive it suffices to show that $G(xt)$ is decreasing for $t\in[0,1]$ for any $x\in[0,\pi/2]$. The $t$ derivative of $G(xt)$ is $2x(\cos(xt))^{2}(4tx(\cos(xt))^{2}-\sin(xt)\cos(xt)-3xt)$. Thus it suffices to show that $4z(\cos z)^{2}-\sin z\cos z-3z\leq 0$ for $z\in[0,\pi/2]$. The $z$ derivative of this expression is $-2\sin z(4z\cos z+\sin z)$ which is clearly negative in this range.

Finally, the small $x$ behaviour of Eq. (11) is obtained from a Taylor expansion in x, with $t$ a parameter. Start with the Taylor expansion $(\cos x)^{4}=\sum_{n=0}^{\infty}a_{n}x^{n}$, where $a_{0}=1$, $a_{2}=-2$, $a_{4}=5/3$ and $a_{n}=0$ for $n$ odd. Using this we may write

[TABLE]

The coefficients $a_{n}$ arising from the expansion of an entire function are exponentially decreasing in $n$ so it is not hard to see that $\sum_{n=6}^{\infty}a_{n}\dfrac{t^{n}-1}{1-t^{2}}x^{n-2}=O(x^{4})$, uniformly for $t\in[0,1]$. Raising both sides to the power $-1/2$ and using the expansion $(1+x)^{-1/2}=1-(1/2)x+O(x^{2})$ leads to

[TABLE]

The result now follows from integrating with respect to $t$.

Appendix D Proof of a bound on $\Delta_{L}$

Here we prove inequality (5) which stated that $|\Delta_{L}(\varphi_{0},n)|\leq\frac{L^{2}}{n}$. To show this, we note that the Euler-Lagrange equation (3) immediately implies that $|\dot{p}|\leq 2$ and hence

[TABLE]

where we used that $p(\ell)=0$. By the dynamical equation again, we thus obtain the uniform bound

[TABLE]

Hence, the total angular change is given by

[TABLE]

where we recalled that $\ell=L/n$. This establishes the bound we had set out to show.

Acknowledgments

The authors would like to thank Geoff Wasteneys for useful feedback on the manuscript. This work was supported by the National Science and Engineering Research Council of Canada.