## Abstract

We describe a method for the simulation of the growth of elongated plant organs, such as seedling roots. By combining a midline representation of the organ on a tissue scale and a vertex-based representation on the cell scale, we obtain a multiscale method, which is able to both simulate organ growth and incorporate cell-scale processes. Equations for the evolution of the midline are obtained, which depend on the cell-wall properties of individual cells through appropriate averages over the vertex-based representation. The evolution of the organ midline is used to deform the cellular-scale representation. This permits the investigation of the regulation of organ growth through the cell-scale transport of the plant hormone auxin. The utility of this method is demonstrated in simulating the early stages of the response of a root to gravity, using a vertex-based template acquired from confocal imaging. Asymmetries in the concentrations of auxin between the upper and lower sides of the root lead to bending of the root midline, reflecting a gravitropic response.

## 1. Introduction

Growth is one of the main responses plants possess with which they adapt to their environment. To extract resources from their surroundings efficiently, plants change their growth patterns in response to external stimuli such as light, touch and the direction of gravity. As a result, the simulation of many processes in plant development requires an integrative model for organ growth. Plant organs usually have a multicellular structure; for example, the primary root of *Arabidopsis thaliana* contains multiple concentric cylindrical layers of cells. As adjacent cells are tightly bound to each other, the overall growth of the organ must be continuous (symplastic). However, different tissue layers have different responses to plant hormones such as auxin and giberellic acid; this facilitates integrated responses to simultaneous environmental and intrinsic stimuli [1].

The rate at which the plant hormone auxin accumulates in a cell is regulated by auxin efflux carriers, such as those of the PIN-FORMED (PIN) auxin efflux carrier family; these have an asymmetric, polar distribution on the plasma membrane of the cell that leads to directed auxin transport on a tissue scale. Simulations of auxin transport in the tip of the primary root [2–5], particularly in the complex region containing the columella cells and lateral root cap, require a geometrical representation that accurately identifies individual cells and their spatial relationship with other cells. Individual cell geometries must be represented with a sufficient level of detail to (at least) permit the calculation of the volumes (areas in two dimensions) of cells and the surface areas (lengths) of the interfaces between neighbouring cells. Early models for auxin transport in the root tip [4,6] used idealized templates to represent the geometries of individual cells. With recent advances in imaging and image analysis, it has become possible to generate templates for simulations from actual cellular geometries, both in two [2] and three dimensions [7,8].

Current mechanical models for plant organs generally fall into one of two categories: those that treat the organ as a continuum, and those that consider individual cells. Treating the organ as a continuum has computational and conceptual advantages, but fails to capture the precise multicellular structure. A range of models that consider individual cells have been developed for plant organs. These include, but are not limited to, cellular Potts, vertex [9,10], vertex-element [11] or finite-element [7,12–14] formulations. Those models which explicitly represent individual cell walls (i.e. vertex-element and finite-element models) can be sensitive to the precise positioning of cell walls [11], which is difficult to obtain with sufficient accuracy for a whole plant organ from microscopy images.

Here, we present a hybrid model, which resolves the geometries of individual cells through a vertex-based representation, but uses a coarser-scale midline representation to track organ growth. The evolution of the organ in its midline representation depends upon the properties of individual cells through averaged constitutive relationships, thereby allowing cell-scale changes in cell-wall mechanical properties, e.g. caused by changes in the concentration of the plant hormone auxin in a particular tissue layer, to affect growth on the whole organ scale. In turn, growth and deformation of the organ midline affects the cell-scale geometry, changing the lengths of cell walls and volumes of cells, thereby modifying hormone transport. It also changes the distance of cells from the organ tip, which in the model used here modifies the mechanical properties of cell walls. This serves as a simple, but concrete, example of a multiresolution method for multicellular organs. The effectiveness of this approach is illustrated by its application to the early stages of the bending of a growing root in response to gravity.

The model approach is in many ways similar to those adopted in implementing multiscale homogenization methods, such as classical techniques [15] and methods such as representative volume elements [16]. Such methods typically make the assumption that the microscale problem is spatially periodic, or can be represented as a realization of a random process, and do not explicitly track all the changes to the microscale structure. While plant tissue is often highly organized, the specific geometries that we wish to use do not satisfy these assumptions. Numerical methods based upon representative volume elements, in which mechanical properties of the material are obtained from simulations of small regions around the quadrature points of a coarse finite-element discretization, have been applied to plant organs [17,18]. However, this technique requires numerical simulation of the microscale model.

