The role of mechanics during brain development.

Convolutions are a classical hallmark of most mammalian brains. Brain surface morphology is often associated with intelligence and closely correlated to neurological dysfunction. Yet, we know surprisingly little about the underlying mechanisms of cortical folding. Here we identify the role of the key anatomic players during the folding process: cortical thickness, stiffness, and growth. To establish estimates for the critical time, pressure, and the wavelength at the onset of folding, we derive an analytical model using the Föppl-von-Kármán theory. Analytical modeling provides a quick first insight into the critical conditions at the onset of folding, yet it fails to predict the evolution of complex instability patterns in the post-critical regime. To predict realistic surface morphologies, we establish a computational model using the continuum theory of finite growth. Computational modeling not only confirms our analytical estimates, but is also capable of predicting the formation of complex surface morphologies with asymmetric patterns and secondary folds. Taken together, our analytical and computational models explain why larger mammalian brains tend to be more convoluted than smaller brains. Both models provide mechanistic interpretations of the classical malformations of lissencephaly and polymicrogyria. Understanding the process of cortical folding in the mammalian brain has direct implications on the diagnostics of neurological disorders including severe retardation, epilepsy, schizophrenia, and autism.


Introduction
For more than a century, the unique surface morphology of the mammalian brain has fascinated scientists across all disciplines (Le Gros Clark, 1945): why does the brain have this complex convoluted structure, and, more importantly, to which extent is brain structure correlated with brain function (Welker, 1990)? From a mechanics point of view, these questions naturally translate into the quest for a basic understanding of brain morphology (Bayly et al., 2014): what are the underlying mechanisms of brain folding?
The mammalian brain is composed of an outer cortical layer of gray matter, consisting primarily of cell bodies, and an inner subcortical core of white matter, consisting primarily of axons. Within the limited space inside the skull, gyrification, the folding of the cortical layer, is viewed as a process to maximize the number of cell bodies and minimize the distance between them (Zilles et al., 2013). The total number of neurons, the number of connections, and the signaling speed are directly correlated with the capacity of information processing. It is thus not surprising that not only the total brain volume, but also the brain surface area, is viewed as strong indicators of intelligence (Roth and Dicke, 2005). The ratio between brain surface area and brain volume, and with it the degree of gyrification, can vary significantly between species (Hofman, 1989). With 86 billion neurons, 0.15 quadrillion connections, and a mass of 1500 g, the human brain is often considered as the most developed mammalian brain (Herculano-Houzel, 2009). Fig. 1 illustrates the average surface-to-volume relation of shrew, hedgehog, cat, monkey, human, and dolphin brains (Hofman, 1989). The graph confirms our intuition that larger mammals generally tend to have larger brains (Welker et al.). However, surprisingly, the surface area of the mammalian brain increases disproportionally faster than its volume (Roth and Dicke, 2005). In the log-log plot, the surface-to-volume ratio scales with a slope of 0.9, which is significantly larger than the slope of 2/3 for isometric scaling: larger mammals not only have larger, but also more convoluted brains. This naturally raises the question (Geschwind and Rakic, 2013): what is the evolutionary advantage of a folded brain? Fig. 2 shows three explanted mammalian brains, which reconfirm that larger mammals tend to have larger brains (Sun and Hevner, 2014): the cow brain is larger than the pig brain, which in turn is larger than the brain of the sheep. Frontal coronal sections of the three brains illustrate that the degree of gyrification increases with brain size: the surface of the cow brain seems more folded than the surface of the pig brain, which seems more folded than the brain of the sheep (Zilles et al., 2013). The cortical thickness, however, seems to vary only marginally between the different species (Welker, 1990). Studies of mammalian brains indicate that brain size is not the only contributing factor to gyrification (Zilles et al., 1988). This motivates our hypothesis that rather than brain size, other anatomical features like cortical thickness, cortical stiffness, and cortical growth during brain development play a crucial role in pattern formation.
The development of the mammalian brain takes place in two distinct stages, which are crucial for cortical folding (Roth and Dicke, 2005): first, progenitor cells located around the ventricles divide symmetrically into two new progenitor cells to increase the total number of cells. Then, these newly created cells divide asymmetrically into a progenitor cell and a neuron (Sun and Hevner, 2014). Neurons migrate toward the surface along radial glial cells to form the cortical layer (Hatten, 1999). According to the radial unit hypothesis, all neurons of the same progenitor cell stack up on top of one another to form a cortical column (Rakic, 1988). Symmetric division is therefore primarily responsible for an increase in brain surface area, while asymmetric division is primarily responsible for an increase in cortical thickness (Roth and Dicke, 2005). Irregularities in cell division or cell migration can evoke abnormalities in surface area or thickness (Hatten, 1999). Those range from polymicrogyria, a malformation associated with an increased surface area and an excessive number of small folds (Tortori-Donati et al., 2005), to lissencephaly, a malformation associated with an increased cortical thickness and a reduced number of shallow folds (Landrieu et al., 1998). Severe malformations are often correlated with neurological disorders, including developmental delay, epilepsy, schizophrenia, and autism (Raybaud and Widjaja, 2011). Despite its pathophysiological importance, the phenomenon of cortical folding remains barely understood (Ronan et al., 2014).
Cytological studies alone fail to explain the process of cortical folding (Schwartzkroin and Walsh, 2000) and seem to indicate that mechanical factors could play a crucial role . Two competing hypotheses suggest that cortical folding is either driven by axonal tension (Van Essen, 1997) or by differential growth (Richman et al., 1975). There is no direct evidence for either of these theories. Axonal tension, a mechanism to bring functionally related units topographically closer Fig. 1. Surface-to-volume ratio of the mammalian brain. Larger mammals have larger brains (Welker at al.). The dashed line with a slope of 2/3 indicates isometric scaling for which brain surface area would scale proportionally with brain volume. The solid line with a slope of 0.9 indicates that the mammalian brain surface area increases disproportionally faster than brain volume (Hofman, 1989). The degree of gyrification increases with brain size. together, can explain folding irrespective of tissue stiffness, but disagrees with dissection experiments, which reveal significant tangential tension in the outer layers but not inside the developing gyri (Xu et al., 2010). Differential growth, a mechanism to release residual stresses by surface buckling, agrees well with dissection experiments, but relies on unrealistic stiffness differences between cortex and subcortex (Bayly et al., 2014). From other biological systems we know that differential growth is capable of generating sufficient compressive stresses to induce structural instabilities (Moulton and Goriely, 2011). Different geometric constraints, stiffness ratios, and growth rates may evoke different types of instabilities like buckles, wrinkles, creases, or folds (Li et al., 2012). Soft materials are especially susceptible to surface folding because of their low material stiffness (Li et al., 2011). Unfortunately, experiments to characterize the stiffness of the brain are rare and measured stiffness values span several orders of magnitude (Franceschini, 2006). Virtual experiments and systematic parameter studies provide a powerful alternative to explore the developing mammalian brain.
Here we model brain development using the continuum theory of finite growth (Rodriguez et al., 1994). We model the cortex as a morphogenetically growing outer layer of cell bodies and the subcortex as a strain-driven growing inner core of axons. Motivated by experiments of axon elongation (Bray, 1984), we assume that chronic axonal overstretch activates mechanotransduction pathways, which collectively result in a gradual increase in axonal length (Dennerll et al., 1989). This approach combines the two popular hypotheses of cortical folding, axonal tension and differential growth . Using this model, we explore the effect of three key players on cortical folding: cortical thickness, stiffness, and growth. Since the absolute cortical stiffness and growth rate are poorly characterized , we explore the role of the relative cortical stiffness, and growth rate with respect to the subcortical properties.
The remainder of this paper is organized as follows: in Section 2, we present our analytical model for cortical folding to establish analytical estimates for the critical time, the critical pressure, and the critical wavelength at the onset of folding. In Section 3, we introduce our continuum model for finite growth to predict brain surface morphologies beyond the onset of folding. In Section 4, we summarize its computational solution within a nonlinear finite element setting. In Section 5, we utilize our model and perform systematic sensitivity studies of cortical thickness, stiffness and growth to understand the origin of pathological malformations. We close by a critical discussion of our results and their potential impact on understanding brain development in Section 6.

