## Abstract

The shape of tissues arises from a subtle interplay between biochemical driving forces, leading to cell growth, division and extracellular matrix formation, and the physical constraints of the surrounding environment, giving rise to mechanical signals for the cells. Despite the inherent complexity of such systems, much can still be learnt by treating tissues that constantly remodel as simple fluids. In this approach, remodelling relaxes all internal stresses except for the pressure which is counterbalanced by the surface stress. Our model is used to investigate how wettable substrates influence the stability of tissue nodules. It turns out for a growing tissue nodule in free space, the model predicts only two states: either the tissue shrinks and disappears, or it keeps growing indefinitely. However, as soon as the tissue wets a substrate, stable equilibrium configurations become possible. Furthermore, by investigating more complex substrate geometries, such as tissue growing at the end of a hollow cylinder, we see features reminiscent of healing processes in long bones, such as the existence of a critical gap size above which healing does not occur. Despite its simplicity, the model may be useful in describing various aspects related to tissue growth, including biofilm formation and cancer metastases.

## 1. Introduction

Morphogenesis is controlled via the complex interaction between biochemical and mechanical signals arising from cell-generated tensions [1,2]. Such signalling is mediated through the presence of boundaries, whose shapes determine the formation of morphogen gradients [3] and the local stress state, which in turn controls cell proliferation [4], differentiation [5] and apoptosis [6,7]. Depending on the physical environment, these boundaries may either be static, as is the case for solid substrates, or they may move as new tissue is formed [8]. One important simplification that can be made to help understand this problem is based on the observation that tissues, or at least cell agglomerates, can behave like viscous fluids with measureable surface tensions when observed for sufficiently long timescales [9–12]. If one describes tissues as fluids, then the equilibrium shape of their boundaries will be determined on one hand by the wettability of any substrates upon which they are sitting [13] and on the other hand by the Laplace–Young equation giving a link between interfacial curvature and tissue pressure [14]. Interestingly, it has been shown that even simple shapes of wettable regions on flat substrates display a rich variety of equilibrium liquid droplet morphologies at constant volume when surface tension plays an important role in the energy of the system [15–17]. This suggests the possibility of a similar richness for fluid-like tissues growing under the influence of large surface tensions. As the stress state of cells is coupled to growth [18,19], an exploration of the consequence of geometric boundary conditions on surface tension mediated growth may help explain the observed influence of shape on tissue growth *in vitro* [20,21] and *in vivo* [22,23], the goal of this paper.

Many physical models for tissue growth have been studied previously (e.g. [8,24,25]) and they often use rather subtle material models, including second gradient approaches (see e.g. [26,27] and [28] for a recent overview). To simplify our model, however, we simplify the material behaviour to that of an ideal fluid (i.e. constant pressure), which relaxes local stress concentrations by viscous flow at a rate much faster than the growth that causes them. As in our case, growth dominates over elastic deformations, we can also avoid treating potential coupling between the two (e.g. [29–31]). We start with a simple approach fully compatible with thermodynamics [32–36] which considers both a biochemical driving force (in the form of a chemical potential) and a mechanical one (including volume stresses as well as surface stresses). Then we assume that local tissue remodelling is sufficiently rapid compared with the growth rate, so that the tissue can essentially be regarded as an ideal fluid on the timescale of growth. Tissue remodelling is the process by which cells can rearrange their orientation and position with respect to each other as well as to move, reorient and degrade the extracellular matrix in which they are sitting. The driving force for remodelling has its origin in the thermodynamics of the irreversible process controlled by the dissipation; for details, we refer to [36] and the contributions by Ambrosi *et al.* [24,33,34]. With these assumptions, we show that the problem reduces to a single free parameter interpreted as a critical curvature for tissue growth. This critical curvature is given by the ratio of the chemical potential, describing the tendency of the tissue to grow as a result of biochemical signals (e.g. growth factors), and the surface energy, which leads to a tissue pressure that may inhibit or accelerate growth. The appearance of surface energy in the critical curvature means that geometric boundary conditions may play an important role in growth, a major focus of this paper. For the following, we assume that the tissue fully wets a given surface and explore the role of surface shape on the resultant growth. Our problem is thus related to the evolution of droplets interacting with partially wettable surfaces [15–17], where growth corresponds to an increase of the droplet volume. We first investigate the growth of an unsupported volume of tissue (or nodule) and then its growth upon a circular disc, this analysis is then extended to tissue growth on a circular ring. Despite the simplicity of the geometries analysed, tissues grow into shapes with exquisite complexity with many aspects emulating what can be observed in callus formation and evolution during bone healing. We compare some of the predictions of the model to the kinetics of bone healing in various mammalian species and to the existence of a critical defect size above which healing is prevented. The surprising level of coincidence suggest that during bone healing, cells might indeed behave like liquids and thus profit from the self-organization capacities of wetting and surface tension.

