Data-driven Normative Reference of Pediatric Cranial Bone Development.

Available normative references of cranial bone development and suture fusion are incomplete or based on simplified assumptions due to the lack of large datasets. We present a fully data-driven normative model that represents the age- and sex-specific variability of bone shape, thickness, and density between birth and 10 years of age at every location of the calvaria.
Methods: The model was built using a cross-sectional and multi-institutional pediatric computed tomography image dataset with 2068 subjects without cranial pathology (age 0-10 years). We combined principal component analysis and temporal regression to build a statistical model of cranial bone development at every location of the calvaria. We studied the influences of sex on cranial bone growth, and our bone density model allowed quantifying for the first time suture fusion as a continuous temporal process. We evaluated the predictive accuracy of our model using an independent longitudinal image dataset of 51 subjects.
Results: Our model achieved temporal predictive errors of 2.98 ± 0.69 mm, 0.27 ± 0.29 mm, and 76.72 ± 91.50 HU in cranial bone shape, thickness, and mineral density changes, respectively. Significant sex differences were found in intracranial volume and bone surface areas (P < 0.01). No significant differences were found in cephalic index, bone thickness, mineral density, or suture fusion. Conclusions: We presented the first pediatric age- and sex-specific statistical reference for local cranial bone shape, thickness, and mineral density changes. We showed its predictive accuracy using an independent longitudinal dataset, we studied developmental differences associated with sex, and we quantified suture fusion as a continuous process.

Takeaways

Question: What are the normative local distributions of cranial shape, thickness, and mineral density, and how are they modulated by age and sex in children?

Findings: We present the first data-driven normative reference of local cranial bone development between birth and 10 years built using a CT image dataset of 2068 subjects. We studied the local sex-specific continuous evolution of cranial shape, thickness, and bone and suture mineral density.

Meaning: Our data-driven model is a necessary tool to accurately identify, quantify, and study the temporal evolution of cranial development anomalies to objectify and personalize diagnoses and treatments.

INTRODUCTION

The calvaria is the upper part of the cranium that encloses the brain and is primarily formed by five major bones: left and right frontal, left and right parietal, and occipital bones. These bones are separated by the sutures, which allow the bones to rapidly expand and accommodate the space needed for a fast-growing brain. Since cranial bone growth is primarily induced by the developing brain, brain growth anomalies usually translate into cranial bone abnormalities and vice versa. Hence, normative quantitative models of local cranial bone anomalies are essential to study not only the phenotypes of cranial pathology but also how they may affect brain development.

The traditional clinical assessment of cranial bone anomalies relies on subjective evaluation because of the lack of quantitative normative references of cranial development. Simple clinical metrics of shape and size, such as cephalic index or head circumference, are used clinically even though they have shown low sensitivity identifying pathology. Therefore, a high variability in incidence reports of pathology, treatment approaches, and outcomes have been reported. Moreover, there are no widespread clinical metrics of essential aspects of cranial development such as bone thickness or bone and suture mineral density. Although longitudinal computed tomography (CT) images could be used to quantify cranial bone changes during childhood, they are not widely available because they require the use of harmful radiation and sedation in young children.

State of the Art

Average normative references have been used to guide cranial bone surgical reconstruction but they were not personalized to the specific anatomy of a patient. Our team created a statistical cranial shape model in the pediatric population using CT images to quantify shape malformation in individual patients. We also used similar models to calculate the optimal personalized surgical plan and to evaluate postsurgical head shapes. However, these models did not consider the normal anatomical changes that occur as children grow to evaluate and predict development. Partially overcoming this limitation, an age-dependent statistical model of cranial bone shape and thickness was built for children between 0 and 3 years old. However, this model built from a small dataset (N = 56) did not account for sex differences and only modeled a simplified representation of the cranium using radial basis functions. With a more detailed cranial representation, another age-dependent statistical model of cranial bone shape and thickness was built for children between 3 and 10 years old. However, the dataset was limited (N = 42), and the model assumed linear temporal changes, did not address sex differences, and did not model growth during the first three years of life when most cranial growth occurs.