Analytical model
To establish analytical estimates for the brain surface morphology, we approximate cortical folding as the instability problem of a confined, layered medium subjected to growth-induced compression. We adopt the Föppl-von Kármán theory (Föppl, 1907;von Kármán, 1910), and model the cortical deflection w using the classical fourth order plate equation (Dervaux et al., 2009), 1cm sheep sheep pig pig cow cow 1cm 1cm 1cm 1cm 1cm Fig. 2. Surface morphology of the mammalian brain. Larger mammals have larger brains: the cow brain, right, is larger than the pig brain, middle, which is larger than the sheep brain, left. Photographs of the entire brain, upper row, and a frontal coronal brain section, lower row, illustrate the characteristic folding pattern: the surface of the cow brain, right, is more folded than the surface of the pig brain, middle, which is more folded than the sheep brain, left. The degree of gyrification increases with brain size. Fig. 3 illustrates the analytical model with the cortical deflection w, Young's modulus of the cortex E c , Poisson's ratio of the cortex ν c , the cortical thickness t c , the cortical pressure P, and the deflection-induced transverse force of the subcortical foundation q. We adopt a sinusoidal ansatz for the cortical deflection, where the amplitude w 0 represents the sulcal depth, the wavenumber n represents the number of gyri and sulci, and the wavelength λ crit represents the distance between two neighboring gyri. With this ansatz, the fourth order equation for cortical folding (1) takes the general form (Biot, 1957), where the transverse force q depends on the nature of the subcortical foundation. In the following, we illustrate analytical estimates for the critical pressure P crit and the critical wavelength λ crit for both an elastic foundation (Biot, 1937) as shown in Fig. 3, left, and a growing foundation (Biot, 1957) as shown in Fig. 3, right.

Growing cortex on elastic subcortex
To establish analytical estimates for cortical folding on an elastic subcortical foundation, we interpret the subcortex as an infinite half-space and impose a sinusoidal deflection wðxÞ ¼ w 0 cos ðnxÞ with a wavelength n and an amplitude w 0 on the upper boundary. This deflection evokes a transverse force qðxÞ with the same wavelength n, Its amplitude q 0 depends on the amplitude of the deflection w 0 , and, through the analytical solution of the elastic half space in terms of the Airy stress function, on Young's modulus of the subcortex E s , on Poisson's ratio of subcortex ν s , and on the wavenumber n (Biot, 1937). Inserting this ansatz (4) into the Föppl-von Kármán plate equation (3) yields the following equation for the cortical pressure P, Fig . 4 illustrates the growth-induced cortical pressure P as a function of the wavenumber n for varying cortical and subcortical stiffnesses. The minimum of each curve corresponds to the critical pressure P crit at which folding occurs. The corresponding value of n is the critical wavenumber. The dotted line illustrates the critical pressure for varying wavenumbers. Growing layer on an elastic foundation, left, and on a growing foundation, right. We model cortical folding using the classical fourth order Föppl-von Kármán plate theory and adopt a sinusoidal ansatz for the deflection w, which generates a sinusoidal transverse force q. This provides analytical estimates for the critical cortical pressure P crit and for the wavelength λ crit parameterized in terms of the cortical thickness tc, the cortical and subcortical Young's moduli Ec and Es, and the cortical and subcortical growth rates Gc and Gs.
To determine the critical wavenumber n, we elaborate the minimization problem PðnÞ-min, or, equivalently, dP λ crit p t c , and to the third root of the ratio of the film-to-substrate stiffness, λ crit p ffiffiffiffiffiffiffiffiffiffiffi ffi E c =E s 3 p .

