Elastic instabilities in a layered cerebral cortex: A revised axonal tension model for cortex folding

We model the elasticity of the cerebral cortex as a layered material with bending energy along the layers and elastic energy between them in both planar and polar geometries. The cortex is also subjected to axons pulling from the underlying white matter. Above a critical threshold force, a"flat"cortex configuration becomes unstable and periodic unduluations emerge, i.e. a buckling instability occurs. These undulations may indeed initiate folds in the cortex. We identify analytically the critical force and the critical wavelength of the undulations. Both quantities are physiologically relevant values. Our model is a revised version of the axonal tension model for cortex folding, with our version taking into account the layered structure of the cortex. Moreover, our model draws a connection with another competing model for cortex folding, namely the differential growth-induced buckling model. For the polar geometry, we study the relationship between brain size and the critical force and wavelength to understand why small mice brains exhibit no folds, while larger human brains do, for example. Finally, an estimate of the bending rigidity constant for the cortex can be made based on the critical wavelength.


I. INTRODUCTION
The cerebral cortex, or grey matter, is the outermost layer of nerve tissue covering the cerebrum and plays a key role in high-level cognitive functions, such as decision-making. The nerve cells in the cerebral cortex contain nonmyelinated axons, and the cortex is distinguished from the underlying nerve tissue consisting of nerve cells with myelinated axons, otherwise known as white matter. The geometry of the cortex varies across mammals [1]. See Figure 1. In mice, the cortex is smooth, while in larger mammals, the cortex develops folds. These folds allow for greater surface area of the cortex such that more neurons can participate in, and, therefore, presumably enhance higher-level cognitive functions. To date, there are two competing mechanisms proposed to drive cortex folding. The first proposal claims the folds are driven by axonal tension from the underlying white matter drawing the sides of gyri (outward folds) together [2]. See Figure 2a. This mechanism is appealing because it can be related to the efficient wiring of neurons via the minimization of distances. Moreover, this model does not invoke any elastic instabilities, i.e. buckling. The second, competing proposal suggests that the folds are driven by buckling [3]. See Figure 2b. More specifically, fast growth of the outermost layer of the cortex produces compressive stress that leads to buckling of this layer as modulated by the stiffness of the underlying foundation (comprised of the remaining layers of the cortex and the white matter). Interestingly, this type of buckling model can also be invoked to study many shape changes in nature ranging from plant growth to geological folds [4].
FIG. 2: Schematic of (a) the axonal tension model, which distinguishes between the cortex (denoted by the grey shading) and the underlying white matter, (b) the differential growth model, which distinguishes between the top layer of the cortex and the rest of the brain matter, (c) the model presented here. The red lines denote axons, the black arrows denote the direction of the force. Only three layers of the cortex are drawn for simplicity.
Indirect evidence for each mechanism exists. For instance, in fetal brains where most of the tissue below the cortex is surgically ablated prior to folds developing, folds eventually do develop [5]. This observation suggests that the intracortical buckling drives folding and not axonal tension from the underlying white matter. In addition, a quantitative model of buckling of an elastic plate (the top layer of the cortex) supported by an elastic foundation (the white matter) yields a critical wavelength for buckling that agrees with the typical distance between folds, provided the Young's modulus of the white matter is 10 times less than the grey matter [3]. This assumption, however, has been called into question [6].
For the axonal tension model, as originally formulated, neuronal pathways connecting gyri should be denser than those connecting sulci (inward folds). Some data supports this notion, though the results may be a matter of defining which surrounding regions belong to gyri and which belong to sulci [7]. Moreover, cortical folds generated by linking different areas of the brain via axonal tension means that denser neuronal pathways should exhibit straighter whitematter trajectories. There exists data to support this notion [8]. On the other hand, cuts in feret brain tissue indicate that the tension does not run between gyri but radially outward [9]. Quantitative data for the axonal tension model at the same level of the buckling mechanism is currently lacking.
Indeed, it could very well be that both mechanisms are at play in the folding of the cortex. If so, can we distinguish between the two? To begin to do so, we develop a new model of the elasticity of the cortex that takes into account (1) the elongated, or rod-like, structure of nerve cells sitting in a "background" of softer, glial and progenitor cells, (2) the layered structure of the cortex [10] (See Figure 3), and (3) brain tissue (the cortex included) is a viscoelastic material with both a non-zero storage modulus and non-zero loss modulus [11]. The combination of these three ingredients may make it reasonable, at some time scale, to model the cortex as a layered liquid crystal with the neurons representing the liquid crystal molecules.
With this "cortex as a liquid crystal structure" in what will turn out to be the smectic phase, we can revisit the axonal tension model and investigate the effect of pulling forces on the cortex. We will do this in both a planar geometry and a polar geometry and demonstrate that "vertical pulling" of the axons in the underlying white matter can lead to buckling in a layered structure. Our analysis allows for an updated version of the axonal tension hypothesis that is more consistent with the data. Prior to this work, all buckling models for cortex folding are based on "horizontal compression".
The paper is organized as follows. The next section details our "cortex" as a smectic" approach in a planar geometry to estimate the critical force needed to generate cortical folds with some critical wavelength. Section III presents results for the polar geometry. Section IV summarizes our results and their implications.