Moreover, unlike classical homogenization problems in continuum mechanics, the microscale problem here is not necessarily well characterized; we have an incomplete understanding of the material properties of cell walls, and it is challenging to measure accurately the geometries of individual cells in a multicellular tissue. As a result, the principal aim is not to solve the microscale problem accurately (as may be the case for homogenization or multigrid methods, for example), but instead to generate effective macroscale equations which capture the qualitative and, to an appropriate degree, the quantitative behaviour of the organ. Our approach exploits the fact that, for symplastic growth, strain fields typically vary smoothly over the tissue scale. By contrast, stress fields can share the heterogeneity of the microstructure [19], supporting integration of stress up to the tissue level.

The method used here is loosely based upon the Cauchy–Born rule for the mechanical properties of periodic structures: while the cell-scale structure here is not periodic, we assume that the deformation on the cell-scale is a smooth interpolation of the deformation on the tissue scale. The Cauchy–Born rule has been shown to have deficiencies when applied to models for animal tissues [18,20,21], largely because of the potential for lattice-scale deformations. However, plant tissues are highly ordered and grow in a symplastic fashion, which we expect to mitigate such deficiencies. Related methods have been developed in computer graphics, where the simulation of the mechanical behaviour of objects with highly complex (and sometimes topologically defective) geometries is of interest. One approach, e.g. used in [22], is to generate a simpler, coarse geometry; the mechanical behaviour of this coarse geometry is simulated and, through suitable interpolation, used to deform the detailed geometry.

The paper is organized as follows. Section 2 briefly describes the cellular-scale (vertex-based) representation and details the midline approach and computational implementation. Section 3 describes illustrative numerical results, and we summarize the advantages and limitations of the results in §4. Details of the auxin transport model, the coordinate transformations used for the midline model and the precise form of the growth profile used are described in the appendices.

## 2. Material and methods

Here, in order to describe the method in a simple context, we consider a two-dimensional model of a plant organ, where all points under consideration lie in a plane; therefore, individual cells occupy two-dimensional regions, and the organ midline is a curve lying in this two-dimensional plane.

### 2.1. Vertex-based organ representation

We use a vertex-based representation of the cellular-scale geometry of the organ (figure 1), as has previously been employed in a large number of other simulations (e.g. [23]). Cells occupy polygonal regions, described by the positions of their vertices; neighbouring cells have common vertices, which enforce the symplastic nature of organ growth. In the simulations of this paper, the implementation of a vertex-based tissue from the OpenAlea software platform [24] was used. Using such a representation, it is straightforward to calculate geometrical quantities such as cell areas and the lengths of the interfaces between neighbouring cells, which are needed for models of auxin transport and other cell-scale processes.

While the method described here can be used without triangulating the polygonal regions occupied by each cell, for the purposes of the simulations described below it proves very useful to do so. One reason for this is that the growth of organs varies over length scales that are comparable with the size of cells [25], which suggests that it is likely that the mechanical properties of cell walls are not uniform on each cell. Another reason is that it simplifies the calculation of the integrals appearing in the discretized versions of the midline equations. Triangularization of the vertex-based tissue representation is perfomed using *gmsh* [26], which uses a Delaunay algorithm. Additional vertices are inserted within cell walls and in the interior of cells to enforce a maximum triangle area.

### 2.2. Midline organ representation

Here, we outline the mathematical details of our midline description. A more general kinematic framework for three-dimensional deformations of elongating rods can be found in [27], for example. Compared with that study, we consider the midline of a two-dimensional object (rather than a three-dimensional rod) restricted to deform in a plane, and as we consider a viscous (rather than elastic) constitutive model we will be concerned with the strain rates (rather than strains) experienced by the material.

The geometry is described on a larger scale by a curved line through the midline of the organ, as can be seen in figure 2. For the purposes of simulation, the midline of the organ is parametrized using a Lagrangian representation, i.e. 2.1where is constant following a material point; corresponds to a fixed point on the organ (the base) and corresponds to the tip or apex of the organ. Using such a representation, common measures of organ geometry can be easily calculated; distance along the midline is given by the integral of the stretch 2.2so the total length of the organ is given by 2.3

Note also that . Some plant organs, such as roots, grow mostly in a region near the tip, so for the purposes of specifying growth profiles we also introduce the distance from the tip . The local unit tangent, ** t**, and normal,

**, to the midline are given by 2.4where**

*n***is a unit vector in the positive**

*k**z*-direction, normal to the plane of interest. The angle , between the positive

*x*-axis and the tangent to the midline, is chosen such that 2.5

