Buckling without bending morphogenesis: Nonlinearities, spatial confinement, and branching hierarchies

During morphogenesis, a featureless convex cerebellum develops folds. As it does so, the cortex thickness is thinnest at the crest (gyri) and thickest at the trough (sulci) of the folds. This observation cannot be simply explained by elastic theories of buckling. A recent minimal model explained this phenomenon by modeling the developing cortex as a growing fluid under the constraints of radially spanning elastic fibers, a plia membrane and a nongrowing sub-cortex (Engstrom, et. al., PRX 2019). In this minimal buckling without bending morphogenesis (BWBM) model, the elastic fibers were assumed to act linearly with strain. Here, we explore how nonlinear elasticity influences shape development within BWBM. The nonlinear elasticity generates a quadratic nonlinearity in the differential equation governing the system's shape and leads to sharper troughs and wider crests, which is an identifying characteristic of cerebellar folds at later stages in development. As developing organs are typically not in isolation, we also explore the effects of steric confinement, and observe flattening of the crests. Finally, as a paradigmatic example, we propose a hierarchical version of BWBM from which a novel mechanism of branching morphogenesis naturally emerges to qualitatively predict later stages of the morphology of the developing cerebellum.


I. INTRODUCTION
When morphogenesis is viewed through the lens of a physicist, one of the typical mechanistic routes to take is morphoelasticity induced by varying internal stresses [1]. Shape change in response to internal stresses, within elastic objects which tend to retain shape, provides us with a large collection of shapes. Differential growth can be a source of such an internal stress. This reasonable viewpoint is made manifest in the widespread use of the Euler buckling instability in purely elastic materials to explain the onset of folds in morphogenesis problems [2] as diverse as cerebra [3][4][5][6][7][8], intestinal crypts and villi [9,10], airway mucus wrinkles [11,12], tooth ridges [13] and hair-follicle patterns [14]. Recent work accounting for cerebellar foliation falls under the same framework [15,16]. These models typically predict a characteristic length scale between folds and state quantitative agreement between prediction and observation to validate such an approach.
On the other hand, tissue fluidity in early stages of developing embryos has emerged as a driver of shape change in both zebrafish and the insect Tribolium castaneum [17,18]. Fluids, unlike elastic solids, do not retain their shape. Then how does fluidity -the antithesis of elasticity -affect the shape of the developing organ? Presumably biological systems have figured out ways to combine elements of elasticity and fluidity to bring about an even more complex collection of shapes at later stages of development. The buckling without bending morphogenesis (BWBM) model is one such model affirming the * mgandiko@syr.edu † jschwarz@physics.syr.edu cleverness of biological systems [19]. Its origin is rooted in explaining recent quantitative observations of the developing murine cerebellum, a much less studied component of the brain as compared to the cerebrum. Vertebrate brains are typically divided into three different sectors -the forebrain, the midbrain, and the hindbrain. While the forebrain consists of the cerebrum, the hypothalamus, and the thalamus, the hindbrain contains the spinal cord, medulla oblongata, and the cerebellum. Since the shape of the forebrain varies more across species than the midbrain and hindbrain, its development has been the focus of much study. On the other hand, by studying conserved patterns in the hindbrain, we gain complementary insight into the morphogenetic processes of the brain.
To be specific, let us characterize the difference in shape between the cerebrum and cerebellum. First, the shape of the cerebellum has an approximate cylindrical symmetry, whereas such a symmetry is absent in the cerebrum. Second, mammalian cerebellums of all sizes produce folds [20], whereas small mammalian cerebrums do not show folds [21]. Given that the overall size of both organs increases with increasing species size, it is possible that the differences in shape arises via different physical shape change mechanisms. If so, do different physical shape change mechanisms lead to different emergent functionalities of the two? Moreover, the conservation of the number of the 8-10 primary lobes of the cerebellum across species demands a treatment on its own footing as it suggests a scale invariant shape change mechanism. Number of secondary folds however, can differ across species.
Recent experimental observations of the developing murine cerebellum found that the proliferating cells in the cerebellar cortex are motile with neighbor exchanges on the order of minutes [22]. It was also observed that the developing cerebellar cortex varies in thickness with the trend being that the cortex is thinnest at the crest (gyri) and thickest at the trough (sulci). Moreover, the developing cerebellum is under tension and not under compression as evidenced by both radial and circumferential cuts. While these observations cannot be explained by elastic wrinkling theories, all three of these findings can be explained by the linear BWBM model in which a growing cerebellar cortex is fluid-like and the sub-cortex is a nongrowing core [19]. The cortex is under tension due to the presence of radial glial cells spanning the cerebellum as well as the pial membrane. Finally, Bergmann glial cells spanning the cortex attempt to maintain constant thickness of the cortex. See Fig. 1. Cell growth in the presence of such constraints drives a featureless convex cortex to form folds. In addition, the BWBM model also offers an explanation for the length scale invariance of the formation of cerebellum folds to understand the conservation of 8-10 primary lobes across vertebrate species spanning a range of sizes [19]. The linear BWBM model provides a quantitative framework for the onset of shape change that manifest as smooth cortex oscillations in the developing cerebellum. However, as the shape of the cerebellum continues to evolve, we observe cusped sulci and wide gyri (see Fig. 2). The linear BWBM model, which can be mapped to a forced, simple harmonic oscillator, cannot account for such nonlinear phenomena. These observations necessitate the exploration of nonlinear elasticity within the BWBM model, particularly in the context of tensioned radial glial cells given that robust nonlinearities have been observed in stretched cells [23]. In addition to exploring the effect of nonlinear springs in the BWBM model, we will also explore the impact of spatial confinement on shape change, which imposes a different form of nonlinearity in the BWBM model. Finally, since we are exploring stages of shape development beyond the onset of shape change, we note the hierarchy of folds that emerges in the developing cerebellum. This hierarchy manifests itself as folds within folds. As the developing cerebellum grows in size, this leads to a type of branching morphogenesis. While branching morphogenesis has been thought of in the context of developing lungs, kidneys, and other organs [24,25], a hierarchical version of BWBM, given its scale invariant solutions, naturally emerges as a potential candidate to explain branching morphogenesis within the cerebellum.
The outline of the manuscript is as follows. In Section II, we review the linear BWBM model focusing on the early stages of shape change in the developing cerebellum. Section III discusses the effect of nonlinear springs that describe the radial fibrous glial cells and its role in sharpening the sinusoidal sulci of the linear model. Section IV discusses the effect of steric hindrances on the flattening of gyri and Section V discusses the hierarchy of subsequent folds emerging at even later stages of cerebellar development. Section VI summarizes the work and addresses its implications. Within a timespan of a day, the featureless, convex developing cerebellum begins to develop folds. The BWBM model provides a physical basis for the onset of these folds. Unlike the cerebrum, the cylindrical symmetry of the cerebellum allows for two-dimensional modeling of its sagittal cross section. The two-dimensional BWBM model [19] offers a quasi-static description of the foliation i.e., it does not describe the dynamics of growth but determines the final equilibrium shape for a given set of parameters. As the parameters change with time, so does the shape. In polar coordinates, the radius of the cerebellum and the thickness of the cortex, otherwise known as the external granular layer (EGL), are represented with r, t respectively (see Fig. 1). The energy functional for the bi-layer model with θ parametrizing the two degrees of freedom is, Here, all the constants k r0 , k t , β, r 0 , t 0 are positive. The first term represents the elastic contribution of the radial glial fibrous cells spanning the cerebellum and the elastic pial membrane surrounding the cerebellum. Moreover, k r is the elasic modulus for this term with r 0 as its rest length radius. The second term is a growth potential for the cortex with k t controlling the contribution of the term and t 0 setting the rest length thickness. The third term represents the Bergmann glial fibrous cells resisting thickness variations of the cortex with β being the accompanying elastic modulus. Note that the third term is not a bending energy term. A bending energy term would involve the curvature as given by the second derivative with θ, unlike the first derivative used here. Finally, given the slower growth of the sub-cortex/core, as compared to the cortex, we demand that the area of the sub-cortex does not change, at least over the time scale of the onset of the foliation, i.e. a day. In other words, It is the area conservation that sets up the competition between the radial elastic energy of the glial fibers and pial membrane with the growth potential term. The extremum of the variational problem subject to this constraint yields a pair of coupled Euler-Lagrange equations for r, t. The solution to these equations are linear sinusoidal functions of the form, where The dimensionless constants are c := k r0 /k t , := µ/k r0 and ρ := k t /β. The linear BWBM model successfully explains the cortex thickness oscillations being out of phase with the subcortex deformation for a developing cerebellum which leads to a thinner cortex at the crest and a thicker cortex at the trough. This model also allows the amplitude of the cortex thickness oscillation to be the same size or even greater than the amplitude of the substrate oscillations. Such phenomena is in contrast to elastic wrinkling models which fail to wrinkle for thicker cortices at the trough [26].

