Are Homeostatic States Stable? Dynamical Stability in Morphoelasticity

Biological growth is often driven by mechanical cues, such as changes in external pressure or tensile loading. Moreover, it is well known that many living tissues actively maintain a preferred level of mechanical internal stress, called the mechanical homeostasis. The tissue-level feedback mechanism by which changes in the local mechanical stresses affect growth is called a growth law within the theory of morphoelasticity, a theory for understanding the coupling between mechanics and geometry in growing and evolving biological materials. This coupling between growth and mechanics occurs naturally in macroscopic tubular structures, which are common in biology (e.g., arteries, plant stems, airways). We study a continuous tubular system with spatially heterogeneous residual stress via a novel discretization approach which allows us to obtain precise results about the stability of equilibrium states of the homeostasis-driven growing dynamical system. This method allows us to show explicitly that the stability of the homeostatic state depends nontrivially on the anisotropy of the growth response. The key role of anisotropy may provide a foundation for experimental testing of homeostasis-driven growth laws.


Introduction
Biological tissues exhibit a wide range of mechanical properties and active behavior. A striking example is biological growth in response to the tissues mechanical environment. Artery walls thicken in response to increased pressure (Berry 1974;, axons can be grown by applying tension (Lamoureux et al. 1992;Recho et al. 2016), and plant growth is driven by various mechanical cues (Goriely et al. 2008;Boudaoud 2010). The general idea underlying these phenomena is that the internal stress state is a stimulus for growth. As stress is rarely uniform, mechanically induced growth often coincides with differential growth, in which mass increase occurs non-uniformly or in an anisotropic fashion. In turn, differential growth produces residual stress, an internal stress that remains when all external loads are removed, appearing due to geometric incompatibility induced by the differential growth. Residual stress has been observed in a number of physiological tissues, such as the brain (Budday et al. 2014), the developing embryo (Beloussov and Grabovsky 2006), arteries (Fung and Liu 1989), blood vessels (Fung 1991), solid tumors (MacArthur and Please 2004), and in a wealth of examples from the plant kingdom (Goriely 2017). In many cases, residual stress has been found to serve a clear mechanical function; for instance in regulating size and mechanical properties.
Many living tissues actively grow in order to maintain a preferred level of internal residual stress, termed mechanical homeostasis. This phenomenon is characterized by growth being induced by any difference between the current stress in the tissue and the preferred homeostatic stress. Mechanically driven growth toward homeostasis poses several interesting and important questions, at the biological, mechanical, and mathematical level. For instance, what determines the homeostatic stress state? At the cellular level, the growth response may be genetically encoded, with a homeostatic state manifest by differential cellular response to mechanical stimuli. From a continuum mechanics point of view, a residually stressed configuration is typically thought of as corresponding to a deformation from an unstressed configuration; however, it is not clear that such a deformation should exist to define a homeostatic state. Connected to this is a question of compatibility: is it actually possible for a system to reach mechanical homeostasis? For example, the boundary of an unconstrained tissue will always be traction free, and thus, if the homeostatic stress for those boundary cells is nonzero, then the system can never completely reach homeostasis. From a dynamics point of view, there is a natural question of stability: is the homeostatic state stable, i.e., if the system is perturbed from its homeostatic equilibrium, is it able to grow in such a way to return to this state? There is also a practical issue of connecting experiment to theory: how does one quantify the homeostatic state and form of growth response?
Mathematical modeling can be of significant value in addressing such questions and in suggesting potential experimental measures to quantify the properties of homeostasis. In the simplest and most widely used form, the mathematical description involves a growth law of the form Here overdot represents time derivative, G is a growth tensor, characterizing the increase or decrease in mass as a local property, T is the Cauchy stress tensor, T * is the preferred homeostatic stress tensor, and K is a fourth-order tensor characterizing the growth response rate due to differences in current and preferred stress. Laws of the form (1), or slight variations thereof, in which growth is coupled to Cauchy stress, have been examined by a number of authors (Vandiver and Goriely 2009;Bowden et al. 2016;Ramasubramanian 2008;Taber 2008), though the most appropriate form of growth law is a much-debated issue (Taber 1995;Ambrosi et al. 2011;Jones and Chapman 2012;Goriely 2017). An alternative but related approach involves coupling growth and Eshelby stress (Ambrosi and Guana 2005) based on thermodynamical arguments (Epstein 2000;Erlich et al. 2015). Attempts to restrict the form of growth laws through thermodynamical considerations such as the Coleman-Noll procedure (Coleman and Noll 1963) have been of limited success due to the inherent thermodynamical openness and non-equilibrium nature of biological systems (Maugin 1999;Lebon et al. 2008). The integration of micro-mechanical models with tissue-level modeling has also been difficult, partly because the lack of periodicity and crystal symmetry in biological tissues makes the application of homogenization techniques difficult (Chenchiah and Shipman 2014). Growth dynamics that depend on the current stress state are inherently challenging to study analytically. Both stress and growth will tend to be spatially dependent, with stress being determined through the solution of a force balance boundary-value problem, and thus, any model will by nature involve a partial differential equation system. The situation is simplified somewhat by the slow-growth assumption, which states that growth occurs on a much longer timescale than the elastic timescale and hence the system is always in a quasi-static mechanical equilibrium.
In this paper, we study mechanically driven growth in the context of growing tubular structures. One motivation for a cylindrical geometry is that such structures are ubiquitous in the biological world, from plant stems ) to axons and airways (Moulton and Goriely 2011a, b), and exhibit diverse mechanical behavior. Working within a constrained geometry will also enable us to gain qualitative insight into the dynamics of structures with growth driven by mechanical homeostasis and to formulate a basic framework for studying the stability of a homeostatic state. Even in an idealized geometry, the full growth dynamics still consists of a set of partial differential equations, with mechanical equilibrium requiring the solution of a boundary-value problem at each time step, and a highly nonlinear growth evolution for components of the growth tensor. There is no mathematical theory, yet, that allows for such an analysis. Our approach is therefore to devise a discretization through a spatial averaging scheme that converts the system to a much more manageable initial-value problem, to which we can apply standard techniques from dynamical systems. The discretization we propose consists of defining annular layers of the tubular structure, such that growth is uniform in each layer, driven by averaged values of the stress components in a law of the form (1). While this approach enables us to study efficiently properties of the continuous (non-discretized) system as the number of layers increases, for a smaller number of layers it is also a useful model of a multilayered tube commonly found in many biological systems. This paper is structured as follows. In Sect. 2, we discuss the general deformation and growth dynamics for a tubular structure that is homogeneous in the axial direction. In Sect. 3, we focus on a tubular system made of two layers, illustrating the main ideas of our discretization approach and illustrating the rich dynamics of this system. In Sect. 4, we generalize from two to N layers. Here we find a rapid convergence of behavior as the number of layers increases and investigate how the anisotropy of the growth affects the stability.

