Post-buckling behaviour of a growing elastic rod

We consider mechanically-induced pattern formation within the framework of a growing, planar, elastic rod attached to an elastic foundation. Through a combination of weakly nonlinear analysis and numerical methods, we identify how the shape and type of buckling (super- or subcritical) depend on material parameters, and a complex phase-space of transition from super- to subcritical is uncovered. We then examine the effect of heterogeneity on buckling and post-buckling behaviour, in the context of a heterogeneous substrate adhesion, elastic stiffness, or growth. We show how the same functional form of heterogeneity in different properties is manifest in a vastly differing post-buckled shape. Finally, a fourth form of heterogeneity, an imperfect foundation, is incorporated and shown to have a more dramatic impact on the buckling instability, a difference that can be qualitatively understood via the weakly nonlinear analysis.


Introduction
Mechanically-induced pattern formation is a phenomenon prevalent in the morphogenesis of many biological structures, from airway wall remodelling [27], to wrinkling of skin [17], to blades of grass [13].The prevailing feature in such systems is the deformation from a 'trivial' base state to a more complex geometry, with buckling induced by mechanical stress.This feature is not unique to biological systems; indeed the same basic wrinkling pattern found in an elephant's skin can be induced by compressing a sheet of rubber.What is unique to the biological world is that such patterns tend to form without any external influence, rather the stress needed for mechanical instability is produced internally.The primary origin for this stress is differential growth, i.e. different parts of a tissue growing at different rates.Perhaps the simplest example is a tissue layer that grows relative to an underlying substrate to which it is adhered.The growth induces a compressive stress in the growing layer, and at some critical threshold the tissue buckles, exchanging compression energy for bending energy.This situation underlies the formation of numerous biological patterns, including brain development-induced by the differential growth between the cortex and subcortex [8]; intestinal crypt fission, in which epithelial tissue grows but is tethered to underlying tissue stroma [38]; even seashell ornamentation, characterised by the adhesion of the growing mantle tissue to the rigid shell that it secretes [10].Aside from differential growth, inherent geometrical constraints can also induce buckling, for instance a row of cells growing uniformly but within a closed space will similarly develop compressive stress and ultimately buckle.
A mathematical description of growth-induced mechanical buckling can take a variety of forms.A number of discrete cell-based models have been devised, in which mechanical interactions between individual cells, coupled with cell growth or proliferation, can lead to deformation in the form of folds [15], invaginations [31], or protrusions [9,24].This approach is more amenable to the inclusion of cell-level biological detail.Continuum modelling, while less amenable to this level of detail, allows one to utilise analytical tools for differential equations, which may improve insight and reveal parametric relationships that are more difficult to attain from discrete models.Within 3D elasticity, growth is naturally incorporated via decomposition of the deformation gradient tensor into a growth tensor describing the local change of mass and an elastic tensor accounting for the elastic response [34].However, beyond simple geometries such as the buckling of a sphere [5], a 3D description of buckling typically requires fully computational techniques such as finite element methods [3].In many cases, the geometry under consideration is well-suited for a reduced dimensional analysis.This is clearly true in the case of filaments.Filaments by definition have one length scale much longer than the other two and hence are well suited to a 1D description.Kirchhoff theory for elastic rods has been applied to a diverse range of filamentary systems, such as DNA coiling [36], neurite motility [33], plant tendril twisting [20], and many more.A planar rod description may also be relevant even in an inherently 2D sys-tem: for instance when a sheet of tissue deforms approximately uniformly in a transverse direction, or when a cross-section of tissue deforms such as in the circumferential wrinkling of a tube [26,6,4].
Within a continuum formalism, a commonly-considered problem is a planar rod on an elastic foundation or substrate, typically with the ends fixed.The basic premise is that growth (axial elongation) of the rod generates compressive stress, creating buckling from a flat to a curved state, but this deformation is resisted by elastic tethering to a fixed substrate.This basic setup, or very similar, has been studied for many years in an engineering context (where 'thermal expansion' typically plays the role of 'growth'), dating back to classic works of Timoshenko [37] and Biot [7], and continues to find new interest and applications.
It is only more recently that the relation to biological pattern formation has become clear and similar systems have been specialised for biological problems.An important aspect that distinguishes biological systems from the above examples is the various ways in which growth can occur; a key challenge here is connecting a continuum level description of growth to underlying cell-level processes.This has stimulated extensive mathematical modelling development and has created the need for a systematic framework.There are three such works of particular relevance for the present paper.In Moulton et al. [28], the theory of Kirchhoff rods was extended by incorporating growth in a manner inspired by the morphoelastic decomposition.This is the framework upon which our analysis is built.Also of note are recent descriptions of buckling in the context of intestinal crypt formation, invaginations that are present throughout the intestines.Edwards and Chapman [16] applied a continuum mechanics approach to the formation of a single crypt.They modelled the crypt epithelium and its underlying tissue stroma as a beam upon a viscoelastic foundation.By performing a linear stability and eigenvalue analysis of buckling, they examined the effect of changes to proliferation, cell death, adhesion, or motility.Nelson et al. [29] complemented this analysis with a 'bilayer' model representing an epithelial layer connected to a flexible substrate.Nelson et al. also conducted an eigenvalue analysis similar to Edwards and Chapman and combined this with a numerical analysis of the full, nonlinear model, demonstrating the influence of heterogeneity in both growth and bending stiffness on the resulting buckled crypt shape.Similar systems, but in 2D using plate theory, have also explored pattern formation due to growth instabilities.Hannezo et al. [21] characterised transitions from crypt-like to herringbone and labryinth patterns; Nelson et al. extended their 1D models in [30] and found that crypt patterning could be most strongly controlled through heterogeneity in growth.
Our objective in this paper is to use the morphoelastic rod framework of Moulton et al. [28] to extend the results of [16,29] and analyse unexplored features that are of general relevance.A focal point for our analysis is the behaviour of the system beyond the initial buckling.In an engineering context, buckling will typically signify failure, and so the threshold value to induce buckling may be the most relevant quantity.For pattern formation in biology, on the other hand, the shape evolution well beyond the initial instability is often critical to the final pattern (and its biological functionality), and hence analysis only of the onset of instability is insufficient.
Also of relevance in many biological systems is understanding the role of heterogeneity in mechanical pattern formation.Heterogeneity can arise in three main forms: growth1 , the mechanical properties of the rod, or substrate adhesion.Here it is important to distinguish between growth and remodelling.Growth refers to an increase (or decrease) in mass, i.e. a change in size of a tissue layer without any change in its material properties.Remodelling, on the other hand, refers to a change in material properties without any change in mass, e.g.due to fibre reorientation or cell differentiation.In a growing tissue, both of these processes occur and will commonly occur non-uniformly.The crypt, for instance, is not a layer of uniform cells, but rather consists of a clear proliferative hierarchy of cells with varying rates of division as one moves up the crypt axis [39].Heterogeneity in adhesion may occur due to non-uniform changes in the substrate layer, or in a biological context due to changes in the cells, or may occur due to buckling itself, for instance due to viscoelastic effects.
Both of these features-large deformation beyond buckling and heterogeneitypose significant mathematical challenges.To capture post-buckling behaviour requires analysis of a nonlinear system of equations, as opposed to the linear stability analysis that can be used to detect buckling.Furthermore, including heterogeneity complicates the use of many analytical tools, either rendering the system analytically intractable or complicating attempts to unfold bifurcations.Here, rather than rely fully on computational techniques, our approach is to analyse post-buckling behaviour and the effect of heterogeneity through a combination of a weakly nonlinear analysis and numerical solution.This approach yields a broad understanding of the role of heterogeneities in growth, material properties, and adhesion, and reveals features of post-buckling pattern formation not described in previous analyses.
We consider a 1D model system of a growing planar rod on an elastic foundation, serving both as an extension of the classical setup and as an abstracted framework for several of the aforementioned biological systems.The rod is subject to growth in the axial direction and clamped boundary conditions, which drive buckling at a critical growth.The goal of this paper is to understand the factors driving the onset of buckling and the post-buckling behaviour.In particular, we investigate how the buckling and post-buckling behaviour changes in the presence of spatial heterogeneity in material properties, obtaining explicit relations for how the pitchfork bifurcation that arises is impacted by heterogeneity, and exploring the shape evolution in the nonlinear post-buckled regime.
The remainder of this paper is structured as follows.In Section 2, we outline the general theory for Kirchoff rods and incorporation of growth, as developed in [28].Then, in Section 3 we summarise the linear stability analysis before extending to a weakly nonlinear analysis.The results of the weakly nonlinear analysis and numerical analysis of the full nonlinear model are presented in Section 4, first in a homogeneous setting, then with the addition of different material heterogeneities.Finally, we close by discussing the implications of our results and directions for future model extensions.