## 2. Formulation of the model hypotheses

We use the formalism for growth outlined in [35,36] but assume that the growing body behaves like an ideal fluid, with a defined surface stress, *γ* (for details, see [14]). This relies on the assumption that there are three separable timescales, that of growth, remodelling and elasticity. The typical time for growth is taken to be much longer than the relaxation of stresses by tissue remodelling processes. If the body behaves like an ideal fluid (at long timescales) and is subjected to a constant surface stress, then the shape of the body will have a constant surface mean curvature when in mechanical equilibrium. This is formulated by the Laplace–Young equation, *p* = *γK*, where *p* is the pressure acting in the body, *γ* is the surface stress (assumed to have the same value as the surface energy) and *K* is the surface curvature (defined as the absolute value of the trace of the curvature tensor or twice the surface mean curvature; for details, see, e.g. [14]). For simple geometric boundary conditions, the equilibrium fluid shape can be calculated analytically, giving curves that relate pressure (or curvature) to volume, for more complex geometries numerical solutions such as ‘Surface Evolver’ [37] can be used.

Once the relation between pressure (or surface curvature) and volume is known for a given tissue geometry in mechanical equilibrium, then we can analyse its stability under the conditions of growth. We do this in the following using the equation of growth from [35]. Here, **g** is the growth eigenstrain tensor, ** σ** the stress tensor which is given as −

*p*

**I**, where

**I**is the unity tensor, the chemical potential tensor and

*f*a mobility coefficient. The chemical driving force for growth is a ‘chemical potential’ , with a sign chosen such that

*µ*is generally positive. The growth equation may be simplified with these assumptions to an equation depending on a chemical driving force

*μ*≥ 0, 2.1where we have defined

*V*as the volume and as its time rate. Note that the convention is that tensile stresses are positive, while the pressure

*p*is taken positive when compressive (thus the inversion of the signs in [1]). Using the Laplace–Young relation, equation (2.1) can be rewritten as 2.2It is clear from this equation that the tissue nodule will grow whenever its surface curvature

*K*is smaller than a critical curvature

*K*

_{C}=

*μ*/

*γ*. This critical value depends on the ratio of chemical potential for growth (generated by biochemical signals such as growth factors) and on the surface stress on the outer boundaries of the nodule and is interpreted as a kinetic control parameter for growth.

## 3. Growth of an unsupported tissue nodule

As a first and most simple case, we consider an unsupported tissue nodule growing in suspension. Local mechanical equilibrium implies that the nodule must be spherical with radius *R* = 2/*K*. Under these conditions, equation (2.2) can be rewritten in a form reminiscent of a nucleation problem (e.g. [38]).
3.1where *R*_{C} = 2/*K*_{C} = 2*γ*/*μ* is the critical radius. From this equation, the tissue nodule will be in equilibrium when , i.e. *R* = *R*_{C}. Whether this is stable or not, depends on how changes for small variations in *R* close to *R*_{C}. As for *R* > *R*_{C} and for *R* < *R*_{C}, the unsupported nodule will keep on growing when larger than the critical radius (at least as long as the supply of nutrients and therefore the critical curvature *K*_{C} is maintained at the same level) but will disappear when it is smaller than *R*_{C}, meaning the equilibrium configuration is indeed unstable.

## 4. Growth of a tissue nodule perfectly wetting a disc

We now consider a nodule that perfectly wets the circular cap on top of a full cylinder for which the walls are not wettable (figure 1*c*). We suppose that there is a constant contribution of the surface stress at the contact between the support (light grey in figure 1*c*) and the nodule (perfect wetting). For this boundary condition, the nodule will satisfy mechanical equilibrium when it is a segment of a sphere or calotte with a radius *R* = 2/*K*. Figure 1*c* shows several examples for spherical nodules sitting on top of a full cylinder. Calling *H* the height of the calotte, measured from the surface of the support, the volume of the nodule is given as
4.1where *R*_{O} is the outer radius of the circular support. From now on, we express all linear dimensions in units of the outer radius, *R*_{O}, but keep to the same symbols for the sake of simplicity. Time *t* is also renormalized by introducing *τ* = *f γt*/

*R*

_{O}. Inserting this into equation (2.2) yields 4.2where the dot now indicates the derivation with respect to normalized time

*τ*. In addition to equation (4.2), we can also write a relation between volume

*V*and

*H*as well as

*H*and the surface curvature