One key measurement of the local shape of the organ is its curvature, defined by
2.6In the simulations, we will fix the position and angle at the end of the midline where , so and for all *t*. The midline path is then fully specified, and can be found from integrating (2.4)–(2.6) with respect to .

A number of measures have been proposed to describe the local kinematics of growth. The (Lagrangian) velocities of material points on the midline (relative to ) are given by 2.7The speed of a point along the midline is similarly given by 2.8We also define the spin (angular velocity) of a point on the midline to be 2.9The most common measure for the rate of growth along the organ axis is the relative elemental growth rate or relative elongation rate (RER), given by 2.10

Here, we choose to describe bending of the root midline through the rate of change of curvature (following a material point), .

#### 2.2.1. Relating midline and Cartesian coordinates

The midline representation describes the locations and trajectories of material points lying on the midline of the organ. By using curvilinear coordinates defined relative to the midline, we may extend this description to the whole of the organ in the form 2.11

For each point ** x** in the organ, we find the nearest point

**on the midline; is then the midline coordinate of**

*r***, and is the signed distance from the point**

*r***to the midline at**

*x***. Assuming that the tissue does not undergo shear deformation, these coordinates are orthogonal, i.e. everywhere, and uniquely defined in a local neighbourhood of the midline (provided is less than the maximum value of ), assuming that the midline does not turn through more than 180° (i.e. does not double back upon itself). The vertices in the vertex-based representation will be assumed to remain at fixed positions in the midline coordinate system.**

*r*As described in more detail in appendix A, the rate of strain tensor is given by 2.12i.e. it is diagonal in these curvilinear coordinates, with zero strain rate normal to the midline, and strain rate parallel to the midline, where this is given by 2.13

### 2.3. Viscous and elastic growth models

Here, we will consider plant organs over sufficiently long time scales and sufficiently short length scales that the effects of inertia are negligible. As a consequence, each material element in the tissue is in mechanical quasi-equilibrium. This will require a model for the constitutive law, which relates the deformation in the tissue to material stresses.

Multicellular models, which explicitly represent individual cell walls, have tended to use one of two main classes of constitutive laws: elastic growth models and viscous growth models.

Elastic growth models treat the cell wall on short time scales as an elastic material. The strain generated by cell-wall deformation relaxes over longer time scales through growth. Various different models for strain relaxation have been proposed [12,13,28].

Viscous growth models treat the cell wall as a viscous material, with cell-wall stresses depending on the strain rate in the wall rather than the strain [29,30]. Such a model is appropriate to capture growth processes, which occur over time scales longer than those on which elastic stresses relax, and has the advantages of being computationally and conceptually simpler than elastic growth models, as the viscous constitutive law in some respects (see below) combines the elastic constitutive law and the growth law. However, viscous models do have the weaknesses that they are not appropriate to study deformations occurring over short time scales (such as indentation experiments).

#### 2.3.1. Relationship between viscous and elastic growth models

A number of different elastic growth models are possible. The simplest possible one-dimensional model consists of a linearly elastic spring, with undeformed length that grows at a rate proportional to the strain (or stress) in the spring, for which
2.14where *l* is the length of the spring and the spring constant. We assume that *l* varies over time scales much longer than the relaxation time scale , which we take to be small.

The second part of (2.14) gives us that, to leading order in , . To next order in , this equation becomes . Substituting these approximations into the first part of equation (2.14) gives 2.15so the forces in the spring can be seen to be approximated by a viscous element (or dashpot) with viscosity , provided that the time scale of deformation is much longer than the relaxation time for viscous stresses. In this sense, viscoelastic growth models of the form (2.14) can be approximated by viscous growth models over long time scales. In modelling plant tissues, it is common for viscous relationships such as (2.15) to be supplemented with a threshold stress (or strain) above which the material yields irreversibly [31], but we do not incorporate this feature here.

### 2.4. Midline evolution equations

The midline evolution equations that will be obtained are similar to those in [29]. However, in that work the equations were derived from consideration of the forces acting on a cross-section normal to the axis of the root, which requires that both the cross-sectional cellular geometry and the mechanical properties of cells vary slowly along the organ axis. Here, we wish to derive the midline equations in a variational form, as this makes it simpler to deduce the discrete forms of the midline equations that we use for simulation. This is made more complicated by the fact that variational principles for dissipative systems (away from equilibrium) may not be as well grounded as those for elastic solids, for example. Below, we sketch out a plausible approach to the derivation of the governing system that we then implement computationally; given our focus on the practical implementation of the hybrid model, we here leave more rigorous considerations on one side.

#### 2.4.1. Energy dissipation for viscous models

