Buckling without Bending: A New Paradigm in Morphogenesis

A curious feature of organ and organoid morphogenesis is that in certain cases, spatial oscillations in the thickness of the growing “film” are out of phase with the deformation of the slower-growing “substrate,” while in other cases, the oscillations are in phase. The former cannot be explained by elastic bilayer instability, and contradict the notion that there is a universal mechanism by which brains, intestines, teeth, and other organs develop surface wrinkles and folds. Inspired by the microstructure of the embryonic cerebellum, we develop a new model of 2D morphogenesis in which system-spanning elastic fibers endow the organ with a preferred radius, while a separate fiber network resides in the otherwise fluidlike film at the outer edge of the organ and resists thickness gradients thereof. The tendency of the film to uniformly thicken or thin is described via a “growth potential.” Several features of cerebellum, +blebbistatin organoid, and retinal fovea morphogenesis, including out-of-phase behavior and a film thickness amplitude that is comparable to the radius amplitude, are readily explained by our simple analytical model, as may be an observed scale invariance in the number of folds in the cerebellum. We also study a nonlinear variant of the model, propose further biological and bioinspired applications, and address how our model is and is not unique to the developing nervous system.


I. INTRODUCTION
An elastic instability, driven by differential growth, is thought to be broadly responsible for many of the motifs seen in organ morphogenesis [1]. Indeed, this mechanism has been studied in the context of brain folds [2][3][4][5][6][7][8][9][10], intestinal crypts and villi [11,12], airway mucus wrinkles [13,14], tooth ridges [15], and hair-follicle patterns [16,17], among others. Wrinkling or buckling provides a means for these organs' shapes to emerge reliably from their respective starting geometries, without appealing to spatial variation in gene expression or other biochemical prepatterning of folds. Apart from geometry, all that is required is a competition between the bending energy of a uniformly growing film (cortex, mucus, or epithelium) and the energy to stretch and compress a slower-growing substrate (subcortex, submucus, or mesenchyme).
Biological tests of wrinkling predictions have largely focused on pattern formation and wavelength, the latter scaling with film thickness and power 1/3, 1/4, or 1/6 of the stiffness contrast, depending on the substrate model [18]. (Stiffness contrast is defined as Young's modulus of the film divided by that of the substrate.) The wavelength test is typically not a conclusive test for two reasons. First, measuring the stiffness contrast in vivo, at the appropriate developmental period, and on the timescale relevant to morphogenesis, is technically challenging. To our knowledge it has not been done, and the best available data come from ex vivo measurements on freshly harvested embryonic organs [12,19]. Second, the weak dependence of wavelength on stiffness contrast exacerbates the uncertainty of the latter: a large range of stiffness contrasts might be said to "agree" with an observed wavelength. Considering these difficulties, it is surprising that other wave properties predicted by wrinkling theory have not been given more attention. In particular, the phase and amplitude behaviors associated with elastic wrinkling are qualitatively distinct and constitute a stringent pass-or-fail test with regard to experimental biology observations. Let us see how this works.
In simple wrinkling analyses, e.g., Euler buckling of a film adhered to a bed of springs which yields the power 1/4 mentioned earlier, the film thickness is typically assumed to be spatially uniform. In reality, the film thickness is modulated by the substrate deformation. A crude consideration of the forces acting on the quasistatic interface indicates the film thickness oscillations should be in phase with the substrate deformation. That is, thick spots in the film (regions of the film under internal tension) should be matched up with thick spots in the substrate (regions of the substrate under internal tension), and thin spots with thin spots. Finite element simulations confirm the essential correctness of this argument for wrinkling in planar geometry, and also circular geometry, provided the film is thin compared to the substrate radius and the stiffness contrast is not too large (see Fig. 1). In general, however, simulations reveal two film thickness minima per wrinkle: one coincides with the wrinkle valley, and the other coincides with the wrinkle crest, as can be seen in the lower panels of Fig. 1. The latter has small depth and width compared to the former (i.e., the thickness profile is overall in phase) except when the modulus ratio is large and/or the film thickness is a significant fraction of the substrate radius. Generically, the amplitude of thickness oscillations is much smaller than the wrinkling amplitude. The large thickness limit that one might expect to yield comparable amplitudes tends not to undergo wrinkling, but rather a global buckling mode [20]. Full details of these simulations, as well as additional examples (one of which shows a wrinkling + global buckling mixed mode), are given in the Appendix.
The simple quasistatics argument for in-phase behavior straightforwardly extends to growth and elastic-modulus profiles (as a function of radius) that are more complicated than a simple bilayer profile, including even the possibility of nonmonotonicity. First notice that any smooth, continuous growth and/or modulus profile can be represented by many discrete layers. Specializing to those layers having thicknesses that are small compared to their radii, we observe that any two adjacent such elastic layers must exert equal and opposite normal forces at every point along their shared interface. As before, normal tension and normal compression correspond to thicker than average, and thinner than average, respectively, assuming initial axisymmetry. Thus, where any one layer is thicker than its average thickness, all the other layers must be relatively thick, and any arbitrary contiguous grouping must be relatively thick. This generalization serves to emphasize that in-phase thickness behavior is a generic consequence of differentially growing elastic multilayers (assuming near-planar geometry), and is independent of the precise details of the growth and/or modulus profiles.
So while elastic instabilities driven by differential growth exhibit a number of interesting features, it follows that if a differentially growing biological system were found with quasistatic out-of-phase behavior (film thickness maxima coinciding with substrate valleys, and thickness minima with hills), elastic wrinkling could potentially be ruled out as the mechanism of shape change. Moreover, a thickness amplitude that is not small compared to the surface height amplitude would be at odds with elastic wrinkling. And yet it is these very "antiwrinkling" behaviors that show up in several motifs in organ morphogenesis: the cerebellum, certain brainlike organoids grown in vitro, and the retinal fovea [see Figs. 2(a)-2(c)]. Note that all three examples are pertinent to nervous-system development. The first two of these have been previously considered as elastic wrinkling problems [9,10], as has brain folding more generally [2][3][4][5][6][7][8]. Fovea formation has not previously been credited to elastic instability, as far as we know, but there does appear to be differential growth between the constituent layers of the retina in the vicinity of the developing fovea, e.g., Fig. 2 in Ref. [21]. In effect, the first two of these biological systems, and tentatively the third about which less is known, are counterexamples to the elastic instability paradigm for morphogenesis. They call for a new paradigm.
In constructing this new paradigm for shape change in developing organs, we go beyond modeling the constituent layers as elastic materials with different elastic constants. If one looks at microstructures within the developing nervous system, one typically finds globular cells, such as granular cell precursors in the cerebellum or neural precursors in the cerebrum, and fiberlike cells, such as radial glial. At some length scale both types of cells, if they are interconnected enough, can presumably be modeled as an elastic continuum. But at what length scale and degree of interconnectedness is such modeling justified? The microstructure of the developing organ hands us a clue as to when such an approximation breaks down and so we should look to it in developing this new paradigm. Moreover, it has recently been shown that in a model for confluent cell tissue, the introduction of cell division drives the tissue from an elastic solid to a fluid [25]. As with any developing mammalian organ, cell division is tantamount to the process, so perhaps we must further relax the notion that the different components involved are all elastic and begin to consider that at least some of the components are fluidlike. window just mentioned, the mouse cerebellum transitions from having a featureless convex surface, to developing smooth, sinusoidal radius oscillations, to then forming cusped invaginations, termed "anchoring centers." The anchoring centers delineate smooth outward protuberances called lobules [22]. Later in development, some of the first-generation lobules subdivide into second-generation lobules [28]. (Human cerebella appear to have several generations of subdivisions.) In the following, we will exploit the quasi-2D nature of the cerebellum (its rows of parallel folds indicate that most of the biomechanical action is in the parasaggital plane), which makes it a simpler system to analyze than the cerebrum with its 3D folds. We note that earlier work on modeling the cerebral cortex as a smectic liquid crystal (as opposed to a purely elastic solid), being pulled on by axonal tension, is a first step in the direction of incorporating more of the microstructure into cerebral development [5].
As for additional applications of the model, a recent in vitro experiment demonstrates shape change in a human embryonic stem-cell-derived aggregate inserted into a quasi-2D microfabricated compartment filled with hydrogel [10]. After several days, elongated fiberlike cells are present along the periphery of the aggregate, while globularlike cells remain in the interior. The former are representative of a cortex and the latter, an inner core. During its second week in development, the surface of the organoid invaginates in a manner reminiscent of brain folds. The introduction of blebbistatin, a myosin inhibitor, produces qualitatively different shape changes as compared to the untreated case. We will apply our prototype model to the +blebbistatin organoid, as opposed to the untreated case, since the developing shapes in +blebbistatin organoid more closely resemble the shapes obtained in the model. We will discuss the untreated case as a potentially more sophisticated model. Finally, we consider retinal development as a third possible application of (a variant of) our prototype model. The developing retina is encased in a rigid shell, with several distinct layers of cells supported by the inner part of the shell, followed by vitreous, a clear gel, filling the space between these layers and the lens of the eye [23]. Müller glial fibers span some of the layers in the macular region, where the center of the field of vision is focused. A depression known as the foveal pit begins to form in this region around human gestational week 25, and develops over a timescale comparable to that of brain folding [21].

II. MODEL
Let us model the growing cerebellar cortex as a 2D annuluslike region having outer radius r and thickness t, which are scalar functions of an angular coordinate 6. In other words, t is defined to be measured in the radial direction (see Fig. 3). This simple parametrization is valid only for weak deviations from an annulus. Consider, e.g., that a deep or overhanging surface fold could generate multivalued r(θ), as well as lead to a t(θ) that violates one's sensibilities around the usual notion of thickness. Under this restriction, we introduce the quasistatic energy functional E r, t, dt dθ to be minimized subject to a constraint on the area of the nongrowing subcortex, i.e., The variational problem at hand is thus where μ is a Lagrange multiplier whose value will be determined upon simultaneous solution of Eqs. (2) and (3).
In Eq. (1), k r , k t , and β are all positive constants. The first term encodes a preferred radius r 0 , or more generally, a preferred shape r 0 (θ). Because of its negative contribution, the second term favors thickening (or thinning) of the annulus with respect to a reference thickness t 0 , and for simplicity, we will take t 0 to be a constant. Thus, while k r is a modulus, k t can be regarded as a "growth potential," and the corresponding terms compete with one another because of the subcortex incompressibility. This competition tends to drive the system away from its preferred shape. The third term in Eq. (1) penalizes spatial variations in thicknessit is not a bending term as its appearance may suggest to some. (A bending term would involve a squared second derivative of the film deflection.) In fact, the absence of a cortex bending modulus is a key feature of the present model that distinguishes it from elastic bilayer models of brain folding. Here, the cortex resembles a mixture of fluid + fibrous scaffolding held inside a container. Both container and scaffolding are flexible in bending, but the scaffolding resists gradients in the container's thickness. As for the physical nature of this scaffolding in an actual cerebellum, we suggest its main component is the cortexspanning Bergmann glia, depicted schematically in Fig. 3. The nongrowing subcortex, in contrast, is treated as an elastic solid. Together with radially oriented fibers that span the system (radial glia), this is the origin of the preferred radius term in Eq. (1). We suggest that radial glia "tether" the surface to the subcortex, using their washerlike pial end feet to distribute the load over the flexible surface, while their other ends remain securely anchored in the solid subcortex. Implicit in this anchoring mechanism is the subcortex's finite shear modulus, however, the actual shear deformation energy of the subcortex is neglected.
So that we may make progress analytically, let us take the preferred shape to be an ellipse having semimajor axis a and (small) eccentricity e, i.e., With this choice, decoupling the Euler-Lagrange equations results in an unconventionally driven oscillator equation for the thickness }(e 2 /2), ∊ = μ/k r , and c = k r /k t are all constants (within the quasistatic description). In the regions of the ∊ -c plane where q is real (see Fig. 4), Eq. (5) has the general solution in terms of which the radius is Combining Eqs. (2), (6), and (7), we find the thickness amplitude is given by and evidently we must impose a lower bound on A 0 to ensure this amplitude is real valued: The physical significance of the Lagrange multiplier's sign is illuminated by a relationship between angle-averaged quantities which is exact to all orders in e. Since physically viable solutions require rt > 0 for all θ, Eq. (10) says that ∊ > 0 corresponds to growth of the cortex while ∊ < 0 corresponds to shrinkage (not in a dynamical sense, necessarily, but with respect to the quasistatic value of t 0 ). Note 〈r -t〉 > 0 also implies that oscillatory solutions having ∊ > 1 are presumably unphysical in the sense at 0 < 0. Consequently, we will restrict our focus to ∊ < 1.
Apart from its sign, the Lagrange multiplier may be thought of as a pressurelike quantity, or a "chemical potential" for changing the core area by a unit amount (the symbol μ was deliberately chosen to suggest this analogy). It can also be argued that ∊ should set the amplitude T. Consider that as ∊ → 0, the competition between preferred radius and film growth disappears, as nothing prevents the film from growing uniformly inward and consuming the core. Absent this competition, there is no driving force for film thickness oscillations, one might argue. The simplest way to make the model consistent with this argument would be to set T ~ ∊, i.e., make the film thickness amplitude linearly go to zero as ∊ goes to zero (the radius amplitude would go to zero quadratically). This assignment would also establish ∊ as the key parameter governing all phase and amplitude aspects of shape change. Endowing ∊ with time dependence would therefore be a natural (and minimal) starting point for dynamical morphogenesis problems, within the current modeling framework. For the remainder of this paper, however, we shall mostly confine ourselves to a statics approach, as well as regard T and ∊ as independent fitting parameters. Time dependence of the model parameters is further investigated elsewhere [27].
The simplicity of the above-described model belies its novelty, in that several of the predicted behaviors are opposite those predicted by differentially growing, elastic bilayer models. For example, t and rt oscillations are always out of phase when ∊ < 1 (compare τ and r s oscillations, respectively, in Fig. 1), while t and r oscillations can be either in phase, or out of phase, depending upon the position in the ∊ -c plane (see Fig. 4). Likewise, the amplitude of t oscillations can be either greater than, or less than, the amplitude of r oscillations (for ∊ < 1/2 and ∊ > 1/2, respectively). As another example, notice that as e ∊ → 0, the wave number q depends only on a ratio of energy scales associated with microscopic mechanisms, i.e., k t /β. In contrast, elastic wrinkling predicts that the number of waves depends on a ratio of length scales, roughly equivalent to our a/t 0 , which may or may not be scale invariant.
In Figs. 2(d) and 2(e), we compare plots of r and rt with images of an embryonic mouse cerebellum and a brainlike organoid treated with blebbistatin. (The region in between the two plotted curves, representing the cortex, is filled with dark purple and dark green, respectively.) Some comments specific to the organoid application are now in order. As previously noted, the inner core of this structure consists of globularlike cells, while the periphery consists of fibrous cells as well as motile cells that divide and move along the fibrous cells in such a way that the periphery grows faster than the core. The fibrous cells and motile cells are reminiscent of Bergmann glia and granular precursor cells in the cerebellum, respectively. There are assumed to be radial glial cells or some other means of transmitting radial tension throughout the organoid, such as that discussed in Ref. [27]. From this description, one could argue that the model at hand may be applicable to the untreated organoid. However, the model does not take into account the active contractility of the core, due to the myosin that is present there. This activity is likely to be mechanosensitive such that assumptions beyond those we have already made would be required. While this is an interesting avenue to pursue, we will restrict our application of the model to the blebbistatin treated organoid with a less active core.
Clearly, two of the fit parameters in Figs. 2(d) and 2(e) are associated with an elliptical preferred shape (e and ϕ). That leaves five dimensionless parameters for a circular preferred shape: ∊, c, k t /β, T/a, and t 0 /a. [It is convenient to regard T/a, rather than A 0 /a, as a fitting parameter, as one may then guarantee Eq. (9) is satisfied.] The first of these five is constrained by an experimental image from which one can measure (or at least estimate) the ratio of the t and r amplitudes. One can also measure T and 〈t〉 in units of 〈r〉, as well as count the number of invaginations q. In general, then, there are four independent constraints on these five model parameters. Suppose, however, the first measurement yields ∊ ≪ 1, such as appears to be the case at the onset of shape change (around embryonic day 16.5 in mice [22]). In this limit, neither the leading-order AC terms nor leading-order DC offsets involve c, because it only appears in Eqs. (6) and (7) as a product with ∊. Absent c, the remaining model parameters are all directly measurable: k t /β ≈ q 2 , T/a ≈ T/〈r〉, t 0 /a ≈ 〈t〉/〈r〉.
The scaling behavior q k t / β can be understood in a simple way, as follows. As ∊ → 0, the mechanical constraints on the outer surface of the fluid layer become relatively severe compared to those on the inner surface, because either the system-spanning radial springs are being turned into rigid rods (k r → ∞), or the core area constraint is being removed (μ → 0). Thus, the degrees of freedom representing different configurations of the outer surface are effectively frozen out of the inner surface problem. Only two terms now contribute significantly to the energy: the growth term, which scales as E grow ~ −k t T 2 , and the gradient term, which scales as E grad ~ βT 2 r 2 /λ 2 . Here, λ is the wavelength of film thickness oscillations. Comparing these two energy scales yields λ r β/k t . Again, this scaling of the wavelength with system size is in contrast to elastic bilayer wrinkling, where the wavelength scales with film thickness. Higher-order growth instabilities lead to different scaling behavior just as different substrate models generate different scaling behavior in elastic bilayer wrinkling [18]. A lower order, hence linear, growth instability leads to unphysical results in the small e limit, which is why we did not implement it here.