Model setup
As a model system to investigate post-buckling in growing slender structures, we consider in this paper an extensible and unshearable planar rod in quasi-static mechanical equilibrium.The rod is constrained geometrically by clamped-clamped boundary conditions, and is also adhered to an elastic (Winkler) foundation.Growth of the rod is modelled under the morphoelastic rod framework of [28].In this framework, one identifies three distinct configurations: an initial reference configuration parametrised by the arc length S 0 ; a grown virtual configuration, parameterised by S and referred to as the 'reference' configuration; and the current configuration, parametrised by s.This framework is summarised in Figure 1.In the initial, reference, and current configurations, the total rod lengths are L 0 , L, and l, respectively.The rod arc length is assumed to evolve through an axial growth process, described by the growth stretch γ(S 0 ) = ∂S/∂S 0 , followed by an elastic response, encapsulated by the elastic stretch α.The total stretch of the rod λ from initial to current configuration is then given by This is the 1D analogue of the multiplicative decomposition of the deformation gradient tensor employed in 3D morphoelasticity [34,5,28].
Let the rod's centerline be given by (x, y), and let θ denote the angle between the tangent vector, d 3 = cos θe x + sin θe y , and the x-axis.This is expressed by (The scale factor αγ accounts for the fact that we parametrise the system by the initial arclength parameter S 0 , as opposed to the current arclength parameter s.) Defining the resultant force and moment in the rod by n = F e x + Ge y and m = me z , the balance of linear and angular momentum give Here f = f e x + ge y is the external body force, which is assumed to be solely due to the underlying foundation.
Fig. 1: The different morphoelastic rod configurations.A rod that is initially confined to a (finite) interval grows in a virtual, unstressed reference configuration, before being mapped to the current configuration, where it is subject to boundary conditions and loads.The rod is parametrised by a different arc length in each configuration.The respective rod lengths have been indicated.The parameters γ, α and λ denote the growth, elastic, and total stretches, respectively.
To this system we add a constitutive equation relating moment to curvature: The parameter E is the Young's modulus of the rod, and is taken for now to be constant, and I denotes the (second) moment of inertia, respectively.For an extensible rod, we take the axial stress to be related to the elastic stretch α via a linear constitutive relation where A is the cross-sectional area of the rod.The foundation is assumed to be a linearly elastic medium occupying an interval along the x-axis, as is the rod.Initially the foundation is a distance y 0 from the rod centreline and the rod is glued to the x-axis; that is, a point (S 0 , 0) along the x-axis is attached to a point (S 0 , y 0 ) on the rod.We can set y 0 = 0 without loss of generality.No remodelling takes place, so these two points are still connected in the reference configuration, and are now at (S/γ, 0) and (x, y) respectively.Therefore the body force acting on the rod (as a force per initial length) is where we assume the foundation stiffness is proportional to that of the rod, with the dimensionless positive parameter k comparing the stiffness of the foundation to that of the rod.The factor of 1/γ indicates that the body force is parametrised with respect to the initial configuration, and that no remodelling occurs after growth.
The system is closed with the clamped boundary conditions at S 0 = 0 and L 0 :

Non-dimensionalisation
Next, we non-dimensionalise the system using the standard Kirchhoff scaling [12,19,18] and circumflexes to denote non-dimensional quantities: Dropping the circumflexes of independent and dependent variables for notational convenience, Equations ( 2)-( 4) simplify to The constitutive law for extensibility now reads The dimensionless boundary conditions are where the remaining (non-dimensional) parameters are Note that the non-dimensional rod length L 0 depends on the ratio of two length characteristics of the rod: its initial total length L 0 and the thickness, characterised by (I/A) 1/2 .For instance, for a rod with circular cross-section of radius r, we have I = πr 4 /4, A = πr 2 , and hence L 0 = 2L 0 /r 1.

Stability Analysis
In this section, we present the analytical tools that we will use to investigate the buckling and post-buckling behaviour of the morphoelastic rod.We first adapt and summarise the linear stability analysis from [28], used to calculate the growth bifurcation value, γ * , before unfolding the bifurcation with a weakly nonlinear analysis.