*K*as 4.3which allows

*K*to be expressed as a function of

*V*. Insertion of this expression into equation (4.2) yields the rate of increase of volume as a function of the nodule volume. A similar approach also gives the rate of increase of the nodule height as a function of its current height, expressed analytically as 4.4Equations (4.2) and (4.4) now predict two equilibrium states for

*H*, when and , provided that

*K*

_{C}< 2. The stability of these states can be visualized with the help of figure 1

*a*,

*b*which shows plots of equations (4.2) and (4.4) for different values of

*K*

_{C}. Here, a short remark concerning the term ‘stable’ may be useful. We consider a kinetic process and a kinetic stability criterion. We observe a stable behaviour of the system, if a small increase of

*H*from a starting position leads to negative values of and consequently a tendency of the system to move back to the starting position—and vice/versa for a small decrease of

*H*and positive values of . Such stable configurations are highlighted by full circles in figure 1

*a*,

*b*. The situation is exactly opposite for the unstable equilibria (open circles). Stable equilibria correspond to nodules smaller than a half-sphere (full lines in figure 1

*c*), while nodules larger than a half-sphere are in unstable equilibria (broken line in figure 1

*c*). Hence, for a given biochemical driving force expressed by the critical curvature

*K*

_{C}, there are stable nodule configurations described by the corresponding curves in figure 1

*c*and with characteristic volumes and heights as given by figure 1

*a*,

*b*.

Figure 1*a*,*b* shows that growth of the nodule will not stop when the critical curvature *K*_{C} exceeds the value of 2. Interestingly (and in contrast to the unsupported nodule discussed in the previous section), the unstable equilibrium nodule (open circles in figure 1*a*,*b*) does not disappear, when its curvature is smaller than *K*_{C}. It runs for *K*_{C} < 2 along the curves in figure 1*a*,*b* from the unstable (open circles) to the stable state (full circles). In this sense, the nodule at equilibrium is dormant, but can reach a state of unconstrained growth as soon as the critical curvature *K*_{C} driving force increases for some reason beyond 2 (for example by a burst of growth factors). This is summarized in figure 1*d*, where *V*_{C}, the nodule volume corresponding to *K*_{C}, is plotted as a function of *K*_{C}. Portions of the curve with positive slope correspond to stable branches. Conversely, negative slopes correspond to unstable branches. Indeed, if a decrease in thermodynamic driving force (corresponding to a decrease in critical curvature) would lead to an increase in volume, the configuration cannot be stable.

## 5. Growth of a tissue nodule perfectly wetting an annulus

We now consider the growth of tissue pinned to the top of an open thick-walled cylinder (that is a flat circular ring, or annulus). We again normalize the linear space dimensions such that the outer radius of the open thick-walled cylinder is *R*_{O} = 1, and assume the walls of the cylinder are non-wetting. The dimension of the inner radius is kept as a parameter *R*_{I}.

The analysis is done in a completely analogous way to the previous sections. Outside the contact to the annulus (the top of the open thick-walled cylinder), the nodule will have a surface with constant curvature. Such surfaces are in general nodoids [39]. To describe this axisymmetric surface, we introduce a radial coordinate *r*, such that *R*_{I} ≤ *r* ≤ *R*_{O} = 1, and an axial coordinate *z* parallel to the cylinder axis of the support. The function linking *z* with *r* in the nodoid depends on two parameters, the (constant) surface curvature *K* and the position *R*_{H} (*R*_{I} < *R*_{H} < *R*_{O}), where there is a horizontal tangent to the profile (that is, the position of the maximum of the nodoid, sitting on top of the open thick-walled cylinder), and can be written as
5.1

with 5.2

The quantity *θ* represents the polar angle as parametrization with the starting value *θ*_{I} ≤ *θ*. The current formulation ensures that *z*(*θ*_{I}) = *θ*. The functions E and F are the incomplete elliptic integrals of the first and the second kind. Both signs in front of the bracket point to the fact that *z* may describe both the contour and its mirror image with respect to the plane *z* = 0. For details, see the electronic supplementary material and the references therein.

It is interesting to note that several topologically different shapes (as previously described by [15,17,40]) may fulfil the requirements of constant mean curvature and of meeting the boundary condition (pinning) on the top of the open thick-walled cylinder. These are sketched in figure 2. The first topology (figure 2*a*) is a spherical cap covering the central channel, the second topology (figure 2*b*,*c*) corresponds to various nodoids and the third (figure 2*d*) is a nodoid with a spherical droplet of identical curvature attached to it. This last (much simplified) configuration was considered here, because Lipowsky and co-workers [15,17,40] showed that nodoidal objects could bifurcate into conformations with a significant spherical bulge. Figure 2*e*–*g* was calculated using Brakke's ‘Surface Evolver’ software [37] (for details, see the electronic supplementary material).

