Regional growth and atlasing of the developing human brain

Detailed morphometric analysis of the neonatal brain is required to characterise brain development and de ﬁ ne neuroimaging biomarkers related to impaired brain growth. Accurate automatic segmentation of neonatal brain MRI is a prerequisite to analyse large datasets. We have previously presented an accurate and robust auto- maticsegmentationtechniqueforparcellatingtheneonatalbrainintomultiplecorticalandsubcorticalregions.In thisstudy,wefurtherextendoursegmentationmethodtodetectcorticalsulciandprovideadetaileddelineation ofthecorticalribbon.Thesedetailed segmentations are usedto builda4-dimensionalspatio-temporalstructural atlasofthebrainfor82corticalandsubcorticalstructuresthroughoutthisdevelopmentalperiod.Weemploythe algorithmtosegmentanextensivedatabaseof420MRimagesofthedevelopingbrain,from27to45weekspost-menstrualageatimaging.Regionalvolumetricandcorticalsurfacemeasurementsarederivedandusedtoinves- tigate brain growth and development during this critical period and to assess the impact of immaturity at birth. creased with increasing age at scan. Relative volumes of cortical grey matter, cerebellum and cerebrospinal ﬂ uid increased with age at scan, while relative volumes of white matter, ventricles, brainstem and basal ganglia and thalami decreased. Preterm infants at term had smaller whole brain volumes, reduced regional white matter and cortical and subcortical grey matter volumes, and reduced cortical surface area compared with term born controls, while ventricular volumewasgreater inthe preterm group. Increasing prematurity atbirth wasassoci-atedwithareductionintotalandregionalwhitematter,corticalandsubcorticalgreymattervolume,anincrease in ventricular volume, and reduced cortical surface area.