Kinematics
We consider a cylindrical tube, consisting of an incompressible isotropic hyperelastic material, the inner wall of which is attached to a fixed solid nucleus, with the outer wall unconstrained (see Fig. 1). We restrict to growth and deformations only in the cross section, such that the cylindrical geometry is always maintained and there is no axial strain. Moreover, we assume that there are no external forces, so that any deformation is caused purely by growth and the elastic response.
Geometrically, we work in a planar polar coordinate basis e R , e θ (the same basis vectors apply to both initial and current configurations), in which the deformation can be described by the map x : B 0 → B t given by: For this map, the deformation gradient is The elastic deformation gradient takes the form Incompressibility requires det A = 1; we thus define α := α θ , so that α −1 = α r . We assume a diagonal growth tensor where the difference between radial growth (γ R > 1) and circumferential growth (γ θ > 1) is shown schematically in Fig. 2. In matrix form (with the basis e R , e θ implied), we have In the initial (stress-free) reference configuration B 0 , the inner cylinder wall is located at R 0 = A 0 and the outer wall is located at R 0 = B 0 . From the morphoelastic decomposition F = AG, we find r = γ R /α and r /R 0 = αγ θ . By eliminating α, we obtain Imposing the boundary condition r (A 0 ) = A 0 , due to the unmoving solid nucleus, we can integrate (7) as The mathematical and biological significance of our solid nucleus setup is discussed in Sect. 4.5.

