Brain folding is initiated by mechanical constraints without a cellular pre-pattern

During human brain development the cerebellum and cerebral cortex fold into robust patterns that increase and compartmentalize neural circuits. Although differential expansion of elastic materials has been proposed to explain brain folding, the cellular and physical processes responsible at the time of folding have not been defined. Here we used the murine cerebellum, with 8-10 folds, as a tractable model to study brain folding. At folding initiation we considered the cerebellum as a bilayer system with a fluidlike outer layer of proliferating precursors and an incompressible core. We discovered that there is no obvious cellular pre-pattern for folding, since when folding initiates, the precursors within the outer layer have uniform sizes, shapes and proliferation, as well as a distribution of glial fibers. Furthermore, although differential expansion is created by the outer layer expanding faster than the core at folding initiation, thickness variations arise in the outer layer that are inconsistent with elastic material models. A multiphase model was applied that includes radial and circumferential tension and mechanical constraints derived from in vivo measurements. Our results demonstrate that cerebellar folding emerges from mechanical forces generated by uniform cell behaviors. We discuss how our findings apply to human cerebral cortex folding.


Introduction 30
Recent work in brain folding has primarily focused on the cerebral cortex and involved models of 31 differential expansion of elastic materials to drive mechanical instabilities that lead to folding 1-7 . Current 32 models are able to create three dimensional shapes strikingly similar to the folds seen in the human 33 cortex 8 . However, the cell and tissue level mechanics actually present at the initiation of folding have not 34 been considered or defined. The simple alignment of folds along the anterior-posterior axis of the murine 35 cerebellum and the genetic tools available in mouse allow for precise developmental interrogation, to 36 identify and analyze the cellular behaviors driving growth that could create differential expansion, and to 37 reveal additional critical forces. 38

39
The developing cerebellum is distinct from the cerebral cortex, as it has an external granule cell layer 40 (EGL) of proliferating precursors that covers the surface and generates growth primary in the anterior-41 posterior (AP) direction 9-11 . A bilayer system is therefore an excellent starting approximation for 42 cerebellar folding. Here we show that cerebellar folding emerges from differential expansion between an 43 un-patterned uniformly expanding fluid-like EGL and an incompressible underlying core. Additionally, 44 we demonstrate that thickness variation in the EGL at the base of the forming fissures, termed anchoring 45 centers (AC) 12 , are inconsistent with traditional elastic wrinkling models driven by differential growth. 46 We constrained a recent model 13 that takes into consideration that tension could play an import role in the 47 developing brain and used our in depth developmental data to capture the initiation of cerebellar folding. 48 49