Growing cortex on growing subcortex
To establish analytical estimates for cortical folding on a growing foundation, we again interpret the subcortex as an infinite half-space and impose a sinusoidal deflection wðxÞ ¼ w 0 cos ðnxÞ on its upper boundary. Yet, now we assume that this deflection is the sum of an elastic subcortical deflection w e 0 and subcortical growth w g 0 , such that w 0 ¼ w e 0 þ w g 0 . Similarly, we can additively decompose the deflection rate _ w 0 into an elastic part _ w e 0 and a growth part _ w g 0 (Lubarda, 2004), The elastic deflection, w e 0 ¼ À2½1 À ν 2 s =½E s nq 0 , and its rate, _ w e 0 ¼ À2½1 À ν 2 s =½E s n _ q 0 , follow from inverting the elastic halfspace relation in Eq. (4). The growth deflection rate _ w g 0 represents the growing subcortex. We assume that subcortical growth is stretch-induced, proportional to the elastic subcortical deflection w e 0 , scaled by the subcortical growth rate G s . Combining all three equations yields an equation for the transverse force q 0 as a Maxwell-type viscoelastic response to the deflection w 0 , We introduce a sinusoidal representation of the transverse force q(x) and adopt a convolution type solution for the amplitude q 0 , We choose an exponential ansatz for the deflection amplitude, w 0 ðtÞ ¼ W 0 expðGtÞ, such that dw 0 =ds ¼ Gw 0 ðtÞ, where G is the characteristic time constant of cortical folding. We insert Eq. (10) into the Föppl-von Kármán plate equation (3) to eventually obtain the equation for the cortical pressure, Similar to the elastic foundation, we evaluate the minimization problem PðnÞ-min to determine the critical wavenumber n, The estimates for the critical pressure P crit and the wavelength λ crit then follow immediately as functions of the cortical thickness t c , the cortical and subcortical Young's moduli E c and E s , the cortical and subcortical Poisson's ratios ν c and ν s , the subcortical growth rate G s , and the characteristic time constant of cortical folding G, Again, we can assume that Poisson's ratios of the cortex and subcortex are of the same order, ν c $ ν s . From the remaining parameters, we conclude that the critical wavelength, the distance between two neighboring gyri, is directly proportional to the cortical thickness, λ crit p t c , to the third root of the ratio of the cortical and subcortical stiffnesses, λ crit p ffiffiffiffiffiffiffiffiffiffiffi ffi E c =E s 3 p , and to the subcortical growth rate, λ crit p . At the onset of folding, we assume that the growth-induced pressure in the cortical layer is equivalent to the elastic modulus of the cortex, E c =½1 Àν 2 c , multiplied by the amount of cortical growth. In the simplest case, we can represent cortical growth as the product of the cortical growth rate G c and the critical folding time t crit , By combining the critical pressure P crit at the onset of folding (13.1) with the cortical pressure generated by growth (14), we obtain a critical condition at the onset of folding, expressed in terms of the critical folding time t crit and the characteristic time constant for cortical folding G.
To further evaluate this condition, we choose G ¼ 1=t crit , to obtain a quintic equation for the critical folding time, For given cortical and subcortical Young's moduli E c and E s , cortical and subcortical Poisson's ratios ν c and ν s , and cortical and subcortical growth rates G c and G s , we solve this equation for the folding time t crit and then determine the critical pressure P crit from Eq. (14) and the critical wavelength λ crit from Eq. (13.2).
Fig . 5 illustrates the critical folding pressure P crit , left, and the critical folding time t crit , right, for varying growth ratios G c =G s and stiffness ratios E c =E s . For simplicity, we have assumed that cortex and subcortex have identical Poisson's ratios, ν c ¼ ν s . The graphs agree with our intuition: the larger the cortical stiffness, the larger the required folding pressure P crit and the smaller the folding time t crit . Fig. 6 illustrates the critical wavelength λ crit for varying stiffness ratios E c =E s and varying cortical thicknesses t c , left, and for varying growth ratios G c =G s and varying cortical thicknesses t c , right. The graphs visualize our analytical estimates: the wavelength increases linearly with increasing cortical thickness t c , increases with increasing stiffness ratio E c =E s , and decreases with increasing growth ratio G c =G s .
Remark 2 (Special case of incompressibility). If we assume that the cortex and the subcortex are incompressible, ν c ¼ 0:5 and ν s ¼ 0:5, the analytical estimates for the critical pressure and the wavelength reduce to the following expressions, This simplification agrees with findings in the literature : the wavelength of cortical folding is proportional to the cortical thickness, λ crit p t c , to the third root of the ratio of the cortical-to-subcortical stiffness, . Critical wavelength λ crit for varying stiffness ratios Ec=Es and varying cortical thicknesses tc, left, and for varying growth ratios Gc=Gs and varying cortical thicknesses tc, right. The wavelength λ crit increases with increasing cortical stiffness Ec, with increasing subcortical growth rate Gs, and with increasing cortical thickness tc. Remark 3 (Special cases of solid-and fluid-like subcortex). As indicated in Fig. 3, we have modeled the growing subcortex as Maxwell type viscoelastic solid with _ q 0 þ G s q 0 ¼ ÀE s =½2½1 Àν 2 s n _ w 0 in Eq. (9). This implies that the subcortex will behave solid-like and fluid-like in the two extreme cases, For small subcortical growth rates G s 5G, we immediately recover the solution for the elastic subcortex from Section 2.1. In Fig. 5, left and right, and Fig. 6, right, small subcortical growth rates correspond to large growth ratios G c =G s . We recover the special case of a solid-like subcortex as the asymptotic behavior for G c =G s -1. Increasing the subcortical growth rate increases the wavelength λ crit . For large enough subcortical growth rates G s bG, we can suppress folding entirely. We recover the special case of a fluid-like subcortex as the asymptotic behavior for G c =G s -0 with t crit -1 and λ crit -1. These extreme cases agree with the literature of a wrinkling layer on a viscous substrate (Huang, 2005): the solid-like subcortex corresponds to the extreme case of a glassy substrate; the fluid-like subcortex corresponds to the extreme case of a rubbery substrate.