Mechanics
Given that all deformations are diagonal in the coordinate basis considered here, the Cauchy stress is also diagonal Let W α R , α θ be the strain-energy density, which relates to the Cauchy stress tensor by T = AW A − p1, where p is the Lagrange multiplier enforcing incompressibility.
In components, this reads With no external loads, mechanical equilibrium requires div T = 0, which takes the form Defining W (α) := W α −1 , α , we have The above formulation is valid for all isotropic, incompressible material. To make progress in our analysis, we further restrict it to the simplest possible material model of a neo-Hookean quadratic strain-energy density given by: where μ in small deformations can be identified with the shear modulus of the material. In this case, for which (11) becomes Along with (15), we impose T R R (B 0 ) = 0, i.e., the outer edge is stress-free. Equations (7) and (15), along with boundary condition T R R (B 0 ) = 0, completely determine the deformation and stress state. Due to the fixed inner boundary condition, for a given growth tensor (7) can be integrated separately, i.e., the deformation is determined independently from the stress, and the radial Cauchy stress is then determined by integrating (15). Once the radial stress component T R R is determined, the circumferential component satisfies Note also that for constant γ R and γ θ , these integrals may be performed analytically, giving explicit expressions for the stress and deformation in terms of the growth. As we show later, the same holds when extending from one layer to multiple layers; if the growth in each layer is constant, the stress components may be written explicitly. It is this fact that we exploit below in formulating a discretized growth dynamics. This is the main motivating reason for the fixed core geometry we consider. Under different boundary conditions, the deformation and stress would be coupled, requiring for instance a root finding exercise to determine the outer radius for which the stress boundary condition is satisfied. In such a case, the framework below applies at the expense of added computational complexity.

Growth Law
We now impose a homeostasis-driven growth law of the form (1). In the plane polar geometry, this takes the forṁ Here K R R := K R R R R , K Rθ := K R Rθθ , K θ R := K θθ R R , K θθ := K θθθθ are the only non-vanishing components of the fourth-order tensor K , and are assumed to be constant in space and time.

Discretization Approach
For given homeostatic stress values and components of K , the growth dynamics is fully defined, with the growth components evolving according to (17). Even in the simplified cylindrical geometry, this comprises a system of nonlinear partial differential equations. Moreover, viewing the dynamics as a discrete process is still complicated by the fact that at each time step updating the growth requires knowing the stress components, which requires integration of (15), which requires integration of (7), which cannot be done analytically for general spatially dependent γ R and γ θ .
However, as stated above, for constant γ R and γ θ , the integrals determining stress may be computed analytically. This suggests a discretization process whereby the annular domain is divided into discrete layers, each with constant growth, and such that the growth in each layer evolves according to averaged values of the stress. In this way, analytical expressions may be determined for both the stress and the average stress, and hence, the dynamics is reduced to a set of ordinary differential equations for the growth components. The inhomogeneity of the full model is replaced by a piecewise homogeneous model. This preserves the key idea of inhomogeneity (allowing, for instance, circumferential growth to be higher near the nucleus than away from it), but is more analytically tractable and allows for precise statements about the long-term dynamics, stability, and qualitative investigation such as the influence of radial versus circumferential stress to the growth dynamics.

Kinematics
We first consider two elastic layers attached to a solid nucleus and in perfect mechanical contact at their interface. In the initial reference configuration B 0 , the inner wall has the radial coordinate R 0 = A 0 , the middle wall at R 0 = A 1 and the outer wall at R 0 = A 2 . In the current configuration B t , the same material points have coordinates that are r (A 0 ) = A 0 , r (A 1 ) = a 1 and r (A 2 ) = a 2 (see Fig. 3).
We impose that in the reference configuration the two annular layers enclose the same area πΔ 2 . The initial reference radii of the two rings thus satisfy The deformation follows the same equations formulated in Sect. 2.1, but with piecewise homogeneous growth where γ 1 and γ 2 are constant. Note that our convention is to use subscript to denote different layers and superscripts for the coordinate basis index. Here, we have imposed Fig. 3 Kinematic setup for the two-layer system. The innermost layer is attached to an unmoving nucleus (a 0 = A 0 ) and the boundary condition at the outer layer is no pressure T R R (A 2 ) = 0 isotropic growth, i.e., γ R 1 = γ θ 1 = γ 1 and γ R 2 = γ θ 2 = γ 2 . The same ideas apply for anisotropic growth, but this simplification reduces the dynamics to a 2D phase space for γ 1 , γ 2 . In principle, one could also have piecewise material properties and piecewise K values; however, our objective is to consider the dynamics in a reduced parameter space, hence the only distinction between the layers is the different growth rates.
The deformation in each layer comes from integrating (8), subject to r (A 0 ) = A 0 and r (A 1 ) = a 1 . We obtain Note that at R 0 = A 1 , r is continuous but not differentiable.

