## Abstract

Mechanical forces are closely involved in the construction of an embryo. Experiments have suggested that mechanical feedback plays a role in regulating these forces, but the nature of this feedback is poorly understood. Here, we propose a general principle for the mechanics of morphogenesis, as governed by a pair of evolution equations based on feedback from tissue stress. In one equation, the rate of growth (or contraction) depends on the difference between the current tissue stress and a target (homeostatic) stress. In the other equation, the target stress changes at a rate that depends on the same stress difference. The parameters in these morphomechanical laws are assumed to depend on stress rate. Computational models are used to illustrate how these equations can capture a relatively wide range of behaviours observed in developing embryos, as well as show the limitations of this theory. Specific applications include growth of pressure vessels (e.g. the heart, arteries and brain), wound healing and sea urchin gastrulation. Understanding the fundamental principles of tissue construction can help engineers design new strategies for creating replacement tissues and organs *in vitro*.

## 1. Introduction

Physical scientists strive to express the fundamental laws of nature in mathematical form. Doing this enables them to test these laws using precise quantitative measurements, as well as to predict behaviour outside the boundaries of the laboratory. Without mathematics, it would be difficult to confirm, for example, the validity of Einstein’s general theory of relativity.

In contrast, the laws of biology rarely are written in the language of mathematics. This is due, in part, to a general feeling (among scientists from all fields) that biology is primarily a qualitative science. More importantly, living systems possess characteristics that set them apart from non-living systems. For example, biological systems respond actively to external stimuli, and their behaviour is influenced by the need to survive natural selection (Forgacs & Newman 2005; Nanjundiah 2005). Evolution also imparts a certain stochastic quality to biology via random genetic mutations that may or may not be passed on to future generations.

Randomness, however, should not, in itself, rule out the existence of quantitative laws for biology. In physics, for example, quantum mechanics provides only a probability that a particle travels with a certain velocity at a certain place at any given instant. Yet the accuracy of this theory has been confirmed by numerous experiments. Another example comes from materials science, where theories of failure are subject to the random appearance of cracks in a material.

In fact, some aspects of the theory of evolution may be analogous to behaviour encountered in the stability analysis of nonlinear inanimate systems. Near a bifurcation point, for instance, small random perturbations can send a system down one of two, possibly widely disparate, paths. Correspondingly, evolutionary changes may be viewed as a series of bifurcations biased by genetic mutations and developmental constraints (Alberch 1980; Oster & Alberch 1982). Moreover, the governing equations or parameters in a biological theory may evolve in time, with the form of the equations determined by natural selection.

Within this context, the present study explores a possible unified theory for the mechanics of morphogenesis, which is the creation of biological form. The focus here is primarily on epithelia (cell sheets), which play important roles in embryonic development. In an earlier study (Taber 2008), we used mathematical modelling to examine a theory for epithelial morphogenesis based on a mathematical formulation of Beloussov’s hyper-restoration (HR) hypothesis (Beloussov *et al.* 1994; Beloussov 1998). According to this hypothesis, epithelia *actively* elongate in response to increased stress and shorten in response to decreased stress. This response tends to restore, but overshoot, the original stress to the other side. For a number of illustrative problems, our models showed that the stress overshoot can drive and sustain morphogenesis.

The present study is motivated by thought experiments that reveal contradictions in some predictions of the HR hypothesis when applied to general developmental processes. One example is growth of a pressure vessel (e.g. the heart, arteries, bladder and brain tube) in the developing embryo. While the internal pressure rises, the circumferential stress in the vessel wall increases. As shown later, the HR hypothesis predicts that the vessel radius would grow unrealistically without bound, even after the pressure reaches a steady state, while the target stress becomes negative at maturity. A second example is epithelial wound healing, which is considered by some to be a good model for normal epithelial fusion in the embryo (Davies 2005). A wound in an embryonic epithelium closes (without scarring) by the ‘purse-string’ action of an actomyosin ring that forms and contracts around the edge of the wound (Jacinto *et al.* 2001). However, when a hole is cut into a membrane under tension, the tangential stress increases near the wound, and an HR response would cause circumferential growth to reduce this stress, thereby enlarging the wound further.

One difference between these problems is how fast changes in load occur. Blood pressure increases relatively slowly (order of days), morphogenesis takes place over a period of hours, and a wound can occur in seconds or less. From these and other considerations, we speculate that different responses are elicited by differences in the rate of change in tissue stress. Accordingly, we modify the HR hypothesis and propose a unified quantitative law for morphomechanics. Computational models are used to show that this law is consistent with the observed behaviour of growing pressure vessels and healing wounds. Other examples, including stretching of an axon and sea urchin gastrulation, then are used to illustrate the generality and some of the limitations of our theory.

## 2. Preliminaries

Unquestionably, genes play a key role in regulating development. It is not clear, however, how much of an embryo’s life they control directly. Do genes dictate every action of every cell at each instant of time? Or do cells have some degree of autonomy as they respond and adapt to environmental cues? Here, we leave this debate for another time and assume that cells respond, either directly or through genetic signals, to environmental stimuli according to a set of feedback rules. In the present study, we consider only the mechanical environment and refer to these rules as ‘morphomechanical laws’. Our task is to determine the mathematical form of these laws.

An embryo does not simply grow into its final form. In fact, many morphogenetic processes involve *deformation* of tissues by mechanical and chemical forces, with growth (changes in tissue volume) not contributing substantially to these changes in shape. On the mechanical side, both passive and active forces cause changes in cell and tissue shape. Passive forces include pressure via swelling of extracellular matrix, as well as pushing and pulling between neighbouring tissues. Active forces include those generated by differential growth and contraction of the cytoskeleton.

