Modeling the effect of brain growth on cranial bones using finite‑element analysis and geometric morphometrics

Purpose Brain expansion during ontogeny has been identified as a key factor for explaining the growth pattern of neurocranial bones. However, the dynamics of this relation are only partially understood and a detailed characterization of integrated morphological changes of the brain and the neurocranium along ontogeny is still lacking. The aim of this study was to model the effect of brain growth on cranial bones by means of finite-element analysis (FEA) and geometric morphometric techniques. Methods First, we described the postnatal changes in brain size and shape by digitizing coordinates of 3D semilandmarks on cranial endocasts, as a proxy of brain, segmented from CT-scans of an ontogenetic sample. Then, two scenarios of brain growth were simulated: one in which brain volume increases with the same magnitude in all directions, and other that includes the information on the relative expansion of brain regions obtained from morphometric analysis. Results Results indicate that in the first model, in which a uniform pressure is applied, the largest displacements were localized in the sutures, especially in the anterior and posterior fontanels, as well as the metopic suture. When information of brain relative growth was introduced into the model, displacements were also concentrated in the lambda region although the values along both sides of the neurocranium (parietal and temporal bones) were larger than under the first scenario. Conclusion In sum, we propose a realistic approach to the use of FEA based on morphometric data that offered different results to more simplified models.


Introduction
The growth of vault bones is closely related to the increase of brain volume during early ontogeny [13,24].The pressure exerted by the brain would be, at least partially, responsible for this interaction with the skull [19].Although the specific mechanisms are still not well understood, different studies suggest that mechanical pressure resulting from brain expansion stimulates bone cells to proliferate and differentiate, as well as to produce extracellular matrix [22].This is specially the case of intramembranous bone formation in the sutures of vault bones.The influence of soft tissues on neurocranium becomes evident in the context of different pathological malformations, such as microcephaly and macrocephaly, which primarily affect brain development but also produce morphological alterations on skull bones [19].Moreover, some features of the human skull are thought to be a result of the expansion of the brain along hominin evolution [2].Consequently, the study of the structural and functional relations between the brain and skull along ontogeny is relevant in biomedical as well as evolutionary contexts.
The responses of cranial bones to changes in brain size and shape can be approached by models that simulate the interactions between these two structures under different conditions.Such models can contribute to broaden the current knowledge on morphological changes of the skull along human evolution but also to the diagnosis and treatment of pathologies affecting craniofacial development, such as the premature fusion of sutures (craniosynostosis) and the alteration of normal brain growth (micro, macrocephaly), among others [1,13,19].In this sense, finite-element analysis (FEA) is a promising alternative as was shown in studies that model infant cranial bones with the aim of simulating the response of the growing skull to mechanical pressure and hits in utero [11].However, the use of FEA to assess the effect of brain growth on the skull has not been fully explored yet [10,12,15].
The aim of this study was to assess by means of the finite-element method the influence of brain expansion on cranial bones during early postnatal ontogeny under two different scenarios of brain growth.In the first scenario, brain volume increases with the same magnitude in all directions, while the second one takes into account the relative growth of the different brain regions.Here, geometric morphometric techniques were used to quantify morphological variation in endocasts of a postnatal ontogenetic sample.Endocasts have been widely used in paleontological studies and contexts, where soft tissues are not preserved, and they proved to be a good proxy of overall brain size and shape [5].The captured variation was then used as an input for generating a finite-element model of cranial bone displacements under the effect of brain growth.

Sample composition
We analyzed computed tomography (CT) images of the skull of 23 subadult individuals.They belong to the Bosma Collection [26].The skulls were scanned in the Johns Hopkins University of Baltimore (USA) using a GE Medical Systems Genesis Jupiter scanner and the following parameters: 512 × 512 matrix size, 0.039 mm voxel size and 1.5 mm slice thickness [18].Images were generously shared by Dr. Joan Richtsmeier (Department of Anthropology, College of the Liberal Arts, Pennsylvania State University).
Age was estimated in each case by examining the patterns of dental formation and eruption [6].According to these estimations, the age at death of these individuals range between 0 and 12 years (Table 1).