Mechanics
The stress balance (15) determines the radial stress as From (16), we then obtain the circumferential stress T θθ R 0 : The expressions T R R 1 and T R R 2 as well as T θθ 1 and T θθ 2 can be determined analytically as functions of A 0 , A 1 , A 2 , μ, γ 1 and γ 2 , though the exact expressions are long and have been suppressed here.
We note that the radial stress component is continuous at the interface between layers. This can be seen by evaluating The circumferential stress, however, is discontinuous at the interface unless the growth rates of adjacent layers are equal, γ 1 = γ 2 . See also Fig. 4.
Sample stress profiles for varying values of γ 1 (with γ 2 = 1) are given in Fig. 4. With γ 1 > 1, the inner layer grows uniformly; hence, its reference state is a uniformly expanded annulus; however, it is constrained by attachment to the core and to the ungrowing outer layer. Thus the inside of the inner layer is in radial tension (the inner

bottom) components of Cauchy stress for
edge is "stretched" radially to match the core), the outside is in radial compression, and the entire layer is in compression in the hoop direction. The outer layer, on the other hand, is forced to expand circumferentially to accommodate the growing inner layer and is in circumferential compression; this is balanced by a compression in the radial direction. The inverse effect occurs with γ 1 < 1.

Growth Law
We define the average stresses T 1 and T 2 , for both radial and circumferential stress components, as Our approach is to modify the growth dynamics so that the (constant) growth in each layer evolves according to the averaged stress values. That is, we study the systeṁ Note that the isotropic growth enforces K R R = K θ R and K θθ = K Rθ ; hence, there are only two (rather than four) growth rate constants K R R and K θθ . To further reduce the parameter space, we make the additional assumption that the homeostatic stress values are equivalent in layers 1 and 2, that is We emphasize that while T R R i and T θθ i for i = 1, 2 are averages over actual stresses according to (23), the homeostatic values T R R i * and T θθ i * for i = 1, 2 are prescribed values that may, but need not, correspond to averages of physically realizable stresses.
To facilitate the analysis, we rescale all stress quantities by a characteristic value σ , e.g.,T R R = T R R /σ , and rescale time ast = tσ K θθ . We also introducẽ The parameterK is a measure of anisotropy of the mechanical feedback, i.e., a weighting of the contribution of radial versus circumferential stress to the (isotropic) growth response. The rescaled growth law is theṅ Here we have re-defined the overdot as derivative with respect to the rescaled time, and we have dropped all hats for notational convenience. Note that all stress averages depend nonlinearly on γ 1 and γ 2 , but not on the spatial coordinate R 0 , which has been integrated out.

Stability Analysis
To investigate the behavior of the growth dynamics, we can now apply standard techniques of dynamical systems to (27); that is we seek equilibria satisfyingγ 1 = 0 anḋ γ 2 = 0 and compute their stability. Let γ eq 1 , γ eq 2 denote an equilibrium state. The nonlinear nature of the dependence of T R R 1 , T R R 2 , T θθ 1 and T θθ 2 on γ 1 , γ 2 makes it difficult to compute analytically the number and location of equilibrium states as a function of the parametersK and T * and we shall use numerical methods to this end.
For a given equilibrium state, we then perform a linear stability analysis. Let 0 < ε 1 and expand as Introducing γ = (γ 1 , γ 2 ) to describe the state of the system (27), its linearly expanded version (to order ε) takes the formγ where the Jacobian matrix has entries Stability is determined in the usual way by the form of eigenvalues of J, which are the roots of the characteristic equation

