1. Introduction {#s1}
===============

Through advances in molecular biology, much information has been gathered about plant hormone response networks (Santner et al., [@B63]) and the mathematical modeling of small regulatory networks has led to increased understanding of their function (Middleton et al., [@B44]). New reporters for hormone levels (Band et al., [@B5]; Brunoud et al., [@B13]), combined with image analysis tools (Pound et al., [@B56]), now permit the acquisition of data for individual cells within a whole organ. However, a current challenge is to understand how cell-scale processes influence behavior at organ scales; in particular, how variations in the mechanical properties of individual plant cell walls, affected by cell-wall-remodeling enzymes under the control of various hormones and regulatory networks, in turn control growth on the whole-organ scale. The complexity of this process, involving hydraulics, turgor pressure regulation, and a large number of biochemical changes in the constituents of the cell wall, necessitates a modeling approach in which these factors can be integrated. Multicellular models which incorporate cell mechanical properties are necessary to understand how the growth of individual cells is organized to drive organ-scale growth and development (Band et al., [@B4]). For instance, spatially one-dimensional models which treat the root as a line of cells have been used to examine the regulation of growth by hormones (Chavarría-Krauser and Schurr, [@B17]; Chavarría-Krauser, [@B14]; Muraro et al., [@B48]).

Vertex-based models (Weliky and Oster, [@B76]; Nagai and Honda, [@B49]) are an attractive framework for examining the effects of cell-scale processes on behavior at the whole-organ scale, as they explicitly treat each cell as an individual object. Symplastic growth (Priestley, [@B58]), where neighboring cell walls adhere and do not slide across each other, can be imposed automatically by the sharing of vertices and walls between adjacent cells, and many of the complications in simulations of animal tissues (caused by cells separating and sliding over each other) are considerably less significant for plant organs. However, standard implementations of vertex-based models can be inadequate to represent elongating organs such as roots because of the large aspect ratio of cells, which makes it difficult to prevent the transverse swelling of cells (normal to the primary axis of the organ).

Simple mass-spring-based models \[e.g., vertex-vertex systems (Smith et al., [@B65])\] and vertex-based models (Rudge and Haseloff, [@B60]; Dupuy et al., [@B27]; Jönsson et al., [@B39]; Smith et al., [@B66]; Dupuy et al., [@B26]; Hamant et al., [@B33]; Merks et al., [@B43]; Abera et al., [@B1]; Alim et al., [@B2]; Uyttewaal et al., [@B71]) have previously been used to examine plant organ growth. These models represent three-dimensional (3D) plant tissues as a collection of 2D polygons \[possibly on a 3D surface, e.g., in Smith et al. ([@B66]) and Hamant et al. ([@B33])\] in order to simplify computation. They capture mechanical processes such as cellulose microfibril reorientation and anisotropic viscous properties of the cell walls to only a limited extent, as they do not explicitly consider cell walls in the plane of the 2D simulation. None of these models has been applied to a growing primary root, for example. More sophisticated finite-element models for plant organs exist (Jönsson et al., [@B39]; Huang et al., [@B36]), but implementing cell division within these models, whilst possible, is a complex task. To address these problems, we present here a new 2D "vertex-element" model, combining vertex-based simulations with a finite-element discretization of in-plane walls, which is computationally efficient and capable of simulating problems on the scale of a whole organ, which in the present case is the elongating section of the root of *Arabidopis thaliana*. This method has similarities to the cell-level finite-element methods for animal cells developed by the Brodland group (Chen and Brodland, [@B18]; Brodland et al., [@B12]), but is tailored to the distinct mechanical properties of plant cells.

The composite structure of the primary plant cell wall endows it with complex properties that can be independently regulated by a variety of enzymatic and metabolic processes (Geitmann, [@B31]). Oriented cellulose microfibrils promote anisotropic elongation in the direction orthogonal to the fibrils (Baskin, [@B7]). The microfibrils are embedded in a pectin matrix and are crosslinked by xyloglucan polymers (Cosgrove and Jarvis, [@B20]). A variety of constitutive models have been proposed for the different wall components at different levels of organization. At the level of a whole cell or tissue, elongation is commonly modeled as a viscoplastic process, parameterized by a yield stress and an extensibility (Lockhart, [@B42]). Continuum models at the level of the wall have demonstrated how fiber reorientation (in addition to matrix stiffening) can suppress cell elongation (Dyson and Jensen, [@B29]) and have shown how crosslinks and matrix properties might independently contribute to yield and extensibility properties (Dyson et al., [@B28]). Macromolecular simulations (Kha et al., [@B40]; Yi and Puri, [@B78]) give insights into how detailed aspects of chemical structure influence mechanical properties. This hierarchy of complementary modeling approaches must be integrated in order to understand plant tissue growth fully. Here, we lay out a computational framework for connecting properties of individual cell walls to a multicellular tissue (i.e., spanning scales from cell walls to the organ level), using the Arabidopsis root as a template.

To illustrate the potential of our approach, to understand basic properties of the multicellular simulation tool (particularly its capacity to characterize the biomechanics of a growing tissue) and to discuss some of the potential pitfalls in its computational implementation, we show how differences in mechanical properties between cell files in an elongated organ can affect growth at the organ level. In particular, we increase the extensibility of a single cell file in a simple organ geometry, causing the growing organ to bend. We also explore the retardation of cell elongation though reorientation of microfibrils. Then, to show how mechanical properties can be readily controlled by hormone levels in such a framework, we consider a simple model in which cell mechanical properties (specifically, yield stress) are regulated by a substance (a growth inhibitor) which is produced in the root apex and undergoes passive diffusion between cell compartments. Cell division is also included in this simulation. The geometry used to initialize the simulation is based upon confocal imaging. We use the model to illustrate the impact of detailed cellular geometries on organ-scale morphology.

2. Materials and methods {#s2}
========================

2.1. Vertex-element model
-------------------------

Simulations are performed in a 2D plane; for the current problem this is a 2D longitudinal section through the midline of the root, as illustrated in Figure [1](#F1){ref-type="fig"}. The 2D model captures aspects of the 3D structure by distinguishing three classes of wall: those in the plane of the simulation; those perpendicular to the plane of the simulation and aligned with the direction of elongation ("axial" walls); and those perpendicular to the plane of the simulation and orthogonal to the direction of elongation ("cross" walls). Cross walls represent the approximately circular end plates of elongating cells, whilst the remaining wall types represent the curved walls parallel to the direction of elongation. In the 2D representation, cells occupy polygonal regions, the boundaries of which are described by a list of edges, each edge being a line joining two vertices. Adjacent cells share edges. Whilst the common edge is treated as a single entity for the representation of the organ geometry, for mechanical purposes we consider the physical cell walls associated with each of the two cells separately. This structure automatically ensures that the organ remains contiguous, as is the case during symplastic growth.

![**Axial and cross walls. (A)** Classification of cell walls perpendicular to the plane of the simulation as axial (red) and cross (blue); walls in the plane of the simulation are shaded gray. **(B)** Corresponding walls in the 3D organ.](fpls-04-00233-g0001){#F1}

The organ grows primarily through the motion of its vertices. Other changes, such as refining an edge into two smaller edges or the introduction of new walls during cell division, may (where necessary) be performed during the simulation. Topological transitions, which are important in animal epithelial tissues and soap bubble rafts, do not occur here because of the inability of neighboring plant cell walls to slide across each other or of an edge section to shrink to zero length, except during a small number of special events such as lateral root emergence, abscission, cell death etc.

### 2.1.1. Forces

Mechanical simulation using such a model involves consideration of the forces acting on each of the vertices. These are viscous and elastic forces from walls perpendicular and parallel to the plane of the simulation, and pressure forces acting on the boundaries of each cell. Owing to the small length scale of the cells, inertia is taken to be negligible, so these forces are always in balance. The forces acting on each vertex are calculated as functions both of the vertex positions and of the velocities, and equated to zero. This gives a system of differential equations relating the vertex positions and velocities, which are integrated numerically. The specific forms of the forces used are described in more detail below; in order for the system to be numerically well-behaved (more precisely, to be a system of differential, rather than differential-algebraic, equations) and to avoid transverse cell expansion, the mechanical properties of cell walls perpendicular to the plane of the simulation are considered. In order to prevent solid-body translation and rotation of the growing organ, the positions of vertices on the shootward (here corresponding to the right-hand) end of the organ are fixed.

***2.1.1.1. Turgor pressure***. Plant cells maintain their shape though an internal turgor pressure, which is generated by osmosis. Turgor pressure is taken to be constant and uniform in all cells; whilst experimental measurements (Pritchard et al., [@B59]) and theoretical considerations (Chavarría-Krauser et al., [@B15]) indicate that the turgor pressure is likely to decrease along the axis, in both cases this change was relatively small (\<10%). There is some evidence that it varies between cell types (and that for roots it is higher in cortical than epidermal cells) (Javot et al., [@B38]; Passioura and Boyer, [@B52]).

This pressure exerts an outwards force on each cell edge (per unit length normal to the plane of the simulation), normal to the edge and with magnitude proportional to its length, as shown in Figure [2A](#F2){ref-type="fig"}. Consider the edge *e* between the vertices with indices *i* and *j*, with positions **x**~*i*~ = (*x*~*i*~, *y*~*i*~) and **x**~*j*~ = (*x*~*j*~, *y*~*j*~). If these describe a wall of the cell with index *m* (and internal pressure *p*~*m*~) in the positive (anti-clockwise) direction, the force on the wall due to the pressure in cell *m* is $$\left. p_{m} \middle| \mathbf{x}_{j} - \mathbf{x}_{i} \middle| \mathbf{\hat{n}}_{ij} \equiv p_{m}\text{R}(\mathbf{x}_{j} - \mathbf{x}_{i}),\text{~~~~}\mathbf{\hat{n}}_{ij} \equiv \frac{(y_{j} - y_{i},x_{i} - x_{j})}{\left| \mathbf{x}_{j} - \mathbf{x}_{i} \right|}. \right.$$

![**Turgor pressure forces. (A)** Forces from the pressure *p*~*m*~ in cell *m* acting on the wall between vertices **x**~*i*~ and **x**~*j*~. **(B)** Total force acting on the vertex **x**~*i*~ from the pressure in cell *m*.](fpls-04-00233-g0002){#F2}

Here R(**x**) ≡ **x** × **k**, where **k** is the unit vector pointing out of the plane, corresponding to in-plane rotation by −π/2, and × denotes the vector product. The turgor force is taken to be equally distributed between the two vertices, giving the contributions $$\frac{p_{m}}{2}\text{R}(\mathbf{x}_{j} - \mathbf{x}_{i})$$ to the total forces acting at **x**~*i*~ and **x**~*j*~. The total pressure force on each vertex is the sum of the contributions from all the cells of which it is on the boundary, i.e., the force on vertex *i* can be expressed in the form $${\sum\limits_{m|i\mspace{2mu} \in B(m)}\frac{p_{m}}{2}}\text{R}\left( {\mathbf{x}_{m_{+}(i)} - \mathbf{x}_{m_{-}(i)}} \right)$$ where the sum runs over all cells *m* for which the vertex with index *i* is within its set of boundary vertices *B*(*m*), and *m*~−~(*i*) and *m*~+~(*i*) denote the indices of the previous and next vertices to *i* around cell *m* in an anticlockwise direction (see Figure [2B](#F2){ref-type="fig"}).

***2.1.1.2. Constitutive assumptions***. Turgor pressure forces are resisted by the cell walls of the tissue, which are assumed here to be under tension. The cell walls have a strongly anisotropic distribution of cellulose microfibrils embedded in a viscous matrix; these microfibrils increase the resistance to deformation of the cell wall in the fiber directions, thereby promoting elongation perpendicular to these directions. In simulating growth, over timescales of hours, axial cell elongation is determined by irreversible creep of the axial cell walls, which can be modeled using a viscoplastic constitutive law (Passioura and Fry, [@B53]; Cosgrove, [@B19]).

Wall stresses in the current configuration are obtained from a non-linear anisotropic constitutive law for the Cauchy stress resultant **σ**, which we assume takes the form $$\mathbf{\sigma} = \mathbf{\sigma}_{y} + \mathbf{\sigma}_{a} + \mathbf{\sigma}_{e}.$$

Here **σ**~*y*~ is a non-linear isotropic viscous contribution (incorporating yielding effects), **σ**~*a*~ is an additional anisotropic viscous term which includes the effects of the cellulose microfibrils and **σ**~*e*~ is the elastic contribution to the stress tensor. We attribute different combinations of these components to the three classes of walls in the simulation. We also assume that walls are sufficiently thin for their bending resistance to be negligible. We do not take explicit account of wall thickness but instead lump such factors into material parameters appearing below.

We model the non-linear isotropic component by $$\mathbf{\sigma}_{y} = 2\left( {\mu_{1} + \frac{\tau_{w}}{\varepsilon^{*}}\left( {1 - \exp\left( {- \frac{\text{ε}^{*}}{\varepsilon}} \right)} \right)} \right)\mathbf{Ε},$$ where $$\varepsilon \equiv \sqrt{\sum\limits_{i = 1}^{2}{\sum\limits_{j = 1}^{2}E_{ij}^{2}}}.$$

Here Ε is the rate-of-strain tensor, μ~1~ is an isotropic viscosity associated with the pectin matrix (plus embedded hemicellulose crosslinks), and τ~*w*~ and ε^\*^ are the yield stress and yield strain rate of the wall. This expression is adapted from the (one-dimensional) model of Dyson et al. ([@B28]), and has the advantage of being related directly to the properties of cross-links between cellulose microfibrils. If we assume τ~*w*~/ε^\*^ » μ~1~, the wall is then more extensible in the post-yielded state (ε » ε^\*^) and stiffer in the pre-yielded state (ε « ε^\*^). This non-linear model is similar to those used to regularize numerical simulations of visco-plastic (Bingham) fluids (Papanastasiou, [@B51]) and generalizes the traditional Lockhart equation (Lockhart, [@B42]). For simplicity, as in Dyson and Jensen ([@B29]), we assume that the cell wall remains of constant thickness (through the addition of new material to the cell wall) in the current configuration even as the wall is stretched. For simulations on a synthetic geometry (Figures [5](#F5){ref-type="fig"}--[7](#F7){ref-type="fig"} below), we set τ~*w*~ = 0.

We include the effect of microfibrils by considering a pair of microfibril directions, **a**~1~ and **a**~2~, that are embedded in the material of the cell wall and are advected with it. The resulting anisotropic viscous stresses are given by $$\mathbf{\sigma}_{a} = {\sum\limits_{k = \mspace{2mu} 1}^{2}\left\{ {\mu_{2}(\mathbf{a}_{k} \cdot \ \mathbf{Ε}\mathbf{a}_{k})\mathbf{a}_{k} \otimes \mathbf{a}_{k} + \mu_{3}(\mathbf{a}_{k} \otimes \ \mathbf{Ε}\mathbf{a}_{k} + \ \mathbf{Ε}\mathbf{a}_{k} \otimes \mathbf{a}_{k})} \right\}}$$ where the symbol ⊗ denotes the outer product of two vectors (**a** ⊗ **b** = **a b**^*T*^). The viscosity μ~2~ represents the additional resistance to extension in the fiber directions; **a**~*k*~ · Ε**a**~*k*~ gives the strain rate in the direction parallel to the fibers and **a**~*k*~ ⊗ **a**~*k*~ is a tensor which represents the direction of the fiber family. μ~3~ is an additional viscosity which resists shear of the cell wall parallel to the fiber direction. When μ~2~ is chosen to be much larger than μ~1~, the cell wall element has low extensibility in the direction of the fibers, and this leads to anisotropic cell expansion (Dyson and Jensen, [@B29]). We note that (Equation 7) is readily extended to include a distribution of fiber angles; for simplicity we restrict attention here to just two primary fiber directions.

While growth is a primarily viscous process, it is helpful to allow non-yielded walls in particular to sustain an elastic stress **σ**~*e*~. We adopt a simple linear relationship between stress and elongational strain with stiffness λ, as shown in more detail below.

***2.1.1.3. Cell wall discretization***. Cell walls in the plane of the simulation play an important role in preventing transverse swelling and in determining the rate of axial elongation. For such walls we make the constitutive assumption **σ** = **σ**~*a*~ + **σ**~*y*~. To simulate their properties, each cell wall is sub-divided into triangular elements (see Figure [3](#F3){ref-type="fig"}). We now describe how to calculate the forces exerted by a particular element on its vertices, without including the element label to simplify our notation. Such calculations are common in implementations of finite element methods---cf. Zienkiewicz and Taylor ([@B79]) or (Bonet and Wood, [@B11]). Most vertices will be associated with multiple triangular elements (see Figure [3](#F3){ref-type="fig"}), and so the total force acting on each vertex will be the sum of the contributions from all triangles for which it is a corner. Each element has a pair of fiber directions **A**~1~, **A**~2~ in the reference configuration. These fiber directions map to **a**~1~, **a**~2~ in the current configuration, which are calculated using $$\mathbf{a}_{k} = \frac{\text{F}\mathbf{A}_{k}}{\left| \text{F}\mathbf{A}_{k} \right|},\ k = 1,2,\text{~~F} = \text{~F}_{2}\text{F}_{1}^{- 1},$$ where **F**~1~ and **F**~2~ are the gradients of the maps from the unit triangle \[with vertices at (0, 0), (1, 0), (0, 1)\] to the reference and current configurations, respectively (see Figure [4A](#F4){ref-type="fig"}). These maps are given by linear interpolation, so that position vectors in the two configurations can be written as $$\mathbf{X}(\mathbf{\xi}) = {\sum\limits_{I = 1}^{3}\mathbf{X}_{I}}N_{I}(\mathbf{\xi}),\text{~~}\mathbf{x}(\mathbf{\xi}) = {\sum\limits_{I = 1}^{3}\mathbf{x}_{I}}N_{I}(\mathbf{\xi}),$$ where **ξ** = (ξ~1~, ξ~2~) parameterizes points in the unit triangle, the shape functions *N*~1~, *N*~2~, *N*~3~ are defined by $$N_{1} = 1 - \xi_{1} - \xi_{2},\text{~~}N_{2} = \xi_{1},\text{~~}N_{3} = \xi_{2}$$ and **x**~*I*~(*t*) (*I* = 1, 2, 3), **X**~*I*~(*t*) (*I* = 1, 2, 3) are the positions of the vertices of the triangles in the current and reference configurations (see Figure [4A](#F4){ref-type="fig"}). The gradients **F**~1~ and **F**~2~ are therefore given by $$\begin{array}{l}
{\left( F_{1} \right)_{\alpha B} = \frac{\partial X_{\alpha}}{\partial\xi_{B}} = {\sum\limits_{I = 1}^{3}X_{I\alpha}}N_{I,\mspace{2mu} B}\left( \xi \right),} \\
{\left( F_{2} \right)_{\alpha B} = \frac{\partial x_{\alpha}}{\partial\xi_{B}} = {\sum\limits_{I = 1}^{3}x_{I\alpha}}N_{I,B}\left( \xi \right),} \\
\end{array}$$ where $$N_{I,\mspace{2mu} B}(\mathbf{\xi}) \equiv \frac{\partial N_{I}}{\partial\xi_{B}}$$ and *X*~*I*α~ denotes the α component of **X**~*I*~ (α = 1, 2;*B* = 1, 2).

![**Triangulated cell walls in the plane of the simulation**. Edges of cells in the plane (axial or cross walls) are shown as heavy black lines; vertices are shown as black or red circles. Note that vertices may lie in the interior of cells (red circles); these are required to avoid using triangles with large aspect ratio in the initial configuration. As the tissue elongates, triangle aspect ratios increase, and either a more highly refined triangular mesh or periodic remeshing is required to accurately simulate organ bending (see Figure [6](#F6){ref-type="fig"}).](fpls-04-00233-g0003){#F3}

![**Coordinate systems and forces for triangular elements lying within walls in the plane of the simulation. (A)** Vertices of a triangular element have positions **X**~1~, **X**~2~ and **X**~3~ in the reference configuration, and **x**~1~, **x**~2~ and **x**~3~ in the current configuration. The gradient of the deformation **F** is calculated using **F** = **F**~2~ **F**^−1^~1~, where **F**~1~ and **F**~2~ are the gradients of the maps, defined by linear interpolation, from the unit triangle to the reference and current configurations, respectively. The elements each have a pair of microfibril directions, **A**~1~ and **A**~2~, in the reference configuration; these are embedded in the material of the cell wall, so are deformed with it to directions **a**~1~ and **a**~2~ in the current configuration. **(B)** Notation used in Equations (16) and (17) for the forces on and normals to the edges of each triangular element.](fpls-04-00233-g0004){#F4}

The velocity **v** of a material point is given by $$\mathbf{v} = {\sum\limits_{I = \mspace{2mu} 1}^{3}\mathbf{v}_{I}}N_{I}(\mathbf{\xi}),$$ where **v**~*I*~(*t*) = ∂**x**~*I*~/∂*t* (the Lagrangian derivative; *I* = 1, 2, 3) is the velocity of the *I*th vertex of the triangular element in which the point sits. The rate of strain tensor Ε appearing in Equation (5)--(7) is then $$\mathbf{Ε} = \frac{1}{2}\left( {\ \nabla\mathbf{v} + {(\nabla\mathbf{v})}^{T}} \right),$$ where the spatial velocity gradient ∇**v** is given by $${(\nabla\mathbf{v})}_{\alpha\beta} = {\sum\limits_{I = 1}^{3}{\sum\limits_{B = 1}^{2}v_{I\alpha}}}N_{I,\mspace{2mu} B}(\mathbf{\xi}){(F_{2}^{- 1})}_{B\beta},$$

We calculate the forces from a triangular element on each of its vertices in a similar manner to the pressure forces, as illustrated in Figure [4B](#F4){ref-type="fig"}. The forces **f**~12~, **f**~23~, **f**~31~ on each of the edges of the triangle are given by $$\begin{array}{l}
\left. \mathbf{f}_{12} = - \sigma \cdot \mathbf{n}_{12} \middle| \mathbf{x}_{2} - \mathbf{x}_{1} \middle| \text{~~}\mathbf{f}_{23} = - \sigma \cdot \mathbf{n}_{23} \middle| \mathbf{x}_{3} - \mathbf{x}_{2} \right| \\
\left. \mathbf{f}_{31} = - \sigma \cdot \mathbf{n}_{31} \middle| \mathbf{x}_{1} - \mathbf{x}_{3} \right| \\
\end{array}$$ where **σ** is given by Equation (4) with no elastic component (**σ**~*e*~ ≡ 0); the outward-pointing normals to each of the edges of the triangle are $$\begin{array}{l}
{\mathbf{n}_{12} = \frac{\text{R}(\mathbf{x}_{2} - \mathbf{x}_{1})}{\left| \mathbf{x}_{2} - \mathbf{x}_{1} \right|},\text{~~}\mathbf{n}_{23} = \frac{\text{R}(\mathbf{x}_{3} - \mathbf{x}_{2})}{\left| \mathbf{x}_{3} - \mathbf{x}_{2} \right|},} \\
{\mathbf{n}_{31} = \frac{\text{R}(\mathbf{x}_{1} - \mathbf{x}_{3})}{\left| \mathbf{x}_{1} - \mathbf{x}_{3} \right|}.} \\
\end{array}$$

Distributing the force on each wall equally to the vertices at either end, the forces **f**~1~, **f**~2~, **f**~3~ associated with stretching of cell walls on the three corners of the element are $$\mathbf{f}_{1} = \frac{1}{2}(\mathbf{f}_{12} + \mathbf{f}_{31}),\text{~~}\mathbf{f}_{2} = \frac{1}{2}(\mathbf{f}_{23} + \mathbf{f}_{12})\text{~~}\mathbf{f}_{3} = \frac{1}{2}(\mathbf{f}_{31} + \mathbf{f}_{23}).$$

The contributions from all elements are combined to give to total force from the walls in the plane on each vertex.

Cell walls perpendicular to the plane of the simulation are considered as viscoelastic; whilst the model is similar to that for cell walls in the simulation plane, for simplicity we do not include the effects of microfibril orientation, setting **σ**~*a*~ = 0 in Equation (4). Furthermore, because in our 2D simulation the walls perpendicular to the plane between adjacent vertices can only undergo expansion, we can restrict attention to a single elongational stress component in the wall direction. However, it is helpful in simulations to retain an elastic element in each "cross" wall, as these grow more slowly and are not clearly in a yielded state. For the purposes of mechanical calculations, we treat an edge "viewed from" each of the cells that border it as two distinct wall segments, and reference them by the pair (*m, e*), where *m* is the index of the cell and *e* the index of the edge. For symplastic growth, the two wall segments corresponding to a given edge have the same length and hence strain-rate.

From the constitutive law (Equation 4), the tensile stress resultant in the wall (*m, e*) is assumed to be $$\sigma_{(m,\mspace{2mu} e)} = 2\left( {\frac{\tau_{(m,e)}}{\varepsilon_{(m,\mspace{2mu} e)}^{*}}\left( {1 - \exp\left( {- \frac{\varepsilon_{(m,e)}^{*}}{\varepsilon_{e}}} \right)} \right) + \mu_{(m,e)}} \right)\varepsilon_{e} + \lambda_{(m,e)}\left( {\frac{l_{e}}{l_{0}{}_{(m,\mspace{2mu} e)}} - 1} \right)$$ where ε~*e*~ is the rate of strain of the edge, μ~(*m,\ e*)~ is the viscosity of the wall, ε^\*^~(*m,\ e*)~ is the yield strain-rate, τ~(*m,\ e*)~ is the yield stress, *l*~*e*~ is the current length of the edge, *l*~0(*m,\ e*)~ is the rest length of the wall and λ~(*m,\ e*)~ is the spring constant of the wall. The resulting contributions to the forces on the vertices are $$\pm \sigma_{(m,e)}\frac{\mathbf{x}_{j} - \mathbf{x}_{i}}{\left| \mathbf{x}_{j} - \mathbf{x}_{i} \right|},$$ where the positive and negative signs correspond to the forces on **x**~*i*~ and **x**~*j*~, respectively.

The non-linear model (Equation 19) reduces to a linear model when τ~(*m,\ e*)~ = 0; as for the in-plane walls, we use the linear model for simulations on artificial geometries (Figures [5](#F5){ref-type="fig"}--[7](#F7){ref-type="fig"}), and the non-linear model (τ~(*m,\ e*)~ \> 0) for simulations on realistic root geometries (Figure [8](#F8){ref-type="fig"}). For simulations on an artificial geometry, cross walls are considered to be viscoelastic (λ~(*m,\ e*)~ = λ~cross~ = const, μ~(*m,\ e*)~ = μ~cross~ = const in Equation (19)), but axial walls are treated as purely viscous (λ~(*m,\ e*)~ = 0, μ~(*m,\ e*)~ = μ~axial~); for the simulations of Figures [5](#F5){ref-type="fig"}--[7](#F7){ref-type="fig"}, cell wall properties are assumed to be constant in time, but we permit differences in axial cell wall extensibilities between different cells. In the simulation in Figure [8](#F8){ref-type="fig"} below, cross walls have constant and uniform properties (τ~(*m,\ e*)~ = τ~cross~ and ε^\*^~(*m,\ e*)~ = ε^\*^~cross~), axial walls have constant and uniform yield strain rate ε^\*^~(*m,\ e*)~ = ε^\*^~axial~, but the axial wall yield stress τ~(*m,\ e*)~ varies between cells and with time in a manner explained in section 2.1.2 below.

![**Viscous cell walls in the plane of the simulation avoid transverse swelling of cells**. The synthetic geometry **(A)** was used to initialize simulations, with results at *t* = 1 hr shown in **(B)** and **(C)**. In **(B)**, with large anisotropic wall viscosity, cells retain their rectangular shape and elongate in the *x*-direction, whereas in **(C)**, with no anisotropic viscosity and only small isotropic wall viscosity (μ~1~ = 0.002 MPa hr, μ~2~ = μ~3~ = 0 MPa hr), cells become rounded. Scale bars indicate 100 μm. Parameters are as listed in Table [1](#T1){ref-type="table"}, except that μ~axial~ = 10.0 MPa μm hr in **(C)** to compensate for loss of cell-wall viscosity in the plane of simulation. 7772 triangular elements were used for the cell walls in the plane of the simulation; these elements are visible in the magnified images.](fpls-04-00233-g0005){#F5}

![**Variations in cell mechanical properties lead to organ-scale bending. (B,C)** Simulation results at *t* = 2 hr with initial synthetic geometry **(A)**. The gray cells in the uppermost file of the organ in **(C)** have slightly lower isotropic viscosity (μ~1~ = 0.13 MPa hr) than the other cells in **(B)** and **(C)** (μ~1~ = 0.15 MPa hr). **(D,E)** The total length of the midline of the organ *L* and the angle θ made between the tip and the *x*-axis, respectively. Initially, both organs grow roughly exponentially; however, the bending rate may decrease after a while if triangular elements become highly elongated **(F)**, in which case linear interpolation over them gives a poor approximation to the strain rate of the cell walls. This computational artifact can be avoided either by generating a more refined initial mesh ("refine", **(G)**) or by periodically regenerating the triangulation of the tissue ("remesh", **(H)**), as shown in **(D)** and **(E)**. Crosses marked **(B,C)** in (**D, E**) correspond to solutions **(B,C)** at *t* = 2 hr. Parameters are as listed in Table [1](#T1){ref-type="table"}. Scale bars indicate 100 μm.](fpls-04-00233-g0006){#F6}

![**Cell wall microfibril reorientation can slow organ elongation. (B,C)** Simulation results after 3 hr with initial tissue geometry displayed in **(A)**. In **(B)**, microfibrils (indicated by red and green lines) are initially aligned with the *y*-axis, whilst those in **(C)** are oriented at angles of ± 0.2 to the *y*-axis. **(D)** shows the lengths of the organs, and **(E)** shows the mean angle ϕ between the microfibril orientations and the *y*-axis. The organ in **(B)** grows (roughly) exponentially, whilst the elongation rate of the organ in **(C)** slows as the microfibrils become more closely aligned with the direction of expansion. Triangulation of geometry has 7722 elements; results with refined initial mesh or periodic remeshing (not shown) are very similar to those shown here (\<1% relative difference). Markers on **(D,E)** show solutions (**B,C**) at *t* = 3 hr. Parameters are as listed in Table [1](#T1){ref-type="table"}. Scale bars indicate 100 μm.](fpls-04-00233-g0007){#F7}

### 2.1.2. Diffusible growth inhibitor

For simulations in a realistic geometry, it is necessary to specify the mechanical parameters of individual cell walls to capture observed patterns of growth. In order to demonstrate the capabilities of this method, we consider a growth regulator which inhibits growth, with concentration *b*~*m*~ (where *m* runs over the indices of all cells). This is a relatively simplistic model for the regulation of cell mechanical properties by hormones, but it provides a framework that is readily generalized. Chavarría-Krauser et al. ([@B15]) proposed a related model for the regulation of growth in a line of cells, where cytokinin produced in the cell at the tip of the root played similar role to the growth inhibitor here. Muraro et al. ([@B48]) also proposed a model for the regulation of growth in the root involving auxin and cytokinin signaling.

The concentration of the representative growth regulator is assumed constant within each cell and transported between cells only by diffusion through cell walls. The cellular concentrations *b*~*m*~ evolve according to $$\begin{array}{l}
{\frac{\text{d}b_{m}}{\text{d}t} = \frac{1}{A_{m}}{\sum\limits_{n \in \mathcal{B}_{m}}\left( {P_{b}S_{m,n}(b_{n} - b_{m})} \right)} - \lambda b_{m} + \alpha_{m},} \\
{\text{~~~}m = 1,\ldots,N_{c}} \\
\end{array}$$ where $\mathcal{B}_{m}$ is the set of indices of cells adjacent to the cell with index *m, A*~*m*~ is the area of the cell with index *m, P*~*b*~ is a permeability coefficient, *S*~*m,\ n*~ is the length of the wall between the two cells with indices *m* and *n*, α~*m*~ is the production rate of substance *b* in cell *m*, λ is the decay rate of substance *b* and *N*~*c*~ is the total number of cells. In our realistic root geometry, each cell has an additional attribute labeling its tissue type. We consider production of *b* to occur only in those cells which are part of the quiescent center (QC). The growth inhibitor concentration serves as a proxy for distance from the QC, and controls growth through modulation of the yield stress of walls in the plane of the simulation, τ~*w*~ (see Equation 7), and the yield stress of axial walls perpendicular to the plane of the simulation, τ~axial~ (see Equation 19). Specifically, we assume that $$\tau_{w} = \frac{\tau_{w}^{0}}{\left( {1 + \left( \frac{k_{b}}{b_{m}} \right)^{n_{b}}} \right)},\text{~~}\tau_{axial} = \frac{\tau_{axial}^{0}}{\left( {1 + \left( \frac{k_{b}}{b_{m}} \right)^{n_{b}}} \right)}.$$

The half-saturation constant, *k*~*b*~, is chosen to control the position at which rapid cell elongation begins and the Hill constant, *n*~*b*~, controls the sharpness of this transition.

### 2.1.3. Cell division

In order to further demonstrate the capabilities of this computational approach, we integrated a simple model for cell division within this mechanical framework. Both the timing of cell division (Dissmeyer et al., [@B22]) and the choice of cell division plane (Sahlin and Jönsson, [@B61]; Besson and Dumais, [@B10]; Alim et al., [@B2]) have been modeled previously. Here we employ a simple model (related to that of Chavarría-Krauser and Schurr, [@B17]) to specify the division time of cells. Each cell has an associated variable *t*~*m*~, which is initially taken to be a random number drawn from a uniform distribution on \[0, 1). This increases with time in all cells at a rate β. A cell will divide provided *t*~*m*~ \> 1 and the growth inhibitor concentration is sufficiently large (*b*~*m*~ \> *k*~*b*~). Following cell division, *t*~*m*~ for both daughter cells is set to zero. Such a model results in a rate of cell division which is roughly constant in the meristem and zero outside it, as observed by Beemster and Baskin ([@B9]). It also requires cells to wait a period 1/β inbetween successive divisions.

We choose to divide cells perpendicular to the direction of maximum strain rate (Hejnowicz, [@B34]; Nakielski, [@B50]). This direction is obtained by calculating the area-weighted average of the rate of strain tensor Ε (Equation 14) over the triangles representing the wall of the cell in the plane of the simulation. This is a symmetric tensor, and the direction of maximum growth is given by the eigenvector associated with the largest eigenvalue. The plane of division is chosen to be that which passes through the geometric center (centroid) of the cell, perpendicular to the direction of maximum growth.

When a cell divides, the polygon representing the walls of the cell perpendicular to the plane of the simulation is divided into two along the plane of division. This requires subdividing two edges of the parent cell into two; a new triangulation is generated for the cells which share these edges and the two daughter cells. This new triangular mesh is obtained from a constrained Delaunay triangulation of the region. This mesh may contain new vertices, and the positions of these in the reference configuration (**X**~*I*~) are calculated by barycentric interpolation. The microfibril directions in the reference configuration, **A**~1,\ 2~, are constant and uniform within each cell in the simulations of this paper, and the fiber directions for the daughter cells are set to be the same as those of the parent cell.

### 2.1.4. Computational implementation

The majority of the simulations was implemented in Python, using the OpenAlea (Pradal et al., [@B57]) simulation environment, which provides data structures for representing vertex-based tissues and routines to manuipulate them. Computationally intensive parts of the simulations (such as solution of mechanical equations) were implemented in C++. Constrained Delaunay triangulations of cell walls in the plane of the simulation were obtained using Triangle (Shewchuk, [@B64]). All simulations were performed on a standard desktop workstation (16 GiB RAM, 3.4 GHz quad core Intel processor).

An implicit system of differential equations for the vertex positions is given by equating the total forces on each vertex to zero, obtained by adding the pressure forces (Equation 3) from each cell, the forces (Equation 20) from each out-of-plane wall segment, and the in-plane forces (Equation 18) from each triangular element. This system (with c. 12,000 degrees of freedom for the simulations of Figure [8](#F8){ref-type="fig"}), is solved numerically using the backward Euler method or the differential-algebraic solver (IDA) from the SUNDIALS suite (Hindmarsh et al., [@B35]). Both schemes use the Newton--Raphson method for the solution of non-linear systems of equations, see e.g., Iserles ([@B37]). The Jacobian matrix is calculated numerically using finite differences, with the known sparsity pattern used to group terms such that multiple columns can be calculated simultaneously. The direct sparse solver UMFPACK (Davis and Duff, [@B21]) is used for linear algebra. IDA is used for the simulations of Figures [5](#F5){ref-type="fig"}--[7](#F7){ref-type="fig"}, as it uses an adaptive stepsize method to control the error in the solution. However, we use the backward Euler method for the simulations of Figure [8](#F8){ref-type="fig"}, as this proves more robust for simulations on complicated geometries.

In simulations on a realistic geometry, the equations for the growth inhibitor (Equation 21) form an additional system of ordinary differential equations which are coupled to the mechanical model. For simplicity of implementation, we use a split-timestep method, in which we simulate each of the chemical and mechanical components of the model in turn over a fixed timestep. The calculated concentrations of growth inhibitor in each cell are then considered to be constant for the simulation of the mechanical parts of the model over one timestep. Cell division is considered at the end of each timestep.

### 2.1.5. Synthetic geometry

In order to investigate basic mechanical features of the model, a rectangular "root" was generated, with *N*~*y*~ = 5 files of cells of length *L*~*x,\ tot*~ = 2000 μm in the *x*-direction and height *L*~*y*~ = 20 μm in the *y*-direction. This height is comparable to the diameter of an Arabidopsis cortical cell (Dolan et al., [@B23]), although it should be noted that the other cell files in the real root are of smaller diameters. Each file of cells was further subdivided in the *x*-direction (starting at the wall with smallest *x*-coordinate, corresponding to the rootward end) with each cell length drawn randomly from the uniform distribution on \[*L*~*x*~, 2*L*~*x*~\], where *L*~*x*~ = 50 μm.

### 2.1.6. Realistic geometry

Longitudinal sections of Arabidopsis roots, with cell walls highlighted by propidium iodide staining, were obtained using confocal microscopy. These images were traced by hand using a vector graphics package (Inkscape). The OpenAlea environment (Pradal et al., [@B57]) provides routines to import tissues from Scalable Vector Graphics files. The root was aligned such that the long axis of the root was parallel to the *x*-axis. Edges were classified as axial or cross according to the angle between the edge and the *x*-axis (those edges within π/4 being axial and the others cross).

2.2. Parameter estimates
------------------------

Turgor pressures within Arabidopis have been measured to be of the order of 0.3 MPa (e.g., see Javot et al., [@B38]). The remaining parameters in the model are difficult to measure experimentally; the viscosity of an elongating cell wall is an emergent property arising from the dynamic breaking of hemicellulose crosslinks and the viscosity of the pectin matrix. However, cell elongation rates within the root have been measured by Peters and Baskin ([@B55]) to be about 0.1 hr^−1^ near the QC, and 0.5 hr^−1^ in the center of the elongation zone, and this guides our choices of parameter values, as described below. In the absence of independent information, the main criterion for these choices was that they gave biologically sensible results. Note that the growth profile of roots has been measured to vary significantly between different seedlings (Walter et al., [@B75]); a number of other authors have also measured growth profiles in roots using modern techniques (Walter et al., [@B74]; van der Weele et al., [@B72]; Basu et al., [@B8]; Chavarría-Krauser et al., [@B16]).

For an isolated rectangular cell of length *L*(*t*) and uniform width *H*, with fibers oriented perpendicular to its primary axis, we have from Equation (7) and (19), neglecting yield, that $$Hp = \left( {2\mu_{1}H + 4\mu_{axial}} \right)\varepsilon,\text{~~}\varepsilon = \frac{1}{L}\frac{\text{d}L}{\text{d}t},$$ showing how in-plane and side walls together regulate elongation. Taking $$\mu_{1} \approx \frac{p}{4\varepsilon},\text{~~}\mu_{axial} \approx \frac{Hp}{8\varepsilon},$$ with ε = 0.5 hr^−1^, the axial cellular elongation rate will be of the correct order, and the walls parallel and perpendicular to the plane of the simulation will make similar contributions to the forces resisting elongation. Yield strain rates (ε^\*^, ε^\*^~axial~) are chosen such that all axial walls and cell walls in the plane of the simulation (in the meristem and elongation zones) are in a yielded state. Parameters for cross walls and the anisotropic components of the viscosity in the plane of the simulation (μ~2~ and μ~3~) were chosen to give minimal cell growth in the lateral direction. For the growth inhibitor, we took the permeability coefficient to be *P*~*b*~ = 2 × 10^3^ μm hr^−1^ (following the estimates for auxin from Kramer, [@B41]); the decay and production rates were chosen to give a plausible elongation rate profile (Walter et al., [@B75]; van der Weele et al., [@B72]; Peters and Baskin, [@B55]; Chavarría-Krauser et al., [@B16]) along the root axis. A full list of parameter values is given in Table [1](#T1){ref-type="table"}.

###### 

**Table of default parameter values**.

  **Symbol**     **Description**                                 **Default value**   **Units**
  -------------- ----------------------------------------------- ------------------- -----------
  *p*~*m*~       Pressure                                        0.3                 MPa
  μ~cross~       Cross wall viscosity                            0.05                MPa μm hr
  λ~cross~       Cross wall spring constant                      100                 MPa μm
  τ~cross~       Cross wall yield stress                         0                   MPa μm
  ε^\*^~cross~   Cross wall yield strain rate                    0.05                hr^−1^
  μ~axial~       Axial wall viscosity                            1.5                 MPa μm hr
  τ~axial~       Axial wall yield stress                         0                   MPa μm
  ε^\*^~axial~   Axial wall yield strain rate                    0.05                hr^−1^
  μ~1~           Wall isotropic viscosity                        0.15                MPa hr
  μ~2~           Wall anisotropic (extensional) viscosity        20.0                MPa hr
  μ~3~           Wall anisotropic (shear) viscosity              0                   MPa hr
  ε              Wall yield strain rate                          0.05                hr^−1^
  τ~*w*~         Wall yield stress                               0                   MPa
  *P*~*b*~       Wall permeability for growth inhibitor          2 × 10^3^           μm hr^−1^
  λ~*b*~         Growth inhibitor decay rate                     16                  hr^−1^
  α~*b*~         Growth inhibitor production rate in QC cells    2                   nM hr^−1^
  *k*~*b*~       Half saturation constant for growth inhibitor   1 × 10^−5^          nM
  *n*~*b*~       Hill coefficient for growth inhibitor action    2                   n/a
  τ^0^~axial~    Axial wall yield stress                         0.1                 MPa μm
  τ^0^~*w*~      Basal wall yield stress                         0.1                 MPa
  β              Cell division rate                              0.1                 hr^−1^
  *N*~*y*~       Number of cell files in synthetic geometry      5                   n/a
  *L*~*y*~       Radial size of cell files                       20                  μm
  *L*~*x*~       Minimum (and half maximum) axial cell length    50                  μm

3. Results {#s3}
==========

We initially (Figures [5](#F5){ref-type="fig"}--[7](#F7){ref-type="fig"}) consider simulations on a synthetic geometry, with constant mechanical properties for individual cells prescribed at the start of the simulation and not including growth regulation by diffusible substances. We use models with linear viscous terms (τ = 0, τ~*w*~ = 0), allowing us to explore the properties of the mechanical model in more detail. Later (Figure [8](#F8){ref-type="fig"}), we will consider simulations on a realistic geometry, and include growth controlled by a diffusible growth inhibitor and non-linear viscous terms in order to represent the yielding behavior of cell walls. Unless otherwise noted, all parameters are as listed in Table [1](#T1){ref-type="table"}.

Figure [5](#F5){ref-type="fig"} shows the importance of the mechanical properties of cell walls in the plane of the simulation in promoting elongation. We consider two simulations with identical synthetic geometries as their initial configuration (shown in Figure [5A](#F5){ref-type="fig"}). In Figure [5B](#F5){ref-type="fig"}, cell walls in the plane of the simulation have significant viscosity, both isotropic (μ~1~) and parallel to the microfibrils (μ~2~), which initially lie in the (transverse) *y*-direction. We find that the organ elongates preferentially in the (axial) *x*-direction, with cells maintaining a rectangular shape, demonstrating how fiber reinforcement allows anisotropic expansion. Note that this is a consequence of our choice of substantial viscosity μ~2~ in the microfibril direction. In Figure [5C](#F5){ref-type="fig"}, the viscosity of the cell walls in the plane is reduced to a very small value (μ~1~ = 0.002 MPa hr), but remains non-zero in order to regularize the numerical solution, and the anisotropic components of the viscosity (μ~2~, μ~3~) are set to zero. We increase the viscosity of axial walls (μ~axial~ = 10 MPa μm hr) in order to compensate for the loss of viscosity in the plane of the simulation. The cells bulge and become rounded, as there are no forces acting to hold together their axial sides. This behavior is similar to the transverse expansion of cortical cells arising when the elongation of endodermal cells is inhibited through blocking their response to gibberellic acid (Ubeda-Tomas et al., [@B70]).

Figure [6](#F6){ref-type="fig"} shows how cell-scale variations in mechanical properties can lead to organ-scale behavior. All the cells in the "uniform" case, shown in Figure [6B](#F6){ref-type="fig"}, have the same isotropic viscosity in the plane of the simulation (μ~1~ = 0.15 MPa hr), along with significant viscosity parallel to the microfibrils (μ~2~ = 20.0 MPa hr), which causes the organ to elongate in the *x*-direction. The "non-uniform" case, shown in Figure [6C](#F6){ref-type="fig"}, is identical to the uniform case, except that cells in the uppermost file (gray) have lower isotropic viscosity (μ~1~ = 0.13 MPa hr). In both cases, the initial configuration is as in Figure [6A](#F6){ref-type="fig"}. The graph of organ length against time in Figure [6D](#F6){ref-type="fig"} shows that the difference in the length of the midline of the organ between the uniform and non-uniform cases is small. However, the angle θ made between the left-hand end of the growing organ and the *x*-axis, shown in Figure [6E](#F6){ref-type="fig"}, indicates a clear bending response (at *t* = 2 hr) in the non-uniform case, as shown in Figure [6C](#F6){ref-type="fig"}. In the absence of adhesion to their neighbors, cells with lower viscosity would extend more rapidly; as the cells are tightly adherent to each other, the organ bends. The tip angle would be expected to be an increasing (and accelerating) function of time, but simulation results in Figure [6E](#F6){ref-type="fig"} do not always show this. This is because, for an organ undergoing bending, the finite element discretization provides a less accurate approximation to the strain rate tensor in the plane of the simulation as the aspect ratio of the triangular elements increases (during elongation of the organ), as shown in Figure [6F](#F6){ref-type="fig"}. In Figures [6G,H](#F6){ref-type="fig"} we illustrate two approaches to avoid this effect. In the first one ("refine"), we use a finer spatial discretization (77,653 triangular elements, compared with 7722 in the "uniform" and "non-uniform" cases). With the finer mesh, the tip angle is an accelerating function of time, but the simulation is considerably slower (c. 60 min compared with c. 3 min for the non-uniform case). In the second approach ("remesh"), we regenerate a new constrained Delaunay triangulation of the geometry every 5 timesteps (0.5 hr). This introduces new vertices within the wall of each cell; the positions of these vertices in the reference configuration (**X**) are calculated by barycentric interpolation. The tip angle in this case is in agreement with the results from simulations on the finer mesh, but this method is computationally more efficient (c. 8 min). In summary, care must be taken to avoid erroneous predictions associated with distortion of elements, but practical steps can be taken that do not make computation time excessive.

Figure [7](#F7){ref-type="fig"} illustrates the effect upon organ growth of anisotropic viscosity caused by oriented cellulose microfibrils. In particular, it shows how reorientation of microfibrils can slow organ elongation (cf. the multinet model of Green, [@B32]). In Figure [7B](#F7){ref-type="fig"}, the microfibrils are initially oriented transversely (parallel to the *y*-axis), whereas in Figure [7C](#F7){ref-type="fig"}, the microfibrils are initially oriented in a pair of directions making angles ± 0.2 radians to the *y*-axis. Again, the initial configuration is as in Figure [7A](#F7){ref-type="fig"}. As can be seen from Figures [7D,E](#F7){ref-type="fig"}, the organ with microfibrils initially parallel to the *y*-axis extends exponentially with time, and there is little reorientation of the microfibrils. In simulations with the microfibrils starting at a small angle to the *y*-axis, initially they do not increase the resistance to elongation significantly, as can be observed in Figure [7D](#F7){ref-type="fig"}. However, as the organ extends, the angle ϕ made by the microfibrils with the *y*-axis increases (as they are embedded in the material of the wall), as shown in Figure [7E](#F7){ref-type="fig"}. This increases their contribution to the viscosity in the direction of elongation, thereby slowing the growth of the organ. (The situation in Figure [7B](#F7){ref-type="fig"}, where microfibrils are perpendicular to the direction of elongation, is not stable to small numerical perturbations; slight asymmetries may cause the microfibril orientations to rotate.) This demonstrates how fiber reorientation alone is sufficient to suppress tissue growth, supporting predictions for individual cells by Dyson and Jensen ([@B29]). Note that, in an actual root, we expect fiber angles to be spatially non-uniform, being aligned more closely with the axis of elongation as we move further away from the root apex.

Figure [8](#F8){ref-type="fig"} shows that the computational model can be used to simulate growth on the scale of the root apex of an Arabidopsis seedling. In the present model, we simulate growth by cell expansion alone, ignoring cell division. The geometry used (Figure [8A](#F8){ref-type="fig"}) consists of the whole root tip and about half of the rapidly growing region (elongation zone), but not the region further away from the tip within which growth slows (the "growth terminating zone") (Verbelen et al., [@B73]). In Figure [8](#F8){ref-type="fig"}, levels of the diffusible growth inhibitor are high near the QC, but decrease along the main axis of the root. This modifies the yield stress of the cell walls \[according to Equation (22)\], and leads to cells undergoing a rapid transition between slowly growing and rapidly elongating states, as is found in experimental measurements of growing roots (Chavarría-Krauser, [@B14]; Peters and Baskin, [@B55]). (The simulations of Figure [8](#F8){ref-type="fig"} have 6265 triangles and 12,860 degrees of freedom; the simulation (Figures [8A--C](#F8){ref-type="fig"}) for 2 hr took less than 2 min.) The local cellular pattern near the root apex at *t* = 1 hr and *t* = 2 hr is shown in Figures [8D,E](#F8){ref-type="fig"}. Cell divisions can be observed in a number of cells, showing the capability of the model. However, it is clear that additional information about the regulation of growth and division in the region near the QC is necessary for the model to be used to study this region in detail.

![**Simulation of a growing Arabidopsis root apex. (B,C)** Simulation results at *t* = 1 hr and *t* = 2 hr with initial geometry **(A)**. Shading indicates cell growth inhibitor concentration (on a logarithmic scale); graphs below (**B,C**) show the expansion rate of the cells (hr^−1^), plotted at the *x*-position of the cell centroids. **(D,E)** The regions highlighted in gray in **(B,C)** in more detail. Triangulation of geometry has 8937 elements (see **(A)**). Simulation parameters are as in Table [1](#T1){ref-type="table"}, except μ~axial~ = 0.05 MPa μm hr, τ^0^~axial~ = 0.1 MPa μm, μ~1~ = 0.3 MPa hr, τ^0^ = 0.1 MPa for (**B,C**) and μ~axial~ = 0.5 MPa μm hr, τ^0^~axial~ = 0.1 MPa μm, μ~1~ = 0.2 MPa hr, τ^0^ = 0.1 MPa for **(D)**.](fpls-04-00233-g0008){#F8}

In order to obtain simulations which qualitatively agree with observations of growing roots (Figures [8A--C](#F8){ref-type="fig"}), it proved necessary to modify the simulation parameters, setting μ~axial~ = 0.05 MPa μm hr, τ^0^ = 0.1 MPa, τ^0^~axial~ = 0.1 MPa μm, μ~1~ = 0.3 MPa hr. These changes were made to reduce the importance of the walls perpendicular to the plane of the simulation. However, with greater axial wall viscosity (μ~axial~ = 0.5 MPa μm hr; the isotropic viscosity μ~1~ was decreased to 0.2 MPa hr to increase the overall expansion rate of the organ), it was found that the organ undergoes substantial bending, as shown in Figure [9A](#F9){ref-type="fig"}. This is a result of the cells on the lower side of the geometry being slightly narrower than those on the upper side; such differences are pronounced in 2D as cell widths depend upon the plane used to generate the 2D section. This effect is illustrated in Figure [9B](#F9){ref-type="fig"}; this simulation is identical to Figure [6B](#F6){ref-type="fig"}, except that the bottom file of cells is made to be slightly narrower (18 μm compared with 20 μm), and this results in the organ bending. Despite the limitations of the simulation, this result illustrates the sensitivity of the overall organ morphology to the shape of its constituent cells when one takes full account of the mechanical and kinematic constraints of symplastic growth.

![**Bending caused by small asymmetries. (A)** Simulation results at *t* = 2 hr with greater axial wall viscosity μ~axial~ than in Figure [8](#F8){ref-type="fig"}; bending is a consequence of cells on the lower side of the root being narrower than those on the upper side. **(B)** The result of a simulation on the artificial geometry of Figure [6A](#F6){ref-type="fig"} at *t* = 2 hr, where the lowest file of cells is slightly thinner (18 μm) than the other files (20 μm); cf. Figure [6C](#F6){ref-type="fig"}. Simulation parameters are as in Table [1](#T1){ref-type="table"}, except μ~axial~ = 0.5 MPa μm hr, τ^0^~axial~ = 0.1 MPa μm, μ~1~ = 0.2 MPa hr, τ^0^ = 0.1 MPa for **(A)**.](fpls-04-00233-g0009){#F9}

4. Discussion {#s4}
=============

We have presented a 2D vertex-element model for growing plant tissues that allows us to connect detailed mechanical properties of individual cell walls to the shape and growth rate of a whole organ. We showed for example that consideration of cell walls in the plane of the simulation is necessary when cells have a large aspect ratio, as is the case for cells in the elongation zone of *A. thaliana*, otherwise cells bulge and become rounded (Figure [5C](#F5){ref-type="fig"}). Note that angular springs could also be used to prevent this bulging, and are likely to be less computationally expensive; alternatively, beam elements could be used (Dupuy et al., [@B25]). However, our model directly includes the wall mechanical property that is believed to be responsible for the anisotropy of growth (Baskin, [@B7]). The constitutive law used for cell walls in the plane of the simulation incorporates anisotropy generated by cellulose microfibrils, and the yielding of cell walls is modeled using a non-linear viscous law, derived from a model of crosslink turnover (Dyson et al., [@B28]). While we have not here presented a complete parametric study, our model can nevertheless be used to assess the effect of enzymes that target different components of the cell wall on organ-level morphologies.

Simulations on a synthetic geometry showed how variations in cell scale properties could lead to organ-scale behaviors. These simulations also indicated the difficulty in simulating highly elongated tissues, as the elements in the plane of the simulation also become elongated and so less able to approximate the strain rate in the tissue. We demonstrated a solution to this computational difficulty through periodic remeshing of the tissue (Figure [6E](#F6){ref-type="fig"}). Our results highlight the importance of careful assessment of numerical errors that can easily arise in simulations of this complexity.

The effects of cellulose microfibrils on plant cell growth were investigated. Microfibrils rotate as the organ elongates, reducing the angle between their direction and the midline of the organ (Anderson et al., [@B3]). This increases the effective viscosity of the organ along the midline, thereby reducing its growth rate (Figure [7](#F7){ref-type="fig"}). Such observations support the hypothesis that microfibril reorientation may play a role alongside stiffening of the pectin matrix in determining when cells stop expanding as they leave the elongation zone.

The mechanical model was then applied to simulate a growing root apex (Figure [8](#F8){ref-type="fig"}), using a geometry acquired from confocal imaging. The yield stresses of cell walls were regulated by a diffusible growth inhibitor. Such a model is intentionally artificial. It is thought that a number of different hormones regulate growth in Arabadopsis roots (Ubeda-Tomas et al., [@B69]). Whilst the details of some hormone response networks have become sufficiently clear to permit the development of mathematical models (Middleton et al., [@B46], [@B45]; Dupeux et al., [@B24]), there are gaps in our understanding of how the levels of these hormones are controlled through synthesis and transport, and the downstream effect of these hormones upon cell wall biochemistry needs further investigation (Sanchez-Rodríguez et al., [@B62]). There is also only partial information about the effects of the numerous cell wall remodeling enzymes on cell wall mechanical properties, and we did not include spatial variation in the initial microfibril orientations. Moreover, our model does not at present include a number of processes, including the shedding of lateral root cap cells and forces imposed by the external medium, which are important for the maintenance of the structure and shape of the root apex. As the understanding of these effects becomes more well developed, these can be readily included in the current framework.

For simulations of the whole root apex, simulation parameters had to be chosen carefully to avoid bending because of small asymmetries in the geometry (Figures [8C,D](#F8){ref-type="fig"}). These asymmetries may be partly a consequence of the 2D representation of a 3D tissue, and also because of inaccuracies in the acquisition of the tissue geometry. By stiffening epidermal tissues and softening inner tissues (simulations not shown), we find that the tendency to grow straight is promoted. This may be a benefit of residual stress patterns in elongating organs, or it may compensate for the fact that a 2D projection of a 2D root underestimates the proportion of cells in peripheral tissues.

For roots to grow straight, in the absence of other external stimuli, may require additional mechanisms to coordinate growth, for instance sensing the stress or strain within individual cell walls. Alternatively, this bending may benefit the plant by encouraging the root to explore its local environment. These questions may need 3D models to be properly resolved. Nevertheless, it is notable that our method can be implemented to produce uniform elongation of a structure that has internal asymmetries in its architecture. It will be of significant interest to understand in more detail the interplay between geometric asymmetries and gradients in mechanical properties in processes such as gravitropic bending, extending existing models for root gravitropism (Barlow et al., [@B6]; Stočkus and Moore, [@B68]; Chavarría-Krauser, [@B14]; Moulia and Fournier, [@B47]) to include more cell-scale detail. Our results show that small differences in the mechanical properties in the outer tissue layers can lead to substantial curvature generation.

This vertex-element framework can be extended to include a number of biological processes which are important in regulating growth. Water transport between individual cells is regulated by the osmotic potential of the cells and the permeability of the cellular membranes. The latter can be controlled by aquaporins; certain aquaporins were observed to be specifically expressed in certain regions during the development of lateral root primordia (Péret et al., [@B54]), and a mathematical model which lumped tissues into different compartments elucidated the role of this process in controlling lateral root emergence. The current framework may be readily extended to incorporate non-uniform turgor pressure and water transport between cells, along with water transport in the apoplast and the supply of water to the xylem and from the phloem (Steudle and Peterson, [@B67]; Passioura and Boyer, [@B52]; Wiegers et al., [@B77]). The more detailed representation presented here will be of use in interpreting the importance of cell-scale spatial and temporal variations in aquaporin expression during lateral root emergence.

We have shown that, even in 2D, it is necessary to simulate both the edges and faces of individual cells in order to predict realistic organ-level growth. In extending the present model to 3D, there are a number of key difficulties which need to be overcome. Acquiring accurate 3D geometries for a whole root apex is currently beyond the capabilities of confocal imaging, although much progress has been made for other organs such as the shoot apical meristem (Fernandez et al., [@B30]). As the orientation of cell walls with respect to the direction of growth is critical, these geometric representations need to be of a very high quality. Finite-element discretization of cell walls in 3D is also more complicated, particularly if the cell walls are represented by shell elements: although many finite element packages implement this, cell division proves to be difficult in such a context. Such considerations illustrate the many challenges that remain in further developing such approaches and reinforce the current value of 2D formulations, both in their own right and in laying crucial foundations for further 3D studies.

Conflict of interest statement
------------------------------

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

This work was conducted in the Centre for Plant Integrative Biology, University of Nottingham, UK, which is jointly funded by the BBSRC/EPSRC (BB/D0196131) as part of their Systems Biology Initiative. John A. Fozard and Oliver E. Jensen were also supported by a BBSRC responsive mode grant (BB/J009717/1). John R. King would also like to acknowledge funding from the Royal Society and Wolfson Foundation.

[^1]: Edited by: Sanna Sevanto, Los Alamos National Laboratory, USA

[^2]: Reviewed by: Pierre Barbier De Reuille, University of Bern, Switzerland; Andres Chavarria-Krauser, Heidelberg University, Germany; Tim Rudge, University of Cambridge, UK

[^3]: This article was submitted to Frontiers in Plant Biophysics and Modeling, a specialty of Frontiers in Plant Science.