No cellular pre-pattern at initiation of folding 50
To systematically determine the cellular behaviors underlying the initiation of cerebellar folding, we first 51 tested whether regional differences in EGL proliferation rates are present that could influence the folding 52 pattern of the cerebellum, since the EGL drives the majority of cerebellar growth 9-11 . Proliferation rates (S 53 phase index) were measured in the EGL during folding initiation (E16.5 and E17.5) in the inbred FVB/N 54 strain to reduce variation between samples. First we asked if the regions that will give rise to distinct sets 55 geometry, with the stiffness of the external layer defined as E o, the stiffness of the core as E i , and the 133 thickness of the external layer denoted as t , the folding wavelength λ is given by 27 134 If the length of the system is , then the number of folds is 135 In other words, the number of folds is inversely proportional to the thickness of the EGL. 136 137 To more closely model the geometry of the cerebellum, we explored a standard elastic bilayer model in a 138 circular geometry using the observed ratio of thickness of the EGL to radius of the cerebellum near the 139 onset of shape change (smooth to wrinkled after E16.5) and invoking a neo-Hookian elastic solid for both 140 layers 28 . Previous work applied a tri-layer elastic model to the cerebellum, incorporating the molecular 141 layer 7 . We have not included the molecular layer in our model since it is not present when folding 142 initiates. The resulting shape change was studied as a function of the ratio of the layer stiffness values 143 (Fig. 4a). We found that to induce the observed number of folds at initiation of folding through wrinkling 144 based models constrained by our measurements of the embryonic cerebellum, a large stiffness ratio was 145 required of around 50. To map the stiffness contrast in the cerebellum we used scanning acoustic 146 microscopy to measure the bulk modulus of the cerebellum daily from E16.5 to P18.5 ( Fig. 4b-c, 147 Supplemental Fig. 2). For small deformations, we expect the instantaneous bulk modulus to be linearly 148 related to the stiffness and, therefore, the ratio of the instantaneous bulk moduli should scale similarly to 149 the ratio of stiffnesses (assuming the same Poisson's ratio). We found that although the EGL has a higher 150 instantaneous bulk modulus than the core at all stages measured, the ratio was not sufficient to produce a 151 folding wavelength similar to that in the cerebellum (Fig. 4d). Small modulus contrasts have been 152 reported for other brain regions with other test modes. 29 153 154 Wrinkling models predict compressive forces in the outer layer. Consistent with this prediction, 155 simulations of cuts through the outer layer and into the inner layer predict that the outer layer should not 156 open (Fig. 4e). Using surgical dissection blades we made radial cuts across the meninges, EGL, and into 157 the core of live E16.5 tissue slices. Time-lapse imaging revealed that unlike the prediction, the EGL opens 158 as well as part of the underlying cut in the core (Fig. 4f-h, Supplemental Fig. 4a-c, and Supplemental 159 Movie 1). This result is consistent with there being circumferential tension within the outer layers of the 160 cerebellum. 161

162
An additional aspect of the model that does not fit the biology of the developing cerebellum is that it 163 requires the EGL to be thinnest at the base of each AC, which are the lowest parts of the cerebellar 164 surface, and thus to have an "in-phase" thickness variation 13 . Without this feature, a purely elastic model 165 cannot be in mechanical equilibrium (in the quasistatic limit). However, we previously noted that the EGL 166 is thickest in the ACs when folding initiates, i.e., it has an "out-of-phase" thickness variation 12 . To 167 validate this observation, we quantified the thickness of the EGL and found it to be around 1.2 -1.4 times 168 thicker in the ACs than in the surrounding EGL at E16.5 and 17.5 when folding initiates (Fig. 5a-d and  169 Supplemental Fig. 3), and the thickness ratio increased to 1.7 times at E18.5 (Fig 5d). It is of interest to 170 note that the final thicknesses of the layers of the cortex of the mature cerebellum, just as in the cerebral 171 cortex, are in-phase. These results further show that traditional wrinkling models cannot capture the 172 initiation of cerebellum folding, and highlight the importance of making measurements at the time of 173 folding rather than when it is complete. 174 175

A multiphase model approximates folding 176
We recently developed a model for folding from differentially expanding bi-layer brain tissues that takes 177 into account the possible contribution of mechanical constraints present in the developing brain 13 . We 178 made five primary assumptions based on the measurements presented above. First, the core is an 179 incompressible material (µ), as indicated by our instantaneous bulk modulus measurements. Second we 180 utilize a uniform expansion of the EGL (k t ) as indicated by the proliferation rate. Third, we do not assume 181 that the EGL is an elastic material. Forth, we assume an elastic component radially and circumferentially 182 to the entire cerebellum (k r ), possibly mediated by fibers spanning the cerebellum including those of the 183 Bergmann glia and radial glia as well as the pial surface and meninges covering the cerebellum. And fifth, 184 we posit that the EGL is constrained towards a uniform thickness by fibers spanning the EGL (β), 185 possibly the Bergmann glia. Given the interplay between incompressible material, compressible fibrous 186 material, and a proliferating fluid-like EGL, this model is multiphase. 187

188
We constructed an energy functional parameterized by both the inner and outer boundary of the EGL and 189 incorporating the above five assumptions into three dimensionless parameters (µ/k r , k r /k t , k t /β). We also 190 constrained the parameters with our developmental data. Minimization of the energy functional yields an 191 equation for a driven harmonic oscillator yielding sinusoidal shapes for both the inner and outer boundary 192 of the EGL. In contrast with the elastic bilayer wrinkling model, EGL thickness oscillations are found to 193 be out-of-phase with the surface height (radius) oscillations when 0 < µ/k r < 1. The ratio of these two 194 quantities, the measured surface height amplitude (A r ) and the EGL thickness amplitude (A t ) is given by 195 which need not be ≫ 1 as is typical of elastic bilayer wrinkling, and the number of initial folds at E16.5 196 is determined by 197 Note that in contrast with elastic wrinkling, the number of initial folds does not depend on the thickness 199 To test the model, we used the observed EGL thickness amplitude, average EGL thickness, average 202 cerebellum radius, and the number of initial folds at E16.5 to constrain 3 of the 5 model parameters. A 4th 203 parameter (u/k r , denoted as ε) was assumed to scale linearly with time and thus generate predictions of 204 cerebellum shapes at later developmental stages (E17.5 and E18.5). We indeed found that the model well-205 approximates the phase and amplitude behavior of EGL thickness and radius oscillations, during these 206 stages. (Fig. 5e-g, Supplemental Fig. 6). 207 208 As the model requires radial tension in addition to the circumferential tension demonstrated above, we 209 examined evidence of radial tension between the EGL and the ventricular zone (VZ) at the initiation of 210 folding. Horizontal cuts were made across the core of live E16.5 tissue slices below the EGL and above 211 the VZ that would cut across anterior radial fibers (Fig 5h). As predicted, after cutting, the tissue relaxed 212 revealing tension directed radially within the cerebellum ( The cerebellum has hierarchical folding in which the initial folds become subdivided. Given that ACs 219 hold their position during development and compartmentalize granule cells within lobules of the EGL 10 220 we reasoned that the ACs could be acting as mechanical boundaries enabling similar mechanics to drive 221 the secondary folding. To test this possibility we measured the expansion of the EGL and the core of the 222 individual lobule regions from E18.5 to P3. We found that indeed in the lobule regions that undergo 223 folding there is a temporal correlation between the onset of sub-folding and differential expansion occurs 224 ( Fig. 6a-d). In contrast, the region (L45) that does not fold during the same time period has a different, 225 more rectangular shape, and the ratio of EGL growth to core growth is proportional for a rectangle during 226 the time measured (Supplemental Fig. 5). We propose that ACs create mechanical boundaries that divide 227 the cerebellum into fractal-like domains with similar physical properties to the initial unfolded 228 cerebellum. Additionally, since ACs compartmentalize granule cells within the EGL lobule regions 10 , 229 once separated the lobule regions can develop distinct characteristics, like the observed differential 230 proliferation rates. We speculate, therefore, that the folding patterns seen across cerebella in different 231 species evolved by adjustment of global as well as regional levels of differential expansion and tension 232 which ultimately mold the functionality of the cerebellum. 233 234

Discussion 235
Here we have provided the first evidence that brain folding emerges from a uniform expansion of cells 236 without obvious pre-patterning. Thus, traditional morphometric cellular behaviors such as changes in cell 237 shape, size and proliferation do not direct where cerebellar folding initiates. Furthermore, our 238 developmental interrogation revealed that thickness variations arise in the EGL that are fundamentally 239 inconsistent with traditional elastic bilayer wrinkling models. By applying a novel multiphase model 240 based on the cellular structure of the cerebellum at the time of folding initiation we were able to capture 241 the shape variations with the correct number of folds. The model accounts for: 1) a rapidly expanding 242 fluid-like EGL, whose thickness is regulated by Bergmann glial fibers, 2) a slower growing 243 incompressible core, and 3) fibrous material in the form of glial fibers (and possibly axons) as well as the 244 meninges that provide radial and circumferential tension (Supplemental Fig. 6). One prediction of the 245 model is that adjusting the amount of tension spanning the cerebellum will change the degree of folding. 246 Indeed, alterations of the cells that likely create tension-based forces could explain the dramatically 247 disrupted folding seen in mouse mutants in which radial glia do not produce Bergmann glia 30 . Without 248 Bergmann glia, the EGL would be expected to not form a layer with regular thickness and it should be 249 more sensitive to variations in radial glial tension. Consistent with this prediction, mutants without 250 Bergmann glia have more localized and less regular folds 30 . Our combination of experimental studies and 251 modeling thus provide new insights into brain folding, including an underappreciated role for tension. 252 253 For our multiphase model to predict the observed shape changes in the murine cerebellum from E17.5 to 254 E18.5 the ratio of the core stiffness over the radial tension must increase. As the measured bulk modulus 255 of the core shows no increase during development, our model predicts that the radial tension must 256 decrease during development. While the cerebellum is crossed by many fibers at folding initiation, radial 257 glial fibers are an attractive candidate to mediate this change in radial tension 31,32 . First, they span from 258 the VZ to the surface of the cerebellum at E16.5. Additionally, during folding initiation the radial glia 259 undergo a transition into Bergmann glia where they release their basal connection to the VZ and the cell 260 body migrates towards the surface 19 . This transition could lead to a reduction in the global radial tension 261 and thus would be consistent with our model prediction.