The present study deals only with tissue-level behaviour, and the analysis is phenomenological in nature. One disadvantage with this approach is that molecular mechanisms cannot be studied directly. Clearly, gene expression, molecular signalling and detailed cytoskeletal activity provide the biological foundation for any macro-scale law. Another disadvantage is that different morphogenetic processes can produce similar macroscopic effects (Brodland & Clausi 1994). For example, multiple mechanisms (e.g. growth, cell intercalation, active cell shape change and contraction/relaxation) can elongate or shorten tissues. Clearly, some macroscopic features can help sort out the various effects, e.g. the stiffening that accompanies contraction and the volume changes that are associated with growth. For simplicity, however, we herein lump these mechanisms together and treat them mathematically as special cases of volumetric growth, e.g. contraction is modelled as negative growth along contractile fibres. But, for convenience, we distinguish between active lengthening and shortening by referring to elongation of the zero-stress length of a material element as ‘growth’ and shortening as ‘contraction’.

## 3. Continuum theory for morphomechanics

Embryonic development involves growth (change in volume), remodelling (change in material properties) and morphogenesis (change in shape) (Taber 1995). While a comprehensive theory for morphomechanics must include all of these processes, the present study focuses primarily on morphogenesis, with all intrinsic force-generating mechanisms modelled by anisotropic growth.

More than 25 years ago, Skalak (1981) laid out the fundamental ideas that underlie many of the modern continuum theories for soft-tissue growth. These ideas have since been formalized mathematically by Tozeren & Skalak (1988) for small deformation and by Rodriguez *et al.* (1994) for large deformation. The present analysis is based on the latter theory. Here, we consider pseudoelastic behaviour, including large strains and nonlinear material properties, with viscoelastic effects neglected.

The basic theory is described briefly below. Further details can be found in previous reports on the mechanics of growth and morphogenesis (Rodriguez *et al.* 1994; Taber 2001). A new aspect of the present work is the form of the morphomechanical laws.

### (a) Kinematic and equilibrium equations

As suggested by Rodriguez *et al.* (1994), the kinematic description of growth is perhaps best visualized by considering a series of virtual configurations (figure 1). At the onset of development, we assume that the tissue is in a stress-free reference configuration *B*. Suppose now that *B* is cut into a set of infinitesimal elements, which then grow according to the growth tensor **G**. With **G** being uniform in each element, the elements remain stress free, and this growth yields the current zero-stress state *B**. Next, the elements are re-assembled into the configuration *B*_{R}. After growth, the elements generally are no longer geometrically compatible, and this assembly requires deformation that produces residual stress. Finally, external loads applied to *B*_{R} give the current loaded state *b*.

The total deformation gradient tensor **F** defines the mapping from the reference configuration *B* to the current configuration *b*. The elastic deformation gradient tensor **F*** defines the deformation of *b* relative to the current zero-stress state *B**. These tensors are related by the equation (Taber 2001),
3.1

Special circumstances can constrain the components of **F*** and **G**. For instance, material incompressibility constrains the elastic deformation through the relation . But if an incompressible material grows, then for added volume and for subtracted volume. On the other hand, if a tissue changes shape with no change in volume, then . Unless stated otherwise, all of the models in this paper are based on morphogenetic shape changes in an incompressible material, i.e. .

With inertial effects being negligible, morphogenesis can be treated as quasi-static. In terms of the Cauchy stress tensor **σ**, the equilibrium equation is
3.2
where **∇** is the gradient operator defined in *b*.

### (b) Mechanical constitutive relations

In this theory, stress is associated only with the elastic deformation given by **F***. For an incompressible material, the constitutive relation can be written in the form (Taber 2001)
3.3
where *p* is a Lagrange multiplier used to enforce incompressibility, *W*(**E***) is the strain-energy density function, **I** is the identity tensor and is the Lagrangian strain tensor for the elastic part of the deformation. For simplicity, the strain-energy density function is taken in the neo-Hookean form
3.4
where is the first invariant of the right Cauchy–Green deformation tensor defined by **C***=**F***^{T}⋅**F*** and *C* is a material constant.

### (c) Morphomechanical laws

In the present formulation, growth is a tensor quantity defined relative to the reference configuration *B*. In two or three dimensions, it is convenient to assume that growth occurs only along locally orthogonal directions, which we take as the principal material directions in *B* or, for isotropic materials, the principal directions of stress. Choosing principal material directions as a basis for growth is consistent, for example, with muscle fibres growing longer and thicker within an isotropic matrix.

Accordingly, growth occurs locally along three mutually orthogonal unit vectors **N**_{i} (*i*=1,2,3) in *B*. In addition, it is convenient to write the morphomechanical laws in terms of the rate-of-growth tensor defined by
3.5
where dot denotes time differentiation. Thus, we write
3.6
With written in the same form, these equations give , which represents the rate of change in length (due to growth) per unit length of an infinitesimal element in the current zero-stress state.

According to Beloussov (1998), embryonic tissues respond to perturbations in stress, rather than strain or some other mechanical quantity. Here, we adopt this as a key assumption (for supporting arguments, see Taber 2008). So, the problem boils down to determining the functional form of a morphomechanical law in which the rate-of-growth tensor depends on the Cauchy stress tensor, i.e. **G**° = **G**°(**σ**).

Finding such a law is analogous to determining the form of a mechanical constitutive law for a nonlinear material. Hence, it is possible that the law for each type of tissue may have a different form. One strategy is to assume a functional form, determine the parameters experimentally, and then test the law using independent experiments (Humphrey 2002). As the laws of physics generally have relatively simple forms, it is best to begin with the simplest form possible, to be modified later as data dictate. Unfortunately, published quantitative data on morphomechanics are sparse, and the need for temporal data makes this problem particularly challenging.

