|Home | About | Journals | Submit | Contact Us | Français|
Folding of the human cerebral cortex has intrigued many people for many years. Quantitative description of cortical folding pattern and understanding of the underlying mechanisms have emerged as an important research goal. This paper presents a computational 3D geometric model of cerebral cortex folding that is initialized by MRI data of human fetus brain and deformed under the governance of partial differential equations modeling the cortical growth. The simulations of this 3D geometric model provide computational experiment support to the following hypotheses: 1) Mechanical constraints of the brain skull influence the cortical folding process. 2) The cortical folding pattern is dependent on the global cell growth rate in the whole cortex. 3) The cortical folding pattern is dependent on relative degrees of tethering of different cortical areas and the initial geometry.
Anatomy of the human cerebral cortex is extremely variable across individuals in terms of its size, shape and structure patterning . The fact that folding pattern of the cortex is a good predictor of its function  has intrigued many people for many years [3~5]. Recently, understanding of the underlying folding mechanisms [3~5] have emerged as important research goals. Though each human brain grows from a similar shape of neuronal tube in the very beginning, the major cortical folding variation develops from 8 months of fetus brain. Since many developmental processes are involved in this cortex folding development, including neuronal proliferation, migration and differentiation, glial cell proliferation, programmed cell death, axon development and synaptogenesis, how these processes interact with each other and dynamically accomplish the cortical folding is still largely unknown .
In the neuroscience community, several hypotheses have been proposed to explain the gyrification or folding of the cerebral cortex from different perspectives [3, 6]. Since the cortex area becomes almost three times larger than the cranial area in human brain after evolution, the mechanical constraint was firstly considered as the major factor that determined the cortical folding pattern . Then, the areal difference, especially the cytoarchitectonic difference, which causes regional mechanical property variation, is considered as the direct reason that causes gyrification . The axongenesis process, which might cause areal differentiation [3, 6], is also considered as an important factor that determines the folding pattern.
Due to the complexity of the neurobiology involved in gyrification of the cerebral cortex, computational modeling of such a process is very challenging. In the literature, there have been a couple of attempts to develop computational models to understand the cortex folding process. For example, in , the authors proposed a 2D continuum mechanics based model of growth to synthesize brain cortex shapes by using physical laws. Recently, a computational morphogenetic model was proposed to study the fundamental mechanisms of cortical folding in . In this model, the 2D annular cortex and radial glial fibres were modeled by finite-elements and the growth of cortex was modeled as the expansion of finite elements.
In this paper, a 3D morphogenetic model is proposed to study the mechanisms of cortical gyrification. Since the human cortex can be considered as a highly convoluted thin shell , the triangulated surface of the developing cortex is adopted as its representation. The geometric model is driven by mechanical forces occurring in the growing brain for simulation of shape dynamics of the developing cortex. Parametric models for cerebral cortex are initialized by the structural fetus MRI data  from 22 to 36 week gestation. The morphogenetic model is deformed under the governance of partial differential equations. The mechanical forces and boundary conditions are guided by MRI data. After applying mechanical and growth properties of cortex on each triangular element, the time-varying system is solved by the Newmark scheme.
T2 MRI data of fetus brain (SE sequence TR=1329.88ms and TE=88.256ms)  from 22 to 36 gestation weeks are used for the cortical surface reconstruction and folding model initialization. Fetus brain is firstly extracted from the structural MRI data and then the skull, cortex plate, white matter zone (including subplate and intermediate zones), ventricle zone, basal ganglia and thalami are distinguished from the brain regions (Figure 1a). The inner surface of the skull and outer surface of cortex plate are reconstructed from the volumes (Figure 1b).
The proposed computational model of folding of the human cerebral cortex is composed of four key components: 1) Deformable model of the cortex. This model provides a geometric representation of the cortex that can be deformed into dynamic shapes via mechanical forces. 2) Driving forces of the folding. The mechanical driving forces are inferred from the cortex growth model and deform the cortical surface under partial differential equations. 3) Geometric constraints. The growth of the cortex is constrained by the geometric constraints of the brain structures, as well as boundary conditions. 4) Model solvers. The numerical solution to the partial differential equations provides the dynamic evolution of the geometric surface model of the cortex. Figure 2 provides an overview of the proposed computational model.
To model the mechanical properties of developing cortex, the elasto-plasticity model is adopted in our methods, akin to that in . The elasticity property enforces the surface model to restore to the original shape and at the same time, the plasticity property tends to maintain the deformed shape permanently.
The developing cortex is represented as tens of thousands of triangle elements (21,136 in our experiments), each of which presents a small cortical region. During the development or growth of cortex, the elastic and plastic properties could be obtained from the deformed ( ) and reference ( ) coordinates of the triangle corner i . The elastic force of each triangle element is defined along each edge of the triangle:
where i and j are the triangle vertex indices, is the rest length of edge ij , and is the current distance between vertices i and j , Kc is the elastic constant. As plasticity plays an important role in soft tissue deformation, in our model, the plasticity of cortex is modeled as an adaptation of the reference configuration to the deformed configuration:
where τce is the time constant for the plasticity, which is similar to that in .
Since the cortex is modeled as a zero thickness surface, the rigidity of cortex is introduced as bending energy . As computing mean curvature on the triangulated surface is computationally expensive, the simplified version of bending stress is defined on each edge of the surface triangle:
where i and j are the two vertices of edge e , p and q are the two vertices that share the same triangle with edge e as illustrated, is the rest distance between vertices p and q , is the current distance between vertices p and q , and Kb is the elastic constant. Also, the bending energy decreases with the plasticity of :
where τcb is the time constant for the plasticity of bending energy.
Since modeling and quantification of average size and number of neurons at the cellular level is still an open problem in the computational neuroscience field, we adopt the classic logistic-growth function  to describe the growth of cortical tissues as that in :
where m is known as the Malthusian parameter and k is the carrying capacity of the system . By changing the rest area of triangular element, the growth of cortex will generate mechanical stress that partly drives the deformation of the cerebral cortex surface.
The development of cortex is limited by some mechanical or boundary conditions such as cranial volume and self collision . To prevent the cortex from developing into the skull or other brain tissues, a volumetric constraint model is maintained during the simulated folding of the cortex. Voxels from the skull, basal ganglia/thalami or ventricular zone are extracted from the scanned 3D MRI image and grouped into a new image as a mask, outside which volumes of other brain regions are set to zero. And then, the cortex model can only be deformed in the zero value space.
Also, we applied similar techniques in  to prevent the self-intersection of deformed cortical surfaces. Specifically, current deformed cortex surface is rasterized to a volumetric model. When any vertex of the surface is being deformed to a new position x′ , the new position would be checked if it's inside the developing skull or other brain tissue volume, or cause its neighboring triangles to intersect with the other rasterized triangles. Otherwise another new position for this vertex would be found.
The proposed system can be formulated as a time-varying partial differential equation, which is commonly adopted in many deformable model approaches . Specifically, the dynamics of each cortical surface vertex can be simply formulated as i=fi/mi, where i is the acceleration of the vertex i , fi is the force that combines all forces affecting vertex i and mi is the mass of vertex i , which is usually defined as the sum of one thirds of the masses of all triangles around vertex i . By combining all equations on the cortical surface together, we have the discrete form of developing cerebral cortex as:
where x , and are the vertex's position, velocity and acceleration respectively, M is a 3n×3n ( n is the number of vertices) diagonal mass matrix on vertices where diag(M) = (m1, m1, m1, m2, m2, m2, ..., mn, mn, mn ) and F(x, ) is the net force vector that combines all forces on cortical surface vertices. The Newmark scheme, which is widely used for solving ODE (ordinary differential equation) integration problem, is adopted in our implementation.
We believe that the combination of computational simulation of cortex folding and invivo neuroimaging of fetus brain development will be very helpful to understand the mechanisms of cortical folding, as well as interactions between different mechanisms. Also, in the long term, the combination of computational modeling and neuroimaging might be helpful to understand the mechanisms of many neurological disorders that result in abnormal cortical folding.