Bifurcation Diagram
The number of equilibrium states and their stability depend on the values ofK and T * . In Fig. 5a, we present a phase diagram that shows four regions with distinct dynamical behavior. These can be summarized as follows: -Region I has four equilibrium states, of which one is a stable node, two are saddles, and the fourth is either an unstable node or an unstable focus. -Region II has four equilibrium states: two are saddles and the other two are either stable nodes or a stable focus and stable node. A Hopf bifurcation at the interface of Regions I & II transforms the unstable focus into a stable focus. -Region III has two equilibrium states, one of which is a stable node, the other a saddle node. At the interface between Regions II and III, a saddle node bifurcation occurs that annihilates the stable node and saddle node in Region II. -Region IV has no equilibrium states.
In Fig. 5b, we show phase portraits for the selected points P1-P5. Nullclines are plotted as blue and green curves, illustrating the appearance and disappearance of equilibrium states as categorized above.
As is evident in Fig. 5, there is a wealth of possible dynamical behavior exhibited in this system. That an idealized two-layer model with isotropic growth and equivalent homeostatic values in each layer has such a rich structure highlights a more generic complex nature of mechanically driven growth. Our intent is not to fully categorize the behavior; rather this system should be seen as a paradigm to illustrate complex dynamics. Nevertheless, several observations are in order.
One observation from the phase portraits in Fig. 5b is that unbounded growth is not only possible but "common," at least in the sense that many parameter choices and initial conditions lead to trajectories for which γ i → ∞. Perhaps the most natural initial condition is to set γ 1 = γ 2 = 1, which corresponds to letting the system evolve from an initial state with no growth. Examining the trajectories in Fig. 5b shows that . 5 a Bifurcation diagram for two layered actively growing piecewise homogeneous system. b Equilibrium states and their dynamical characterization. Parameter values were A 0 = 1, A 1 = 1.562, A 2 = 1.970 points P1 and P2 would not evolve toward the single stable state, but rather would grow without bound. Another point of interest is that while regions I, II and III contain stable equilibria, the stable states in Regions I and III satisfy γ eq 1 γ eq 2 < 1. These are equilibria for which one of the layers has lost mass (at least one of the γ i < 1). Growth in both layers requires both γ i > 1, and we find that such an equilibrium only exists in a small subset of Region II, shaded dark blue in Fig. 5. We further see that T * < 0 in the dark blue region, andK approximately in the range 10-17. This implies that in order for a stable equilibrium to exist where both layers have grown, the homeostatic stress must be compressive in one or both components, and the system must respond more strongly to radial than to circumferential stress.
Admissible Versus Inadmissible Homeostatic Values. In Fig. 5, we imposed the homeostatic stress T * to be equal in each layer. Moreover, T * could take any value and thus had no direct correspondence to a physically realizable stress state. We now define an admissible homeostatic value as the average over a stress field that can be physically realized with the given geometry and boundary conditions. Such an admissible homeostatic stress state derives from a homeostatic growth, i.e., a given growth field γ * = γ * 1 , γ * 2 T defines a spatially dependent stress, and averaging according to (23) then gives admissible values for the homeostatic stress: An inadmissible homeostatic value is one that cannot be expressed as an average over an actual stress, i.e., there exists no γ * defining T * .

Growth Law with Admissible Homeostatic Values.
To conclude our analysis of the two-layer system, we return to the same growth law, but for admissible homeostatic values. Due to the spatial inhomogeneity of the stress profile in the two-layer cylinder (see for instance Fig. 4), it is not possible to have equal homeostatic values in each layer 1 and 2. The growth law with admissible homeostatic values readṡ The phase space for this system is now inherently three dimensional, as the homeostatic stress values are defined by the two choices γ * i as opposed to the single value T * . Here we restrict our analysis to a single example, with γ * 1 = 5.867, γ * 2 = 3, andK = 23.5, thus representing a preferred state defined by significant growth in each layer, and with strongly anisotropic growth dynamics due to the large value ofK . The dynamics are presented in Fig. 6. The contour plot in Fig. 6a shows that there are in total four equilibrium states. The streamlines and trajectory plots in Fig. 6b, c reveal that the equilibria consist of a stable spiral, two saddles, and one stable node. It is interesting to note that P4, which is the equilibrium state at which both γ eq i = γ * i , is unstable; that Trajectories and layer sizes for highly anisotropic growth law with admissible homeostatic state. a Contours forγ 1 = 0 andγ 2 = 0 for the system 6. As can be confirmed from the stream plots b, c, there is one stable spiral, two saddles, and one stable node. The saddle point P4 in (b) is the homeostatic equilibrium γ * 1 , γ * 2 . Parameters: μ = 2, Δ = √ 3 (A 0 = 1, A 2 = √ 7).K = 23.5. Homeostatic growth: γ * 1 = 5.867, γ * 2 = 3 is, the system does not remain at the equilibrium state through which the homeostatic values were defined. Included in Fig. 6b are three sample trajectories, with the size of each layer shown at different times, and illustrative of the variety of dynamical behavior. The green trajectory quickly settles to a stable state marked by significant resorption (both γ i < 1); the blue and red trajectories sit outside the basin of attraction of P1 and show an initial period of resorption followed by significant growth. The red trajectory is in the basin of attraction of the stable focus and thus oscillates between growth and decay as it approaches the stable point at P3, while the blue trajectory, just outside the basin of attraction, ultimately grows without bound, never reaching an equilibrium state.