III. CONCAVE VARIANT
Suppose, instead of the convex bilayer depicted in Fig. 3, the system of interest is a concave bilayer contained within a rigid circular boundary having radius R. The substrate with conserved area A 0 is in contact with the boundary wall, and the growing (or shrinking) film is interior to that. Both layers are annular in shape, with thicknesses rt and t, respectively.
Such a geometry has been used in prior work on morphogenesis of the gut [11,12], airways [13,14], and other tubular structures [29]. The sole difference in our formulation of this variant is that the variational principle involves an extra term with respect to the convex problem: δ{E -μ ∫ dθ[(r -t) 2 -2R(r -t)]} = 0. One finds, in this case, that R appears only in the DC offset terms of the solution and Thus, the phase diagram is the same as before, save for which regions correspond to growth and shrinkage. This latter difference is because which reduces to Eq. (10) only on the phase boundary c = 1 -∊ −1 . Note also that as R → ∞, growth and shrinkage correspond to c < 1 and c > 1, respectively. The thickness amplitude is given by where ρ = [a -t 0 + ∊(c-1)R]/(1 -∊ + ∊c), which implies an upper bound on A 0 .
In Fig. 2(f), we plot a solution of this concave variant next to the image of the retinal fovea, as its (presumably differentially growing) layered structure is enclosed by a nongrowing rigid shell on one side and vitreous gel on the other, while fibrous cells effectively span its "cortex." The parameters used (see the Fig. 2 caption) are consistent with 〈t -t 0 〉 > 0, i.e., growth of the ganglion cell layer, despite the fact that ∊ is negative, which is a qualitative difference from the convex model.
One may point out that our modeling of foveal pit morphogenesis is unrealistic because many pits are generated instead of one or two. But just as localized growth leads to localized buckling in elastic bilayer models, we can potentially resolve the discrepancy at hand by introducing a spatially inhomogeneous growth potential k t = k t (θ). Two analytic cases are deserving of mention. For the first case, suppose k t (θ) is a piecewise periodic function, having a constant value in an arbitrarily small interval [θ 1 , θ 2 ] and a different constant value everywhere else. In the variational problem, we are free to take the bounds of integration as θ 1 and θ 2 (because there are no long-range interactions in the circumferential direction) and solve two independent problems: one for each growth potential region. Both solutions are, of course, given by Eqs. (11) and (12), and it remains to apply matching boundary conditions. Unfortunately, taking k t → 0 regionally leads to unphysical behavior in the analytic limits of the DC offsets (the same is true of the convex variant), so one must settle for a small but nonzero value of k t . Alternatively, one might be interested in a slowly varying growth potential, and in this case, we expect a local density approximation to be valid. For small ∊, this would read q(θ) ≈ k t (θ)/ β. These possibilities for localized buckling also apply to the convex variant.