In analogy with Stokes' law for viscous flow, we will use an energy dissipation principle to derive the evolution equations of the organ midline. Introducing a wall viscosity , the energy dissipation rate per unit area is , so the overall rate of energy dissipation is given by
2.16where this integral is over the region in two dimensions *R* occupied by the organ. We expect that the evolution of the midline will be an extremum of *D*, subject to the constraint that the work done by the internal pressure forces (turgor) be equal to the viscous dissipation, i.e. that
2.17

Thus, using the method of Lagrange multipliers, we find that the midline evolution must be a stationary point of 2.18

#### 2.4.2. Midline equations from variational principle

The evolution of the organ midline is a stationary point of (2.18) at each time *t*, with respect to variations in the RER and curvature generation rate . Using (2.13) gives
2.19

Here, for clarity we show the explicit dependence of all quantities on *s*, and *t*; in the remainder of this section we will omit the dependence on *t* for conciseness. We desire that the midline evolution be an extremum of *H*. Writing , , we have
2.20

Splitting this into two gives 2.21

To be a stationary point, these two integrals must independently be zero for all variations and . When we change the integration variables to midline coordinates , for which , this condition becomes
2.22*a*and
2.22*b*for all , where these integrals at each *s* are over the portion of the line normal to the midline that lies within the region occupied by the organ, (i.e. these amount to integrated versions of the Euler–Lagrange equations arising from (2.19)). Identifying
2.23and
2.24as the net tension and moment acting on each cross-section normal to the root midline, these equations are
2.25for all . Using the expressions for and , we can find and at each point *s* as the solution of the two-dimensional linear system
2.26and
2.27

When the organ curvature is small relative to its thickness, i.e. , we have
2.28if the midline is chosen such that , where here
2.29denotes the average of a quantity over the line normal to the midline at each *s*. These correspond to the two-dimensional version of the midline equations from [29].

### 2.5. Midline discretization

Simulation of the midline equations for complex geometries generally requires discretization of the midline geometry. The simplest approach is to treat the midline as a piecewise circular curve, as shown in figure 3; such a choice is natural as it leads to a piecewise constant curvature , although other discretizations with sufficient smoothness such as piecewise polynomial splines would be feasible. This consists of *N* sections with length and curvature , for . We then replace the midline coordinates for each vertex, by , where *j* is the index of the piecewise circular section containing the vertex and is the proportion of the distance along the circular section that the nearest point on the midline lies.

Substituting this discretization into (2.20) gives 2.30 2.31

where is the region , occupied by the organ within the *n*th midline section, where
2.32

On transforming the integration to midline coordinates, this can be written as 2.33where 2.34and 2.35

Thus, we find that the net tension and moment vanish over each section. This gives equations for the evolution of and , namely that for each , and are the solutions of the two-dimensional linear system
2.36*a*and
2.36*b*Here
2.37

For the purposes of simulation, it is easier to compute all integrals in terms of spatial coordinates, in particular this is because edges between vertices will be treated as straight lines in Cartesian coordinates, rather than in the midline coordinate system. Transforming back into Cartesian coordinates for the integration variables, and for simplicity subtracting times (2.36*b*) from (2.36*a*), the equations for each section become
2.38*a*and
2.38*b*

Here, 2.39

These integrals are approximated by numerical quadrature using the triangularization of the vertex-based tissue representation. Integrands are evaluated at the three vertices of each triangle; each vertex contributes to the integral corresponding to the midline section that contains it, where *A* is the area of the triangle and *f* the value of the integrand at the vertex.

### 2.6. Regulation of cell-wall mechanical properties

One key feature of this computational approach is that it can incorporate regulation of cell-wall mechanical properties through cell-scale processes. The regulation of growth in the primary root tip is likely to be the consequence of the combination of many different signals, some of which are typically diffusible hormones. Modelling the complex interactions between multiple hormones is beyond the scope of this study. We instead choose a simple model for the regulation of cell-wall properties. The hormone auxin is thought to be the key signal generating growth differences between the upper and lower sides of the root during gravitropic bending; evidence indicates that the epidermal (outermost) cell tissue plays a primary role in responding to this auxin asymmetry and generating bending [6]. All cells are taken to have uniform constant turgor pressure *P*, and the viscosities of the cell walls (apart from those in the epidermis) are chosen such that the RER of the root adopts the Peters-Baskin step-stool function ([25], and here described in appendix D)
2.40where is the distance from the root tip. However, as noted above, we wish to let the cell-wall properties of epidermal cells outside the meristem be regulated by auxin. As a simple method, we choose epidermal cell viscosities to be
2.41where *a* is the cellular auxin concentration (in arb. units, scaled such that the concentration in the central vasculature at the base of the organ is 2), is the half-saturation constant, controls the amplitude of the response and is chosen to be a smooth function which is approximately zero within the meristem and one outside it, namely
2.42where , , and erfc denotes the complementary error function. Cellular auxin concentrations are computed using a model described in appendix C.