II. CORTEX AS A LIQUID CRYSTAL SMECTIC: PLANAR GEOMETRY
The cortex consists of neurons, glial cells, and progenitor cells. The glial cells provide nutrients for the nerve cells they surround. Progenitor cells eventually become nerve cells (since nerve cells do not divide). The shape of each nerve cell is rod-like with a cross-sectional diameter of order a micron and a length ranging from several hundred of microns to approximately a millimeter. The mechanical rigidity along the axon of the nerve cell is provided for by microtubules. Microtubules are semiflexible polymers with a persistence length of approximately 1 mm [12]. Therefore, the nerve cells are rather rigid "molecules". The surrounding glial cells are softer [13]. Moreover, they presumably provide for most of the viscosity observed in indentation experiments [13]. Since the glial cells are softer, we will assume that the elasticity of the cortex is dominated by the rigid, rod-like nerve cells.
How are these rigid, rod-like nerve cells arranged in the cortex? As indicated in Figure 3, they are predominantly oriented perpendicular to the outer surface of the cortex [1]. Moreover, the nerve cells in the cortex arrange themselves into six layers with the morphology differing slightly between layers. For simplicity, we assume that all the layers are equivalent in thickness and in elastic properties. Given this extra spatial structure, we model the elasticity of the cortex as a smectic liquid crystal and then ask the following: What are the consequences of axons from the underlying white matter pulling vertically on the cortex in this planar geometry? The pulling of axons (nerve cells) has been well-established [14] and given the orientation of axon highways in the underlying white matter [15], vertical pulling is in keeping with observations. So, as with the original version of the axonal tension model, here, the white matter enters the model solely via an applied strain and via boundary conditions. We will also include the effect of uniform cortical growth, as opposed to differential cortical growth, to begin to look for potential interplay between axonal pulling and growth driving cortex folding.
To quantify the effect of the applied vertical stain on the cortex smectic, we consider the set of surfaces ω(x, k) = x · n − kl = 0, where x denotes the position in Cartesian coordinates (x, y, z), n is the unit normal to the layer, k ∈ Z, and l is the interlayer spacing [16]. We will assume translational invariance in the y-coordinate such that the model is two-dimensional with the layers described by curves. The deformation of layers (curves), from the initial configuration described by x 0 = xe x + ze z with layer normal n 0 = e z along z-axis, to the current configuration, x, is characterised by the deformation gradient F = ∂x/∂x 0 . In this planar geometry, we assume the following mapping x → (1 + α)x and z → z + U (x, z), where α characterizes the lateral growth of the cortex and U (x, z) is a displacement field due to vertical pulling of the axons. Then, the deformation gradient matrix and its inverse transposed in {e x , e z } basis are Thus, the spatial gradient of isosurfaces in the current configuration can be computed as [16] The thickness of the deformed layers corresponds to l 0 /|∇ω|, where l 0 is the thickness of undeformed layers. The elastic free energy density, accounting for the finite deformation such as compression and bending of layers, consists of two terms, respectively where B is compression modulus and K is the bending rigidity. The ratio K/B defines the characteristic length scale of the order of the layer thickness ≃ 1 mm. The displacement field, U (x, z), can be further decomposed into a uniform dilation of layers along z-direction and an inhomogeneous displacement, u(x, z), so that U (x, z) = γz + u(x, z). Assuming u ≪ 1 and expanding the layer dilation and the unit normal we find The regularity of the wavelength of brain folds allows us to assume the selection of the certain wavelength and thus look for the periodic solution, u(x, z) = φ(z) cos(qx). Replacing this ansatz into Eq. 3 and integrating over the period 2π/q we arrive at the free energy, where we introduced the dimensionless variables This free energy has not previously been studied beyond the harmonic limit [16]. As for parameters, the thickness of the cortex scales does not vary much among the mammals. More specifically, h ≃ 2 − 5 mm [10]. Note that the uniform growth factor α is absorbed in the coefficients λ and dimensionless wavelength χ. Similar to the Helfrich-Hurault instability in nonliving smectic liquid crystals [16][17][18][19] we expect, that above some critical threshold γ cr , undulations of layers are energetically favoured to minimize the compression energy in expense of the bending. More precisely it follows from Eq. 6 that if γ(1 + γ) > λ 2 χ 2 , a non-trivial solution φ = 0 exists [17]. This is necessary but not sufficient condition for the buckling profile to be favoured. Below we derive the stability criterion, which relates the control parameter γ and the wavelength L x = 2π/q with thickness h and elastic moduli K/B, which are intrinsic parameters of the system.
In equilibrium, we require the vanishing of the first variation of the free energy δF = 0, yielding the Euler-Lagrange equation associated with Eq. 6 and the boundary conditions. We assume that at the lower interface grey-white matter (z = 0) the displacement field vanishes, φ|z =0 = 0. The upper interface with surrounding fluid (z = 1) is free, and the condition ∂zφ|z =1 = 0 should be satisfied to guarantee δF = 0.
The first integral of the Euler-Lagrange equation associated with Eq. 6 is Rewriting the potential on the RHS as V(φ) = 3 16 (1 + γ) 2 χ 4 (φ 2 − φ 1 )(φ 2 − φ 2 ) and intergrating, we find the general solution in terms of the Jacobi elliptic function, or which satisfies the boundary conditionφ|z =0 = 0. The period of Eq. 9 is 4K/ν, where K is the complete elliptic integral of the first kind. The condition at the upper interface, corresponding to the maximum of the displacement field, determines the threshold criterion which is with µ = χ 2 (γ + γ 2 − λ 2 χ 2 ). In the limiting case C → 0, shown in Fig. 5a, the above condition falls into √ µ = π/2 with the threshold being written explicitly as This result coincides with the harmonic case studied in [16], modulo a factor of two in γ 0 c . Note that in this case, the wavelength of the instability, L 0 x = 2π/q c = 2 2πK/B (Eq. 7), is given by the ratio of elastic constants and neither depends on the thickness h nor on the growth rate α. On the contrary, the threshold γ 0 c (Eq. 11), the minimum of the curve in Fig. 5a and Fig. 5c for C = 0, is lower in presence of the growth α = 0 and for higher thickness h of the cortex (small λ). This means that for thicker cortices the instability is more likely to happen, assuming the same elastic constants K/B and growth rate α. Now that we know γ 0 c , is this a physiologically accessible value? To answer this question, we can estimate the vertical stress exerted by axons at the white-grey matter interface, z = 0, required for the instability to happen since σ zz ≃ Bγ 0 c . The typical value of elastic modulus of cortical neurons is approximately 200 P a [6], which is related to B. The value of bending rigidity K is not available in the literature, however, based on the analogy with smectic liquid crystals we assume K/B ∼ 0.1 − 1 mm is of the order of the layer thickness. Thus, for a human brain with h ≃ 4 mm we can estimate λ ∼ 0.025 − 0.25 (assuming no growth), yielding the threshold γ 0 c ≃ 0.07 − 0.5 with necessary stress σ zz ∼ 10 − 100 P a, or the force of 10 − 100 pN per the unit area of 1 µm 2 . This prediction is consistent with tension measurements of neurons [14].
In general case, with (∂zφ) 2 |z =0 = C = 0, we are above the threshold in Eq. 11, and the axons should exert higher stresses than σ zz Bγ 0 c . Since the cortex is bounded below by the white matter, with its own inherent size (as opposed to being a pure elastic foundation of infinite size) and whose elastic properties are similar to the cortex [6], we do not expect ∂zφ|z =0 to vanish, but that its contribution be small such that C ≃ 0.1 − 0.2 potentially. In Fig. 5b, we plot threshold curves (Eq. 10) in the γ-χ plane for λ = 0.1 and different values of C. We cannot to identify the threshold (the minimum) in closed form, however, we can compute the values {γ c , χ c } as shown in Fig. 5c as function of the dimensionless variable, λ. As expected, γ c > γ 0 c , and the wavelength L x ∼ 1/χ c increases above the threshold. Note that Ref. [17] obtains the opposite trend. This difference is due to the additional (1 + γ) factors in our free energy, where we have not assumed that γ is small.
Before concluding this section, let us address this instability for different mammalian species. Assuming that elastic constants are of the same order for different species [20], we expect that for smaller h, λ increases, thus, the threshold, γ c , also increases. See Fig. 5c. Plausibly, in small species (small h) not enough force is generated by axons to overcome the threshold such that no folds emerge. On the other hand, the cortex thickness scales logarithmically with brain size [21] so that we should analyse the role of geometry and confinement on the instability threshold for layer buckling before drawing conclusions across species.