Continuum model
To explore the folding pattern beyond the onset of folding, we model growth using the nonlinear field theories of mechanics supplemented by the theory of finite growth. This results in a set of five equations, which define the kinematics, the constitutive behavior, the mechanical equilibrium, the growth kinematics, and the growth kinetics. Fig. 7 illustrates our multiscale continuum model for cortical and subcortical growth in the developing mammalian brain.
To characterize the kinematics of finite deformation, we introduce the deformation map φ, which maps points X from the undeformed configuration to their new positions x ¼ φðX; tÞ in the deformed configuration. We then introduce the deformation gradient, which we decompose multiplicatively into an elastic part F e and a growth part F g (Garikipati, 2009), A similar multiplicative decomposition holds for the Jacobian J, which we decompose into an elastic part J e and a growth part J g . To define the growth kinematics, for simplicity, we assume that growth is purely isotropic, parameterized in terms of a single scalar-valued growth multiplier ϑ, This implies that the grown volume J g is identical to the growth multiplier cubed ϑ 3 . In the initial ungrown state, the growth multiplier is one, ϑ ¼ 1, such that ϑ41 and ϑo1 characterize volume growth and shrinkage. The elastic tensor F e and its Jacobian J e then take the following explicit forms, For simplicity, we assume that both cortex and subcortex are isotropic and elastic. We characterize their constitutive behavior through the following Neo-Hookean free energy parameterized exclusively in terms of the elastic tensor F e and its Fig. 7. Multiscale continuum model for cortical and subcortical growth. The cortex, the gray matter, grows morphogenetically at a constant rate Gc. Cortical growth induces subcortical deformation, which triggers subcortical growth. The subcortex, the white matter, grows at a stretch-dependent rate as Gs〈J e À J 0 〉, where Gs mimics the axon elongation rate and 〈J e À J 0 〉 activates growth only, if the elastic volume stretch J e exceeds its baseline value J 0 .
Jacobian J e , for example as where λ and μ are the Lamé constants. This implies that only the elastic part of the deformation induces stress. For the following considerations, it proves convenient to reparameterize the free energy (20) in terms of the total deformation gradient F and the growth factor ϑ, Following standard arguments of thermodynamics, the Piola stress P follows as energetically conjugate to the deformation gradient F, The Piola stress enters the standard balance of linear momentum, the equation of mechanical equilibrium. In the absence of volume forces, the balance of linear momentum reduces to the vanishing divergence of the Piola stress, It remains to define the kinetics of growth, the equations that characterize the evolution of the growth multiplier in time (Menzel and Kuhl, 2012). Since the cortex consists primarily of cell nuclei whereas the subcortex consists primarily of axons, we assume different growth kinetics for the cortex and subcortex.

Cortical growth
For the cortex, we assume that growth is purely morphogenetic (BenAmar and Goriely, 2005), independent of mechanical stress or strain (Ambrosi and Mollica, 2002), characterized exclusively by the growth rate G c , The cortical growth rate G c may vary in time and space, depending on the current stage of development and on the regional location. For simplicity, here we let the cortex grow linearly in time and homogeneously in space, G c ¼ const.

Subcortical growth
For the subcortex, we assume that growth is stretch-induced (Kuhl, 2014). Cortical growth induces extreme deformations in the subcortex. The subcortex is primarily populated by axons, which lengthen gradually when subject to chronic stretch (Bray, 1984). We make the following ansatz, where G s is the subcortical growth rate. With 〈J e À J 0 〉 ¼ J e À J 0 for J e 4 J 0 and 〈J e ÀJ 0 〉 ¼ 0 otherwise, the term in the Macaulay brackets activates growth only if the elastic volume stretch J e exceeds its baseline value J 0 , i.e., when axons are stretched beyond their physiological limit (Dennerll et al., 1989).
Remark 4 (Stresses). To provide a more intuitive illustration of the Piola stress, we re-evaluate the thermodynamic stress definition P ¼ ∂ψ =∂F of Eq. (22), yet now, formulated in terms of the elastic deformation F e , using the Neo-Hookean free energy (20), This indicates that the Piola stress in Eq. (22) is nothing but the growth-weighted classic elastic Piola stress, For the special case of isotropic growth, growth-weighting simplifies to dividing the elastic Piola stress by the growth factor ϑ.

Computational model
To solve the nonlinear equations for brain development, we adopt a finite element discretization in space and a finite difference discretization in time. We introduce the growth multiplier ϑ as an internal variable, and solve its evolution equation for cortical growth (24) and subcortical growth (25) locally at the integration point level. We approximate the growth rate _ ϑ through a finite difference ansatz, where ϑ is the current growth multiplier, ϑ n is the growth multiplier of the previous time step, and Δt ¼ t À t n denotes the current time increment.

Cortical growth
For the cortex, we determine the new growth multiplier ϑ explicitly through a linear update of the previous growth multiplier ϑ n using Eqs. (24) and (28) To solve the global set of equations, we determine the tangent moduli of the cortex through the total derivative of the Piola stress P from Eq. (22) with respect to the deformation gradient F and fix the current growth multiplier ϑ, Since the growth multiplier is independent of the current deformation, the tangent moduli simply consist of the growthweighted classical elastic tangent moduli (Papastavrou et al., 2013), Here, we have used the following abbreviations, f ○g ijkl ¼ fg ik f○g jl and f ○g ijkl ¼ fg il f○g jk , for the non-standard fourth order products.

Subcortical growth
For the subcortex, we apply an implicit time integration scheme and reformulate the evolution equation (25) with the help of the finite difference ansatz (28). This introduces the discrete residual R in terms of the unknown growth multiplier ϑ, To solve this nonlinear equation, we adopt a local Newton iteration. For each iteration step, we calculate the linearization of the residual R with respect to the current growth multiplier ϑ, The Heaviside step function H is one during growth, J=ϑ 3 À J 0 40, and zero otherwise. Within each iteration step, we update the unknown growth multiplier, ϑ'ϑÀR=K; until we achieve local convergence, i.e., until the absolute value of the growth update jR=Kj is smaller than a user-defined convergence threshold. To solve the global set of equations, we determine the tangent moduli of the subcortex through the total derivative of the Piola stress P from Eq. (22) with respect to the deformation gradient F, The first term, the Hessian of the free energy function for constant growth, ϑ¼const., defines the growth-weighted classical elastic tangent moduli similar to the case of cortical growth (31), The second term is specific to the constitutive equation (20) and characterizes the sensitivity of the Piola stress P with respect to the growth multiplier ϑ, The third term is specific to both the kinetic equation for growth (25) and its algorithmic solution (33), where ∂ _ ϑ=∂ϑ ¼ K and ∂ _ ϑ=∂F ¼ G s JF À t =ϑ 3 . The Hessian of the free energy function for constant deformation, F ¼ const., defines the correction of the constitutive moduli due to growth, The fourth-order tangent moduli for the cortex (30) and for the subcortex (35) enter the iteration matrix for the global Newton iteration. Upon convergence of the global Newton iteration, we store the current growth multiplier ϑ locally at the integration point level.
Remark 5 (Tangent moduli). Again, we can reformulate the definition of the elastic tangent operator A e ¼ ∂P=∂F in terms of the elastic deformation F e to obtain a more intuitive interpretation of its terms, This indicates that the elastic tangents in (31) and (36) are nothing but the growth-weighted standard tangent operator, For the special case of isotropic growth, this growth-weighting simplifies to dividing the tangent by the growth factor squared ϑ 2 .