Author Manuscript
Author Manuscript Author Manuscript

Author Manuscript
Finally, we suggest that it would also be interesting to apply this concave variant to the gut, e.g., if the out-ofphase motif were discovered at the appropriate time of development.

IV. NORMAL THICKNESS VARIANT
When there are strong deviations from circularity, the "radial thickness" t is no longer a good measure of film thickness. A more natural measure is the "normal thick-ness" τ = t(r ⋅ n), where r = cos(θ)x + sin(θ)y is the unit where N and M are shorthand for the components of r in the direction of the surface normal and surface tangent, respectively, and M/N = (tan ϕ + cot θ)/(1 -tan ϕ cot θ). After simultaneously solving Eqs. (18) and (19), the radial thickness may be recovered per its definition, t = τ/N.
Though complicated and nonlinear, Eqs. (18) and (19) are readily converted into a system of first-order difference equations, and solved numerically. We use multidimensional Newton-Raphson iteration on a periodic grid, with Krylov approximation of the Jacobian. The analytic solution of the corresponding linear problem provides a convenient initial guess for relaxation. While these numerical solutions show some sensitivity to the number of grid points used, we find that 10 000 grid points gives rapid convergence and reproduces the linear model's behavior at small e, i.e., where the radial thickness and normal thickness definitions coincide. Furthermore, at radial vector and n = − sin(ϕ)x + cos(ϕ)y is the unit surface normal vector, with We now modify the variational problem [Eqs. (1) and (3))], by replacing t with τ in the last two terms of the energy functional E* r, τ, τ′ = ∫ dθ k r r − r 0 2 − k t τ − t 0 2 + β dτ dθ 2 , (16) but leaving the variational principle otherwise unaltered, Choosing r and τ as the independent variables, as we have done here, reduces the order of the Euler-Lagrange (EL) equations. Had we chosen instead to keep r and t as the independent variables, a higher derivative r″ would be introduced into the energy functional by the (dτ/dθ) 2 term, and fourth-order EL equations would result (see, e.g., Ref. [30]). The reduced-order EL equations are given by And modest values of e and c, the observed deviation of the shapes from those of the linear model can be rationalized as follows. A sharp bend in a film having uniform normal thickness incurs no penalty from a (dτ/dθ) 2 term, but it does incur a penalty from a (dt/dθ) 2 term (although not in the same sense as that for bending an elastic film). Therefore, we might expect the normal thickness variant to exploit the extra degrees of freedom afforded by these low-energy, sharp bends in negotiating the competition between growth or shrinkage and the preferred radius. Figure 5 shows these sharp bends occur preferentially where the film is thin. One interesting consequence of this behavior is that additional minima in the substrate radius can be introduced, e.g., in the right panel of Fig. 5. In contrast, where the film is thick, its shape is not significantly distorted from that in the corresponding linear model.

