Volumetric finite-element modelling of biological growth

Differential growth is the driver of tissue morphogenesis in plants, and also plays a fundamental role in animal development. Although the contributions of growth to shape change have been captured through modelling tissue sheets or isotropic volumes, a framework for modelling both isotropic and anisotropic volumetric growth in three dimensions over large changes in size and shape has been lacking. Here, we describe an approach based on finite-element modelling of continuous volumetric structures, and apply it to a range of forms and growth patterns, providing mathematical validation for examples that admit analytic solution. We show that a major difference between sheet and bulk tissues is that the growth of bulk tissue is more constrained, reducing the possibility of tissue conflict resolution through deformations such as buckling. Tissue sheets or cylinders may be generated from bulk shapes through anisotropic specified growth, oriented by a polarity field. A second polarity field, orthogonal to the first, allows sheets with varying lengths and widths to be generated, as illustrated by the wide range of leaf shapes observed in nature. The framework we describe thus provides a key tool for developing hypotheses for plant morphogenesis and is also applicable to other tissues that deform through differential growth or contraction.

The subscripts range over 1 − 3, the spatial dimensions. D is the elasticity tensor, a fourth rank tensor with 81 elements, with certain symmetry properties implying that it has only 15 independent parameters. While GFtbox supports the use of arbitrary elasticity tensors, which may also vary over the material, the examples in this paper and in most of our use of GFtbox assume an isotropic and uniform elasticity tensor, for which the elasticity tensor simplifies to the form given in equation 2 below, with only two independent parameters. The equilibrium conformation is that which minimises the total deformation energy 1 2 ij σ ij ij . The transient diffusion-decay equation for a morphogen p is dp dt = K∇ 2 p − Ap, where K is its diffusion constant and A its decay rate. As with elasticity, GFtbox supports K and A that vary in time and across the tissue, but in all of the examples these are constant and uniform.
Boundary conditions for the diffusion models take the form of fixing p to certain values at certain places.
GFtbox also supports setting rates of production of morphogens and interactions between morphogens.

Definition of the finite elements
The finite elements used in all the volumetric models presented here are first-order tetrahedra, with quadratic Gaussian integration. GFtbox also supports first-order pentahedra and hexahedra, also with quadratic Gaussian integration, but does not support remeshing for these.

Visualisation
For a solid object only the outside can be seen. Additional tools must be implemented to see inside, using cutaway sections or partial transparency. A field of polarisation arrows or tensor indicatrices on a surface is legible to the eye, but a cloud of such decorations in a spatial volume may be more difficult to interpret. Tools were developed for creating legible 2D representations of these 3D objects, and for exporting to 3D formats either for postprocessing in 3D graphics tools or for presenting directly in web pages or augmented/virtual reality applications.
Grid lines are embedded in the tissue and move with it, to illustrate the deformation arising through growth. (The grid is independent of the finite element decomposition, which is generally much finer.) Poisson's ratio is zero unless otherwise stated. (The behaviour of all the examples is qualitatively not much affected by the value.) Within each figure, images are drawn to the same scale, except when scale bars indicate otherwise.
Distance and time are measured in arbitrary units. Strain is a dimensionless quantity. Growth rate and strain rate have dimension 1/time. In the examples, the tissues generally have a radius or semi-diameter of 1 unit, the maximum growth rate is 1 per time unit, and the examples are run for a time suitable for making the resulting form clear. The visualisations of residual strain rate should be considered to be of the ratio of residual strain rate to maximum growth rate. Since the latter is always 1 this is numerically the same, but is a dimensionless quantity with a physical meaning.

Description of models
Here we give a systematic summary, figure by figure, of the GFtbox models illustrated.
When growth is isotropic (k par = k par2 = k per ) we write k for the common value of the growth rates.
The organisers referred to in the text as +ORG, −ORG, +ORG2, and −ORG2 are in the GFtbox models named ID SOURCE, ID SINK, ID SOURCE2, and ID SINK2. Each of these takes a value of either 1 or 0 at each vertex of the mesh. When one of these defined to be 1 "near" a certain point of the surface and 0 elsewhere, the neighbourhood in which it is 1 is defined as all points subtending an angle of no more than 15 • with the given point, and within a distance of 0.3 of the surface.
Except for Figure 16, the coordinate frame in all models, as depicted, has the x axis from front left to rear right, the y axis from front right to rear left, and the z axis upwards. Figure 16 has the x axis from rear left to front right and the y axis from front left to rear right.
The centre of the initial state of every model is at the origin of coordinates.