Results
To illustrate the features of our computational model, we expand the analytical study in Section 2 beyond the onset of folding, which is difficult to assess analytically and has only been addressed recently by a few groups Cao and Hutchinson, 2012). During brain development, constrained growth of the cortical layer induces compressive stresses P. Once these stresses reach a critical value P crit , the brain surface buckles into a wavy pattern to partially release the growth-induced stress. In analogous to the analytical estimates in Section 2 and Fig. 6, we explore the role of the three main contributors to the folding pattern: the cortical thickness t c , the stiffness ratio between cortex and subcortex E c =E s , and the growth ratio between cortex and subcortex G c =G s .

Sensitivity of surface morphology with respect to cortical thickness
To explore the effect of the initial cortical thickness t c on the folding pattern, we explore constrained growth in a regular rectangular slice of 2 Â 1 Â 0.05 cm 3 of a transverse brain section. We discretize the slice with 80 Â 40 Â 1 ¼3200 tri-linear Q1 elements and 19,926 degrees of freedom and assume a plane strain state. To constrain growth, we fix the left, bottom, and right boundary nodes orthogonal to the boundary, but allow them to slide freely along the edge. We model the cortex as Neo Hookean elastic with Lamé constants λ c ¼ 34:2 kPa and μ c ¼ 3:3 kPa (Soza et al., 2005) and assume that the subcortex is three times softer with λ s ¼ 11:4 kPa and μ s ¼ 1:1 kPa (Budday et al., 2014). We fix the cortical growth rate to G c ¼ 2:0, the subcortical growth rate to G s ¼ 0:003, and the physiological limit for axonal growth to J 0 ¼ 1:0 (Chada et al., 1997). Since a perfectly regular rectangular domain would not fold in the computational simulation, we trigger an initial imperfection by selectively increasing the subcortical growth rate G s by 10% in a 0.05 cm thin vertical band in the center of the rectangle . We gradually increase the initial cortical thickness from t crit ¼ 0:125 mm to t crit ¼ 1:000 mm in eight equal steps of Δt crit ¼ 0:125 mm. Fig. 8 illustrates the sensitivity of the wavelength λ crit with respect to the initial cortical thickness t c . The eight dots indicate the computationally predicted wavelengths for the eight different cortical thicknesses. The eight figures inside the graph additionally illustrate the corresponding folding patterns. Simulations with coarser and finer meshes and with smaller and larger perturbations predicted similar folding patterns and similar wavelengths. This indicates that the computationally predicted surface morphology is relatively insensitive to the underlying discretization and to the imposed mode of perturbation. The dashed line shows the analytical wavelength-thickness relation for a growing cortical layer on a growing subcortical foundation according to Section 2. According to Eq. (8), the analytical solution is based on the additive decomposition of the rate of deformation into elastic and growth parts. It follows from evaluating equation (13), the equation to estimate the critical wavelength, λ crit ¼ 2πt c , as a function of the characteristic time scale of folding G, which we obtain from iteratively solving the quintic equation (16) as G 5 À64G 3 c ½G s þ G 2 6 0. The solid line shows the averaged computational wavelength-thickness relation for a morphogenetically growing cortex on a stretch-driven growing subcortex according to Sections 3 and 4. According to Eq. (17), the computational prediction is based on the multiplicative decomposition of the deformation gradient into elastic and growth parts. Fig. 8 suggests that the analytical estimate with the additive decomposition and the computational prediction with the multiplicative decomposition agree well in the small deformation limit (Li et al., 2011;Lubarda, 2004). Their direct comparison confirms that the wavelength increases linearly with the initial cortical thickness, λ crit p t c . The slope of m ¼10.20 of the computational model is slightly higher than the slope of m ¼8.92 of the analytical model, which indicates that the computational prediction is slightly stiffer than the analytical estimate. This discrepancy is in agreement with the well-known overly stiff response of tri-linear finite elements, in particular in bending-dominated problems. In addition, the chosen homogeneous Neumann boundary conditions at the lateral sides enforce symmetric folding patterns for which the computationally predicted wavenumbers are always multiples of one half. Aside from these limitations, the analytically and computationally predicted surface morphologies are in excellent quantitative and qualitative agreement. Yet, additional discrepancies might arise between the wavelengthto-thickness relation on initially flat geometries as studied here and initially curved geometries (Li et al., 2011) of real mammalian brains.