Because of the lack of large pediatric datasets, other works have used simulation models to predict temporal cranial shape changes. Although finite-element methods have been extensively used to simulate cranial growth for surgical planning, their clinical translation has been hindered by the simplified assumptions about the biophysical processes driving cranial bone development.

In our previous work, we showed that data-driven cranial shape models can identify common temporal cranial growth patterns and enable a personalized prediction of development. However, we only modeled shape changes during the first two years of life; our method was not computationally feasible for large datasets at high temporal resolutions; and our approach did not consider sex, bone thickness, or mineral density.

Indeed, a few works have quantified cranial bone thickness and mineral density changes. Although thickness was modeled in, it suffered from limited datasets used for their constructions, and it assumed unrealistic linear temporal changes. Because of the absence of large CT image datasets, we also built a dense cranial bone thickness model during early childhood using magnetic resonance images. However, accurate bone segmentation from magnetic resonance images is still an ongoing area of research. To our knowledge, there are no quantitative data-driven models of cranial bone mineral density.

We present a fully data-driven normative statistical model of pediatric cranial bone development between birth and 10 years of age. Unlike other studies, our model was constructed without making assumptions about development; incorporates sex variability; and quantifies cranial bone shape, thickness, and mineral density at every location of the calvarial surface.

MATERIALS AND METHODS

Dataset

Retrospective cross-sectional CT images of 2068 subjects (1103 boys, 965 girls, age 3.12 ± 3.05 years, see Fig. 1 for age distribution) were used for model construction with IRB approval at University of Colorado Anschutz Medical Campus (protocol 20-1563) and Children’s National Hospital (protocol 3792). This dataset was screened to exclude patients with malformations, hydrocephalus, intracranial tumor, hemorrhage, hardware (eg, shunts), and prior history of surgery. Average in-plane resolution was 0.36 ± 0.06 mm with slice thickness of 1.97 ± 1.28 mm. In addition, longitudinal CT images of 51 additional independent subjects (28 boys, 23 girls, age at first scan 2.24 ± 2.22 years, age at second scan 3.55 ± 2.71 years, and time between scans 1.29 ± 1.71 years) were used to evaluate the predictive accuracy of the model.

Fig. 1.

Age distribution of the cross-sectional dataset, stratified by sex.

Automatic Image Processing

We used previous methods to automatically segment the cranial bones from the CT images. In summary, adaptive thresholding algorithms were used to extract the cranial bones (Fig. 2A). Image registration to a reference template was used to isolate the calvaria from the rest for the skull, and graph-cut-based segmentation was used to automatically label each cranial bone (Fig. 2B). All results were manually screened for errors in cranial base identification or bone labeling. After converting the segmented images to volumetric meshes using the marching cubes algorithm (Fig. 2C), we created a simplified representation of the calvaria using spherical mapping. Hence, three 2D spherical maps containing the shape (Euclidean coordinates in mm), bone thickness (in mm), and the Hounsfield units (HUs) from the CT images as surrogate of bone mineral density were generated, as shown in Figure 3. To establish local anatomical correspondences between the spherical maps of every subject guided by the location of the cranial sutures, we further aligned the spherical maps with the bone labels of all subjects using the diffeomorphic demons algorithm (Fig. 3).

Fig. 2.

CT image preprocessing pipeline. A, Cranium extracted from a CT image using adaptive thresholding. B, Cranial base landmarks represented as gray spheres at the glabella, clinoid processes of the dorsum sellae, and opisthion, together with the bone labels color coded on the cranium. C, Average Hounsfield unit values from the CT image color coded on the cranium.

Fig. 3.

Standard maps generated using spherical sampling. A, Thickness map (in mm). B, Image intensity map (in Hounsfiled units). C, Magnitude of the Euclidean coordinates map (in mm). D, Bone labels (LF for the left frontal bone, RF for the right frontal bone, LP for the left parietal bone, RP for the right parietal bone, and O for the occipital bone).

Model Construction

The aligned spherical maps were used to train a statistical model of cranial bone development. We used principal component analysis (PCA) to represent information in the spherical maps as

where represents all spherical maps in our dataset, the average spherical map, the principal components, and the PCA coefficients.

Given the orthogonality of the principal components, we used regression to model each PCA coefficient similar to the previous methods. We evaluated the performance of different nonlinear temporal regression functions summarized in Table 1‚ accounting for sex unlike other works. We estimated the average PCA coefficient fitting error for each regression function as