Following our previous study of the HR hypothesis (Taber 2008), we assume that the morphomechanical law has the relatively simple form
3.7
where **A** is a fourth-order tensor of morphomechanical coefficients, and **σ**_{0} the homeostatic or ‘target stress’ tensor. (Morphogenetic equilibrium occurs when **σ**=**σ**_{0}.) In general, it is convenient to express the stress tensors in terms of convected coordinates, which rotate with the material (Taber 2004).

Growth laws of the form of equation (3.7) have been used in prior models for the cardiovascular system (Rachev *et al.* 1996; Taber 1998; Rodriguez *et al.* 2007; Alford *et al.* 2008). Usually, the target stress is taken as constant, which seems to work reasonably well for functional adaptation of mature tissues. For developing tissues, however, this assumption can be problematic. For example, consider a pressure vessel whose growth is driven by wall stress that begins at zero and increases with the pressure. If the target stress begins at the mature value and the growth coefficients are positive, then the vessel would atrophy initially. One possible solution is that the target stress also begins at zero and increases during development (Taber 1998).

A wandering target stress is one way to generate a stress overshoot (HR) in a feedback control system. Another way is for the system to contain a delayed response. In this case, the target stress would remain constant, and the tissue stress would exhibit damped oscillations about the target stress before finally settling in on its equilibrium value. In this case, morphogenesis would then cease and await a new perturbation. In contrast, a permanent or extended stress overshoot due to a changing target stress would further perturb adjacent tissue, which then would respond in turn to propagate morphogenesis (Taber 2008). Beloussov (1998) has speculated that such behaviour causes a succession of developmental events to create form. Here, we adopt this idea and, as in a previous study (Taber 2008), speculate that the target stress evolves with time.

To determine an evolution equation for the target stress, we assume that this stress changes with the zero-stress configuration, as defined by **G**. Hence, , and we write
3.8
which has the same form as equation (3.7), with **B** being a fourth-order tensor of morphogenetic coefficients.

Together, equations (3.7) and (3.8) represent our first-approximation morphomechanical laws. As shown later, the tensors **A** and **B** play key roles in extending the HR hypothesis to a more general theory.

It is important to mention that, although **G** can be anisotropic, the material properties are taken as isotropic for simplicity. There really is no reason that an initially isotropic tissue cannot grow different amounts in different directions in response to an anisotropic stress field. Although this simplification should not influence the fundamental behaviour of our models, future work should address the effects of material anisotropy.

### (d) Simplified equations in principal coordinates

For problems with appropriate symmetry, the principal directions of stress and strain align with the principal material directions. In this case, we can write
3.9
and so on, with the λ’s being stretch ratios. The strain invariant in equation (3.4) becomes
3.10
relative to the zero-stress state. In addition, we assume for simplicity that all of the off-diagonal terms of **A** and **B** are zero. With these simplifications, the governing equations are specialized for each example as needed.

## 4. Fundamental behaviour in one dimension

Previously, we used specialized forms of equations (3.7) and (3.8) to investigate the HR response of some basic problems in morphogenesis (Taber 2008). The main idea of the HR hypothesis is illustrated by the following one-dimensional example. Suppose an elastic bar is initially in a stress-free homeostatic state, i.e. the axial stress σ and target stress σ* are zero everywhere. Then, the bar is stretched and held at a fixed length. According to the HR hypothesis, the bar responds by growing longer to lower the tension until a new equilibrium is achieved with the bar in compression. Hence, the bar grows more than that needed to restore a stress-free state, and the response overshoots the original stress of zero to the other side, with the initial tension becoming compression. Similarly, if the bar is shortened and held, the bar contracts, transforming compression into tension. This example is used below to demonstrate the implications of equations (3.7) and (3.8).

### (a) Basic equations for stretching of a bar

For uniaxial stretching of a bar composed of neo-Hookean material, incompressibility, symmetry, the boundary conditions and equations (3.3) and (3.4) yield (Taber 2004)
4.1
where σ and λ* are the axial stress and elastic stretch ratio, respectively. In addition, equation (3.1) yields the total stretch ratio
4.2
in which *G* is a growth ‘stretch ratio’ to be found from the morphomechanical laws
4.3
and
4.4
as provided by the one-dimensional forms of equations (3.7) and (3.8).

The solution to this problem was computed using Matlab for both isometric and isotonic loading. In the simulations, the initial conditions are taken as *G*=1 and σ_{0}=0. For the isometric case (λ is a specified constant), equations (4.1) and (4.2) provide the stress, and equations (4.3) and (4.4) are integrated using forward differences to obtain *G* and σ_{0} for the next time step. Then, σ is computed at the next time step, and so on.

For isotonic loading (the first Piola–Kirchhoff stress *P* is a specified constant), the Cauchy stress is given by σ=λ*P*=*G*λ**P*, where constant tissue volume is assumed With *G* computed as above at each time step, this expression and equation (4.1) are solved for λ*, and then equation (4.2) gives the total stretch ratio.

### (b) Unified morphomechanical law in one dimension

The signs of the morphomechanical coefficients *A* and *B* define the qualitative behaviour of models based on this theory. For HR, increased stress induces growth, while decreased stress stimulates contraction. This response requires *A*>0. For an oppositely directed overshoot, the target stress must move in the opposite direction, giving *B*<0.

The HR response for this and other problems (with *A*>0 and *B*<0) was explored previously (Taber 2008). As mentioned in §1, however, some problems in development apparently are not consistent with the HR hypothesis. One example cited was growth of a pressure vessel, where HR predicts unbounded growth and an unrealistic negative target stress at maturity. This contradiction can be resolved by taking *B*>0 (along with *A*>0). The other example mentioned was wound healing, for which an HR response would open the wound further. The correct behaviour requires *A*<0, so that the elevated tension parallel to the edge of the wound elicits a contractile response. Because the contraction would further elevate the tension, a new equilibrium would require *B*>0 also in this case.