V. DISCUSSION
We have demonstrated that the out-of-phase behavior observed in certain morphogenesis contexts is not only at odds with an elastic bilayer wrinkling mechanism, but indicative that microstructural details affect morphogenesis in heretofore unappreciated ways. Simply treating the problem as one of differential growth between two homogeneous elastic materials appears insufficient to capture the unique shapes of the developing cerebellum, +blebbistatin organoid, and the retinal fovea. We have constructed a minimal model that captures these shapes (at least qualitatively), and in our model radially oriented fibers play a key role.
It is reassuring, then, to notice that all three of the above-mentioned biological systems contain one or more types of radial fibers. One is, therefore, led to speculate whether our model is solely applicable to the central nervous system, where requisite fibers such as glia are present during development of all regions, and certain of whose organs exhibit the telltale out-of-phase behavior. At the time of writing, we are not aware of definitive out-ofphase behavior in other morphogenesis contexts, but we are actively searching. The developing cerebrum is a natural place to look, but interestingly, in adult humans, the cerebral cortex thickness appears to be in phase with its surface height (see Fig. 1 in the Supplemental Material of Ref. [31]). Another place to look is where there is smooth muscle.
If the smooth muscle tissue is very thin and minimally connected, presumably the fibrous nature of the cells making up the tissue call for a microscale mechanical description like the one we have developed here. A novel modeling description for such a system could lead to additional novel shape-changing mechanisms for developing organs.
Application of the model to other developing organs may also potentially call for an extension into the third dimension. One can then define a spherical inner core and a spherical outer shell of proliferating cells as well as fibrous cells radially extending themselves throughout the system. The developing undulations along the perimeter of the shape in two dimensions now becomes developing undulations over a surface and one can ask whether or not the undulations on the outer surface are in phase or out of phase with the inner surface, depending on the parameters. Since the key new aspects of the model are all contained within the two-dimensional version, we focus on its results for now.
The possibility of constructing a synthetic device that embodies the Hamiltonian given by Eqs. (1) and (2) is an intriguing one. Such a device would consist of (1) an incompressible core, (2) a growing film that has the curious combination of flexibility in bending and stiffness against thickness gradients, and (3) system-spanning, radially oriented, elastic fibers. (See the discussion in Ref. [27] for an alternative, potentially simpler way to realize the ~k r term.) The film component would likely have to be some sort of fiber-matrix composite material, as it also appears to be in the cerebellum. The nontrivial mechanical properties of this film component could perhaps be reverse engineered by trying different composite formulations and watching for out-of-phase behavior.
Potential applications of the model outside of molecular and tissue-scale biology may also be worth investigating. Morphogenesis of large cities, e.g., shares several of the key ingredients of our model (at least in a qualitative sense). Cities are 2D structures, and there is in certain cases an "incompressible" urban core, surrounded by a growing suburban belt. The preferred radius concept is also not unreasonable: people relocating to or within the city might be expected to strike a balance between commuting time to the city center and housing prices, which tend to fall off with distance from the center. Chengdu, the capital of China's Sichuan province, is one example of a large and rapidly growing city that is free of significant geological constraints on its shape changes (as evidenced by its nearly circular shape). It would be interesting to try and interpret Chengdu's shape changes within the past few decades [32] from within a differential growth framework.
Model variants beyond the two described here are, of course, possible, as well as numerous. Already mentioned was the notion of putting spatial dependence into one or more parameters. One could also, e.g., introduce a curvature-dependent growth potential, k t = k t (r ″), which might be expected to control the degree of cuspiness of the invaginations. An anharmonic correction to the preferred radius term, or a tight-packing constraint, might generate self-contacting folds by squashing or squeezing the lobules (i.e., wrinkle crests) together. Neither the linear model nor the normal thickness variant that we have investigated here appear capable of producing cusped invaginations or self-contacting folds, such as those that occur in mouse cerebellum development after about 18 embryonic days [22]. Extending the model beyond just the onset of shape changes, to the large growth and The preferred radius r 0 naturally suggests a mechanism for subdivisions, and as previously mentioned, subdivisions occur in the mammalian cerebellum during later development. Time dependence could be introduced into our quasistatic model or one of its variants in such a way that lobules are dynamically growing in area. When a lobule's area reaches a value r 0 2 , the lobule may in fact constitute a subsystem that resembles the full system at an earlier point in time. Shape change of the subsystem would be akin to a folding hierarchy, and one might imagine this going on for several generations. An implementation in which subdivisions emerge spontaneously, as opposed to being put in "by hand," is a very interesting direction for future work.
Finally, we point out that an apparent scale invariance in the wave number q, found here in the limit ∊ → 0, is reminiscent of "Larsell's criterion for the vermis" [28]. Adult mammalian cerebella span 2 orders of magnitude in size, ranging from ~1 mm in mouse to ~10 cm in humans, and Larsell observed and argued that these can all be considered as having an underlying ten fold motif. In his book, Larsell notes that other researchers before him including Bradley, Bolk, and Riley conducted comparative anatomy studies of the cerebellum (an embryological study spanning six species in the case of Bradley) and drew similar conclusions about its scale invariance [33]. While Larsell's criterion is not universally accepted by biologists, the fact remains that the cerebellum of all mammal species is highly folded, whereas the cerebrum is unfolded (i.e., lissencephalic) in the smallest mammals, thus there is at least a hint of scale invariance in cerebellar morphogenesis, which our model may capture.