III. NONLINEAR ELASTICITY
While the linear BWBM model addresses the onset of shape change, a next follow up question to ask is -what are the limitations of the linear BWBM model in explaining the more dramatic shape changes in the developing cerebellum at later stages of development? As the radial glial cells and the pial membrane become more stretched by the developing crests, the enhanced stretching may lead to detachment of the radial glial cells from the pial membrane or may lead to nonlinear elastic effects. It has long been known that stretched cells act as nonlinear springs [23] and that collagen, which is one of dominant fibrous proteins constituting the pial membrane exhibits nonlinear elasticity as strain increases [28]. Such effects are not addressed in the linear BWBM model. It would therefore be prudent to examine the role of nonlinear elasticity in the BWBM model. To do so, we promote the radial spring constant k r0 to k r (r), a radially dependent spring 'constant' - The above three terms correspond to quadratic, cubic and quartic energy terms of the radial glial spring potential energy. The cubic energy term brings an asymmetry in the radial spring energy across r 0 and the quartic energy term counteracts the destabilizing effect of the cubic term. The more general form of the uncoupled differential equation where every spring constant is allowed to be nonlinear is explored in Appendix VII A. We nondimensionalize the problem asẼ := E/(k r0 r 2 0 ),r := r/r 0 ,t := t/r 0 ,t 0 := t 0 /r 0 and wherek r1 = k r1 r 0 /k r0 andk r2 = k r2 r 2 0 /k r0 . The nondimensional energy functional with dimensionless variables and coefficients is, whose minimization is subject to the constraint, The variational problem at hand is, The corresponding Euler-Lagrange equations arẽ 1 ρc For the purpose of numerically solving these coupled differential equations, we may as well stop here. However, the exercise of uncoupling the differential equations points to an important source of nonlinearity. Towards finding it, we differentiate Eq. 9 twice with θ to arrive at Here we see the nonlinear term (dr/dθ) 2 . This is a consequence of the inhomogeneous radial spring constant and the coupling between the Euler-Lagrange equations brought about by the Lagrange multiplier in Eq. 8. Solving Eq. 9 fort, we have, Substituting Eq. 12 in Eq. 10 leads to, Substituting, Eq. 13 in Eq. 11, we find the shape equation of the system to be, The importance of the term (dr/dθ) 2 , a source of nonlinearity in Eq. 14, lies in the robustness of its appearance. Irrespective of the form of nonlinearity in k r (r), we retain this nonlinear term whereas the coefficients in Eq. 14 correspondingly vary (see Appendix VII A). Such a nonlinearity is also seen when one attempts to uncouple the Lotka-Volterra equations [29]. We use the RK45 method of scipy.integrate package in python for numerical integration of all differential equations in this paper. For generating the phaseportraits in Fig. 7, we use XPPAUT [30].
Results of numerical integration in Fig. 3a,b, show asymmetric oscillations -sharper sulci and smooth gyri. This holds true even withk r1 = 0, i.e. with no explicit imposed asymmetry in the energy functional (see Eq. 6). This points to the possible role of the (dr/dθ) 2 nonlinearity in bringing about asymmetric oscillations. In the context of explaining the sharp sulci in the cerebral cortex of larger mammals, the sulcification instability in nonlinear elastic materials has been used [7]. Here we observe, sans an elastic instability, sharper sulci in comparison to their sinusoidal counterparts. Even under the influence of these nonlinearities, the system robustly retains the thicker cortical troughs and thinner cortical crests observed in the linear BWBM model. For smallk r1 ,k r2 , the number of folds do not change appreciably in comparison to the linear BWBM model.

