The Emergence of Complexity from a Simple Model for Tissue Growth

The growth of living tissue is known to be modulated by mechanical as well as biochemical signals. We study a simple numerical model where the tissue growth rate depends on a chemical potential describing biochemical and mechanical driving forces in the material. In addition, the growing tissue is able to adhere to a three-dimensional surface and is subjected to surface tension where not adhering. We first show that this model belongs to a wider class of models describing particle growth during phase separation. We then analyse the predicted tissue shapes growing on a solid support corresponding to a cut hollow cylinder, which could be imagined as an idealized description of a broken long bone. We demonstrate the appearance of complex shapes described by Delauney surfaces and reminiscent of the shapes of callus appearing during bone healing. This complexity of shapes arises despite the extreme simplicity of the growth model, as a consequence of the three-dimensional boundary conditions imposed by the solid support.

properties of single cells and that can be described by physical laws [7,8]. In addition, tissues are active and can grow in response to physical, chemical and biological stimuli from the environment. When this growth is slow compared to local remodelling rates of the extracellular tissue, shear stresses in the matrix are expected to relax fast enough to reach the mechanical equilibrium status of a fluid before the volume significantly increases. Under such conditions, the shape of the tissue "droplet" corresponds at any moment in time to the equilibrium shape governed by the surface stress and corresponding energy state [9]. Despite the extreme simplicity of such a model, fluid drops with remarkably complex shapes can emerge simply by interacting of the surface with three-dimensional substrates of complex geometry, see e.g. [9][10][11].
In order to get better insight into the development of geometric complexity, we further investigate a simple model introduced in earlier work [9], where growth is described as the time derivative of the tissue volume, V , being proportional to − K − K c V, with K being the total surface curvature (i. e. the trace of the surface curvature tensor) of the tissue aggregate. K is twice the mean curvature and should not be confused with the Gaussian curvature which plays no role in the following (we use the symbol K to be consistent with our previous papers). The parameters and K c correspond, respectively, to the specific surface energy and a critical curvature related to the driving force for growth. This model was derived based on earlier work [8,[12][13][14], whereby tissue is described as a fluid. The surface stress state causes a surface-curvature dependent pressure inside the tissue. This internal pressure is proportional to the curvature K according to the Laplace-Young law. It counteracts the pressure induced by tissue growth and compensates it exactly at a curvature K c .
When writing this relationship for a spherical droplet of radius R, it becomes apparent that the tissue growth model belongs to a wider class of models that has been intensively studied in the context of phase separation [15][16][17]. Indeed, the equation Ṙ = A n R∕R c − 1 ∕R n describes the growth of a particle from a medium that provides the necessary source of matter, whereby R c is the critical radius and A n a kinetic coefficient. The integer n is = 1 when particle growth rate is limited by the accretion reaction on the surface of the particle and = 2, 3, 4 when its growth is limited by diffusion through the bulk of the medium (such as a solution or an alloy), diffusion along planar structures (such as grain boundaries) or along linear paths (such as dislocation lines), respectively [18]. Interestingly, the tissue growth model considered here corresponds to the same class of models, with n = 0. All these models have in common that a particle will grow when its radius R is larger than R c , but shrinks when it is smaller. In phase separation models, R c depends on the chemical potential difference between the phases and, thus, on the supersaturation of the solution from which the particle grows [19]. In our biological tissue growth model, the meaning is analogous and R c describes the driving force for tissue growth [9]. In general, R c may depend on time through the interaction with other growing or shrinking particles. The kinetics of phase separation processes and the associated growth and coarsening of particles has been analysed in great detail both by statistical physics and lattice methods [20][21][22][23] as well as through particle models [24][25][26][27][28].
A significant increase in structural complexity occurs when phase separation models are combined with wetting of a substrate to which the growing new phase droplet adheres. Surface-directed phase separation has been studied in detail since the 1990s [29][30][31][32]. The wetting surface then represents a boundary which confines the growth in certain directions, leading to a deviation from spherical droplet shapes. In a previous paper we investigated this for the tissue growth model assuming that the tissue wets planar annulus-like substrates [9], an idealized geometry motivated by healing of the fracture gap in long bones [33,34]. Here we expand this approach, exploring the increased structural complexity based on wetting of non-planar surfaces. Inspired again by callus formation during bone healing [33,34] we particularly chose geometries where the tissue is wetting parts of the outer shell of a cylinder.
Two stability criteria need to be considered when analysing droplet shapes. The first one is mechanical stability for a constant droplet volume that we assume to be reached nearly instantaneously. Mechanical stability corresponds to a minimum of the total (Gibbs) energy, accounting for the surface energy and the internal pressure that fulfils the Laplace-Young equation. In addition, there is a second stability criterion with respect to unlimited growth or resorption of the tissue droplet. Kinetic stability implies that the equation of growth has a stable solution with constant volume ( V = 0, which means that K = K c ). Given that V ∼ − K − K c V, such a solution is stable if a fluctuation corresponding to a small increase in volume leads to an increase in curvature, so that V < 0 by the kinetic equation [9]. This is not the case for spherical droplets, for example, because the curvature of a sphere decreases when its volume increases. This is the well-known instability of droplets at the critical radius. Larger particles have a tendency to grow, while smaller ones shrink. We showed in our previous work that this situation changes for droplets that partially wet a substrate, so that kinetic stability exists under certain conditions [9].
In this paper we explore the mechanical and kinetic stability of tissues growing on nonplanar surfaces inspired by simple "bone-like" geometries. The analyses are not meant to accurately describe the growth of any biological material that has usually hierarchical structure and is controlled by a variety of signals [35]. Rather we are exploring the predictions of the model when the tissue droplets adhere to non-planar substrates. A surprisingly big manifold of shapes of bodies are shown to be possible, some of which are kinetically stable. These observations may have implications for our general understanding of the development of biological form, as it shows that surprising complexity can arise from extremely simple physical principles.

A Simple Growth Model
What follows is a short outline of our growth model based on a continuum mechanics framework. For a description of growth using continuum mechanics we refer to the overview article by Ambrosi et al. [7] and the book by Goriely [36]. Our approach follows the concept of surface growth as discussed by Skalak, see e.g. [37]. However, growth is not considered as the time-evolution of a finite displacement field. The model description includes continuum thermodynamics as introduced by Ambrosi and co-workers [8,12] and thermodynamic forces. Such an approach is followed by several groups, see [38] for a recent example. Our model was introduced to study tissue growth in confined geometries [13,14] involving growing volumes together with moving surfaces. Recently Swain and Gupta [39] extended the thermodynamic concept including an incoherent interface. Details of our model can be found in Sect. 2 of [9], only the key aspects are outlined here: The growing tissue with volume V is assumed to behave like an ideal fluid forming a surface A with constant curvature K (trace of the curvature tensor, e.g., 2/R for a sphere, 1/R for a cylinder, with radius R, respectively), a surface energy and a surface stress s .
According to the Laplace-Young equation (from 1805, for details see [40]) the internal pressure p follows from the local mechanical equilibrium of a surface element as [40] In the following context we follow the standard assumption that the surface energy γ exhibits the same value as the (mechanical) surface stress. The justification of this assumption goes back to seminal works from the seventies of the last century [41,42].
Any possible singularities, incompatibilities and residual stresses are assumed as compensated by relaxation. Following the thermodynamic procedure by Ambrosi and co-workers, see above and applying Eq. (1), a growth law for the volume V of the growing tissue is derived in [9] as with V the time rate of V, f a kinetic constant and K c a critical curvature, which is defined as the ratio of a chemical potential μ and the surface stress s . With = s and a linear transformation of time, = f t, we obtain The presence of K c meets the recently shown necessity to include a critical stress in a growth concept (note that K c is proportional to a critical pressure p c according to Eq. (1) to obtain, e.g., the homeostasis [43,44]. The growth law (Eqs. 2, 3) has been successfully applied to simulate rather intricate topologies of grown tissue structures that would develop during bone healing on a flat circular ring acting as a substrate [9]. One limitation of such a growth law is that it assumes growth is uniform throughout a tissue. It is clear that cells deep inside a growing tissue may change their behavior (see e.g. [45].) due to hypoxia, differentiation or even formation of extracellular matrix and further work will be required to incorporate the process of tissue maturation into such a model. Cell-based models such as (e.g. [46,47].) also show similar behaviors in cases when nutrients are sufficiently supplied and no cell-maturation occurs. Our surface-stress based approach does not take into account the relative roles of cell-cell adhesion or cortical tension [6] which may require (local) elastic contributions to the surface energy. In the following we explore solely how such a growth model functions for a constant surface energy.

Stability of Equilibrium and of Growth Process
The term "stability" is interpreted differently in biology. For example, Finean wrote on "chemical stability" of a membrane already more than five decades ago [48]. Islam and Dimov reported on the "mechanical stability" of membranes nearly three decades ago [49]. The term "mechanical stability" has often been used for high stiffness or high deformation resistance of membrane-like structures in colloquial language. We consider here "stability" in the sense of both the classical thermodynamics and the kinetics of continua.

Continuum Description
We consider a tissue body with the volume V and the surface A = A S + A W . The wetted surface A W is considered to be spatially fixed to a solid substrate. The tissue surface A S is described by a constant trace of the curvature tensor K, owing to the fact that a free surface of a fluid has a constant mean curvature at mechanical equilibrium.

Mechanical Stability
Our, starting point is the total Gibbs energy G = G(K) of the system with the only free variable K given as see, e.g., Sect. 2.2 in [50]. Note that the Gibbs energy must also include the external load term, which has often been ignored.
Equilibrium exists, if the first variation of G becomes δG = 0, yielding The equilibrium state attains (thermodynamic or static) stability, if G obtains a minimum, which means that the second variation δ 2 G becomes > 0, yielding However, this second variation δ 2 G must include Eq. (5), as a necessary condition for an extremum, together with Eq. (1) and γ = γ s as side conditions. Therefore, the derivative A s ∕ K, see Eq. (5), is expressed by using these side conditions as A s ∕ K = K V∕ K. Inserting this relation in Eq. (6) yields concluding that the equilibrium state is stable, if V∕ K > 0.
The problem can also be formulated as a variational problem with the Lagrangian describing the tissue through looking for a minimum of A s for a fixed volume V = V 0 with p taking the role of a Lagrange multiplier, see, e.g., Sect. 2.3 in [50]. However this variational formulation deals with geometrical quantities and not mechanical quantities.
It shall be mentioned that Gillette and Dyson [51] and later Lowry and Steen [52] investigated surfaces of revolution within a simplified variational framework. Both groups completely describe all surfaces of revolution with a constant trace of the curvature tensor K. Nodoids, undoloids and catenoids (in addition to spheres and cylinders) meet the condition, Eq. (8) (see Appendix). The variational framework as mentioned above can be extended to general surfaces. However, in most cases, only a computational algorithm allows finding constant mean curvature surfaces; here we refer to the Surface Evolver developed by Brakke, see [53][54][55].
A possible extension of our approach would be to take both stretching and bending energy terms of the surface energy into account (as thin shells instead of membranes) in G (in our case G(K)), see the early work by Helfrich [56], the book by Safran [50] and the recent paper [57]. Also a continuum mechanical framework has been provided for this energy minimization problem [58,59]. By including bending and stretching deformation terms in the total elastic strain energy one enters the field of structural stability. This very prominent field of the (structural) stability of equilibrium is being applied to biological systems as, e.g., Sect. 15.3 in [36], and recent physically oriented publications on the stability of shells [60], sometimes in the context of bioinspired materials [61].

Kinetic Stability
Concerning the growth law, Eq. (2), we previously formulated a stability criterion [9] which can be summarized as following: (i) stable behaviour of the system is observed, if a small increase of V from a starting position leads to negative values of V , and consequently a tendency of the system to move back to the starting position; (ii) unstable behaviour of the system exists, if a small increase of V from a starting position leads to positive values of V , and drives, therefore, the system further away from the starting position.
This criterion is based on Lyapunov's stability theory framework, see e.g. [62]. The system is subjected to a small perturbation of the (static) equilibrium configuration. If the time evolution of the perturbation is decaying to zero, then process stability occurs. Recently this criterion has been denoted as "mechanobiological stability" criterion, see, e.g., the papers by Cyron et al. [63][64][65]. This "mechanobiological" aspect is extended to a system approaching to a "homeostatic" state, which is attained by remodelling, for details see Cyron and Aydin [66].

Conclusions About Stability
In Sect. 4 of our previous paper [9] the stability criterion from above was reformulated by applying the functional relation between the volume V and the trace of the curvature tensor K as follows: • for branches in the K-V-diagram with positive slope the growth process is stable; • for branches in the K-V-diagram with negative slope the growth is unstable.
Considering the stability of equilibrium (Sect. 2.2) and the stability of the growth process (Sect. 2.3), leads to the conclusion that a single criterion ensures both static and kinetic stability of the growing object.

Admissible Surfaces with Constant Curvature K
We investigate axisymmetric tissues described as fluid droplets that completely wet the top and part of the sides of a hollow cylinder (Fig. 1a, b). The wettable surface geometry consists of the outer wall (with a wetting length of WL) and the top annulus of the cylinder.
To be mechanically stable (i.e. to satisfy the Laplace-Young equation, Eq. (1)) the surfaces must have a constant trace of the curvature tensor, denoted as the curvature K. In general, the set of axisymmetric constant mean curvature surfaces are the Delaunay surfaces consisting of nodoids, unduloids, catenoids, cylinders and spheres [51,67]. Due to the sharp interface on the outer radius of the annular top of the wettable surface, there are eight different axisymmetric surface configurations which satisfy the boundary conditions (Fig. 1c). It is assumed that tissues will only be pinned to the top outer edge of the cylinder, if the resultant membrane force on this edge is directed inwards the substrate. The first three configurations ( Fig. 1c-1a, 2a and 3a) are the same surface configuration as reported in our previous paper [9] together with a monolayer on the outer substrate surface. The first two surface configurations (1a and 2a) are nodoids, the third consists of two joined spherical caps with identical curvature, one of them convex up attached to the outer edge of the top annulus, the other convex down attached to the inner edge of the annulus. The next three configurations ( Fig. 1c-1b, 2b, 3b) consist of surfaces in which the outer substrate surface is completely wetted by a fluid, however pinned to the outer edge of the annulus. The axisymmetric solutions of tissues pinned to the outer part of the cylinder are, in order of increasing volume, described by cylinders (for zero volume), unduloids, spheres, and then nodoids.
The assumption is here that the growing tissue/fluid has the same pressure (curvature) in both regions however remains pinned at the outer edge of the annulus resulting in a kink in the surface along this line. If the surface becomes depinned on this line, additional two axisymmetric surface configurations can be found: nodoid surfaces (Fig. 1c, configuration 4) and spherical surfaces (Fig. 1c, configuration 5). Surface configurations that intersect with the substrate are not considered (Fig. 1c i); neither are those in which the surface curvature of the tissue on the top is different from that of the tissue on the sides (Fig. 1c ii). Any other (similar) combination of constant mean curvature surfaces which introduce a kink (that points inwards) have no reason to be pinned at the ridge (Fig. 1c iii) and are thus excluded from our analysis.

Results and Discussion
Using the equations of Gillette and Dyson ( [51], see Appendix) the volumes of all axisymmetric configurations that satisfy the boundary condition were calculated and plotted as a function of curvature (Fig. 2). A large number of solutions is found with configurations providing volume-curvature relationships having regions with both positive and negative slopes as well as having different branches with different volumes for the same curvature.
To better understand the differences between the individual configurations it helps to follow the three trajectories in K-V space indicated by the black/grey, blue and red lines in Fig. 2 starting at zero volume. Let us consider the first trajectory given by the black (1a), and grey (1b) curves. This trajectory starts with a monolayer of tissue growing on the top and side surfaces of our cylinder. For the curve (1a) only the tissue on the upper surface will grow with solutions being nodoid surfaces pinned to the upper annulus [9]. Tissue will then only start to grow on the outer sides of the substrate when the upper curvature of the tissue matches the curvature of the cylindrical outer walls (K = 1). From this point on the tissue will grow on both the sides and the top (with identical curvatures), with nodoid solutions on the top annulus and unduloid solutions on the sides. At a curvature of 1.7 the outer side surface of the tissue transits from an unduloidal through a spherical to a nodoidal surface. The two nodoid solutions, despite having the same curvature, do not have a common tangent, and are still pinned to the corner. This arises as there are two nodoid solutions for surfaces attached to the upper annulus, and the one with smaller volume is described in this case. The system comes to a point of maximum curvature (point D in Fig. 2) after further growth. At this point the slope of the curve becomes negative and the system becomes unstable to further growth.
The second trajectory (Fig. 2 in blue 2a, 2b and 4) and the third trajectory (Fig. 2 in red 3a, 3b and 5), both start with a thin tissue layer covering the hollow part of the inner cylinder like a lid. If this tissue layer breaks down or is pierced, then nodoid solutions (from the higher-volume branch) are stable. These would then grow following line 2a, again until a curvature of 1 is reached, at which point growth is triggered also on the outer cylinder going from unduloid to sphere to nodoid solutions. When the curvature reaches a value of around 2 (point B in Fig. 2), both the upper nodoidal surface and the outer nodoidal surface are tangent and the surface can unpin. At this point growth becomes unstable and will follow the branch 4. The resulting configurations have a small region of stability with a maximum volume at point C (Fig. 2). The restriction to axisymmetric solutions adds additional stability to the configurations. We explored this in more detail by performing numerical simulations of these surfaces using Brakke's Surface Evolver [53,54] and demonstrated that small perturbations in the surface geometry around point C (Fig. 2) lead to asymmetric bulging similar to what was observed in our previous paper [9].
For the third trajectory ( Fig. 2 in red 3a, 3b and 5) the thin tissue layer starts growing on the top as two connected spherical caps (see configuration 3a in Fig. 1). These tissues again grow up to a curvature of 1, at which point tissue starts growing on the outer side surface. Growth continues on the top and the sides along trajectory 3b up to point A, which is a spherical surface just touching the outer edge of the cylinder. At this point the surface is unstable and will then follow branch 5 (red-dashed line in Fig. 2).
The set of these growth trajectories suggests that only certain configurations are attainable by growth from other configurations. More information about the relative stability of the different configurations can be seen in Fig. 3, which illustrates the relationship between area and total Gibbs energy ∫ tissue body ( dA s − pdV) vs curvature and volume, see also Sect. 2.2. of this manuscript. The plot of area vs curvature (Fig. 3a) is very similar in appearance shape to Fig. 2. This similarity can be understood from Fig. 3B, since the area vs volume curves have similar trajectories. Configurations 1a and 1b have the lowest area until they cross the configurations for the spherical cap (5 in Fig. 1c), which has the lowest surface per volume at large volumes. Although the area vs volume curve (Fig. 3b) may give an indication of which systems are likely to be more stable (in terms of minimum surface area), it is more informative to plot the total energy of the system, i.e. the sum of the surface energy and the potential energy of the load (i.e. the pressure p). One can see clearly in Fig. 3c, d that a growth trajectory following the black (1a and 1b), blue (2a and 2b) or red (3a and 3b) curves all result in a decrease in energy up until the point of instability highlighted by the black dots in all figures. The lower trajectory in black and grey has a much lower energy than the other two branches. A consequence of this is that transitions between the branches 1a or 1b to the other branches (2a, 2b, 3a 3c etc.) are unlikely as they require a large spontaneous change in the energy of the system. In contrast, the small differences between the branches 2 and 3, not only in terms of geometry, but also in terms of the total energy, makes transitions between the different configurations on these two branches more likely upon perturbation. Once the system becomes unstable with respect to growth (black points in Fig. 3) the total energy of the growing tissue increases. These plots again highlight the complexity of potential configurations arising from relatively simple boundary conditions and growth laws.

Conclusions
A phenomenological model of tissue growth was studied to highlight the importance of the geometry of the substrate on the resulting form of the tissue. First of all the substrate is essential for the existence of non-trivial equilibrium forms since without substrate the tissue would either shrink or grow infinitely [9]. Secondly, the substrate geometry gives rise to a multitude of tissue forms characterized by a constant mean curvature (Fig. 1). Thirdly, these equilibrium forms are not equivalent in terms of their stability and possible transitions between forms can be predicted. The consideration of a wettable substrate that is non-planar results in new features of the tissue shape. As shown in the studied example the surface can get pinned at edges of the substrate leading to a kink in the surface. This pinning can be stable (due to the restriction that the interface cannot penetrate the substrate like in the configuration i) in Fig. 1C) or unstable giving rise to a change in configurations. The multitude of resulting tissue shapes is remarkable having in mind the simplicity of the growth law (Eq. 3) including only one crucial parameter K c defined as the ratio between the chemical potential and the surface stress, see [9]. While a large chemical potential favors growth, a large surface tension retains growth. Consequently, K c can be described as the propensity of the tissue to grow. Biochemically this propensity can be enhanced by either increasing the chemical potential by growth factors, or by decreasing the surface tension, e.g., by a slackening of the cytoskeleton of the cells in the tissue [45].
The important influence of geometrical constraints on the evolution and the equilibrium configurations is well accepted in the context of physical systems. The role of geometrical constraints to influence or even to guide biological processes is also starting to emerge [68][69][70][71][72]. In a biological context the role of geometrical constraints goes beyond what was explored for physical systems. Indeed, a particular fascination of biological systems is their display of intricate forms. The seminal book by D'Arcy Thompson "On Growth and Form" [73] is after its centennial anniversary still a source of inspiration for the richness of form in biology. The model investigated here is just a just another example for the emergence of complex three-dimensions form from a simple growth law and, therefore, pertains to the more general field of pattern formation in biology [74].
where F( , t) and E( , t) are respectively the first and second kinds of elliptic integrals given by and The parameter shifts the curve up and down parallel to the z-axis, and the parameters and are scale and shape parameters respectively.
Using the first and second elliptic integrals, it is possible to calculate the mean curvature of a surface of revolution defined by these equations which is The Gauss curvature For the nodoid and unduloid solutions on the outer wall, the volume is given by and the area is given by where t lim is the value of the parametrization at the top edge. The volumes of the nodoid portions (attached to the top) are calculated in a piecewise manner, by first calculating the volume of the outer part of the nodoid to the parameter value of the maximum axial height and second by subtracting from this the volume of the portion that curves inwards. Area is calculated in a similar manner.  cos ) .