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 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.

Figure 1.

Longitudinal T1, T2, and FA images of the same subject along with its inner cortical surfaces reconstructed from tissue segmentation results.

Figure 1.

Longitudinal T1, T2, and FA images of the same subject along with its inner cortical surfaces reconstructed from tissue segmentation results.

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,1,2,3,4} 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).

Figure 2.

Illustration of the elastoplasticity properties defined on each edge (i, j), and also the bending properties for the edge (i, j) defined on the virtual edge (p, q). (a) The original shape of edge (i, j) and its 2 neighboring triangles Δpij and Δqij. (b) Edge (i, j) and its 2 neighboring triangles are deformed by, respectively, contracting the edge (i, j) and rotating the 2 triangles. (c) The plasticity of the original lengths l0i,j and l0p,q under the current deformation. (d) The elastic forces on the edge (i, j) and the virtual edge (p, q) under the current deformation.

Figure 2.

Illustration of the elastoplasticity properties defined on each edge (i, j), and also the bending properties for the edge (i, j) defined on the virtual edge (p, q). (a) The original shape of edge (i, j) and its 2 neighboring triangles Δpij and Δqij. (b) Edge (i, j) and its 2 neighboring triangles are deformed by, respectively, contracting the edge (i, j) and rotating the 2 triangles. (c) The plasticity of the original lengths l0i,j and l0p,q under the current deformation. (d) The elastic forces on the edge (i, j) and the virtual edge (p, q) under the current deformation.

Elastic property and plastic property. At the developmental stage t, the elastic property fe on each edge (i, j) of the triangulated inner surface is modeled as: 

(1)
fei,j(t)=Ke(l0i,j(t)lci,j(t))l0i,j(t)(xci(t)xcj(t))lci,j(t),
where 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)=xci(t)xcj(t) 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. Ke is the elastic constant that controls the elasticity of the edge (i, j). At the initial stage (t = 0), we have l0i,j(t)=xci(t=0)xcj(t=0), 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): 

(2)
dl0i,j(t)dt=1τe(lci,j(t)l0i,j(t)),
where τ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 fb on each virtual edge (p, q) is modeled as: 

(3)
fbp,q(t)=Kb(l0p,q(t)lcp,q(t))l0p,q(t)(xcp(t)xcq(t))lcp,q(t),
where l0p,q(t)is the rest length of the virtual edge (p, q), lcp,q(t)=xcp(t)xcq(t) is the current length of the virtual edge (p, q), and Kb is the elastic constant that controls the bending elasticity of the cortex. At the initial stage (t = 0), we have l0p,q(t)=xcp(t=0)xcq(t=0), 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): 
(4)
dl0p,q(t)dt=1τb(lcp,q(t)l0p,q(t)),
where τ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): 

(5)
dAc(t)dt=Ac(t)·g·(1Ac(t)AStk),tk1<t<tk,k≥1,
where Ac(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, k1}, 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 St0is 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) can be found on the target cortical surface Stk. Then, the guidance force fgi can be modeled as: 