Table 1.

Regression Functions Evaluated to Model Temporal Development of the PCA Coefficients

where is the number of the principal components, is the number of subjects, is the PCA coefficient of components representing subject , represents age, represents sex (0 for girls and 1 for boys), is the regressed PCA coefficient, and represents the regression parameters. In addition, we calculated the fitting error in mm (for shape and thickness) and HUs (for bone mineral density) as

where represents the number of sampled data points in the calvaria, are the regressed PCA coefficients at the age and sex of every training subject, and is the Frobenius norm. We solved the following weighted mean-squared-error minimization problem to estimate the regression parameters :

The introduction of the weighting term allowed compensating for the uneven number of subjects at different ages, which was modeled as

The introduction of the weighting term allowed compensating for the uneven number of subjects at different ages, which was modeled as

where is a Gaussian kernel.

To study the age- and sex-specific variability of each coefficient representing the training subjects on each principal component , we created a regression model of their standard deviation as

where is the standard deviation regression function for the principal component and its parameters. Polynomial functions of degrees 1–3 were evaluated.

Normative Metrics

After model construction, we reconstructed the average cranial shape, thickness, and bone mineral density between birth and 10 years using Equation 1. In addition to quantifying the normative ranges for these three metrics at every location, we calculated the cephalic index, cranial bone surface areas, and intracranial volume. Since our methods quantify for the first time the variability in bone mineral density at every location of the cranial bones and sutures, we proposed a suture fusion index that quantifies the age- and sex-specific degree of fusion of every suture as

where represents every point within suture , represents the mineral density at , represents the minimum mineral density of suture , and represents the median mineral density of the bone plates surrounding suture at location , which we implemented by sampling bone mineral density values in the perpendicular direction of the suture at location . Note that is a normalization factor to keep the fusion index between 0 (open suture) and 1 (fused suture).

Temporal Predictions

We evaluated two different approaches to predict temporal anatomical changes using our model. First, given the CT image of a subject with sex and age , we predicted cranial anatomy at age using an additive scheme:

where represents the predicted value at time .

Our second approach incorporates the regressed age-specific standard deviation of the PCA coefficients to preserve population variability using a multiplicative scheme:

RESULTS