Analyses of morphological variation of the endocasts
With the aim of estimating the magnitude and pattern of morphological variation during postnatal growth in the brain, we manually segmented the endocranial cavity of each specimen and obtained a 3D endocast of the neurocranium (Fig. 1).This procedure was performed manually using 3D Slicer [8].
Size changes in the brain along ontogeny were analyzed using the volume of the endocasts, while shape changes were described by a set of landmarks and semilandmarks.First, 21 landmarks were digitized on each virtual endocast, following definitions summarized in Table 2 and illustrated in Fig. 1b [3].Since the endocranium is characterized by extensive surfaces lacking identifiable homologous discrete landmarks, a complementary set of surface semilandmarks was used to improve the coverage of the structure of interest.For this  purpose, 2500 semilandmarks were obtained from one of the meshes using the Poisson Disk Sampling and Element Mesh Sampling functions in MeshLab (Fig. 1c).These semilandmarks were projected from the reference to the target surfaces-aligned by their homologous landmarks-using the placePatch function in Morpho package [25].Finally, the set of 2521 points were aligned by a Generalized Procrustes superimposition and the semilandmarks slid by the minimum Procrustes distance criteria.In sum, the set of 21 landmarks was first used to align the semilandmarks and, then, landmarks and semilandmarks were used jointly to describe morphological changes in the endocranium.Since we expect brain effect on cranial bones to be due to size and shape changes, we analyzed variation in form by multiplying superimposed coordinates (describing shape) by centroid size (describing size).On these coordinates of form, we performed a principal component analysis (PCA).These analyses were carried out using geomorph package for R.

Finite-element analysis
A structural static analysis using the FEA Package ANSYS 17.1 in a Dell Precision™ Workstation T5500 with 48 GB and 5.33 GHz was performed in the youngest specimen (Bosma 2).To study the deformation of the neurocranial bones of a neonate individual we simulated the effect of brain growth using two different loading scenarios (see skull reconstruction in Online Resource 1).First, the segmented model was converted to a CAD model.During this step, the irregularities in the surface that were due to the scanning process were repaired using refinement and smoothing tools from Geomagic (3D Systems).The sutures were also defined in the cranial vault (Fig. 2).
Second, we created the FEA model to compute the displacements in the two different neurocranial bones loading conditions, which reflect two different simulated scenarios.In the first scenario, the growth was simulated by applying a pressure of 1 MPa perpendicular to the inner surface.In this scenario, a non-homogeneous effect is expected due to the geometry of the neurocranium and to the differences in the properties of bone and sutures.In the second scenario, different perpendicular pressures were applied in relation to the relative growth of brain regions: temporal, parietal (superior and inferior), occipital, frontal and cerebellum (Fig. 2).To obtain growth vectors for each brain region, we performed a superimposition of the mean configuration of points of neonates with a 1.5-year-old individual.The coordinates were separated into six regions: temporal lobe, frontal lobe, occipital lobe, superior and lateral parietal lobe and cerebellum.A superposition was performed through Protest without correction by size and the mean of the residuals by region was estimated to describe the magnitude of change in form, or relative growth, between 0 and 1.5 years.In the post-processing stage, we obtained displacement measures from the values in each node of the mesh.
For the two scenarios, boundary conditions were defined to represent the loads and fixed displacements.We established fixed displacements in the occipital condyles needed to avoid movements in the model when the pressures are applied (Fig. 2).Elastic, linear and homogeneous material properties were assumed for the bone using values from previous studies [10,17].Finally, the neurocranial skull was meshed using an adaptive mesh of hexahedral elements [16] with 123.338 elements describing vault and cranial base bones as well as sutures was generated.

Cortical bone Sutures Fixed displacements
Fig. 2 Model created for the finite-element analysis including the fixed displacements and the regions with bone and sutures