Introduction
The incidence of preterm birth continues to rise, with an estimated 14.9 million infants (representing 11.1% of all births) delivered worldwide each year (Blencowe et al., 2012). Insights into impaired neurodevelopment in these vulnerable infants have been gained from magnetic resonance (MR) imaging studies assessing brain development during the period between preterm birth and the normal time of birth. Brain volumetric and cortical surface measurements provide important information regarding brain development and have the potential to predict long-term neurodevelopmental performance (Peterson et al., 2003;Counsell et al., 2008;Thompson et al., 2008;Rathbone et al., 2011;Boardman et al., 2010). Previous studies in preterm infants have demonstrated reduced brain volume (Peterson et al., 2003;Inder et al., 2005;Thompson et al., 2007;Ball et al., 2012) and decreased cortical surface area (Ajayi-Obe et al., 2000;Kapellou et al., 2006), which are related to subsequent adverse neurodevelopmental outcome (Kapellou et al., 2006;Rathbone et al., 2011). However, studies to date have focused on characterising brain tissue volumes or volumes of large subcortical structures. In addition, sample sizes have usually been small, over a limited age range and detailed regional brain growth has not been studied (Hüppi et al., 1998;Murphy et al., 2001;Peterson et al., 2003;Inder et al., 2005;Prastawa et al., 2005;Mewes et al., 2006;Nishida et al., 2006;Zacharia et al., 2006;Gilmore et al., 2007;Song et al., 2007;Thompson et al., 2007;Xue et al., 2007;Anbeek et al., 2008;Dubois et al., 2008a,b;Pienaar et al., 2008;Rodriguez-Carranza et al., 2008;Yu et al., 2010;Cardoso et al., 2013;Wang et al., 2012;Moeskops et al., 2013Moeskops et al., , 2015. Due to the lack of manually segmented atlases, previous automatic methods (Peterson et al., 2003;Mewes et al., 2006;Gilmore et al., 2007;Thompson et al., 2007) did not segment deep grey matter structures and parcellated cortical grey matter (CGM) and white matter (WM) regions according to arbitrary linear parcellations which did not reflect regional anatomy. The first regional atlases of the brain were manually delineated by Gousias et al. (2012) that define 50 brain regions in 20 term-equivalent neonatal brains. In Makropoulos et al. (2014) we presented the first study to automatically segment the neonatal brain from the early preterm period to term-equivalent age, into 50 structures with the use of these atlases (82 regions when the WM/CGM regions defined in the atlases are further subdivided into the corresponding WM and CGM parts).
Accurate delineation of the cortex in the neonatal brain is challenging due to partial volume effects and limits in MR imaging resolution. The interior cortical boundary is difficult to delineate as CGM-WM partial volume (PV) can lead to overestimation of the segmented CGM. Furthermore, accurate delineation of the exterior cortical boundary is challenging as the complexity of the cortical surface in conjunction with limits in the MR imaging resolution renders the sulci delineation problematic. Only a few segmentation approaches focus on delineating the cortical ribbon in terms of morphology in the neonatal population (Xue et al., 2007;Wang et al., 2011Wang et al., , 2013. However, surface measurements are affected from undetected sulci, leading to increased cortical thickness estimates. These technical challenges have resulted in a limited number of studies assessing cortical surface measurements in the neonatal brain. Cortical surface area measurements have been reported in Kapellou et al. (2006), Xue et al. (2007), Dubois et al. (2008a,b), Pienaar et al. (2008), Rodriguez-Carranza et al. (2008), Rathbone et al. (2011), Lefevre et al. (2015), Moeskops et al. (2015) and curvature measurements/gyrification indices in Xue et al. (2007), Pienaar et al. (2008), Rodriguez-Carranza et al. (2008), Moeskops et al. (2013), Wang et al. (2013), Lefevre et al. (2015), Moeskops et al. (2015).
This study: • proposes a novel methodology for delineating the cortical ribbon. • assesses regional brain growth in the developing preterm brain.
• constructs a 4-dimensional spatio-temporal structural atlas with 82 labelled structures of the developing brain. • investigates the effect of prematurity on regional brain growth and cortical development.

Subjects
MPRAGE and T2-weighted MR imaging was performed on a Philips 3 Tesla Achieva system sited on the neonatal intensive care unit using an 8 channel phased array head coil. A cohort of 380 infants recruited from the Neonatal Intensive Care Unit at Queen Charlotte's and Hammersmith Hospitals was used in this study. Exclusion criteria for this study were focal abnormalities on MR imaging as determined by an expert perinatal neuro-radiologist. Out of the 380 infants in the cohort (467 images), 14 infants (19 images) were excluded due to abnormality: 4 infants with cystic periventricular leukomalacia (PVL) (7 images), 5 with Hemorrhagic Parenchymal Infarction (HPI) (6 images), 2 with pseudocysts (3 images), 1 with cerebellar haemorrhage, 1 with multiple white matter infractions and 1 with multiple cystic lesions following meningitis. 28 images were further excluded due to motion artefacts. The resulting dataset was used for the analysis in this work.
We studied 338 infants (298 preterm infants and 40 healthy term born infants) who were born at a median (range) gestational age (GA) of 30 (23 +2 -42) weeks. 49 preterm infants had chronic lung disease (CLD), 49 patent ductus arteriosus (PDA) and 32 culture positive sepsis. 66 preterm infants were scanned on 2 occasions and 8 on 3 occasions, resulting in 420 MR imaging data-sets. The interval between scans ranges from 6 days to 16 weeks, with a median of 8 +1 weeks. The term born infants included in this study were scanned only once. The median PMA at imaging was 39 +1 (27 +1 -44 +6 ) weeks and postnatal age at scan 5 +3 (0-19 +5 ) weeks. There were no significant differences in GA at birth, incidence of sepsis, PDA or CLD between the study group and those excluded due to focal lesions or motion corrupt images. There was no significant difference in age at scan between the study group and those with focal lesions, however the motion corrupted images were acquired at a significantly younger age at scan (p = 0.0016). The characteristics of the term-born and preterm infants are presented in Table 1.

MR imaging acquisition
All examinations were supervised by a paediatrician experienced in MR imaging procedures. Preterm infants at term age were sedated with oral chloral hydrate (25-50 mg/kg) prior to scanning and pulse oximetry, temperature and electrocardiography were monitored. Ear protection was used, comprising earplugs moulded from a silicone-based putty (President Putty, Coltene Whaledent, Mahwah, NJ, USA) placed in the external auditory meatus and neonatal earmuffs (MiniMuffs, Natus Medical Inc., San Carlos, CA, USA). T2-weighted images (TR 8670 ms; TE 160 ms; flip angle 90°; slice thickness 2 mm acquired with an overlap of 1 mm; in plane resolution 0.86 × 0.86 mm) were used for segmentation and cortical analysis.

Image segmentation
The T2 images were segmented with the pipeline presented in Makropoulos et al. (2014). The algorithm is based on an Expectation-Maximization (EM) scheme similar to Van Leemput et al. (1999) with a spatial prior term and an intensity model of the image. Spatial priors of the structures are obtained by averaging the warped labels from the 20 manually segmented atlases of Gousias et al. (2012). Image intensities are modelled with a Gaussian Mixture Model (GMM). Proposed adaptations of the segmentation model limit the influence of intensity in the delineation of structures with very similar intensity profiles. An extensive validation was performed to verify the robustness of the algorithm at different ages of the developing neonatal brain. This method allows regional brain growth in neonates to be assessed in an automatic and reproducible way. We refer the reader to Makropoulos et al. (2014) for more details on the individual parts of the pipeline. Table 2 presents the automatically parcellated regions.
In the next sections we present two extensions for the detailed delineation of the cortical ribbon in the neonatal brain. Examples of the proposed segmentation technique applied to an early preterm and term-equivalent brain can be seen in Figs. 1 and 2 respectively.

CGM-WM partial volume correction
Due to the partial volume between WM and CGM at the interface between the two tissues, automatic techniques tend to overestimate the CGM volume . Fig. 3 depicts this effect. The voxels between WM and CGM have an intermediate intensity and it is difficult to attribute them to either tissue. A Gaussian Mixture Model (GMM) with single classes for WM and CGM tends to overestimate the CGM extent. To account for this effect in Makropoulos et al. (2012) we implemented a partial volume correction for the CGM-WM boundary. In this study, an additional class is added to the GMM as a partial volume class between CGM and WM. Once the EM scheme has converged, the PV class is merged with the WM class to reduce the CGM overestimation and enhance the WM tissue estimate.

Sulci detection and enhancement
Delineation of the sulci in the neonatal cortex is difficult due to the limited resolution that leads to partial volume effects in the exterior cortical surface (see Fig. 4). Here, we present a novel approach for sulci enhancement for the neonatal images, based on the assumption that cortical thickness at sulci locations should be similar over a local neighbourhood of the cortical ribbon.
Sulci detection is performed in a way similar to Han et al. (2004). The CGM-WM interface (interior cortical surface) is iteratively propagated with the fast marching method with a speed function derived from the cerebrospinal fluid (CSF) posterior obtained with EM. Shock points are then detected where two sulcal banks merge in the propagated surface as in Han et al. (2004). Due to the limited resolution, the cortical ribbon in the two hemispheres may appear to be connected in different parts of the midsection of the brain (an example can be seen in Fig. 5.A). A second type of shock points is also added here for neighbouring CGM voxels from different hemispheres. Having detected the shock points, Han et al. (2004) perform morphological thinning to create a thin layer, one voxel thick, of CSF that splits the two sulcal banks apart. However, defining the CSF inside the sulcus to be one voxel wide might lead to incorrect cortical thickness estimates at points in the sulcal banks. As can be seen in Fig. 5.B, since the cortical thickness is estimated between the CGM-WM interface and the CGM-CSF interface, the width of the layer chosen for the shock points has a direct effect on the cortical thickness in the sulcal regions.
In this work the shock voxels are attributed to CSF only if their distance to the CGM-WM interface-the potential thickness if they belong to CGM-is larger than the cortical thickness of neighbouring parts of the cortical ribbon. Voxelwise cortical thickness is estimated as described in Jones et al. (2000). The sulcal points are thus prevented from thickness inconsistent with the rest of the cortex, since their cortical thickness is approximated based on the thickness of neighbouring parts in the cortical ribbon. Further details of the proposed sulci correction method can be found in Appendix A. An illustration of the method is presented in Fig. 6.

Cortical surface reconstruction
Cortical surface meshes were obtained by triangulation of the CGM-WM isosurface with the marching cubes algorithm (Lorensen and Cline, 1987). The CGM binary mask was initially blurred with a Gaussian kernel of 1 mm standard deviation to avoid a blocky surface due to the limited resolution. Laplacian smoothing was applied to the surfaces (Herrmann, 1976) to improve the mesh quality. The surface region that belongs to the boundary between WM and deep GM was excluded from the cortical surface (see Fig. 7). Example cortical surfaces are presented in Fig. 8.
Topological correction of the surfaces was not addressed. The interior cortical surface does not suffer from the problem of merging gyri and therefore does not present major topological defects.

Spatiotemporal structural atlas construction
Structural information of probabilistic brain atlases is constructed by averaging segmentations of different brains in the same coordinate space, in order to account for the anatomical variability in the brain. This section describes the construction of the first regional spatiotemporal structural atlas for the neonatal brain with a methodology similar to Kuklisova-Murgasova et al. (2011), andSerag et al. (2012). As in Kuklisova-Murgasova et al. (2011), andSerag et al. (2012), the segmentations are averaged with a non-parametric kernel regression according to the age at scan of the subjects.
The spatio-temporal template of Serag et al. (2012) is used as the coordinate space of the atlas. This template defines mean brain images for 28 to 44 weeks post-menstrual age at scan (with a week interval). The derived atlas is defined in the same age range. The segmentations of the 420 T2 images are warped to the space of the template according to the age at scan of each subject. In order to enforce consistency of the atlas in the time domain (different ages of the template), each segmentation is warped to the mean images in the range [a − 3,a + 3], where a is the rounded age at scan of the subject in weeks. The transformations are calculated with non-rigid registration of the subject's T2 image to the corresponding mean images of the template. The nonrigid registration is carried out using free-form deformations with control point grid spacings of 20 mm, 10 mm, 5 mm and 2.5 mm and normalized mutual information (NMI) as similarity measure Rueckert et al. (1999). Having computed the transformations for all the subjects, the age-dependent probability map P k,t of each structure k at time t of the atlas can be computed as: where s denotes the different subjects (S = 420 in total). M s,k is the binary mask of structure k from the segmentation of subject s. M s,k is warped (∘) under the transformation T s,t of s to the mean image of the template at age t. w(ts,t) is an age-dependent weight of subject s according to how closely the age t s of the subject matches the age t of the atlas. The weight is defined according to a Gaussian kernel: where σ w is set to 1 week. The probabilistic atlas at each age t is then defined as the union of the probability maps at the corresponding age: P t = {P k,t }. Fig. 4. Axial slice of a T2-weighted MRI (A) and magnified region of the cortex (B). Due to PV effects, the CSF inside the cortical sulci is often hard to discriminate, and consequently delineate with intensity-based segmentation techniques. Especially in areas where cortical gyri "touch" each other there is often very little evidence, in terms of intensity, of the CSF inside the sulcus. A maximum-probability version of the atlas at time t is further constructed by assigning the structure with the maximum probability to each voxel at age t: The constructed structural atlas incorporates the 82 structures of the brain, the CSF, the intra-cranial and the extra-cranial background. Illustrations of the probabilistic and maximum-probability versions of the atlas are presented in Figs. 9, 10. The atlas is publicly available from http://brain-development.org/.

Volumetric measurements
Absolute and relative volumes of the tissues and 82 structures of the brain were measured directly from the segmentations. Relative volume was determined as the ratio of the structure volume to the total brain volume (excluding the CSF and ventricles).

Cortical surface measurements
Surface area and curvature measures of the cortex were computed from the cortical surface meshes. A number of curvature measures were adopted from Rodriguez-Carranza et al. (2008) with Tnormalization, effectively normalized according to T ¼ 3 volume surface area , that are invariant to the surface area of the brain. This allows comparison of brains with different sizes, as is the case of the developing neonatal brain. The curvature measures included in this study are: global curvedness, mean curvature L 2 norm and Gaussian curvature L 2 norm. Their formulation is presented in Table 3. Regional cortical surface measurements were measured based on the segmented CGM structures. The cortical labels were propagated to the surface meshes. Each vertex of the mesh was labelled with the closest CGM structure in the 3-dimensional space.

Statistical analysis
We determined the centiles of the volumetric and surface measurements with increasing age at scan for the preterm datasets. Correlations with the age at scan were calculated to investigate premature brain development. The measurements were assessed individually according to the Pearson correlation coefficient and adjusted R 2 values from fitting a linear model to the data. Since the relationship between regional brain volumes and brain development may not always be characterised Fig. 5. A: T2 with the cortical segmentation overlaid. The arrows show parts of the cortical ribbon connected across the two hemispheres in the midsection of the brain. B: Example shock points (in pink) detected for the cortical segmentation (in red). Shock voxels are labelled as CSF if their distance D WM,i to the WM is larger than D allowed . D allowed is estimated from neighbouring parts of the cortical ribbon with streamlines that do not cross shock points (yellow lines).  with a linear model, we additionally use a Gompertz-like function to model this relationship, similar to Wright et al. (2014). We compare the fit of the two models according to the sum of squares error (SSE). The Gompertz function used in this study allows for both growth and decline to be studied and is detailed in the following section. The contribution of each measurement to brain development was further assessed with multiple linear regression of all the measurements of the same type (e.g. volume, surface area) against the age at scan.
The effect of preterm birth was assessed by comparing the group of term controls with an equal-sized group of early preterm subjects, born before 30 weeks, with equivalent ages at scan. Only a single scan per subject was included in the group analysis. There was no significant difference at the age at scan of the two groups (unpaired t-test, p = 0.98), however 10 of the preterm subjects had chronic lung disease. Group comparison was performed with unpaired t-tests when the data were normally distributed. Data that did not approximate a normal distribution were log transformed and checked again for normality. If the data were normally distributed in the log domain, comparison was again performed with an unpaired t-test in the log transformation. Nonnormal distributions were compared with the non-parametric  Wilcoxon rank sum test. Normality was assessed with the Lilliefor's test. Due to the large number of regions (82 in total) and the small number of term controls used in this study, we assess the effect of preterm birth using only univariate statistics.
The impact of increasing prematurity was further explored for all the preterm subjects by assessing correlations of the volume and surface measurements with GA at birth, correcting for the PMA at scan (correlation coefficient and adjusted R 2 values). The measurements were further entered in a multiple regression model with GA at birth (correcting for the PMA at scan) to explore their relative importance in a combined model.
We assessed the effect of including multiple scans of preterm infants in the analysis with a linear mixed-effects model. Linear mixedeffects have been used in the literature for analysis of longitudinal data (Bernal-Rusiel et al., 2013;Sadeghi et al., 2013). We examined the relationship with brain development and increasing prematurity by entering the age at scan and age at birth (correcting for age at scan) respectively as fixed effects into the model. In order to account for the correlation of measurements within the same subject, we included intercepts for the subjects as random effects. After fitting the model, we examined the significance of the random effect with a likelihood ratio test of the linear mixed model and a linear model (the equivalent of the linear mixed model without the random effect). Marginal R 2 values (Nakagawa and Schielzeth, 2013) of the fixed effects are reported and compared with the adjusted R 2 values from the linear model.
Univariate statistics in all cases were adjusted with Bonferroni correction. Significance was assumed at p b 0.05. Tables 7 and 8 present a summary of the statistically significant measurements in this study based on the univariate and multivariate statistics respectively.  Table 3 Area-independent curvature measures. Notation: H = mean curvature, G = Gaussian curvature, c = curvedness, A = surface area, T = 3*volume/A.  Tables presenting all the results of the analysis in this study can be found in the supplementary material.

Gompertz model
Gompertz functions are sigmoid functions often used to model growth, where growth is initially low, then accelerates to reach peak growth, and then decelerates to low growth again at the end (see Fig. 11). As in Wright et al. (2014), we examine the following Gompertz-like function to model the relation of brain descriptors with age: where β 1 is the initial value of f(t) and β 1 + β 2 the ending value of f(t). β 3 is the parameter that adjusts the growth rate and β 4 the t value where f(t) reaches peak growth. In this study we deviate from Wright et al. (2014) by allowing the parameters of the Gompertz function to take negative values. This allows us to further model decline of a descriptor in relation to age (see Fig. 11).

Volumetric measurements
Absolute and relative volumes of the brain tissues with increasing age at scan are illustrated in Figs The absolute volume of each of the 82 structures was significantly correlated with the age at scan. In addition, most of the structures have a significant linear correlation of their relative volumes to age at scan. These correlations were positive for relative volume of the CGM, CSF, cerebellum and corpus callosum and negative for relative volume of the ventricles, the majority of WM regions and the basal ganglia Fig. 11. Example of a positive Gompertz function (blue line) that models growth and a negative Gompertz function (red line) that models decline as a function of time t. The left plot displays the function f(t) and the right plot the gradient of the function df(t). The peak growth/decline occurs at time t = β t displayed with a dotted line. and thalami. The Gompertz model fitted more accurately the relation of both absolute and relative volumetric measurements with age at scan for the majority of regions. However, the reduction in the sum of squares error was smaller than 5% for most of the measurements. This suggests a significant linearity in the structure development. The structures that their volume was better fitted with the Gompertz function were: the anterior cingulate gyrus WM bilaterally, insula WM and CGM bilaterally, posterior medial and inferior temporal gyrus CGM bilaterally and the left posterior superior temporal gyrus WM, CGM. The WM and the WM parts reach a peak growth at around 33 weeks age at scan, while the CSF and CGM and its parcellations present peak growth later at around 38-39 weeks. The subcortical structures have the largest growth around the 31 weeks, with the exception of cerebellum and ventricles that present the largest increase in volume later at 37 and 40 weeks respectively. The total brain growth reaches peak growth at the age of 35 weeks. When entered in a multiple regression model with age at scan, the structures that had significant effect on the model were: brainstem, corpus callosum, cerebellum right, subthalamic nucleus left, the WM and CGM parts of the left medial anterior temporal lobe and left anterior fusiform gyrus, and the WM parts of right anterior cingulate gyrus, left posterior superior temporal gyrus, anterior gyri parahippocampalis bilaterally. The combined model was highly correlated with age at scan with adjusted R 2 = 0.935. Total brain volume of the preterm group was significantly smaller than the term controls. Reduction in total brain volume was significantly associated with increasing prematurity in the preterm population.
The preterm infants had significantly reduced WM volume and more specifically in the lateral anterior temporal lobe bilaterally, anterior medial and inferior temporal gyrus bilaterally and the left hemispheric parts of the parietal lobe, medial anterior temporal lobe, posterior cingulate gyrus, middle superior temporal gyrus and anterior gyri parahippocampalis. Increasing prematurity was associated with a reduction in total and regional WM volume in the majority of WM parts.
The CGM was less affected overall with significant group differences localised in the lateral anterior temporal lobe bilaterally and the left medial anterior temporal lobe. Decreasing GA at birth was significantly associated with reduced total CGM volume and was negatively correlated with CGM volume in the medial and lateral anterior temporal lobe bilaterally, medial and inferior temporal gyrus bilaterally, anterior gyri parahippocampalis bilaterally, anterior fusiform gyrus bilaterally left frontal lobe, left anterior cingulate gyrus and the right middle superior temporal gyrus.
The volume of subcortical structures was also affected by prematurity. The preterm subjects had reduced volume in the areas of the corpus callosum, subthalamic nucleus, left amygdala and right caudate nucleus. Volumes of all the subcortical structures were significantly reduced with increasing prematurity. The absolute volume of the left ventricle and relative volume of CSF and both ventricles were significantly increased for the preterm subjects. This increase of relative volume was additionally correlated with increasing prematurity. Prematurity was further associated with relative volume changes in multiple brain regions.
The tissues that were significantly associated with increasing prematurity according to multiple linear regressions were the WM, brainstem, CSF and ventricles. The multiple regression model including all the 82 structures resulted in the following significant predictors of increasing prematurity: the WM parts of left medial anterior temporal lobe, right lateral anterior temporal lobe, left posterior cingulate gyrus, left anterior and right posterior medial and inferior temporal gyrus, left anterior fusiform gyrus and the anterior gyri parahippocampalis, the CGM parts of left frontal lobe, right occipital lobe, right anterior cingulate gyrus, right middle superior temporal gyrus, left anterior and posterior gyri parahippocampalis, and the brainstem, corpus callosum and left lateral ventricle.
Accounting for the repeated measurements with the linear mixed model yields statistically different results from the linear model as assessed by the likelihood ratio. However, the obtained R 2 values are very similar to the linear model's R 2 values. The vast majority of R 2 values (98%) computed to quantify brain development and effect of prematurity had difference between the 2 models less than 0.02, which can be considered negligible.

Cortical surface measurements
The different surface measures of the cortex with respect to age at scan are illustrated in Fig. 14 and correlations in Table 4. An initial experiment was performed to look into the effect of the Laplacian smoothing of the cortical surfaces. Cortical surfaces reconstructed with and without Laplacian smoothing present very similar results in correlation with ages at birth and scan. When compared with an unpaired t-test, the mean curvature L 2 norm and global curvedness were significantly different between the two reconstructions, while the surface area and

Medial and inferior temporal g A l gm
Medial and inferior temporal g P r gm Medial and inferior temporal g P l gm Gaussian curvature L 2 norm did not reach significance. In the following, we present the analysis based on the smoothed surfaces. The curvature measures and surface area were positively related to the PMA at scan for the whole cortex and for almost all the cortical regions after Bonferroni correction (with the exception of the mean curvature L 2 norm and global curvedness of the left anterior gyri parahippocampalis). The regions whose surface area was significantly related to age at scan in a multivariate model were: right parietal lobe, left occipital lobe, right anterior cingulate gyrus, middle superior temporal gyrus bilaterally, anterior medial and inferior temporal gyrus bilaterally, left anterior and right posterior fusiform gyrus. Similarly for the curvature measures, the common regions that presented significant associations were: left medial anterior temporal lobe, left posterior superior temporal gyrus, anterior medial and inferior temporal gyrus bilaterally and left anterior gyri parahippocampalis. The adjusted R 2 values of the multivariate model for surface area, mean curvature L 2 norm, Gaussian curvature L 2 norm and global curvedness were respectively: R 2 = 0.896, R 2 = 0.921, R 2 = 0.906, R 2 = 0.916. The relative surface area, the surface area of a region normalized to the total surface area, presents both regional increases and decreases with increasing age at scan that are significantly associated in the majority of cortical regions.
The Gompertz model provided a better approximation to the surface measures with respect to age at scan than the linear model in the majority of regions. However, in all the measures apart from the global curvedness the improvement in the SSE for most regions was not considerable (less than 5%). Exception is the surface area of the anterior gyri parahippocampalis and the curvatures of the posterior parts of the medial and inferior temporal gyrus and fusiform gyrus, and the medial and lateral parts of the anterior temporal lobe. The global curvedness of more than half the cortical regions presented an improvement in SSE larger than 5% with decreases up to 31%. The structures that presented the largest gain from the Gompertz fitting were the posterior part of the medial and inferior temporal gyrus and the frontal, parietal and occipital lobe. The increase in surface area is maximised at the ages (at scan) of 34-38 weeks. Most of the regional curvature measures present peak growth at either the earliest or latest ages at scan.
The cortical surface area was found to be significantly reduced in the preterm subjects and more specifically in the lateral anterior temporal lobe bilaterally and the left hemispheric parts of the medial anterior temporal lobe and anterior gyri parahippocampalis. Increasing prematurity was further associated with decreasing surface area in the whole cortex and most of the regions. From a multiple linear regression model, the structures that were significantly related to the age at scan were: right frontal lobe, right parietal lobe, right occipital lobe, left medial anterior temporal lobe, right anterior temporal lobe, right middle superior temporal gyrus, left anterior medial and inferior temporal gyrus, posterior medial and inferior temporal bilaterally and left anterior gyri parahippocampalis.
The curvature measurements were not associated with age at birth in the whole cortex and the majority of cortical regions. A  ✓ Cingulate g A l wm ✓ Cingulate g P r wm Cingulate g P l wm ✓ Superior temporal g middle r wm Superior temporal g middle l wm Superior temporal g P r wm Superior temporal g P l wm ✓ Medial and inferior temporal g A r wm Medial and inferior temporal g A l wm ✓ ✓ (continued on next page) notable exception is the anterior part of the temporal lobe that consistently presented a positive correlation with increasing prematurity in all of the curvature measures. The preterm infants further demonstrated increased curvature measurements compared to the term controls, especially in the mean curvature L 2 norm and global curvedness. The right parts of the parietal lobe, posterior medial and inferior temporal gyrus, posterior superior temporal gyrus were significantly higher in the preterm subjects in all the curvature measures. These structures additionally presented increased relative volume in the preterm subjects and this increase was correlated with increasing prematurity. Increasing prematurity was significantly associated according to a multiple regression model with: the right occipital lobe, right anterior medial and inferior temporal gyrus, left anterior and right posterior gyri parahippocampalis and left anterior and right posterior fusiform gyrus. Similar to the volumetric results, R 2 values computed with the linear mixed model were essentially the same as that of the linear model with a maximum difference of 0.014 across all the surface measurements with respect to age at scan and age at birth.

Comparison with manual measurements
The tissue volumes obtained after the proposed corrections are similar to volumes in the literature evaluated using manual segmentation approaches. Anbeek et al. (2008) provide average tissue volumes (mL) of 13 subjects around term who were born over a wide age range of gestations (gestational age 25.9-42.9 weeks, corrected age at test -3.6-5.1 weeks): CSF 51.4, CGM 101.2, WM 146.4, BGT 20, Brain 319. Corresponding CGM and WM volumes are obtained here over the scan ages of 36-40 weeks (see Fig. B.3 of Appendix B). The CGM in Anbeek et al. (2008) represents about 32% of the brain volume and the WM around 46% of the brain volume. The relative volumes obtained here are 34% for the CGM and 48% for the WM around term (see Table 5). It should be noted that the relative volumes prior to correction for the CGM-WM partial volume and sulci correction were 43% and 40% for the CGM and WM respectively. This overestimation of CGM obtained prior to the corrections is consistent with previous automatic segmentation studies (see Table 6).
Similar volumes to manual results are further obtained for the early preterm period after the proposed corrections (see Table 5). Moeskops et al. (2013) present a relative CGM volume of 18% and relative WM volume of 70% for 10 neonates scanned at 30.8 ± 0.7 weeks age at scan. The relative volumes in our study are 25% for the CGM and 57% for the WM around 30 weeks age at scan (prior to the corrections the relative CGM and WM were 31% and 52% respectively). The CGM oversegmentation caused by a Gaussian Mixture Model that assumes one class for WM and one class for CGM can be observed in Fig. 15. Without the introduction of the CGM-WM partial volume correction, the segmentation tends to attribute a larger proportion of the brain to the CGM.
Median thickness across the subjects in the cohort is presented in Fig. 15. The cortical thickness estimated using the cortical segmentations without the sulci correction produces an increasing thickness with age at scan. The thickness of the uncorrected segmentations correlates significantly with the age at scan (p b 10 −36 ). However, with the introduction of the sulci correction, the cortical thickness measured over the subjects remains unaffected by the age at scan (p = 0.07) of the neonate. The cortical thickness of the neonatal brain has a median value of 1.59 ± 0.09 mm across the database (the 25th and 75th percentiles are 1.54 and 1.65 mm respectively). A Gompertz fit to the thickness data decreases the mean squared error by a factor of 2% and displays a gradual increase of thickness after the age of 38 weeks from 1.59 to 1.65 mm.

Discussion
Quantitative measurements of the developing neonatal brain are required to study normal brain growth and potentially aid the early diagnosis of later neurological impairments. In this study, we employed the automatic segmentation algorithm proposed in Makropoulos et al. (2014) to delineate 82 regions of the brain in a large group of infants and derive a number of volumetric and cortical surface measurements. Two corrections are incorporated for the detailed delineation of the cortical ribbon in the neonatal brain. The first correction estimates a partial volume class between the CGM and WM which is consequently relabelled as WM in order to limit the over-inclusion of voxels in the CGM tissue. The second correction detects and delineates the cortical sulci that are hard to segment with intensity-based segmentation techniques. We initially detect the cortical sulci from the expansion of the interior cortical surface as areas of the surface that collapse to each other similar to Han et al. (2004). The thickness of the detected sulcal areas is then approximated from neighbouring parts of the cortical ribbon where the thickness can be accurately measured. It should be noted that the sulci correction proposed assumes a similar cortical thickness in the local neighbourhood of the cortical ribbon. This could potentially introduce errors in regions where this assumption does not hold. However since the cortical thickness is typically less than 2 voxels, we would expect only negligible errors from this assumption. Derived volumetric and thickness measures after the application of the method presented here are similar to measurements obtained from manually segmented data. A structural atlas is constructed for different ages of the neonatal brain for all the segmented brain structures and is made publicly available. The atlas defines the structure probability and average segmentation respectively of each structure in the spatio-temporal space of Serag et al. (2012).
Cortical surface area measurements have been previously presented for the neonatal brain with a range of 150-1500 cm 2 between 27 and 44 weeks PMA at scan (Kapellou et al., 2006;Xue et al., 2007;Dubois et al., 2008b;Pienaar et al., 2008;Rodriguez-Carranza et al., 2008;Moeskops et al., 2013Moeskops et al., , 2015. The surface area of the cortex in this study was around 120-1100 cm 2 in the corresponding ages at scan. Curvature measurements have been reported in a limited number of Medial and inferior temporal g P r wm ✓ Medial and inferior temporal g P l wm Gyri parahippocampalis A r wm ✓ ✓ Gyri parahippocampalis A l wm ✓ ✓✓ Gyri parahippocampalis P r wm Gyri parahippocampalis P l wm Fusiform g A r wm Fusiform g A l wm ✓ ✓ ✓ Fusiform g P r wm Fusiform g P l wm studies (Xue et al., 2007;Pienaar et al., 2008;Rodriguez-Carranza et al., 2008;Moeskops et al., 2013;Moeskops et al., 2015) which used different definitions of curvature measures and included only small numbers of subjects. Here, the curvature measures from Rodriguez-Carranza et al. (2008) were adopted that are invariant to the surface area. Similar positive correlations of cortical curvature with age at scan to Rodriguez-Carranza et al. (2008) are derived in this analysis. Xue et al. (2007), Dubois et al. (2008b), Moeskops et al. (2015) likewise presented increasing mean curvature and gyrification with increasing age at scan. We should note here that the cortical surfaces were smoothed prior to obtaining the surface measurements which might influence the results. Laplacian smoothing did not have a major impact on consequent measurements. The cortical surface was reconstructed after blurring the cortical mask to avoid a blocky appearance due to the limited resolution, which would lead to extreme curvature values. Although this might affect the measurements, we do not expect a significant effect due to the small amount of blurring introduced (Gaussian kernel of 1 mm standard deviation) and due to the fact that median values are only reported in this study.
Cortical thickness measurements in the neonatal population have been previously presented in (Xue et al., 2007;Moeskops et al., 2013;Moeskops et al., 2015) for limited datasets. These studies obtain a median cortical thickness of around 1-1.4 mm for the neonatal brain. Moeskops et al. (2015) further demonstrate an increase in cortical thickness from around 1 mm to 1.4 mm between 30 and 40 weeks age at scan in the preterm brain. The cortical thickness estimated here remains almost constant in the neonatal brain, similar to Xue et al. (2007), with a value around 1.6 mm, from the early preterm period to term-equivalent age. Differences in thickness values can be attributed to the different in-plane resolution of the MRI (Moeskops et al. (2013 have a highly anisotropic resolution, 0.34 × 0.34 × 2 and 0.35 × 0.35 × 1.2 mm, while the analysed data have a near isotropic resolution of 0.86 × 0.86 × 1 mm), and different thickness measurement methods (Xue et al., 2007).
We also obtained absolute and relative volumes of brain tissues from the early preterm period to term-equivalent age. Surface area and curvature measures of the whole cortex and regional cortical parts were estimated based on the segmentations. Our results show that, with the exception of cortical thickness, regional brain and cortical growth is significantly associated with brain maturation. A Gompertz function presents a better approximation than a linear model for the relation of the volumetric and surface measurements with age at scan. This is expected given the increased degrees of freedom of the Gompertz function (4 variables to adjust instead of 2 for the linear model). However, the reduction in the sum of squares error with respect to the linear model is less than 5% in the majority of measurements, suggesting that the linear model can still capture the relationship with age at scan. A similar relationship of cortical folding with age is exhibited in the fetal brain (Wright et al., 2014) over the ages of 27-39 weeks. However, the fetal data in Wright et al. (2014) additionally cover an earlier spectrum of ages starting from 22 weeks gestational age. Over the age range of 22-27 weeks the fetal brain presents large, non-linear, increases in cortical folding with peak growth at 30-32 weeks. Consequently, the Gompertz model in the fetal data provides a better approximation with a fitting error approximately half of the linear model (Wright et al., 2014). A recent study comparing cortical folding between preterm newborns and fetuses (Lefevre et al., 2015) presented similar results. The cortical curvature in both the fetal and neonatal brain demonstrated a linear relationship with age over the range of around 28-36 weeks, although of different extent in the two groups. The fetal brains prior to that age range demonstrated a non-linear increase with age, although this effect was not specifically explored (Lefevre et al., 2015). The neonatal brain presents peak volumetric growth at 35 weeks age at scan which is presented progressively at the subcortical regions except for the cerebellum, the WM regions, the cerebellum, the CSF, the CGM parts and the ventricles. The growth of the regional surface area is maximised at the ages of 34-38 weeks. Multivariate analysis of the volumetric and surface measures including all the regions result in adjusted R 2 values ranging from 0.896 to 0.935.
We used the volumetric, surface area and curvature measures to characterise the effect of prematurity in the neonatal brain, and to compare preterm brain development with that of healthy term born controls. Total brain volume in the preterm infants was reduced compared to term controls. CLD has been previously associated with reduced total brain volume compared to preterm infants without CLD who do not seem to have reduced brain volumes . Almost half the infants in the preterm group (17 out of 40) had CLD and this may be reflected in the reduced brain volume in this study. Total brain volume was also highly negatively correlated to increasing prematurity, as has been reported previously (Peterson et al., 2003;Inder et al., 2005;Thompson et al., 2007;Ball et al., 2012). Prematurity is related to reductions in the volume of the WM, CGM and subcortical structures, and increases in the relative CSF and ventricular volume and these alterations are significantly associated with decreasing age at birth. Regional decreases in the WM have been previously described in Mewes et al. (2006), and Thompson et al. (2007). These studies presented both reductions and increases in the regional CGM volumes. Similar volumetric associations in the subcortical GM and more specifically in the amygdala, thalamus, hippocampus and lentiform nucleus have been reported in Peterson et al. (2000), Srinivasan et al. (2007), and Ball et al. (2012). Larger volumes of CSF and ventricles in the preterm subjects have been found in previous studies (Peterson et al., 2000;Mewes et al., 2006;Thompson et al., 2007).
Preterm infants had reduced surface area compared to term controls in only a few parts of the cortex. Surface area reductions were however significantly correlated with increasing prematurity in the majority of the cortical regions. Ajayi-Obe et al. (2000) similarly presented reduced cortical surface area in preterm infants compared to term controls. Kapellou et al. (2006) further demonstrated a decreasing surface area in the cortex with increasing prematurity. Cortical curvature was largely not associated with the age at birth of the infants. An exception is the anterior temporal lobe that presents a positive correlation with increasing prematurity. Kesler et al. (2006) demonstrated similar results in prematurely-born children, where the temporal lobe was shown to be specifically disrupted by preterm delivery with increased gyrification in the preterm population. Kesler et al. (2006) suggested that increased gyrification may be due to abnormal growth of the inner cortical layers. Here, these alterations are specifically localised in the anterior part of the temporal lobe. The preterm infants further demonstrated increased curvature measurements compared to the term controls. Specifically, the curvature of the right parts of the parietal lobe, posterior medial and inferior temporal gyrus and posterior superior temporal gyrus was significantly higher in the preterm subjects for all the curvature measures. The relative volume of these structures was also significantly increased compared to term controls and was correlated with increasing prematurity. Future studies with the inclusion of clinical variables and neurodevelopmental outcome will help to further elucidate the effect of prematurity in the neonatal brain. Appendix A. Sulci correction 'Shock' points are detected as in Han et al. (2004): The CGM-WM interface is propagated with the fast marching method. The distance D FMM to the CGM-WM interface is defined by solving the following Eikonal equation: where i indexes the voxels of the image. D FMM (i) = 0 for voxels in the CGM-WM interface and F(i) = 1 − 0.9p i,CSF according to the CSF posterior p CSF of the EM algorithm. The gradient ∇D FMM can subsequently be used to identify 'shock' points where the spatial derivative is not well defined: Here T is set to 0.8 as in Han et al. (2004). At these points two sulcal banks will merge in the propagated surface and ∇D FMM will be small.
Sulcal points are then selected among the 'shock' points based on the assumption that the thickness of the cortical ribbon is locally consistent. A distance D allowed,i is measured for each 'shock' voxel as the locally weighted mean cortical thickness measured on other points of the cortical ribbon. D allowed,i is only averaged over points with streamlines that do not cross a 'shock' point, in essence points that are not inside a cortical sulcus. Connected components of 'shock' points, points that belong to the same sulcus, are identified with connected component labelling. Points that belong to the same sulcus are then defined to have the same D allowed , which is the mean D allowed,i of the 'shock' points in the component. Finally, 'shock' points are labelled as CSF only if their distance to the CGM-WM interface D WM,i is larger than D allowed . Relative volume (ratio of structure's volume to the total brain volume) centiles of tissues with increasing scan age. The yellow line for each region represents the 25% centile over the subjects in the database, the red line the median value, and the blue line the 75% centile.