Figure 6
The initial tissue is a sphere of unit radius.
The growth pattern is set up in the initial state and thereafter moves with the tissue.
Growth is isotropic everywhere.

Figure 7
The initial tissue is a thin disc of unit radius, with an initial small perturbation to its flatness.
The growth pattern is set up in the initial state and thereafter moves with the tissue.
Growth is isotropic everywhere.
The initial tissue is a cube of diameter 2.
The growth pattern is set up in the initial state and thereafter moves with the tissue.
Growth is isotropic everywhere.
Growth rate is a function of the distance r from the centre: for r > 1 (and is therefore 1 at the corners). (a,b,e,f) The initial tissue is a cube of diameter 2.
The growth pattern is set up in the initial state and thereafter moves with the tissue.
Growth is isotropic everywhere. k = (x + 1)/2, hence 0 on the left-hand face and 1 on the right-hand face.
(c,d,g,h) The initial tissue is a thin square parallel to the xz plane, of diameter 2. The growth pattern is defined in the same way as for the cube.

Figure 10
The initial tissue is a sphere of unit radius.
The growth pattern is set up in the initial state and thereafter moves with the tissue.
Growth is isotropic everywhere.
POL is fixed at 1 where ID SOURCE = 1, and has diffusion constant 0.1 and decay rate 0.5.

Figure 12
The initial tissue is a cube of diameter 2, whose centre is at the coordinate origin. k par = 1, k par2 = k per = 0.

Figure 13
The initial tissue is a sphere of unit radius. These models differ only in how the polarity is organised. For all images in the third column, k par = 1, k par2 = 0, k per = 0. For all images in the fourth column, k par = 0, k par2 = 1, k per = 1. POL is fixed at 1 where ID SOURCE = 1, and has diffusion constant 0.1 and decay rate 0.5.

Figure 15
The initial tissue is a cube of diameter 2, whose centre is at the coordinate origin. It diffuses by the steady state diffusion equation.
POL2 is similarly determined by ID SOURCE and ID SINK.

Figure 16
The initial tissue is a sphere of unit radius.
ID SOURCE = 1 at all vertices on the equatorial plane and 0 elsewhere.
POL is fixed at 1 where ID SOURCE = 1, and has diffusion constant 0.1 and decay rate 0.5.

Complications with elisions
If a tetrahedron has any of the thin forms of Figure 1 (a-c), it would be desirable to merge the sets of vertices that lie very close to each other, eliminating this tetrahedron, some or all of its faces, and every other tetrahedron that shares any of those faces, as in Figure 1(d-f). However, in some circumstances, merging adjacent vertices in this manner can create an invalid mesh. This is because to be compatible with the finite element method, all internal edges of the tetrahedral mesh should belong to at least three tetrahedra, and this condition may be violated. For example,consider the three tetrahedra shown in the upper left of Figure 2 (a) and in exploded view (b). The internal edge shown in blue is shared by all three tetrahedra. If we elide the edge marked in red in Figure 2 (a), we obtain tetrahedra which have all of their vertices in common (c) and the blue edge is now shared by only two tetrahedra, which is illegitimate. Figure 2(d) shows the exploded view. Either the two remaining tetrahedra are both flat, or one has been everted. Neither case is compatible with the finite element method. Thus before carrying out elisions, there must be a check as to whether they would generate a legitimate mesh. Alternatively, the two coincident tetrahedra could be eliminated by a transformation similar to that of Figure 3(iii), but again, the validity of the transformation would have to be checked. Due to these complications, we have not yet implemented elision transformations.