Sensitivity of surface morphology with respect to stiffness ratio
To explore the effect of the stiffness ratio E c =E s on the folding pattern, we simulate growth of an elliptic brain slice to mimic an idealized transverse brain section during early development. While the cortical thickness t c is a parameter that is easy to measure experimentally, the cortical and subcortical stiffnesses E c and E s are relatively difficult to determine. On one hand, in vivo experiments on living brain tissue seem virtually impossible. On the other hand, it remains questionable to which extent ex vivo experiments can provide useful estimates for the material properties of the living brain in vivo. Some effort has been made to identify the elastic material parameters of the human brain , yet, the reported values deviate considerably: Young's modulus was found to vary four orders of magnitude, from 0.5 kPa to 500 kPa, and even Poisson's ratio was reported to range from 0.2 to 0.5 (Franceschini, 2006). Due to the structural difference between neuronal cell bodies and neuronal axons, it seems reasonable to assume that the stiffness is not even uniform across the brain, and that cortical and subcortical stiffnesses are inherently different.
Although there is no general agreement of the absolute stiffness values of cortex and subcortex, we can still perform a sensitivity study of the cortical-to-subcortical stiffness ratio E c =E s . To this end, we simulate growth of an elliptic brain slice of an area of 4 cm 2 , an ellipticity ratio of 1.15, and a thickness of 0.005 cm. We discretize the ellipse with 3328 tri-linear Q1 elements and 20,358 degrees of freedom and assume a plane strain state. This discretization introduces 128 nodes on the outer boundary, which implies that it can capture a folding pattern with 16-folds at a resolution of eight notes per wavelength. We model the cortex as Neo Hookean elastic with Lamé constants λ c ¼ 34:2 kPa and μ c ¼ 3:3 kPa and systematically double the cortical-to-subcortical stiffness ratio E c =E s from 2 2 to 2 5 in four subsequent steps by adjusting the subcortical Lamé constants λ s and μ s . We fix the cortical growth rate to G c ¼ 2:0, the subcortical growth rate to G s ¼ 0:003, and the physiological limit for axonal growth to J 0 ¼ 1:0 (Chada et al., 1997). In addition, we gradually vary the initial cortical thickness from t crit ¼ 0:25 mm to t crit ¼ 1:00 mm in four equal steps of Δt crit ¼ 0:25 mm. In contrast to the perfectly regular rectangular domain, the elliptic domain possesses an inherent imperfection because of its varying curvature, and we do not need to impose additional artificial imperfections to trigger folding. Fig. 9 illustrates the sensitivity of the wavelength λ crit with respect to the initial cortical thickness t c and the stiffness ratio E c =E s . Snapshots of each column have the same stiffness ratio E c =E s and are displayed at the same stage of cortical growth ϑ, the stage of first self-contact within the corresponding column. Snapshots of each row have the same cortical thickness t c Fig. 8. Sensitivity of surface morphology with respect to initial cortical thickness for constrained growth in a rectangular domain. The dots illustrate the computationally predicted wavelengths λ crit for varying cortical thicknesses tc. The solid line shows the averaged computational wavelength-thickness relation for a morphogenetically growing cortex on a stretch-driven growing subcortex. The dashed line shows the analytical wavelength-thickness relation for a growing cortical layer on a growing subcortical foundation. and are displayed for subsequent stages with increasing cortical growth ϑ as the stiffness ratio increases. The computational simulation agrees well with the analytical estimates in Section 2: the wavelength λ crit increases with increasing cortical thickness t c , from top to bottom, and with increasing stiffness ratio, E c =E s , from left to right. In all cases, folding started first in the region of lowest curvature, on the shorter symmetry axis, and gradually propagated outwards to the regions of highest curvature. While cortical growth is identical in all 16 cases, and homogeneously distributed across the entire cortex, subcortical growth varies significantly across the 16 simulations and displays pronounced regional heterogeneities. In general, larger cortical wavelengths induce larger subcortical stretch resulting in larger subcortical growth. As the wavelength increases, the individual folds become deeper. As a consequence, subcortical growth is largest in the gyri and smallest in the sulci. Fig. 10 summarizes the computationally predicted average wavelength λ crit , i.e., the elliptical circumference divided by the number of folds n, for the four varying cortical thicknesses and the four varying stiffness ratios. In agreement with the analytical estimates in Section 2, the average wavelength increases linearly with increasing cortical thickness t c , from left to right. Also in agreement with the analytical estimates, the average wavelength increases with increasing stiffness ratio E c =E s , from lower blue dots to upper red dots.

Sensitivity of surface morphology with respect to growth ratio
To explore the effect of the growth ratio G c =G s on the folding pattern, we simulate the same idealized elliptic transverse brain section as in Section 5.2 with an elliptic area of 4 cm 2 , an ellipticity ratio of 1.15, and a thickness of 0.005 cm. Again, we discretize the ellipse with 3328 tri-linear Q1 elements and 20,358 degrees of freedom and assume a plane strain state.
We model the cortex as Neo Hookean elastic with Lamé constants λ c ¼ 34:2 kPa and μ c ¼ 3:3 kPa (Soza et al., 2005) and assume that the subcortex is three times softer with λ s ¼ 11:4 kPa and μ s ¼ 1:1 kPa (Budday et al., 2014). We fix the subcortical growth rate to G s ¼ 0:003 and fix the physiological limit for axonal growth to J 0 ¼ 1:0. We systematically increase the cortical-to-subcortical growth ratio G c =G s from 10 À 1 via 10 0 and 10 1 to 10 2 . Similar to the previous example, we also vary the initial cortical thickness from t crit ¼ 0:25 mm to t crit ¼ 1:00 mm in four equal steps of Δt crit ¼ 0:25 mm. Fig. 11 illustrates the sensitivity of the wavelength λ crit with respect to the initial cortical thickness t c and the growth ratio G c =G s . All snapshots correspond to the stage of cortical growth, at which the final folding pattern of all 16 ellipses had fully developed. Again, folding started first in the region of lowest curvature, on the shorter symmetry axis, and gradually propagated outwards to the regions of highest curvature. Cortical growth is identical in all 16 cases and homogeneously distributed across each slice, whereas subcortical growth varies significantly across the 16 simulations and displays Fig. 9. Sensitivity of surface morphology with respect to initial cortical thickness and stiffness ratio for elliptic geometry. The wavelength λ crit increases with increasing cortical thickness t c , from top to bottom, and with increasing stiffness ratio, Ec=Es, from left to right. Larger wavelengths induce larger subcortical stretch resulting in larger subcortical growth.
pronounced regional heterogeneities. The left column indicates that slow cortical growth rates G c allow for balanced cortical and subcortical growth, which keeps the wavelength uniform and generates simple sinusoidal folding patterns. In contrast, large cortical growth rates G c scale down subcortical growth. This generates higher compression in the cortex, which initiates the formation of secondary folds. According to the analytical estimates from Section 2, shorter critical wavelengths require a larger critical pressure before buckling is induced. For instance, in the bottom, left corner of Fig. 11, growth has not Fig. 11. Sensitivity of surface morphology with respect to initial cortical thickness and growth ratio for elliptic geometry. The wavelength of primary folding λ crit increases with increasing cortical thickness t c , from top to bottom. The overall wavelength increases with increasing growth ratios Gc=Gs, from left to right. Fig. 10. Sensitivity of surface morphology with respect to initial cortical thickness and stiffness ratio for elliptic geometry. The dots illustrate the computationally predicted average wavelengths λ crit for varying cortical thicknesses tc and varying stiffness ratios Ec=Es. The average wavelength increases with increasing cortical thickness t c , from left to right, and with increasing stiffness ratio, Ec=Es, from lower blue dots to upper red dots. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.) yet generated sufficient compression to induce folding. If the cortex keeps on growing, this ellipse will buckle at a shorter wavelength than any of the other shown ellipses. The concurrence of slow growth rates and thick cortices prevents reaching a load high enough to initiate buckling.