III. CORTEX AS A LIQUID CRYSTAL: POLAR GEOMETRY
Now, we model the cortex as 2D set of curves in polar coordinates (ϕ, r), confined between two radii r ∈ [R 0 ; R 0 + h] and ϕ ∈ [0; 2π]. Here, R 0 is the lateral size of an idealised circular brain, which varies among mammalian species from 5 mm to 20 cm [21]. The layers in the ground state are concentric circles radial with the position x 0 = re r , and the normal to the layers n = e r . Neglecting the effect of growth in this section, we consider a deformation map ϕ → ϕ and r → r + v(r, ϕ) with the corresponding deformation gradient Note that we do not incorporate growth here since it does not drive the instability found in the previous section. Then, the spatial gradient of isosurfaces in the current/target configuration can be computed as ∇ω = F −T e r , or The thickness of the deformed layers corresponds to l 0 /|∇ω|, where l 0 is the thickness of undeformed layers. The normal to the layers and its divergence, i.e. curvature, are given by where we have linearized n r ≃ 1 and n ϕ ≃ −∂ ϕ v/r for small deformations v ≪ 1. The free energy density in polar coordinates can now be obtained by inserting the above expressions into Eq.( 3). Looking for solutions in the form v(r, ϕ) = γr + ψ(r) cos(q ϕ ϕ) and integrating out the ϕ-dependence, keeping only terms quadratic in ψ, where y = log(r/R 0 ),ψ = ψ/R 0 , ξ 2 = K/(BR 2 0 ), and η = h/R 0 . The equilibrium equation forψ is a Schrödinger-type equation, describing a particle with energy E in the potential well V (y) (see Fig.6a), written as Assuming the WKB approximation for the classical region E > V (y) [22], we obtaiñ Satisfying the boundary conditionsψ WKB | y=0 = 0 and ∂ yψWKB | y=log(1+η) = 0, we find a relationship similar to Eq. 10, which reads as In fact, the first term can be integrated and cast in the closed form using More importantly, the curves defined by Eq. 18 have minima in the γ-q plane. These minima determine the threshold γ c . For the following analysis, we assume a logarithmic dependence of h ∈ [0.5 : 5] mm on the size of the brain R 0 ∈ [0.5 : 25] cm [21]. We observe that γ c decreases with the increasing system size R 0 , while the number of undulations q ϕc increases as shown in Fig. 6b. In other words, the smaller the brain, the less likely it is for cortex folds to develop given the increased threshold. Moreover, even if the threshold were met, the distance between folds would be larger such that there would be fewer folds. Our results may explain why mice brains do not exhibit folds, while human brains do. To be more precise, for a typical human brain with R 0 ∼ 10 cm, we obtain the typical dimensionful wavelength L x = 2πR 0 /q ϕc of the order of 1 cm for K/B ≃ 1 mm (ξ = 0.01) and 4 mm for K/B = 0.2 mm (ξ = 0.002, or 25 times smaller K) as shown in Fig. 6b with the colored sidebar. For a typical mouse brain, R 0 ∼ 1 cm such that the typical wavelength is about 6 mm (for ξ = 0.1). However, the critical strain required to initiate the instability exceeds unity such that the instability would not be accessible. We should also note that the value of γ c related to the stress exerted by the axons, σ rr ∼ Bγ c , and depends strongly on the thickness h of the cortex.

