Mechanics of invagination and folding: hybridized instabilities when one soft tissue grows on another

We address the folding induced by differential growth in soft layered solids via an elementary model that consists of a soft growing neo-Hookean elastic layer adhered to a deep elastic substrate. As the layer/substrate modulus ratio is varied from above unity towards zero we find a first transition from supercritical smooth folding followed by cusping of the valleys to direct subcritical cusped folding, then another to supercritical cusped folding. Beyond threshold the high amplitude fold spacing converges to about four layer thicknesses for many modulus ratios. In three dimensions the instability gives rise to a wide variety of morphologies, including almost degenerate zigzag and triple-junction patterns that can coexist when the layer and substrate are of comparable softness. Our study unifies these results providing understanding for the complex and diverse fold morphologies found in biology, including the zigzag precursors to intestinal villi, and disordered zigzags and triple-junctions in mammalian cortex.

We address the morphogenetic folding induced by differential growth in soft layered solids via an elementary model that consists of a soft growing neo-Hookean elastic layer adhered to a deep elastic substrate. We find a transition from supercritical smooth folding followed by cusping of the valleys to direct subcritical cusped folding as the layer/substrate modulus ratio is varied from above unity toward zero. Beyond threshold the high amplitude fold spacing converges to about four layer thicknesses across a range of modulus ratios. In three dimensions the instability gives rise to a wide variety of patterns of cusped folds, including almost degenerate zigzag and triple-junction patterns that can coexist when the layer and substrate are of comparable softness, providing a quantitative understanding for the complex fold morphologies found in soft tissues.
PACS numbers: 46.32.+x, 62.20.mq, 87. 19.lx Some biological tissues have evolved into folded morphologies to accommodate large surfaces in small volumes; folding of the mammalian brain is associated with expansion of its outer cerebral cortex, and the folded epithelium in the gut provides large surface area for absorption of nutrients. Pattern formation in morphogenesis is typically attributed to chemical instabilities within Turing's reaction-diffusion paradigm [1], but increasing experimental evidence shows that these folded morphologies actually emerge from smooth precursors via mechanical instabilities driven by differential growth of layered tissues [2,3]. Recent studies have considered mechanical descriptions of brain folding [4][5][6][7][8], epithelial folding in the gut [3,9], fingerprint formation [10], and circumferential buckling in tubular organs and tumors [11][12][13].
We examine the stability of a layer occupying 0 < z < a atop a substrate occupying z < 0. The layer has undergone growth, swelling, or pre-strain such that, without adhesion, it would undergo a relaxing deformation G = diag(g x , g y , g z ), becoming a flat slab with thickness h = a/g z . We model both tissues as incompressible neo-Hookean elastic solids so that, if this reference state is subject to a displacement u and corresponding deformation gradient F = I + ∇u, the elastic energy per unit volume is 1 2 µ 2 Tr F · F T in the substrate and 1 2 µ 1 Tr F · G −1 · G −T · F T in the layer, while incompressibility requires Det(F ) = Det(G) = 1. We first consider uniaxial growth, G = diag(g, 1, 1/g), and consider the layer's linear stability to infinitesimal isochoric perturbations, u = f (z) sin(kx)ẑ + (f (z)/k) cos(kx)x. The full calculation is in the supplementary material [31]; here we restrict ourselves to an overview. The linearized Euler-Lagrange equations for the solids require g 4 f (4) − k 2 (g 4 + 1)f (2) + k 4 f = 0 (with g = 1 in the substrate), solved by f (z) = a 4 e kz +a 3 e −kz + a 2 e kz/g 2 +a 1 e −kz/g 2 in the layer and f (z) = e kz (b 1 +zkb 2 ) in the substrate. Imposing material continuity at the interface and natural boundary conditions at the free surface and interface we conclude the layer becomes linearly unstable when g satisfies Fig 1a plots this threshold as a function of wavelength for different modulus ratios. Minimizing this threshold over wavelengths gives the wavelength and threshold of the first unstable mode as a function of η, plotted in Fig 1b-c. As expected the instability reproduces the long-wavelength wrinkling limit for However, the zero wavelength Biot instability is never observed: it is always preceded by a point of nonlinear instability to infinitesimal cusped sulci [29]. In incompressible neo-Hookean materials this occurs at a uniaxial stretch λ s ≈ 0.647 [29], corresponding to g s ≈ 1.55. Thus we expect finite wavelength wrinkling as predicted by linear stability for η > η s ≈ 1.14, where the linear threshold intersects g s (see Fig. 1c), and direct sulcification at g s for η < η s . We verify these predictions and explore folding beyond threshold with numerical simulations [31], starting with uniaxial growth in two dimensions. When η > η s the surface indeed buckles smoothly at the linear instability threshold. This buckling amplifies surface compression in the valleys, producing the sulcification strain there shortly thereafter. The valleys then sulcify in a secondary bifurcation which occurs before g s up to η ≈ 10 ( Fig.  2a). Cusps form first on the surface and then on the interface between the layer and substrate at high g (Fig.  2c). However, the interface in real systems may not be strictly sharp, suppressing interfacial cusping.
When η < η s sulcification sets in at g s , before the linear instability. For η near η s the energy of the system at g = g s is minimized by a cusped fold that has a finite amplitude (see the layer with η = 1 and g = 1.55 ≈ g s in Fig. 2c), implying that sulcification is a nonlinear subcritical instability when η is near η s (Fig. S1). In the soft limit, η → 0, the substrate provides a rigid constraint for the layer and the situation is similar to the pure sulcification [29,30]; past g s sulcification proceeds as a nonlinear supercritical instability, forming sulci with a wavelength comparable to the thickness of the layer. Simulations indicate supercritical sulcification for η < ∼ 0.7.
We optimized the wavelengths of the folds for minimum energy (Fig. 2b) finding that in the smooth-folding regime, η > η s , the wavelength l decreases past threshold (in Eulerian coordinates), which is similar to folding of stiff layers [14]. In contrast, when η < η s the optimal wavelength increases past threshold. Consequently the optimal wavelength settles to about four times the undeformed layer thickness over a wide range of η when g becomes large. This agrees with the folding patterns in the brain and the gut, and is a hallmark of mechanical compressive folding. The near-threshold fold energy depends weakly on the wavelength for direct sulcification, η < η s , whereas there is a well defined wavelength in the smooth buckling regime. Consequently the folding patterns are expected to be more sensitive to imperfections when η < η s . The low η folds are also notably asymmetric, with a deeply cusped upper surface and a wavy bottom surface, in contrast to nearly top-bottom symmetric folds at η > η s (Fig. 2c).
We now turn to consider folding in three dimensions. If G = diag(g x , g y , 1/(g x g y )) is not uniaxial, but we nevertheless consider x − z plane-strain deformations, the elastic energy depends on two dimensionless parameters, η/g 2 x and ηg 2 x g 2 y . Our plane-strain calculations thus hold with the identifications η → g y η and g → g x √ g y . In particular, for transversely isotropic growth, g = g x = g y , the sulcification threshold is at g x √ g y = g 3/2 ≈ 1.55, i.e, g s ≈ 1.34, and the linear threshold intersects g s at η s ≈ 0.86.
To simulate three dimensional folded states we start with transversely isotropic growth of the layer on a large domain (width 30h), and allow the layer to fold and relax spontaneously (Figs. 3a and S2). These simulations show trends similar to the uniaxial case; when η > η s ≈ 0.86 the layer buckles smoothly, forming a pattern of hexagonal dents near the threshold. When η < ∼ η s , a pattern of line-like sulci with alternating directions forms (Fig. S2). These patterns have previously been identified as the near threshold patterns of stiff layers (η η s ) [18] and Folding of a uniaxially growing layer in two dimensions. (a) Morphologies as a function of modulus ratio η and growth ratio g = gx. The diagram includes the linear instability curve (blue line, dotted for η < ηs ≈ 1.14 where the sulcification threshold is met before the linear instability) and numerically determined thresholds for cusp formation (black line). (b) Energy minimizing wavelengths as a function of g for a range of moduli ratios η. Hollow and solid circles correspond to smooth and cusped folds, respectively. The blue curve is the wavelength of linear instability (the dotted end is above the sulcification threshold gs ≈ 1.55). (c) Simulated minimum energy states near threshold (left column) and far from threshold (right column) for modulus ratios η = 0.5, 1 and 4. Color indicates isotropic stress: dark red ≤ −5µ, green = 0, dark blue ≥ 5µ.
clamped soft layers (η = 0) [30], respectively (although analyses based on the plate theory predict a checkerboard pattern in the stiff film limit [17,18]). In contrast to these comparatively ordered near-threshold patterns, the far-from-threshold patterns that result from unaided relaxation are disordered, reminiscent of the complex fold patterns on the surface of the brain when η ≈ 1 (Figs. 3a,b).
To understand the emergence of the disordered farfrom-threshold patterns, we recall some previous results: for stiff films the zigzag (herringbone) pattern has been identified as an energetically preferred far-from-threshold pattern [18] and for a clamped soft layer we previously identified a hexagonal array of triple-junction sulci as an optimal far-from-threshold pattern [30]. There should thus be a transition from zigzags to triple-junctions as the layer stiffness is reduced from large toward zero. To investigate this we simulate periodic patterns on domains that contain one unit cell of the pattern with energetically optimized wavelengths. We use external guiding forces to inititate the patterns [31] that, once formed, are stable in the neighborhood of η,g-space where they are the energy minimizing patterns. Optimal near-threshold patterns, including a rectangular array of line-like sulci at η ≈ η s and smooth hexagonal dents at η > η s , are shown on the left of Fig. 3c. These patterns minimize the energy only at / c close to unity, where we define the relative growth strain / c = g−1 gc−1 to conveniently quantify distance from the threshold (g c is the growth ratio at the threshold). For far-from-threshold states ( / c ≥ 2) we construct zigzag patterns as well as patterns of hexagonal triplejunctions, discovering that these two patterns are almost equally good energy minimizers (Fig. 3e). However, the triple-junctions are slightly better at low η and/or high / c , while zigzags become relatively better toward large η, as expected. Fig. 3d shows a phase diagram summarizing these findings.
Given the nearly equal energies of the zigzag and triplejunction folds, it is not surprising that systems with modest modulus ratios, such as the folded cerebral cortex or its gel layer mimics [8], show disordered mixtures of zigzags and triple junctions (in contrast, epithelium of the gut shows ordered zigzags [3] owing to a higher modulus ratio and anisotropy). The disordered patterns include some T-shaped triple-junction folds that lack the rotational symmetry of the Y-shaped ones, see Figs. 3a,b. However, T-shaped folds do not form a simple periodic pattern and they may thus be considered as variants of the ideal Y-shaped folds. We note that the disordered patterns resulting from unaided folding, such as that in Fig. 3a, have slightly higher energies than the corresponding periodic zigzag or triple-junction patterns.
The optimal wavelengths of the zigzags and triplejunctions are shown in Fig. 3f and 3g, respectively. As in the uniaxial case, the wavelengths decrease with g when η > η s , but remain nearly steady when η ≈ η s , resulting in an approximate collapse of the wavelengths toward large g such that the width of the folds is again roughly four times the undeformed layer thickness. This means that further expansion of the layer is accomodated by increasing amplitude. The optimal x, y-aspect ratio of zigzags decreases with decreasing η such that the zigzags in soft layers are relatively shorter than their stiff-layer counterparts [18].
Our study demonstrates that complex fold patterns emerge from even the simplest systems incorporating elasticity and differential growth, providing a guide for more specific applications incorporating, for example, curved geometries, anisotropy, prestressed substrates, compressibility and stress relaxation. We show a transition as a function of modulus ratio from supercritical smooth finite-wavelength buckling followed by sulcification of valleys to direct strongly subcritical sulcification at 0.86 < η s < 1.14, depending on growth anistropy. This is accompanied by a transition from contraction to elongation of past threshold wavelengths, driving a convergence to a spacing of about 4h over a range of modulus ratios in both two and three dimensions. In equibiaxial growth there is also a beyond threshold transition from zigzags to triple-junction folds within the typical range of modulus ratios in soft layered tissues, underpinning the morphogenetic diversity of mechanical folding.
Supplementary Material for "The mechanics of invagination and folding: a soft layer growing on a soft substrate"