Growth of Discrete N Layer System
Next, we generalize the dynamical system of the previous section from two to N layers where growth and stresses are constant throughout each layer. If N is sufficiently large, a system of N layers can be used as a suitable spatial discretization of a continuous growth profile on which precise statements can be obtained. In this case, we can generalize Eq. (33) to N coupled ODEs. We will analyze the stability of this system near a homeostatic equilibrium and show to what extent the results obtained for N = 2 remain unchanged as the discretization is refined (N increases), which informs the stability of the continuous (N → ∞) system.
A major difference compared to the two-layer model is the method to obtain homeostatic values. Previously, homeostatic values were prescribed via the homeostatic growth values γ * 1 , γ * 2 . In the present model, homeostatic values are obtained by assuming the existence of a prescribed continuous homeostatic growth profile γ * R 0 . The homeostatic values γ * i are then obtained through local averaging of the prescribed profile γ * R 0 over an interval by generalizing Eq. (23). These values are admissible by construction.
Since growth is taken as constant in each layer, the stresses can be determined fully analytically and a stability analysis can then be performed. The stability analysis will inform under which conditions the dynamical system will either relax to a homeostatic state after a small perturbation or lead to an instability.

Kinematics
We consider N perfectly connected annuli, separated by N + 1 interfaces, which in the initial reference configuration have the radial coordinate values {A 0 , A 1 , . . . , A N } as sketched in Fig. 7. The K -th annulus is defined by A K −1 ≤ R ≤ A K for K ∈ {1, . . . , N }. We choose a particular discretization so that the area between layers, πΔ 2 , is constant: We can write A K explicitly as Given a continuous curve γ R 0 , we define the piecewise constant growth profile by taking the average The growth value γ K is constant for all K . We demonstrate the construction of the discrete profile {γ K } from the continuous profile γ R 0 in Fig. 8, in which we consider Fig. 7 Kinematic setup for an isotropically growing N layered system. Note that the discretization is chosen such that the areas of each layer are equal as an example the continuous function Once {γ K } are obtained, we compute the radial map r K R 0 from the discrete profile {γ K }. Note that while γ K is a constant throughout the K -th layer, the radial map r K is a function of the radial coordinate R 0 : Explicitly, this implies Notice that the recursive expression (38) and the explicit expression (39) are consistent with the requirement which means that r is continuous at each interface A K −1 (Fig. 9).

Mechanics
Stress Components. In the continuous version, the radial stress T R R is obtained from (15). The discrete version reads Fig. 9 Radial function r K R 0 for the case of discrete growth γ i , computed according to (39). The dashed line represents the case of no deformation r = R 0 ; everything below the dashed line is resorption ("shrinking"), everything above this line is growth (Parameters as in Fig. 8) Traction continuity at the interfaces implies We define τ R R R 0 as the indefinite integral over the right hand side of (41) (dropping the integration constant), from which, we express the radial stress in the K-th layer as The circumferential stress T θθ is related to the radial stress T R R through (16). The discrete version of the relationship between T R R and T θθ is given by where Stress profiles corresponding to the growth law (37) are depicted in Fig. 10a (radial) and Fig. 10b (circumferential).
Average Stress. As in the two-layer case, average values for the radial and circumferential stress can be computed exactly. The average radial stress in the K -th layer T R R K is Fig. 10 Stress profile and stress averages for the growth profile (37). a Radial stress profile T R R and average stress profile T R R . The analytical curve was obtained from (44) and the numerical curve (for validation) was obtained from (15). In both the numerical and analytical case, the piecewise growth profile γ i according to (36) was used. The average stress was computed according to (47) with the same growth profile as the other curves. b Circumferential stress profile T θθ and average stress profile T θθ . The analytical curve was obtained from (45) and the numerical curve (for validation) was obtained from (16). The average stress was computed according to 49. All other parameters are as in Fig. 8, with Young's modulus μ = 1 where ν K R 0 is defined as We have seen in (45) how the circumferential stress T θθ relates to the radial stress T R R . The average over that expression is We have presented an expression for κ K in (46). The average over κ K is According to (49), the expression for T θθ K is the sum of κ K (see 50) and T R R K (see 47). The average radial and circumferential stress components are depicted as horizontal lines in the respective layers in Fig. 10a (radial) and Fig. 10b (circumferential).