Discussion
Despite its tremendous significance, little is known about the origin of cortical folding in the developing mammalian brain (Bayly et al., 2014). Two popular but competing hypotheses suggest that cortical folding originates either in the subcortex, driven by axonal tension (Van Essen, 1997), or in the cortex, driven by differential growth (Richman et al., 1975).
Here we have combined both hypotheses into a bilayered material model for cortical folding, in which we represent the cortex as a morphogenetically growing outer layer (Holland et al., 2013), and the subcortex as a strain-driven growing inner core (Budday et al., 2014).
To gain first insight into these competing mechanisms , we have established analytical estimates for the critical time, pressure, and wavelength at the onset of folding by modeling cortical folding by means of the instability problem of a confined, layered medium under growth-induced compression (Biot, 1957). We have shown that the critical wavelength λ crit , the distance between two neighboring gyri, is directly proportional to the cortical thickness, λ crit p t c , proportional to the third root of the cortical-to-subcortical stiffness ratio, λ crit p , and proportional to the third root of the subcortical growth rate, λ crit p ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðG þG s Þ=G 3 p . To explore the evolution of surface morphologies beyond the onset of folding, we have proposed a continuum model for finite growth, which we have solved computationally within a nonlinear finite element setting (Ambrosi et al., 2011). In regular rectangular geometries, we had to apply a small perturbation to trigger the formation of instabilities (Papastavrou et al., 2013). In agreement with the literature , we observed that the computationally predicted surface morphology was relatively insensitive to the imposed mode of perturbation. In elliptic geometries, the heterogeneity in curvature was sufficient to initiate folding (Eskandari et al., 2013). The instability originated at the center of the long axis and then spread symmetrically outward.
We have systematically varied cortical thickness, stiffness, and growth and predicted folding patterns that were in excellent agreement with our analytical estimates. As expected, our computational model predicted a much wider variety of surface morphologies than the analytical solutions. In some cases, it even predicted the formation of secondary folds . In general, folding patterns deviated from the symmetric sinusoidal ansatz towards morphologies with larger gyri and smaller sulci. This asymmetry reflects the impact of chronic axon elongation on gyral regions with positive stretch (Bray, 1984), which induces subcortical growth. As a natural consequence, the subcortical growth multiplier ϑ displays significant regional variations with maxima of ϑ42:0 in the gyral centers and minima of ϑ ¼ 1:0 indicating no growth at the sulcal base. This confirms that analytical modeling can provide valuable first insight into regular folding patterns (Biot, 1957), but computational modeling is mandatory to explore irregular brain surface morphologies (Budday et al., 2014).

Sensitivity of surface morphology with respect to cortical thickness
Of all parameters, our model seems to be most sensitive to variations in cortical thickness. Our simulations suggest that the intersulcal distance increases linearly with increasing cortical thickness (Biot, 1937). A considerably thickened cortex can even suppress the formation of folds entirely (Raybaud and Widjaja, 2011). This tendency is consistent with clinical pictures of diseased human brains: Lissencephaly, a malformation with a markedly thickened cortex, is characterized by a smooth brain surface (Landrieu et al., 1998); Polymicrogyria, a malformation with a regionally thinned cortex, is characterized by a highly convoluted brain surfaces with many small and superficial folds (Tortori-Donati et al., 2005). Thin cortices and decreased gyrification are associated with epilepsy, attention deficit hyperactivity disorder, dementia, mental retardation, and dyslexia; thick cortices and increased gyrification are associated with Williams syndrome, autism, and schizophrenia (Zilles et al., 2013). We conclude that the cortical thickness directly influences the gyral wavelength and is a key parameter to control surface morphology and primary folding.

Sensitivity of surface morphology with respect to stiffness ratio
The stiffness ratio between cortex and subcortex has a similar effect as the cortical thickness, however, only when scaled by its third root. Controlling surface morphology through the stiffness ratio has been discussed intensely in the materials sciences community (Cai et al., 2011), where thin stiff films on compliant substrates with stiffness ratios of up to four orders of magnitude play a major role . As the cortex with its neuronal cell bodies and synapses is much denser packed than the subcortex with its myelinated axons, it is intuitive that it may have a larger mechanical stiffness. Yet, experiments have shown that the cortical stiffness is less than an order of magnitude larger than the subcortical stiffness (vanDommelen et al.,). Some studies only found a stiffness difference of 50% (Christ et al., 2010). This has, in fact, been the major criticism of the first mechanical model for cortical folding based on the hypothesis of differential growth (Richman et al., 1975). We conclude that the stiffness ratio may influence surface morphology, but because of its small variation, it cannot be the single main driving force to explain cortical folding and morphological abnormalities.
6.3. Sensitivity of surface morphology with respect to growth ratio Of all parameters, the growth ratio between cortex and subcortex seems to be the least well understood. Yet, it is probably the most important parameter to control the formation of secondary folds (Zang et al., 2012). In the continuum model, cortical and subcortical growth are introduced constitutively through the kinetics of growth in Eqs. (24) and (25). Eventually, we hope to tie these equations to cellular mechanisms such as axon elongation (Bray, 1984). In the analytical model, we have made a critical assumption to evaluate the relation between the critical wavelength, the time constant of growth, and the subcortical growth rate, namely that G ¼ 1=t crit . At this point, this is a plain assumption, yet it provides some insight into the two extreme cases of abnormally slowly and abnormally fast growing cortices. Abnormally slowly growing cortices create an almost fluid-like behavior of the subcortex: axons are capable of responding almost instantaneously to growth-induced subcortical deformation and, in extreme cases, folding is suppressed entirely. Abnormally fast growing cortices create an elastic solid-like behavior of the subcortex: axons are incapable of responding to stretch, the pressure in the cortical layer raises quickly, and provokes the formation of secondary folds (Li et al., 2012). In the human brain, for example, primary folding begins at 22 weeks gestation and secondary folding takes place between weeks 25 and 30 (Raybaud and Widjaja, 2011). In the materials science community, these two types of folds are associated with kinematically induced instabilities, controlled by thickness and stiffness, and dynamically induced instabilities, controlled by growth rates (Huang, 2005). The interplay of kinematic and dynamic instabilities generates a wide variety of surface morphologies, in which secondary folding serves as a mechanism to release large compressive stresses in the outer layer . We conclude that the growth ratio is a key parameter to control irregular surface morphologies and secondary folding.