Linear stability analysis
We first determine the critical growth stretch γ * and corresponding buckling mode using a linear stability analysis.The calculations in this section are also present in Moulton et al. [28], but are summarised here for completeness and to motivate the weakly nonlinear analysis.Inspecting the system (10)-( 14), for all γ > 1 there exists a base solution corresponding to a straight, compressed rod; that is, with θ ≡ 0 and the total stretch λ = 1: Expanding each variable about the base solution so that, for example, x = x (0) + δx (1) + O(δ 2 ), where δ is an arbitrary small parameter, and considering O(δ) terms, leads to the linearised system F (1) − kx (1) = 0, G − ky (1) = 0, ( θ (1) − γm (1) = 0, m where = ∂/∂S 0 .We observe that {x (1) , F (1) } are decoupled from {y (1) , G (1) , m (1) , θ (1) }, and that ( 17)-( 19) can be expressed as two ordinary differential equations for x (1) and y (1) , M y (1) := y (1) + 2ay (1) where the coefficients a and b are defined by Solving Equation (20) subject to x (1) (0) = x (1) ( L 0 ) = 0 leads to the trivial solution x (1) ≡ 0 (and subsequently F (1) ≡ 0).Turning to Equation (21), our linearised boundary conditions are Seeking solutions of the form y (1) ∼ e iωS0 (valid on an infinite domain) leads to the following oscillation frequencies We remark that in order for non-damped oscillations to exist over the (finite) domain, we require a ≥ b.Applying the four boundary conditions for y (1) , specified by Equation ( 23), yields the solution where the constants C 2 and C 3 are given by and the critical growth value γ * must satisfy the relation If ( 27) is satisfied, then the buckled solution is given by ( 25), but with arbitrary constant C 1 .For later convenience, we thus define the function ŷ to be the determined part of this function, i.e.
The smallest value of γ > 1 that satisfies (27) occurs when a = b (and, hence, where k = kI/A 2 .Although it appears that γ * inf does not vary with any length scale, this is not entirely true.For example, for a rod with circular cross-section of radius r, then from (15), k ∝ k.However, for a rectangular cross-section with height h and width w, k ∝ kh/w.Additionally, this value of γ only leads to oscillations if the rod length L 0 is infinite, and results in the trivial solution over a finite domain.Therefore the critical growth value γ * is the first value of γ > γ * inf that solves (27).

Weakly nonlinear analysis
Having summarised the above results from [28] that are key for us, we now carry out a weakly nonlinear analysis.For a given root of ( 27), the buckled solution is only determined to within the arbitrary constant C 1 .In order to understand the behaviour of the buckled rod as it continues to grow, it is necessary to determine how the buckling amplitude C 1 depends on γ, which we accomplish through a weakly nonlinear analysis.We unfold the bifurcation by introducing the ansatz where ε is a fixed small parameter (whereas δ is small, but arbitrary) and γ (1) = O( 1) is a control parameter describing the proximity to the growth bifurcation point γ * .Substituting ( 30) into ( 12), re-expanding θ and m and retaining higher-order terms in θ reveals that the nonlinearities will be balanced by growth if ε = O(δ 2 ).Setting ε = δ 2 and re-expanding our variables about the trivial solution ( 16) as a power series in δ leads to the following system of differential equations for each order O(δ n ), n ≥ 1, Here, the functions h , and h m (n) denote inhomogeneities due to lower order terms.As was the case in Section 3.1, the system decouples into two linear operators acting on x (n) and y (n) , For general n ≥ 1, the boundary conditions for y (n) are now given by Observe that the homogeneous problems, ( 20) and ( 21), along with the boundary conditions-x (1) (0) = x (1) ( L 0 ) = 0 and ( 23), respectively-are self-adjoint.
When n = 1, we recover the linearised system described by Equations ( 17)- (19), i.e.H x (1) = H θ (1) ≡ 0. At n = 2 we find that H θ (2) ≡ 0, giving us no further information on the buckling amplitude C 1 ; however, H x (2) is nonlinear in θ (1) and thus a non-trivial solution for x (2) exists, as x (1) is trivial and the Fredholm Alternative Theorem is immediately satisfied.Hence, we must consider O(δ 3 ) terms to obtain the amplitude equation for C 1 .This leads us to consider the inhomogeneities H y (3) , which can be expressed solely in terms of x (2) and y (1) and their derivatives.Then, by the Fredholm Alternative Theorem, a solution for y (3) exists if and only if the following solvability condition is satisfied: where ŷ solves the homogeneous problem, defined by (28).The only unknown quantity in (37) is the coefficient C 1 ; hence this solvability condition for y (3)  provides a relation between the buckling amplitude C 1 and the distance from the critical buckling growth parameter, expressed by γ (1) .Simplifying this condition (see Appendix A) leads us to deduce that The constants K 1 , K 2 are tedious to compute analytically, but nevertheless only depend only on the material parameters k and L 0 , and through them the critical growth γ * .It can also be shown (Appendix A) that K 2 > 0 for all parameter choices.Equation (38) shows that the system exhibits a pitchfork bifurcation, a known property of similar systems [29,23].The three branches of the pitchfork are given by This shows that the bifurcation will be supercritical if K 1 < 0, and subcritical if K 1 > 0. In the next section, we explore the dependence of K 1 on k and L 0 and its effect on the buckling and post-buckling behaviour.

Buckling and post-buckling behaviour
Having established a relationship for the post-buckling amplitude, we now explore the effect of material parameters and heterogeneity on the buckling and post-buckling behaviour.First we examine the form of bifurcation in a homogeneous setting, effectively by analysing how the critical buckling growth γ * , the buckling mode, and the pitchfork constants K 1 and K 2 vary with the two free parameters in the non-dimensional system, L 0 and k.We then adapt the weakly nonlinear analysis to incorporate heterogeneity and investigate the impact of non-uniformity in foundation stiffness, rod stiffness, and growth.In each case, we complement the analytical work by solving the full nonlinear system (10)-( 14), using the numerical package AUTO-07p [14].AUTO-07p uses pseudo-arclength continuation to trace solution families and solves the system with an adaptive polynomial collocation method.

Effect of length and foundation stiffness
In Figures 2-3, we plot bifurcation diagrams for varying values of the dimensionless foundation stiffness k and the dimensionless rod length L 0 , respectively, for increasing growth.The horizontal axis in each case is the growth parameter γ, and the vertical axis plots the non-trivial branches, ± y := ±max S0 |y(S 0 )|, which is closely related to the value of the constant C 1 but more representative of the post-buckling amplitude.The solid lines are determined from the weakly nonlinear analysis, while the dashed lines are numerical results.We also plot the buckled shape (x(S 0 ), y(S 0 )) at the specified points for each branch.
As k increases, the buckling occurs for increased mode number, reflecting the energy trade-off that as the foundation stiffness is increased, a large amplitude is penalised more by a high foundation energy, and hence higher bending energy is sacrificed to have a lower amplitude.An increased value of k also leads to an increase in γ * , which shows that the foundation is serving to stabilise the rod against buckling.This can again be understood in terms of an energy trade-off, but in this case it is the compressive energy in the grown but unbuckled state that is being sacrificed.It is important to note that this feature could not occur in an inextensible rod, for which the trivial state does not exist for any γ > 1.
Considering length L 0 leads to similar changes in both critical buckling growth and mode number.However, while the buckling mode increases for increasing length, the critical buckling growth decreases.Recalling the scaling L 0 = L 0 (A/I) 1/2 , this reflects the notion that a short or thick rod can endure more growth before buckling, and will buckle at lower mode.We note also that as L 0 → ∞, γ * → γ * inf , the critical growth value for buckling on an infinite domain (29).
Perhaps most notable is the transition from supercritical to subcritical bifurcation evident in both diagrams.We find that subcritical bifurcations occur for large enough k or small enough L 0 .We find through numerical continuation that the subcritical branches then fold back, a feature not captured by the weakly nonlinear analysis at order O(δ 3 ).A linear stability analysis (see Appendix B) confirms that the portion of the subcritical branch before folding back is unstable, while the portion after the fold-back is stable.(As would be expected, the curved branches are stable in the supercritical case.)This implies that a subcritical bifurcation signifies a discontinuous jump from the trivial flat state to the finite amplitude stable branch, as well as the presence of a hysteresis loop if γ is subsequently decreased.

Locating the pitchfork transition
For given parameters k and L 0 , the weakly nonlinear analysis enables us to determine the type of pitchfork bifurcation simply by computing the sign of K 1 .Figure 4 shows the regions in k-L 0 space where supercritical and subcritical pitchforks occur.Note that despite having an explicit expression for K 1 , the actual computation of its value was done numerically as it requires root finding for the eigenvalue γ * .Hence, to produce Figure 4 we computed K 1 over a discrete grid in the k-L 0 plane.The transition boundary was then verified at several points through numerical path following in AUTO-07p.
Unexpectedly, we do not find a simple monotonic transition boundary, as seen in similar studies [22], but rather an intricate pattern with an oscillatory structure that, to our knowledge, has not been reported before.This structure implies that multiple transitions between super-and subcritical buckling can occur for a fixed k and varying L 0 (or vice versa); that is, simply increasing the length of the rod monotonically can create repeated transitions between super-and subcritical bifurcation.For an infinite rod, the transition can be computed as k ≈ 0.38196 (see B); this point is included as a dashed, horizontal line in Figure 4, and it appears that as L 0 → ∞, the oscillations dampen and the transition boundary approaches this constant value.By contrast, as L 0 decreases, the oscillation amplitude increases, although the slenderness assumption of the rod breaks down as L 0 ∼ O(1).Whether the non-monotonic structure in the intermediate region, which still has oscillations of moderate amplitude, has an intuitive explanation, or whether it can be observed physically, are interesting open questions, but ones that we postpone for future study.

Parameter heterogeneity
Thus far, we have assumed spatial homogeneity in model parameters.However, in many biological systems, heterogeneities are inherent in the system.This raises the question of how a given heterogeneity is manifest in the buckled pattern.To explore this, we initially consider three distinct forms of heterogeneity in: the foundation stiffness, the rod stiffness, and the growth.Each heterogeneity has a clear biological interpretation.For example, in the intestinal crypts, these heterogeneities would correspond to: the different types of extracellular matrix secreted by the cells comprising the underlying tissue stroma (founda- tion heterogeneity), the mechanical properties of epithelial cells in the crypt (stiffness heterogeneity), and the variations in the proliferative capacity of these cells (growth heterogeneity).In order to understand the effect of each type of spatial heterogeneity, we examine heterogeneity for each parameter in isolation.
With parameter heterogeneity, it becomes increasingly difficult to obtain analytically tractable results with a weakly nonlinear analysis, especially if the amplitude of the heterogeneity is pronounced.Nevertheless, when the parameters are close to homogeneous, we can extend the weakly nonlinear analysis and, in particular, ask how the heterogeneity impacts the pitchfork bifurcations observed in the homogeneous case.We complement this analysis with numerical solutions of the full system defined by Equations ( 10)- (14).For computational convenience, heterogeneity is incorporated via a sequence of numerical continuations in the growth and heterogeneity parameters.
We model heterogeneity as a spatial deviation from a baseline homogeneous state.In general, for an arbitrary parameter µ, we consider Here, the constant µ 0 corresponds to the baseline homogeneous value, and the function ξ(S 0 ) captures the spatial variation, modulated by the amplitude factor and constrained only by the requirement that µ ≥ 0. With heterogeneity introduced, the weakly nonlinear analysis can be viewed as a two-parameter unfolding, both in the distance from the critical buckling growth, via γ = γ * + εγ (1) , and in the distance from homogeneity, characterised by .We are thus faced with balancing three small parameters: ε, , and the order of the expanded variables, which we denoted by δ, e.g. as in y = δy (1) + O(δ2 ).In the homogeneous case the correct balance is given by = δ 2 .With > 0, numerous balances could be sought, and a full analysis of the two-parameter unfolding is beyond the scope of this paper.Our approach involves starting from homogeneity, and increasing the order of to see when and how it first impacts on buckling.Hence, we again take = δ 2 , and consider = δ β η, with β > 1 and η an O(1) control parameter.
The three cases we wish to consider for heterogeneity are: -Foundation heterogeneity, for which µ = k; -Rod stiffness heterogeneity, for which µ = E, the Young's modulus 2 ; Perturbing each of these parameters via (40) has a similar effect on the weakly nonlinear analysis.In each case, it is easy to show that for β > 2, the heterogeneity does not affect the weakly nonlinear analysis up to O(δ 3 ), and therefore does not affect the buckling amplitude C 1 .Consequently, the bifurcation relation (39) is unaffected.When β = 2, the heterogeneity first has an impact (up to O(δ 3 )) and, hence, it is for this case that we adapt the analysis.At O(δ 3 ), the corrective term y (3) now satisfies The first term on the right hand side, H old y (3) , corresponds to the inhomogeneities in the homogeneous case, while H new y (3) describes the effects of the heterogeneity (40).
For each parameter considered, evaluating the solvability condition (37) leads to a new equation for C 1 : The constants K 1 and K 2 are identical those in Equation (38) (see A).The heterogeneity is fully encapsulated in the term K 3 , defined straightforwardly by The heterogeneous model, hence, undergoes a translated pitchfork bifurcation, where the branches are given by Observe that since ξ appears in H new y (3) only, and hence in K 3 only, it does not affect the type of pitchfork that occurs, but merely translates it.That is, setting C 1 = 0 in the non-trivial branch gives γ (1) = −(K 3 /K 2 )η; recalling (30), the critical growth γ * is now shifted to where γ * 0 is the critical growth stretch for the rod in a homogeneous setting, as determined from the linear stability analysis in Section 3.1.Hence, we see that the material heterogeneity (40) results in an O(δ 2 ) shift in γ * .Since K 2 > 0, the direction and degree of the shift is determined by K 3 .Generally, K 3 provides the 'metric' for whether the heterogeneity has a net effect of strengthening or weakening the effect of the material parameter.
In order to investigate greater amplitudes of heterogeneity and the postbuckling shape evolution, we perform numerical continuation on the full model.
As an illustrative example, we apply the same form of heterogeneity for each of the three parameters: ξ(S 0 ) = cos(2πS 0 / L 0 ) and = 0.9, characterising a significant decrease in the middle region and increase in the outer regions.Figure 5 depicts the resultant rod shapes.As evident in Figure 5, the heterogeneity has a markedly different effect for each material parameter.
Foundation stiffness heterogeneity Here, the modified foundation is softer in the middle and stiffer near the endpoints, causing a significant increase in amplitude in the middle of the rod, where the resistance to deformation is weaker.This phenomenon can be generally understood and quantified by applying the weakly nonlinear analysis.Note from Equation ( 21) that the operator M on the left hand side of Equation ( 41) is the linearised (beam) equation, i.e. the vertical force balance for an extensible rod upon a foundation.Consequently, the term H new y (3) captures additional forces due to the imposed heterogeneity.In the case k(S 0 ) = k 0 (1 + δ 2 ηξ(S 0 )), at O(δ 3 ) this term takes the particularly simple and instructive form: The heterogeneity thus acts as an amplifying force where ξ(S 0 ) < 0, and a resistive force where ξ(S 0 ) > 0. This is apparent in Figure 5(b): the magnitude of y (1) is largest in the middle, with ξ(S 0 ) < 0, reducing the effects of H old y (3)   and leading to an increase in amplitude.
Rod stiffness heterogeneity In the case of rod stiffness, the dominant trend is compression in the middle region, leading to a significant decrease in amplitude and arclength, and formation of a near cusp-like point, reflecting the reduced energy cost of both bending and stretching in the middle region.We have also examined the competing energies within the system: bending versus stretching versus foundation (defined in E).We find that both the bending and foundation energy are reduced as increases, despite the cusp-like formation, while the stretching energy increases (see Figure 10).The total energy remains roughly constant through most of this tradeoff, but eventually, at large values of , the stretching penalty outweighs the benefit to the bending and foundation energies and a sharp rise in the total energy occurs for 0.7.

Growth heterogeneity
In the case of growth heterogeneity, note that for the form of heterogeneity considered, ξ(S 0 ) = cos(2πS 0 / L 0 ), the net growth, defined by is unchanged from the homogeneous case, γ(S 0 ) = γ 0 .Thus, for varying , there is no change in net growth, merely a redistribution of material from for each of (b) foundation stiffness heterogeneity k(S 0 ), (c) rod stiffness heterogeneity E(S 0 ), and (d) growth heterogeneity γ(S 0 ).The heterogeneity amplitude has been continued to = 0.9 from the homogeneous state ( = 0).
the middle region to the sides.Accordingly, Figure 5(d) shows a significant change in shape: the middle region flattens while the left and right regions, with increased material, show an increase in amplitude and curvature.The loss of material from the middle also has the effect of increasing the elastic stretch α, resulting in a transition from compression to tension.An intuitive explanation for this can be seen by examining the flat solution ( 16): a growth stretch of γ < 1 implies that the horizontal force F (0) > 0, i.e. the rod is in a state of tension.For growth heterogeneity, we observe behaviour in energy that is qualitatively opposite from rod stiffness heterogeneity: as heterogeneity is increased, both bending and foundation energies increase, while stretching and total energies decrease.The flattening of the middle region caused by loss of material to the edges reduces the compressive energy locally, while the redistribution of material to the sides leads to a net increase in bending and foundation energy (Fig. 11).

Inverse problem
Figure 5 shows that applying the same form of heterogeneity to different material properties can result in radically different shapes.Two natural and related questions follow from this: (i) can one tailor the heteregeneities to achieve a desired shape?And (ii) given a particular shape, can one infer the form and type of any material heterogeneity present?These questions, with significant relevance both from morphogenetic and tissue engineering perspectives, are related to the mathematical inverse problem.Such an inverse problem is inherently complex, as the shapes considered are (partial) solutions of a high order nonlinear boundary value problem, only achieved in the forward direction through numerical path continuation.In order to develop some intuition, here we consider a simple but illustrative example.Our approach is as follows: we take a candidate shape with embedded heterogeneity-foundation heterogeneity following Fig 5(b)-and we try to match that shape by varying the heterogeneity in the rod stiffness and growth (as well as the net growth); we then consider characteristics other than the shape itself and seek distinctive differences, i.e. signatures of the heterogeneities.
We then examined all other features (e.g.moment m, force components F , G, etc), seeking distinguishing characteristics, from which we identified two salient features.Figure 6(b) shows the axial stress n 3 in the outer region.For growth heterogeneity (dashed-dot blue curve), we find a marked increase in n 3 compared to the other two heterogeneities.To produce the desired shape with heterogeneous growth required redistribution of material from the outside regions to the middle, thus creating a region of axial tension.In contrast, the rod remains in compression for foundation and rod stiffness.In principle, this difference could be detected in a cutting experiment, as for instance used to determine residual stress in arteries [11] and solid tumours [35].Only the situation with heterogeneous growth would spring open in such a region, as tension is released.
In distinguishing foundation stiffness heterogeneity, the foundation energy density U F , defined explicitly in E, provided the clearest indicator.Figure 6(c) plots U F in the middle section of the rod, where the shapes are qualitatively most similar, and we find a significant decrease in the case of foundation heterogeneity.This is perhaps not surprising, reflecting the fact that a locally softened foundation results in lower foundation energy, but there is a clear dif-Fig.6: Characterisation of different material heterogeneities.The baseline foundation stiffness k 0 and rod length L 0 have been set to k 0 = 0.04 and L 0 = 20, respectively.(a) Different values of and different choices of ξ(S 0 ) have been set for foundation stiffness heterogeneity (red, solid), rod stiffness heterogeneity (brown, dashed), and growth heterogeneity (blue, dot-dashed) respectively to produce similar post-buckled shapes (x(S 0 ), y(S 0 )).(b) Despite the similar shape and amplitude in this region, growth heterogeneity results in regions of low growth resulting in tension (n 3 > 0).(c) Foundation stiffness heterogeneity decreases the foundation energy density U F where the foundation has been weakened, even though the different heterogeneity cases exhibit similar shapes and amplitude.ference in magnitude despite the similarity in rod amplitudes.In principle, this could also be detected in a different cutting experiment: separating the growing cells (rod) from substrate (foundation) would generate a less pronounced substrate recoil in the case of foundation heterogeneity.However, it must be noted that there is only a comparative difference, and one that may not exceed experimental error in a real biological context.Moreover, we observed no features that clearly distinguished the case of rod stiffness heterogeneity from the other two.Hence while we have identified potential signatures, this simple example also illustrates the inherent challenge in resolving the inverse problem.

The role of extensibility.
It is important to note that many of the above trends are reliant on the assumption of rod extensibility.In an inextensible model, axial compression is not permitted, as the arclength is fixed, which is equivalent to the geometric constraint α ≡ 1.For explicit comparison, we consider the same stiffness heterogeneity in an inextensible rod.Consequently, only bending is affected by the heterogeneity (40).In Figure 7, we compare the shape evolution with increasing in both inextensible and extensible models, for ξ(S 0 ) = cos(πS 0 / L 0 ) and ξ(S 0 ) = cos(2πS 0 / L 0 ).For the inextensible models, we take an equivalent foundation stiffness, k = 0.04, but set γ = 1.1 to obtain a similar initial amplitude.In an inextensible rod, the arclength is fixed and thus the response to heterogeneity is to alter the shape towards aligning points of minimal and maximal curvature with material points of maximal and minimal stiffness, respectively.Hence in Figure 7(a) the inextensible rod shifts to have maximal amplitude on the soft region on the right side, whereas the extensible rod (Fig. 7(b)) compresses on the right side, thus producing a completely different morphology with minimal amplitude.In Figure 7(c), the inextensibility leads to a localisation of curvature in the soft middle region, as opposed to the strong compression in the extensible case, shown in Figure 7(d).These simulations illustrate the dramatic effect that extensibility can have on shape morphology and the response to material heterogeneity.

Foundation imperfection
The heterogeneities we have considered thus far, while having significant impact on the post-buckling shape evolution, have had a relatively minor effect on the bifurcation itself, only serving to translate the pitchfork, and by modest amounts.This is in contrast to typical results in a Koiter imperfection sensitivity analysis [23,2], in which material imperfections may shift the bifurcation to occur at significantly reduced loads and in an imperfect fashion (a 'broken pitchfork').Here, the small change in bifurcation can be understood within the framework of our model by considering the form of heterogeneity imposed.As derived in section 3.2, the base (homogeneous) equation for the pitchfork bifurcation is 1) = 0. Perturbations to the system in the form of heterogeneities have the potential to change this to The cases we have examined lead to K 4 = 0 and K 3 = 0, which merely translates the pitchfork (as C 1 = 0 is still a solution branch).Breaking the pitchfork would require K 4 = 0.The reason that the additional term obtained has a factor of C 1 is that we have only considered multiplicative heterogeneity, i.e. we The resulting forms of E(S 0 ) (brown, dashed line) have also been plotted.Arrows in the plots indicate the evolution of the rod shape (blue, sold lines) in the increasing direction of the continuation parameter, while darker blue lines correspond to higher values of .When rod stiffness is asymmetric, competition in curvature causes the inextensible rod to redistribute its material more so than the extensible rod.In the symmetric case, extensibility leads to compression at the locations of maximal curvature, which is not seen for the inextensible rod.
have imposed heterogeneity in terms that multiply dependent system variables: foundation stiffness k multiplies x and y in the force balance, stiffness E multiplies θ as well as α, and growth γ appears in the system multiplicatively in multiple places (as evident in Equations ( 2)-( 4)).Thus, at the relevant order in an asymptotic expansion, a perturbation to these parameters always appears multiplicatively with the base solution y 1 = C 1 ŷ, and thus the additional term in the solvability condition that provides the bifurcation condition is of the form K 3 C 1 .
In order to produce a non-zero added term K 4 , independent of C 1 , we must consider additive heterogeneity.One possible type of additive heterogeneity, commonly considered in imperfection analyses, is in the shape of the foundation; that is, we consider the foundation to have spatially-varying imperfections present.That is, we modify the force balance equations (11) to where ξ(S 0 ) is the shape of the imperfection and captures the magnitude.This form of heterogeneity first affects the weakly nonlinear analysis when = O(δ3 ).Setting = δ 3 η, where η acts as a control parameter away from homogeneity, the solvability condition ( 37) is shifted by a factor independent of C 1 , and the bifurcation now satisfies The constant K 4 is defined by (The constants K 1 and K 2 are the same as in the homogeneous case (38).)Note the loss of both the trivial amplitude branch and the symmetric nature of the non-trivial amplitude branches.Therefore, with this underlying imperfection, the model undergoes an asymmetric (or imperfect) pitchfork bifurcation, with the branch selected determined by the sign of K 4 .As η is increased, so is the deviation from the homogeneous amplitude equation (38), and hence the splitting of the initially-symmetric non-trivial branches is amplified.Note also that K 4 involves a simple inner product with ξ and the buckling mode y.
The heterogeneity thus has maximum effect when the imperfection to the foundation is of the same shape as the buckling mode, i.e. when ξ ∝ y. Figure 8 displays the bifurcation diagram for ± y against γ for various values of η, both from the weakly nonlinear analysis (dashed curves) and numerical solution of the full system (solid curves).We have taken the heterogeneity ξ(S 0 ) = y(S 0 ).As expected, increasing η further splits the branches and increases the predisposition to the upper branch 3 .Due to the scaling = δ 3 , large values of η are needed to observe a noticeable difference in the bifurcation.In Figure 8 we have taken η = O(10 5 ), which grossly violates Fig. 8: Bifurcation diagram for additive foundation heterogeneity (48).The foundation stiffness and rod length have been set to k = 0.04 and L 0 = 20, respectively, while the foundation heterogeneity ξ(S 0 ) has been prescribed to ξ(S 0 ) = y(S 0 ), as defined by Equation (28).The amplitudes ± y = ± max S0 |y(S 0 )| against γ from weakly nonlinear analysis (blue, dotted lines) and numerical continuation (red, solid lines) have been plotted for different values of η.For numerical calculations, the small parameter δ is set to δ = 0.01 for all cases.Increasing η further biases the rod to γ = 1 to the upper non-trivial amplitude branch.the asymptotic assumption that η = O(1); since in this plot δ = 10 −2 , even with η = 10 5 , = 10 −1 , i.e. the perturbation is still small, and we find that the weakly nonlinear analysis matches the full model reasonably well.Also evident is that as heterogeneity is increased, a larger deflection is observed for a given γ.This is consistent with typical results that imperfection deforms load-deflection curves so that higher deflections occur under smaller loads [23].
As a final point of interest, we wish to measure which form of heterogeneity has the greatest impact on the bifurcation.As a 'metric' for comparison, following the typical engineering analysis of load-deflection, here we consider the compressive force as a function of growth.To examine this, we define the net axial stress Fig. 9: Effect of the considered heterogeneities on compressive stress.
The net axial load n 3 is plotted against γ.As before, the foundation stiffness and rod length have been set to k = 0.04 and L 0 = 20, respectively, while the foundation heterogeneity ξ(S 0 ) has been prescribed to ξ(S 0 ) = y(S 0 ) (see Eq. ( 28)).The homogeneous case (blue, solid line) is compared against the foundation stiffness heterogeneity (orange, dashed lines); rod stiffness heterogeneity (pink, open circles); growth heterogeneity (green, dash-dotted lines); and foundation imperfection (48) (brown dots).For numerical calculations, the small parameter δ is set to δ = 0.01 and η is set such that = 0.1 for all cases.The additive foundation heterogeneity is quickest to relieve the axial stress induced by rod growth.
In Figure 9, we compare the net axial stress with increasing growth for each of the four heterogeneities considered.In each case we have imposed ξ = ŷ, and due to the different nature of the perturbation schemes, we have chosen the scale factors such that the total perturbation from the uniform state is equivalent across the four cases.The perfect buckling case appears as the solid blue line, with the sharp cusp appearing at γ * and signifying that buckling occurs at a critical compressive stress, which is relieved partially through the buckling.Each of the heterogeneities produces a similar curve, though we see that the additive heterogeneity in foundation shape just considered has the most significant effect, followed by growth heterogeneity.Both foundation and stiffness heterogeneity follow the perfect case very closely; zooming in on the cusp region (see inset) shows that these forms lead to a delayed bifurcation, and, in the case of foundation stiffness, the bifurcation occurs at slightly larger stress before subsequently compensating and dipping below the perfect case.

Discussion
We have investigated the buckling and post-buckling behaviour of a planar morphoelastic rod attached to an elastic foundation.We extended the original linear stability analysis by Moulton et al. [28] by conducting a weakly nonlinear analysis, complemented with numerical solutions of the full, nonlinear model.We first considered a homogeneous setting, and then explored the effect of heterogeneity in material parameters.
In the homogeneous case, we obtained a classic pitchfork bifurcation, with buckling occurring at a critical growth.The nature of the bifurcation (its location and type-supercritical or subcritical) could be characterised via two dimensionless parameters, one ( L 0 ) relating to length of the finite rod, and another ( k) comparing the relative stiffness of foundation and rod.Increasing length was found to destabilise the rod, causing bifurcation at a smaller value of growth and with increased mode number.Increasing the foundation stiffness, on the other hand, stabilises the rod, increasing the critical growth and the mode number.These results are consistent with previous similar analyses [32].The type of bifurcation, however, was non-standard; the boundary between supercritical and subcritical bifurcations exhibited an unexpected complexity.In a biological context, where monotonically increasing growth is a natural driver of the formation and subsequent evolution of spatial patterns, this transition has critical importance, signifying where a smooth shape evolution (supercritical) can be expected.Here the effect of a finite domain is also apparent, as the complexity of the transition becomes less pronounced as L 0 increases.
Multiplicative heterogeneity with respect to three different material properties was then considered: the foundation stiffness, the rod stiffness, and growth.A modified weakly nonlinear analysis showed that in each case the heterogeneity served to translate the bifurcation point, but did not alter its nature.Explicit relations for the shift in bifurcation allowed us to determine how the form of the heterogeneity influences the direction and degree of the translation.For example, the simplest relation appeared with heterogeneity in the foundation stiffness, in which the greatest effect occurs when the heterogeneity is aligned the square of the buckling mode.This reflects the intuitive notion that weakening the foundation attachment in regions where the uniform rod deforms maximally has the strongest impact.
To complement the weakly nonlinear analysis, the full nonlinear system was solved with numerical continuation; this enabled us to investigate the postbuckling for more pronounced heterogeneity and at large growth values.A common feature was an induced 'asymmetry' of the buckled shape.With heterogeneous foundation stiffness, softer (stiffer) parts of the foundation give rise to increased (decreased) rod amplitudes, as might be expected.With heterogeneity in rod stiffness, the situation is less straightforward.Softer parts of the rod are more easily curved, and thus it might be reasonable to expect such regions to correspond to higher amplitude; however, compression is also less costly in the soft regions.In all cases we have examined, the rod flattens through compression in the soft regions, a deformation that increases stretching energy, but is compensated by a decrease in both bending and foundation energies.Here, the assumption of extensibility is crucial, as compression is not permitted in an inextensible model.Indeed, a direct comparison of an inextensible and extensible model post-buckling revealed significant morphological differences, highlighting the importance of a critical assessment of when the inextensible assumption is warranted.
In the case of non-homogeneous growth, we showed that even with zero net growth, heterogeneity, interpreted as a redistribution of rod material from spatial regions with decreased γ to those with increased γ, can significantly impact the post-buckling behaviour.The general trend is not surprising: the rod flattens in regions where material is lost.What is perhaps surprising is that the distribution of material seems to play as important a role in the shape evolution as the total amount of material added through growth.
Having examined the effect of heterogeneity on post-buckled shape, we investigated the extent to which these heterogeneities can be characterised.Growth heterogeneity produces disparate regions of tension (where growth is reduced) and compression (where growth is increased).These regions, and in principle the heterogeneity itself, are amenable to detection through cutting experiments that release residual stresses.Heterogeneity in foundation stiffness is manifest through regions of increased amplitude and decreased foundation energy.Again, such regions may in principle be experimentally detectable through release of compression, in this case by separating the rod from the foundation.
While these thought experiments suggest the possibility of distinguishing between forms of heterogeneity and using heterogeneity to tailor properties, it is clear that this is not a straightforward problem, even in the idealised planar rod setting we have considered.In a 3D setting, more measurable quantities are available, for example, stress in transverse directions, which could potentially yield measurable differences in behaviour.On the other hand, the general complexity of the inverse problem will increase as the number of variables increases.
In the final section, we have examined a fourth type of heterogeneity, with a view to establishing why the impact of heterogeneity on the bifurcation itself was relatively minor in the previous scenarios.Here we made the key distinction between multiplicative and additive heterogeneity.A multiplicative heterogeneity appears in a term that multiplies dependent variables in the system; due to the nature of the perturbative expansion, such terms only serve to shift the pitchfork bifurcation.An additive heterogeneity, for which a perturbation is applied to a term that does not multiply dependent variables, can have a significant effect, creating an imperfect bifurcation (broken pitchfork) and creating a larger deviation from the perfect, homogeneous, case.Here we considered an imperfectly-straight foundation, and showed that the effect is maximal when the form of the imperfection matches that of the buckling mode.
It is worth discussing how these results compare to similar studies.In Nelson et al. [29], bending stiffness heterogeneity and growth heterogeneity were considered.In the case of the former, the authors found that localised regions of softened bending stiffness leads to a localisation of buckling.This result is similar to Figure 7 for the inextensible case, where the rod shape shifts towards points of softened rod stiffness.Here, it is important to note the significance of the extensibility assumption; for an extensible rod, we found that compression is the dominant factor in altering rod shape.In the case of growth heterogeneity, Nelson et al. concluded that net growth affects the post-buckling behaviour more than heterogeneity.This discrepancy can be attributed to the authors' prescription of growth as an increase of the total rod length at the boundaries, averaging over any spatial variations in local growth.Here, we considered growth as a local increase of arc length, showing that heterogeneity (without change to net growth) leads to drastically different post-buckled shapes.Interestingly, these results on growth are in qualitative agreement with a later study by Nelson et al. [30], where they extended the model using 2D plate theory to consider a monolayer tethered to a substrate and foundation.They concluded that in this 2D setting, growth heterogeneity is capable of generating more complex buckled shapes than in the earlier 1D model.Furthermore, while substrate heterogeneity is capable of quickening the onset of buckling, it is subservient to growth heterogeneity in pattern formation.From our weakly nonlinear analysis, we have shown that buckling can actually be quickened or delayed by material heterogeneity, and have provided generalised metrics that quantify the resulting shift in bifurcation, a particular benefit of the weakly nonlinear analysis.Moreover, considering the inverse problem has revealed that it is possible to produce similar post-buckled shapes through different heterogeneity types.
In this paper, we have assumed each model parameter to be independent from the others and, for the sake of clarity, varied each parameter in isolation.A natural extension would be to consider the combined effect of several simultaneously-varying parameters via inter-parameter coupling.For example, one could consider a rod with non-uniform stiffness and a growth evolution law that depends on axial stress.This would naturally induce heterogeneity in multiple parameters, and the resultant competing effects would likely produce a complex solution space.Another natural extension is to consider nonlinear constitutive effects.There would certainly be benefit to considering a nonlinear constitutive relation between axial stress and the elastic stretch α, in particular because many of our simulations featured significant compression potentially beyond the threshold for quantitative validity of the Hooke's law considered here.Another useful extension is to incorporate nonlinearity in the response of the foundation, a phenomenon that has been studied in great detail in systems without growth [22], but whose role in the context of growth remains unclear.
Finally, we note that many of our modelling choices were motivated by observations on the intestinal crypt (and other, physiologically similar structures).Thus while the model represents an idealised version of a crypt, there are several extensions that would render it biologically realistic.For instance, our growth parameter contains no information about the timescale of growth, which is a fundamental aspect of many biological systems, particularly the crypt.Therefore, one could introduce time-dependent growth or time-dependent mechanical relaxation (for example in the foundation), allowing remodelling to occur over time.Alternatively, the proliferative structure within a crypt suggests the spatial form of the growth stretch should be bimodal [1,39].In the context of mechanosensitive growth, Miyoshi et al. [25] showed that a specific subset of stromal cells is activated during wound healing to increase stem cell proliferation in the crypt, as one example.The crypt also provides a natural setting to investigate possible feedback mechanisms between growth and the underlying foundation; this work is currently underway.
Using the Fredholm Alternative Theorem and considering (67) in powers of C 1 , the constants K 1 and K 2 in Equation ( 38) are obtained by evaluating the following integrals where x := C −2 1 x (2) and y := C −1 1 y (1) .Applying integration by parts to (69) yields Therefore K 2 > 0 for all parameter values, provided that y is non-trivial.Hence, the sign of constant K 1 determines the nature of the pitchfork.

B Pitchfork bifurcation on an infinite domain
We show that in the case of an infinitely-long rod, the value of k for which the system transitions from a supercritical pitchfork bifurcation to a subcritical pitchfork bifurcation can be calculated exactly.Moreover, we show that the transition occurs only once.
At O(δ 2 ), we substitute (72) into (58) and, after solving the equation and imposing boundedness, obtain At O(δ 3 ), we obtain the amplitude equation for C 1 by substituting (72) and ( 74) into (67) and impose that secular terms vanish.This leads to Equation ( 38), K 1 C 3 1 + K 2 C 1 γ (1) = 0, where the constants K 1 and K 2 are given by K 1 = k 2 7 + 100 k + k 2 We note that we have substituted Equation (29) to express K 1 and K 2 in terms of k only.As K 2 > 0 ∀ k > 0, the sign of K 1 completely determines the type of pitchfork bifurcation the system undergoes, as confirmed in Appendix A. Therefore, the value of k > 0 where the transition occurs is found by solving K 1 = 0, yielding It can be verified easily that K 1 < 0 when k < k * and K 1 > 0 when k > k * , indicating a transition from a supercritical pitchfork bifurcation to a subcritical pitchfork bifurcation.Moreover, there is only one transition when the rod is of infinite length, in contrast to the multiple transitions that can occur when the rod length is finite. ( dyn θ (3) dyn A solution to Equation (101) exists if and only if the determinant of the left hand side matrix is equal to zero for a given σ 2 .Therefore, if there is a value of σ 2 < 0 such that the determinant vanishes, then the buckled solution xeq is unstable, as the time perturbation grows exponentially as T → ∞.

E Energies
Here, we give the forms of the energy functionals that are considered in Sections 4.3.We then show how these defined energies change for heterogeneity in rod stiffness and growth, helping to elucidate the observed changes to rod shape in Figure 5.The total energy of the system, E, comprises contributions from bending, stretching, and the underlying foundation.
In the reference configuration, these dimensional individual energy densities may be written respectively as: After nondimensionalisation, we write the energy densities as where E and k represent the dimensionless rod and foundation stiffness repectively (note that in the homogeneous case, E = 1).Therefore, the total energy density is given by The latter two terms correspond to the geometric constraints placed on the rod, with the horizontal and vertical forces F and G acting as Lagrange multipliers.However, these contributions will vanish after energy minimisation.Hence, the total energy E total is given by where Figures 10-11 illustrate how the different energetic contributions change for increasing heterogeneity amplitude , in rod stiffness heterogeneity and growth heterogeneity, respectively.As increases, we observe tradeoffs between bending and foundation energy, and stretching energy.In Figure 10, the stretching energy increases, while bending and foundation energy decrease.This tradeoff helps to explain why compression dominates in the case of rod stiffness heterogeneity.In the case of growth heterogeneity, the behaviour is qualitatively opposite, due to the transition to tension in the middle region of the rod.Fig. 10: Energy plot for heterogeneous rod stiffness.The bending energy E bend (blue, solid), stretch energy E stretch (green, dashed), E foundation (brown, dotted), and total energy E total (pink, dot-dashed) have been plotted.The parameters are the same as those in Figure 5.We see that E stretch increases while E bend and E foundation decrease for increasing heterogeneity. .The bending energy E bend (blue, solid), stretch energy E stretch (green, dashed), E foundation (brown, dotted), and total energy E total (pink, dot-dashed) have been plotted.The parameters are the same as those in Figure 5.As growth heterogeneity increases, E stretch decreases drastically, with E bend and E foundation increasing slightly.

Fig. 2 :
Fig. 2: Bifurcation diagram for varying foundation stiffness.Bifurcation diagram plotting ± y = ±max S0 |y(S 0 )| against γ from weakly nonlinear analysis (blue, dotted lines) and numerical continuation (red) of the full model (10)-(14); stable solutions are marked with solid lines, while dashed lines represent unstable solutions.Each non-trivial branch denotes a continuation in γ for a fixed value of k.The buckled shape (x(S 0 ), y(S 0 )) is plotted at the point y = 1 (orange dot) for (a) k = 10 −3 , (b) k = 0.15, (c) k = 0.375, and (d) k = 0.625.The dimensionless rod length is fixed to be L 0 = 20 for all cases.

Fig. 4 :
Fig.4: Phase diagram of pitchfork bifurcations.The regions have been determined by computing the sign of K 1 at each point on a discretised grid of k and L 0 .Subcritical pitchfork bifurcations (K 1 > 0) have been labelled with dark blue crosses, while supercritical pitchfork bifurcations are labelled with an orange dot (K 1 < 0).The dashed line corresponds to the transition value of k in the infinite-length case, k ≈ 0.38196 (see B).

Fig. 5 :
Fig. 5: The effect of heterogeneity on the post-buckled shape (x(S 0 ), y(S 0 )) and the underlying foundation.The baseline foundation stiffness k 0 and rod length L 0 have been set to k = 0.04 and L 0 = 20, respectively.(a) The growth stretch γ (γ 0 for γ(S 0 )) has been continued until y = 2.75 for the homogeneous case.The heterogeneity function ξ(S 0 ) has been set to ξ(S 0 ) = cos 2πS0 L0

Fig. 7 :
Fig. 7: The effect of extensibility on rod shape.The specified heterogeneities are (a)-(b) ξ(S 0 ) = cos πS0 L0 and (c)-(d) ξ(S 0 ) = cos 2πS0 L0 .The dimensionless foundation stiffness and rod length have been set to k = 0.04 and L 0 = 20, respectively.The growth parameter γ was set to γ = 1.1 and γ = 1.8 for the inextensible and extensible cases respectively.Continuation in is over the interval ∈ [0, 0.85].The resulting forms of E(S 0 ) (brown, dashed line) have also been plotted.Arrows in the plots indicate the evolution of the rod shape (blue, sold lines) in the increasing direction of the continuation parameter, while darker blue lines correspond to higher values of .When rod stiffness is asymmetric, competition in curvature causes the inextensible rod to redistribute its material more so than the extensible rod.In the symmetric case, extensibility leads to compression at the locations of maximal curvature, which is not seen for the inextensible rod.