A. Assisting-dampening oscillator
In Sec. III, the nonharmonic radial glial springs generate a shape equation 14 which has several sources of nonlinearity. Given the assured presence of the quadratic nonlinearity (dr/dθ) 2 in the shape equation irrespective of the nonlinearity introduced, we seek to isolate the effect it has in determining the shape of the system. Towards that end, we study the simple harmonic oscillator with a velocity dependent force in this section.
The differential equation of motion for a one dimensional, simple harmonic oscillator of mass m with an external nonlinear forcing term, −γ(dx/dt) 2 is written as, The quadratic term assists/dampens the negative/positive velocity respectively bringing an asymmetry in the problem. We refer to this unphysical oscillator as the 'assisting-dampening oscillator' (ADO). The parameters in this equation and the length scale x 0 from the initial conditionx(0) =x 0 can be used to introduce nondimensional variables, x :=x/x 0 and t :=t/ m/k. The equation of motion effectively determined by a single tuning parameter Γ := γx 0 /m is, The nonlinear forcing term is nonconservative since it cannot be written down as the gradient of a position dependent potential. This problem can also not be formulated within the framework of Lagrangian mechanics as this velocity dependent force can not be represented by a function U (x,ẋ) such that, [31]. However, the time reversal symmetry of Eq. 15 ensures all trajectories sufficiently close to the equilibrium point, (x 0 , y 0 ) = (0, 0) to be closed orbits in the phase space [32]. This translates to periodic motion in the position-time space (see Fig. 7).
The second order differential equation 15 can be rep-resented as two coupled first order differential equations, Dividing the set of equations, we have, In contrast to the differential equation of motion in position-time space (see Eq. 15), the phase-space orbit is represented by a first-order differential equation. This makes finding an exact solution to the orbit in phasespace simpler. The variable transformation, u = y 2 linearizes the above differential equation as, The solution to the homogeneous differential equation, is u(x) = A e −2Γx where A is a constant that needs to be fixed by the initial condition. The particular solution to Eq. 18 is u(x) = −x/Γ + 1/(2Γ 2 ). Reverting to y, the total solution satisfying the initial condition y(−1) = 0 is, In the limit of Γ → 0, we have y(x) ∼ ± √ 1 − x 2 which are semi-circular arcs as expected for a simple harmonic oscillator.
The numerical integration of the differential equation is presented in Fig. 3. The sharp troughs in the positiontime space are solely due to the effect of the quadratic nonlinear term (dx/dt) 2 .
The red lines in Figs. 3e, 7 are velocity nullclines which are the locus of points in phase-space where the acceleration, dy/dt = d 2 x/dt 2 = 0. From Eq. 15, it follows that the nullcline is a quadratic function in x, It is interesting to note that the quadratic nullcline (see Fig. 7b) bears resemblance to a similarly curved nullcline of the BWBM with nonharmonic springs (see Fig.  7d). The nullcline's curvature is dictated by the same quadratic nonlinear term (dx/dt) 2 in both cases. For comparison with this unphysical oscillator, we also study a conservative Hamiltonian system and observe that it is also capable of exhibiting oscillations with sharp troughs. The simplest such system is obtained by adding cubic and quartic potential energy terms to the simple harmonic oscillator. The former term provides an explicit asymmetry in the energy functional, and the latter term ensures stability about the equilibrium point. In Appendix VII B, this system is explored numerically and analytically using a perturbation series constructed via the Lindstedt-Poincaré method [33].

