Buckling without bending morphogenesis: Nonlinearities, spatial confinement, and branching hierarchies
BBuckling without bending morphogenesis:Nonlinearities, spatial confinement, and branching hierarchies
M. C. Gandikota ∗ and J. M. Schwarz
1, 2, † Physics Department and BioInspired Institute, Syracuse University, Syracuse, NY 13244 USA Indian Creek Farm, Ithaca, NY 14850 USA (Dated: February 18, 2021)During morphogenesis, a featureless convex cerebellum develops folds. As it does so, the cortexthickness is thinnest at the crest (gyri) and thickest at the trough (sulci) of the folds. This observa-tion cannot be simply explained by elastic theories of buckling. A recent minimal model explainedthis phenomenon by modeling the developing cortex as a growing fluid under the constraints ofradially spanning elastic fibers, a plia membrane and a nongrowing sub-cortex (Engstrom, et. al.,PRX
I. INTRODUCTION
When morphogenesis is viewed through the lens of aphysicist, one of the typical mechanistic routes to take ismorphoelasticity induced by varying internal stresses [1].Shape change in response to internal stresses, within elas-tic objects which tend to retain shape, provides us witha large collection of shapes. Differential growth can be asource of such an internal stress. This reasonable view-point is made manifest in the widespread use of the Eu-ler buckling instability in purely elastic materials to ex-plain the onset of folds in morphogenesis problems [2]as diverse as cerebra [3–8], intestinal crypts and villi[9, 10], airway mucus wrinkles [11, 12], tooth ridges [13]and hair-follicle patterns [14]. Recent work accountingfor cerebellar foliation falls under the same framework[15, 16]. These models typically predict a characteristiclength scale between folds and state quantitative agree-ment between prediction and observation to validate suchan approach.On the other hand, tissue fluidity in early stages ofdeveloping embryos has emerged as a driver of shapechange in both zebrafish and the insect
Tribolium cas-taneum [17, 18]. Fluids, unlike elastic solids, do not re-tain their shape. Then how does fluidity - the antithesisof elasticity - affect the shape of the developing organ?Presumably biology has figured out ways to combine el-ements of elasticity and fluidity to bring about an evenmore complex collection of shapes at later stages of devel-opment. The buckling without bending morphogenesis(BWBM) model is one such model affirming the clev-erness of biology [19]. Its origin is rooted in explaining ∗ [email protected] † [email protected] recent quantitative observations of the developing murinecerebellum, a much less studied component of the brainas compared to the cerebrum.Vertebrate brains are typically divided into three dif-ferent sectors – the forebrain, the midbrain, and the hind-brain. While the forebrain consists of the cerebrum, thehypothalamus, and the thalamus, the hindbrain containsthe spinal cord, medulla oblongata, and the cerebellum.Since the shape of the forebrain varies more across speciesthan the midbrain and hindbrain, its development hasbeen the focus of much study. On the other hand, bystudying conserved patterns in the hindbrain, we gaincomplementary insight into the morphogenetic processesof the brain.To be specific, let us characterize the difference inshape between the cerebrum and cerebellum. First, theshape of the cerebellum has an approximate cylindricalsymmetry, whereas such a symmetry is absent in the cere-brum. Second, mammalian cerebellums of all sizes pro-duce folds [20], whereas small mammalian cerebrums donot show folds [21]. Given that the overall size of both or-gans increases with increasing species size, it is possiblethat the differences in shape arises via different physi-cal shape change mechanisms. If so, do different physi-cal shape change mechanisms lead to different emergentfunctionalities of the two? Moreover, the conservationof the number of the 8-10 primary lobes of the cerebel-lum across species demands a treatment on its own foot-ing as it suggests a scale invariant shape change mech-anism. Number of secondary folds however, can differacross species.Recent experimental observations of the developingmurine cerebellum found that the proliferating cells inthe cerebellar cortex are motile with neighbor exchangeson the order of minutes [22]. It was also observed that thedeveloping cerebellar cortex varies in thickness with the a r X i v : . [ c ond - m a t . s o f t ] F e b trend being that the cortex is thinnest at the crest (gyri)and thickest at the trough (sulci). Moreover, the devel-oping cerebellum is under tension and not under com-pression as evidenced by both radial and circumferentialcuts. While these observations cannot be explained byelastic wrinkling theories, all three of these findings canbe explained by the linear BWBM model in which a grow-ing cerebellar cortex is fluid-like and the sub-cortex is anongrowing core [19]. The cortex is under tension due tothe presence of radial glial cells spanning the cerebellumas well as the pial membrane. Finally, Bergmann glialcells spanning the cortex attempt to maintain constantthickness of the cortex. See Fig. 1. Cell growth in thepresence of such constraints drives a featureless convexcortex to form folds. In addition, the BWBM model alsooffers an explanation for the length scale invariance ofthe formation of cerebellum folds to understand the con-servation of 8-10 primary lobes across vertebrate speciesspanning a range of sizes [19].The linear BWBM model provides a quantitativeframework for the onset of shape change that manifestas smooth cortex oscillations in the developing cerebel-lum. However, as the shape of the cerebellum contin-ues to evolve, one observes cusped sulci and wide gyri(see Fig. 2). Robust nonlinearities have been earlier ob-served in stretched cells [23]. The linear BWBM model,which can be mapped to a forced, simple harmonic os-cillator, cannot account for such nonlinear phenomena.These observations necessitate the exploration of nonlin-ear elasticity within the BWBM model, particularly inthe context of the tensioned radial glial cells. In additionto exploring the effect of nonlinear springs in the BWBMmodel, we will also explore the impact of spatial confine-ment on shape change, which imposes a different form ofnonlinearity in the BWBM model. Finally, since we areexploring stages of shape development beyond the on-set of shape change, we note the hierarchy of folds thatemerges in the developing cerebellum. This hierarchymanifests itself as folds within folds. As the developingcerebellum grows in size, this leads to a type of branch-ing morphogenesis. While branching morphogenesis hasbeen thought of in the context of developing lungs, kid-neys, and other organs [24, 25], a hierarchical versionof BWBM, given its scale invariant solutions, naturallyemerges as a potential candidate to explain branchingmorphogenesis within the cerebellum.The outline of the manuscript is as follows. In SectionII, we review the linear BWBM model focusing on theearly stages of shape change in the developing cerebel-lum. Section III discusses the effect of nonlinear springsthat describe the radial fibrous glial cells and its role insharpening the sinusoidal sulci of the linear model. Sec-tion IV discusses the effect of steric hindrances on theflattening of gyri and Section V discusses the hierarchyof subsequent folds emerging at even later stages of cere-bellar development. Section VI summarizes the work andaddresses its implications. II. LINEAR BUCKLING WITHOUT BENDINGMORPHOGENESIS MODEL
FIG. 1.
Schematic of the buckling without bending morpho-genesis (BWBM) model.
The radial glial cells span the cere-bellum and Bergmann glial cells span the cortex and resistthickness variations in the cortex. The area of the sub-cortexdoes not grow as fast as the cortex and so it is approximatedas nongrowing over some time scale.
Within a timespan of a day, the featureless, convex de-veloping cerebellum begins to develop folds. The BWBMmodel provides a physical basis for the onset of thesefolds. Unlike the cerebrum, the cylindrical symmetry ofthe cerebellum allows for two-dimensional modeling ofits sagittal cross section. The two-dimensional BWBMmodel [19] offers a quasi-static description of the folia-tion i.e., it does not describe the dynamics of growth butdetermines the final equilibrium shape for a given set ofparameters. As the parameters change with time, so doesthe shape. In polar coordinates, the radius of the cerebel-lum and the thickness of the cortex, otherwise known asthe external granular layer (EGL), are represented with r, t respectively (see Fig. 1). The energy functional forthe bi-layer model with θ parametrizing the two degreesof freedom is, E (cid:20) r, t, dtdθ (cid:21) = (cid:90) dθ (cid:40) k r ( r − r ) − k t ( t − t ) + β (cid:18) dtdθ (cid:19) (cid:41) . (1) Here, the first term represents the elastic contribution ofthe radial glial fibrous cells spanning the cerebellum andthe elastic pial membrane surrounding the cerebellum.The second term is a growth potential for the cortex andthe third term represents the Bergmann glial fibrous cellsresisting thickness variations of the cortex. Note thatthe third term is not a bending energy term. A bendingenergy term would involve the curvature as given by thesecond derivative with θ , unlike the first derivative usedhere. All coefficients in the energy functional are positive.Finally, given the slower growth of the sub-cortex/core,as compared to the cortex, we demand that the area ofthe sub-cortex does not change, at least over the timescale of the onset of the foliation, i.e. a day. In otherwords, 12 (cid:90) dθ ( r − t ) = A = constant . (2)It is the area conservation that sets up the competitionbetween the radial elastic energy of the glial fibers andpial membrane with the growth potential term.The extremum of the variational problem subject to thisconstraint yields a pair of coupled Euler-Lagrange equa-tions for r, t . The solution to these equations are linearsinusoidal functions of the form, t ( θ ) = T sin( nθ + φ ) ,r ( θ ) = (cid:18) − (cid:15) (cid:19) r ( θ ) − (cid:18) (cid:15) − (cid:15) (cid:19) t ( θ ) , (3)where n = (cid:112) ρ (1 + (cid:15)c/ (1 − (cid:15) )) ,T = √ (cid:15) (cid:115) A π − (cid:18) a − t − (cid:15) + (cid:15)c (cid:19) . The dimensionless constants are c := k r /k t , (cid:15) := µ/k r and ρ := k t /β .The linear BWBM model successfully explains the cor-tex thickness oscillations being out of phase with the sub-cortex deformation for a developing cerebellum whichleads to a thinner cortex at the crest and a thicker cor-tex at the trough. This model also allows the amplitudeof the cortex thickness oscillation to be the same size oreven greater than the amplitude of the substrate oscilla-tions. Such phenomena is in contrast to elastic wrinklingmodels which fail to wrinkle for thicker cortices at thetrough [26]. Another feature of this model is that it offersan explanation for the length scale independence of thefolding morphogenesis in the cerebellum. Mammaliancerebellums are seen to develop folds irrespective of thesize of the organ [20]. This is in contrast to small mam-malian cerebrums which do not develop folds [27, 28]. Inthe limit of small (cid:15) of the linear BWBM model, the num-ber of primary folds N scales as √ ρ . Thus, in this limit,the number of folds developed in the cerebellum is inde-pendent of the size of the cerebellum and is determinedsolely by the material elastic constants. If the materialelastic constants do not dramatically differ across mam-malian species, the linear BWBM model can potentiallyexplain the conservation in the number of primary folds. III. NONLINEAR ELASTICITY
While the linear BWBM model addresses the onset ofshape change, a next follow up question to ask is – whatare the limitations of the linear BWBM model in explain-ing the more dramatic shape changes in the developingcerebellum at later stages of development? As the radialglial cells and the pial membrane become more stretchedby the developing crests, the enhanced stretching maylead to detachment of the radial glial cells from the pialmembrane or may lead to nonlinear elastic effects. Ithas long been known that stretched cells act as nonlinear
FIG. 2.
Visual comparison of experiment and model.
Themorphogenesis of cerebellum involves the development of foldsin an initially featureless convex cerebellum. The top rowshows midsagittal cross-sections of mouse cerebellum prior tobirth. Scale bar is 200 µ m. The arrowheads are from the orig-inal figures in Ref. [29] and indicate fissures. Reprinted withpermission from Springer Nature. The linear BWBM modelreproduces the smooth sinusoidal like folds that occur at theonset of folding. The BWBM model with nonharmonic radialglial springs reproduces the sharp cusped sulci that we see atthe later stages of cerebellar development. Note that time isnot explicit in the quasi-static BWBM. We start from sys-tem parameters describing the unstable state and minimizethe energy while obeying constraints to arrive at the nonlin-ear BWBM configuration. We use brain folding terminology- gyri (sulci) and oscillator terminology - crests (troughs) in-terchangeably to refer to hills (valleys) of oscillations respec-tively. Parameters: c = 0 . , . , (cid:15) = 0 . , ρ = 31 . , ˜ t =0 . , . , ˜ k r = 0 , − . , ˜ k r = 0 , .
05 for the linear and non-linear BWBM respectively. [˜ r, d ˜ r/dθ ] θ =0 is [2 . , − . springs [23] and that collagen, which is one of dominantfibrous proteins constituting the pial membrane exhibitsnonlinear elasticity as strain increases [30]. Such effectsare not addressed in the linear BWBM model. It wouldtherefore be prudent to examine the role of nonlinearelasticity in the BWBM model. To do so, we promotethe radial spring constant k r to k r ( r ), a radially depen-dent spring ‘constant’ - k r ( r ) := k r + k r ( r − r ) + k r ( r − r ) . (4)The above three terms correspond to quadratic, cubicand quartic energy terms of the radial glial spring poten-tial energy. The cubic energy term brings an asymmetryin the radial spring energy across r and the quartic en-ergy term counteracts the destabilizing effect of the cubicterm. The more general form of the uncoupled differen-tial equation where every spring constant is allowed to benonlinear is explored in Appendix VII A. We nondimen-sionalize the problem as ˜ E := E/ ( k r r ) , ˜ r := r/r , ˜ t := t/r , ˜ t := t /r and˜ K r ( r ) := k r ( r ) k r = 1 + ˜ k r (¯ r −
1) + ˜ k r (¯ r − , (5)where ˜ k r = k r r /k r and ˜ k r = k r r /k r . The nondi-mensional energy functional with dimensionless variablesand coefficients is,˜ E (cid:20) ˜ r, ˜ t, d ˜ tdθ (cid:21) = (cid:90) dθ (cid:26) ˜ K r ( r )(˜ r − − c (˜ t − ˜ t ) + 1 ρc (cid:18) d ˜ tdθ (cid:19) (cid:27) , (6)whose minimization is subject to the constraint,12 (cid:90) dθ (˜ r − ˜ t ) = A r = dimensionless constant . (7)The variational problem at hand is, δ (cid:20) ˜ E − (cid:15) (cid:90) dθ (˜ r − ˜ t ) (cid:21) = 0 (8)The corresponding Euler-Lagrange equations are˜ K r ( r )(¯ r − − (cid:15) (¯ r − ˜ t ) + ˜ K (cid:48) r ( r )2 (¯ r − = 0 , (9)1 ρc d ˜ tdθ + 1 c (˜ t − ˜ t ) − (cid:15) (˜ r − ˜ t ) = 0 . (10)For the purpose of numerically solving these coupled dif-ferential equations, we may as well stop here. How-ever, the exercise of uncoupling the differential equationspoints to an important source of nonlinearity. Towardsfinding it, we differentiate Eq. 9 twice with θ to arrive at (cid:20) ˜ K r ( r ) + 2(¯ r −
1) ˜ K (cid:48) r ( r ) + (¯ r − K (cid:48)(cid:48)(cid:48) r ( r ) − (cid:15) (cid:21) (cid:18) d ¯ rdθ (cid:19) + (cid:20) K (cid:48) r ( r ) + 3(¯ r −
1) ˜ K (cid:48)(cid:48) r ( r ) + (¯ r − K (cid:48)(cid:48)(cid:48) r ( r ) (cid:21) (cid:18) d ¯ rdθ (cid:19) = − (cid:15) d ˜ tdθ . (11)Here we see the nonlinear term ( d ˜ r/dθ ) . This is a con-sequence of the inhomogeneous radial spring constantand the coupling between the Euler-Lagrange equationsbrought about by the Lagrange multiplier in Eq. 8. Solv-ing Eq. 9 for ˜ t ,˜ t = ¯ r − (cid:15) (cid:34) (¯ r −
1) ˜ K r ( r ) + ˜ K (cid:48) r ( r )2 (¯ r − (cid:35) . (12)Substituting Eq. 12 in Eq. 10 leads to d ˜ tdθ = ρ (cid:20)(cid:18)(cid:26) (¯ r −
1) ˜ K r ( r ) + ˜ K (cid:48) r ( r )2 (¯ r − (cid:27) − ¯ r (cid:19) × (cid:18) (cid:15)c(cid:15) (cid:19) + (cid:15)c ¯ r + ˜ t (cid:21) . (13) Substituting, Eq. 13 in Eq. 11, we find the shape equa-tion of the system to be, (cid:20) ˜ K r ( r ) + 2(¯ r −
1) ˜ K (cid:48) r ( r ) + (¯ r − K (cid:48)(cid:48) r ( r ) − (cid:15) (cid:21) (cid:18) d ˜ rdθ (cid:19) + (cid:20) K (cid:48) r ( r ) + 3(¯ r −
1) ˜ K (cid:48)(cid:48) r ( r ) + (¯ r − K (cid:48)(cid:48)(cid:48) r ( r ) (cid:21) (cid:18) d ˜ rdθ (cid:19) − ρ(cid:15) ¯ r + ρ (1 + (cid:15)c ) (cid:20) (¯ r −
1) ˜ K r ( r ) + (¯ r − K (cid:48) r ( r ) (cid:21) + ρ(cid:15)t = 0 . (14)The importance of the term ( d ˜ r/dθ ) , a source of non-linearity above lies in the robustness of its appearance.Irrespective of the form of nonlinearity in k r ( r ), we re-tain this nonlinear term whereas the coefficients in Eq.14 correspondingly vary (see Appendix VII A). Such anonlinearity is also seen when one attempts to uncouplethe Lotka-Volterra equations [31].We use the RK45 method of scipy.integrate pack-age in python for numerical integration of all differen-tial equations in this paper. For generating the phase-portraits in Fig. 7, we use XPPAUT [32].Results of numerical integration in Fig. 3a,b, showasymmetric oscillations - sharper sulci and smooth gyri.This holds true even with ˜ k r = 0, i.e. with no explicitimposed asymmetry in the energy functional (see Eq. 6).This points to the possible role of the ( d ˜ r/dθ ) nonlinear-ity in bringing about asymmetric oscillations. In the con-text of explaining the sharp sulci in the cerebral cortexof larger mammals, the sulcification instability in nonlin-ear elastic materials has been used [7]. Here we observe,sans an elastic instability, sharper sulci in comparison totheir sinusoidal counterparts. Even under the influenceof these nonlinearities, the system robustly retains thethicker cortical troughs and thinner cortical crests ob-served in the linear BWBM model. For small ˜ k r , ˜ k r ,the number of folds do not change appreciably in com-parison to the linear BWBM model. A. Assisting-dampening oscillator
In Sec. III, the nonharmonic radial glial springs gen-erate a shape equation 14 which has several sources ofnonlinearity. Given the assured presence of the quadraticnonlinearity ( d ˜ r/dθ ) in the shape equation irrespectiveof the nonlinearity introduced, we seek to isolate the ef-fect it has in determining the shape of the system. To-wards that end, we study the simple harmonic oscillatorwith a velocity dependent force in this section.The differential equation of motion for a one dimen-sional, simple harmonic oscillator of mass m with an ex- FIG. 3.
Asymmetry in the widths of gyri (crests) and sulci (troughs). a,b) Numerical solutions in polar and Cartesiancoordinates for BWBM with nonharmonic radial glial springs. Nonsinusoidal sharp sulci can be seen for systems with nonzero˜ k r , ˜ k r (green, blue, yellow). c) Phase space orbits of solutions in (a,b). The left-right asymmetry in the orbits seen as asteeper drop on the left translates to sharper troughs in the coordinate space. d) Numerical solution in Cartesian coordinatesfor the ADO oscillator. Sharper troughs are observed due to the addition of the nonlinear force term ( dx/dt ) to the simpleharmonic oscillator. The acceleration of the oscillator is zero at the points marked with the asterisk. The width of crest andtrough is shown. e) Numerical and analytical solution of the phase-space orbit for the ADO. The velocity nullcline, in red, isthe locus of points with zero acceleration. The position coordinate (marked in asterisk) at which the nullcline intersects theorbit decides the widths of crests and troughs in (d). The shape of the orbit is seen to be especially similar to the orbits in (c)for the cases of ˜ k r = 0 , .
19. f) Crest - trough asymmetry for the ADO and the BWBM model as a function of their respectivetuning parameters: Γ and − ˜ k r = ˜ k r . The dashed black line is the upper limit of the asymmetry measure. Parameters c = 0 . , (cid:15) = 0 . , ρ = 15 . , ˜ t = 0 .
7. For periodicity in (a), parameter c is appropriately tuned forchoices of ˜ k r , . 6-lobe BWBM in (f): c = 0 . , (cid:15) = 0 . , ρ = 31 . , ˜ t = 0 .
5. [˜ r, d ˜ r/dθ ] θ =0 is [1.1,0],[1.25,0] for 5-/6- lobe. ternal nonlinear forcing term, − γ ( d ˜ x/dt ) is written as, m d ˜ xd ˜ t = − k ˜ x − γ (cid:18) d ˜ xd ˜ t (cid:19) . The quadratic term assists/dampens the nega-tive/positive velocity respectively bringing an asymmetryin the problem. We refer to this unphysical oscillatoras the ‘assisting-dampening oscillator’ (ADO). Theparameters in this equation and the length scale x fromthe initial condition x (0) = x can be used to introducenondimensional variables, x := ˜ x/x and t := ˜ t/ (cid:112) m/k .The equation of motion effectively determined by asingle tuning parameter Γ := γx /m is, d xdt = − x − Γ (cid:18) dxdt (cid:19) . (15)The nonlinear forcing term is nonconservative since itcannot be written down as the gradient of a position dependent potential. This problem can also not be for-mulated within the framework of Lagrangian mechanicsas this velocity dependent force can not be representedby a function U ( x, ˙ x ) such that, − Γ ˙ x = − ∂U∂x + ddt (cid:18) ∂U∂ ˙ x (cid:19) , where ˙ x = dx/dt [33]. However, the time reversal sym-metry of Eq. 15 ensures all trajectories sufficiently closeto the equilibrium point, ( x , y ) = (0 ,
0) to be closedorbits in the phase space [34]. This translates to periodicmotion 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, dxdt = y,dydt = − x − Γ y . (16)Dividing the set of equations, we have, y dydx = − x − Γ y . (17)In contrast to the differential equation of motion inposition-time space (see Eq. 15), the phase-space orbitis represented by a first-order differential equation. Thismakes finding an exact solution to the orbit in phase-space simpler. The variable transformation, u = y lin-earizes the above differential equation as, dudx + 2Γ u = − x. (18)The solution to the homogeneous differential equation, dudx + 2Γ u = 0 , (19)is u ( x ) = A e − x where A is a constant that needs tobe fixed by the initial condition. The particular solutionto Eq. 18 is u ( x ) = − x/ Γ + 1 / (2Γ ). Reverting to y , thetotal solution satisfying the initial condition y (1) = 0 is, y ( x ) = 1 √ (cid:113) − x − (1 − e − x − . (20)In the limit of Γ →
0, we have y ( x ) ∼ ±√ − x whichare semi-circular arcs as one expects for a simple har-monic oscillator.The numerical integration of the differential equationis presented in Fig. 3. The sharp troughs in the position-time space are solely due to the effect of the quadraticnonlinear term ( dx/dt ) .The red lines in Figs. 3e, 7 are velocity nullclines whichare the locus of points in phase-space where the acceler-ation, dy/dt = d x/dt = 0. From Eq. 15, it follows thatthe nullcline is a quadratic function in x , x = − Γ (cid:18) dxdt (cid:19) . (21)It is interesting to note that the quadratic nullcline (seeFig. 7b) bears resemblance to a similarly curved null-cline of the BWBM with nonharmonic springs (see Fig.7d). The nullcline’s curvature is dictated by the samequadratic nonlinear term ( dx/dt ) in both cases.For comparison with this unphysical oscillator, we alsostudy a conservative Hamiltonian system and observethat it is also capable of exhibiting oscillations with sharptroughs. The simplest such system is obtained by addingcubic and quartic potential energy terms to the simpleharmonic oscillator. The former term provides an ex-plicit asymmetry in the energy functional, and the latterterm ensures stability about the equilibrium point. InAppendix VII B, this system is explored numerically andanalytically using a perturbation series constructed viathe Lindstedt-Poincar´e method [35]. B. A measure for crest - trough asymmetry
Stokes, in his study of propagating waves approachingthe shore, discussed the development of narrow crestsand wider troughs [36]. In this context measures for ve-locity and acceleration skewness have been proposed andstudied. In the folds of cerebellum, which is a stationaryspatial oscillation, a similar asymmetry presents itself inthe form of wide gyri/crest and sharp sulci/trough. Wepropose the following measure for quantifying the asym-metry in widths -crest - trough asymmetry := t crest − t trough t crest + t trough . (22)Here, the length the parameter θ traverses between con-secutive pair of points at which d x/dt = 0 on the climb-ing (falling) wave, and the immediately following falling(climbing) wave defines t crest ( t trough ), respectively (seeFig. 3d). These points are the points of intersection ofthe phase-space orbit with the velocity nullcline (see Fig.3e). The horizontal parabolic-shaped nullclines, thus in-fluence the position of these points of intersection andplay a prominent role in bringing the asymmetry betweenthe widths of the crests and troughs. For the ADO, thevelocity nullcline is exactly parabolic (see Eq. 21 andFig. 7). The parabolic shape of the nullcline is due tothe quadratic nonlinearity ( dx/dt ) in Eq. 15. For theBWBM model with nonlinear radial glial springs, thenullcline will not be exactly parabolic due to the pres-ence of other nonlinearities in Eq. 14.The measure in Eq. 22 is bounded within ( − , θ in lieu of t in Eq. 22. Studying this asymmetry forthe ADO and the nonlinear BWBM model for a rangeof tuning parameters, we see that the asymmetry scalesas 1.04 ± ± ± IV. SPATIAL CONFINEMENT
A cerebellum does not grow in isolation. It encoun-ters steric effects from the cerebrum, brain-stem, and theskull. The effect of the skull on a growing cerebrum hasearlier been studied computationally [37]. We model thesteric effects on the developing cerebellum by incorpo-rating a logistic function (1+tanh( x ))/2 into the energyfunctional, E (cid:20) r, t, dtdθ (cid:21) = (cid:90) dθ (cid:26) k r ( r − r ) − k t ( t − t ) + β (cid:18) dtdθ (cid:19) + K c q ( r − r c )]) (cid:27) . (23) FIG. 4.
Space constraint locally flattens the lobes of the linear BWBM model. a,b) Polar and cartesian representation of thesolutions to the differential equation of the linear BWBM model under a tanh(˜ r − ˜ r c ) space constraint (see Eq. 26). Theblack confining circles in (a) are the confining walls and are represented by dashed lines in (b). c) Phase-space orbits, unlikethe nonlinear BWBM model and the ADO model, remain symmetric. Parameters: c = 0 . , (cid:15) = 0 . , ρ = 15 . , ˜ K c = 30 , q =10 , ˜ t = 7 q, ˜ r = 1 q . [˜ r, d ˜ rdθ ] θ =0 is [˜ r c , The constants of the logistic function are omitted as theyvanish in the resulting Euler-Lagrange equations. Here, K c is the coupling constant and q − is the width of thestep size of the tanh function. Since K c q has the dimen-sions of k r , it is a measure of the strength of the stericinteraction. Renormalizing the parameters and variablesby q − avoids singularities in the Euler-Lagrange equa-tions in the limit of q →
0. Dividing Eq. 23 by k r q − ,we have,˜ E (cid:20) ˜ r, ˜ t, d ˜ tdθ (cid:21) = (cid:90) dθ (cid:26) (˜ r − ˜ r ) − c (˜ t − ˜ t ) + 1 ρ c (cid:18) d ˜ tdθ (cid:19) + ˜ K c tanh[˜ r − ˜ r c ] (cid:27) , (24)where ˜ E := E/ ( k r q − ), ˜ r := qr , ˜ t := qt , ˜ r c := qr c and˜ K c = K c q /k r . All parameters and variables are nowrendered dimensionless. Given that there is no explicitdependence on q , we can be assured that it won’t show upin the Euler-Lagrange equation either. The correspond-ing area conserving constraint is,12 (cid:90) dθ (˜ r − ˜ t ) = A q = dimensionless constant . (25)The Euler-Lagrange equation is { (1 − (cid:15) ) − ˜ K c sech (˜ r − ˜ r c ) tanh(˜ r − ˜ r c ) } (cid:18) d ˜ rdθ (cid:19) + ˜ K c sech (˜ r − ˜ r c ) { − (˜ r − ˜ r c ) } (cid:18) d ˜ rdθ (cid:19) + (cid:8) ρ(cid:15) c + (1 − (cid:15) ) γ (cid:9) ˜ r + ˜ K c γ sech (˜ r − ˜ r c )+ ρ(cid:15) ˜ t − γ ˜ r = 0 (26)Here γ = (cid:112) ρ ( (cid:15)c + 1). We, again, observe the ( dr/dθ ) term. However, its coefficient is zero for all r (cid:54) = r c , thus effectively localizing its effect. Results of numeri-cal integration in Fig. 4a,b shows local flattening of thecrests/gyri at contact with the confining wall. There ap-pears to be no nonlocal effects of the confining wall onthe lobes. V. A BRANCHING HIERARCHY
As the gyri/crests of the cerebellar cortex continueto grow, the sulci/troughs sharpen and become an-chored [38]. The anchoring is due to a combination ofthe radial glial cells and the pial basement membrane,a thin sheet of extracellular matrix (ECM) made up ofcollagen, laminin, and other ECM components [39]. Itis known that basement membranes stick together [40],so as the troughs sharpen, the pial membrane from eachside of the trough begins to come into contact and stick toreinforce the anchoring. The resulting anchoring centersdelineate the petal-like projections called lobules [29].We hypothesize that when the anchoring centers serveas effective boundaries between the lobules, each lobe be-comes its own subsystem. This subsystem then consistsof its own subcortex/subcore, all within the encompass-ing primary cortex/core geometry. Some of these feature-less subsystems go on to develop folds of their own to, inturn, generate another generation of subsystems and soon. See Fig. 5a. In other words, as the cerebellum con-tinues to develop, a branching hierarchy of subsystemsemerges. This branching hierarchy is yet another distin-guishing feature of the cerebellum that sets it apart fromthe cerebrum.The hierarchical generation of lobules within lobulespoints to a scale-invariant branching process. The lin-ear BWBM model, in the limit of (cid:15) →
0, also demon-strates that the number of folds N , depends on the ma-terial properties as N ≈ √ ρ and the size of the systemis rendered irrelevant. Since the formation of crests and FIG. 5.
Branching morphogenesis in the cerebellum and length-scale invariance in BWBM foliation. a) Midline saggital cross-section of mouse cerebellum 3 days post birth. Scale bar is 200 µ m. Reprinted from Fig. 7b of Ref. [22] published underthe 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 systemswith elliptical r of (d). Each lobe becomes a sub-system of its own and spawns its own folds. Here, the first and each of thesecond generation systems are generated numerically with the second generation overlaying the first after a rescaling factorof 0.28. d) Folds develop transversely to the major axis of the elliptical r (blue). Parameters: c = 0 . , (cid:15) = 0 . , ρ =15 . , ˜ t = 0 . , ˜ k r = ˜ k r and eccentricity e = 0 ,
06 respectively for the purple and green 4 lobed systems respectively. Theinitial condition at t = 0 is [˜ r, d ˜ rdθ ] = [0 . , e = 0 . troughs do not depend on system size, the BWBM frame-work is a natural description for the hierarchy. Morespecifically, as the size of the subsystems decreases witheach successive generation, crests and troughs can stillform as long as the material properties do not change.As an idealized purely geometric example of the hierar-chy, we consider an initial zeroth generation circle. Alongits perimeter, six first generation circles are generated.This process can proceed ad infinitum, to generate a frac-tal structure with a fractal dimension of log(3) / log(6 /π ).We show four generations of this hierarchy in Fig. 5b. Inthe same vein, in Fig. 5c, we represent the hierarchy offolds in the cerebellum, using the linear BWBM model.The ‘zeroth’ generation is the circular r which generatessix first generation foliations. Each lobe within consecu-tive folds is now considered as an independent sublobe.Each of these sublobes then branch into two second gen-eration lobes. As we see in Fig. 5a, the free part of thefirst generation lobe i.e. the part which is not sticking toits neighbors appears elliptical with the major axis of theellipse being parallel to the outermost edge of the lobe.Thus, we use an elliptical r while retaining the employ-ment of the linear BWBM model to generate the secondgeneration folds. See top portion of Fig. 5c. The twosecond generation lobes are generated via the method ofimages i.e. we generate a 4 lobed BWBM system as inthe top portion of Fig. 5d and use only the top half ofthis solution to represent the second generation folds inFig. 5c.It is also interesting to note that the most prominentfold occurs transversely to the horizontal major axis ofthe elliptical r (see top portion of Fig. 5d). This couldsimply be the consequence of the system trying to mini-mize its energy contribution from the radial glial springs.For comparison, we show on the bottom portion of Fig. 5d, a four lobed system generated from a circular r . VI. DISCUSSION
Inspired by cerebellar shape development, we studythe effects of nonlinear elasticity, steric confinement, anda branching hierarchy within the BWBM model. Explor-ing the effects of nonlinear elasticity of the fibrous radialglial cells, the interplay between geometry and nonlinear-ity is seen to give rise to troughs sharper than the troughsobtained in the linear BWBM model and we arrive atan asymmetry between the crests (gyri) and the troughs(sulci). This asymmetry can be understood thus: therelatively slow growth of the sub-cortex is taken into ac-count by demanding area conservation of the sub-cortexin the two-dimensional BWBM model. The associatedLagrange multiplier couples the radius of the cerebel-lum and the thickness of the cortex. Nonlinear radialsprings, in association with this coupling, results in therobust quadratic nonlinear term of the form ( dr/dθ ) inthe shape equation. To illustrate the role of the quadraticnonlinear term in sharpening the sulci, we study the sim-ple harmonic oscillator with the same form of nonlinear-ity and observe its sufficiency in achieving sharp sulci.Several other nonlinearities emerge in the shape equa-tion including a spatially-varying effective mass.The perspective of cerebellum foliation as the actionof a nonlinear oscillator can be a useful one given the ex-tensive theoretical studies of such oscillators [41, 42]. ForBWBM of the cerebellum, the linear model with constant r maps to a forced harmonic oscillator and, for smalleccentricities of r , maps to an unconventional Duffingoscillator. For nonlinear ˜ K r ( r ), we attempt to under-stand the corresponding nonlinearity in the context ofthe assisting-dampening oscillator. We hope the studyof cerebellar foliation as a nonlinear oscillator problemcontinues to be fruitful. Interestingly, the existence of anew morphological instability in confined nonlinear elas-tic sheets was found in the context of a period-doublingbifurcation, exhibiting an analogy with parametric reso-nance in another nonlinear oscillator [43].The period-doubling hierarchy found in the Ref. [43]is very different from the new hierarchy found here. Thehierarchy found here is one due to boundary conditionsin the form of anchoring centers to create sub-regions, orsub-systems, from which the same type of scale-invariantfoliation emerges, at least in the limit of small (cid:15) . Hadthe foliation mechanism not been scale-invariant in anylimit, such as with a purely elastic system, the smallersub-lobes would soon become featureless as the numberof foliations depend linearly on the perimeter of the sub-system. Within BWBM, therefore, we have identifieda new scale-invariant branching morphogenesis mecha-nism. It is not yet clear how generic this new branchinghierarchy mechanism is in terms of moving beyond thecerebellum. Ref. [19] addressed potential applications ofBWBM to two-dimensional brain organoids [44] and thedeveloping retina [45].Referring again to Fig. 5(a), one of the second genera-tion sublobes labelled L45 does not branch. Perhaps thematerial properties are altered in this sublobe so thatfeatures do not form. For sublobe L678, the two newsub-sublobes are not similar in size. This could be dueto changes in curvature of the confinement from grow-ing, surrounding tissue. So far, we have only addressedsteric, static confinement. Certainly, such variabilitiesfrom sublobe to sublobe can be explored in less idealizedconditions.We note that within the context of purely elasticitytheory, an explanation for the thicker-sulci/thinner-gyriof the developing cerebellar cortex was recently achievedby the addition of surface tension [46]. While more ener- getic contributions can certainly be added to either thepurely elastic model or to the BWBM model, the recentexperimental observations needs to be incorporated inthe modeling – the cells in the cortex are motile with cel-lular rearrangements on the time scale of minutes [22]and the cerebellum is under tension as it develops (asopposed to compression). These observations render apurely elastic model suspect. However, the differentialgrowth between the cortex and sub-cortex remains cen-tral in both classes of models. To make progress, we needfurther experimental falsification tests to rule out classesof models.Finally, given our focus on the cerebellum here, one isnaturally led to wonder whether some form of BWBM isapplicable to the cerebrum. As mentioned previously, thecerebrum and the cerebellum have rather different mor-phologies. In terms of development in the cerebellum,the predominant growth of the cells is in the cortex [38].Many such cells migrate inward to become part of thecore of the cerebellum. In the cerebrum, much of thecell proliferation is in the core and the progenitor cellsmigrate outward to become part of the cortex [47, 48] .In this sense, the two organs are inverses of each other.Given the presence of motile cells in the developing cere-brum, one may wonder whether a purely elastic approachto the developing brain is reasonable. Without a doubt,the time scales of cell migration decide the contest be-tween elasticity and fluidity. Under the framework ofliquid crystals, earlier work on the developing cerebrumarrived at a prediction for the bending modulus of thecortex [27]. This approach was based on a revised viewof the axonal tension model for the developing cerebrum[49, 50]. These novel approaches have the potential tobuild new inroads in quantitative biology.MCG acknowledges useful discussions with Manu Man-nattil. JMS acknowledges financial support from NSF-DMR-1832002 and an Isaac Newton Award from theDoD. [1] M. Ben Amar, P. Nassoy, and L. LeGoff, Philosophi-cal Transactions of the Royal Society A , 20180070(2019).[2] C. M. Nelson, Journal of Biomechanical Engineering (2016).[3] S. Budday, P. Steinmann III, and E. Kuhl, Frontiers inCellular Neuroscience , 257 (2015).[4] D. P. Richman, R. M. Stewart, J. W. Hutchinson, andV. S. Caviness, Science , 18 (1975).[5] R. Raghavan, W. Lawton, S. Ranjan, andR. Viswanathan, Journal of Theoretical Biology , 285 (1997).[6] P. Bayly, R. Okamoto, G. Xu, Y. Shi, and L. Taber,Physical Biology , 016005 (2013).[7] T. Tallinen, J. Y. Chung, J. S. Biggins, and L. Mahade-van, Proceedings of the National Academy of Sciences , 12667 (2014).[8] T. Tallinen, J. Y. Chung, F. Rousseau, N. Girard, J. Lef`evre, and L. Mahadevan, Nature Physics , 588(2016).[9] E. Hannezo, J. Prost, and J.-F. Joanny, Physical ReviewLetters , 078104 (2011).[10] A. E. Shyer, T. Tallinen, N. L. Nerurkar, Z. Wei, E. S. Gil,D. L. Kaplan, C. J. Tabin, and L. Mahadevan, Science , 212 (2013).[11] B. R. Wiggs, C. A. Hrousis, J. M. Drazen, and R. D.Kamm, Journal of Applied Physiology , 1814 (1997).[12] B. Li, Y.-P. Cao, X.-Q. Feng, and H. Gao, Journal ofthe Mechanics and Physics of Solids , 758 (2011).[13] J. W. Osborn, Journal of Theoretical Biology , 338(2008).[14] A. E. Shyer, A. R. Rodrigues, G. G. Schroeder, E. Kas-sianidou, S. Kumar, and R. M. Harland, Science ,811 (2017).[15] E. Lejeune, A. Javili, J. Weickenmeier, E. Kuhl, andC. Linder, Soft Matter , 5613 (2016). [16] E. Lejeune, B. Dortdivanlioglu, E. Kuhl, and C. Linder,Soft Matter , 2204 (2019).[17] A. Mongera and et al., Nature , 401 (2018).[18] A. Jain and et al, Nature Communications , 5604(2020).[19] T. Engstrom, T. Zhang, A. Lawton, A. Joyner, and J. M.Schwarz, Physical Review X , 041053 (2018).[20] O. Larsell, The Comparative Anatomy and Histology ofthe Cerebellum (University of Minnesota Press, 1967).[21] S. Herculano-Houzel, Front. Hum. Neurosci. , 31 (2009).[22] A. K. Lawton, T. Engstrom, D. Rohrbach, M. Omura,D. H. Turnbull, J. Mamou, T. Zhang, J. M. Schwarz,and A. L. Joyner, Elife , e45019 (2019).[23] P. Fernandez, P. A. Pullarkar, and A. Ott, Biophys. J. , 3796 (2006).[24] V. D. Varner and C. M. Nelson, Development , 2750(2014).[25] E. Hannezo, C. L. G. J. Scheele, M. Moad, N. Drogo,R. Heer, R. V. Sampongna, J. van Rheenen, and B. D.Simons, Cell , 242 (2017).[26] R. Lagrange, F. L. Jim´enez, D. Terwagne, M. Brojan,and P. Reis, Journal of the Mechanics and Physics ofSolids , 77 (2016).[27] O. Manyuhina, D. Mayett, and J. Schwarz, New Journalof Physics , 123058 (2014).[28] T. Tallinen, J. Y. Chung, J. S. Biggins, and L. Mahade-van, Proceedings of the National Academy of Sciences , 12667 (2014).[29] H. Marzban, Development of the Cerebellum from Molec-ular Aspects to Diseases (Springer, 2017).[30] C. Storm, J. J. Pastore, F. C. MacKintosh, T. C. Luben-sky, and J. P. A, Nature , 195 (2005).[31] H. T. Davis,
Introduction to Nonlinear Differential andIntegral Equations (US Government Printing Office,1961).[32] B. Ermentrout,
Simulating, Analyzing, and AnimatingDynamical Systems: A Guide to XPPAUT for Re-searchers and Students (SIAM, 2002).[33] H. Goldstein, C. Poole, and J. Safko, “Classical mechan-ics,” (2002).[34] S. H. Strogatz,
Nonlinear Dynamics and Chaos: WithApplications to Physics, Biology, Chemistry, and Engi-neering (CRC Press, 2018).[35] D. Jordan, P. Smith, P. Smith, et al. , Nonlinear Ordi-nary Differential Equations: An Introduction for Scien-tists and engineers , Vol. 10 (Oxford University Press onDemand, 2007).[36] G. G. Stokes, Transactions of the Cambridge Philosoph-ical Society (1880).[37] J. Nie, L. Guo, G. Li, C. Faraco, L. S. Miller, and T. Liu,Journal of Theoretical Biology , 467 (2010).[38] A. Sudarov and A. L. Joyner, Neural Development , 26(2007).[39] W. Halfter, S. Dong, Y.-P. Yip, M. Willem, andU. Mayer, J. Neuroscience , 6029 (2002).[40] K. Baumann, Nature Reviews Molecular Cell Biology ,767 (2014).[41] A. H. Nayfeh and D. T. Mook, Nonlinear Oscillations (John Wiley & Sons, 2008).[42] C. Hayashi,
Nonlinear Oscillations in Physical Systems (Princeton University Press, 2014).[43] F. Brau, H. Vandeparre, A. Sabbah, C. Poulard,A. Boudaoud, and P. Damman, Nature Physics , 56(2011). [44] E. Karzbrun, A. Kshirsagar, S. R. Cohen, J. H. Hanna,and O. Reiner, Nature Physics , 515 (2018).[45] H. Kolb, in Webvision: The Organization of the Retinaand Visual System , edited by H. Kolb, R. Nelson, E. Fer-nandez, and B. Jones (Webvision, 2018) Chap. 1.[46] D. Riccobelli and G. Bevilacqua, Journal of the Mechan-ics and Physics of Solids , 103745 (2020).[47] P. Rakic, Nature Reviews Neuroscience , 724 (2009).[48] Z. Moln´ar, G. J. Clowry, N. ˇSestan, A. Alzu’bi,T. Bakken, R. F. Hevner, P. S. H¨uppi, I. Kostovi´c,P. Rakic, E. Anton, et al. , Journal of Anatomy , 432(2019).[49] D. C. Van Essen, Nature , 313 (1997).[50] D. C. Van Essen, Proceedings of the National Academyof Sciences , 32868 (2020). VII. APPENDIXA. Uncoupling general Euler-Lagrangeequations for r, t
The uncoupled nonlinear differential equation in r forthe case of having k r := k r ( r ) is shown in Eq. 14.The form of such an equation is, however, dependent onthe choice of nonlinear coefficients. For general coupledEuler-Lagrange equations, we have f ( r, t ) = 0 , (27) t (cid:48)(cid:48) = g ( r, t ) . (28)Here the prime superscript indicates a derivative withrespect to θ . The derivative of Eq. 27 with respect to θ gives, ∂f ( r, t ) ∂r r (cid:48) + ∂f ( r, t ) ∂t t (cid:48) = 0 . (29)Another derivative yields, ∂f ( r, t ) ∂r r (cid:48)(cid:48) + ∂f ( r, t ) ∂t t (cid:48)(cid:48) + ∂ f ( r, t ) ∂r r (cid:48) + 2 ∂ f ( r, t ) ∂r∂t r (cid:48) t (cid:48) + ∂ f ( r, t ) ∂t t (cid:48) = 0 . (30)Eq. 27,29 can be solved to find t := t ( r ) , t (cid:48) = t (cid:48) ( r, r (cid:48) ).The uncoupled differential equation in r would then be, ∂f ( r, t ( r )) ∂r r (cid:48)(cid:48) + ∂f ( r, t ( r )) ∂t t (cid:48)(cid:48) ( r ) + ∂ f ( r, t ( r )) ∂r r (cid:48) + 2 ∂ f ( r, t ( r )) ∂r∂t r (cid:48) t (cid:48) ( r ) + ∂ f ( r, t ( r )) ∂t t (cid:48) ( r ) = 0 . (31)The last two terms for Eq. 14 is absent and this particularcase of nonlinear differential equation is of the form, ∂f ( r, t ( r )) ∂r r (cid:48)(cid:48) + ∂ f ( r, t ( r )) ∂r r (cid:48) + ∂f ( r, t ( r )) ∂t g ( r, t ( r )) = 0 , (32)1where the implicit dependence - t ( r ) is to be taken intoaccount after evaluating the partial derivatives in t, r .This renders an uncoupled differential equation of motionfor r . B. Higher order corrections to the SimpleHarmonic Oscillator
With no loss of generality, we can study the simpleharmonic oscillator with cubic and quartic corrections tothe potential energy as,
FIG. 6.
Lindstedt-Poincar´e perturbation method for the cubicand quartic energy corrections to the SHO. O ( (cid:15) ) curve followsthe leading edge of the numerical solution but doesn’t capturethe offset in troughs and the relative difference in widths ofcrests and troughs. Here (cid:15) = 0 . , Γ (cid:48) = 0 . d xdt = − x − x − Γ (cid:48) x , (33)a single parameter equation where all variables are nondi-mensionalized. To evolve a perturbative scheme, we treatthe last two terms of the above equation as a perturba-tion to the simple harmonic oscillator. To that effect,with (cid:15) > d xdt = − x − (cid:15) ( x − Γ (cid:48) x ) . (34)The Lindstedt-Poincar´e method [35] introduces angu-lar frequency ω by change of variable τ = ωt . This allowsfor coupling between the frequency and amplitude. Eq.34 is now written as, ω d xdτ + x + (cid:15)x − (cid:15) Γ (cid:48) x = 0 . (35)The series expansions for x and ω are, x ( t ) = x ( t ) + (cid:15)x ( t ) + (cid:15) x ( t ) + ...ω = 1 + (cid:15)ω + (cid:15) ω + ... (36)The boundary conditions are, x (0) = A , ˙ x (0) = 0 ,x i (0) = 0 , ˙ x i (0) = 0 i > . (37) O ( (cid:15) ) We substitute the series expansions of Eq. 36 in Eq. 35.All the terms at every given order of (cid:15) must add up tozero for Eq. 35 to hold. The response equation at thezeroth order of (cid:15) then is, d x dτ + x = 0 , (38)which is just the equation of motion of the simple har-monic oscillator. The zeroth order response that obeysthe boundary conditions 37 is, x ( τ ) = A cos( τ ) . (39) O ( (cid:15) ) At the first order of (cid:15) we have, d x dτ + x = x + Γ (cid:48) x − ω d x dτ . (40)Solutions of lower order response equations are used tofind solutions of higher order response equations. ω isfixed by the demand that the solution to Eq. 40 be pe-riodic. For this, upon substituting the zeroth order solu-tion (see Eq. 39) in Eq. 40 we set the coefficient of cos( τ )in the resulting RHS of Eq. 40 to vanish. We obtain ω = − (cid:48) a . (41)The amplitude-frequency dependence can be seen here.The first order correction to x which obeys the boundaryconditions is, x ( τ ) = A (cid:18) Γ (cid:48) A −
13 cos( τ ) (cid:19) + A − A τ ) − Γ (cid:48) A
32 cos(3 τ ) . (42) O ( (cid:15) ) At the second order of (cid:15) , we have, d x dτ + x = 3Γ (cid:48) x x +2 x x − ( ω +2 ω ) d x dτ − ω d x dτ . (43)Carrying out the same procedure at this order, we have, ω = 1384 (cid:18) − a Γ (cid:48) + 144 a Γ (cid:48) − a Γ (cid:48) ω − a + 128 a ω − ω (cid:19) , (44)and x ( τ ) = − a [ − τ )(3 a Γ (cid:48) − a − ω )+ a (cos(3 τ )(9 a Γ (cid:48) + 96 a Γ (cid:48) + 216Γ (cid:48) ω + 64)+ 60 a Γ (cid:48) cos(4 τ ) + 9 a Γ (cid:48) cos(5 τ ) − a Γ (cid:48) + 128)] . (45)2 FIG. 7.
Phase portraits of BWBM and ADO systems . Blue lines are position nullclines: drdt = 0 and red lines are velocitynullclines: d rdt = 0. The tangent to the flow lines are vertical when they pass the position nullcline and horizontal whenthey pass the velocity nullcline. a)Linear BWBM b) ADO c) Cubic and quartic energy corrections to the SHO d) BWBMwith nonharmonic springs. In (b,c,d) left-right asymmetric orbits where there is a steeper fall on the left translates to sharpertroughs in position-time space. We see left-right asymmetric orbits with or without bent nullclines (see b,d & c). The perturbative expansions are compared with numer-ical solution in Fig. 6.
C. Phase portraits
Closed orbits in phase-space translate to periodic mo-tion in the position-time space. Plotting orbits in phasespace for different initial conditions builds the phase por-trait of the system. The x ( dx/dt ) nullcline are the locus of points where dx/dt = 0 ( d x/dt = 0). The nullclinesare generally used to divide the phase-portrait into re-gions where the tangent of the orbit points in the samegeneral direction (NW, NE, SE or SW). The orbit has avertical (horizontal) tangent when it passes through the x ( dx/dθdx/dθ