A FINITE ELEMENT METHOD FOR GROWTH IN BIOLOGICAL DEVELOPMENT

We describe finite element simulations of limb growth based on Stokes flow models with a nonzero divergence representing growth due to nutrients in the early stages of limb bud development. We introduce a “tissue pressure” whose spatial derivatives yield the growth velocity in the limb and our explicit time advancing algorithm for such tissue flows is described in detail. The limb boundary is approached by spline functions to compute the curvature and the unit outward normal vector. At each time step, a mixedhybrid finite element problem is solved, where the condition that the velocity is strictly normal to the limb boundary is treated by a Lagrange multiplier technique. Numerical results are presented.

1. Introduction.The subject of limb development has generated much recent interest among biologists and physicists.The reasons for this interest are its importance as an example of well defined organogenesis during embryological development and that the biological and physical process underlying skeletogenesis are still far from clear.Among the many open questions are those related to how the overall limb shape develops.There exists a vast literature on the molecular biology involved in limb development [19], [21], [26], [27], [4], but how this molecular biology translates into patterning and growth is less clear.New experiments, however, suggest that the time is now ripe to investigate computationally how overall limb shape develops during vertebrate limb growth.As we are dealing with a complex free boundary problem, we will need to develop new algorithmic approaches to biological fluid flows in nonconvex domains, and this paper contributes to this area of biocomplexity.Such work is not only important from a conceptual viewpoint.The work's health implications are significant, because this research can be expected to influence several pharmacutical and bioengineering technologies.
In section II we describe the basic biology required to understand how the overall limb develops its complex asymmetric form.This is a rich source of biocomplexity, for concurrent with the development of overall growth and form, internal 340 C. M. MUREA AND H. G. E. HENTSCHEL spatiotemporal distributions of morphogens, activators, inhibitors and associated gene products occur that both depend on and control limb growth and form.
In section III we describe the types of free boundary problems associated with creeping flows in nonconvex growing domains that occur during organogenesis.Specifically, in such developing domains, we need to solve for the growth velocity in the limb using a Stokes flow with a nonzero divergence representing local nutrients generating the observed growth.We introduce a "tissue pressure" whose gradient yields the growth velocity and calculate the resulting scalar field using biologically plausible boundary conditions including expressions for the tissue pressure at the limb boundary formulated in terms of the instantaneous limb curvature, and the imposition on the internal epithelial surface of the limb the biologically plausible boundary condition that the tangential velocity field is zero.At the boundary joining the limb with the main trunk of the vertebrate embryo, we impose a less restrictive slip condition for the growth velocity.The growth rate of the limb is then given by the normal velocity of the fluid at this moving boundary.
In section IV we describe in detail a new finite element algorithm for studying such flows.Mathematically, several general frameworks for solving Stokes equations in moving domain have been developed.These include the arbitrary Lagrangian Eulerian together with the finite element method [16], the level set method [25], the immersed boundary method [22], and the particle method [5], [17].The approach we develop here is an explicit time advancing scheme belongs to the framework called "front-tracking methods." In section V we apply our algorithm to track the free boundary and internal growth velocity field in both initially semicircular (in the very early vertebrate embryo the limb bud is approximately semicircular, see Figure 1) and in nonconvex domains.Finally we discuss our results in section VI.
2. Biology underlying vertebrate limb development.Studies of limb development involve many interconnected questions, from what are the mechanisms controlling overall limb shape to how internal structure in the growing limb bud develops.There exists a large literature on the molecular biology involved in limb development [19], [21], [26], [27], [4], but how this molecular biology translates into growth and form is less clear.An examination of limb physiology shows that this is a complex process.Clearly defined axes exist-proximal-distal, anterior-posterior, and dorsal-ventral.Different sizes and shapes for the stylopod (one bone in the upper arm or thigh), zeugopod (two bones in the forearm or calf) and autopod (different numbers of nonidentical segmented digits) are observed in the tetrapod limb.From a computational viewpoint, the situation is equally challenging.We need to understand how the overall limb develops its complex asymmetric form and, concurrent with this process, how internal asymmetric spatiotemporal distributions of morphogens, activators, inhibitors and associated gene products result in the skeletal limb forms created by evolution.
Recently there have been new insights into skeletal development.Much recent evidence suggests that the early stages of skeletal pattern formation in the developing vertebrate limb depend on complex dynamics involving several growth factors and differentiation of cells with receptors that allow response to these factors.We have shown that this biology is indeed sufficient to generate the basic patterning of the generic vertebrate limb [14].Computational work in three dimensions [18] [2] has both confirmed and extended this mechanism [14].It incorporates a core set of cellular-biochemical processes known to occur in limb bud mesenchyme and is capable of generating wrist and ankle spot-like elements [1] in addition to the longer stripe-like bone elements.
But nearly all studies described above were carried out in growing rectangular or parallelpiped domains.Real biological development, however, involves both growth and changes of form of free-moving boundaries [29].This is certainly the case of the limb bud (see Figure 1).Therefore in this paper we describe tissue flows and their associated algorithmic implementation that both help shape the embryonic limb and help convect internal morphogens and gene products vital to the development of internal form.Because we are concentrating on external epithelial domain grown and form, we suppose it can be described mathematically as a free-moving boundary problem controlled by internal Stokes flow attributable to internal tissue growth fed by a continuous source of nutrients S(x).This boundary value problem is similar in some respects to other two-fluid flow interfaces in Hele-Shaw cells with surface tension.Such flows are known to give rise to nontrivial interfacial structure because of the existence of the Mullins-Sekerka instability.More generally they fall into the general category of creeping flows in the presence of moving boundaries.The fact that creeping flows are involved can easily be seen by estimating the associated growth Reynolds number.When the typical length scales in the developing limb are L ∼ 10 −1 cm, while typical convective velocites induced by growth are V ∼ 10 −6 cm/sec, while the kinematic viscosity of water ν ∼ 10 −2 cm 2 /sec, the resulting flows typically have Reynolds numbers in the range Re = LV /ν ∼ 10 −5 .Because mesodermal cellular flows will have effective viscosities significantly larger than water, the Reynolds numbers will be even smaller [8].Another similarity is that, at least in the first approximation, growth is two-dimensional.The developing limb fibres connect the dorsal and ventral walls of the limb bud [3] leading to two-dimensional flows.Also, the phenomenon of convergent extension [28], in which flattened cells tend to develop in the two-dimensional plane defined by the proximal-distal (shoulder to digit tip) and anterior-posterior (thumb to little finger) axes, will help justify a growth description in terms of two-dimensional flows.
There are, however, several differences from the usually studied incompressible Stokes flows.As mentioned above, growth due to mitosis and nutrients ensures that material is constantly being added (and sometimes removed when cell death or apoptosis occurs).In addition, the surface tension in the developing limb embryo is heterogeneous, because the epithelial cell layer is weaker near the AER, and consequently boundary conditions will result in a more complex boundary value problem than those studied in Hele Shaw cells.
There exists a previous integration of the influence of growth on form in the context of the avian limb bud [6].There it was assumed that the flow was a Navier-Stokes flow in the presence of homogeneous boundary conditions.In addition, the growth was assumed to be strongly dependent on the local FGF concentration released by the apical ectodermal ridge (AER) near the tip of the limb bud.More recent evidence suggests that mitosis and therefore growth are not strongly influenced by the local FGF concentration, but rather the main influence of mitosis is on cell differentiation.
The need to develop general methods for integrating creeping flows in biology has motivated our search to develop finite element algorithms for creeping flows in the context of avian limb development.We also believe that similar finite element algorithms will be useful in the more general context of organogenesis.