Figure 3 shows relations between volume and curvature for all these cases and for three values of the inner radius (i.e. the ratios of the inner radius *R*_{I} and the outer radius *R*_{O}) of the support cylinder. For relatively thick-walled cylinders (figure 3*a*,*b*), there are three stable branches (labelled *S*), two of which correspond to nodoids and there are three unstable branches (labelled *U*). Branch *S*_{1} corresponds to a wetting of the ring-like cylinder top and the gradual development of a nodoid with increasing volume (at increasing critical curvature *K*_{C}, shapes 1, 2 and *G* in figure 2*b*). Beyond some value of the critical curvature *K*_{C} (close to 2.2 and to 4 in figure 3*a* and *b*, respectively), the branch ends and growth becomes unstable, see branch *U*_{1}. At the end of this branch, denoted by *R* in figure 3*c*, another stable branch *S*_{2} is reached, where a further increase of the volume is not possible. This is the largest possible nodoid that satisfies the boundary conditions. Nodoids along *S* will respond to a decrease in critical curvature by following path *S*_{2} (with shapes of the type shown in figure 2*c*). Interestingly, branch *S*_{2} is tangent (and nearly identical at small volumes) to another stable branch *S*_{3} that corresponds to a closed spherical cap (figure 3*a*) as a second possible topology. Hence, bifurcation from the topology shown in figure 2*c* to *a* is possible and very likely. When calculating the volume–curvature relations for the third topology (figure 2*d*), the corresponding branches (*U*_{2} and *U*_{3}) are barely visible in the upper right corner of figure 3*b* and completely outside the field of view in figure 3*a*. The situation is rather different for relatively thin-walled cylinders (figure 3*c*). It turns out that the curves for nodoids of type shown in figure 2*b*,*c*, with a drop-like spherical bulge (*U*_{2} and *U*_{3}) are actually in the same range as the ones for nodoids and bifurcations are possible. Most notably, the curves cross at a value of the critical curvature close to 8, where it is possible that the nodoid (of the type depicted in figure 2*c*) bifurcates to another topology (of the type depicted in figure 2*d*). This consists of a nodoid (of the type depicted in figure 2*b*) of smaller volume but with the same curvature by generating an additional spherical bulge with exactly this curvature taking up the extra volume. After this bifurcation, growth of the spherical bulge becomes unstable (branch *U*_{2} in figure 3*c*), in a manner equivalent to the unstable branch in figure 1*d*. The curves shown in figure 3 were also verified with numerical calculations using ‘Surface Evolver’ [37] that also showed a bulging bifurcation (figure 2*g*). Despite the strong simplification of the analytical model (as shown in figure 2*d*), there is excellent qualitative agreement between the analytical and numerical solutions.

This analysis of potential bifurcations merits a more thorough analysis. For this, we plot the total free surface area of the tissue nodule as a function of critical curvature for the relatively thin-walled cylinder with an inner radius of 0.75 (figure 4). Given that surface energy governs the development of shape, the total area, which is proportional to total surface energy, will tend to be minimal. A bifurcation from one topology to another is only possible without a jump in volume and at given curvature. Indeed, curvature is imposed on the system by the biochemical driving force expressed by *K*_{C} and changes in volume require the synthesis of additional tissue and cannot, therefore, be instantaneous. In figure 4, two possible bifurcations can be observed, namely one between the paths *U*_{1} and *U*_{2} crossing close to a critical curvature of *K*_{C} = 8 and the other one between the paths *S*_{2} and *S*_{3} approaching asymptotically for small curvatures. These two points are marked by arrows in figure 4. It turns out that configurations due to *S*_{2} and *S*_{3} have nearly the same surface area at small curvatures, which means that a bifurcation is likely but not necessary. However, if biochemical driving forces increase the critical curvature *K*_{C} again, it is more likely that the system follows the path *S*_{2}, which has the lower surface area and consequently surface energy at equal critical curvature. The bifurcation between the paths *U*_{1} and *U*_{2} is of a totally different nature. When curvature increases along the path *S*_{1} and exceeds the maximum value at *G*, so that further growth becomes unstable along path *U*_{1}, the formation of a spherical bulge significantly lowers the surface energy, so that a bifurcation from path *U*_{1} to path *U*_{2} becomes highly probable. The bulge itself is unstable, and will continue to grow in a similar way to the unstable branch of the calotte (figure 1*d*). It should be mentioned again that for thicker cylinder walls, *U*_{1} and *U*_{2} do not intersect (figure 3*a*,*b*) so that the bifurcation towards uncontrolled growth is not expected above a certain wall thickness. Without the bifurcation towards a spherical bulge (which is the case for inner radii significantly smaller than 75% of the outer radius), the configurations labelled *G* and *R* in figure 3 correspond to those with the largest curvature and the largest volume, respectively. Quantitative data for these two configurations are given in table 1.

