Biological backgroundDuring human brain development, four fundamental cell types play an essential role, radial glial cells (RGCs), intermediate progenitor cells (IPCs), Outer radial glial cells (ORGCs), and neurons, as illustrated in Fig. 1, left. At the very beginning of neurogenesis, the original stem cells, which are called neuroepithelial cells, transform into RGCs in the brain’s ventricular zone around the cerebral ventricles15. Between gestational weeks 5 and 7, the RGCs amplify their number through symmetric division that leads to a significant increase in both the thickness and surface area of the ventricular zone (VZ)16,17,18. Thereafter, the newly generated RGCs switch from symmetric to asymmetric division and become the initial source of IPCs19. The IPCs migrate to populate the subventricular zone (SVZ), located above the ventricular zone (VZ). In the SVZ, IPCs undergo multiple rounds of division before they ultimately produce cortical neurons20,21. Importantly, apart from the limited quantity of neurons generated by RGCs during the initial stages of neurogenesis, IPCs serve as the exclusive source of cortical neurons3,22. In gyrencephalic species, especially in humans, by around gestational week 11, RGCs undergo a transition from generating IPCs to generating ORGCs23,24,25. The latter exhibit a distinctive division behavior, swiftly translocating upon generation to the outer layer of the SVZ, forming a distinct zone known as the outer subventricular zone (OSVZ). This unique division behavior is reffed to as ’mitotic somal translocation’ (MST)24,26,27. As a result of this behavior, the boundary of the OSVZ expands outward in radial direction over time, which increases its capacity to produce new cells21. During gestational weeks 11 to 24, ORGCs undergo asymmetric division to produce IPCs, making ORGCs the second source of these cells22. The abundance of IPCs in gyrencephalic species plays a crucial role in doubling the neurons generated. Consequently, these species exhibit a significantly greater number of cortical neurons compared to lissencephalic species24. The newly generated neurons migrate towards the cortical layer, arranging themselves in a sequence from inside to outside to form the six-layered cortex28. According to the radial unit hypothesis, the fibers of both RGCs and ORGCs form a scaffold for neuronal migration29. Therefore, both types of radial glial cells function as central units for both cell generation and the organization of neuronal migration. After gestational week 23, brain folds start to emerge as the outer brain surface significantly expands circumferentially5,30, as indicated in Fig. 1, right. Intriguingly, this expansion of the brain’s surface coincides with the development of neuronal connectivity31.Fig. 1Kinematics of human brain growth model between gestational weeks (GW) 18 and 38. The deformation map \({\varvec{{x}}}=\varvec{\varphi }({\varvec{{X}}}, t)\) maps the material point \({\varvec{{X}}}\) from the stress-free reference configuration \(\mathcal {B}_\text{0}\) with outward pointing normal vector \({\varvec{{N}}}\) to the spatial configuration \(\mathcal {B}_\text{t}\) with outward pointing normal vector \(\varvec{n}\) at a specific time t. Inserting the intermediate growth configuration \(\mathcal {B}_\text{g}\) leads to a multiplicative decomposition of the deformation gradient \({\varvec{{F}}} = \varvec{\nabla }_{{\varvec{{X}}}} \varvec{\varphi }\) into an elastic component, \(\varvec{F}^\text{e}\), and a growth component, \({\varvec{{F}}}^\text{g}\). The Jacobian \(J= \text {det} \, {\varvec{{F}}}\) maps the cell density c from the reference to spatial configuration, i.e., \(c_0 = Jc\). The cellular processes on the microscopic scale drive the cortical \(\vartheta ^\parallel\),\(\vartheta ^\bot\), and subcortical \(\vartheta _{\text{s}}\) growth multipliers in \(\mathcal {B}_\text{g}\).Computational modelTo numerically study the role of the different cell types mentioned in the previous section on the resulting folding patterns during human brain development, we develop a multi-field computational model coupling growth and mechanical deformation with a cell-density field. The cell density field is characterized by a system of advection-diffusion equations, each formulated to mimic the distinct behavior of a specific cell type. The mechanical deformation is coupled to the cell density field as growth is controlled by the cell density. In the subsequent sections, we present the fundamental equations that outline this coupled problem. Initially, we provide a broad overview of the kinematics of finite growth, the mechanical balance, and the constitutive equations. Following this, we delve into the specifics of the cell-density problem. Next, we introduce the discretization process in space and time, and finally, we discuss the selection and calibration of model parameters.KinematicsIn our approach, we use the principles of nonlinear continuum mechanics, enriched with concepts from the theory of finite growth. On the macroscopic (tissue) scale, we define the deformation map \({\varvec{{x}}} = \varvec{\varphi } ({\varvec{{X}}}, t)\) to depict the transition of a material point from the brain’s initial (reference) state, \(\mathcal {B}_{0} \subset \mathbb {R}^3\), to its spatial (deformed) state, \(\mathcal {B}_\text{t} \subset \mathbb {R}^3\), at a subsequent time \(t \in \mathbb {R}_{+}\) during its development, with \({\varvec{{x}}} \in \mathcal {B}_\text{t}\) and \({\varvec{{X}}} \in \mathcal {B}_0\). Additionally, we introduce the deformation gradient, \({\varvec{{F}}} = \varvec{\nabla }_{{\varvec{{X}}}} \varvec{\varphi }\), a derivative of the deformation map with respect to the position vector of the reference point, to transform a line element from its original to its new position. The determinant of the deformation gradient, the Jacobian \(J= \text {det} \, {\varvec{{F}}}\), is used to measure the resulting local volumetric change, as illustrated in Fig. 1. Consistent with the theory of finite growth, an intermediate (stress-free) configuration, \(\mathcal {B}_\text{g}\), is inserted between the brain’s two primary configurations. Consequently, the deformation gradient is multiplicatively decomposed into an elastic component, \({\varvec{{F}}}^\text{e}\), and a growth component, \({\varvec{{F}}}^\text{g}\)32,33, such that$$\varvec{F} = \varvec{F}^{{\text{e}}} \cdot \varvec{F}^{{\text{g}}} .$$
(1)
Likewise, the Jacobian is split into \(J = J^\text{e} \cdot J^\text{g}\), where \(J^\text{e} = \text {det} \, {\varvec{{F}}}^\text{e}\), and \(J^\text{g} = \text {det} \, {\varvec{{F}}}^\text{g}\). It should be highlighted that the elastic deformation tensor captures the brain’s purely elastic response to both external and internal forces, which preserves the continuity of the tissue. In contrast, the growth tensor specifies the extent and direction of free expansion led by the cell-density field. Notably, the elastic deformation is reversible, whereas the growth part leads to an irreversible change in the system.To represent the microscopic (cell) scale, we introduce the independent scalar field of material cell density, \(c_0({\varvec{{X}}})\), which represents the total number of cells per unit volume of the brain in its reference state. The spatial cell density, \(c({\varvec{{x}}},t)\), defined for a specific time t in the deformed configuration, is linked to the material cell density through the Jacobian as \(c_0 = Jc\).Mechanics problemFor the mechanics problem, we solve the balance of linear momentum in the absence of external forces, which is represented in the spatial configuration as,$$\begin{aligned} \text {div} \,\varvec{\sigma } = {\textbf {0}}, \end{aligned}$$
(2)
where \(\text {div}\) is the divergence operator with respect to spatial position \({\varvec{{x}}}\), and \(\varvec{\sigma }\) is the Cauchy stress tensor. The stress tensor is derived from the elastic deformation tensor since only elastic deformations contribute to stress within the material, i.e., \(\varvec{\sigma } = \varvec{\sigma }({\varvec{{F}}}^\text{e})\).We assume that the response of brain tissue during development is isotropic and hyperelastic, as viscous effects, notable at higher strain rates, diminish in significance over the developmental stages that span months of gestation. Moreover, within the specific regions of interest, i.e., the subcortex and cortex, the tissues have shown isotropic properties in experiments2,34. Mathematically, we introduce the strain energy function \(\psi _0\) in the spatial configuration \(\mathcal {B}_{0}\), that is corresponding to function \(\psi _{\text{g}}\) in the growth configuration \(\mathcal {B}_{g}\), where \(\psi _0 = J^\text{g} \psi _\text{g}\). Based on our previous work2, we adopt a neo-Hookean strain energy function to describe the material behavior of brain tissue during cortical folding. This is formulated as follows,$$\begin{aligned} \psi _{\text{g}}({\varvec{{F}}}^\text{e}) =&\frac{1}{2} \, \lambda \, \text {ln}^{2}(J^\text{e}) + \frac{1}{2} \, \mu \, \left[ \, {\varvec{{F}}}^\text{e}:{\varvec{{F}}}^\text{e} -3 -2 \text {ln}(J^\text{e}) \right] , \end{aligned}$$
(3)
where \(\mu\) and \(\lambda\) are the Lamé parameters that relate to the Poisson’s ratio \(\nu\) in the linear regime as \(\nu = \lambda / 2 [\lambda + \mu ]\).According to our previous findings35, the human brain tissue tends to have a varying stiffness during development due to the changes in its local microstructure36,37. Consequently, we propose that the cortical shear modulus \(\mu _{\text{c}}\) increases linearly from the subcortical shear modulus (\(\mu _{\text{s}} = \text {constant}\)) to that of a fully developed human cortex, \(\mu _{\infty }\), such that,$$\begin{aligned}&\displaystyle {\mu _{\text{c}}(c) = \text {min} \left( \, \mu _{\infty } \, , \,\mu _{s} + m^{\mu }_{\text{c}} \, \langle (c-c^{\mu }_{\text {min}} \,) \rangle \right) } \\&\text {with} \qquad \displaystyle {m^{\mu }_{\text{c}} = \frac{\mu _{\infty } – \mu _{\text{s}}}{c^{\mu }_{\text {max}} – c^{\mu }_{\text {min}}}}.\nonumber \end{aligned}$$
(4)
where \(c^{\mu }_{\text {min}}\) is the minimum threshold, \(c^{\mu }_{\text {max}}\) is the maximum threshold, and \(\langle \rangle\) are Macaulay brackets that return its argument if the argument is positive and zero otherwise as shown in Fig. 2.Fig. 2An idealized evolution of the cell-density- dependent cortical shear modulus \(\mu _{\text{c}}\) during brain development, assumed to increase linearly from subcortical shear modulus \(\mu _{\text{s}}\) to its final value corresponding to the fully developed cortex \(\mu _{\infty }\).Adhering to the Clausius-Planck inequality, which presumes no internal energy dissipation within the hyperelastic materials38, we deduce the Cauchy stress as follows:$$\begin{aligned} \varvec{\sigma }&= \frac{1}{J^\text{e}} \frac{\partial \psi _\text{g}({\varvec{{F}}}^\text{e})}{ \partial {\varvec{{F}}}^\text{e}} \cdot {{\varvec{{F}}}^\text{e}}^{T} = \frac{1}{J^\text{e}} \left[ \left[ \lambda \, \text {ln}(J^\text{e}) – \mu \right] {\varvec{{I}}} + \mu \, {\varvec{{b}}}^\text{e} \right] , \end{aligned}$$
(5)
where \({\varvec{{b}}}^\text{e} = {{\varvec{{F}}}^\text{e}} \cdot {{\varvec{{F}}}^\text{e}}^{T}\) is the elastic left Cauchy-Green tensor, and \({\varvec{{I}}}\) is the second order unit tensor.Mechanical growth problemThe growth tensor \({\varvec{{F}}}^\text{g}\) is the key feature of the model, linking the mechanics problem with the cell-density problem. It thus plays a crucial role in controlling tissue deformation. We establish the mathematical formulation of the growth tensor based on the hypothesis of differential growth and guided by our understanding of cellular mechanisms5,39,40. During the early stages of human brain development, cellular processes, such as cell division, translocation of progenitor cells, and neuronal migration, all contribute to the subcortex’s isotropic expansion. As development continues, neurons settle in the cortex and start synapse formation, the development of neural connections results in more pronounced circumferential growth compared to radial growth14,31. This observation, which aligns with the differential growth hypothesis, is demonstrated in Fig. 1. Accordingly, we introduce independent growth multipliers \(\vartheta ^\bot\) and \(\vartheta ^\parallel\) in the circumferential and radial direction, respectively, as:$$\begin{aligned} \vartheta ^\bot = \left[ 1+ \kappa ^\bot \, c \right] ^\alpha \quad \text {and} \quad \vartheta ^\parallel = \left[ 1+ \kappa ^\parallel \, c \right] ^\alpha , \end{aligned}$$
(6)
where the growth factors \(\kappa ^\bot\), \(\kappa ^\parallel\), and the growth exponent \(\alpha\) determine the strength and shape of the coupling between the cell-density field and the mechanical field. Those multipliers are identical in the subcortex, denoted as \(\vartheta ^\bot = \vartheta ^\parallel = \vartheta _{s}\), but differ in the cortex, where \(\vartheta ^\bot> \vartheta _{s} > \vartheta ^\parallel\), as explained in more detail in Sect. “Model parameters”. Consequently, the growth tensor is written as,$$\begin{aligned} {\varvec{{F}}}^\text{g} = \vartheta ^\bot \; \left[ {\varvec{{I}}} – {\varvec{{N}}} \otimes {\varvec{{N}}} \right] + \vartheta ^\parallel \; {\varvec{{N}}} \otimes {\varvec{{N}}}, \end{aligned}$$
(7)
where \({\varvec{{N}}}\) is the outward pointing normal vector in the reference configuration \(\mathcal {B}_{0}\), linked to the spatial normal vector through \({\varvec{{N}}} = {\varvec{{F}}}^{-1} \cdot {\varvec{{n}}}\)12,35.Cell-density problemTo describe the behavior of each cell type involved in brain development, we additively decompose the cell-density field into four separate components. These components correspond to radial glial cells (RG), intermediate progenitor cells (IP), outer radial glial cells (ORG), and neurons (N), respectively, such that,$$\begin{aligned} c = c_{\,\text {\tiny RG}} + c_{\,\text {\tiny IP}} + c_{\,\text {\tiny ORG}} + c_{\,\text {\tiny N}}. \end{aligned}$$
(8)
To cover the cell-density problem, we introduce a system of mass balance equations in the spatial configuration \(\mathcal {B}_{t}\). Each equation is formulated to mimic the distinct behavior of a particular cell type. This system of equations is effectively summarized in a single mathematical expression as follows,$$\begin{aligned} \frac{\dot{J}}{J} \, c_{ \, \bullet } + \dot{c}_{\, \bullet } = \text {div} \, {\varvec{{q}}}_{\, \bullet } ({\varvec{{x}}}, \, c_{\, \bullet }) + f_{\, \bullet } (t, \, c_{\, \circ }), \end{aligned}$$
(9)
where \(\bullet\) serves as a placeholder for the specific cell type in question, \({\varvec{{q}}}_{\, \bullet }\) is the flux term, and \(f_{\, \bullet }\) is the source term in the spatial configuration, with \(\circ \in \{\text { RG, IP, ORG, N}\}\).The flux term characterizes cell motion and is constructed as the sum of two distinct components: the advection term and the diffusion term, which are generally expressed as,$$\begin{aligned} {\varvec{{q}}}_{\, \bullet }=&-\underbrace{c_{\, \bullet } \hat{{\textbf {v}}}_{\, \bullet }({\varvec{{x}}}, c_{\, \bullet })}_{\text {advection term}} +\underbrace{{\varvec{{d}}}^{\, \mathrm c}_{\, \bullet } ({\varvec{{x}}}, c_{\, \bullet }) \cdot \varvec{\nabla }_{{\varvec{{x}}}} \,c_{\, \bullet }}_{\text {diffusion term}}. \end{aligned}$$
(10)
In this formulation, the advection term describes cell transport, including translocation (such as for ORGCs), migration (as seen with neurons), or more regular movements (notably for IPCs). Meanwhile, the diffusion term describes cell dispersion, occurring either in the subcortex (for progenitor cells) or in the cortex (for neurons). The dynamics of transport and dispersion are cell-type specific and can be adjusted through certain parameters, which will be elaborated upon in Sect. “Model parameters”. In the advection term, the transport velocity vector \(\hat{{\textbf {v}}}\) describes cell movement in radial direction, following the path of radial glial cell fibers towards the increasingly distant outer brain surface due to the growth41,$$\begin{aligned} \hat{{\textbf {v}}}_{\, \bullet }({\varvec{{x}}}, c_{\, \bullet }) = \mathcal {H}(c_{\, \bullet } – c_{0}; \gamma ^\text{c} ) \, v_{\, \bullet }({\varvec{{x}}}) \, \frac{{\varvec{{n}}}}{\Vert {\varvec{{n}}} \Vert }, \end{aligned}$$
(11)
where \(v_{\, \bullet }\) defines the speed of transport, and the current norm vector \({\varvec{{n}}}\) determines the direction of transport. The latter depends on the mechanical deformation through the push forward operator, i.e., \({\varvec{{n}}} = {\varvec{{F}}} \cdot {\varvec{{N}}}\). Here, the term \(1/{\Vert {\varvec{{n}}} \Vert }\) serves to normalize the direction vector. Finally, we use the nonlinear regularized Heaviside function \(\mathcal {H}(c_{\, \bullet } – c_{0}; \gamma ^c )\), with the threshold \(c_{0}\) and exponent \(\gamma ^c\), to model the impact of adaptive transport based on type-specific cell density. It ensures a smooth increase in the amount of transport as the cell density rises. The general formulation of this function is given by,$$\begin{aligned} \mathcal {H}(\circ ; \gamma ) =\frac{e^{\gamma \, \circ }}{ 1 + e^{\gamma \, \circ }}. \end{aligned}$$
(12)
Assuming isotropic diffusion, the diffusion tensor may be written as,$$\begin{aligned} {\varvec{{d}}}^{\, \mathrm c}_{\, \bullet }({\varvec{{x}}}, c_{\, \bullet }) = \left[ {d}^{\, \mathrm c}_{\, \bullet } ({\varvec{{x}}})+ \nu _{\, \bullet }^{\, \mathrm c}(c_{\, \bullet }) \right] \, {\varvec{{I}}} \end{aligned}$$
(13)
where \({d}^{\, c}_{\, \bullet }\) is the type-specific diffusivity that is a function of the spatial position, and \(\nu _{\, \bullet }^{\, c}\) is an artificial (numerical) diffusivity. Although this artificial diffusivity has no physical meaning, it is intended to ensure numerical stability within the advection-diffusion equation, addressing the disparity between the advection and diffusion terms in the transition area between the cortex and subcortex at the onset of transition. To compute this value, we follow the streamline upwind Petrov-Galerkin (SUPG) method. For more computational details, we refer to42.Based on the lineage relationships between the different cell types under consideration, we define the source term in Eq. (9) as follows,$$\begin{aligned} f_{\, \bullet } (t, \, c_{\, \circ }) = \sum _{\begin{array}{c} \circ \in \{ \text {\tiny RG, IP, ORG, N} \} \end{array}} G_{ \, \bullet } (\circ , \text {P}(t)) \, c_{\, \circ }, \end{aligned}$$
(14)
where the summation \(\sum\) extends over the four cell types, \(G_{ \, \bullet } (\circ , \text {P})\) is a division ratio between the cell types \(\bullet\) and \(\circ\) in the particular division phase \(\text {P}\), and \(c_{\, \circ }\) denotes the density of cell type \(\circ\). More details regarding the source term will be given in the model parameter Sect. “Model parameters”.Discretization in timeSince the system of cell density balance equations given in 9 includes two time-dependent components, namely the Jacobian and the type-specific cell density, time discretization is essential. In our approach, we have, in addition to the cell density field, the nonlinear deformation field that requires a nonlinear solver. Therefore, we can employ the nonlinear solver to implicitly calculate the time-dependent terms using an Euler discretization scheme13. We subdivide the time domain \(\mathcal {T}\) into \(\mathrm n_{\, \mathrm st}\) intervals such that,$$\begin{aligned} \mathcal {T} = \bigcup _{\mathrm{n =1}}^{{\text{n}}_{\, \text{st}}} \left[ t_\mathrm{n-1} , t_{\text{n}}\right] , \end{aligned}$$
(15)
and obtain,$$\begin{aligned} \dot{c}_{ \, \bullet }= \displaystyle {\frac{c_{{\, \bullet }, \, \text{n}}- c_{{\, \bullet }, \, \mathrm{n-1}}}{\Delta t}}, \qquad \text {and} \qquad \dot{J} = \displaystyle {\frac{J_{\text{n}}-J_{\mathrm{n-1}}}{\Delta t}}, \end{aligned}$$
(16)
where \(c_{{\, \bullet }, \, \text{n}}\) and \(J_{\text{n}}\) are the unknown values at the current time step, \(c_{{\, \bullet }, \, \mathrm{n-1}}\) and \(J_{\mathrm{n-1}}\) are the known values from the previous time step, and \(\Delta t = t_{\text{n}}-t_{\mathrm{n-1}}\) is the time increment. As a result, we derive the fully time-discretized system of equations for the cell density balance as follows,$$\begin{aligned}&\displaystyle {\frac{J_{\text{n}}-J_{\mathrm{n-1}}}{J \, \Delta t}} \, c_{ \, \bullet , \, \mathrm n} + \displaystyle {\frac{c_{{\, \bullet }, \, \text{n}}- c_{{\, \bullet }, \, \mathrm{n-1}}}{\Delta t}} = \text {div} \, {\varvec{{q}}}_{\, \bullet }({\varvec{{x}}}, \, c_{{\, \bullet }, \, \text{n}}) + f_{ \, \bullet }(t_{\text{n}}, \, c_{{\, \bullet }, \, \mathrm{n-1}}). \end{aligned}$$
(17)
The source term is approximated explicitly as,$$\begin{aligned}&f_{\, \bullet } (t_{\text{n}}, \,c_{{\, \circ }, \, \mathrm{n-1}}) = \sum _{\begin{array}{c} \circ \in \{ \text {\tiny RG, IP, ORG, N} \} \end{array}} \Delta t \, G_{ \, \bullet } (\circ , \text {P}(t_{\text{n}})) \, c_{\, \circ , \, \mathrm{n-1}}. \end{aligned}$$
(18)
Geometry and discretization in spaceTo systematically study the interaction between the proliferation of different cell types in the developing human brain and cortical folding, we chose to use a two-dimensional domain that represents an exemplary part of a coronal slice at gestational week 24, as shown in Fig. 3. The selected domain facilitates comparison with magnetic resonance (MR) images and histologically stained sections of fetal brains (HBS)35, while simultaneously lowering computational cost, making it an optimal choice for our application. The following Dirichlet and Neumann boundary conditions are applied on the domain,$$\begin{aligned}&\varvec{\varphi }({\varvec{{x}}}) = {\varvec{{0}}}&\text {for} \hspace{25pt} {\varvec{{x}}}\in \partial \varvec{\Omega }^{\varphi } \nonumber \\&\varphi _2 ({\varvec{{x}}})= 0&\text {for} \hspace{25pt} {\varvec{{x}}}\in \partial \varvec{\Omega }^\text{t}_1 \nonumber \\&\varphi _1 ({\varvec{{x}}})= 0&\text {for} \hspace{25pt} {\varvec{{x}}}\in \partial \varvec{\Omega }^\text{t}_2 \nonumber \\&\varvec{\sigma }\cdot {\varvec{{n}}}={\varvec{{0}}}&\text {for} \hspace{25pt} {\varvec{{x}}}\in \partial \varvec{\Omega }^\text{t}_3 \nonumber \\&{\varvec{{q}}}_{\, \bullet }\cdot {\varvec{{n}}}={\varvec{{0}}}&\text {for} \hspace{25pt} {\varvec{{x}}}\in \partial \varvec{\Omega }^\text{t}_3 \end{aligned}$$
(19)
where \(\partial \varvec{\Omega }^\varphi\) is the Dirichlet boundary and \(\partial \varvec{\Omega }^\text{t}\) are the Neumann boundaries satisfying \(\partial \varvec{\Omega }^{ \mathrm t} \cup \partial \varvec{\Omega }^\varphi = \partial \varvec{\Omega }\) and \(\partial \varvec{\Omega }^\text{t}\cap \partial \varvec{\Omega }^\varphi = \emptyset\). Here, \(\partial \varvec{\Omega }^\text{t}_1 = \{ {\varvec{{x}}} \in \partial \varvec{\Omega }^\text{t} \, \big | \, x_2 = 0 \}\) and \(\partial \varvec{\Omega }^\text{t}_2 = \{ {\varvec{{x}}} \in \partial \varvec{\Omega }^\text{t} \, \big | \, x_1 = 0 \}\), where \(\varphi _1\) and \(\varphi _2\) are the first and second component of the deformation map. In addition, we segment the simulation domain into five neighboring sections, representing the germinal zones within the brain (see Fig. 3). Mathematically, each zone is distinguished by a specific radius; \(r_{\text {vz}}\) represents the upper limit of the ventricular zone, \(r_{\text {isvz}}\) represents the upper boundary of the inner subventricular zone, \(r_{\text {osvz}}\) represents the upper and final boundary of the outer subventricular zone, while \(r_{\text {cp}}\) represents the lower boundary of the cortex. Finally, we discretize the geometry into a regular finite element mesh and linearize the nonlinear coupled problem using the Newton-Raphson method. Accordingly, we end up with the following global system of algebraic equations that can be solved numerically,$$\begin{aligned} \underbrace{\begin{bmatrix} \textbf{K}^{\varphi \, \varphi } & \textbf{K}^{\varphi \, \text {\tiny RG}} & \textbf{K}^{\varphi \, \text {\tiny IP}} & \textbf{K}^{\varphi \, \text {\tiny ORG}} & \textbf{K}^{\varphi \, \text {\tiny N}}\\ \textbf{K}^{ \text {\tiny RG} \, \varphi }& \textbf{K}^{\text {\tiny RG} \, \text {\tiny RG}} & 0 & 0& 0\\ \textbf{K}^{ \text {\tiny IP} \, \varphi } & 0& \textbf{K}^{{\text {\tiny IP}} \,{\text {\tiny IP}}} & 0& 0\\ \textbf{K}^{ {\,\text {\tiny ORG}}\, \varphi } & 0 & 0& \textbf{K}^{{\text {\tiny ORG}}{\,\text {\tiny ORG}}} & 0\\ \textbf{K}^{ {\text {\tiny N}} \, \varphi } & 0 & 0& 0& \textbf{K}^{{\text {\tiny N}}\,{\text {\tiny N}}} \\ \end{bmatrix}}_{\textbf{K}} \cdot \underbrace{\begin{bmatrix} \Delta \, \varvec{\varphi }\\ \Delta \, c_{\,\text {\tiny RG}} \\ \Delta \, c_{\,\text {\tiny IP}} \\ \Delta \, c_{\,\text {\tiny ORG}} \\ \Delta \, c_{\,\text {\tiny N}} \end{bmatrix}}_{\Delta \, \varvec{\phi }} = \underbrace{\begin{bmatrix} \textbf{R}^{\varphi }\\ \text {R}^{c_{\,\text {\tiny RG}}}\\ \text {R}^{c_{\,\text {\tiny IP}}}\\ \text {R}^{c_{\,\text {\tiny ORG}}}\\ \text {R}^{c_{\,\text {\tiny N}}} \end{bmatrix}}_{\textbf{R}}, \end{aligned}$$
(20)
where \(\textbf{K}\), \(\Delta \, \varvec{\phi }\), and \(\textbf{R}\) are the total tangent matrix, total incremental solution vector, and right-hand side residual vector, respectively. Importantly, the zeros in the tangent matrix indicate the absence of direct coupling (implicit) between the different cell density components. Rather, the coupling is addressed indirectly (explicitly) via the source terms, which appear on the right-hand side in the residual vector.Fig. 3The simulation domain represents a 2D exemplary part of a coronal slice at gestational week 24 with the outer radius R and ventricular radius r. The schematic provides various views of the domain through representations. The continuum representation demonstrates the Dirichlet \(\partial \varvec{\Omega }^\varphi\) and Neumann boundaries \(\partial \varvec{\Omega }^\text{t}\). The germinal zones representation highlights the spatial locations of these zones within the simulation domain. The mesh representation illustrates the initial mesh of the domain.Model parametersThe developmental dynamics of the human brain and its heterogeneous microstructure lead to a non-uniform distribution of model parameters throughout the simulation domain, i.e., the shear modulus, growth factors, amount of transport, and diffusivity. These parameters vary significantly across the different brain zones2, introduced in Fig. 1. Furthermore, the real human brain does not show distinct boundaries between its zones. To address this issue, on the one hand, and ensure numerical stability, on the other hand, we use the Heaviside function, introduced in equation (12) to model a smooth transition at the interfaces between zones. Here, we define the Heaviside function as a function of the radius \(r_{i}\) of a material point \({\varvec{{X}}}\) relative to the geometric center \({\varvec{{X}}}_0\), calculated as \(r_{i} = \parallel {\varvec{{X}}} – {\varvec{{X}}}_0 \parallel _{2}\).Given that the cortex is \(\beta _{\mu }\) times stiffer than the subcortex, we define the transition of the shear modulus at the cortex-subcortex boundary \(r_{\text {cp}}\), as,$$\begin{aligned} \mu (r_{\text{i}}) = \mu _{\text{s}}+\left[ \mu _{\text{s}} \left[ \beta _{\mu}-1 \right] \times \mathcal {H} (r_{\text{i}} – r_{\text {cp}} ; \gamma ^\text{r}_\text {cp})\right] , \end{aligned}$$
(21)
where \(\beta _{\mu } = \mu _{\infty } / \mu _{\text{s}} \ge 1\) denotes the stiffness ratio between cortex and subcortex, and \(\gamma ^\text{r}_{\text {cp}}\) represents the Heaviside spatial exponent at the cortex. The distribution of the shear modulus along the radial direction of the geometry is illustrated in Fig. 4 (top). In the same manner, to satisfy the differential growth discussed in Sect. “Mechanical growth problem”, the growth factors in the circumferential and radial directions \(\kappa ^\bot\) and \(\kappa ^\parallel\), respectively, are modeled as,$$\begin{aligned}&\kappa ^\bot (r_{\text{i}}) = \kappa _{\text{s}}+\left[ \kappa _{\text{s}} \left[ \beta _{\kappa }-1 \right] \mathcal {H} (r_{i} – r_{\text {cp}};\gamma ^\text{r}_\text {cp} ) \right] \quad \text {and} \nonumber \\&\kappa ^\parallel (r_{\text{i}}) = \kappa _{\text{s}}+\left[ \kappa _\text{s} \left[ \frac{1}{\beta _{\kappa }}-1 \right] \mathcal {H} (r_{\text{i}} -r_{\text {cp}};\gamma ^\text{r}_\text {cp} ) \right] , \end{aligned}$$
(22)
with the growth ratio \(\beta _{\kappa }\) given as,$$\begin{aligned} \beta _{\kappa }= \frac{ \kappa ^\bot }{ \kappa _{\text{s}}} = \frac{ \kappa _{\text{s}}}{\kappa ^\parallel } \ge 1. \end{aligned}$$
(23)
Again, Fig. 4 shows how the growth factor differentiates in the cortex into circumferential and radial factors. Finally, we introduce the parameters for the advection-diffusion equation in Sect. “Cell-density problem” to capture the distinct behaviors of each cell type. IPCs are characterized by a notably slow movement to occupy the space between the ventricular zone and the outer subventricular zone,$$\begin{aligned} v_{\, \text {\tiny IP}}(r_{\text{i}}) = v_{\, \text {\tiny IP}} \left[ 1- \mathcal {H} (r_{\text{i}} – r^\text{t}_{\text {osvz}}; \gamma ^\text{r}) \right] , \end{aligned}$$
(24)
while ORGCs are translocated directly after they are born to the outer subventricular zone,$$\begin{aligned} v_{\, \text {\tiny ORG}}(r_{\text{i}}) = v_{\, \text {\tiny ORG}} \left[ 1- \mathcal {H} (r_{\text{i}} – r^\text{t}_{\text {osvz}}; \gamma ^\text{r}) \right] . \end{aligned}$$
(25)
On the other hand, neurons migrate towards the cortex,$$\begin{aligned} v_{\, \text {\tiny N}}(r_{\text{i}}) = v_{\, \text {\tiny N}} \left[ 1- \mathcal {H} (r_{\text{i}} – r_{\text {cp}}; \gamma _{\text {cp}}^\text{r}) \right] . \end{aligned}$$
(26)
The second parameter that needs to be adjusted is the diffusivity. RGCs undergo diffusion after division within the ventricular zone,$$\begin{aligned} d_{\, \text {\tiny RG}}^{\, \mathrm c}(r_{\text{i}})&= d^{\, \mathrm c} \, \left[ 1- \mathcal {H} (r_{\text{i}} – r_{\text {vz}}; \gamma ^\text{r}) \right] . \end{aligned}$$
(27)
Similarly, IPCs and ORGCs disperse within their respective effective zones, namely between the ventricular zone and the outer subventricular zone,$$\begin{aligned} d_{\, \text {\tiny ORG}}^{\, \mathrm c}(r_{\text{i}}) = d_{\, \text {\tiny IP}}^{\,\mathrm c}(r_{\text{i}})&= d^{\, \mathrm c} \, \left[ 1- \mathcal {H} (r_{\text{i}} – r^\text{t}_{\text {osvz}}; \gamma ^\text{r}) \right] . \end{aligned}$$
(28)
In contrast, neurons exhibit a pure migration behavior within the subcortex. Once they reach the cortex, they diffuse to occupy the cortical layer,$$\begin{aligned} d_{\, \text {\tiny N}}^{\, \mathrm c}(r_{\text{i}})&= d_{\, \text {\tiny N}}^{\, \mathrm c} \, \mathcal {H} (r_{\text{i}} – r_{\text {cp}};\gamma _{\text {cp}}^\text{r} ). \end{aligned}$$
(29)
The radial distributions of the type-specific amount of transport and diffusivity are demonstrated in Fig. 4 (bottom).Fig. 4Distribution of model parameters along the domain’s radial direction from the ventricular zone to the outer cortical surface: shear modulus \(\mu\); tangential and radial growth factors \(\kappa ^\bot \, ,\, \kappa ^\parallel\); diffusivity \(d^\text{c}\) of neurons, radial glial cells (RGCs), and outer radial glial cells (ORGCs), intermediate progenitor cells (IPCs); migration speed v of neurons and translocation speed of outer radial glial cells and intermediate progenitor cells. The corresponding model parameters are summarized in Table 2.To account for the impact of the mitotic small translocation behavior of ORGCs, which significantly pushes the outer boundary of the subventricular zone between gestational weeks 11 and 2423, we define the time-dependent radius \(r^\text{t}_{\text {osvz}}\) in Eqs. 24 , 25, and 26. This radius is formulated as a linear function of time within the specified period, increasing from the inner subventricular zone upper boundary until reaching its maximum value, such that$$\begin{aligned} r^\text{t}_{\text {osvz}} = \text {min} \{ r_{\text {osvz}} \, , \, r_{\text {isvz}} + m_{\text {mst}} \, \langle t – t_{\text {P2}}\rangle \}, \end{aligned}$$
(30)
where \(r_{\text {osvz}}\) is the upper and final boundary of the outer subventricular zone, as mentioned in Sect. “Geometry and discretization in space”, \(m_{\text {mst}}\) is the mitotic small translocation constant and \(t_{\text {P2}}\) is the end time of the second phase. The constant \(m_{\text {mst}}\) is determined based on the observed distance between the inner and outer subventricular zones in histologically stained sections of the human fetal brain from gestational weeks 17–32. This constant is calculated as the ratio of this observed distance over the corresponding time period.In the equations defining the radial distribution of all previous parameters, \(\gamma ^\text{r}\) serves as the general spatial exponent for the Heaviside function, while \(\gamma _{\text {cp}}^\text{r}\) specifically applies at the subcortex-cortex interface. Both exponents are critical for controlling how smoothly parameters change at the boundaries between different zones. To calibrate these exponents, we analyze how well our simulation results match with histologically stained images of the human fetal brain around gestational week 18. Regarding the method used in preparing, staining, and processing these sections, please see14. Figure 5 illustrates a comparison between the observed cell density distribution along a trajectory from the ventricular surface to the outer brain surface for the stained section and simulation results for five different exponent combinations, as listed in Table 1. To enhance comparability and avoid differences in the dimensions between the histologically stained sections and the simulation domain, we normalize the domain’s radius according to the extension from the ventricular to the outer cortical surface in the stained sections. Furthermore, cell density is normalized relative to its average value. The distribution characterized by \(\gamma ^\text{r}=20\) and \(\gamma _{\text {cp}}^\text{r} = 50\) closely replicates the trend observed in the histologically stained sections. The distinct peaks at the ventricular zone, outer subventricular zone, and cortex align precisely with those identified in the stained images. Similarly, the simulation results for the parameter combination 5 and the corresponding distribution accurately depict the observed local minima at the inner subventricular zone and intermediate zone.Fig. 5Evolution of the cell density distribution along the normalized radius from the ventricular surface to the outer brain surface for numerical simulation results corresponding to gestational week 18 across different radial distributions of parameters, as introduced in Table 1, alongside a stained human brain section at gestational week 18 (HBS). The cell density is normalized to its average value. Combinations 1 to 5 display the simulation results, each employing a distinct set of exponents in the parameters’ radial distribution equations. The corresponding model parameters are summarized in Table 2.Table 1 Exponents of the parameters’ radial distribution equations corresponding to the distributions shown in Fig. 5.All model parameters are summarized in Table 2. In general, the mechanical parameters are adopted from our prior studies on human brain tissue2,9,34. Geometric parameters, in contrast, are evaluated from histologically stained sections of the human fetal brain, which reveal the boundaries of the germinal layers. However, an initial radius estimate is made to ensure optimal visualization of the results. In the results Sect. “Results and discussion”, we calibrate our simulation results for convenient comparison with Magnetic resonance images. Growth parameters were originally obtained from12, yet they have been adapted to guarantee algorithmic convergence. Finally, the parameters for the cell density problem have been carefully adjusted experimentally to better align with actual human brain development, while preserving numerical stability.Table 2 Model parameters.To complete the set of model parameters, we need to define the division ratio between cell types for each division phase, i.e., \(G_{ \, \cdot } (\circ , \text {P})\). This step required an initial understanding of the lineage relationships between different cell types. Figure 6 illustrates the lineage tree and progression among the four considered cell types during the various gestational weeks. According to the literature3,5,24,43,44,45,46, the lineage tree can be organized into five primary phases, designated from P1 through P5, advancing from the early to the later stages of development. Additionally, defining a timeline for our simulation is essential; for this purpose, we establish a formula linking gestational weeks \(t_{\text {GW}}\) to computational time t, expressed as$$\begin{aligned} t_{\text {GW}} = 0.3 \times t + 4. \end{aligned}$$
(31)
The division ratios between cells across different division phases are detailed in Supplementary Table 1. Here, we provide an example of how we determined the division ratios. Considering the division of RGCs during the third phase (P3), each RGC divides into one RGC and one ORGC, represented as \(G_{\text {\tiny RG}} (\text { RG}, \text { P3}) = 1\) and \(G_{\text {\tiny ORG}} (\text { RG}, \text { P3}) = 1\), respectively. ORGCs then divide again, each giving rise to another ORGC, thus \(G_{\text {\tiny ORG}} (\text { ORG}, \text { P3}) = 1\). Meanwhile, each IPC divides to form two new IPCs, and one of these IPCs further divides to produce two neurons, leading to \(G_{\text {\tiny IP}} (\text { IP}, \text { P3}) = 1\) and \(G_{\text {\tiny N}} (\text { IP}, \text { P3}) = 2\). Following the calculation of division ratios, the source term can be computed using Eq. (14), as outlined in the Supplementary Material.Fig. 6The lineage tree and progression among the four fundamental neural progenitor cells during various gestational weeks in human brain development. This lineage tree is organized into five primary phases, designated P1 through P5, advancing from the early to the later stages of development. The term “virtual case of the intermediate progenitor cells” describes a scenario that physically transpires within the brain during development, yet is not accounted for in numerical simulations.
Exploring the role of different cell types on cortical folding in the developing human brain through computational modeling