### 2.7. Effect of midline evolution on cell-level vertex-based representation

The equations (2.38) govern the evolution in time of the discrete midline; here, we describe how this is used to update the vertex-based representation of the organ.

At the start of each time step, we calculate the discrete midline coordinates for each of the vertices in the cell-scale representation using the method described in appendix Ba. Note that, while usually constant between time steps, these coordinates may change owing to subdivision of the discrete midline, and new vertices may be introduced following cell division. Cell-level properties (such as auxin concentrations) are updated (e.g. using the auxin transport model described in appendix C) to permit the calculation of cell-wall viscosities.

Then, using the discretized midline equations (2.38), the evolution of the organ midline over one time step is calculated. The new positions of the vertices in Cartesian coordinates are calculated using the discrete midline coordinates and the relationship (B 1). As a result of this step, cell areas and cell-wall lengths may change, and these are recomputed. Concentrations of substances contained by cells (such as hormones) are updated to reflect the effect of dilution. Vertex velocites can also be calculated using the derivatives of (B 1) with respect to the midline variables; these are used for the specification of the axis of cell division.

From this, we see that the evolution of the macroscopic midline affects the cell-level properties of the simulation in a number of different ways. Changes in the areas of cells and the lengths of the walls between neighbouring cells modify the transport of auxin between cells, thereby affecting the auxin levels experienced by cells. These changes in cell area also change hormone concentrations through dilution, although for auxin this effect is small compared with transport. The evolution of the macroscopic midline also changes the distance of cells from the apex of the tissue, and, therefore, cell-wall mechanical properties.

### 2.8. Simulation overview

At each time step, the following operations are performed in order.

— Processes on the vertex-based mesh, such as auxin transport, are evolved for one time step.

— The viscosity for each (triangle of the) cell wall is calculated.

— The mean viscosity , the corresponding moments , , and the mean and moment of the pressure are calculated for each section of the midline.

— The lengths and curvatures of each of the midline sections are updated from (2.38) using a simple first-order forward Euler method.

— The positions of each of the vertices in Cartesian coordinates is updated using the coordinate transformation (B 1).

— Cell-level quantities, such as auxin concentration, are updated based on the change in cell size over the growth step.

— Cell division is performed.

— Midline sections, which have exceeded a specified threshold in length, are divided into two, and midline coordinates for the vertices in this section updated to reflect this change.

We note here that the numerical approximations made in this simulation are numerical integration the system of the ordinary differential equations for auxin transport [32], quadrature over the triangles to approximate the integrals occurring in (2.36), the piecewise constant approximation of the midline, and the forward Euler method used for the evolution of the lengths and curvatures of the midline sections.

## 3. Results

### 3.1. Organ bending

We first investigate the behaviour of this method when applied to the growth of an organ with a relatively simple, artificial geometry. The initial shape of the organ was taken to be a rectangle of height and length . This was subdivided in the *y*-direction into files of cells (of height ). Each file of cells was further subdivided in the *x*-direction into cells, starting at the left-hand end of the tissue (smallest *x*-coordinate), with each cell length being drawn from a uniform distribution on , where . Cell-wall viscosities were chosen to be , except for the uppermost layer of cells, which we took to be easier to extend, with . Cellular turgor pressures were taken to be . (Note that similar simulations have been carried out in the context of a vertex-element method; see [11, fig. 6].) The midline was initially discretized into sections of equal length.

Results of these simulations are shown in figure 4. From figure 4*d*, we see that the organ extends, and from figure 4*e* that the tip angle of the organ increases with time. The number of midline sections, *N*, increases adaptively with time, through midline sections being divided into two when they exceed . Simulations with increased number of initial midline sections () and smaller threshold for division () gave results which were in close numerical agreement; the relative difference in both midline length and tip angle after were less than .

We also explored the effect of the choice of the organ midline. The midline of an extended body is not well defined; one candidate for the organ midline, namely the medial axis [33], tends to branch towards the tip of the organ. Using a midline offset by a distance from the centre of the organ figure 4*b*, it can be seen that the length figure 4*d* and tip angle (figure 4*e*) of the resulting simulation are very close to those obtained using the centreline of the organ. Note that, in the simulations of figs 4–7, the small curvature correction terms (1 − *κη*) were omitted from the calculations of the averaged quantities in (2.38); including them makes the agreement even better.