B. A measure for crest -trough asymmetry
Stokes, in his study of propagating waves approaching the shore, discussed the development of narrow crests and wider troughs [34]. In this context, measures for velocity and acceleration skewness have been proposed and studied. In the folds of cerebellum, which is a stationary spatial oscillation, a similar asymmetry presents itself in the form of wide gyri/crest and sharp sulci/trough. We propose the following measure for quantifying the asymmetry in widthscrest -trough asymmetry := t crest − t trough t crest + t trough .
Here, the length the parameter θ traverses between consecutive pair of points at which d 2 x/dt 2 = 0 on the climbing (falling) wave, and the immediately following falling (climbing) wave defines t crest (t trough ), respectively (see Fig. 3d). These points are the points of intersection of the phase-space orbit with the velocity nullcline (see Fig.  3e). The horizontal mouthed parabolic-shaped nullclines, thus influence the position of these points of intersection and play a prominent role in bringing the asymmetry between the widths of the crests and troughs. For the ADO, the velocity nullcline is exactly parabolic (see Eq. 21 and Fig. 7). The parabolic shape of the nullcline is due to the quadratic nonlinearity (dx/dt) 2 in Eq. 15. For the BWBM model with nonlinear radial glial springs, the nullcline will not be exactly parabolic due to the presence of other nonlinearities in Eq. 14.
The measure in Eq. 22 is bounded within (−1, 1). For a sinusoidal wave, which is perfectly symmetric, the measure equals zero. For the case of studying this asymmetry in the lobes of the nonlinear BWBM model, we use θ in lieu of t in Eq. 22. Studying this asymmetry for the ADO and the nonlinear BWBM model for a range of tuning parameters, we see that the asymmetry scales as 1.04 ± 0.05 and the nonlinear BWBM model shows a parameter dependent scaling of 0.61 ± 0.02 and 0.79 ± 0.08 for two chosen sets of parameters of 5 lobed and 6 lobed systems respectively (see Fig. 3f).