VI. BIOLOGICAL PERSPECTIVE
A biologist might naturally ask, if buckling is a mechanical phenomenon that requires no genetic prepatterning, then why does the specific buckling mechanism (with versus without bending) matter for biological contexts? What can I infer about the biology of a tissue from the phase and amplitude behavior of its wrinkles?
In our view, these kinds of questions target the notion of emergent behavior in bulk material systems. To use a classic example, some bulk materials superconduct at low temperatures while others do not-these differences in emergent behavior can teach us something about the basic atomic building blocks of the materials and the way those blocks are arranged via, e.g., the electron-phonon inter-action. This is "backwards" from the reductionist view-point, that would have us start from the basic atomic building blocks to try to understand superconductivity.
Likewise, in organ morphogenesis, the macroscopic phase, amplitude, and wavelength behaviors are associated with the emergent phenomenon of buckling instability, and they can teach us something about the microbiological building blocks. In-phase behavior (and other behavior consistent with the rightmost column of present, they should be mixed together on a fine length scale in every direction. Furthermore, there must be little to no stress-relieving rearrangements and neighbor exchanges of cells (i.e., fluidlike behavior) on the approximately 12-h timescale of shape changes; all significant stress reduction happens over the much larger length scale of the wrinkles. Out-of-phase behavior teaches us that certain types of fiberlike cells are present and span relatively long length scales, and further suggests what will happen when we add or remove these fibers genetically (see the consequences of tuning k r and β summarized in Table I). It also provides some clues as to what is happening to the fibers as the organ grows over time (see our experimental companion paper, Ref. [27], which further tests both models against mouse cerebellum development, including predictions beyond those listed in Table  I). Finally, and perhaps most importantly, out-of-phase behavior is indicative of fluidlike rearrangements and neighbor exchanges of at least one cell type in the film component, presumably as a result of both cell divisions and cell motility. Such lessons from emergence may ultimately help us understand and correct what is going wrong, microbiologically, in brain-folding-related diseases such as lissencephaly and polymicrogyria, and in certain malformations of the retina such as the bifoveality shown in Fig. 2(c). and outer surface of the film, as shown in Fig. 6. We observe some noise in these thickness data, which can be attributed to the finite size of the FE mesh. However, this noise is a second-order perturbation and will not change the overall thickness variation of the film. The film thickness variation is compared to the substrate deformation in Fig. 7.

FIG. 3.
Idealized, parasaggital section of an embryonic cerebellum. The cortex is modeled as a growing fluidlike "film." Bergmann glia fibers span this film, while radial glia fibers span the structure. Periodic boundary conditions are likely not relevant to a real cerebellum, whose cortex is discontinuous owing to a ventricular zone and attachment to the brain stem, but they will be adopted here for simplicity.  Additional bilayer wrinkling simulations, similar to those shown in Fig. 1, but with different geometries and modulus ratios. Summary of main predictions for the growth regime (0 < ∊ < 1) of the convex variant, with comparison to the corresponding predictions of a standard morphogenesis model. The ∊ ≪ 1 limit pertains to early development, i.e., just after the onset of shape change. The five dimensionless model parameters discussed above are written in a slightly different (but equivalent) form here in order to emphasize their physical meaning.

Conventional elastic bilayer model
Phase relationship between film thickness and substrate deformation (the planar limit t/r