Generating a Homeostatic State from a Prescribed Growth Profile.
The discretization and averaging process described above enables for a concise framework for studying growth dynamics. As the homeostatic state is defined by a growth profile-a function that is only constrained to be positive-a generic classification of dynamic behavior is likely untractable. Our intent, rather, is to briefly investigate stability and the rate of convergence in terms of the number of layers. For this, we restrict attention to a linear homeostatic growth profile γ * R 0 , characterized by a single parameter, C 1 , Note that this growth profile satisfies γ * (A 0 ) = 1, i.e., no growth at the inner boundary. We obtain the discrete homeostatic stress profile γ * i from the continuous profile γ * R 0 by computing the average according to (36). The homeostatic stress T (γ * ) is computed from the discrete homeostatic stress profile γ * i according to (44) and (45). The homeostatic values T (γ * ) are obtained as averages according to (47) and (49). It is important to note that the homeostatic stress is generated by prescribing a growth profile (51), which by definition ensures that the homeostatic stress is admissible.

Growth Dynamics
We consider a growth law that generalizes (33) to N layers. The main difference with (33) is that the values for homeostatic stress are obtained by the linear growth profile.
The growth law readṡ In order to consider the stability of (52) in the neighborhood of the homeostatic state, we expand growth around its equilibrium values: To linear order in ε, the dynamical system simplifies tȯ The eigenvalues of the Jacobian matrix J characterize the stability of (52) near the homeostatic state. The components of the N × N matrix J are We characterize the stability in the neighborhood of the homeostatic state as a function of two non-dimensional parameters: The mechanical feedback anisotropy parameter K and the slope of the homeostatic growth profile C 1 . The latter appears in (55) through γ * (see Sect. 4.3). Figure 11a shows a bifurcation diagram of the stability of the dynamical system (52) as a function ofK −1 and C 1 for N = 9 layers (note that unlike in Fig. 5, here we use the inverse ofK to focus on large circumferential stress). The regions are colored according to the largest real part of the eigenvalues λ i of J, that is λ = Max(Re λ 1 , Re λ 2 , . . . Re λ N ). There are three parameter regions: an unstable region (orange), a stable region (blue), and an undecidable region (green) for which λ is within a small tolerance of zero. This last region is included as it is typically within numerical error and its inclusion allows to make precise statements about stability. This relatively shallow region of λ is further explored in Fig. 12 and allows us to identify the clearly stable and clearly unstable regions of the diagram. Figure 11b-e shows that for increasing values of N (that is, a refinement of the discretization), the regions are practically unchanged (b-d), and that the largest eigenvalue of four selected points converges reliably to a finite positive (P1 and P2) or negative (P3 and P4) eigenvalue.
The green shallow region is more explicitly visualized in Fig. 12. This plot shows in the vertical axis the value of the largest real eigenvalue computed at K −1 , C 1 from the Jacobian matrix (55). The planes λ = tol and λ = −tol are shown in dark gray, and eigenvalues between are assumed to be in the shallow (green) region in which stability cannot be decided from an expansion of γ according to (53) to first order in ε.
Thus, we see that there exist a region of stability, and a region of instability, which both persist (for large enough N ) independently of the discretization. A strongly anisotropic growth law (K −1 close to zero or negative) is required for the system to be unstable. We also considered the convergence as N increases for a representative sample of points in the stable and unstable regions and confirmed that there was no significant change in λ. We expect that the stable and unstable regions represent the true behavior of the full (inhomogeneous) system discussed in Sect. 2. The intermediate (green) shallow region of eigenvalues has a more complicated structure due to the discretization that is not expected in the full system.