We speculate that these different responses are caused by differences in the rate of change in stress. In the embryo, growth occurs in response to loads (pressure, flow, etc.) that normally change relatively slowly, over a period of days. Moreover, embryonic growth occurs primarily by cell division, which may take a day or so per cycle. Morphogenesis is generally a faster process, with significant changes in form occurring on a scale of hours. Finally, wound healing is a response to a stress perturbation that occurs in seconds or less.

From these considerations, we assume that the same morphomechanical law applies to all of these problems, but the coefficients in the law depend on stress rate. The functional form of this dependence is left for future work. Here, as a first approximation, the following response classes are defined:

### (c) Results for stretching of a bar

Before showing how this theory resolves the issues raised above, we illustrate the basic behaviour for uniaxial loading of a bar. For this purpose, it is useful to write the governing equations in non-dimensional form. With *T* being a characteristic time, the following dimensionless variables are introduced:
4.5
With these definitions, the basic forms of the governing equations are essentially unchanged if bars are placed over all variables, and dot indicates differentiation with respect to . For convenience, overbars are omitted throughout the remainder of the paper, with the understanding that all referenced variables are non-dimensional, including those contained in plots.

With either the length of the bar or the applied load fixed, the only time scales are associated with the rates of change of *G* and σ_{0}, as characterized by *A* and *B*, respectively. Hence, the qualitative temporal response depends only on the ratio *B*/*A*, and, without loss of generality, we take *A*=1 or −1 and change the value of *B*. The specific values of *B* (given in the figure captions) are chosen for clarity of presentation.

For a step stretch (λ=1.2), illustrative results are compared for the three -dependent response classes defined above (figure 2). For morphogenetic HR, the bar grows following the stretch, while both the stress σ and the target stress σ_{0} drop from their initial values. Eventually, σ catches σ_{0}, and growth stops with the bar in compression, in accordance with the basic behaviour described above for the HR hypothesis. For growth (GR), both *G* and σ_{0} increase while σ decreases but remains tensile. This behaviour resembles stress relaxation for a viscoelastic solid, with the bar achieving a new equilibrium at a lower tension. Finally, stretch activation (SA) induces contraction which increases the tension. In this example, the target stress increases fast enough to eventually catch the stress, yielding a new equilibrium at an elevated tension.

A step shortening of the bar (λ=0.8) induces oppositely directed responses, with the bar ending up in tension for HR and compression for GR and SA (figure 3). Notably, the time course of the HR stress response differs qualitatively from the case of initial tension (compare figures 2*b* and 3*b*), suggesting a way to test these ideas experimentally.

Clearly, these results depend on the ratio *B*/*A*. Illustrative results for different values of *B* (with *A*=1) are shown for the HR case (figure 4). The magnitude of the overshoot increases with the magnitude of *B*, but a new equilibrium always is established.

Under isotonic loading, equilibrium demands that the bar remains in tension or compression, depending on the sign of the load. Hence, overshoot to the opposite direction cannot occur, and the HR response is inherently unstable, i.e. the bar grows without bound (figure 5). In contrast, the GR and SA responses reach new equilibria at increased and decreased lengths of the bar, respectively.

## 5. Resolution of problems inconsistent with the HR hypothesis

Now, we turn to the two aforementioned problems that apparently are inconsistent with the HR hypothesis—growth of a pressure vessel and epithelial wound healing. Finite element models for these problems were developed using Comsol Multiphysics (v. 3.4; Comsol AB, Burlington, MA, USA), with the morphomechanical laws included as auxiliary differential equations. Details on the implementation of volumetric growth in Comsol can be found in Taber (2008).

### (a) Growing pressure vessel

Consider a thick-walled neo-Hookean cylinder of fixed length with an outer radius 1 and inner radius 0.8 (figure 6, inset). The cylinder is subjected to an internal pressure *p*_{i} that increases linearly in time from zero at *t*=0 to a peak value *p*_{i}/*C*=0.1 at *t*=0.1 and remains constant thereafter. In the heart, for example, *t*=0 would represent the time that the heart tube first forms in the embryo, and *t*=0.1 would be the time of birth. The pressure could be the mean, peak systolic, or pulse pressure. For *t*≥0, the cylinder grows according to equations (4.3) and (4.4), where *G*, σ and σ_{0} represent quantities in the circumferential direction. Radial and axial growth, as well as axial stretch, are not included in this illustrative example. The initial conditions are *G*=1 and σ_{0}=0 everywhere.

For the slowly increasing pressure that occurs during development, the response is treated as stress-induced growth with *A*>0 and *B*>0. Results are shown for *A*=1 and different values of *B* (figure 6). As the pressure increases for *t*≤0.1, the lumen grows at an accelerating pace; after the pressure plateaus, the behaviour depends on the value of *B* (figure 6*a*). If *B* is large enough, the radius continues to increase to a new homeostatic value. At equilibrium, although the stress is not uniform, σ=σ_{0} throughout the wall (figure 6*b*). For *B*=1, however, unbounded growth occurs.

This behaviour can be explained by considering the average circumferential stress , as given by Laplace’s law, where *r*_{i} and *h* are the deformed inner radius and wall thickness, respectively. With the pressure fixed (*t*≥0.1), circumferential growth continues as long as σ>σ_{0}, increasing the radius and . If *B* is large enough, the target stress increases fast enough to eventually catch the rising stress, and growth stops. If *B* is too small, however, the tube continues to grow without bound.

These results show that, in contrast to the HR response, realistic stable growth occurs for values of *B*/*A* large enough to establish a homeostatic radius at maturity. Confirmation of this prediction awaits experimental measurements of *A* and *B*.

