Orchestrated neuronal migration and cortical folding: A computational and experimental study

Brain development involves precisely orchestrated genetic, biochemical, and mechanical events. At the cellular level, neuronal proliferation in the innermost zone of the brain followed by migration towards the outermost layer results in a rapid increase in brain surface area, outpacing the volumetric growth of the brain, and forming the highly folded cortex. This work aims to provide mechanistic insights into the process of brain development and cortical folding using a biomechanical model that couples cell division and migration with volumetric growth. Unlike phenomenological growth models, our model tracks the spatio-temporal development of cohorts of neurons born at different times, with each cohort modeled separately as an advection-diffusion process and the total cell density determining the extent of volume growth. We numerically implement our model in Abaqus/Standard (2020) by writing user-defined element (UEL) subroutines. For model calibration, we apply in utero electroporation (IUE) to ferret brains to visualize and track cohorts of neurons born at different stages of embryonic development. Our calibrated simulations of cortical folding align qualitatively with the ferret experiments. We have made our experimental data and finite-element implementation available online to offer other researchers a modeling platform for future study of neurological disorders associated with atypical neurodevelopment and cortical malformations.


Finite element implementation
First, we briefly summarize the strong form of the governing equations in Section 4.4. The strong form of the coupled PDEs are given by u =ȗ on S u , Tn =t on S t .
Balance of momentum (S.1) Next, we let w 1 , w 2 denote two test fields which vanish on S c and S u , respectively. After multiplying the PDEs with two test fields and performing integration by parts, the weak form could be expressed as As is routine, the body is discretized into finite elements such that B t = ∪B e t , and the nodal variables are taken to be the cell densities and displacement, which are interpolated inside each element by with the index A = 1, 2, ... denoting the nodes of the element, c A i and u A are the nodal cell densities and displacements, respectively, and N A the shape functions. Next, the test fields w 1 and w 2 are interpolated by the same shape functions in the Galerkin approach, viz.
After substitution of Eq (S.3) and Eq (S.4) into the weak forms in Eq (S.2), we arrive at elementlevel of residuals. Since the expressions are similar for different neuron types, we only demonstrate the residual for c 1 (x, t), The element level residual for the displacement field u(x, t) is given by Since we consider cohorts of electroporated neurons at three different times in this work, four element-level residuals in total are assembled into a global residual, which, when set to zero, represents a system of non-linear equations for the nodal degrees of freedom.
The following tangents are required by the iterative Newton-Raphson procedure for convergence. More specifically, the tangent K AB uu is given by where the spatial tangent modulus A is related to referential tangent modulus A R through and the referential tangent modulus is given by Piola stress [1]. After using Eq (15), the Piola stress is now given by The three diagonal tangents for neuron density are similar, and here we only demonstrate K AB c 1 c 1 , which is given by where the derivative ∂H ∂c 1 is given by We note the current work does not address coupling between cohorts of neurons, thus we set The displacement-cell density tangents are similar among the three types of neurons; we only demonstrate K AB u i c 1 , which is given by where the derivative ∂T ∂c 1 ij could be further expressed as (S.14) The derivative ∂T ∂F g ijkl in Eq (S.14) is given by The derivative ∂F g ∂c 1 kl in Eq (S.14) is given by The cell density-displacement tangents among the three neuron types are similar, and we only demonstrate K AB c 1 u k , which is given by After implementing a two-noded linear element in Matlab for verification of our numerical implementation (see below), we then implement four-noded quadrilateral and eight-noded brick elements, (UPE4 and U3D8, respectively) in [2] by writing user-defined element (UEL) subroutines, following our previous work [1,3,4,5,6].

Verification of our numerical implementation
In order to verify our finite element implementations, we decouple the problem and consider each phenomenon (spatiotemporal cell density and mechanical deformation) separately. For the former, we compare our independent numerical solutions in Matlab and Abaqus to each other; in the latter case, we compare our numerical solutions in Abaqus with an analytically tractable solution.
First, we consider cell migration in a rigid slender body with a length of l = 1 m. Here, we only restrict our attention on one cohort of cells, c 1 , and allow them to migrate along x direction. Under these assumptions, the PDE in Eq (5) is now reduced to 1-D, The cell density flux q 1 in Eq (S.18) is given by where the velocity fieldv 1 (x) is expressed aŝ The source term f c 1 in Eq (S.18) is given by the following function,  Fig 1a). For the mechanical problem, we prescribe a simple shear motion with an angle of θ to a unit cube ( SI Fig 1b), such that the associated deformation gradient is given by [7] where γ = tan θ denotes the amount of shear. Next, two assumptions are made make the analytical solution tractable: 1) perfect incompressibility (i.e. J e = 1), and 2) no growth (i.e., F g = 1). Under these assumptions, the Cauchy stress in Eq (15) is now given by where P is the constitutively indeterminate pressure that is introduced to satisfy the incompressibility constraint. In our simulation, we prescribe the same deformation to a single U3D8 element in Abaqus. Note that in our numerical treatment, we take L = 10 3 µ to approximate a nearly incompressible material in the simulation. We compare the analytical solutions for shear stress and normal stress difference, given by T 12 = µγ and T 11 − T 33 = µγ 2 , (S.24) respectively, against the numerical solutions (SI Fig 1c). Once again, the excellent agreement between analytical and numerical results indicate the mechanical portion of our finite element implementation is verified.

Robustness of genetic algorithm
To investigate the robustness of the genetic algorithm used in our model calibration, we focus on exploring the evolution of some of the key material parameters. More specifically, we keep track of four key material parameters (G c , v i , D, and k s ) and plot them as a function of genome up to 40 generations (SI Fig 2). We note that, under different ranges of initial guesses used in the calibration, most parameters can converge to steady values after roughly about 10 generations/100 genomes.
Parameter G c appears to fall into a local minimum around 100 genomes, before converging to a lower value after approximately 300 genomes. This study reveals the robustness of the genetic algorithm.