Mathematical derivations
Radially symmetric growth in a sphere The symmetry of the spherical models of Figure 6 allows the equations of elasticity to be reduced to an ordinary differential equation that can be solved analytically. Comparing the result with that of the finite element calculation provides a validation test. We outline the calculation here, generalising the setting to a possibly hollow sphere in any number of spatial dimensions (i.e. including the case of a disc or annulus confined to the plane).
Suppose we have a possibly hollow sphere whose internal and external radii are r 0 ≥ 0 and r 1 > r 0 . Let it have a growth rate at each point of g r in the radial direction and g t in all tangential directions, each of these being a function only of r. We want to calculate the equilibrium configuration, assuming that it remains spherically symmetrical.
For generality, we solve the problem in an arbitrary number of dimensions k + 1 (k being the number of independent tangential directions). The practical cases are k = 1 and k = 2. The first task is to determine the stiffness tensor of an isotropic elastic substance. The simplest expression of this tensor, valid in any number of dimensions, is in terms of Lamé's first and second parameters λ and µ (see e.g. [3], eqn. 5.4.2): C pqrs = λδ pq δ rs + µ(δ pr δ qs + δ ps δ qr ) (p, q, r, s : 1 . . . k + 1) C has three families of non-zero elements: C pppp = λ + 2µ, C ppqq = λ (p = q), and C pqpq = µ (p = q).
By considering the effect of either an isotropic compression or a uniaxial compression we can express the Lamé parameters in terms of the bulk modulus K and Poisson's ratio ν: µ is the shear modulus, which must be positive for a physically realistic material. Therefore ν < 1/k. For convenience we define = 1 − kν.
For a strain whose principal components are v = (v 1 , ..., v k+1 ), the energy density, up to a constant factor is ν(Σ i v i ) 2 + Σ i v 2 i . (Notice that the bulk modulus K is part of that constant factor, as there are no external forces, only the material interacting with itself.) Now let a spherically symmetrical deformation of the sphere be described by a function ζ(r) which maps the radius of a point in the original sphere to its radial displacement in the deformed sphere. Its new radius is thus r + ζ(r). The radial strain is ρ = ζ − g r , and the tangential strain τ = ζ/r − g t (on each of k axes). Substituting these in the expression for energy density, multiplying by r k (which is proportional to the area of the surface at radius r), and ignoring constant factors, we obtain the energy density per unit radius, which we write as: The equilibrium conformation is obtained by minimising the integral of this from r 0 to r 1 . By the calculus of variations, the integral is minimised if L satisfies the Euler-Lagrange equation: This works out to the basic equation: where α = k +ν . At a free surface the boundary condition is: This is a quantity that we might call the dual Poisson's ratio, hence the notation, being the expansion in one direction resulting from a uniform compression in all the perpendicular directions. (For k = 1 this is identical to Poisson's ratio. For larger k it is larger.) For a solid sphere (r 0 = 0) the boundary condition at the centre is dictated by continuity: ζ(0) = 0.
If g t and g r have the forms a t r n and a r r n for constants a t and a r , the equation can be solved analytically. The solution has the form: for constants B, C, and D, where P n (r) = r 1 r n−1 dr. This is (r n − 1)/n for n = 0 and log r for n = 0. Defining the function P n in this way avoids having to treat n = 0 as a special case. B is determined by substituting the general solution into Eq 4 (in which C and D cancel out), and C and D are obtained from the two boundary conditions. Calculating these constants is complicated by several special cases that must be treated separately. In particular, if = 0 the equation has no unique solution, as any volume-preserving flow can be added without violating the equations. In that case we can take the limit of the solution as tends to zero, which exists and has the same general form. The different boundary condition for r = 0 also implies a different calculation of the constants (C must be zero in that case). There is a further special case where n = −k − 1, for which the solution has a different form: the constants then being calculated in the same way.

Some special cases of the exact solution Zero residual strain
For the radial and tangential growth rates to be such that there is no residual strain, we must have ρ = τ = 0. This implies g r = g t + rg t . If g t = ar n , this implies g r = a(n + 1)r n = (n + 1)g t . From τ = 0 we obtain ζ = ar n+1 (which is the integral of g r from 0 to r).

Uniform isotropic growth
Uniform isotropic growth is obtained when n = 0 and a r = a t . This implies that C = 0 and D = a, as we expect.

Isotropic non-uniform growth
The examples of Figure 6 have isotropic growth rate depending on radius of (a) 1, (b) r 2 , and (c) 1 − r 2 . The initial sphere has radius 1 and Poisson's ratio zero. The exact solutions that are plotted in the figure are r, 0.4r + 0.2r 3 , and 0.6r − 0.2r 3 .