IV. SPATIAL CONFINEMENT
A cerebellum does not grow in isolation. It encounters steric effects from the cerebrum, brain-stem, and the skull. The effect of the skull on a growing cerebrum has earlier been studied computationally [35]. We, therefore, are compelled to explore the effects of steric confinement with the BWBM model, particularly since there may be interplay between the area conservation of the sub-cortex and the imposed steric effects at a radial boundary. More specifically, if the lobes are not allowed to grow radially, they may grow tangentially, making way for sharper folds in the cerebellum. To check for this, we model the steric effects on the developing cerebellum by incorporating a logistic function (1+tanh(x))/2 into the energy functional, The constants of the logistic function are omitted as they vanish in the resulting Euler-Lagrange equations. Here, K c is the coupling constant, q −1 is the width of the step size of the tanh function and r c is the radial position of the step. The strength of the steric interaction is taken to be K c q 2 , a quantity whose dimensions matches those of k r . Renormalizing the parameters and variables by q −1 avoids singularities in the Euler-Lagrange equations in the limit of q → 0. Dividing Eq. 23 by k r q −2 , we have, whereẼ := E/(k r q −2 ),r := qr,t := qt,r c := qr c and K c = K c q 2 /k r . All parameters and variables are now rendered dimensionless. Given that there is no explicit dependence on q, we can be assured that it won't show up in the Euler-Lagrange equation either. The corresponding area conserving constraint is, The Euler-Lagrange equation is, Here γ = ρ( c + 1). We, again, observe the (dr/dθ) 2 term. However, its coefficient is zero for all r = r c , thus effectively localizing its effect. Results of numerical integration in Fig. 4a,b shows local flattening of the crests/gyri at contact with the confining wall. Interestingly, there are no nonlocal effects of the confining wall on the lobes, i.e., steric effects do not contribute to the sharpness of the troughs/sulci, at least for the parameters we study. The area conservation on the sub-cortex is not a strong enough constraint to effect the shape of the lobes near the troughs as the sub-cortex can still change its shape.