## 6. Comparison to callus formation and bone healing

The model described above has only one free parameter, the critical curvature *K*_{C}, and is based on three main biological hypotheses. To recap, the model relies on the following conditions:

(1) Tissue is remodelled on a timescale that is much shorter than growth, which justifies the description of the tissue as an effective ideal fluid (on the timescale of months) [10].

(2) There is a surface stress, which results in a surface curvature-dependent internal pressure inside the tissue [9,11,12,21,41].

(3) Growth of a tissue, by cell proliferation and extracellular matrix formation, depends on its pressure through the growth equation. Such a dependency on pressure is known from the growth of tumour spheroids [18,42] and, in principle, can be controlled by the organism (e.g. via growth factors).

Despite the simplicity of the above model, the coupling of growth with the external surface geometry via the internal pressure leads to surprisingly complex behaviour as a function of the boundary conditions. Interestingly, this behaviour is reminiscent of many features seen in biological tissues. While an unsupported nodule is unstable (disappearing at small driving forces and showing unlimited growth at large driving forces), there are equilibrium sizes for the nodule when adhering to a substrate. The stabilizing effect of substrates on growing tissues has been observed in cancer metastases [18,43] and clearly plays an important role in biofilm formation [44].

Our model shows that growth of a tissue is also dependent on the shape of the surface upon which it adheres. When the tissue grows on top of a hollow cylinder, which is the same geometry as the free ends of bone in the simulated fracture gap after an osteotomy [45], there is a maximum height the tissue will grow to, for a given value of critical curvature. Assuming that similar growth occurs from both sides of the fracture gap in an osteotomy, our model thus predicts the existence of a maximum bridgeable gap size. This is reminiscent of the so-called ‘critical defect size’ known from bone healing [46]. Such calculations of course assume no wetting of the external bone surface which may occur during bone fracture healing due to damage or removal of the periosteum (e.g. [47]). The prediction of the model is that the critical gap size is 70%, 78% and 61% of the outer diameter of the cylindrical tube, when the inner diameter is 75%, 50% and 25% of the outer diameter, respectively (table 1). Most land mammals have a ratio of inner to outer diameter lower than 75% [48], and typically have critical defect sizes of one to two times the outer bone diameter [49]. For example, the ratios of critical defect sizes to bone diameters are reported to be: 1.3–1.5 times in feline tibias [49], 1.5–2 times in dog ulnas [49] and 2–2.5 times in sheep tibias [46]. These ratios are of a similar order of magnitude to the model predictions. For thin-walled cylinders (wall thickness less than 25% of the radius), the model predicts a bifurcation towards uncontrolled growth of spherical bulges, which could have implications in the fracture healing of thin-walled bones such as those found in birds. Our theory predicts that if the fracture gap is too large and then tissue starts to resorb, a cap will form that will gradually close the tubes at each side of the gap (figure 2*e*). Remarkably, such closure is observed in pseudarthrosis or non-union of bone fractures (e.g. [50]). In our theory, this arises due to a kinetic instability of the growth process, which is described by the evolution equation of the shape of the actual configuration leading to different paths for resorption and growth (figure 3).

In conclusion, if the above conditions 1–3 are fulfilled in a growing tissue, then the organism needs to control only one parameter, namely the proliferation rate (that leads to the growth pressure or the critical curvature as used in our model), in order to achieve a complex sequence of structures surprisingly similar to *in vivo* observations. In addition to bone fracture healing, our model may be interesting for biofilm formation and cancer research as it shows how the presence of tissue surface tension, and liquid like properties can lead to radically different behaviours of tissues on surfaces.

## Authors' Contributions

P.F. and F.D.F. designed the study and developed the analytical model, G.Z. solved the nodoid equations and J.D. performed the numerical growth simulations. P.F., F.D.F. and J.D. wrote the manuscript and all authors gave final approval for publication.

## Competing Interests

We have no competing interests.

## Funding

We acknowledge the Alexander von Humboldt foundation for supporting the visit of F.D.F to the Max Planck Institute of Colloids and Interfaces through the Humboldt Award.

## Acknowledgements

The authors wish to thank Ron Shahar and Reinhard Lipowsky for stimulating discussions.

- Received February 9, 2015.
- Accepted May 1, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.