3.1.
The basic ingredients.We consider, therefore the following minimal model, which incorporates the key features of this biological growth.Addition of material occurs everywhere either at a constant rate S or more generally the rate of growth can be assumed to be S(x, t) as addition of material could be both spatially varying and have a temporal dependence because of genetic switching mechanisms.This means that that the tissue flow in the limb will include a continuous distribution of sources and therefore obey where v is the fluid velocity.We treat growth of the limb as due to a creeping flow, because of the very low Reynolds numbers involved [8].Therefore we can expect the flow to obey the Stokes equation with volume source where p is a pseudopressure field defined by p = P − p air , where P is the pressure of the fluid and µ is the viscosity of the fluid.Finally we need boundary conditions.As it appears there is no flow of material into the main body of the organism, we shall take slip boundary conditions at the boundary of the limb connected to the main body where ν is the outer unit normal vector to the boundary, while the elastic properties of the epithelial layer of cells forming the skin layer will result a pressure at the free growing boundary that obeying where γ(s) is the effective surface tension of the limb at a point s on the free boundary [9] while κ(s) is the limb curvature at point s.
The equation of the normal velocity of the free boundary is 3.2.Detailed structure of the dynamics.We now describe in detail the structure of the creeping flow dynamics we wish to integrate.We study the evolution of a bounded connected open domain Ω(t) of R 2 with boundary where Γ 1 (t) and Γ 2 (t) are two nonempty subsets of ∂Ω(t).Here t ≥ 0 is the time.We assume that Γ 1 (t) is a nonclosed curve of class C 2 and its ends evolve on the Oy axis (see Figure 2).The boundary Γ 2 (t) is the segment which has the same ends as Γ 1 (t).Let ν = (ν 1 , ν 2 ) denote the unit outward normal vector and τ = (−ν 2 , ν 1 ) the unit tangential vector to the boundary.In the moving domain Ω(t), we have to find the velocity v( and the pressure p(•, t) : Ω(t) → R of the fluid, such that where µ > 0 is the viscosity of the fluid, f = (f 1 , f 2 ) are the applied volume forces, S > 0 is the rate of growth, γ > 0 is the effective surface tension and κ is the curvature of Γ 1 (t).We use the sign convention that convex domains have positive curvature of the boundary.
Remark 1.In a previous work [20], we studied a similar problem, where the fluid velocity was not necessarily normal to the boundary Γ 1 (t).In [15], Guermond and Quartapelle suggest the importance to prescribe either the normal or the tangential component of the fluid velocity, to establish the stability of the finite element approximations.In the present paper, the fluid velocity is supposed to be normal to the boundary Γ 1 (t) (see equation ( 8)).Even if the assumption concerning the fluid velocity on free boundary is not the most appropriate from the biological point of view [7], this constraint is required by the boundary condition concerning the pressure (9).
Equation ( 11) is a natural boundary condition associated with the essential boundary condition (10).
The boundary Γ 1 (t) evolves according to the law where V ν is the normal velocity of Γ 1 (t).
The combination of ( 8) and (10) requires that the the boundaries Γ 1 (t) and Γ 2 (t) are normal at the intersections, more precisely where ν |Γ 1 (t) and ν |Γ 2 (t) are the outer unit normal vectors to Γ 1 (t) and Γ 2 (t), respectively.Without this condition, the fluid velocity will not be continuous in the two corners.This requirement is biologically plausible in terms of flow, though of course the boundaries of real limb domains are not exactly normal to each other in geometrically.
We know the initial domain We consider the problem ( 6)-( 14) of determining the evolution of Ω(t) and to find the velocity v(x, y, t) and the pressure p(x, y, t) for t ∈ [0, T ], where T > 0 is a given real constant.
4. The finite element algorithm.We develop a key algorithm to study such creeping flows in the presence of moving boundary conditions.
To evaluate the curvature terms, the boundary is approached by cubic spline interpolation, which gives a curve twice continuously differentiable.The curvature is computed using the parametrization of the splines.
At each time step, a new mesh is generated, but the generation of the current mesh is independent from the previous one.
To numerically solve these equations we use finite element methods.We introduce for each t ∈ [0, T ] the following Hilbert spaces: The weak form of the problem ( 6)-( 11) is to find where Remark 2. Some authors use the notations ∂x − ∂v1 ∂y for two-dimensional vectors.We have to note that, in this case, v × ν and ∇ × v are scalars, not two-dimensional vectors.The bilinear applications (21) and ( 23) could be rewritten more concisely as For the weak form of Stokes equations with Dirichlet boundary condition on the velocity, the standard bilinear form is The boundary condition on the pressure requires the use of the alternative bilinear form (21).
From Green's formula and the following identity we can prove that if v, p is a strong solution of ( 6)-(11), then v, p, ω = ∇ × v is a solution of ( 18)-(20).
The system ( 18)-( 20) is a mixed-hybrid-like problem in that the trial spaces W(t), Q(t) and Λ(t) are independent and some trial functions are defined on the physical domain, while others are defined only on the boundary.The main advantage of this framework is the treatment of the constraints ( 7) and ( 8) by the Lagrange multiplier technique, consequently, we are not forced to use finite elements which verify (7) and (8).
The finite element approximation of Stokes equation with boundary condition on pressure was studied in [23], [24] and [12], but the condition (8) was treated in a strong way.
We denote by ∆t the time step and by N = T /∆t the number of time steps.We approximate Γ 1 (n∆t) by a polygonal line , where the vertices A n i have the coordinates (x n i , y n i ) for i = 0, . . ., M .It is assumed that x n 0 = x n M = 0 for all n, which implies that A n 0 and A n M evolve on the Oy axis.We denote by Ω n h the polygonal domain bounded by Γ n 1,h and the segment [A n M , A n 0 ].Let T n h a triangular mesh of Ω n h .For the approximation of the fluid velocity v, we have used the triangular finite elements P 1 + bubble, also called MINI elements, introduced by Arnold, Brezzi and Fortin (see the general reference [10]), for the fluid pressure p the triangular finite elements P 1 and for the Lagrange multiplier ω, the segment finite elements P 1 .
We denote by v n h , p n h , ω n h the finite element approximations of v(n∆t), p(n∆t), ω(n∆t), respectively.
Let ν n h denote the unit outward normal vector to the boundary Γ n 1,h , which is constant on each edge A n i , A n i+1 .To compute in (23) the integral term containing ν, we have used the approximation In a similar way, we can approximate in (24) the term Γ1(t) S (w • ν) ds.
We have for i = 0, . . ., M − 1 where We set m x 0 = m x M = 0 and m x 1 , . . ., m x M −1 is the solution of the linear system For y(ξ), we have similar formulas.For i = 0, . . ., M − 1 where The linear system to solve is The curve S n has a continuous curvature given by In the sequel, κ n i = κ(ξ i ) stands for the curvature in the vertex A n i using the spline functions.
Do compute in (24) the integral term containing the curvature, we have used the approximation 2. An explicit time-advancing scheme.From ( 12), a point on the boundary Γ 1 (t) moves along the normal to the boundary with the velocity v • ν.
In a previous section, we introduced ν n h the unit outward normal vector to the boundary Γ n 1,h , constant on each edge A n i , A n i+1 , which, however, is not well defined on the vertices A n i .
To compute the position of the vertices A n+1 i of the polygonal line Γ n+1 1,h , we will use the normal to the spline function S n given by More precisely, we set which is the forward Euler's scheme for the numerical approximation of (12).
For each n from 0 to N − 1 Step 1: Generate T n h a triangular mesh of Ω n h .Knowing the boundary points A n i , i = 0, . . ., M , the mesh can be generated automatically, using FreeFem++ [13].
Step 2: Compute the spline function . ., M .The details were presented in section 4.1 "Treatment of the curvature terms." Step 3: Compute κ n i = κ(ξ i ) the curvature at each vertex A n i , using the formula (26).
Step 4: Find v n h , p n h , ω n h the finite element solution of ( 18)- (20).After the finite element discretization, the problem to solve is a symmetric linear system, not positive-definite, of the form Step 5: Compute ν n S the unit outward normal vector to the spline function S n using the formula Step 6: Move the vertices of the boundary using the forward schema (27).
5. Numerical tests.We have tested the algorithm presented in this paper for two kind of geometrical shapes: a semicircle and a nonconvex domain.Both examples were also discussed in [20] using Darcy's law for the fluid flow, while the actual model is based on the Stokes equations with volume source.
5.1.The initial domain is a semicircle.First let us consider the academic case in which the initial domain is a semicircle of ray R 0 .We assume that the rate of growth S and the surface tension γ are constants, and we set the applied volume forces f = (0, 0).Of course, these assumptions are not biologically reasonable, but in this case, we know an exact solution of the free boundary problem ( 6)-( 14).
If we set the parametric representation of Γ 1 (0) as the evolution of the boundary Γ 1 (t) is described by The velocity and the pressure have the form .
The algorithm was implemented using the language FreeFem++ [13], and the numerical results were displayed using gnuplot [11].
The number of the segments of the polygonal line Γ n 1,h is M = 32, and the number of edges on the vertical boundary [A n M , A n 0 ] is 20 for all n.The initial mesh has 200 vertices, 346 triangles, and after 100 time steps we use a mesh of 207 vertices, 360 triangles (see Figure 4).
The computed velocity (see Figure 5) is radial as the theoretical solution.

5.2.
A nonconvex initial domain.Let now consider a case when the initial domain is nonconvex.The boundary Γ 1 (0) has two flat parts (one on the bottom and on the top), and three semicircles of rays r1 = 0.6, r2 = 0.2 and r3 = 0.2, respectively.
The simulation was performed for: µ = 1, f = (0, 0), S = 2, γ = 0.5, ∆t = 0.0001 and N = 1800.The number of segments of the polygonal line Γ n 1,h is M = 48 and the number of edges on the vertical boundary [A n M , A n 0 ] is 20 for all n.The initial mesh has 302 vertices and 534 triangles (see Figure 6).In [20], the flat parts of the boundary Γ 1 (0) do not move; consequently, the growth is only in the Ox direction.As we see in Figures 8 and 7, the domain growths also appear along the Oy axis.
In view of Figure 8, we can suppose that the domain will evolve to a domain with a cut, while in [20] the same domain seems to evolve to a convex domain.
Using the same time step ∆t = 0.0001, we have performed more than 1800 iterations, while in [20] we have obtained a self-intersecting moving boundary after 182 iterations.6. Discussion.The finite element algorithm described above contains the basic ingredients required for the study of the development of biological form as a consequence of growth via creeping flows due to nutrient addition to a closed bounded domain surrounded by epithelial cell walls.Perhaps the most interesting results of these simulations can be seen in Figure 7.The internal growth velocity field has clearly developed a nontrivial spatial and temporal structure that will affect not only the external shape of the limb bud but also internal processes, such as the spatiotemporal distribution of gene products, and consequently the development of internal structure.Our approach of course represents a minimal model, in that it incorporates many crucial biological processes but at the same time leaves out many critical elements that need to be considered in future algorithms.Such critical elements include internal domains created by skeletal elements that can both channel and hinder fluid flow.Also there exist a core set of cellular-biochemical processes known to occur in limb bud mesenchyme that will further sculpt the developing form of the organism.For example, nonuniformly distributed gene products, such as Sonic Finally, of course, explicit three-dimensional simulations need to be undertaken to investigate how spatially varying surface tension caused by heterogeneous epithelial properties affects growth and form.This is more complex algorithmically, since in the two-dimensional case powerful tools exist that can automatically generate moving meshes when the parametric description of the boundary is provided.We are now studying the three-dimensional version of these algorithms and are especially interested in how heterogeneous properties of the epithelial layer will affect growth and form.
We believe, however, that all such considerations can be incorporated into developments of the basic algorithm proposed here and need to be studied further.

Figure 1 .
Figure 1.A schematic drawing of the early stages of both external and internal growth of the embryonic verebrate limb bud.

Figure 2 .
Figure 2. Schematic illustration of the free boundary problem

Figure 3 .
Figure 3. Spline approximation of the initial boundary (left), after 50 time steps (middle) and after 100 time steps (right)

Figure 4 .
Figure 4.The initial mesh (left) and the mesh after 100 time steps (right).

Figure 5 .
Figure 5.The initial velocity (left) and after 100 time steps (right).The scaling factor is 0.2.

A-Figure 6 .
Figure 6.The initial mesh (left) and the mesh after 1800 time steps (right).

Figure 7 .
Figure 7.The velocity at different time steps.The scaling factor is 0.2.

Figure 8 .
Figure 8. Spline approximation of the boundary at different time steps