V. A BRANCHING HIERARCHY
Mammalian cerebellums are seen to develop folds irrespective of the size of the organ [20]. This is in contrast to small mammalian cerebrums which do not develop folds [7,36]. A feature of the linear BWBM model is that it offers an explanation for the length scale independence of the folding morphogenesis in the cerebellum. In the limit of small of the linear BWBM model, the number of primary folds N scales as √ ρ. Thus, in this limit, the number of folds developed in the cerebellum is independent of the size of the cerebellum and is determined solely by the material elastic constants. If the material elastic constants do not dramatically differ across mammalian species, the linear BWBM model can potentially explain the conservation in the number of primary folds. In this section, we discuss how the linear BWBM can be used to describe the branching hierarchy within a given cerebellum.
As the gyri/crests of the cerebellar cortex continue to grow, the sulci/troughs sharpen and become anchored [37]. The anchoring is due to a combination of the radial glial cells and the pial basement membrane, a thin sheet of extracellular matrix (ECM) made up of collagen, laminin, and other ECM components [38]. It is known that basement membranes stick together [39], so as the troughs sharpen, the pial membrane from each side of the trough begins to come into contact and stick to reinforce the anchoring. The resulting anchoring centers delineate the petal-like projections called lobules [27].
We hypothesize that when the anchoring centers serve  Fig. 7b of Ref. [22] published under the terms of the Creative Commons Attribution. b) Idealized geometrical example of hierarchy manifests as a fractal. c) Representation of branching morphogenesis implemented using the geometrical idea of (b) and the 4 lobe BWBM systems with elliptical r0 of (d). Each lobe becomes a sub-system of its own and spawns its own folds. Here, the first and each of the second generation systems are generated numerically with the second generation overlaying the first after a rescaling factor of 0.28. as effective boundaries between the lobules, each lobe becomes its own subsystem. This subsystem then consists of its own subcortex/subcore, all within the encompassing primary cortex/core geometry. Some of these featureless subsystems go on to develop folds of their own to, in turn, generate another generation of subsystems and so on. See Fig. 5a. In other words, as the cerebellum continues to develop, a branching hierarchy of subsystems emerges. This branching hierarchy is yet another distinguishing feature of the cerebellum that sets it apart from the cerebrum. It is to be noted that even though cerebellums of all sizes demonstrate folds, smaller cerebellums have fewer hierarchial branchings.
The hierarchical generation of lobules within lobules points to a scale-invariant branching process. Given that within the framework of BWBM, the formation of crests and troughs do not depend on system size, this framework offers a natural description for the hierarchy. As the size of the subsystems decreases with each successive generation, crests and troughs can still form as long as the material properties do not change. The validity of the two dimensional model to describe growth in three dimensions remains valid since the cerebellum retains its cylindrical symmetry during development [22].
As an idealized, purely geometric example of the hierarchy, we consider an initial zeroth generation circle. Along its perimeter, six first generation circles are generated. This process can proceed ad infinitum, to generate a fractal structure with a fractal dimension of log(3)/log (2). We show four generations of this hierarchy in Fig. 5b.
In the same vein, using the linear BWBM model, we represent the hierarchy of folds in the cerebellum in Fig.  5c. The 'zeroth' generation is the circular r 0 (not shown in figure) which generates six first generation foliations.
Each lobe formed within consecutive folds is now considered an independent subsystem and sets the length scale for the r 0 of the following, second generation of lobes. Fig. 5a suggests the free part of the first generation lobe i.e. the part that is not sticking to its neighbors, to be 'elliptical' with the major axis of the ellipse being parallel to the outermost exposed edge of the lobe. This is especially evident for the L678 lobe. To accomodate this visual observation in generating second generation folds, we employ a geometric non-linearity in the form of an elliptical r 0 with an eccentricity e within the BWBM model. We assume there are no residual stresses and generate a 4 lobed BWBM system as in Fig. 5d. We then use only the top half of this solution to represent the second generation lobes in Fig. 5c.
It is also interesting to note that in Fig. 5d, the most prominent fold occurs transversely to the horizontal major axis of the elliptical r 0 . This is simply the consequence of the system trying to minimize its energy contribution from the radial glial springs. For comparison, we show in Fig. 5e, a four lobed system generated formed from a circular r 0 .