Results
Endocast volume increased strikingly during the first years of life, while it stabilized around 3 years, growing from ~ 300 to 1300 cc 3 (Fig. 3a).The first axis of the PC1 represents variation in form associated with age (Fig. 3b).
Along PC1, two groups can be distinguished: one is composed by neonates and the other comprises the rest of the sample (~ from 1.5 to 13 years) (Fig. 3b).Brain expansion with age, as is summarized along PC1, is characterized by a large anterior expansion of frontal lobes, which acquire a more rounded shape, as well as a lateral enlargement of parietal lobes (Fig. 3c).The average residuals that summarize the distance between nearest points of the mean configuration of neonate endocasts and a 1.5-year-old individual are shown in Table 3.The cerebellum is the region with the largest average residual, indicating that this region went under a larger relative growth compared with the other lobes.The temporal lobe, in contrast, displayed the smallest residuals during this period.The average residuals were transformed into proportions by assigning a value of 1 to the cerebellum and proportional values to the rest of the lobes, which were then used as relative pressures under scenario 2 for the FEA (Table 3).
Displacements of neurocranial bones under the pressure of the brain were simulated by finite-element models using the two different scenarios for brain pressures.Under the first scenario, in which uniform pressures were applied throughout the neurocranium, we found that the largest displacements were located around the superior portion of the parietal bone and the lambda point (Fig. 4).In the second scenario, in which pressures were proportional to the relative growth of brain lobes, displacements were also concentrated in the lambda region but the values along the lateral walls of the lateral portions of the neurocranium were larger than under a homogeneous pressure (Fig. 4).In both cases, the lowest values of displacements were in the cranial base.Accordingly, the index of these displacements, estimated at the points used to measure the cranial width and length, is close to 0.8 mm when a non-uniform pressure is applied, while in the first scenario is 0.6 mm.