The number of principal components was selected as the minimum number at which the total explained variance converged. We used 35, 326, and 448 principal components to model cranial bone shape, thickness, and mineral density, respectively. (See figure, Supplemental Digital Content 1, which shows the cumulated percentage of variance as a function of the number of principal components to model (d) cranial bone shape, (e) thickness, and (f) HU values, http://links.lww.com/PRSGO/C124.) We use the inverse hyperbolic sine regression function to build our average model and linear polynomials to represent standard deviations. (See Supplemental Digital Content 1 (a)-(c), which shows the first PCA coefficient fitting using the Arcsinh formulation, http://links.lww.com/PRSGO/C124. See table, Supplemental Digital Content 2, which shows average model fitting error for the different PCA coefficient regression functions calculated using Eq. 3 and Eq. 2, http://links.lww.com/PRSGO/C125.) (See table, Supplemental Digital Content 3, which shows average squared fitting error of the temporal standard deviation model, http://links.lww.com/PRSGO/C126.) We achieved minimum fitting errors of 2.36 ± 1.36 mm, 0.01 ± 0.74 mm, and –3.01 ± 218.09 HU for cranial bone shape, thickness, and mineral density, respectively.

Figure 4 shows the evolution of the intracranial volume and cephalic index generated from the average cranial anatomies from our model. Sex differences were evaluated using the Student t test. Significant differences were found in the intracranial volume between boys and girls (from 0.38 0.07 l and 0.33 0.06 l at birth to 1.34 0.12 and 1.21 0.11 l at 10 years, respectively, P < 0.001). No significant differences were found in the cephalic index between boys (84.79% 4.54% at birth and 79.20% 3.94% at 10 years) and girls (83.58% 4.18% at birth and 78.83% 3.47% at 10 years), with P = 0.11 and P = 0.39 at birth and 10 years, respectively.

Fig. 4.

Average cephalic index (A) and intracranial volume (B) for boys and girls between birth and 10 years generated from the constructed model. Dashed lines represent the range of two standard deviations.

Table 2 summarizes the average cranial bone surface area, thickness, and mineral density at birth and 10 years (see their temporal distribution in SDC 4). (See figure, Supplemental Digital Content 4, which shows average evolution of cranial bone surface areas, thickness, and HU values for every age and sex generated from the constructed model, http://links.lww.com/PRSGO/C127.) Significant sex differences were found in all bone surface areas except for the occipital bone. There were no significant sex differences in thickness and mineral density except for the cranial thickness in the parietal bones at 10 years.

Table 2.

Cranial Bone Measurements for the Frontal, Parietal, and Occipital Bones at Birth and at 10 Years of Age

Metrics	Age	Boys	Girls	P
Frontal bones
Surface area (cm2)	Birth	44.91 ± 8.30	40.71 ± 7.25	<0.001*
	10 y	114.45 ± 10.00	105.70 ± 9.28	<0.001*
Thickness (mm)	Birth	1.13 ± 0.25	1.14 ± 0.32	0.07
	10 y	3.24 ± 0.68	3.02 ± 0.57	0.10
Mineral density (HU)	Birth	477.27 ± 92.78	479.84 ± 93.84	0.31
	10 y	1125.27 ± 113.08	1139.52 ± 112.35	0.52
Parietal bones
Surface area (cm2)	Birth	149.75 ± 13.99	139.18 ± 12.80	<0.001*
	10 y	322.45 ± 23.18	298.47 ± 19.66	<0.001*
Thickness (mm)	Birth	0.96 ± 0.21	1.00 ± 0.27	0.36
	10 y	3.31 ± 0.66	3.08 ± 0.57	0.03
Mineral density (HU)	Birth	486.66 ± 94.24	492.53 ± 99.35	0.33
	10 y	1082.72 ± 95.29	1087.74 ± 84.16	0.93
Occipital bone
Surface area (cm2)	Birth	40.34 ± 5.79	36.87 ± 5.07	0.21
	10 y	99.19 ± 8,86	95.44 ± 9.37	0.24
Thickness (mm)	Birth	1.28 ± 0.27	1.27 ± 0.31	0.07
	10 y	3.22 ± 0.63	3.05 ± 0.57	0.18
Mineral density (HU)	Birth	525.70 ± 94.91	540.99 ± 105.17	0.62
	10 y	1131.26 ± 99.66	1131.94 ± 96.35	0.24

Significant differences associated with sex for P < 0.01 are indicated with an asterisk.

Figure 5 shows the average suture fusion index for each suture. The average metopic suture fusion index reaches a value of 95% at 2.8 years and 2.5 years for boys and girls, respectively. The coronal suture fusion index reaches a value of 95% at 8.1 years and 8.7 years for boys and girls, respectively. The sagittal suture fusion index reaches a value of 95% at 7.7 and 8.2 years for boys and girls, respectively. Finally, the lambdoid suture fusion index remains below 95% at 10 years.

Fig. 5.

Evolution of the cranial suture fusion index with age for boys (A) and girls (B).

Our PCA model achieved similar fitting errors representing the first scans in the longitudinal dataset compared with the training dataset (P = 0.17, 0.26, and 0.45 for shape, thickness, and mineral density, respectively). The additive temporal predictive approach (Eq. 8) achieved errors of 2.98 ± 0.69 mm, 0.27 ± 0.29 mm, and 76.72 91.50 HU predicting shape, thickness, and mineral density changes, respectively, at the time of the second CT image. The multiplicative approach (Eq. 9) achieved errors of 2.98 ± 0.70 mm, 0.30 ± 0.33 mm, and 79.00 93.39 HU, respectively. Differences between both approaches were not statistically significant with P = 0.99, P = 0.71, and P = 0.91 for the three outcomes, respectively. We also achieved errors of 0.13 0.15 l and 0.13 0.14 l (P = 0.97) predicting cranial volume growth, and 0.83% ± 1.00% and 0.85% ± 0.98% (P = 0.92) predicting cephalic index changes for the two predictive schemes, respectively.

DISCUSSION

Research studies toward understanding and quantifying cranial pathology have been traditionally hindered by the lack of complete data-driven quantitative normative references of cranial bone development, leading to a relative paucity of studies on abnormal cranial growth compared with brain development and limiting our understanding of the coupled development between the brain and cranium. The model presented in this study provides a quantitative reference for the three most important aspects of cranial bone growth (shape, thickness, and mineral density) at every location of the calvaria during the first 10 years of life.

The intracranial volume and cephalic index derived from our model were aligned with previous clinical findings. They also show the importance of accounting for sex during cranial assessment. We show quantitatively that boys present a significantly higher intracranial volume and larger bone surface areas for the frontal and parietal bones. No statistical sex differences were found in bone mineral density and thickness. This suggests that the difference in cranial capacity may be related to a differential brain volume growth between boys and girls and not to differential bone growth and mineralization processes.

Our results highlight the importance of accounting for the stage of development when evaluating cranial shape anomalies. Although a cephalic index range between 76% and 85% is generally reported as a normative reference, our study shows that the cephalic index changes substantially with age from an average of 84.23% 4.31% at birth to an average of 79.01% 3.72% at 10 years. These quantitative findings agree with clinical reports suggesting a decreasing cephalic index with age. We also observed a high variability in its range in the normative population (Fig. 4), which explains its low-sensitivity detecting pathology and the variable normative ranges reported.

We also provide the first quantitative reference of bone mineral density changes and cranial suture fusion using standard HU values. This is of special importance in the assessment of suture anomalies such as craniosynostosis. Even though suture fusion is a local continuous ossification process, most studies address its evaluation as a binary problem, that is, sutures are reported as either fused or open, either completely or partially, based on a visual assessment from CT images. Because of this subjective evaluation, high discrepancies in the normal age of suture fusion are found in the literature, and our understanding of how sutures fuse normally and how pathology modifies this process is still limited.

As shown in our suture fusion index (Fig. 5) and reported in the literature, the mineral density of the metopic suture increases faster, and this suture fuses earlier than the other sutures. Although it is commonly reported that the metopic suture fuses normally between 9 and 24 months, we observe that this is a continuous process and, although mineral density changes are faster during the first year of life, this suture does not reach the mineral density of the frontal bones until close to 3 years. A similar pattern is observed in the other sutures that fuse later. It can also be observed from Figure 6 that the metopic suture area shows a more pronounced thickening at 10 years (3.45 mm) than the frontal bones (2.94 mm), a pattern that is also observed in the sagittal suture. This may be related to a larger available space at the brain’s longitudinal fissure and suggests the role of local brain pressure as a local modulator of bone development.

Fig. 6.

A, Average cranial bone thickness in the metopic suture and in the surrounding frontal bones. B, CT slice of a 4.58-year-old subject showing a substantial thickening at the metopic suture (indicated with a red arrow) compared with the frontal bones. C, CT slice of the same subject showing a similar thickening pattern in the sagittal suture indicated with a red arrow.

Interestingly, most sutures show a slight mineral density decrease after approximately 5 years. This is not a biophysical process, and this finding is due to partial volume effects in the CT images, which is a limitation of this study. Most of the cranial bone volume is represented by compact bone, and the spongy bone is usually not visible in the images of young children. However, as the cranial bone thickens, the notably less dense spongy bone can be seen in the images, and therefore, the average CT image intensity within the bone decreases (Fig. 6).

The cranial development model built from a unique large pediatric dataset presented in this study represents a necessary statistical reference of normality that can be used to identify and quantify abnormal cranial bone shape, thickness, and mineral density development associated with pathology. Normative references are an essential component to objectify patient evaluation and diagnosis and to quantify subtle anomalies that are hard to identify with subjective visual evaluation. Moreover, our quantitative references of mineral density at the cranial sutures will enable the study of the development of craniosynostosis, the causes of which are currently unknown for 70% of patients. Hence, it could lead to improved patient management using quantitative data to objectify and personalize diagnoses and treatments. Finally, given the generative nature of our models, they can also be used to create synthetic datasets for educational and research purposes.

Our model of cranial development at each location of the calvaria is publicly available at https://github.com/cuMIP/normativeCranialGrowth.