IV. DISCUSSION
Given the two pre-existing mechanisms of cortical folding (axonal tension and buckling), our cortex-as-a-smectic approach represents a novel way to think about the elasticity of the cortex. For the first time, we demonstrate that vertical pulling forces via axonal tension can lead to buckling. All prior buckling models of the cortex are a consequence of horizontal compression due to growth of the outermost layer of cortex. Our revised version of the axonal tension idea is in keeping with the observation that neurons in the white matter just beneath the cortex are oriented perpendicularly to the cortex [15]. While some doubt has been cast on the original version of the axonal tension model since circumferential tension along the axes of gyri (from one side of the "hill" to other side) is not observed, but radial tension is [9]. Our model does not conflict with this observation. The observation of circumferential tension near the bases of sulci (the valleys) presumably sets in a later stage in the folding process [9]. Here, we have focused on the intiation of the folding process.
Moreover, with our simple model, we obtain reasonable estimates of the critical wavelength and strain, in contrast to prior buckling estimates, where a large mismatch in the elastic moduli must be assumed to obtain agreement with observations [3]. The large mismatch (of order ten) does not agree with experiments [6]. A more recent buckling model with stress-dependent growth addresses this mismatch [23]. Since the white matter is modelled only as producing strain on the cortex via the underlying axons and boundary conditions, our result does not conflict with observations. On the other hand, our model does suggest the urgency for direct measurements of the bending rigidity of the cortex since our results suggest that K/B ∼ 1 mm. See Figure 7, where the critical wavelength is plotted as a function of cortex thickness for two different values of K/B for the planar and polar geometries presented here. Knowing the geometry of the human brain, Fig. 7a is more in keeping with observations than Fig. 7b.
By investigating the polar geometry, we address cortex folding across mammals. We find that smaller brains require a larger stress to initiate buckling/folding and that the critical wavelength increases with size. While we did not investigate the effect of growth in this geometry since the instability presented here is not driven by growth, it would be interesting to extend this case to include growth, particularly, differential growth, or stress-dependent growth [23].
Most of the cortical folds are simple folds-a simple "indentation", if you will, though some folds exhibit more structure. For instance, there exist secondary folds deeper inside the brain. We will call these more complicated structures, T-folds. And while such folds are more rare, a complete theory of cortex folding should be able to explain such emergent structures. To obtain these structures, growth, more details of the underlying white matter, and possibly constraints [24] will have to be incorporated into the model.
Finally, we have focused on the material properties of the cerebral cortex to better understand its shape. However, how does such properties affect its function? In developing a more accurate theory of the "brain as a material", can we better understand its function? For example, it would be interesting to couple viscoelasticity with connectivity models of the cortex [25] to determine more precisely the interplay between structure and function in the brain.