### (b) Wound healing

The model for epithelial wound healing is a square neo-Hookean membrane that, in the initial unloaded configuration, has sides 2 units in length and a hole of radius 0.1 at the centre (figure 7, inset). A plane-stress finite-element model was constructed for one quadrant with appropriate symmetry conditions enforced. Quadrilateral elements are used with the mesh being finer near the hole, where a stress concentration develops (figure 7*b*). The membrane is given an initial tension by stretching it 5 per cent equibiaxially. (The amount of stretch has little effect on the qualitative behaviour.) Then, the edges are constrained by roller supports, and the material is assumed to grow/contract according to the equations
5.1
and
5.2
where *r* and θ are cylindrical polar coordinates in the plane of the membrane. For isovolumic shape change, the transverse normal growth is *G*_{z}=1/*G*_{r}*G*_{θ}.

For the growth tensor, the initial conditions are *G*_{r}=*G*_{θ}=1 at all points. For the target stresses, the membrane *without a hole* is subjected to 5 per cent stretch, and the computed stress σ_{r}=σ_{θ} provides the initial values of σ_{r0} and σ_{θ0} in the wounded membrane. (The hole is small enough relative to the size of the membrane so that the disturbance caused by the hole decays to nearly zero by the time it reaches the boundaries of the membrane.) In addition, tissue stiffening that accompanies contraction is simulated by setting the material coefficient *C*=1/*G*_{θ}.

The sudden change in stress caused by the wound is assumed to stimulate a stretch-activation response with *A*<0 and *B*>0. Illustrative results are shown for *A*=−1 and various values of *B* (figure 7). Consider first the solid curves for the hole radius (figure 7*a*), which include growth/contraction in both the *r* and θ directions. As expected, circumferential contraction (*G*_{θ}<1, figure 7*c*) in response to the suddenly elevated circumferential tension (figure 7*b*) closes the wound (figure 7*a*). For all values of *B* examined, the general behaviour is similar qualitatively. In particular, wound closure exhibits a biphasic response, with a relatively slow phase followed by a more rapid phase (figure 7*a*). As *B* increases, the duration of the first phase increases. but the wound eventually heals for all values of this parameter. Finally, radial growth has a relatively small effect on healing (dashed curve in figure 7*a* for *G*_{θ} only and *B*=0).

Temporal measurements of wound closure are needed to test these model predictions, as well as to determine the values of *A* and *B*. Currently, some qualitative observations are available for comparison. First, within minutes after an embryonic epithelium is wounded, a narrow contractile ring (‘purse string’) begins to form around the edge of the wound (Jacinto *et al.* 2001; Davies 2005). The localized drop in *G*_{θ} near the edge of the hole is consistent with this finding (figure 7*c*), although the actin ring may actually be narrower than the model predicts. A more refined model that includes a detailed structure of the actin cytoskeleton could improve the accuracy. Second, the increased membrane thickness (λ_{z}>1), due to incompressibility, near the wound edge also agrees with experimental observations (Jacinto *et al.* 2001). On the other hand, measurements have shown that the similar process of dorsal closure in the *Drosophila* embryo proceeds at a relatively steady rate (Hutson *et al.* 2003; Toyama *et al.* 2008), in contrast to the biphasic response predicted by the model. Interestingly, taking *C* as a constant gives a more uniform rate of closure (dash-dot curve in figure 7*a* for *C*=1). To evaluate our model more fully, careful measurements of material properties are needed during wound closure.

In summary, these results agree with much of the observed behaviour of growing pressure vessels and healing wounds. Thus, a theory of morphomechanics that includes the effects of stress rate appears to resolve contradictions raised by applying the HR hypothesis to these problems.

## 6. Counterexample: rapid shortening of an axon

In searching for laws of nature, it is just as important to seek counterexamples as it is to find supporting data. Here, we consider one problem that appears to be inconsistent with our theory. There probably are others.

As the nervous system develops, neurons send out thin, tubular axons to form connections with other cells. A specialized structure, called the growth cone, crawls away from the cell body, stretching the axon and carrying its distal end to a target. Axons can reach a metre or more in length, indicating that the elongation involves growth induced by axial tension (Lamoureux *et al.* 1989; Gilbert 2003). Using flexible needles, Dennerll *et al.* (1989) measured the tension in axons from dorsal root ganglion neurons of chick embryos. After attaching one end of a needle to the growth cone, they lifted it from the substrate and held the axon at a relatively fixed length until the force equilibrated. (Force was determined from needle deflection.) When the axon was subsequently stretched and held, the force slowly dropped but remained tensile. However, when the axon was rapidly shortened, the force dropped suddenly and then increased to a new equilibrium value that, depending on the neuron, was either below or above its initial value.

This experiment is simulated using the above model for stretching a bar, which represents the axon between the cell body and growth cone. In their experiments, Dennerll *et al.* (1989) failed to measure a number of key variables, e.g. stretch ratios, forcing us to make a number of assumptions. First, assuming that the initial value of σ_{0} is given by the stress in the cell membrane before the axon forms, we obtain a representative value for σ_{0}(0) by stretching the bar by 10 per cent at *t*=0 and computing the corresponding stress. Next, prior to the beginning of the experiment, axon stretching by a growth cone is simulated by stretching the bar at a constant rate until *t*=0.2, after which the length is fixed at λ=1.4 (see solid curve for λ in figure 8*a*). Finally, from a newly achieved equilibrium at *t*=1, the bar is shortened and maintained at λ=1.35 thereafter. During the entire simulation, the bar grows/contracts according to equations (4.3) and (4.4) with *G*(0)=1 and the values of *A* and *B* given below. For this problem, it is instructive to consider two phases—the stretch and hold phase (phase 1, *t*≤1), and the shorten and hold phase (phase 2, *t*>1).

During phase 1, the model must be able to capture two key results from the experiments of Dennerll *et al.* (1989). First, the axon grows longer in response to tension applied by the growth cone. Growth cones can move tens of microns per hour, falling in the range of morphogenesis (HR).1 Because tension-induced growth also is consistent with an HR response with *A*>0 and *B*<0, taking *A*=1 and *B*=−1 for phase 1 gives the correct qualitative behaviour (see dashed curve in figure 8*a* for *t*<1).

Second, after the axon is stretched and held, the stress drops but remains tensile. This finding leads to a relatively subtle but important point. If the target stress begins at zero, then HR would transform the initial tension into compression. However, by assuming that the target stress is initially positive (equal to the tension in the cell membrane), then an HR response lowers the stress, but tension can remain when a new equilibrium state is achieved (see solid curve in figure 8*b* for *t*<1).

Thus far, everything seems consistent with our theory, but phase 2 is a different story. When axon stretch is decreased rapidly in the experiments, the tension increases to a new value *above or below* the original tension. It is important to note here that this phase apparently begins from a homeostatic state with σ=σ_{0}. Hence, the shortening drops σ below σ_{0}, i.e. σ−σ_{0}<0. According to our theory, the sudden change in loading should elicit an SA response with *A*<0. But this would lead to more growth and a further drop in stress (figure 8, phase 2 curves SA), contrary to the experimental finding of increased tension.

On the other hand, HR gives a response closer to reality, with the stress increasing to a new equilibrium (figure 8*b*). However, HR always predicts an overshoot of the stress present at the end of phase 1, whereas Dennerll *et al.* (1989) found an undershoot in some neurons.2

These results force us to re-evaluate our theory. They may indicate, as some have suggested, that searching for a universal law for morphomechanics is a fruitless endeavour (Forgacs & Newman 2005; Nanjundiah 2005). This is certainly possible. Before giving up, however, let us consider two other possibilities. First, it is possible that the axon grows thinner as it responds to a change in tension. (Dimensional changes were not measured in the experiments.) In this case, the force would decrease for a given value of Cauchy stress, potentially reaching a new equilibrium below, rather than above, its value at the end of phase 1. Second, the stress-rate dependence of *A* and *B* was inferred primarily from experiments and observations for embryonic epithelia that generally have not yet differentiated into specific tissues. In contrast, axons grow from cells that already have differentiated into neurons. Moreover, the microstructures of axons and epithelial cells differ considerably. For example, while the actin cytoskeleton dominates the mechanical behaviour of epithelial cells, microtubules play a major role in force generation in axons (Ahmad *et al.* 2000; Baas & Ahmad 2001). Hence, we speculate that the morphomechanical laws of equations (3.7) and (3.8) still apply, but the rate dependence of the morphomechanical coefficients may be tissue specific.

In this regard, it could be that axons are extremely sensitive to loading rate. Small differences in the experimental rate of shortening could lead to different types of behaviour that produce an overshoot or undershoot of the initial stress. (Stress relaxation due to viscoelasticity also is consistent with an undershoot.) Presumably, these effects can be included in the functions that define how the morphomechanical coefficients depend on stress rate. Answering these questions, however, may require experiments and analysis to determine the underlying molecular mechanisms of this behaviour.

## 7. Application to sea urchin gastrulation

In this section, we apply our theory for morphomechanics to the problem of gastrulation, which is a process that establishes the primitive gut and reorganizes the primary germ layers in the early embryo. Owing to its relatively simple geometry, the sea urchin embryo is a popular organism for studying gastrulation. The first series of cell divisions create a fluid-filled spherical shell called the blastula (figure 9a). Soon, a localized dimple (invagination) forms and deepens until it extends the entire diameter of the embryo to form the gut tube (figure 9*b*,*c*). Experiments and models have provided insights into the mechanisms of gastrulation (Davidson *et al.* 1995, 1999), but the main driving forces remain incompletely understood. Here, we examine the possibility that mechanical feedback plays a role in this process.

### (a) Normal embryo

Recently, we proposed a first-approximation shell model for sea urchin gastrulation based on the HR hypothesis (Taber 2008). An initial perturbation consisting of a localized 10 per cent contraction produced a dimple that continued to deepen after the perturbation ended, but the depth of the dimple was quite limited. That model is modified here in two ways.

First, the model is made more realistic by including two layers—an inner layer of epithelial cells and a stiffer outer layer of extracellular matrix representing a combination of the apical lamina and hyaline layer in the sea urchin blastula (see Davidson *et al.* 1995). Studies have shown that the matrix probably plays a role in gastrulation by swelling or deswelling (Lane *et al.* 1993; Davidson *et al.* 1999). Second the initial perturbation is increased in magnitude to a 50 per cent local contraction.

The model for the blastula is a spherical bilayered shell filled with an incompressible fluid (figure 9*d*, open outline). A quasi-static axisymmetric analysis was used to solve this problem in Comsol, with the finite element mesh consisting of triangular elements that are denser in and near the invaginating region. The effects of the fluid are included through a weak constraint condition that maintains a constant cavity volume; a Lagrange multiplier represents the hydrostatic fluid pressure. To prevent rigid-body motion, a point at the centre of the dimple is fixed.

All variables are scaled according to equation (4.5) with *C* being the modulus of the inner layer, and dimensions are scaled by the inner radius. In the reference state at *t*=0, the shell has an inner radius 1 (unit), an outer radius 1.10 and the boundary between layers is located at a radius of 1.07. The material coefficients are taken as *C*=1 for the inner layer and *C*=5 for the outer layer. These relative values are representative of the sea urchin embryo (Davidson *et al.* 1995).

Let (*R*,Φ,Θ) be radial, meridional and circumferential spherical polar coordinates of a material point in the reference configuration. In the deformed configuration (*b* in figure 1), the coordinates become (*r*,ϕ,θ). The active tissue response is assumed to occur parallel to the middle surface of the shell and is governed by the relations
7.1
and
7.2
in which the stresses are written relative to convected coordinates. The parameter α represents a specified rate of growth (or contraction) that is independent of stress. In addition, for cell shape change, *G*_{R}=1/*G*_{Φ}*G*_{Θ}.

For large invaginations, concentrations in stress and strain near the edge of the dimple cause convergence issues. Relaxing the constraint of material incompressibility helps, and thus the strain-energy density function is taken in the modified neo-Hookean form 7.3 where κ is the bulk modulus, and . The value of κ is taken as 10 in both layers. Changing the value of κ somewhat alters the results quantitatively, but not qualitatively (Taber 2008), so the effects on the present qualitative study should not be significant.

The simulation begins with a stress-free state with *G*_{R}=*G*_{Φ}=*G*_{Θ}=1 and σ_{ϕ0}=σ_{θ0}=0 everywhere. Invagination is initiated by imposing an isotropic in-plane contraction in the outer layer of the dimple region3 (,figure 9*d*) by setting α=−5 for 0≤*t*≤0.1 (peak contraction of 50%) and α=0 thereafter. The response of the entire blastula wall is governed by equations (7.1) and (7.2) for all *t*≥0. For this morphogenetic process, an HR response is assumed with *A*=1 and *B*=−1.

The shapes predicted by the model agree reasonably well with experimental observations (compare figure 9*b* and *f*, *c* and *h*). After the prescribed initial contraction reaches a steady state at *t*=0.1, the invagination continues to deepen. Note also the considerable early thickening of the dimple region (figure 9*b*,*f*) and the sharp corner at the edge of the dimple (figure 9*c*,*h*). If the specified initial contraction is smaller, subsequent morphogenesis becomes slower but remains similar qualitatively (figure 10*a*).

During the initial phase of gastrulation, the displaced fluid stretches the shell outside the region of invagination, and the fluid pressure increases (figure 10*b*). The walls then grow in response to the tension, and the HR overshoot eventually produces compressive stresses and a negative cavity pressure. Consistent with this result, blastulae in a number of species are wrinkled in appearance, but the wrinkles depend on the osmolarity of the surrounding medium (Kobayakawa & Satoh 1978). If the blastula wall is under tension before gastrulation begins, then the final wall stress and pressure may remain positive, but the HR hypothesis clearly predicts a drop in pressure as gastrulation proceeds. Experimental pressure measurements during gastrulation are needed to test this prediction.

### (b) Mechanically perturbed embryo

Experimental perturbations of morphogenesis offer data that can be used to test the predictive capability of mathematical models. In an intriguing study, Moore and Burt (1939) divided gastrulating sea urchin embryos into two parts by cutting them circumferentially. They then observed the behaviour of the part containing the invagination. Following the cut, the dimple maintained its shape, but the region outside the dimple first unfurled and became almost straight before reversing course and bending back to form a complete, but smaller, gastrula (figure 11*a*).

Simulating this experiment with our finite element model met additional convergence problems caused by highly deformed elements near the dimple edge. We were able to push the solution a little further by ignoring radial growth, i.e. we set λ_{gR}=1. Because the main concern here is what happens outside the dimple, where relatively little change in wall thickness occurs, this simplification should have little effect on the main results.

To simulate the experiment, the solution for the full gastrulating shell is saved at *t*=1 (figure 11*b*). Next, all equations are deactivated in the top half of the sphere, effectively removing this part and leaving a free edge at the top of the remaining lower hemisphere. Then, a static solution is obtained for the cut shell with all growth held fixed (figure 11*c*). Finally, the time-dependent calculations are restarted with the fluid constraint removed and growth turned on everywhere except inside the dimple. In addition, as cutting the embryo suddenly removes the internal fluid pressure and the wall stresses along the edge of the cut, the morphomechanical coefficients are changed to *A*=−1 and *B*=1 for an SA response (figure 11*d*).

Until the solution procedure diverges, the model exhibits trends that appear to agree with the general behaviour of the cut embryo (figure 11*a*). In particular, the bulk of the embryo outside the dimple bends outwards as the free edge curls inwards (figure 11*d*). While these tendencies resemble the experimental result, the corresponding deformations in the experiments did not occur simultaneously and were considerably more striking (figure 11*a*). This discrepancy could be due to the simplicity of the model or the inadequacy of our theory. Another possibility is indicated by running the model with *A* and *B* remaining at their HR values (*A*=1 and *B*=−1). In this case, the embryo outside the dimple opens substantially more, but the edge does not bend inwards (result not shown).

Taken together, these results suggest that the morphomechanical coefficients may take some time to change even when the stress changes suddenly. In other words, the ‘momentum’ of the HR response takes time to convert into an SA response. To explore this idea, the simulation was run for a gradual transition from HR to SA behaviour After the cut, the HR response first opens the embryo, and then SA closes it again (figure 11*e*). However, this model really is too simple to draw definitive conclusions.

## 8. Discussion

Physical forces drive the deformations that occur during morphogenesis, and researchers have speculated that mechanical feedback plays a central role in regulating these forces (Odell *et al.* 1981; Beloussov *et al.* 1994; Beloussov 1998; Ramasubramanian & Taber 2008). This idea has inspired the development of computational models for a number of processes, including gastrulation (Odell *et al.* 1981; Taber 2008), neurulation (Clausi & Brodland 1993; Brodland & Clausi 1994), cardiac looping (Nerurkar *et al.* 2006; Ramasubramanian *et al.* 2008), and bone development (Carter 1987). The feedback laws used in these models are based on stretch-activated cytoskeletal contraction, HR or stress-induced changes in density and material properties (for bone). While changes in form and structure given by these models agree reasonably well with experimental observations for specific applications, there is no consensus on the specific form that a biomechanical feedback law should take. Beloussov (1998) has spent more than three decades studying the mechanics of morphogenesis. His experiments have revealed relatively consistent patterns in the active response of embryonic tissues to perturbations in stress. Integrating these observations led to the HR hypothesis, which has been used to explain in qualitative terms a number of simple, as well as relatively complex, morphogenetic processes (Beloussov 1998, 2008; Beloussov *et al.*1994, 2006). Moreover, mathematical models based on this idea have yielded some promising results (Belintsev *et al.* 1987; Beloussov & Grabovsky 2006; Taber 2008; Ramasubramanian *et al.* 2008).

As discussed in §1, however, the HR hypothesis is not consistent with all problems in development. We considered two such problems in this report—growth of a pressure vessel and healing of an epithelial wound. To address this issue, we proposed an extension of the HR hypothesis to formulate a more general theory for morphomechanics. In our theory, growth is used to simulate any process that actively changes tissue shape, with evolution equations governing the rates of change in the growth and target stress tensors in response to perturbations in stress. Letting the coefficients in these equations depend on stress rate allows our models to capture observed behaviour that is not consistent with an HR response.

Models based on our theory yield the correct behaviour for a growing pressurized cylinder and a healing epithelial wound, as well as for normal sea urchin gastrulation. However, models for axon shortening and experimentally perturbed gastrulation indicate that there is more work to be done. (A more realistic model for gastrulation based on the research of Davidson *et al.* (1995) may help address some of the outstanding issues.)

The results presented in this report are primarily for illustrative purposes. Our intent is to stimulate further discussion of this topic—we make no claim that the theory embodied in equations (3.7) and (3.8) is correct. Nor do we make such a claim for the proposed morphogenetic mechanisms contained in models for specific problems. It is possible, for example, that growth rate actually is a nonlinear function of stress or strain, and contraction may not initiate gastrulation.

Being phenomenological in nature, the present theory ignores its genetic and molecular underpinnings. Moreover, to avoid an exercise in curve fitting, the number of unknown parameters is kept to a minimum, i.e. the temporal response of each model depends only on the ratio *B*/*A*. It is important to note the relatively robust qualitative response indicated by the parameter studies, and the axon example shows that equations (3.7) and (3.8) are not so general that they can reproduce all possible morphogenetic behaviour.