Animals. 281
The inbred FVB/N stain was used for all proliferation rate, area, length, and expansion rate 282 measurements. Atoh1-CreER 40 , Nestin-CreER 41 , Rosa26 MTMG,42 , were used to quantify cell shape and size 283 as well as fiber distribution and were maintained on the outbred Swiss Webster background. The Swiss 284 Webster strain was used for scanning acoustic microscopy. Both sexes were used for the analysis. For embryonic stages heads were fixed in 4% paraformaldehyde overnight at 4°C. For postnatal animals, 298 the brain was dissected out first before fixation. Tissues were stored in 30% sucrose. For all proliferation, 299 area, length, and thickness measurements brains were embedded in optimal cutting temperature (OCT) 300 compound. Parasagittal sections were collect with a Leica cryostat (CM3050s) at 10µm. 301 302 Prior to IHC, EdU was detected using a commercial kit (Invitrogen, C10340). Following EdU reaction the 303 following primary antibodies were used either overnight at 4°C or 4 hours at room temperature: mouse 304 anti-P27 (BD Pharmingen, 610241), rabbit anti-GFP (Life Technologies, A11122), rat anti-GFP (Nacalai 305 Tesque, 04404-84). All antibodies were diluted to 1:500 in 2% milk (American Bioanalytical) and 0.2% 306 Triton X-100 (Fisher Scientific). Alexa Fluor secondary antibodies (1:500; Invitrogen) were used: Alexa 307 Fluor 488 donkey anti-rabbit, A21206, Alexa Fluor 488 donkey anti-rat, A21208, Alexa Fluor 488 308 donkey anti-mouse, A21202, Alexa Fluor 647 donkey anti-mouse, A31571. EdU was detected using a 309 commercial kit (Invitrogen, C10340). parafinized, hydrated in de-ionized water and scanned on the SAM to generate maps of amplitude, sample 378 thickness, speed of sound, acoustic impedance, attenuation, bulk modulus, and mass density. Co-379 registered histology and SAM amplitude images were used to identify regions-of-interest (ROIs) 380 corresponding to the EGL layer and underlying core of the cerebellum in each sample. Bulk modulus was 381 analyzed as a measure of tissue stiffness: ROI measurements were acquired from 3 sections from 3 382 embryos at each developmental stage. 383

