Human cerebral cortex develops extremely fast in the first year of life. Quantitative measurement of cortical development during this early stage plays an important role in revealing the relationship between cortical structural and high-level functional development. This paper presents a computational growth model to simulate the dynamic development of the cerebral cortex from birth to 1 year old by modeling the cerebral cortex as a deformable elastoplasticity surface driven via a growth model. To achieve a high accuracy, a guidance model is also incorporated to estimate the growth parameters and cortical shapes at later developmental stages. The proposed growth model has been applied to 10 healthy subjects with longitudinal brain MR images acquired at every 3 months from birth to 1 year old. The experimental results show that our proposed method can capture the dynamic developmental process of the cortex, with the average surface distance error smaller than 0.6 mm compared with the ground truth surfaces, and the results also show that 1) the curvedness and sharpness decrease from 2 weeks to 12 months and 2) the frontal lobe shows rapidly increasing cortical folding during this period, with relatively slower increase of the cortical folding in the occipital and parietal lobes.
The human cerebral cortex develops from the smooth neural tube, but its folding pattern is extremely complex and variable across different individuals (Talairach and Tournoux 1988). Although the primary and secondary folding structures have been developed before birth, the tertiary folding structures are still undergoing rapid development (Toro and Burnod 2005), together with the high-level function, after birth. On the other hand, with the acquisition of more and more high-resolution brain MR images from infants (Shaw et al. 2007; Dubois et al. 2008; Kaukola et al. 2009; Awate et al. 2010; Schumann et al. 2010; Shen et al. 2010), it becomes possible to partially reveal the relationship between morphological and functional development of the cerebral cortex in the first year of life (Dubois et al. 2008; Awate et al. 2010) and to characterize certain brain disorders that occur during this stage (Shaw et al. 2007; Kaukola et al. 2009; Schumann et al. 2010) by analysis of longitudinal cortical development.
For early brain development study, all images acquired at different developmental stages need to be warped onto a common space. This is generally achieved by image- or surface-based registration methods (Woods et al. 1998; Shen and Davatzikos 2002; Liu et al. 2004; Vercauteren et al. 2009). However, these registration methods usually use intensity- or feature-based similarity to locate correspondences, which may lead to incorrect correspondence detection especially in the flat cortical surface regions and thus eventually affecting the accurate measurement of cortical development.
Accordingly, several computational mechanical models have been proposed to study the mechanisms of the cortical folding (Toro and Burnod 2005; Geng et al. 2009; Nie et al. 2010). Specifically, in Toro and Burnod (2005), the dynamics of the cortical folding was modeled by a 2D computational mechanical model to study the fundamental mechanisms of the cortical folding. In a recent study (Geng et al. 2009), a finite element method was proposed to model the cortical folding (regulated by tension forces) in which diffusion tensor imaging (DTI) data of the fetal sheep brain was adopted as a cue to model the tension forces that regulate the cortical folding. Recently, the dynamics of the cortical folding was modeled by a deformable computational mechanical model (Nie et al. 2010), in which the cortex is modeled as a deformable surface deforming in the 3D space. The surface-based deformable mechanical model, which we use for cortical development simulation in this paper, is widely adopted for dynamic mechanical simulation, such as cloth and paper (Grinspun et al. 2003), and thus suitable for cortical growth since the cortex (with the thickness between 1 and 5 mm) can also be considered as a thin sheet (Van Essen and Maunsell 1980; Fischl et al. 1999; Xu et al. 1999). By appropriate estimation of mechanical parameters, the real dynamics of the cortical development could be simulated. In particular, during the simulation, each element of the cortex can be easily tracked, and its deformation can be captured by the mechanical model. Thus, these mechanical models can provide not only the dynamics of the cortical development but also more accurate temporal correspondences for the cortex across different brain developmental stages.
In this paper, we present a novel computational growth model for capturing human cortical development in the first year of life. Specifically, a mechanical model is adopted to describe the mechanical properties of the cortex, and also the classic logistic growth function (Murray 1993) is adopted to describe the growth of the cortex. To achieve an accurate simulation, a guidance model is further incorporated into the mechanical model to estimate parameters and guide the simulation. Our proposed method has been tested on 10 healthy subjects with longitudinal brain MR images acquired at every 3 months from birth to 1 year old. The experimental results show that the simulated cortical surfaces are highly similar to the individual’s cortical surfaces reconstructed from MR images at different developmental stages (in terms of both shapes and positions), and also our proposed method can provide a better way of measuring cortical correspondences at any time point in the first year of life.
Materials and Methods
Data Set and Preprocessing
Longitudinal T1, T2, and diffusion-weighted MR images of 10 healthy subjects were acquired at every 3 months from birth to 1 year old, using a 3-T head-only MR scanner. The imaging parameters for T1 images (with 160 axial slices) are as follows: time repetition [TR] = 1900 ms, time echo [TE] = 4.38 ms, flip angle = 7, acquisition matrix = 256 × 192, and voxel resolution = 1 × 1 × 1 mm3; the imaging parameters for T2 images (with 70 axial slices) are as follows: TR = 7380 ms, TE = 119 ms, flip angle = 150, acquisition matrix = 256 × 128, and voxel resolution = 1.25 × 1.25 × 1.95 mm3; and the imaging parameters for diffusion-weighted images (with 60 axial slices) are as follows: TR/TE = 7680/82 ms, matrix size = 128 × 96, 42 noncollinear diffusion gradients, and diffusion weighting b = 1000 s/mm2. T2 image is linearly aligned onto the respective T1 image of the same subject and further resampled to 1 × 1 × 1 mm3. Fractional anisotropy (FA) image is reconstructed from the DTI image and then aligned onto the T1 image of the same subject and further resampled to 1 × 1 × 1 mm3. For each set of the aligned T1, T2, and FA images, the skull stripping is first performed to remove noncerebral tissues, and also the cerebellum and brain stem are removed by the in-house–developed tools (as shown at the first, second, and third columns of Fig. 1). Since the delineation of the gray–white interface at 6 months of age is difficult due to low tissue contrast in the T1 and T2 MR images (as shown in Fig. 1), the DTI images, especially the FA images, are further adopted (as in Liu et al. (2007)) to improve the segmentation of white matter (WM), especially for 6-month images. Specifically, the combined T1, T2, and FA image information is used to segment the brain image into gray matter (GM), WM, and cerebrospinal fluid regions (Wang et al. 2011), as shown at the fourth column of Figure 1. The results in Wang et al. (2011), which compare the automatic segmentation results with manual segmentation results, show that the average Dice overlay ratio of cortex is about 0.9 at 6 months, which is similar to the average Dice overlay ratio of 0.9 at 3, 9, 12 months.
After topology correction of WM volume (Shattuck and Leahy 2001), the inner cortical surface (i.e., the interface between WM and GM) at different developmental stages can be reconstructed and represented by the triangular meshes composed of a set of vertices and triangles (as shown at the fifth column of Fig. 1). Since the transient subplate zone, which is interposed between the immature cortical plate and WM, may still exist at 2 weeks after birth, the inner cortical surface at 2 weeks is defined as the interface between the cortex plate and WM zone (including WM and transient subplate zone). Due to the effect of partial volume in MR images, the outer surface is relatively hard to be accurately reconstructed, especially in the deep sulcal regions (Fischl et al. 1999; Xu et al. 1999). Thus, the inner cortical surface at birth is used as the initial surface in our model to simulate the cortical development in the first year of life, as detailed below.
Computational Growth Model
In our computational growth model, the cerebral cortex is represented by its inner cortical surface and modeled by a deformable sheet with elastoplasticity property. In this way, the cortical surface can be deformed into dynamic shapes via mechanical forces and can also be restored to its original shape (elasticity) or maintain the dynamic shape permanently (plasticity). Specifically, the development of the cortical surface is modeled by a growth model, which generates the driving force by increasing the “rest” area of the cortical surface (here “rest” indicates the steady status of the cortical surface under no external force). To reduce the simulation error, a guidance model is then introduced to estimate both the growth parameters and the cortical shapes at later developmental stages. Also, a volumetric constraint model is adopted to prevent the cortical surface from deforming into other brain tissues and self-collision. These models are finally incorporated into a time-varying system and further solved by the Newmark scheme (Newmark 1959).
In the following, we will use , , , , and to denote the inner cortical surfaces reconstructed from the longitudinal images acquired at 2 weeks, 3, 6, 9, and 12 months old of each subject, respectively. We also use to represent them. Note that our goal is to develop a growth model that can simulate the continuous cortical development under certain mechanical constraints and also approach all those given cortical surfaces at the known time points.
Deformable Mechanical Model
The elastoplasticity model is adopted on each edge of the triangulated inner cortical surface to model the mechanical properties of the developing cortex. For each edge, both “elastic property” and “bending elastic property” are introduced to restore it to its rest shape, and also both plastic property and bending plastic property are introduced to maintain the current shape. An illustration of elastoplasticity property is shown in Figure 2. For each edge element (i, j) and its 2 neighboring triangles Δpij and Δqij, after contracting and bending to a new shape (Fig. 2b) from the original shape (Fig. 2a), the elastic forces are generated on both the edge (i, j) and the virtual edge (p, q) to restore it to its original shape (Fig. 2d), and the rest lengths of the edge (i, j) and the virtual edge (p, q) are adapted to the current lengths by the plastic property (Fig. 2c).
At the same time, the plastic property on each edge (i, j) (as shown in Fig. 2c) for maintaining the developed length of the edge is modeled as the adaptation of the rest length to the deformed length :
Bending elastic property and bending plastic property. To model the rigidity of the cortex, the “bending” elastoplasticity property is also defined on each edge (i, j), through a virtual edge (p, ) as shown in Figure 2c, where p and q are the vertices that share the same 2 triangles with the edge (i, j). When the angle between the triangles Δpij and Δqij changes, as shown in Figure 2b, the distance between vertices p and q will change accordingly. Thus, similar to equation (1), the bending elastic property fb on each virtual edge (p, q) is modeled as:
Cortex Growth Model
Although most neuronal proliferation and migration have been completed before birth (Brown et al. 2002), the volume of the cortex and the size of the brain still keep developing until adolescence because of the further postnatal growth of neuronal dendrites, myelination, and gliogenesis and angiogenesis. Given the later developed surface area of the cortex , the cortex growth model at time t can be described as the classic logistic growth function (Murray 1993; Toro and Burnod 2005):
Since we have reconstructed cortical surfaces from the later developed stages , these cortical surfaces can be used to guide the cortical development in 2 aspects. First, the parameter in equation (5) can be estimated based on the area of the surface . Second, the simulation to the realistic cortical shape and position can be guided by the shape and position of the cortical surface .
The initial cortical surface is first warped into the space of the target cortical surface (i.e., the inner cortical surface at 3, 6, 9, and 12 months) by a high-dimensional nonlinear hybrid (volumetric/surface) registration method (Shen and Davatzikos 2002; Liu et al. 2004). Thus, for each vertex i with the position on the initial cortical surface , a rough corresponding vertex can be found on the target cortical surface . Then, the guidance force can be modeled as:
Since the correspondences in the flat cortical region with low curvature value, such as sulcal wall, determined by the hybrid registration method (Liu et al. 2004) is usually not as accurate as those in the highly bended cortical region with high curvature value, such as gyral crest and sulcal bottom (Li et al. 2009), the weight of the guidance force should be set larger at the highly bended region. Therefore, the weight of the guidance force w is defined as:
To prevent the cortical surface from self-collision and deforming into other brain tissues, such as the skull, basal ganglia/thalami, and cerebellum, these tissues and the current deformed cortical surface are first rasterized into a volumetric model (Nie et al. 2010). Then, the deformation of each cortical vertex is prevented from self-collision and deforming into these tissues in this volumetric model.
For each vertex with the current position , when it’s being deformed to a new position , 2 conditions should be checked: 1) The new position cannot be inside of the developing skull or other brain tissue volume as mentioned before and 2) The deformation of the current vertex should not cause its neighboring triangles to intersect with the rasterized surfaces in other neighborhoods. If any of the 2 conditions is violated, a new deformed position is identified as the closest valid position to on the straight line between and as defined below:
Note that, after updating the position of a vertex, its neighborhood is re-rasterized.
The explicit Newmark scheme (Newmark 1959) is adopted to solve the above dynamic model. Using the position and velocity given in the previous time t, the developed cortical surface at the next time can be estimated as:
Longitudinal Cortical Folding Measurement
With the proposed cortical growth model, we can measure and study the longitudinal changes of the cortical folding. Several quantitative methods for measuring cortical folding have been proposed in the literature. For example, the traditional gyrification index was firstly proposed in Zilles et al. (1988) to measure the cortical folding in a 2D slice and recently extended to the 3D local gyrification by measuring the cortical surface area in a sphere (Schaer et al. 2008; Toro et al. 2008). Meanwhile, curvature-based methods have also been proposed to measure the complexity of the cortical folding especially in the developing brain (Rodriguez-Carranza et al. 2007; Pienaar et al. 2008). Recent comparison on the curvature-based measurement and the gyrification index (Rodriguez-Carranza et al. 2007) also shows that these 2 types of measurements perform similarly on the inner cortical surfaces.
In this paper, 2 curvature-based measures, for example, the curvedness (Koenderink and van Doorn 1992) and sharpness (Pienaar et al. 2008), are adopted to characterize the local change of the cortical folding:
To measure the changes of curvedness and sharpness in cortical regions, the cortical surface is parcellated into a set of regions of interest (ROIs) by an atlas-based warping method (Shen and Davatzikos 2002), and the relative curvedness and sharpness changes in each ROI are defined, respectively, as:
Since the curvature decreases with the increase of the cortical surface area even in the case of linear cortical growth, the change of cortical surface area should be normalized in order to measure the cortical folding relatively. Therefore, the normalized relative curvedness and sharpness changes in each ROI can be defined as below:
The performance of our proposed cortical growth model is evaluated both qualitatively and quantitatively on 10 healthy subjects (each with longitudinal images collected at 2 weeks, 3, 6, 9, and 12 months old) and is further compared with a high-dimensional nonlinear hybrid (volumetric/surface) registration algorithm (Liu et al. 2004). In our method, we use the inner cortical surfaces at 6 and 12 months old to help build our cortical growth model and then use the inner cortical surfaces reconstructed at 3 and 9 months old as the ground truth to evaluate the prediction power of our model. Note that, after building our cortical growth model, given any time in the first year, we can immediately estimate its corresponding cortical surface, that is, at 3 or 9 months old. Similarly, for the hybrid registration algorithm, we can also build its respective brain growth model by registering the 2-week-old image to the 6-month-old image and then to the 12-month-old image. In this case, the brain image at 3 or 9 months old can be estimated by resampling the estimated deformation fields by constraining the respective cortical surface area to be the same as that in the ground truth 3- or 9-month-old cortical surface. It is worth noting that this cortical surface area information is not used in our model when estimating cortical surfaces at 3 and 9 months.
In all experiments, the same parameters are used, that is, Ke = 0.1, τe = 50, Kb = 0.1, τb = 50, g = 0.01, w1 = 0.006, w2 = 0.004, βth = 2.0, Cm = 0.7, α = 3.0, and Δt = 0.05. We use 240 iterations in building the cortical growth model for the first year, thus the real-world unit for time in our implementation is month (each with 20 iterations). Our model takes about 4 h to build the cortical growth model for each subject on a PC with Inter Xeon 2.26 Ghz CPU and 4 GB memory.
Visual Inspection of the Accuracy
Figure 3 gives a comparison between the cortical surfaces predicted by our model (red color) and the ground truth (light blue color) reconstructed at 3, 6, 9, and 12 months old. Note that the ground truth surfaces are set to be transparent for better visual inspection. With the cortical surface guidance at 6 and 12 months old, our predicted cortical surfaces have the similar shapes and positions to the ground truth at the respective time points, as shown in Figure 3b,d, respectively. Also, at other time points without any guidance, our predicted cortical surfaces still have the similar shapes and positions to the ground truth, as shown in Figure 3a,c.
Quantitative Evaluation of the Accuracy
To quantitatively evaluate our model and further compare its performance with the hybrid registration algorithm (Liu et al. 2004), we measure the surface distance error between the predicted surface and the ground truth surface at each time point (from 3 to 12 months old) on all 10 subjects. Note that the surface distance is symmetrically computed between the 2 surfaces under comparison. Figure 4 shows a comparison of surface distance errors between our method and the hybrid registration method on one subject at 4 different time points. It can be observed that the surface distance errors between the cortical surfaces predicted by our model and the ground truth cortical surfaces at 3, 6, 9, and 12 months old, as shown in Figure 4a–d, respectively, are much smaller than the surface distance errors between the cortical surfaces predicted by the hybrid registration algorithm and the ground truth cortical surfaces at 3, 6, 9, and 12 months old, as shown in Figure 4e–h, respectively.
The comparison of the averages and standard deviations of surface distance errors of 10 subjects between our model and the hybrid registration algorithm at 4 time points (from 3 to 12 months old) is shown in Figure 5. With the surface guidance at 6 and 12 months old, the average surface distance for all subjects is smaller than 0.25 mm by both our model and the hybrid registration algorithm, which also use the 6- and 12-month-old images and surfaces to guide the registration. For the time points without surface guidance such as at 3 and 9 months old, the average surface distances for all subjects are 0.48 ± 0.10 and 0.24 ± 0.06 mm by our model, respectively, and 0.68 ± 0.10 and 0.61 ± 0.10 mm by the hybrid registration algorithm, respectively. Considering the small brain size in the first year of life, this improvement is significant for early brain development study. This experiment demonstrates that our proposed model can predict more accurate cortical surfaces with use of the mechanical models compared with the hybrid registration algorithm that considers only the image- and surface-derived information in the registration.
Smoothness and Consistency of Cortex Growth on Vertices
We further visually evaluate the smoothness and consistency of the longitudinal cortical surfaces predicted by our model, by comparison with those by the hybrid registration algorithm. As shown in Figure 6a, the trajectories of vertex growth (red curves) predicted by our model along 5 time points (from 2 weeks to 12 months) are much smoother even in the flat cortical surface regions compared with those (blue curves) produced by the hybrid registration algorithm that registers the 2-week-old image/surface to images/surfaces at other 4 time points. To quantitatively compare the performance, we first perform the linear regression on the trajectory of each vertex and then show the histogram of the residual errors in the whole cortical surface. As shown in Figure 6b, our method gives smaller residual errors than the hybrid registration method.
Smoothness and Consistency of Cortex Growth on ROIs
To further compare the smoothness and consistency in larger cortical regions, an atlas with 90 labeled ROIs (Tzourio-Mazoyer et al. 2002) is warped onto the subject image space to label the subject cortical surface into a set of ROIs. For each ROI, the relative cortical area and the relative edge length are respectively defined as:
To quantify the smoothness and consistency of cortical folding features, the relative curvedness and the relative sharpness are also measured from the cortical surfaces predicted by our method (red) and the hybrid registration (blue) on 20 ROIs in Figure 8c,d. As we can see, both curvedness and sharpness are decreasing smoothly in most ROIs with the growth of the brain according to our model, while they are bumpy by the hybrid registration method.
We also show in Figure 9 the changes of the relative curvedness on 5 major gyri, including the precentral gyrus, postcentral gyrus, superior temporal gyrus, superior occipital gyrus, and superior frontal gyrus, on the 10 subjects by the 2 methods. Although both methods show the consistent decreasing trend of curvedness on all 5 gyral ROIs, the results from our method show much smoother results for all subjects.
Cortical Folding Pattern from 2 Weeks to 12 Months
Since the above results have already demonstrated the accuracy, smoothness, and consistency of our proposed method, we apply it to measure the cortical folding changes from 2 weeks to 12 months. Note that we use the inner cortical surfaces reconstructed from images at all 5 time points to help build our cortical growth model in this experiment.
Since the primary and secondary folding structures have been developed before birth, the simple increasing of the cortical surface area will reduce the averages of curvedness and sharpness from 2 weeks to 12 months as shown in Figure 8c,d, respectively. The normalized changes of curvedness and sharpness from 2 weeks to 12 months on the 10 subjects are illustrated in Figures 10 and 11, respectively.
As we can see, though the normalized changes of curvedness and sharpness are continuously increasing during the first year, the cortical folding shows different changing patterns in different lobes. For example, the frontal lobe shows rapidly increasing cortical folding during this period, which corresponds to high-level function development in this lobe, as indicated by red arrows in both Figures 10 and 11. While in the occipital and parietal lobes, the cortical folding is increasing relatively slower as indicated by blue arrows in both Figures 10 and 11. And in the temporal lobe, the folding changes are even more complex on the 3 major gyri. Specifically, in the first 3 months, the curvedness of the inferior temporal gyrus increases faster than the middle temporal gyrus and increases slower after 3 months as indicated by green arrows in Figure 10. On the contrary, the folding changing speed of the superior temporal gyrus continues to increase during the first year as indicated by purple arrows in Figure 11.
The cortical folding also shows different patterns in different time periods. For example, in the frontal lobe, the curvedness and sharpness increase fast in the first 3 months and then slow down from 3 to 9 months, as indicated by red arrows in both Figures 10 and 11. The average sharpness in both pre- and postcentral gyri increases consistently during the first year after birth, although the precentral gyrus shows stronger increasing folding patterns from 2 weeks to 6 months, as indicated by orange arrows in Figures 10 and 11.
Certain hemispheric asymmetry could also be observed in Figures 10 and 11. For example, the left hemisphere shows faster increase of sharpness on precentral gyri as indicated by light orange arrows in Figure 11 from 2 weeks to 6 months, while the right hemisphere shows faster increase of curvedness on inferior temporal gyrus as indicated by purple arrows in Figure 10 from 2 weeks to 6 months.
Since the number of subjects (only 10 subjects in our results) is limited in our experiment, the statistical significance of curvedness and sharpness changes is measured by the P value in nonparametric permutation tests (Nichols and Holmes 2002) on each ROI with the null hypothesis that the mean values of and are greater than 1, and the results are shown in Figures 12 and 13, respectively. These results show similar pattern as in Figures 10 and 11. However, it is very interesting that during the period from 6 to 9 months, most cortical ROIs show no significant increase of curvedness, especially in the right hemisphere as indicated by red arrows in Figure 12. And during the first 3 months, the increase of sharpness in the most areas of occipital and temporal lobe is not significant as indicated by red arrows in Figure 13.
Discussion and Conclusion
Due to lack of information in the smooth and flat cortical regions, traditional deformable registration method is limited to correctly estimate tissue deforation over time in these regions as shown in our result. Even though longitudinal regulations could be added to partially compensate for this lack of information and also generate longitudinal smooth and consistent deformation to some degree, our model that takes into account the physical properties of the growing brain is expected to perform better than the registration method that simply assumes smoothly varying deformations.
By coupling longitudinal image analysis with physical models of the cortical growth, we have presented a computational growth model for cortical development in the first year of life. Experimental results show its good performance in predicting cortical surface development at unknown time points in the first year, with the help of mechanical modeling. These results also demonstrate that the features of cortical surfaces, such as curvedness and sharpness, are much more smoothly and consistently estimated by our proposed model than those by the hybrid registration method. Certain interesting developmental patterns of cortical folding in the first year of life are also reported in our results: 1) the curvedness and sharpness decrease from 2 weeks to 12 months; 2) the frontal lobe shows rapidly increasing cortical folding during this period, while the cortical folding increases relatively slower in the occipital and parietal lobes; 3) during the first 3 months, the curvedness of the inferior temporal gyrus increases faster than the middle temporal gyrus and increases slower after 3 months, whereas the folding changing speed of the superior temporal gyrus continues to increase during the first year, and 4) the cortical folding of the frontal lobe also seems to develop fastest during the first 3 months and after that slows down. Due to the complex neurobiological processes involved in the cortex development during this stage, such as WM myelination, synapse development, and neuron dendritic projection, our growth model is not complex enough to cover all possible developing factors currently. For example, we currently use the spatially uniform growth function in our growth model since we do not know the spatial distribution of growth rate. And since the axon tension could also have effect on cortical folding (Van Essen 1997), it is possible to incorporate the axon tension into our model by making use of DTI data sets. In the future study, a more advanced and sophisticated model further involving WM and functional development will be developed to simulate this procedure.
Since only the inner cortical surface is adopted in our method to represent the geometry of the cortex currently, certain important features of the cortex could not be easily measured from our results, such as the cortical thickness (Fischl and Dale 2000) and gyrification index (Zilles et al. 1988), which are usually defined on the outer cortical surface. However, our current model can still provide abundant information regarding the cortex development from different levels and aspects. In the future, we will carry out more experiments and measurements to further evaluate our cortical growth model and apply it to a large early brain development study.
National Institutes of Health (EB006733, EB008374, EB009634, MH088520, NS055754, HD053000, and MH070890).
Conflict of Interest: None declared.