The present study also sidesteps the issue of constraining the morphomechanical coefficient tensors **A** and **B** for general three-dimensional problems. Here, we have assumed that growth rates along a particular material direction are affected only by the normal stress in the corresponding direction. It is possible, and perhaps even likely, that stresses in orthogonal directions also have an effect. The signs of the various components of these tensors, as well as other mathematical properties, are left to future theoretical and experimental work.4

In conclusion, this study represents one step in a continuing effort to determine the fundamental laws that govern how an embryo is constructed. As mentioned above, some researchers feel that biological systems do not follow mathematical laws such as those we propose (Forgacs & Newman 2005; Nanjundiah 2005). On the other hand, Beloussov argues rather persuasively that such laws may exist (Beloussov & Grabovsky 2006; Beloussov 2008). In any event, it is difficult to draw any definitive conclusions until more quantitative data are available, especially for the time-dependent *active* response of embryonic tissues to mechanical loads.

Future work should address the following questions:

Do mathematical laws for morphomechanics exist and, if so, do they depend on the specific tissue type or differentiation state?

How does the active response (and morphomechanical coefficients) of embryonic tissue depend functionally on the rates of stress and strain?

What is the value of the target stress when tissues and organs first form in the embryo and does it evolve during development? Or are there multiple stable target stresses as suggested by Fung (1990)? Moreover, is the target stress governed by an optimization principle, as suggested for arteries (Murray 1926)?

