## Abstract

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.

## Introduction

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 mm^{3}; 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 mm^{3}; 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/mm^{2}. T2 image is linearly aligned onto the respective T1 image of the same subject and further resampled to 1 × 1 × 1 mm^{3}. 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 mm^{3}. 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 $St0$, $St1$, $St2$, $St3$, and $St4$ 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 ${Stk,k=0,\u20021,\u20022,\u20023,\u20024}$ 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 Δp*ij* 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).

**Elastic property and plastic property**. At the developmental stage *t*, the elastic property *f*_{e} on each edge (*i*, *j*) of the triangulated inner surface is modeled as:

*i*and

*j*are the 2 ending points on the edge (

*i*,

*j*) as shown in Figure 2a, $l0i,j(t)$is the rest length of the edge (

*i*,

*j*), $lci,j(t)=\u2225xci(t)\u2212xcj(t)\u2225$ is the current length of the edge (

*i*,

*j*), and $xci(t)$ and $xcj(t)$ are the current coordinates of the vertices

*i*and

*j*, respectively.

*K*

_{e}is the elastic constant that controls the elasticity of the edge (

*i*,

*j*). At the initial stage (

*t*= 0), we have $l0i,j(t)=\u2225xci(t=0)\u2212xcj(t=0)\u2225$, where $xci(t=0)$ and $xcj(t=0)$are the original coordinates of vertices

*i*and

*j*on the inner cortical surface $St0$(as described above), respectively. Note that this elastic property will drive each edge to restore to its rest length during the cortical development, as shown in Figure 2d.

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 $l0i,j(t)$to the deformed length $lci,j(t)$:

_{e}is the time constant for the plasticity (Toro and Burnod 2005).

**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*, $q$) 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 *f*_{b} on each virtual edge (*p*, *q*) is modeled as:

*p*,

*q*), $lcp,q(t)=\u2225xcp(t)\u2212xcq(t)\u2225$ is the current length of the virtual edge (

*p*,

*q*), and

*K*

_{b}is the elastic constant that controls the bending elasticity of the cortex. At the initial stage (

*t*= 0), we have $l0p,q(t)=\u2225xcp(t=0)\u2212xcq(t=0)\u2225$, where $xcp(t=0)$and $xcq(t=0)$are the original coordinates of vertices

*p*and

*q*on the cortical surface $St0$. Similar to equation (2), the bending plasticity property can be modeled as the adaptation of the rest length $l0p,q(t)$ to the deformed length $lcp,q(t)$:

_{b}is the time constant for the bending plasticity property (Nie et al. 2010).

#### 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 $AStk$, the cortex growth model at time *t* can be described as the classic logistic growth function (Murray 1993; Toro and Burnod 2005):

*A*

_{c}(

*t*) is the current area of the cortex, and

*g*is known as the Malthusian parameter, and $AStk$ is the carrying capacity of the system (Murray 1993), which is defined as the surface area of the later developed cortex $Stk$. By changing the rest area of each triangle, the rest lengths of the 3 corresponding edges on the triangle will also change accordingly. Thus, the cortical growth will generate a mechanical stress to partially drive the deformation of the cortical surface. At the initial stage, we increase the rest length of each edge element according to the cortical growth function, thus the initial stress could be calculated according to the difference between the current length and rest length of each element.

#### Guidance Model

Since we have reconstructed cortical surfaces from the later developed stages ${Stk, tk-1<t<tk, k\u22651}$, 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 $Stk$. Second, the simulation to the realistic cortical shape and position can be guided by the shape and position of the cortical surface $Stk$.

The initial cortical surface $St0$is first warped into the space of the target cortical surface $Stk$ (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 $xci(t=0)$ on the initial cortical surface $St0$, a rough corresponding vertex $(xki)\u2032$ can be found on the target cortical surface $Stk$. Then, the guidance force $fgi$ can be modeled as:

*w*is the weight of the guidance force varying at different cortical regions, and the setting of

*w*will be detailed later. And $(xi)\u2032(t)$is the attraction point determined by the later developed cortical surfaces and computed as:

*i*at time points

*k*and

*k*+ 1 determined by the registration method, respectively; and $(xi)\u2032$ is a linear combination of these 2 correspondences weighted by β(

*t*) as defined below:

_{th}is a threshold controlling when the surface $Stk+1$ could be involved into the guidance model. Specifically, when comparing to the area difference between the next 2 target surfaces $(AStk+1\u2212AStk)$, if the area of the current surface $(Ac(t))$ is much smaller than the area of the next target surface $(AStk)$, only the next target surface $Stk$ will be used to guide the simulation (β(

*t*) = 1). Otherwise, if the current surface area is close to the next target surface area (i.e., $AStk\u2212Ac(t)AStk+1\u2212AStk\u2264\beta th$), the correspondences on the surface at time

*k*+ 1 could also be involved to further guide the simulation. In this way, the corresponding vertex $(xi)\u2032(t)$ will smoothly change across different time points.

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:

*w*

_{1},

*w*

_{2}, and α are the weight values, $Cmaxi$ and $Cmini$ are the maximum and minimum principal curvatures of the vertex

*i*, and

*C*

_{m}is a curvature constant. Thus, the guidance model would be smooth during the simulation. Since the force $fgi$ is an artificial force that does not truly exist, the weights

*w*

_{1}and

*w*

_{2}are both set as small values during the simulation.

#### Constraints

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 $xci(t)$, when it’s being deformed to a new position $xci(t+\Delta t)$, 2 conditions should be checked: 1) The new position $xci(t+\Delta t)$ 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 $x\u02dcci(t+\Delta t)$ is identified as the closest valid position to $xci(t+\Delta t)$ on the straight line between $xci(t)$ and $xci(t+\Delta t)$ as defined below:

Note that, after updating the position of a vertex, its neighborhood is re-rasterized.

#### Model Solver

Since we are applying a mechanical model on the cortex, the dynamics of the cortical surface can be formulated as Newton’s laws:

**x**, $x\u02d9$, and $x\xa8$ are the position, velocity, and acceleration of all

*n*vertices on the cortical surface, respectively.

**M**is a 3

*n*× 3

*n*diagonal mass matrix with $diag(M)=(m1,m1,m1,m2,m2,m2,\u2026,mn,mn,mn),$ where

*m*

_{i}is the mass of the vertex

*i*, calculated as the sum of one thirds of the masses of all one-ring triangles around the vertex

*i*. $F(x,x\u02d9)$ is the force vector that combines all forces on vertices. For the vertex,

*i*, $Fi(x,x\u02d9)$is calculated as:

*i*, and $\u2211(p=i),qfbp,q(t)$ is the sum of all bending forces involved on the vertex

*i*.

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 $t+\Delta t$ 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:

*i*at time

*t*, respectively. The sharpness $si(t)$ which emphasizes the difference between 2 principle curvatures is suitable to reveal the large primary folds of the cortex, while the curvedness $ci(t)$ can reveal small bumps and ridges of the cortex (Pienaar et al. 2008).

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:

*t*

_{0}(in order to account for the initial differences among subjects).

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:

## Results

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, *K*_{e} = 0.1, τ_{e} = 50, *K*_{b} = 0.1, τ_{b} = 50, *g* = 0.01, *w*_{1} = 0.006, *w*_{2} = 0.004, β_{th} = 2.0, *C*_{m} = 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.

We also measure the smoothness and consistency between neighboring vertices. The relative local area $RAi(t)$ and relative edge length $RLi(t)$ of each vertex *i* at time *t* are defined as:

*T*

^{i}is the one-ring neighboring triangles of vertex

*i*,

*N*

^{i}is the one-ring neighboring vertices of vertex

*i*, and $Atrij(t)$ is the area of triangle

*j*at time

*t*. Thus, $RAi(t)$ is the sum of triangular areas around vertex

*i*normalized by the initial area, and $RLi(t)$ is the sum of vertex distances between vertex

*i*and its neighboring vertices normalized by the initial distance. We compute $RAi(t)$ and $RLi(t)$, respectively, from the cortical surfaces predicted by our model from 2 weeks to 12 months and also by the hybrid registration algorithm that registers the 2-week-old image/surface to images/surfaces at other 4 time points. In Figure 7a,b, the respective results by the 2 methods on 20 randomly selected vertices are shown, with red for our model and blue for the hybrid registration. Compared with the hybrid registration, our method can predict much smoother growth for the relative local area and the relative edge length from 2 weeks to 12 months.

### 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 $RAROI(t)$ and the relative edge length $RLROI(t)$ are respectively defined as:

*T*

_{ROI}and

*L*

_{ROI}are the set of triangles and edges in the ROI, respectively, $ledgej(t)$ is the length of edge

*j*at time

*t*, $RAROI(t)$ is the area of the ROI normalized by its initial area, and $RLROI(t)$ is the sum of all edge lengths in the ROI normalized by its initial value. In Figure 8a,b, the results on 20 randomly selected cortical ROIs (from 90 ROIs) are given for our method (red) and the hybrid registration (blue). As we can see, our model shows much more consistent growth results.

To quantify the smoothness and consistency of cortical folding features, the relative curvedness $RCROI(t)$ and the relative sharpness $RSROI(t)$ 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 $RCROI(t)$ and sharpness $RSROI(t)$ 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 $NRCROI(tk,tk\u22121)$ and sharpness $NRSROI(tk,tk\u22121)$ 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.

## Funding

National Institutes of Health (EB006733, EB008374, EB009634, MH088520, NS055754, HD053000, and MH070890).

*Conflict of Interest*: None declared.