PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Med Image Comput Comput Assist Interv. Author manuscript; available in PMC 2010 July 13.
Published in final edited form as:
Med Image Comput Comput Assist Interv. 2009; 12(Pt 2): 458–465.
PMCID: PMC2902882
NIHMSID: NIHMS135905

A Computational Model of Cerebral Cortex Folding

Abstract

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.

Keywords: cortex folding, simulation

1 Introduction

Anatomy of the human cerebral cortex is extremely variable across individuals in terms of its size, shape and structure patterning [1]. The fact that folding pattern of the cortex is a good predictor of its function [2] 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 [6].

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 [7]. Then, the areal difference, especially the cytoarchitectonic difference, which causes regional mechanical property variation, is considered as the direct reason that causes gyrification [7]. 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 [4], 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 [5]. 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 [2], 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 [8] 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.

2 Methods

2.1 Materials and Pre-processing

T2 MRI data of fetus brain (SE sequence TR=1329.88ms and TE=88.256ms) [8] 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).

Figure 1
Fetus brain structure reconstruction.

2.2 Computational Model

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.

Figure 2
Flowchart of the simulation model.

2.2.1 Deformable model

To model the mechanical properties of developing cortex, the elasto-plasticity model is adopted in our methods, akin to that in [5]. 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 ( xci ) and reference ( xc0i ) coordinates of the triangle corner i . The elastic force of each triangle element is defined along each edge of the triangle:

fcei,j=Kc(lc0i,jlci,j)lc0i,jlci,j(xcixcj),i,j{1,2,3},ij
(1)

where i and j are the triangle vertex indices, lc0i,j=xc0ixc0j is the rest length of edge ij , and lci,j 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:

dlc0i,jdt=1τce(lci,jlc0i,j),i,j{1,2,3},ij
(2)

where τce is the time constant for the plasticity, which is similar to that in [12].

Since the cortex is modeled as a zero thickness surface, the rigidity of cortex is introduced as bending energy [9]. 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:

fbei,j=Kb(lc0p,qlcp,q)lc0p,qlcp,qf(xcpxcq),i,jedgee
(3)

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, lc0p,q=xc0pxc0q is the rest distance between vertices p and q , lcp,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 lcp,q :

dlc0p,qdt=1τcb(lcp,qlc0p,q)
(4)

where τcb is the time constant for the plasticity of bending energy.

2.2.2 Cortex Growth Model

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 [10] to describe the growth of cortical tissues as that in [5]:

dAc0dt=Ac0m(1Ac0k)
(5)

where m is known as the Malthusian parameter and k is the carrying capacity of the system [10]. 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.

2.2.3 Constraints

The development of cortex is limited by some mechanical or boundary conditions such as cranial volume and self collision [7]. 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 [11] 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 x for this vertex would be found.

2.2.4 Model Solver

The proposed system can be formulated as a time-varying partial differential equation, which is commonly adopted in many deformable model approaches [9]. Specifically, the dynamics of each cortical surface vertex can be simply formulated as [x with umlaut]i=fi/mi, where [x with umlaut]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:

x¨=M1F(x,x.)
(6)

where x , x and [x with umlaut] 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, 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.