Answering these questions will require close collaboration among engineers, biophysicists and developmental biologists. A complete understanding of the principles that govern the mechanics of morphogenesis will ultimately require multi-scale analyses that link events at the molecular, cell and tissue levels. Work in this field is just beginning.

## Acknowledgements

This work was supported by grants R01 GM075200 and R01 HL083393 from the National Institutes of Health, as well as grant DMS-0540701 from the National Science Foundation. I am grateful to Pat Alford, who suggested modelling the dissected sea urchin embryo and helped with aspects of the model development.

## Footnotes

One contribution of 12 to a Theme Issue ‘Mechanics in biology: cells and tissues’.

↵1 Clearly, an axon actually grows by adding material as it extends over long distances, but an HR response is more consistent with the stress-rate criteria defined above. How the axon reaches equilibrium at the end of phase 1 is immaterial here, so long as tension is maintained.

↵2 Because force was measured in the experiments, it would be more consistent to examine the first Piola–Kirchhoff stress

*P*, which is proportional to the applied force. Calculations show, however, that the fundamental behaviour is the same, and so Cauchy stress is used for clarity in comparing with the target stress.↵3 For convenience, the term ‘dimple region’ refers only to the region containing the initial contraction, even as the actual invagination continues to spread.

↵4 Theoretical aspects of anisotropy and mechanically based growth laws are discussed in some recent papers (see Maugin & Imatani 2003; Ambrosi & Guana 2007; Menzel 2007).

- © 2009 The Royal Society