(6)
fgi(t)=w((xi)(t)xci(t)),
where w is the weight of the guidance force varying at different cortical regions, and the setting of w will be detailed later. And (xi)(t)is the attraction point determined by the later developed cortical surfaces and computed as: 
(7)
(xi)(t)=β(t)·(xki)+(1β(t))·(xk+1i),tk1<t<tk,
where (xki) and (xk+1i) are the correspondences of vertex i at time points k and k + 1 determined by the registration method, respectively; and (xi) is a linear combination of these 2 correspondences weighted by β(t) as defined below: 
(8)
β(t)={1AStkAc(t)AStk+1AStk>βth1βth·AStkAc(t)AStk+1AStkAStkAc(t)AStk+1AStkβth,
where β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+1AStk), 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., AStkAc(t)AStk+1AStkβ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)(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: 

(9)
w=w1+w2tanh(α·(max(|Cmaxi|,|Cmini|)Cm)),
where w1, w2, and α are the weight values, Cmaxi and Cmini are the maximum and minimum principal curvatures of the vertex i, and Cm 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 w1 and w2 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+Δt), 2 conditions should be checked: 1) The new position xci(t+Δ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˜ci(t+Δt) is identified as the closest valid position to xci(t+Δt) on the straight line between xci(t) and xci(t+Δt) as defined below: 

(10)
x˜ci(t+Δt)=θxci(t)+(1−θ)xci(t+Δt),θ[0,1].

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: 

(11)
x¨=M1F(x,x˙),
where x, x˙, and x¨ are the position, velocity, and acceleration of all n vertices on the cortical surface, respectively. M is a 3n × 3n diagonal mass matrix with diag(M)=(m1,m1,m1,m2,m2,m2,,mn,mn,mn), where mi 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˙) is the force vector that combines all forces on vertices. For the vertex, i, Fi(x,x˙)is calculated as: 
(12)
Fi(x,x˙)=jNifei,j(t)+(p=i),qfbp,q(t)+fgi(t),
where Ni is the one-ring neighboring vertices of the vertex i, and (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+Δt can be estimated as: 

(13)
x˙(t+Δt)=x˙(t)+Δtx¨(t),
 
(14)
x(t+Δt)=x(t)+Δtx˙(t+Δt)+Δt22x¨(t),
where Δt is the time step. Note that, before estimating the cortical surface at new time, the rest lengths of each edge and each virtual edge are updated by equations (3–5).

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: 

(15)
ci(t)=(Cmaxi(t))2+(Cmini(t))222,
 
(16)
si(t)=(Cmaxi(t)Cmini(t))2,
where Cmaxi(t), Cmini(t), ci(t), and si(t) are the maximum principal curvature, minimum principal curvature, curvedness (Koenderink and van Doorn 1992), and sharpness (Pienaar et al. 2008) of vertex 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: 

(17)
RCROI(t)=iROIci(t)iROIci(t0),
 
(18)
RSROI(t)=iROIsi(t)iROIsi(t0),
where t0 is the initial time of the cortex development in the cortical growth model, RCROI(t) and RSROI(t) are, respectively, the average curvedness and sharpness in one ROI, normalized by their initial values at time t0 (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: 

(19)
NRCROI(tk,tk1)=iROIci(tk)iROIci(tk1)·AROI(tk)AROI(tk1),
 
(20)
NRSROI(tk,tk1)=iROIsi(tk)iROIsi(tk1)·AROI(tk)AROI(tk1),
where AROI(tk) is the area of ROI at time tk, NRCROI(tk,tk1) and NRSROI(tk,tk1) are the normalized change rate of curvedness and sharpness between times tk and tk1, respectively.

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, 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.

Figure 3.

Comparison between the cortical surfaces predicted by our model (red colors) and the ground truth cortical surfaces (light blue colors) at 3, 6, 9, and 12 months old, as shown from (a) to (d), respectively.

Figure 3.

Comparison between the cortical surfaces predicted by our model (red colors) and the ground truth cortical surfaces (light blue colors) at 3, 6, 9, and 12 months old, as shown from (a) to (d), respectively.

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.

Figure 4.

Comparison of surface distances errors (mm) of one subject at 4 time points from 3 to 12 months old. (ad) are 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, respectively. (eh) are 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, respectively. The color bars are shown on the right.

Figure 4.

Comparison of surface distances errors (mm) of one subject at 4 time points from 3 to 12 months old. (ad) are 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, respectively. (eh) are 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, respectively. The color bars are shown on the right.

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.

Figure 5.

Comparison of the averages and standard deviations of surface distance errors of 10 subjects between our model (red bars) and the hybrid registration algorithm (blue bars) at 4 time points from 3 to 12 months old.

Figure 5.

Comparison of the averages and standard deviations of surface distance errors of 10 subjects between our model (red bars) and the hybrid registration algorithm (blue bars) at 4 time points from 3 to 12 months old.

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.

Figure 6.

Comparison of the smoothness of vertex growth by our cortical growth model (red colors) and the hybrid registration algorithm (blue colors). (a) Trajectories of vertex growth by 2 methods. (b) The distribution of residual errors after performing linear regression on each trajectory by 2 methods.

Figure 6.

Comparison of the smoothness of vertex growth by our cortical growth model (red colors) and the hybrid registration algorithm (blue colors). (a) Trajectories of vertex growth by 2 methods. (b) The distribution of residual errors after performing linear regression on each trajectory by 2 methods.

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: 

(21)
RAi(t)=jTiAtrij(t)jTiAtrij(t0),
 
(22)
RLi(t)=jNilci,j(t)jNilci,j(t0),
where Ti is the one-ring neighboring triangles of vertex i, Ni 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.

Figure 7.

The growth of the relative local area (a) and the relative edge length (b) on 20 randomly selected vertices from 2 weeks to 12 months, by our model (red) and the hybrid registration method (blue).

Figure 7.

The growth of the relative local area (a) and the relative edge length (b) on 20 randomly selected vertices from 2 weeks to 12 months, by our model (red) and the hybrid registration method (blue).

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: 

(23)
RAROI(t)=jTROIAtrij(t)jTROIAtrij(t0),
 
(24)
RLROI(t)=jLROIledgej(t)jLROIledgej(t0),
where TROI and LROI 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.

Figure 8.

The changes of cortical features on 20 cortical ROIs from 2 weeks to 12 months. (a) Relative ROI area. (b) Relative ROI edge length. (c) Relative ROI curvedness. (d) Relative ROI sharpness. Red and blue curves represent the results by our model and the hybrid registration method, respectively.

Figure 8.

The changes of cortical features on 20 cortical ROIs from 2 weeks to 12 months. (a) Relative ROI area. (b) Relative ROI edge length. (c) Relative ROI curvedness. (d) Relative ROI sharpness. Red and blue curves represent the results by our model and the hybrid registration method, respectively.

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.

Figure 9.

The changes of average relative curvedness on 5 major gyral ROIs of the 10 subjects by our cortical growth model (red curves) and the hybrid registration method (blue curves), from 2 weeks to 12 months.

Figure 9.

The changes of average relative curvedness on 5 major gyral ROIs of the 10 subjects by our cortical growth model (red curves) and the hybrid registration method (blue curves), from 2 weeks to 12 months.

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,tk1) and sharpness NRSROI(tk,tk1) from 2 weeks to 12 months on the 10 subjects are illustrated in Figures 10 and 11, respectively.

Figure 10.

ROI-based average curvedness changes on the 10 subjects from 2 weeks to 12 months old. The changes from 2 weeks to 3 months, 3–6 months, 6–9 months, and 9–12 months are shown from the first to the fourth row, respectively. The color bars are given on the right. In each row, both side views and top view are displayed.

Figure 10.

ROI-based average curvedness changes on the 10 subjects from 2 weeks to 12 months old. The changes from 2 weeks to 3 months, 3–6 months, 6–9 months, and 9–12 months are shown from the first to the fourth row, respectively. The color bars are given on the right. In each row, both side views and top view are displayed.

Figure 11.

ROI-based average sharpness changes on the 10 subjects from 2 weeks to 12 months old. The changes from 2 weeks to 3 months, 3–6 months, 6–9 months, and 9–12 months are shown from the first to the fourth row, respectively. The color bars are given on the right. In each row, both side views and top view are displayed.

Figure 11.

ROI-based average sharpness changes on the 10 subjects from 2 weeks to 12 months old. The changes from 2 weeks to 3 months, 3–6 months, 6–9 months, and 9–12 months are shown from the first to the fourth row, respectively. The color bars are given on the right. In each row, both side views and top view are displayed.

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.

Figure 12.

P values of ROI-based curvedness changes on the 10 subjects, from 2 weeks to 12 months old. The results from 2 weeks to 3 months, 3–6 months, 6–9 months, and 9–12 months are shown from the first to the fourth row, respectively. The color bars are given on the right. In each row, both side views and top view are displayed.

Figure 12.

P values of ROI-based curvedness changes on the 10 subjects, from 2 weeks to 12 months old. The results from 2 weeks to 3 months, 3–6 months, 6–9 months, and 9–12 months are shown from the first to the fourth row, respectively. The color bars are given on the right. In each row, both side views and top view are displayed.

Figure 13.

P values of ROI-based sharpness changes on the 10 subjects, from 2 weeks to 12 months old. The results from 2 weeks to 3 months, 3–6 months, 6–9 months, and 9–12 months are shown from the first to the fourth row, respectively. The color bars are given on the right. In each row, both side views and top view are displayed.

Figure 13.

P values of ROI-based sharpness changes on the 10 subjects, from 2 weeks to 12 months old. The results from 2 weeks to 3 months, 3–6 months, 6–9 months, and 9–12 months are shown from the first to the fourth row, respectively. The color bars are given on the right. In each row, both side views and top view are displayed.

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.

References

Awate
SP
Yushkevich
PA
Song
Z
Licht
DJ
Gee
JC
Cerebral cortical folding analysis with multivariate modeling and testing: studies on gender differences and neonatal development
Neuroimage
 , 
2010
, vol. 
53
 (pg. 
450
-
459
)
Brown
M
Keynes
R
Lumsden
A
The developing brain
 , 
2002
Oxford
Oxford University Press
Dubois
J
Benders
M
Borradori-Tolsa
C
Cachia
A
Lazeyras
F
Ha-Vinh Leuchter
R
Sizonenko
SV
Warfield
SK
Mangin
JF
Huppi
PS
Primary cortical folding in the human newborn: an early marker of later functional development
Brain
 , 
2008
, vol. 
131
 (pg. 
2028
-
2041
)
Fischl
B
Dale
AM
Measuring the thickness of the human cerebral cortex from magnetic resonance images
Proc Natl Acad Sci U S A
 , 
2000
, vol. 
97
 (pg. 
11050
-
11055
)
Fischl
B
Sereno
MI
Dale
AM
Cortical surface-based analysis. II: Inflation, flattening, and a surface-based coordinate system
Neuroimage
 , 
1999
, vol. 
9
 (pg. 
195
-
207
)
Geng
G
Johnston
LA
Yan
E
Britto
JM
Smith
DW
Walker
DW
Egan
GF
Biomechanisms for modelling cerebral cortical folding
Med Image Anal
 , 
2009
, vol. 
13
 (pg. 
920
-
930
)
Grinspun
E
Hirani
AN
Desbrun
M
Schr
P
Discrete shells
Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation
 , 
2003
San Diego (CA)
Eurographics Association
(pg. 
62
-
67
)
Kaukola
T
Kapellou
O
Laroche
S
Counsell
SJ
Dyet
LE
Allsop
JM
Edwards
AD
Severity of perinatal illness and cerebral cortical growth in preterm infants
Acta Paediatr
 , 
2009
, vol. 
98
 (pg. 
990
-
995
)
Koenderink
JJ
van Doorn
AJ
Surface shape and curvature scales
Image Vis Comput
 , 
1992
, vol. 
10
 (pg. 
557
-
564
)
Li
G
Guo
L
Nie
J
Liu
T
Automatic cortical sulcal parcellation based on surface principal direction flow field tracking
Neuroimage
 , 
2009
, vol. 
46
 (pg. 
923
-
937
)
Liu
T
Li
H
Wong
K
Tarokh
A
Guo
L
Wong
S
Brain tissue segmentation based on DTI data
Neuroimage
 , 
2007
, vol. 
38
 
1
(pg. 
114
-
123
)
Liu
T
Shen
D
Davatzikos
C
Deformable registration of cortical structures via hybrid volumetric and surface warping
Neuroimage
 , 
2004
, vol. 
22
 (pg. 
1790
-
1801
)
Murray
JD
Mathematical biology
 , 
1993
2nd ed
Heidelberg (Germany)
Springer-Verlag
Newmark
NM
A method of computation for structural dynamics
J Eng Mech
 , 
1959
, vol. 
85
 (pg. 
67
-
94
)
Nichols
TE
Holmes
AP
Nonparametric permutation tests for functional neuroimaging: a primer with examples
Hum Brain Mapp
 , 
2002
, vol. 
15
 (pg. 
1
-
25
)
Nie
J
Guo
L
Li
G
Faraco
C
Stephen Miller
L
Liu
T
A computational model of cerebral cortex folding
J Theor Biol
 , 
2010
, vol. 
264
 (pg. 
467
-
478
)
Pienaar
R
Fischl
B
Caviness
V
Makris
N
Grant
PE
A methodology for analyzing curvature in the developing brain from preterm to adult
Int J Imaging Syst Technol
 , 
2008
, vol. 
18
 (pg. 
42
-
68
)
Rodriguez-Carranza
CE
Mukherjee
P
Vigneron
D
Barkovich
J
Studholme
C
Comparing 3D gyrification index and area-independent curvature-based measures in quantifying neonatal brain folding
Proc. SPIE
 , 
2007
, vol. 
6512
 pg. 
65120N
 
Schaer
M
Cuadra
MB
Tamarit
L
Lazeyras
F
Eliez
S
Thiran
JP
A surface-based approach to quantify local cortical gyrification
IEEE Trans Med Imaging
 , 
2008
, vol. 
27
 (pg. 
161
-
170
)
Schumann
CM
Bloss
CS
Barnes
CC
Wideman
GM
Carper
RA
Akshoomoff
N
Pierce
K
Hagler
D
Schork
N
Lord
C
, et al.  . 
Longitudinal magnetic resonance imaging study of cortical development through early childhood in autism
J Neurosci
 , 
2010
, vol. 
30
 (pg. 
4419
-
4427
)
Shattuck
DW
Leahy
RM
Automated graph-based analysis and correction of cortical volume topology
IEEE Trans Med Imaging
 , 
2001
, vol. 
20
 (pg. 
1167
-
1177
)
Shaw
P
Eckstrand
K
Sharp
W
Blumenthal
J
Lerch
JP
Greenstein
D
Clasen
L
Evans
A
Giedd
J
Rapoport
JL
Attention-deficit/hyperactivity disorder is characterized by a delay in cortical maturation
Proc Natl Acad Sci U S A
 , 
2007
, vol. 
104
 (pg. 
19649
-
19654
)
Shen
D
Davatzikos
C
HAMMER: hierarchical attribute matching mechanism for elastic registration
IEEE Trans Med Imaging
 , 
2002
, vol. 
21
 (pg. 
1421
-
1439
)
Shen
EY
Wu
KH
Lin
MF
Chen
CY
Study of brain growth in children—a new approach to volume measurements using MRI-reconstructed 3D neuroimaging
Childs Nerv Syst
 , 
2010
, vol. 
26
 (pg. 
1619
-
1623
)
Talairach
J
Tournoux
P
Co-planar stereotaxic atlas of the human brain
 , 
1988
Stuttgart (Germany)
Thieme
Toro
R
Burnod
Y
A morphogenetic model for the development of cortical convolutions
Cereb Cortex
 , 
2005
, vol. 
15
 (pg. 
1900
-
1913
)
Toro
R
Perron
M
Pike
B
Richer
L
Veillette
S
Pausova
Z
Paus
T
Brain size and folding of the human cerebral cortex
Cereb Cortex
 , 
2008
, vol. 
18
 (pg. 
2352
-
2357
)
Tzourio-Mazoyer
N
Landeau
B
Papathanassiou
D
Crivello
F
Etard
O
Delcroix
N
Mazoyer
B
Joliot
M
Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single-subject brain
Neuroimage
 , 
2002
, vol. 
15
 (pg. 
273
-
289
)
Van Essen
DC
A tension-based theory of morphogenesis and compact wiring in the central nervous system
Nature
 , 
1997
, vol. 
385
 (pg. 
313
-
318
)
Van Essen
DC
Maunsell
JH
Two-dimensional maps of the cerebral cortex
J Comp Neurol
 , 
1980
, vol. 
191
 (pg. 
255
-
281
)
Vercauteren
T
Pennec
X
Perchant
A
Ayache
N
Diffeomorphic demons: efficient non-parametric image registration
Neuroimage
 , 
2009
, vol. 
45
 (pg. 
S61
-
S72
)
Wang
L
Shi
F
Yap
P-T
Gilmore
JH
Lin
W
Shen
D
Accurate and Consistent 4D Segmentation of Serial Infant Brain MR Images
LNCS
 , 
2011
, vol. 
7012
 (pg. 
93
-
101
)
Woods
RP
Grafton
ST
Watson
JD
Sicotte
NL
Mazziotta
JC
Automated image registration: II. Intersubject validation of linear and nonlinear models
J Comput Assist Tomogr
 , 
1998
, vol. 
22
 (pg. 
153
-
165
)
Xu
C
Pham
DL
Rettmann
ME
Yu
DN
Prince
JL
Reconstruction of the human cerebral cortex from magnetic resonance images
IEEE Trans Med Imaging
 , 
1999
, vol. 
18
 (pg. 
467
-
480
)
Zilles
K
Armstrong
E
Schleicher
A
Kretschmann
HJ
The human pattern of gyrification in the cerebral cortex
Anat Embryol (Berl)
 , 
1988
, vol. 
179
 (pg. 
173
-
179
)