3. Results

  1. We simulated the computational models of cortical folding introduced in Section 2. The parameters in the deformable model and growth model were same in both simulations with and without skull constraints, which were set as: Kc = 1 , τce = 1000 , Kb = 5 , τcb = 1000 , m = 0.002 , k = 3 , Δt = 0.05 . The folding development results of the 3D morphogenetic model are illustrated in Figure 3. The Figure 3a and Figure 3b show the snapshots of cortex development at the iteration number 0, 40, 80, 120, 160 and 200 with and without skull constraints respectively. By visual evaluation, it is evident that reasonably realistic convolutions are generated during these simulations. To quantitatively evaluate the produced convolutions, we use the curvature as the measurement of cortical folding. Also, we use the total surface area of the developing cortex to measure its growth rate. The differences between the total area of cortex and the average absolute Gaussian curvature of the cortex surface are illustrated in Figure 4a and Figure 4b respectively. Figure 4a shows that the cortical surface area increases almost linearly with the simulation iteration numbers, irrespective of whether or not there is skull constraint. The results in Figure 4a indicate that the cortex increases its surface area and convolutes itself to reduce the fast growing internal tension with or without skull constraint. However, as shown in Figure 4b, it is intriguing that the increase in speed of average curvature in simulation with skull constraint is much higher than the one without skull constraint, meaning that cortex development with skull constraint is much more convoluted than that without skull constraint. This computational simulation result provides direct evidence to the following hypothesis: mechanical constraints of the brain skull influence the cortical folding process. The results in Figure 3 and Figure 4 also tell that the cortex growth pattern is quite different because of the skull constraint. Without skull constraint, the cortex expands with the similar shape of the initial cortex during the first 40 iterations and then small shape changes appear during 40~80 iterations. After that, primary sulci rapidly develope during 80~160 iterations and then, secondary sulci develop between major sulci. In this model simulation, smaller tertiary gyri and sulci are not found to be developed. Since most of the primary sulci developed in this model simulation are long but not deep, the average absolute Gaussian curvature is slowly increased during the cortical development process. It is striking that the cortex growth pattern in the simulation with skull constraint is quite different, that is, the primary and secondary sulci are developed much earlier and faster, e.g., in the iteration 40~120. Also, smaller tertiary gyri and sulci are developed during last iterations (120~200). The sulci developed in this model simulation are deeper, compared to those simulations without skull constraint. Thus, the average absolute Gaussian curvature increases faster as shown in Figure 4b. The quite different cortical growth patterns, due to the existence of skull constraints or not further demonstrate that mechanical constraints of the brain skull significantly influence the cortical folding process.
    Figure 3
    Cortical development with or without skull constraints. (a) and (b) show the snapshots of cortex development at the iteration number 0, 40, 80, 120, 160 and 200 with or without skull constraints respectively. The parameters are set as: Kc = 1 , τ ...
    Figure 4
    The differences of the total area of cortex and the average absolute Gaussian curvature of the cortex surface during convolution development.
  2. This experiment investigates the effect of global cell growth rate on the cortical folding process. The Malthusian parameter m in Eq.(5) models the cell growth rate. When the growth parameter m is increased from 0.002 (Figure 5a) to 0.004 (Figure 5b), the cortical growth speed is almost doubled according Eq. (5). As a result, more deeper sulci and smaller gyri can be found in the simulated folding, as shown in Figure 5b. Also, the average absolute mean curvature is significantly increased from 0.31 to 0.35 in the cortex development simulation with the skull constraint. These simulation results provide computational experimental support to the following hypothesis: the cortical folding pattern is dependent on the global cell growth rate in the whole cortex. It is interesting that the effect of global cell growth rate on the cortical folding process is much less significant in the simulations without skull constraint, as shown in the Figure 5c and 5d. For example, even though the parameter m is increased to 0.004 from 0.002, the average absolute mean curvature is only slightly increased. This result further demonstrates that the skull constraint is an important modulator of cortical folding pattern formation. The computational results in this experiment demonstrate that overgrowth of the cortex will increase the cortical folding and convolution. This simulation result might provide theoretical clues to the following two independent studies: 1) In a MRI study of Autism, it was reported that the left frontal cortical folding is significantly increased in autism patient, but there is also a significant decrease in the frontal folding patterns with age in the autistic group [12]; and 2) It was reported that in Autism, the brain overgrows at the beginning of life and slows or arrests growth during early childhood [13]. In the future, it would be very interesting to combine our folding simulation studies with longitudinal MRI studies of Autistic brains and normal controls to further elucidate the relationship between brain growth and cortical folding in both Autism and normal neurodevelopment.
    Figure 5
    (a) and (b) are simulations with skull constraints. (a) The parameters are set as those in Figure 3. (b) The parameter m is changed to 0.004. (c) and (c) are simulations without skull constraints. Parameters are the same as (a) and (b).
  3. It was hypothezied that relative degrees of tethering of different cortical areas or initial geometry of cortex had significant influence on cortical folding pattern [7]. To examine this hypothesis, a small artificial deformation is applied to the originally reconstructed cortex surface as illustrated in Figure 6c. The green region is manually labeled and is deformed around 2mm towards the inside ventricular zone. The cortical folding simulation results in Figure 6b and 6d demonstrate that initial geometry of the cortex has significant influence on the resulted cortical folding pattern. However, if there is no skull constraint in the simulation, less influence of the intial geometry on cortical folding can be seen, compared to the result with original shape of cortex (figure not shown). It is also apparent from the results in Figure 6 that the local shape differences not only cause the folding pattern variation near the locally deformed region, but also influence the cortical folding pattern all over the cortex, as denoted by red arrows in Figure 6. This result demonstrates that cortical folding process is a global and system-level behaviour, but not a localized procedure, and provides computational experiment support to the idea that cortical folding process might be a minimization of global energy function [3].
    Figure 6
    Effect of initial cortex shape on the development and distribution of convolution. Compared to the original cortex (a), the green region in (c) is deformed around 2mm inside. (b) and (d) are the results of 200 iterations of simulation. Other parameters ...

4. Conclusion

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.

References

1. Talairach J, Tournoux P. Co-planar Stereotaxic Atlas of the Human Brain. Thieme; New York: 1988.
2. Fischl B, et al. Cortical Folding Patterns and Predicting Cytoarchitecture. Cereb Cortex. 2008;18(8):1973–1980. [PubMed]
3. Van Essen D. A tension-based theory of morphogenesis and compact wiring in the central nervous system. Nature. 1997;385:313–318. [PubMed]
4. Raghavan R, Lawton W, Ranjan SR, Viswanathan RR. A continuum mechanics-based model for cortical Growth. J. Theor. Biol. 1997;187:285–296.
5. Toro R, Burnod Y. A Morphogenetic Model of the Development of Cortical Convolutions. Cerebral Cortex. 2005;15:1900–1913. [PubMed]
6. Rakic P. A Century of Progress in Corticoneurogenesis: From Silver Impregnation to Genetic Engineering. Cerebral Cortex. 2006;16:i3–i17. [PubMed]
7. Brown M, et al. The developing brain. Oxford University Press; Oxford: 2002.
8. Beth Israel Deaconess Medical Center http://radnet.bidmc.harvard.edu/fetalatlas/
9. Terzopoulos D, Platt J, Barr A, Fleischer K. Elastically Deformable Models. Computer Graphics. 1987;21:205–214.
10. Murray J. Mathematical biology, Heidelberg. Springer-Verlag; Germany: 1993.
11. Fischl B, Sereno MI, Dale AM. Cortical surface-based analysis II: inflation, flattening and a surface-based coordinate system. NeuroImage. 1999;9:195–207. [PubMed]
12. Antonio Y, et al. Increased frontal cortical folding in autism: a preliminary MRI study, Psychiatry Research: Neuroimaging. 2004;131:263–268. [PubMed]
13. Courchesne E, et al. Mapping early brain development in autism. Neuron. 2007 Oct 25;56(2):399–413. [PubMed]