The results were also compared with those generated by the vertex-element model previously described in [11]. This model uses triangular elements with anisotropic viscous mechanical properties to describe cell walls in the plane of the simulation, and viscoelastic elements to describe cell walls perpendicular to the plane of plane of the simulation. Choosing the parameters of the vertex-element model such that axial elongation was resisted only by the isotropic viscosity of cell walls in the plane of the simulation, we found in figure 4*f* that the length of the organ was in close agreement between the two approaches, but that the tip angle figure 4*g* increased more slowly in the vertex-element simulations than in midline simulations.

A second simulation was performed in which the initial organ geometry was curved. This initial geometry was generated using the rectangular geometry used before, with the midline through the centre of the object. The positions of the vertices were found in the discrete midline coordinates. A curvature of was then applied to the midline, and the new positions of the vertices in Cartesian components calculated. Unlike the earlier simulations, we took all cell walls to have the same viscosity . From figure 5*c–e*, it can be seen that while the organ extends, and the tip angle increases, the average curvature of the root midline remains almost constant over the course of the simulations. Such behaviour is consistent with (2.28). Again, the number of midline sections *N* increases as the organ extends ( initially, and at ).

### 3.2. Axial variation of cell-wall properties

To illustrate that the method is applicable to situations in which cell-wall properties vary with time, we consider a simulation in which the properties of cell walls were function of their distance from the tip of the organ. We take the dependence of cell-wall viscosity on distance to be a function of the distance from the tip , namely
3.1where is as shown in figure 6*a*. (Note that we use a different elongation rate function as the organ considered here is smaller than the full-root geometry used later.) We see the results of this simulation in figure 6*c*. The non-uniformity of the elongation rate profile leads to a non-uniform distribution of cell lengths along the organ axis, with the largest cells at the end of the simulation being found roughly from the tip (left-hand end of the organ). This is to be expected, as these cells have spent most of the simulation in the region of rapid expansion (i.e. in the elongation zone, which here lies between and from the tip). Although the model has no variation in cell-wall properties normal to the organ midline, slight bending of the midline can be seen in the simulations. The polygons occupied by each cell are divided into triangles, and cell-wall properties are assigned depending upon the positions of these triangles; this gives a discrete approximation to the smooth growth rate shown in figure 6*a*. Differences in cellular geometries lead to differences in the triangulation, and therefore cell-wall properties are not exactly uniform over each line normal to the root midline. In addition, the simulation currently uses a crude method to handle triangles which overlap the ends of each piecewise circular section, and differences in the degree of overlap across the root also contribute towards bending of the root midline. Such effects become more significant as the axial gradient in cell-wall properties increases, and indicate the importance of both the discretization of the cells and the midline being sufficiently fine in regions where cell-wall properties change rapidly along the root axis.

Cell division can readily be incorporated in this framework. At each time step (here with length ) cells divide with probability , provided that their area is greater than and that the distance from their centroid to the tip of the organ (measured along the centreline) is less than some threshold . The resulting distribution of cell areas is shown in figure 6*d*. Again, the largest cells can be found just beyond the region of rapid expansion. However, there is now a region of small cells near the tip owing to cell division events. This indicates the possibility of applying models of these types to explore the effects of changes in the spatial distribution of growth and division rates upon the resulting distribution of cell lengths in a multicellular tissue.

### 3.3. Whole root geometry with auxin transport

We now apply this method to a more complex scenario, in which changes in the cellular levels of the hormone auxin regulate the mechanical properties of individual cells.

As a cell-scale geometry, we use a template based upon a longitudinal cross-section of an *A. thaliana* root (from [2]). Owing to limitations of the field of view of the microscope, only part of the growing region of the root was acquired. The geometry was therefore extended by stretching the cells furthest from the tip, and then subdividing these cells using a length distribution that depended upon the cell file and the distance from the tip.

Cell-wall viscosity of epidermal cells (shootwards of the meristem at the root tip) was regulated by auxin in a manner described in §2.6. During gravitropism, it has been observed that the PIN efflux proteins (which transport auxin out of cells) change their localization in certain columella cells within the tip of the root [34]. When a root is growing vertically, PIN efflux proteins were found to be equally distributed on all sides of the columella cells (as illustrated in figure 7*a*). However, Kleine-Vehn *et al.* [34] observed that PIN proteins become localized to the lower sides of columella cells soon after gravistimulation (i.e. when a vertical root is rotated by a quarter-turn such that it is oriented in a horizontal direction).

