Neurocranium thickness mapping in early childhood

The neurocranium changes rapidly in early childhood to accommodate the growing brain. Developmental disorders and environmental factors such as sleep position may lead to abnormal neurocranial maturation. Therefore, it is important to understand how this structure develops, in order to provide a baseline for early detection of anomalies. However, its anatomy has not yet been well studied in early childhood due to the lack of available imaging databases. In hospitals, CT is typically used to image the neurocranium when a pathology is suspected, but the presence of ionizing radiation makes it harder to construct databases of healthy subjects. In this study, instead, we use a dataset of MRI data from healthy normal children in the age range of 6 months to 36 months to study the development of the neurocranium. After extracting its outline from the MRI data, we used a conformal geometry-based analysis pipeline to detect local thickness growth throughout this age span. These changes will help us understand cranial bone development with respect to the brain, as well as detect abnormal variations, which will in turn inform better treatment strategies for implicated disorders.

Finally, neurocranial fractures, the leading symptom of injury in infants due to motor accidents, falls and abuse 15 may also lead to abnormalities. A growth map of the neurocranium may help in addressing correct sutures and bone plates positions in these cases. A previous study has shown that while adult cranial sutures and cranial bones have a similar stiffness, pediatric cranial sutures can deform 30 times more than cranial bone before failing and 243 times more than the adult cranial bone 16 . Due to this, fracture may cause the neurocranium to sustain a major shape change, impacting the brain as well. It is thus imperative for clinicians to accurately diagnose neurocranial deformities early, and to accurately monitor progress and development, both on its own and with respect to the brain. The lack of normative data and quantitative imaging tools for these analyses are major obstacles to the early detection and optimal treatment of conditions of infancy. Additionally, treatment approaches and their success vary widely among institutions and clinicians and are usually based on subjective criteria.
Computed tomography (CT) images are typically used to diagnose abnormalities of the neurocranium. While the use of this modality is often necessary when diagnosing pediatric patients, studies have also shown a small but significant increase in the risk of long term adverse outcomes such as developing brain cancer or leukemia 17,18,21 , which suggests that alternative imaging procedures without ionizing radiation should be used whenever possible 18 . It is, therefore, inappropriate to subject healthy children to CT scans as a way to collect data for growth maps of shape changes of the neurocranium. Fortunately, this structure can also be delineated on magnetic resonance images (MRI). An example of such a segmentation for an 18-month-old brain is provided in Fig. 1.
MRI has also recently shown value in diagnosis of craniosynostosis through a novel sequence that enhances the bone-soft tissue boundary 19 , as well as to investigate craniosynostosis prenatally 20 .
In the current study, we map normative neurocranial thickness in children between 6 and 36 months of age using a large, preexisting database of MRI scans of healthy children ranging in age from 0 to 3 years old, acquired at the Baby Imaging Lab (https ://www.babyi magin glab.com). We focus on thickness because it is an important factor in predicting susceptibility to injury, as well in diseases of the neurocranium such as thalassemia and Paget's disease. Thickness measurements also have potential for applications in electroencephalogram measurements, that can be erroneous due to variation in neurocranium and scalp thickness 21,22 . In this paper, we use neurocranium masks derived from high resolution MRI brain scans in young children and introduce a new pipeline to investigate neurocranial thickness on MRI data. We first generate tetrahedral meshes, allowing subvoxel resolution analysis of the thin neurocranium surfaces. We then solve the Laplace equation at each vertex on the neurocranium surface, creating a harmonic field through which we calculate thickness using streamlines. The weighted-SPHARM approach is used to register surfaces, allowing group comparisons. The general flow of processing is shown in Fig. 2 23 . This paper adds to the body of studies involving CT data in ages under 36 months to map neurocranial thickness.

Method.
Data. All the methods were performed in accordance with relevant guidelines and regulations.
The dataset used includes 56 T1 MP-RAGE MRIs (1.4-1.8 mm 3 ) of healthy normal children between the ages of 6-36 months of age at the time of scanning. The number of subjects in each age group is shown in Table 1. The inclusion criteria is as follows: singleton, full term (37-42 weeks at birth) with no abnormalities on fetal ultrasound and no reported history of neurological events or disorders in the infant. Data acquisition details can be found in 24,25 . The study was approved by the Institutional Review Board of Brown University and informed consent was obtained from the guardians of all participants. All data was de-identified before pre-processing and acquired on a Siemens 3 T Tim Trio scanner equipped with a 12 channel head RF array. To minimize intra-scan motion, children were asleep and swaddled with a pediatric MedVac vacuum immobilization bag (CFI Medical Solutions, USA) and foam cushions. Scanner noise was reduced by lessening the peak gradient amplitudes and slew-rates, and using a noise-insulating scanner bore insert (Quiet Barrier HD Composite, UltraBarrier, USA). www.nature.com/scientificreports/ Pre-processing. The data was pre-processed as follows. First, the MRI brain volume was skull-stripped using FSL BET 26,27 , and resampled to a 1 × 1 × 1 mm 3 resolution for consistency throughout processing. Bias correction was performed using the N4 ANTs bias correction tool 28 . The result was then linearly registered to an agematched custom template, using FSL FLIRT 29,30 with 6 degrees of freedom. This custom template was previously generated from the same dataset as described by Remer et al. 31 , and also resampled to 1 × 1 × 1 mm 3 resolution. The transformation that resulted from the registration was saved and then applied to the original dataset with the neurocranium. The resulting registered image (with the neurocranium) was then bias field corrected using the N4 ANTS bias correction tool. Finally, we applied the FSL BET 26,27 tool to extract the outline of the inner and outer neurocranial surfaces (see example in Fig. 3a). The extracted structure does not include the craniofacial skeleton. Each mask was visually inspected to ensure accuracy, and to eliminate datasets with overlapping voxels between the two surfaces, or those which included areas that were not part of the neurocranium. Note that our very first step involved using skull-stripped images. This was done as we needed to register the data for group comparison. Without removing the neurocranium initially, the high intensity signal from the scalp caused the registration to be imprecise.
Thickness calculation. The neurocranium thickness calculation used in this paper involved a similar approach to that used for the corpus callosum in previous studies 32,33 . We used the Iso2mesh toolbox 34 to create volumetric tetrahedral meshes for the inner and outer neurocranium masks generated by FSL BET 27 . An example is shown in Fig. 3b and c. In short, the thickness is computed from the inner boundary (M0; internal patch) to the vertices on the outer surface mesh (M1; external patch) of a mesh M, by optimizing a string energy function (E):  We are looking for the piece-wise linear function f : V → R that minimizes E. f satisfies Laplace's equation with Dirichlet boundary conditions. We use the volumetric Laplace-Beltrami operator on the tetrahedral mesh 31 , and fix the values of the function f on the external boundary (B1) as 1, the internal boundary (B0) as 0: If we define u(s) is a parametric curve with arc length s and x is a point on the surface patch we start from, streamlines 36,37 are constructed to connect the two surfaces by solving the equation below.
The thickness is defined as the total arc length of the streamline that traverses the neurocranium from external to internal patches 32,35,36 .
Surface registration. Next, we used the weighted-SPHARM approach described in 38,39 to register each individual subject to a common template. This enabled us to get a surface with the same number of vertices for each subject, allowing us to compare thickness locally between groups. Each neurocranium's internal and external surfaces were mapped to a sphere using an area preserving and surface flattening algorithm 38,39 . We utilized a template sphere with 40,962 vertices to compute the spherical harmonic components in order to get a SPHARM representation of the surfaces. We utilized non-parametric t-tests with 10,000 random assignments of subjects to groups to compare consecutive age groups. Finally, we conducted a regression analysis to understand neurocranial thickness developmental trajectories. Due to the small number of 6 month old subjects, we did not use it for group comparison. However, this dataset was included in the regression analysis since age in days was used and there was no need to cluster data by age bracket as with the group comparisons. Figure 4 shows the neurocranial thickness at a.12 months, b. 24 months and c. 36 months. The color bar goes from dark blue to red, with the latter indicating the thicker regions. Our results show consistent increase in the thickness of the neurocranium, around the midline of the occipital bone under the posterior fontanelle, indicated with the arrows. Internal regions (z-value < 95), including those closer to the air cavities and areas closer or below the cerebellum, were disregarded in this study due to variability in the data, and were set to 0 in the results below.

Results
Results in Figs Fig. 6 shows the difference between the neurocranium thickness at 12 months and 36 months. We can clearly see that the differences appear to correspond to the occipital bone under the posterior fontanelle.   Figure 7a1 and a2 illustrate the regression coefficients and the p-values of the neurocranium thickness change with respect to the child's age (in days) from 6 to 36 months. Here, again our results show that the largest change is localized to areas corresponding to the occipital bone of the neurocranium. Using the most significant vertex, we plotted the progression of neurocranium thickness change age in children between the ages of 6 months and 24 months in Fig. 7b. The 36 month data was excluded for this figure as the t-tests showed little difference in neurocranium thickness between the age groups of 24 and 36 months.

Discussion
Our analysis shows that the biggest change in thickness observed in the neurocranium in early childhood is around the midline of the occipital bone under the posterior fontanelle. Additionally the largest increase in this area is between the ages of 12 m and 18 m and fewer changes are observed thereafter, indicating major development in the posterior part of the head happens in the first 1.5 years of life.
While most studies focus on global measures such as neurocranial volume or head circumference, our study using non-invasive MRI adds to the body of previous studies by providing information on the local thickness changes of the developing neurocranium. Our results give us insight into the development of the neurocranium, as maps shows a clear increase in neurocranial thickness in the posterior part of the head, likely overlapping with the posterior fontanelle.
In a previous study 4 , authors developed a statistical head model for children aged 0-3 from CT images, allowing them to predict neurocranium geometry. Parameters such as age and head circumference were used to predict features such as neurocranium size/shape, suture geometry and neurocranium thickness. The authors showed a greater thickness of the neurocranium in the occipital region in ages 1.5 years and older. In addition, the model also showed that the lambdoid suture closes between the ages of 1 and 1.5 and neurocranium thickness increased in regions corresponding fontanelle after it closes. In Fig. 5b, the most rapid increase in thickness was near the posterior fontanelle and along the lambdoid suture between the ages of 12 and 18 months. Hence, the thickness differences we found in Fig. 5b, agree with this neurocranium growth model 4 . Additionally, we also see a thickness increase in the parietal bones surrounding the sagittal suture between these age groups, consistent with marginal increase in thickness shown in the CT data generated statistical model 4 . Figure 6 shows the thickness difference between 12 and 36 months; the biggest difference around the posterior fontanelle and regions corresponding to the lambdoid sutures. Again, this may be attributed to neurocranium growth after posterior fontanelle fuses.
The largest difference between the 18 and 24 months old cohorts (Fig. 5c) was seen in the anterior part of the head and in the bones surrounding the occipital and temporal lobes. The concentration of points in the temporal region corresponds to a region in between the parietal and temporal bones. Similar differences were also seen when comparing the 24 and 36 month groups (Fig. 5d). Figure 5 shows that the thickness increases non-uniformly, showing larger changes in the posterior part of the neurocranium earlier compared to the anterior region. This non-uniform pattern may also be associated to the posterior to frontal development of the brain itself, as shown in several studies on brain development 40,41 . Additionally, it has also been observed that the However, it must be noted that as the differences in Fig. 5c and d are not seen bilaterally, results may also be influenced by the variability of the data or inaccuracies that arise during registration. Alternatively, the patterns may indicate that while most of the neurocranium thickness changes in the posterior and occipital bone area happen prior to 2 years, those of the frontal and temporal bones may continue to evolve even after age 2.
Normal healthy development of the neurocranium and its relationship to brain development is also relevant in understanding deformational plagiocephaly, which occurs due to the continuous placement of an infant in one position the supine position only. Studies have shown that children with deformational plagiocephaly have motor delays in early childhood 1 . Additionally, some children may be predisposed to deformational posterior plagiocephaly due to localized cranial flattening, even though this happens in 13% of healthy newborns. A recent study investigating the association of deformational posterior plagiocephaly and visual field abnormalities showed that in 35% of the 40 infants studied, one or both hemifields were constricted 43 . While these are not causal relationships, and it is possible that deformational plagiocephaly occurs in children who have developmental delays, there are several studies implying that these children face a higher risk of developmental difficulties, language disorders and attention span problems 44,45 . Therefore, it will be beneficial for physicians assessing these children to have a growth model of the normally developing neurocranium to compare to as well as to recognize the relationship between neurocranium and brain development in early childhood to prevent developmental delays from persisting. Additionally, such a growth model may also be beneficial to assess the progress of diseases such as thalassemia, Paget's disease, to monitor intracranial pressure, as well as use in electroencephalogram measurements 21,22 .
The CT data generated statistical neurocranium models 4 also show changes in the anterior part of the neurocranium, near the coronal suture, that we were not able to observe in our group comparisons. This may be due to the fact that this change was smaller compared to the posterior part, and we were not able to observe it due to our smaller sample size. This model also shows a small area of the sagittal suture closed after age 2. Thickness www.nature.com/scientificreports/ may increase gradually after the suture closes and therefore, our dataset may not cover the age range to observe a thickness change that occurs after the suture closes. The data size used in this study and the age range is a limitation as we could reliably extract the neurocranium only from a small number of subjects with the highest quality images, due to the fact that the neurocranium is quite thin in early childhood and difficult to automatically extract from MRI data. To overcome this problem, future studies investigating the neurocranium may benefit from new MRI sequences to enhance the bone-soft tissue boundary 19 .

conclusion
In this study, we investigated the growth of the neurocranium across ages 6 to 36 months. Our results showed the biggest change in the posterior part of the neurocranium, where the lambdoid suture lies, between 12 and 18 months, corroborating results of the geometric head model created from CT data 4 . Therefore, it is important in future studies to investigate younger age groups as well as age groups in smaller intervals to build a more detailed timeline of neurocranium growth. The neurocranium growth model will improve over time by inclusion of more subjects to become a reliable template of neurocranium growth. We show in this pilot study, that using MRI head scans, neurocranium thickness growth can be mapped out to show changes in areas that correspond to previous CT studies. Future work will include a correlation of the neurocranium thickness and parameters acquired from the brain, such as cortical thickness, which will allow us to acquire a more complete picture of the development of the brain and neurocranium. The association between plagiocephaly and neurocranium development, particularly in light in of the fact that our results show consistent group differences in the back of the neurocranium, will also be investigated.