Discussion
The pressure exerted by the growing brain represents an important stimulus for the growth of the skull, especially the vault.The morphometric analysis carried out in this study indicated that during the first year of life the endocranial volume undergoes a pronounced increase, which is less evident after 3 years.This result agrees with previous studies based on endocasts obtained from CT scans and MRI of healthy brains [7,21].Differences in form were also more pronounced between neonates and older ages.The main agerelated shape changes consisted in the relative expansion of the corresponding areas of the frontal and occipital lobes and the cerebellum.A similar pattern was described by Neubauer et al. [21] for the ontogeny of the endocranial changes between the neonatal stage and the first year, which was characterized by a high growth rate of the posterior fossa, probably associated with the growth of the cerebellum.Changes in form found here between neonates and older individuals also included a lateral expansion of the parietal lobes producing a more globular endocranial configuration.This phase of brain globularization has been suggested as a singular feature of the human brain ontogeny [9].In sum, our results as well as previous studies point out a localized growth of brain regions throughout postnatal ontogeny.The two finite-element analyses performed here showed differences in the magnitude and location of displacement of the bones and sutures of the neurocranium.In scenario 1, which applied a uniform pressure to the cranial bones and sutures, the largest displacements were observed along the anterior fontanel and the metopic suture, and in the posterior fontanel.Under scenario 2 of differential growth by brain region, the largest displacements were observed in the posterior fontanel, whereas the displacements in the anterior region were much smaller and similar across the frontal bone.In addition, a larger lateral displacement was observed along the temporo-parietal bones.Given the different properties of bones and sutures it is not surprising that most of the changes in a neonate skull are concentrated on the sutures, although this is more evident under the assumption that brain pressures are uniform.In a previous study, Jin and colleagues noted [10] that the brain expansion forces are not expected to be uniform as the brain develops within the non-uniform concavities of the skull.To take the non-uniform effect of the brain into account, they proposed a finite-element model in which the forces on the cranial bones and sutures develop in response to the pressure from the expanding volume of the brain.However, the simulated expansion of the brain volume was uniform across all dimensions.In the same way, a more recent finite-element analysis used an isotropic expansion of the intra-cranial volume to assess the effect of the growing brain on cranial bones [12].As we showed here, in agreement with other studies [9,21], the lobes of the brain have a differential growth rate postnatally, which also contributes to the non-uniform pressure on the cranial bones.Such models might be improved by incorporating the differential expansion of the brain regions.An alternative to explore in that direction is the use of geometric morphometrics, which provides an accurate description of changes in absolute and relative size (i.e., shape).As we showed here, the residuals that quantify differences in brain form between age stages may serve as an input for finite-element analyses.
The displacements of cranial bones under the second scenario modeled here showed a larger increase in cranial length than breadth, although the differences are not as pronounced as under the uniform pressure.This means that a constant pressure would result in a much longer skull.The relative larger increase in length compared to width is consistent with the observed postnatal changes of the cephalic index that drops from an average of 0.8-0.78from birth to 3 years [14].The hybrid model developed by Jin and colleagues [10] resulted in an increase of the cephalic index with the expansion of the brain, contrasting with our results as well as with pediatric data.Such discrepancies might be due to the specific anatomy of the model used given that the simulations were based on a single specimen.Further models would benefit from using an atlas built from a sample of individuals instead of one particular case.The use of atlases is common in neuroimage studies, although is still less developed for the skull probably due to the scarcity of automatic methods for quantifying the cranial morphology, with some exceptions (e.g., [23]).
Models developed in previous studies as well as the proposed here kept constant the material properties, the size and shape of the cranial bones and sutures [10,12].Because such properties change with age, these models could only accurately represent a rather limited period of time, particularly for early postnatal growth when rapid changes take place.Therefore, further models that simulate the effect of brain expansion at different ages are needed to take into account the changes in bone shape and the progressive closure of sutures with age.This will also enable the simulation of changes in the growth pattern of the skull under different types of abnormal closure of sutures (i.e., craniosynostosis).In these cases, skull growth usually results in a reduced development perpendicular to the fused sutures and a compensatory overgrowth along the sutures that are open, as postulated by Virchow's Law [27].Our study could be used as a baseline to develop detailed simulations of different types of craniosynostosis under realistic assumptions.Previous studies have applied FEA to virtually predict surgical outcomes on skull growth in cases of craniosynostosis after surgery (e.g., [4,29,30]).There are few examples, however, analyzing growth dynamics of skulls with craniosynostosis and they either focus on specific craniofacial regions (e.g., [20]) or use simplified neurocranial geometries to explore displacements in different types of craniosynostosis (e.g., [28]).Further work would be essential to link the apportionments of our model, which is based on realistic growth parameters, with recent advances to formally analyze physical and geometrical constraints in patients with craniosynostosis (e.g., [28]).
In sum, the results of this study suggest that the conjoint use of geometric morphometrics and finite-element analysis is a promising approach for developing models aimed at simulating the effect of brain growth on cranial bones.These model approaches can be of great help for medical practice in planning and monitoring surgical interventions of the skull.

Fig. 1 2
Fig.1Steps for segmentation and reconstruction of endocasts from CT images and landmark digitization.a From the CT scan of each skull, the internal cavity of the neurocranium was identified in each For bone we used E (Young modulus) = 1300 MPa and v (Poisson coefficient) = 0.28, and for sutures we used E (Young modulus) = 200 MPa and v (Poisson coefficient) = 0.28.

Fig. 3 a
Fig. 3 a Age against endocast volume.b Age against scores of the first principal component (PC) derived from the principal component analysis (PCA) of form (shape variables plus centroid size).c Morphings depicting the changes in endocast form towards the negative and positive extremes of the PC1 (left).Heatmaps in the right column represent the point distance between the endocast of the negative extreme (youngest specimens) and the positive extreme of the PC1.The distances are displayed on a neonate endocast

Fig. 4
Fig. 4 Displacements under the scenario of equal pressures (scenario 1, top) and different pressures across brain regions (scenario 2, bottom)

Table 1
Sample composition.

Table 3
Average of residuals between 0 and 1.5 years and applied pressures in MPa in the inner part of the FEA model skull