VI. DISCUSSION
Inspired by cerebellar shape development, we study the effects of nonlinear elasticity, steric confinement, and a branching hierarchy within the BWBM model. Exploring the effects of nonlinear elasticity of the fibrous radial glial cells, the interplay between geometry and nonlinearity is seen to give rise to troughs sharper than the troughs obtained in the linear BWBM model and we arrive at an asymmetry between the crests (gyri) and the troughs (sulci). This asymmetry can be understood thus: the relatively slow growth of the sub-cortex is taken into account by demanding area conservation of the sub-cortex in the two-dimensional BWBM model. The associated Lagrange multiplier couples the radius of the cerebellum and the thickness of the cortex. Nonlinear radial springs, in association with this coupling, results in the robust quadratic nonlinear term of the form (dr/dθ) 2 in the shape equation. To illustrate the role of the quadratic nonlinear term in sharpening the sulci, we study the simple harmonic oscillator with the same form of nonlinearity and observe its sufficiency in achieving sharp sulci. Several other nonlinearities emerge in the shape equation including a spatially-varying effective 'mass' coefficient.
The perspective of cerebellum foliation as the action of a nonlinear oscillator can be a useful one given the extensive theoretical studies of such oscillators [40,41]. For BWBM of the cerebellum, the linear model with constant r 0 maps to a forced harmonic oscillator and, for small eccentricities of r 0 , maps to an unconventional Duffing oscillator. For nonlinearK r (r), we attempt to understand the corresponding nonlinearity in the context of the assisting-dampening oscillator. We hope the study of cerebellar foliation as a nonlinear oscillator problem continues to be fruitful. In a related work, the existence of a new morphological instability in confined nonlinear elastic sheets was found in the context of a period-doubling bifurcation, exhibiting an analogy with parametric resonance in another nonlinear oscillator [42].
The period-doubling hierarchy found in the Ref. [42] is very different from the new hierarchy found here. The hierarchy found here is one due to boundary conditions in the form of anchoring centers to create sub-regions, or sub-systems, from which the same type of scale-invariant foliation emerges, at least in the limit of small . Had the foliation mechanism not been scale-invariant in any limit, such as with a purely elastic system, the smaller sub-lobes would soon become featureless as the number of foliations depend linearly on the perimeter of the subsystem. Within BWBM, therefore, we have identified a new scale-invariant branching morphogenesis mechanism. It is not yet clear how generic this new branching hierarchy mechanism is in terms of moving beyond the cerebellum. Ref. [19] addressed potential applications of BWBM to two-dimensional brain organoids [43] and the developing retina [44].
Referring again to Fig. 5(a), one of the second generation sublobes labelled L45 does not branch. Perhaps the material properties are altered in this sublobe so that features do not form. For sublobe L678, the two new sub-sublobes are not similar in size. This could be due to changes in curvature of the confinement from growing, surrounding tissue. So far, we have only addressed steric, static confinement. Certainly, such variabilities from sublobe to sublobe can be explored in less idealized conditions.
We note that within the context of purely elasticity theory, an explanation for the thicker-sulci/thinner-gyri of the developing cerebellar cortex was recently achieved by the addition of surface tension [45]. While more energetic contributions can certainly be added to either the purely elastic model or to the BWBM model, the recent experimental observations needs to be incorporated in the modeling -the cells in the cortex are motile with cellular rearrangements on the time scale of minutes [22] and the cerebellum is under tension as it develops (as opposed to compression). These observations render a purely elastic model suspect. However, the differential growth between the cortex and sub-cortex remains central in both classes of models. To make progress, we need further experimental falsification tests to rule out classes of models.
Finally, given our focus on the cerebellum here, one is led to wonder whether some form of BWBM is applicable to the cerebrum. As mentioned previously, the cerebrum and the cerebellum have rather different morphologies. In terms of development in the cerebellum, the predominant growth of the cells is in the cortex [37]. Many such cells migrate inward to become part of the core of the cerebellum. In the cerebrum, much of the cell proliferation is in the core and the progenitor cells migrate outward to become part of the cortex [46,47] . In this sense, the two organs are inverse to each other. Given the presence of motile cells in the developing cerebrum, one may wonder whether a purely elastic approach to the developing brain is reasonable. Without a doubt, the time scales of cell migration decide the contest between elasticity and fluidity. Under the framework of liquid crystals, earlier work on the developing cerebrum arrived at a prediction for the bending modulus of the cortex [36]. This approach was based on a revised view of the axonal tension model for the developing cerebrum [48,49]. These novel approaches have the potential to build new inroads in quantitative biology.
MCG acknowledges useful discussions with Manu Mannattil. JMS acknowledges financial support from NSF-DMR-1832002 and an Isaac Newton Award from the DoD.
The perturbative expansions are compared with numerical solution in Fig. 6.

C. Phase portraits
Closed orbits in phase-space translate to periodic motion in the position-time space. Plotting orbits in phase space for different initial conditions builds the phase portrait of the system. The x (dx/dt) nullcline are the locus of points where dx/dt = 0 (d 2 x/dt 2 = 0). The nullclines are generally used to divide the phase-portrait into regions where the tangent of the orbit points in the same general direction (NW, NE, SE or SW). The orbit has a vertical (horizontal) tangent when it passes through the x (dx/dθ) nullcline. The intersection of the nullclines gives the equilibrium point where both the 'velocity' and 'acceleration' of the system is zero. The time-reversal symmetry of the governing differential equations (see Eq. 14, 15, 26) guarantees closed orbits at points close enough to the equilibrium point. This is verified in the phaseportraits in Fig. 7.