In our simulation, we first find the steady state of the auxin transport equations on the static initial geometry with equal abundance of PIN proteins on all sides of the columella cells (figure 7*c*). These auxin concentrations were used as the initial conditions for the simulation. We then changed the localization of PINs on columella cells such that, after gravistimulation, they were only found on the lower faces of the columella cells (figure 7*b*). Soon after the start of the simulation, a difference in auxin concentrations between the upper and lower sides of the root tip can be observed (figure 7*d*). These differences are most pronounced in the cells bordering the columella cells (known as lateral root cap cells), and are smaller in magnitude in the epidermal cells, reflecting observations made using auxin sensors [35,36].

As the simulation proceeds, the increase in cell-wall viscosity for epidermal cells on the lower side of the root and the decrease on the upper-side of the root leads to the generation of curvature (figure 7*e*), and thereby, reorients the root tip. We terminate the simulations after the first 2 h; at later times during the gravitropic response, the asymmetry in auxin concentrations between the upper and lower sides of the root is thought to be switched off [35].

## 4. Discussion

Here, we have described a multiscale computational method that can be used for the simulation of multicellular plant organs undergoing axial elongation. This method captures the multicellular geometry of the plant tissue throughout the simulation, which permits simulation of cell-scale processes such as auxin transport. The macroscale equations for the evolution of the organ midline incorporate appropriately weighted averages of cell-level mechanical properties, which allows the simulation to capture the effects of changes in the mechanical properties of individual cells. This framework provides a convenient tool with which to investigate the intricate internal control systems that regulate growth patterns. However, the model is still vulnerable to the effect of geometric imperfections, whereby internal structural asymmetries can induce bending. This raises the possibility of internal regulatory mechanisms correcting such biases, analogous to the proprioception proposed in other models of gravitropic bending [37], which are naturally implemented in terms of the evolving centreline of the organ.

The method has promise for the investigation of a range of biological problems, in particular responses to gravity, and also responses to other stimuli such as water and light. The midline representation is well suited to long slender organs but it presents a restriction on the deformations that the organ can undergo. The model is unable to explain the complex pattern of growth and cell division that occurs near the root tip [38], or during the development of a new lateral root primordium [39]. Such deficiencies require either a different organ-scale description of growth in these regions [40,41], or coupling the method to the solution of a (finite-element) cell-scale simulation [12,13] in a local region. However, the theoretical details and practical implementation of such a coupled method have not yet been studied. The model is also unable to describe situations in which the organ can widen, for example the increases in root diameter that occur when a root is impeded by soil [42]. Nevertheless, the type of hybrid framework that we have described would seem to have promise for the description of a wide variety of multiscale and multicellular processes. It would also be of interest to extend this approach to study branching root structures, in which the root midline is represented by a (graph-theoretic) tree of circular sections rather than a one-dimensional list, in a similar manner to root architecture simulations (see the review of [43]). Additional consideration of the appropriate behaviour at branching points would be necessary in this case.

We note here that the model for gravitropism is deliberately simplistic, and here to illustrate how cell-scale models can be effectively integrated into this framework. There are still many details of the gravitropic response that are yet to be explained [44]; the sensing of auxin by cells is complex [45], and it is also thought that other mechanisms such as ion channels may play key roles [46].

## Competing interests

We declare we have no competing interests.

## Funding

We would like to thank the Biological and Biotechnology Science Research Council (BBSRC) for responsive mode (BB/J009717/1) and CISB awards to the Centre for Plant Integrative Biology. J.R.K. would like to thank the Isaac Newton Institute for Mathematical Sciences for its hospitality during the programme Coupling Geometric PDEs with Physics for Cell Morphology, Motility and Pattern Formation supported by EPSRC grant no. EP/K032208/1 and the Simons Foundation for a Fellowship during his participation in the programme.

## Appendix A. Derivation of rate of strain tensor in midline coordinates

From the definitions (2.4) of ** t** and

**, A 1and also**

*n*A 2

The first part of (A 2) can be written as A 3

Using (2.10) and (A 3) gives A 4

From the chain rule, temporal derivatives following material points and at fixed distances along the midline are related by A 5and similarly A 6

Here, the subscript indicates that the partial derivative with respect to *t* is to be taken with held constant (also known as a material or Lagrangian derivative), while the subscript *s* denotes a derivative with respect to *t* at constant distance *s* along the midline. Using these gives the compatibility conditions
A 7

Elsewhere in this paper, all temporal derivatives are following a material point, and we omit this notation.

The velocity in Cartesian coordinates of a point at fixed midline coordinates is given by A 8

