Decision letter: Exploring the role of the outer subventricular zone during cortical folding through a physics-based model

Abstract

Article Figures and data Abstract Editor's evaluation Introduction Cellular processes during brain development Model Model validation Results and discussion Conclusion Data availability References Decision letter Author response Article and author information Metrics Abstract The human brain has a highly complex structure both on the microscopic and on the macroscopic scales. Increasing evidence has suggested the role of mechanical forces for cortical folding – a classical hallmark of the human brain. However, the link between cellular processes at the microscale and mechanical forces at the macroscale remains insufficiently understood. Recent findings suggest that an additional proliferating zone, the outer subventricular zone (OSVZ), is decisive for the particular size and complexity of the human cortex. To better understand how the OSVZ affects cortical folding, we establish a multifield computational model that couples cell proliferation in different zones and migration at the cell scale with growth and cortical folding at the organ scale by combining an advection-diffusion model with the theory of finite growth. We validate our model based on data from histologically stained sections of the human fetal brain and predict 3D pattern formation. Finally, we address open questions regarding the role of the OSVZ for the formation of cortical folds. The presented framework not only improves our understanding of human brain development, but could eventually help diagnose and treat neuronal disorders arising from disruptions in cellular development and associated malformations of cortical development. Editor's evaluation Through theoretical analysis, the authors argue that the proliferation of neurons in the outer subventricular zone, which is specific to humans, decreases the distance between neighboring sulci in the cerebral cortex and increases cell density in the ventricular zone. Though the exact mechanisms remain to be further elucidated, the compelling data and approach represent a valuable foundation for the study of cortical folding from the underpinning cellular level as well as the coupling role of mechanics and cellular biology. This study will be of particular interest to the large community of scientists studying the mechanisms of brain development and disorder and even possibly beyond. https://doi.org/10.7554/eLife.82925.sa0 Decision letter Reviews on Sciety eLife's review process Introduction The brain is one of the most fascinating organs in the human body. Its complex structure on both micro- and macroscopic scales closely correlates with the unique cognitive abilities of humans. Cortical folding is one of the most important features of the human brain. Still, compared to other mammals, the human brain is neither the largest nor the most folded brain. However, relative to its size, it has the largest number of cortical neurons that connect with billions of neuronal synapses (Herculano-Houzel, 2009). This fact attracted the attention of neuroscientists over the past few years to explore the source of these cells and how they develop in the early stages of brain development. The number of brain cells is determined in utero through the proliferation process. Previous studies on different lissencephalic species, such as mice, have shown that cell division in the brain is confined to a small region near the cerebral ventricles (Hansen et al., 2010). However, in gyrencephalic species, this seems to be different. Recent findings show that the human brain, for example, is characterized by two proliferation zones with two different types of progenitor cells. Both zones produce neurons that later migrate towards the outer brain surface and form the cortex Lui et al., 2011; Pebworth et al., 2021. In rodents, progenitor cells around the ventricular zone (VZ) generate intermediate progenitor cells as their daughters, which accumulate above the VZ and form a new layer called the subventricular zone (Noctor et al., 2002). In humans, there is an additional outer layer of the subventricular zone, often referred to as outer subventricular zone (OSVZ) (Hansen et al., 2010; Lui et al., 2011; Noctor et al., 2007). This zone was first discovered in the monkey brain by Colette Dehay and her colleagues (Smart et al., 2002), and confirmed in the human brain by several following studies (Huttner and Kosodo, 2005). The OSVZ seems to play a significant role in the proliferation process and affects the size and complexity of the human cortex. The evidence for this allegation is the wave of cortical neurogenesis that coincides with the cell division in the OSVZ (Lukaszewicz et al., 2005). At the macroscopic scale, the high proliferation in the OSVZ coincides with a significant tangential expansion of the cortical layers. The latter is an essential factor for the formation of cortical folds (Reillo et al., 2011). Still, it remains unknown, how exactly this proliferation process in the OSVZ affects gyrification of the forming cortex. Different approaches have been used to understand the relation between cellular mechanisms at the microscopic scale and cortical development at the macroscopic scale. Genetic analyses and experimental studies using cell culture models and brain organoids have given first valuable insights concerning the source of cells and their behavior (Hansen et al., 2010). Here, we intend to complement these studies by using a numerical approach to bridge the scales from the behavior of different progenitor cell types at the cell scale to the emergence of cortical folds at the tissue or organ scale. From a mechanics point of view, forces that are generated due to cellular processes may act as a link to understand the underlying mechanisms behind cortical folding (Budday et al., 2015b). Many previous studies tried to explain normal and abnormal cortical folding either from a purely biological or mechanical perspective (Tallinen et al., 2014; Razavi et al., 2015). However, it will not be possible to capture the folding mechanism without considering both perspectives at the same time (de Rooij and Kuhl, 2018; Zarzor et al., 2021; Wang et al., 2022). In other words, to fully understand the physiological and pathological mechanisms underlying cortical folding in the developing human brain, we need to study the coupling between cellular processes and mechanical forces – to eventually assess how disruption of cellular processes affect the folding pattern and lead to malformations of cortical development (Guerrini et al., 2008; Blumcke et al., 2021; Llinares-Benadero and Borrell, 2019). To fill this knowledge gap, we establish a two-field computational model that accounts for both proliferating zones in the human brain, the VZ and OSVZ. The first field in the model describes the growth and deformation of brain tissue based on the theory of finite growth (Rodriguez et al., 1994; Göktepe et al., 2010). The second field describes the cellular processes occurring during human brain development, where we use an advection-diffusion equation to mimic the proliferation and migration in the subcortex and neuronal connectivity in the cortex (de Rooij and Kuhl, 2018; Zarzor et al., 2021). We add two source terms to consider the division in different zones. We validate our model through a comparison of the simulation results with histologically stained sections of the human fetal brain and address unresolved questions regarding the role of the OSVZ for cortical folding. Cellular processes during brain development The central cellular unit that plays a critical role in essential processes of brain development is a type of progenitor cells called radial glial cells. In the early stage of brain development, when neurogenesis begins, neuroepithelial cells transform into radial glial cells (Noctor et al., 2007). Around gestational week 5, these cells locate near cerebral ventricles, in the VZ, where they undergo interkinetic nuclear migration. The associated symmetric division behavior leads to a significant increase in the number of radial glial cells and results in both increased thickness and surface area of the VZ (Blows, 2003; Fish et al., 2008; Bystron et al., 2008). Subsequently, the cells switch to an asymmetric division behavior and generate intermediate progenitor cells (Noctor et al., 2004). The latter migrate to the subventricular zone, where they proliferate and produce neurons, as illustrated in Figure 1 (Noctor et al., 2007; Pebworth et al., 2021). Ultimately, the majority of cortical neurons are produced by intermediate progenitor cells (Lui et al., 2011; Libé-Philippot and Vanderhaeghen, 2021). Figure 1 Download asset Open asset Schematic illustration of human brain development between gestational weeks (GW) 4 and 38 at the cellular scale (top) and the organ scale (bottom). In the early stage of development, the repetitive division of radial glial cells in the ventricular zone (VZ) significantly increases the total number of brain cells. The newly born intermediate progenitor cells accumulate above the VZ and form a new layer called the inner subventricular zone. The outer radial glial cells (ORGCs) that are produced around gestational week 11 form a new layer called the outer subventricular zone (OSVZ). The neurons generated from progenitor cells migrate along radial glial cell fibers towards the cortex. Around gestational week 28, the migration process is almost finished, and the radial glial cells switch to produce different types of glial cells like astrocytes and oligodendrocytes. According to the radial unit hypothesis proposed by Pasko Rakic over 30 years ago, the radial glial cell fibers organize the migration process, which starts around gestational week 6 (Nonaka-Kinoshita et al., 2013). He postulated that these fibers form a scaffold to guide neurons during their migration from the proliferating zones to their final destination in the cortex, which forms the outer brain surface (Rakic, 1988; Lui et al., 2011). In gyrencephalic species, those fibers have a characteristic fan-like distribution (Borrell and Götz, 2014; Nonaka-Kinoshita et al., 2013). However, it is still under debate whether this unique distribution is the cause or rather the result of cortical folding. Our recent computational analyses support the latter, i.e., that it is the result of cortical folding (Zarzor et al., 2021). The migration process synchronizes with a radial expansion of all brain layers. Still, the VZ does not expand remarkably as the intermediate progenitor cells move outwards to the subventricular zone. The migrated neurons finally organize themselves in the six-layered cortex in an inside-out sequence, where the early-born neurons occupy the inner layers (Gilmore and Herrup, 1997). Until gestational week 23, the outer brain surface is still smooth, although the bottom four layers of the cortex are already filled with neurons (Shinmyo et al., 2017). The first folds appear between gestational weeks 20 and 28, as the cortical layer significantly expands tangentially (Budday et al., 2015b; Habas et al., 2012). Importantly, at around gestational week 25, the cortical neuronal connectivity emerges and comes along with the horizontal elongation of neuronal dendrites (Takahashi et al., 2012). While the processes summarized above are common among mammals, the human brain has some specific features that play a significant role in increasing the number of cortical neurons, and enhancing the complexity of cortical folds (Libé-Philippot and Vanderhaeghen, 2021). At the beginning of the second trimester, around gestational week 11, the original radial glial cells switch from producing intermediate progenitor cells to producing a special kind of cells that is found in all gyrencephalic species but is enriched in the human brain. The newly generated cells are similar to the radial glial cells in terms of shape and function, but unlike the original radial glial cells, they migrate to the outer layer of the subventricular zone, referred to as OSVZ, after they are born. Therefore, they are called outer radial glial cells (ORGCs) (Lui et al., 2011; Fietz et al., 2010; Hansen et al., 2010; Reillo et al., 2011; Nonaka-Kinoshita et al., 2013). We would like to note that some literature refers to this type of cells as basal radial glial cells. While the original radial glial cells have a bipolar morphology with two processes – one extending to the cerebral ventricle and one to the outer cortical surface – the ORGCs have a distinct unipolar structure with only a single process extending to the outer cortical surface (Hansen et al., 2010; Betizeau et al., 2013; Reillo et al., 2011; Nonaka-Kinoshita et al., 2013). The OSVZ shows a significantly more pronounced radial expansion compared to the inner subventricular zone and VZ between gestational weeks 11.5 and 32. The immediate reason causing this difference is the characteristic division behavior of ORGCs: they translocate rapidly in radial direction before they divide, which scientists have referred to as ‘mitotic small translocation (MST)’ (Fietz et al., 2010). Importantly, the MST behavior pushes the boundary of the OSVZ outward, which increases its capacity to produce new neurons. The intermediate progenitor cells have enough space to undergo multiple rounds of division before producing neurons, which increases the overall number of generated neurons (Kriegstein et al., 2006; Lui et al., 2011). The ORGCs, like radial glial cells, play an important role in the proliferation process: they divide symmetrically and asymmetrically to produce further ORGCs and intermediate progenitor cells (Libé-Philippot and Vanderhaeghen, 2021). Intermediate progenitor cells divide to generate a pair of neurons (Lui et al., 2011). According to previous studies, 40% of produced neurons are generated by ORGCs at gestational week 13, but this ratio increases to 60% by gestational week 14, and exceeds 75% by gestational week 15.5. Then, after gestational week 17, the ORGCs become the only source of cortical neurons in the upper cortical layers (Hansen et al., 2010). Besides their role in increasing the number of neurons, the ORGCs generate additional scaffolds that elongate to the outer brain surface and serve as paths for neuronal migration (Llinares-Benadero and Borrell, 2019; Nonaka-Kinoshita et al., 2013). Compared to other mammals, the neurogenesis of the human cortex can thus be divided into two main stages. The first stage is characterized by migration along a continuous scaffold consisting of radial glial cell fibers, which run from the ventricular surface to the outer cortical layer around gestational week 15. During the second stage, the migration path switches to a discontinuous form. After gestational week 17, the radial glial cell fibers run from the ventricular surface to the inner subventricular zone, while ORGC fibers run from the OSVZ to the outer brain surface, as indicated in Figure 1 (Nowakowski et al., 2016). Consequently, migrating neurons follow a sinuous path through numerous radial fibers before they reach its final location in the cortex (Lui et al., 2011). The additional scaffolds formed by ORGC fibers are not only important for neuronal migration, but also for the tangential expansion of the cortex. Previous studies show that a reduced number of ORGCs lead to reduced tangential expansion. In these cases, the cortex is less folded or even lissencephalic (Poluch and Juliano, 2015). In contrast, increasing the number of ORGCs leads to more excessive folding (Florio et al., 2017; Borrell, 2018). Still, it remains unknown whether these effects are a result of the specific proliferation behavior of ORGCs or the associated scaffold of ORGC fibers. What is known, though, is that the existence of ORGCs is a necessary but not sufficient condition for cortical folding (Llinares-Benadero and Borrell, 2019). Around gestational week 28, the migrating neurons occupy the first (top) cortical layer, while the migration and proliferation processes come to an end. After finishing their role during the neurogenesis stage, radial glial cells and ORGCs switch to produce different types of glial cells, e.g., astrocytes and oligodendrocytes (Schmechel and Rakic, 1979). Also, radial glial cells may later convert into ependymal cells that locate around the cerebral ventricles. The oligodendrocytes form myelin sheaths around neuronal axons, wherefore the subcortical layer gains its characteristic white color. Model To numerically study the effect of the VZ and the OSVZ on the resulting folding pattern, we simulate human brain development by using the finite element method. The influence of various factors on the emergence of cortical folds can be best shown on a simple two-dimensional (2D) quarter-circular geometry (Darayi et al., 2022), as illustrated in Figure 2A. In addition, we also investigate the folding evolution on a simplified half-sphere three-dimensional (3D) geometry. In the following, we introduce the main equations describing the coupling between cellular mechanisms in different proliferating zones and cortical folding, which we solve numerically. Figure 2 Download asset Open asset Kinematics of the multifield brain growth model. The reference configuration ℬ0 represents the initial state of the brain at gestational week (GW) 11. The spatial configuration ℬt represents the state of the brain at any time t during development. The stress-free (intermediate) growth configuration ℬg is inserted between reference and spatial configurations. (A) Simulation domain representing a part of the human brain’s frontal lobe. (B) Distribution of model parameters (r1c,r2c,v,anddcc) along the brain’s radial direction ri from the ventricular surface to the outer cortical surface. Kinematics To mathematically describe brain growth, we use the theory of nonlinear continuum mechanics supplemented by the theory of finite growth (Rodriguez et al., 1994; Göktepe et al., 2010). The initial state of the brain at an early stage of development, around gestational week 11, is represented by the reference configuration ℬ0. The state of the brain at time t later during development is represented by the spatial configuration Bt. The deformation map x=φ(X,t) maps a reference point X∈B0⊂R3 to its new position x∈Bt⊂R3 at a specific time t, as illustrated in Figure 2. The derivative of the deformation map with respect to reference point position vector is called deformation gradient F=∇Xφ. The local volume change of a volume element is described by the Jacobian J=detF. Following the theory of finite growth (Rodriguez et al., 1994; Göktepe et al., 2010), we introduce a stress-free configuration between the reference and spatial configuration, the growth configuration Bg. Accordingly, the deformation gradient is multiplicatively decomposed into an elastic deformation tensor Fe and a growth tensor Fg, such that, (1) F=Fe⋅Fg. The elastic deformation tensor describes the purely elastic deformation of the brain under the effect of external forces or forces generated internally to preserve tissue continuity. On the other hand, the growth tensor controls the amount and directions of unconstrained expansion. We note that the elastic deformation tensor is reversible, while the growth tensor is not. To not only predict brain growth but also its relation to cellular processes during brain development, we introduce the spatial cell density c(x,t) as an additional scalar independent field that depends on the spatial point position and time. It represents the number of neurons per unit area (de Rooij and Kuhl, 2018) in 2D and per unit volume in 3D. The corresponding balance equation describes cell division – resulting in newborn cells – through appropriate source terms and cell migration – the directed movement of neurons – through appropriate flux terms. For the two unknown fields, the deformation and the cell density, we introduce not only balance but also constitutive equations in the following that then allow us to compute their evolution in space and time through numerical simulations. We explain how we mathematically describe the mechanical (growth) problem and cellular processes as well as how those are linked to capture feedback mechanisms between cellular processes, mechanics, and growth. Mechanical problem To govern the mechanical problem, we use the balance of linear momentum given in the spatial configuration Bt, (2) ∇x⋅σ=0withσ=σ(Fe), where ∇x is the spatial gradient operator and σ is the Cauchy stress tensor formulated in terms of the elastic deformation tensor, as only the elastic deformation induces stresses. The Cauchy stress describes the 3D stress state in the spatial (grown and deformed) configuration and is computed by deriving the strain energy function ψg with respect to elastic deformation tensor, (3) σ(Fe)=1Je∂ψg(Fe)∂Fe⋅FeT, where Je=detFe. The strain energy function describes the material behavior of brain tissue mathematically. In our case, we consider a nonlinear hyperelastic material model as viscous effects, which have been observed for higher strain rates, become less relevant in the case of a slow process like brain development occurring over the course of weeks and months, as discussed in Budday et al., 2020. Our previous studies have shown that the isotropic neo-Hookean constitutive model best represents the material behavior of brain tissue during cortical folding (Budday et al., 2020). The corresponding strain energy function ψg is given as (4) ψg(Fe)=12λIn2(Je)+12μ(ri)[Fe:Fe−3−2In(Je)], where μ and λ are the Lamé parameters. We use the nonlinear Heaviside function that is given in the general form as H(x;γ)=eγx/(1+eγx) to guarantee a smooth transition from the cortex to the subcortical plate with distinct mechanical parameters, (5) μ(ri)=μs+[[μc−μs]×H(ri−rcp;20)], where the Heaviside function exponent γ equals 20. Please note that we will later use a different value of γ to serve numerical and geometrical requirements regarding the nature of the transition: higher values will lead to sharper transitions, lower values to smaller transitions. A more detailed discussion on the role of the value of γ can be found in Zarzor et al., 2021. Our recent numerical simulation study suggested that the cortical stiffness continuously changes during human brain development due to the changes in the local microstructure (Zarzor et al., 2021). Accordingly, we formulate the cortical shear modulus μc as a function of the cell density, (6) μc(c)={μ∞ifc≥cmax,μs+mc(c−cmin)ifcmax>c>cmin,μsifc≤cmin. It increases with increasing cell density in the range μc⁢(c)∈[μs,μ∞], while the subcortical shear modulus μs remains constant. The slope is defined as mc=μ∞-μs/cmax-cmin and the stiffness ratio as βμ=μ∞/μs. Through Equation 6, the cell density problem controls the effective stiffness ratio between cortex and subcortex (as the cortical stiffness changes while the subcortical stiffness remains constant) and thus also the emerging cortical folding pattern (Budday et al., 2014; Zarzor et al., 2021). Mechanical growth problem The growth tensor introduced in the Kinematics section is a key feature in our model that links the cell density problem with the mechanical problem. As it controls the amount and direction of growth, we need to consider how cellular processes affect the physiological growth behavior in order to find an appropriate formulation. During cellular migration, the subcortical layers expand isotropically. Then, under the effect of neuronal connectivity, the cortex grows – more pronounced in circumferential than in radial direction – as illustrated in Figure 2. Thus, we introduce the growth tensor as (7) Fg=ϑ⊥[I−N⊗N]+ϑ∥N⊗N, where N is the normal vector in the reference configuration B0 (it is linked to the spatial normal vector through N=F−1⋅n), while ϑ⊥ and ϑ∥ denote the growth multipliers in circumferential and radial direction, respectively. Those multipliers control the amount of growth as a function of the cell density, (8) ϑ⊥=[1+κ⊥(ri)c]αandϑ∥=[1+κ∥(ri)c]α, where κ⊥ and κ∥ are the growth factors in the circumferential and radial direction, respectively, and α is the growth exponent. To ensure isotropic growth in the subcortical layers, we formulate those factors as a function of the radius ri, such that (9) κ⊥(ri)=κs+[κs[βκ−1]×H(ri−rcp;20)]and (10) κ∥(ri)=κs+[κs[1βκ−1]×H(ri−rcp;20)], where κs is the growth factor in the subcortical layers, and βκ is the growth ratio between κ⊥ and κs (Zarzor et al., 2021). Through Equation 8, the amount of growth is directly related to the cell density – the higher the cell density, the more growth. Cell density problem We formulate the balance equation of the cell density problem in such a way that we can mathematically describe the different cellular processes occurring at the microscopic scale. Temporal changes in the cell density field are kept in balance by source and flux terms. The balance equation given in the spatial configuration ℬt follows as (11) J˙Jc+c˙=−∇x⋅[v^(c,x)c−dcc(x)⋅∇xc]+r1c(x,s)+r2c(x,s), where the first flux term v^(x)c represents the migration in the subcortical plate, the second flux term dcc(x)⋅∇xc represents the neuronal connectivity in the cortex, the first source term r1c(x,s) represents cell proliferation in the VZ, and the second source term r2c(x,s) cell proliferation in the OSVZ. The migration velocity vector v^(x) guides the cells along radial glial cell fibers and controls their speed, (12) v^(x)=H(c−c0;γc)v(ri)n/∥n∥. The vector n represents the normalized orientation of radial glial cell fibers in the spatial configuration and controls the migration direction of neurons. As the brain grows and folds, the fiber direction changes. Through this feedback mechanism, the mechanical growth problem affects how neurons migrate and the cell density evolves locally. The nonlinear regularized Heaviside function ℋ⁢(c-c0;γc) with the Heaviside exponent γc links the migration speed with the cell density field. Accordingly, the cells start to migrate only when their density exceeds the critical threshold c0. The value v specifies the maximum migration speed of each individual cell in the domain. To ensure that this value vanishes smoothly at the cortex boundary rcp, we formulate it as a function of the radial position ri, as shown in Figure 2B, such that (13) v(ri)=v[1−H(ri−rcp;10)]. After the cells reach the cortex, they diffuse isotropically, as described by the diffusion tensor dcc(x)=(dcc(ri)+vc(c)) I with the diffusivity dcc, the artificial viscosity vc(c), and the second order unit tensor I. The artificial viscosity term vc(c) serves as a numerical stabilization to avoid numerical oscillations associated with the advection-diffusion equation. It only acts when the actual cell density does not satisfy the balance equation and ensures more reliable results without having a particular physical meaning. The diffusivity is defined as a function of the radial position ri to act only in the cortex, (14) dcc(ri)=dccH(ri−rcp;10). The first source term r1c represents the radial glial cell proliferation in the VZ, as demonstrated in Figure 2B, and is given as (15) r1c(x,s)=Gvzs(s)[1−H(ri−rvz;50)]with (16) G{∙}s(s)=G{∙}−{(s−1)G{∙}ifs<1.80.8G{∙}else where rvz is the outer radial boundary of the VZ and Gvzs is the division rate in the VZ as a function of the maximum stretch s in the domain. By applying Equation 16 for the VZ, we ensure that the division rate decreases from its initial value Gvz to a smaller value as the maximum stretch value s in the domain increases, i.e., with increasing gestational age. This constitutes an additional feedback mechanism between the mechanical growth problem and the cell density problem: As the maximum stretch and thus the deformation increases due to constrained cortical growth, the division rate in the VZ decreases, resulting in less newborn cells. Besides the proliferation of radial glial cells around the cerebral ventricles in the VZ, the ORGCs proliferate in the OSVZ. To capture this effect, we add a second source term r2c, as demonstrated in Figure 2B. The second source term is given as (17) r2c(x,s)=Gosvzs(s)[H(ri−risvz;50)−H(ri−rosvz(t);50)], where risvz is the outer radial boundary of the inner subventricular zone and Gosvzs is the division rate in the OSVZ that again decreases with increasing maximum stretch s in the domain. To numerically capture the expansion of the OSVZ under the effect of MST of ORGCs, we formulate the outer radial boundary of the OSVZ as a function of time, such that,rosvz=risvz+mmstt, where mmst is introduced as the MST factor. Again, we apply Equation 16 for the OSVZ, but in this case with the initial division rate G{osvz}. Model parameters and boundary conditions In this work, we will consider two different cases regarding the mechanical model: The first case considers a varying cortical stiffness as introduced in the Mechanical problem section, while the second case assumes a constant cortical stiffness, i.e., μc=μ∞=constant. While our previous study had suggested that the simulations with varying cortical stiffness lead to morphologies that better agree with those in the actual human brain (Zarzor et al., 2021), we still consider both cases in the following, varying stiffness and constant stiffness, as the situation might change when including the OSVZ and we aim to investigate corresponding interdependency effects. Table 1 summarizes the model parameters that are used in the simulation for the 2D case. We will refer to the parameters changes in the 3D case later when we present the corresponding results to avoid confusion. The mechanical and diffusion parameters are adapted from the literature (Budday et al., 2020; de Rooij and Kuhl