LINEAR STABILITY
We consider a three-dimensional elastic reference space labeled by x − y − z Cartesian coordinates and containing a soft layer of thickness a and shear modulus µ 1 occupying the region 0 < z < a bound to a substrate with shear modulus µ 2 ≡ µ 1 /η occupying the half-space z < 0. We assume that the layer has undergone growth (or swelling, or pre-strain before being glued to the substrate) such that, if it were cut free, it would undergo a relaxing shape-change described by the growth-matrix G. The material then undergoes a displacement field u(x, y, z), giving a deformation gradient is F = I + ∇u. If we model both materials as incompressible neo-Hookeans, the total elastic energy of the system, up to an irrelevant prefactor, will be We enforce incompressibility, encoded by Det(F ) = 1, by introducing a Lagrange multiplier pressure field P , leading to the effective energy density Minimizing this energy with respect to variations in u and P gives the Euler-Lagrange equations where σ = ∂L/∂∇u is the stress-tensor, given by Energy minimization also requires the free surface to be stress free, and the layer and substrate stresses to match at their interface: We restrict our attention to growth in the x direction, g x = g, g y = 1, for which we expect buckling in the x − z plane so we consider the linear stability of the flat state to infinitesimal sinusoidal perturbations with wavenumber k, writing the fields as Working in an x − z basis, to first order in we thus have At zeroth order in we have F = F −T = I and hence piecewise constant stress σ = µ G −1 · G −T − P 0 I . The bulk equations are thus already satisfied, and the boundary conditions require σ ·ẑ = 0 throughout, solved by Expanding eqn. 5 to first order in reduces it to ∇ · u = 0, requiring f 1 (z) = f 2 (z)/k. In the layer, eqn. 4 expanded to first order requires −g 4 f 2 (z) + k 2 f 2 (z) + g 2 P 1 (z) = 0.
The first of these can be solved algebraically for P 1 , giving and substituting this into the second gives the governing equation for f 2 (z): The same results hold in the substrate, but with g = 1.
The solution for f 2 in the layer and the substrate is thus where the a i and b i are constants of integration, and we have disgarded the unbounded solutions in the substrate. Continuity of z displacement at the interface requires continuity of f 2 at z = 0, giving while continuity of x displacement requires continuity of f 2 at z = 0, giving The first order correction to the free boundary conditions at z = a require a 1 g 4 + 1 + a 2 g 4 + 1 e and the first order correction to the interfacial boundary conditions at z = 0 require a 1 ηg 4 + η + 2 + a 2 ηg 4 + η − 2 + 2a 3 ηg 4 + g 2 + 2a 4 g 2 g 2 η − 1 = 0 (20) These four linear equations can be re-expressed as single matrix equation, 2e ak g 2 −ak g 4 2e which will only admit solutions if the determinant of the matrix is zero. This requirement reduces to the following condition for instability: 2 cosh ak g 2 g 4 − 1 2 3g 4 + 1 η sinh(ak) + 2g 2 g 4 −2 g 4 + 1 2 η 2 + 5g 4 + 2 + 1 cosh(ak) + sinh ak g 2 g 4 g 12 + 20g 8 + 6g 4 + 4 η 2 − 4 g 8 + 6g 4 + 1 + η 2 sinh(ak) − 2g 2 g 4 − 1 3 η cosh(ak) We find the first point of instability by finding the wave-number with the minimum g. Solutions to this eqn are shown for a variety of different moduli ratio in fig. 1a in the main text. Each line shows the strain required for instability at each wavelength for a given slab. The actual point and wavelength of instability are given by the first unstable mode, that is, by the minimum of the curve. Inspecting fig. 1a we see that all layers are unstable to a zero wavelength instability at g = g b = 1.839..., which is the Biot strain. We can extract this from the above condition by considering that the dominant terms at zero wavelength. Breaking the hyperbolic trig functions into exponentials, we see that the above condition contains two terms that diverge at small wavelength, so we may write the condition as where For asymptotically small wavelength the left exponential dominates and, since A = 0 the condition is only satisfied at the Biot strain g b = 1 1a we see that, for sufficiently soft layers (η ≤ 0.525...) this zero wavelength Biot solution is the first unstable mode, whereas for stiffer layers there is a finite wavelength buckling at a lower growth, giving rise to a more conventional buckling. We characterize the instability by minimizing the curves in fig.1a from the main text over wavelength to find the first unstable mode for each modulus ratio, which gives the wavelength of instability and the growth required.
These are plotted against modulus ratio in figs. 1b-c in the making paper, which again clearly show the transition from Biot zero-wavelength behavior to finite wavelength buckling at η ≤ 0.525.... For layers marginally stiffer than η b the threshold for instability is just below g b so we can analyze the instability by setting g = g b + and considering the leading behavior in . The Biot polynomial then becomes which is solved by This small correction to the Biot-threshold for short but non-zero wavelengths is positive, indicating the Biot point is a minima in fig. 1a, if B(g b , η) > 0 and negative, indicating the Biot point is a maxima, if B(g b , η) < 0. The crossover occurs at B(g b , η b ) = 0 which, inspecting the form of B, requires Further analysis indicates that the transition is continuous but logarithmic, in the sense that the unstable wavelength vanishes as λ ∼ 1/ log(η − η b ).