Using (A 4), we have A 9which can be combined with (A 7) to give A 10where we define A 11

(As we will see below, is the strain rate component parallel to the organ midline.) Also, from (A 8) A 12

Differentiating (2.11), we have that A 13and from (2.4) and (A 1) this becomes A 14

Using these, we find that A 15and A 16so A 17

The rate of strain tensor is therefore given by A 18

## Appendix B. Discrete coordinate transformations

Here, we list for reference the explicit form of the transformations between the discrete piecewise circular midline coordinates and Cartesian coordinates.

**(a) Midline to Cartesian coordinates**

There is a relatively simple transformation between discrete midline coordinates and Cartesian coordinates ** x**, given by
B 1

*a*where , and denote the position, tangent and normal to the start of the

*j*th circular section (figure 3), given by B 1

*b*with B 1

*c*and B 1

*d*

Note that, for , (B 1*a*) becomes
B 1*e*

**(b) Cartesian to midline coordinates**

The inverse transformation from Cartesian coordinates ** x** to midline coordinates first involves identifying a circular section for which

**lies in the region between the normals to the circle at either end of the section (figure 8), i.e.**

*x**j*such that B 2If multiple suitable sections exist, the one whose midpoint is closest to the point

**x**is chosen. The values of and are then given by B 3

*a*and B 3

*b*if , B 3

*c*and B 3

*d*if , and B 3

*e*and B 3

*f*if . Note that the expressions here suffer from round-off errors for small , in which case they are instead calculated in the form B 3

*g*for and B 3

*h*for .

## Appendix C. Auxin transport

Auxin transport was simulated using the model of [2]; we summarize this briefly here. All parameter values for auxin transport are as listed there.

Auxin is considered to be uniformly distributed within each of the *N* cells, with concentration , where denotes the index of the cell (figure 9). Apoplastic compartments are associated with the common edges between each pair of adjacent polygons, and have auxin concentrations , where *i* and *j* are the indices of the two neighbouring cells (we identify ). The volume of each cell is taken to be proportional to its area , while the volume of each apoplastic compartment is taken to be proportional to its length times its thickness ( can be calculated from the polygonal network; is a property associated with each apoplastic compartment).

Auxin influx and efflux carriers facilitate the movement of auxin ions through the cell plasma membrane. Auxin fluxes from apoplast compartments (with label *ij*) into adjacent cells (with label *i*), , are given by
C 4, and are constants, , and LAX_{ij} are 1 if AUX/PIN/LAX proteins are present on the plasma membrane domain of cell *i* facing cell *j*, and 0 otherwise. For each cell *i* we have
C 5where is the auxin production rate in each cell, is the auxin degradation rate and the sum runs over all cells *j* that neighbour cell *i*.

Auxin transport between adjacent apoplastic compartments (i.e. compartments that meet a common vertex) occurs by diffusion. A (small) ‘vertex compartment’ is associated with each vertex, and the diffusive flux between a vertex compartment and the neighbouring wall compartment (per unit wall thickness) is taken to be
C 6where is the auxin concentration in the vertex compartment with index *q*. As the vertex compartment is small, the net flux into the compartment must be approximately zero. This provides equations specifying the vertex compartment concentrations
C 7

Substituting this into (B 6), and again considering conservation of auxin, equations for the wall compartment auxin concentrations are given by
C 8where the first sum (*k*) runs over the pair of cells (*i*,*j*) adjacent to the wall compartment, while the second sum (*q*) is over the two vertices at either end of the wall compartment. Auxin production and PIN efflux carrier locations were specified on the root template using the same rules as in [2], except that during our simulation of the response to gravity we modify the PIN efflux carrier locations on specific (columella) cells located within the tip of the root. Boundary conditions were applied at the base of the root; cells in the central files were set to have constant auxin concentration 1, while cells in the outermost three cell files (from the outside in, these are known as the epidermis, cortex and endodermis) were set to have constant zero auxin concentration.

## Appendix D. Peters–Baskin step-stool function

For reference, here we describe the function that we use to specify growth rates as a function of the distance from the organ tip. Using the definitions given in [25], we have
D 1*a*
D 1*b*
D 1*c*
D 1*d*

and
D 9*e*

for some constants (RER in meristem), (RER in elongation zone), (length of meristem), (distance of end of elongation zone from tip), (controls width of meristem-elongation zone transition), (controls width of elongation zone to mature transition).

## Footnotes

One contribution of 12 to a theme issue ‘Coupling geometric partial differential equations with physics for cell morphology, motility and pattern formation’.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.