Finite element simulations 385
The wrinkle of a circular bilayer structure in Fig. 3a was simulated with commercial software ABAQUS. 386 Both film and substrate were modeled as incompressible neo-Hookean materials. The ratio between shear 387 moduli of the film and substrate was 50 and the initial radius of the simulated structure was 16 times that 388 of the film thickness. The differential growth of the EGL and core was modeled by an isotropic expansion 389 of the film in the bilayer structure. 390

391
To test the elastic wrinkling model, we conducted finite element (FE) simulations for bilayer structures 392 with a film bonded on a substrate, which represents the EGL layer and core structure, respectively. The 393 structures were assumed to be under 2D plane strain deformation to mimic the quasi-2D nature of 394 cerebellum wrinkles. Neo-Hookean model was adopted to describe the elastic properties of both film and 395 substrate, whose strain energy can be expressed as 396 where is the shear modulus and ! represents the first invariant of the deformation gradient tensor . 397 The Poisson's ratios for the film and substrate were set to be 0.5, based on experimental observations that 398 the bulk modulus of EGL and core are in the order of GPa, much larger than the shear modulus of soft 399 tissues (~ kPa). 400

401
We carried out FE simulations through commercial software ABAQUS. A second order 6 node hybrid 402 element (CPE6MH) was utilized to discretize the film and substrate. Very fine FE meshes were used to 403 make sure the results independent of mesh size. To incorporate differential growth in real EGL layer and 404 core, an isotropic growth deformation tension was applied to the modeled film by decoupling the 405 deformation tenor into elastic deformation part and growth part G. 406 = • For simplicity, we assume the growth part is isotropic and controlled by a scalar variable 407

Methods for modeling initiation of cerebellar folds 417
We considered the midsagittal section of the cerebellum and, therefore, formulated a two-dimensional 418 model. The distance of the outer edge of the EGL and, hence, the outer edge of the cerebellum from the 419 center of the cerebellum was defined as r(θ) with θ as the angular coordinate. We assumed that r(θ) was 420 single-valued. The thickness of the EGL was defined as t(θ). See model schematic below. 421