DETAILS OF NUMERICAL SIMULATIONS
Our numerical simulations are based on the explicit finite element method and discretizations based on constant strain triangle (2D) or tetrahedron (3D) elements. The elements have a quasi-incompressible nodal pressure formulation with a bulk modulus 10 3 (2D) or 10 2 (3D) times the shear modulus of the material. The simulated substrate has a thickness ten times the layer thickness. Whether the substrate has a clamped or free base has a negligible effect on the results, confirming that it provides a good approximation of an infinite substrate. The simulation codes are available online at users.jyu.fi/˜tutatall/codes/.
In 2D plane-strain simulations the domain includes a half fold with symmetric boundary conditions at sides. The domain is discretized to a rectangular mesh of triangles such that the growing layer includes 200 elements through its thickness. A fold is initiated by applying a small normal force to the surface at the edge of the domain. The surface  FIG. 4. Amplitude of a fold as a function of g, deterimined by continuously decreasing g while keeping the width fixed. a) A fold with η = 0.5 and width 2.50h unfolds continuously at the sulcification threshold g = gs ≈ 1.546, indicating supercritical sulcification. b) A fold with η = 1 and width 4.71h (corresponding to the energy minimum at the sulcification threshold) unfolds discontinuously below gs, with a jump in amplitude, indicating subcritical sulcification. c) A fold with η = 4 and width 7.54h (corresponding to the wavelength of linear instability) unfolds continuously at the threshold of linear instability, indicating supercritical smooth buckling. The systems in (a) and (c) would follow the same path if folded by continuously increasing g, but the system in (b) would remain flat until gs, followed by instant folding to a finite amplitude (assuming there is a sufficient perturbation to initiate sulcification at g = gs). is prevented to cross the vertical line that defines the boundary of the domain to prevent self-intersection of the fold. The energetically optimal wavelength is searched by adjusting the relative width of the domain.
In 3D simulations the domain has periodic boundary conditions at the lateral boundaries. The domain is discretized to tetrahedrons such that the nodes form a hexagonal prism mesh. The growing layer includes 10 elements through its thickness. Self-contacts of the surface are accounted for by preventing nodes penetrating faces at the surface. We perform two kinds of simulations in 3D: 1) Simulations on a large rectangular domain (width 30 times the layer thickness) where the layer is allowed to fold and relax freely (Fig. 5). The surface of the layer in this case has small random perturbations that allow folding to initiate. 2) Simulations on a small domain that contains one unit cell of the rectangular pattern of zigzags or hexagonal pattern of triple-junctions. In these simulations we use point forces in selected locations of the domain to initiate formation of the patterns. When the dimensions of the domain are close to the optimal wavelength of the desired pattern, a stable periodic pattern forms easily. Dimensions of the domain are fine-tuned to find the energy minimizing wavelength.
An important difference between a simulated and true continuum system arises when the surface reaches the sulcification threshold as a result of increasing g: the continuum surface would sulcify as a result of an infinitesimal perturbation, but the simulated surface requires a finite perturbation to cross the barrier set by the discrete mesh. The curve describing the transition between smooth and sulcified folds in Fig. 2A of the main article was thus determined by finding the sulcification strain at the surface, with simulateneous optimization of wavelength. The foded surfaces shown in Fig. S2 are also obtained by unfolding the layers that were first relaxed at g = g x = g y = 2.