Solid Nucleus Versus Pressurized Cylinder
Our choice for boundary conditions has been a solid nucleus at the inner wall of the cylinder, r (A 0 ) = A 0 , and a no-pressure condition at the outer wall, T R R (B 0 ) = 0. These conditions are advantageous for mathematical analysis because the deformation map (8) can be easily integrated for a piecewise constant growth profile, with the result (39). Similarly, the radial stress (15) and circumferential stress (16) can be explicitly integrated for a piecewise constant growth profile, with stresses (44) and (45), respectively, and average stresses (47) and (49), respectively. Explicit forms of the average stresses are very advantageous for our model, in which derivatives of these expressions must be computed to obtain the Jacobian (55). In many biological tubular systems, however, the boundary condition will be a pressure difference between inner and outer cylinder walls, due to hydrostatic pressure maintained inside the tube.  Fig. 11. The large shallow region has a fine structure which is an artifact of discretization and is not expected in the full inhomogeneous system. For this reason, we choose a three color system in Fig. 11a-d, in which the shallow region and its fine structure are merged into one region defined by −tol ≤ λ ≤ tol. The two planes serving as upper and lower bounds of this region are depicted in dark gray. Values above λ = tol are stable, below λ = −tol are unstable Consider the following boundary conditions: where P is the prescribed pressure. Now, instead of (8), the integration of the radial map reads Notice that where previously we had the known nucleus position A 0 , we now have the unknown current inner cylinder coordinate a. Therefore, even for a piecewise constant growth profile, (57) cannot be directly obtained because a is unknown. To determine a, it is necessary to integrate the radial stress T R R while using the boundary condition (56) in the form −P = T R R (B 0 ) − T R R (A 0 ). This requires a numerical root finding algorithm (e.g., Newton's Method) to determine a, as shown in Ben Amar and Goriely (2005) for a spherical shell. As the difference between a solid nucleus boundary condition and a pressure difference boundary condition (56) lies only in a difference for the inner cylinder wall, we expect that the techniques outlined in this article (stress-averaging, N -layered cylinder) will apply to the pressurized cylinder setup, albeit at extra computational cost.

Conclusion
It is now well appreciated that growth can induce mechanical instabilities (Goriely and Ben Amar 2005;Ben Amar and Goriely 2005). The related problem that we have considered in this paper is the stability of a grown state through its slow-growth evolution. The question is not therefore about mechanical instability but about the dynamic stability of a preferred homeostatic state. While the former is characterized by a bifurcation from a base geometry to a more complex buckled geometry, occurring on a fast elastic timescale, the latter involves the system evolving away from a given stress state on the slow growth timescale. In general the homeostatic state is not homogeneous; hence, the issue of stability requires the analysis of partial differential equations defined on multiple configurations with free boundaries. There are no standard mathematical tools available to study this problem even for simple nonhomogeneous systems. An alternative is to consider the stability of states that are piecewise homogeneous (in space). The problem is then to establish the stability of coupled ordinary differential equations describing locally homogenous states through the traditional methods of dynamical systems. Within this framework, we considered two relatively simple problems. First, we considered the dynamical stability of a two-layer tube with different, but constant, growth tensors in each layer. We characterized the dynamics of the full nonlinear system and showed that the number of equilibria and their stability varies greatly and gives rise to highly intricate dynamics which we organized via several bifurcations. We identified a parameter region where the system is stable. We found that the growth dynamics of tubular structures in the neighborhood of the homeostatic equilibrium depends in a nontrivial way on the anisotropy of the growth response and that the equilibrium becomes unstable for highly anisotropic growth laws. This complexity of dynamics naturally raises the question about stability of homeostatic equilibria for more general systems.
Second, we showed that given a continuous law in a cylindrical geometry, we can introduce a suitable discretization of the problem that keeps all the characteristics of the continuous problem. We showed that for a linear growth law, there are clear regions where stability and instability persist independently of the discretization (for sufficiently large N ). We expect that these regions represent the true behavior of the full inhomogeneous system. This result allows us to characterize the stability of a morphoelastic growing cylinder.
As our results (in particular Fig. 11) show, admissible homeostatic states can lead to either mechanically stable or unstable equilibria. This suggests a way to distinguish between physiological (stable) and pathological unbounded (unstable) growth. Indeed, our model also suggests a natural growth termination mechanism. The question of growth termination (what triggers a tissue to stop growing?) is a much-debated question in developmental biology. It has been particularly well studied in the model system of the Drosophila melanogaster wing disk where the morphogen Decapentaplegic (Dpp) has been identified as the main regulator of growth. A number of models propose growth regulation and termination based on a combination of mechanical effects and Dpp concentration. Some of them are continuum models (Aegerter-Wilmsen et al. 2007;Ambrosi et al. 2015), others are vertex models (Shraiman 2005;Hufnagel et al. 2007) but none of them is entirely satisfactory (Vollmer et al. 2017). We hope that the dynamical stability of homeostatic states offers an alternative way of looking at growth termination, which emerges naturally in our model as a stable equilibrium.
While we have only scratched the surface of the complex dynamic behavior that exists in such systems, the framework presented here provides a tool to explore growth dynamics and stability of homeostatic states and finally address some of the fundamental challenges of morphoelasticity (Goriely 2017): What growth laws, in general, would lead to dynamically stable homeostatic states? What is the final size of a growing organism for a given growth law? What are the conditions under which growth dynamics produces oscillatory growth?