Mechanical modeling explains surface morphologies of mammalian brains
From the explanted mammalian brains in Fig. 2, we conclude that brain size increases with increasing animal size (Hofman, 1989) from 140 g in sheep via 180 g in pigs to 450 g in cows (Nieuwenhuys et al., 1997). Yet, the average cortical thickness in the frontal coronal sections in Fig. 2 varies only marginally from 0.23 cm in the sheep via 0.22 cm in the pig to 0.22 cm in the cow brain. This is in line with the common understanding that the average cortical thickness varies marginally across mammalian species and is independent of brain size (Roth and Dicke, 2005). Similarly, the average gyral wavelength in Fig. 2 varies marginally from 0.61 cm in the sheep via 0.53 cm in the pig to 0.90 cm in the cow brain. According to our model, increasing the brain size at a constant cortical thickness does not affect the absolute gyral wavelength; yet, it increases the relative gyral wavelength when scaled by brain size. This is in agreement with a recent review, which reported the gyrification index, the ratio between the total brain surface area and the exposed surface area, to be 1.94 in sheep and 2.18 in pig (Zilles et al., 2013). Our model can thus explain why the surface-to-volume ratio of mammalian brains in Fig. 1 increases disproportionally faster than predicted by isometric scaling and why the degree of gyrification tends to increase with brain size.

Limitations
Both our analytical and our computational model provide valuable insight into the development of the mammalian brain. Yet, they have a few limitations, which could be addressed to make the models more realistic.
First, in our current model, we have neglected the geometric constraint by the skull. Our results indicate that the skull is not necessary to trigger gyrification; yet, it might be an important regulator of cortical folding (Nie et al., 2010). In our model, folding is constrained exclusively by the subcortical layer underneath the growing cortex. Adding a stiff skull above the growing cortex would certainly influence the final folding pattern, for example, by flattening out the gyral ridges. Yet, the impact if the skull on the initial gyral wavelength at the onset of folding might be rather minor .
Second, here, we have assumed the constitutive behavior as quasi-incompressible with a Poisson's ratio of ν ¼ 0:458 in the linear regime (Soza et al., 2005). Other studies suggest that brain tissue is nearly incompressible with ν ¼ 0:496 (Franceschini, 2006) or entirely incompressible (Rashid et al., 2012(Rashid et al., , 2014). Our analytical model in Section 2 is generally valid for both compressible and incompressible materials. Our continuum model in Section 3 would require a different strain energy function in Eq. (20) to capture incompressibility exactly. Imposing incompressibility in the computational model in Section 4 would require additional modifications both on the constitutive level and on the element level, since the overall response would alternate between compressible during growth and incompressible during purely elastic phases (Rausch and Kuhl, 2014;Schmid et al., 2012).
Third, within this study, we have simplified the elastic response of the brain as Neo Hookean isotropic (Soza et al., 2005). Although the microstructure of both cortex and subcortex is clearly anisotropic, their anisotropic material behavior remains poorly characterized. The microstructure of the cortex is relatively regular and closely correlated with the radial direction r 0 , which defines the orientation of the cortical columns (Rakic, 1988). The microstructure of the subcortex is rather irregular and closely correlated with the axon orientation a 0 , which needs to be identified through diffusion tensor imaging (Mori and Zhang, 2006) or other novel imaging techniques (Chung and Deisseroth, 2013). Detailed measurements of the cortical stiffness, along and perpendicular to the cortical columns, supplemented by measurements of the subcortical stiffness, along and perpendicular to the axon orientation, would be tremendously valuable. To improve the elastic module of our model, we are currently performing a series of nano-indentation tests (Zhang et al., 2010) to characterize cortical and subcortical stiffnesses, the degree of cortical and subcortical anisotropy, and the stiffness variation across different species.
Fourth, we have not only simplified the elastic behavior but also the growth response as isotropy. For cortical growth, we are currently working on replacing the isotropic growth tensor, F g ¼ ϑI, in Eq. (18) by an anisotropic growth tensor, F g ¼ ϑ ? I þ½ϑ J À ϑ ? r 0 r 0 , where r 0 characterizes the radial direction (Rausch and Kuhl, 2014). In this setup, the radial growth multiplier ϑ J characterizes cortical thickening along radial direction (Göktepe et al., 2010), and the surface growth multiplier ϑ ? characterizes area growth perpendicular to it (Buganza Tepole et al., 2011). This allows us to replace the single phenomenological evolution of cortical growth, _ ϑ ¼ G c , in Eq. (24) by two independent equations for cortical thickening and surface growth. We can then mechanistically link surface growth, _ ϑ ? ¼ G ? c , to the symmetric cell division of progenitor cells into two new progenitor cells, and thickness growth, _ ϑ J ¼ G J c to the asymmetric cell division into a progenitor cell and a neuron (Roth and Dicke, 2005). For subcortical growth, we could replace the growth tensor F g ¼ ϑI, in Eq. (18) by an anisotropic growth tensor, F g ¼ I þ½ϑ J À1a 0 a 0 , where a 0 characterizes the axon orientation (Zöllner et al., 2012). The growth multiplier ϑ J characterizes the chronic axon elongation in response to overstretch (Dennerll et al., 1989). We could then replace elastic volume change 〈J e À J 0 〉 as the driving force for subcortical growth, _ ϑ J ¼ G s 〈J e À J 0 〉, in Eq. (25) by the elastic axonal stretch 〈λ e À λ 0 〉 with λ e ¼ ½a 0 Á F t Á F Á a 0 1=2 , to correlate the model parameters to experimentally measured axon elongation rates (Bray, 1984).
Finally, growth is neither homogeneous in space nor constant in time. To functionally correlate cellular and molecular events to cortical and subcortical growth (Knutsen et al., 2013), we could turn different growth rates on and off to better represent the sequence of events during gyrogenesis, including neuronal proliferation, differentiation, apoptosis, dendrogenesis, synapsogenesis, glial proliferation, lamination, and cellular rearrangement (Raybaud and Widjaja, 2011).

Concluding remarks
Mechanical modeling of brain development can explain variations in surface morphology of the mammalian brain. Variations in cortical and subcortical thickness, stiffness, or growth can generate variations in pattern formation. A thinner, softer, or slower growing layer of gray matter generally enhances cortical folding and reduces the gyral wavelength. A thicker, stiffer, or faster growing layer of gray matter reduces cortical folding and increases the gyral wavelength. Larger mammals tend to have larger brains, but similar cortical thicknesses. Our model predicts that the absolute gyral wavelength in mammals is almost constant across different species, while the relative gyral wavelength increases with brain size. This explains why the surface-to-volume ratio of mammalian brains increases disproportionally faster than predicted by isometric scaling. Our model can also explain the pathological malformations of polymicrogyria, associated with a thin and overly convoluted cortex, and lissencephaly, associated with a thick and poorly convoluted cortex. Understanding the mechanisms of cortical folding during brain development may have direct implications on the diagnostics of neurological disorders, including severe retardation, epilepsy, schizophrenia, and autism.