422
Taking into account the four assumptions discussed in the main text, we constructed the following energy 423 functional to be minimized 424 , , with k r as the stiffness modulus (a spring constant in one-dimension) of the radial glial fibers and the pial 427 surface contained in the meninges surrounding the cerebellum since the cerebellar radius is proportional 428 to its perimeter, r 0 as the preferred radius of the cerebellum, k t denoting a growth potential due to cell 429 proliferation, t 0 as thickness of the EGL (cortex), and, β quantified the mechanical resistance to changing 430 the thickness of the EGL. Given our first assumption of an incompressible cerebellar core, we imposed The variational analysis yielded the following equation of shape for t(θ); 444 . The solution to the equation of shape was 448 There was 451 an additional equation of shape for r(θ) from the variational principle that depended on t(θ) and so was 452 We used the measured data at E16.5 to set the parameters to make predictions for the shape of both the 455 EGL and core (and so the relationship between the two) at later times. Plots assumed a circular preferred 456 shape, and with other parameters as follows: = ! is shown in Fig. 5f

Statistical analyses 470
Statistical analyses were performed using Matlab software. Significance was determined at P< 0.05. Two-471 way ANOVA was used for proliferation analysis as two variables were tracked, mouse and region. Cell 472 shape, volume, fiber distribution, EGL thickness and bulk modulus were run under a standard ANOVA. 473 After ANOVA analysis a multiple comparison was run with Tukey's honestly significant difference 474 criterion. F-test for variance and two-tailed student's paired t-test were used for slice cutting and 475 relaxation quantifications. The degrees of freedom, test statistics, and P values, are given in the figure  476 legends, where appropriate. All error bars are standard deviations. No statistical methods were used to 477 predetermine the sample sizes. We used sample sizes aligned with the standard in the field. No 478 randomization was used nor was data collection or analysis performed blind.      Fig. 6: Progressive folding occurs during regional differential expansion. a,b, H&E stained midline sagittal sections of FVB/N cerebella at P1 and P3. c, Expansion of lobule length and lobule area approximate the proportional expansion of a semi-circle (curve) at E18.5 and P0. After P0 the EGL expansion in both regions increases more than the underlying area creating differential expansion. d, Folding initiates during regional differential expansion. Scale bars: 200 µm. Only E16.5 cerebella that showed regional thickening in the geometry where ACs normally arise were used for the measurements, and one embryo did not yet have an AC3. a-c, Thickness variation in and surrounding AC1 (anova a, df = 17, P = 0.13, F = 1.55 b, df = 29, P = 7.0e -14 F = 6.82 c, df = 57, P = 9.1e -11 F = 4.05). d-f, Thickness variation in and surrounding AC2 (anova d, df = 17, P = 0.08 F = 1.74 e, df = 29 P = 3.9 e-22 F = 11.88 f, df = 57 P = 2.9e -35 F = 16.35). g-i, Thickness variation in and surrounding AC3 (anova g, df = 17, P = 0.59 F = 0.89 h, df = 29, P = 2.4e -17 F = 9.81 i, df = 57, P = 7.6e -33 F = 14.57 Supplemental Fig. 4: The stress patterns within the cerebellum are different between the EGL and the VZ. a-c, Example of a live cerebellar slice before a and after b a radial cut through the EGL, and still images from a time-lapse c. Time = 0 minutes is at the time it takes to remove the knife and start the imaging, therefore the cut has already begun opening. d-f, Example of a live cerebellar slice before d and after e a horizontal cut through between the EGL and ventricular zone (VZ), and still images from a movie f. g-i, radial cuts through the EGL open more quickly initially than horizontal cuts between the EGL and the Ventricular zone, but the latter continue to relax for longer (g, f-test for unequal varience P = 0.09, two tailed t-test df = 16, p = 0.03, T = -2.43; h, f-test P = 0.04 and unequal varience two-tailed t-test df = 12.8 P = 0.16, T = -1.48; i, f-test P = 0.49 and two tailed t-test df = 16, P = 0.03, T = 2.43). j The degree of opening is tightly related to the length of the opening in horizontal cuts but not in radial cuts (f-test P = 0.02, unequal varience two tailed t-test df = 11.89, P = 0.02, T = -2.80). Stars: statistical differences. Error bars: S.D.