A physico-mechanical model of postnatal craniofacial growth in human

Summary Our fundamental understanding of the physico-mechanical forces that drive the size and shape changes of the cranium during ontogeny are limited. Biomechanical models based on finite element method present a huge opportunity to address this critical gap in our knowledge. Here, we describe a validated computational framework to predict normal craniofacial growth. Our results demonstrated that this approach is capable of predicting the growth of calvaria, face, and skull base. We highlighted the crucial role of skull base in antero-posterior growth of the face and also demonstrated the contribution of the maxillary expansion to the dorsoventral growth of the face and its interplay with the orbits. These findings highlight the importance of physical interactions of different components of the craniofacial system. The computational framework described here serves as a powerful tool to study fundamental questions in developmental biology and to advance treatment of conditions affecting the craniofacial system such as craniosynostosis.


Figures S1 to S12
Tables S1 to S4 SI references [S1-S7] (same with refs [7-9, 26, 27, 31, 33] in the Main Text)  Flowchart of six sensitivity tests including 16 independent analyses carried out in this study, where the main method (left side) and major improvement from each test (right side) are presented.Each sensitivity test was carried out to improve the predictions at one specific region through several case tests (see Figure S3-S9).The validated configurations or parameters from the former test (in yellow) were used as the base case for the next test.The optimized configurations from the sensitivity test 6 -Case 2 were used for final simulations.This figure is related to STAR Methods.

Figure S3
. Overview of sensitivity test 1.(a) This test was to investigate the contribution of pure radial expansion (αxyz) [S1] of four cranial volumes (intracranial volume-ICV, orbital volume-OV, nasal cavity volume-NCV and upper intraoral volume-UIV) in modeling the postnatal craniofacial growth up to 48 months of age.Nodal constraints in all degrees of freedom were placed around the foramen magnum, as in prior studies [S1,S2], and bone formation at the sutures was modelled using the tissue differentiation algorithm described in STAR Methods.(b) 2D morphological comparison of inner skeleton between pre-aligned in silico models and average skulls at 6, 12, 24, 36 and 48 months of age using cross-sections from three planes.(c)-(f) Absolute differences between the predicted four cranial volumes and target values at each age.The results show the well-predicted calvarial growth (see sagittal plane in (b)) and volumetric changes (c), but less satisfactory predictions across the face and skull base, these include: ( 1   ) by applying additional one-directional nodal displacement constraints to facial component covers (Case 1-3) to limit the undesirable volume expansion at covers (see 2D diagrams).Four local coordinate systems (in parallel with the global coordinate system) were created at the center of the external surface of each cover.Nodes on the external surface of each cover were selected and their coordinates were transformed to corresponding local coordinate systems for initialization (see 3D view).Displacement constraints in one specific direction were placed on selected nodes at the beginning of each load-step.The values of maximum nodal displacement applied to different covers in each load-step were characterized based on Method 1 (Figure S4) for Case 1, and Method 2 (Figure S4) for Case 2. In Case 3, the values for one-directional displacement constraints were optimized and defined through iterative tests at each load-step until no additional deformation of covers was observed (less than 10% changes in cover thickness).(b) 2D morphological comparison of inner and outer craniofacial skeleton between pre-aligned in silico models from all cases and average skulls at 6, 12, 24, 36 and 48 months of age shown as plots of cross-sections from three planes.(c)-(f) Absolute differences between predicted four cranial volumes from all cases and target values at each age.The results show that, by applying the nodal displacement constraints to facial component covers (OVC, NCVC and UIVC), there are no bulges observed at the covers due to the radial volume expansion (see (b)).Under the same expansion coefficients of four cranial volumes at each load-step from sensitivity test 1, the estimated values of displacement constraints applied in Case 1 and 2 led to significant changes in the predicted ICV, OV and UIV (see (c)-(e)), while Case 3 limited expansion of upper palate (see (f)).This figure is related to Figure S2 and STAR Methods. to optimize prediction in the spheno-occipital region.This was carried out through four case tests: (1) Case 1 -Nodal constraints in all degrees of freedom (DOF) were placed around the basilar portion of the occipital bone; (2) Case 2 -Nodal constraints in all DOF were placed around the presphenoid region, following a prior study [S3]; (3) Case 3 -Nodal constraints in one DOF were placed around the presphenoid region, while the one-directional displacement constraints (which allow a certain degree of displacement of constrained nodes in one direction) were placed around the foramen magnum.The foramen magnum displacement values (DO(l)) at each load-step were determined through iterative tests based on the measurement of D8 (see Figure S4), following the same method in sensitivity test 2 -Case 3; Case 4 -Nodal constraints in one DOF were placed around the basilar portion of the occipital bone with the same one-directional nodal displacement constraints placed around the foramen magnum.(b) 2D morphological comparison of inner skeleton between pre-aligned in silico models from all cases and average skulls at 6, 12, 24, 36 and 48 months of age using cross-sections from three planes.(αS,xyz ,αT,xyz ,αO(b),xyz and αO(l),xyz) were determined through iterative tests at each load-step, where the radial expansion coefficient was set to 0.1% at the beginning and then increased by 0.05% at each iteration, until the overall predicted shape matched the reference models at the skull base; (2) Case 2 -The sphenoid was assumed to expand directionally, more in the Y direction (antero-posterior direction) than in the X and Z directions (medio-lateral direction and superior-inferior direction).αS,y was set at 0.5% while αS,x and αS,z were both set at 0.1% at the beginning with the same increment rate (0.05%) at each iteration (other configurations were the same as Case 1); (3) Case 3 -The sphenoid was assumed to expand in Y and Z directions (other configurations were the same as Case 1).(b) 2D morphological comparison of inner skeleton between pre-aligned in silico models from all cases and average skulls at 6, 12, 24, 36 and 48 months of age using cross-sections from three planes.(2) Case 1 -Along with the directional facial volume expansion in Case 1, four facial bony components (upper portion of ethmoid bone-E(u), lateral wall of nasal cavity-LN, zygoma-Z and maxilla-M) were further expanded radially or directionally.Similar to sensitivity test 4 (Figure S7), the expansion coefficient of selected facial components was determined through iterative tests at each age (load-step), where radial or directional expansion coefficients were set at 0.1% and increased by 0.05% at each iteration until the predicted OV, NCV and UIV reached target values, and differences were less than 5% between predicted facial dimensions (upper facial height and biorbital breadth) and mean values from the normative dataset (n=217).(b) 2D morphological comparison of inner skeleton between pre-aligned in silico models from all cases and average skulls at 6, 12, 24, 36 and 48 months of age using cross-sections from three planes.(c)-(f) Absolute differences between four predicted cranial volumes from all cases and target values at each age.(g)-(j) Volume of bone components of all cases measured at each age.The results show that, the directional expansion of NCV and UIV in Case 1 led to better predictions of the nasal cavity region (see sagittal plane in (b)), but the size difference in the midfacial complex between predictions and the average models becomes prominent after 12 months of age (see coronal plane in (b)).In Case 2, radial and directional expansion of key facial components optimized the prediction of midfacial growth and the resulting facial shapes presented a better match to the average models up to 48 months of age (see sagittal and coronal planes in (b)).This figure is related to Figure S2 and STAR Methods.This test further expanded on sensitivity test 5 -Case 2 by altering the elastic modulus of NCV and UIV (referring to the volume of each defined space or cavity) to test the effect of material properties on predicted shapes.This was carried out through three case tests: (1) Case 1 -The elastic modulus (E) of NCV and UIV were set as 30MPa at 3 months of age (load-step 0) with an increment rate of 100 MPa/month (same as the material properties of suture components); (2) Case 2 -The E of NCV and UIV were set as 421MPa at 3 months of age (load-step 0) with E increased by 125 MPa/month (same as the material properties of bone components); (3) Case 3 -The E of NCV and UIV to 3000MPa at 3 months of age (load-step 0) with E increased by 125 MPa/month (an arbitrary extreme-large value of elastic modulus).Note that the expansion coefficients were adjusted to accommodate changes in material properties.(b) 2D morphological comparison of inner skeleton between pre-aligned in silico models from all cases and average skulls at 6, 12, 24, 36 and 48 months of age using cross-sections from four planes.(c)-(f) Absolute differences between four predicted cranial volumes from all cases and target values at each age.The results show that the difference in elastic modulus of NCV and UIV leads to non-significant effects on predicted facial shapes at specific ages.As expansion of these two volumes aims to expand the relevant regions to the correct shape and size rather than predicting the growth contributions of soft tissues such as the brain (intracranial volume), the configurations applied in Case 2 were used in final simulations.This figure is related to Figure S2 and STAR Methods.

Figure S12.
Changes in strain pattern across the skeleton and craniofacial sutures at 6, 12, 24, 36 and 48 months of age during growth.The first principal strains (1st µε) were derived from one load-step between the former age to the target age and mapped onto the in silico model at the target age.Skeleton and suture components were scaled individually with different ranges of strain before 12 months representing both bone formation at sutures and tissue differentiation, then the same strain range was shared by skeleton and suture components up to 48 months.This figure is related to Figure 5. Supplemental Tables Table S1 The anatomical structure of cribriform plate, crista galli, ethmoidal labyrinths and superior concha are simplified.Ethmoidal sinuses (air cells) were filled and assigned the same material properties as bone.
Lateral wall of nasal cavity b 21,22 LN a The anatomical structures of orbital plate and middle nasal concha were simplified and combined in current component.

Nasal septum b 23 NS
The anatomical structures of the perpendicular plate of ethmoid bone, vomer bone and septal cartilage were simplified and combined in current component.Segmented as one plate which covers the upper intraoral region through the alveolar margin.The edges or the contact interfaces of all above components have been flattened, and the width of suture, synchondrosis and fontanelle components has been increased by 2 -4 mm.Common craniofacial components and anatomical abbreviations are sourced from [S4-S7] unless otherwise noted: a Abbreviations defined by authors.b Component contain more than one sutures or bones.
Figure S1.Visualization of the final in silico baseline model and initial reconstruction of in vivo skull from a female individual at around 3 months of age: (a) 3D views.(b) 2D cross-sections in four planes.Note that details of the separated components of in silico model and the planes used for cross-sections are described in Figure 1 in the main text.This figure is related to Figure 1.

Figure S2 .
Figure S2.Flowchart of six sensitivity tests including 16 independent analyses carried out in this study, where the main method (left side) and major improvement from each test (right side) are presented.Each sensitivity test was carried out to improve the predictions at one specific region through several case tests (see FigureS3-S9).The validated configurations or parameters from the former test (in yellow) were used as the base case for the next test.The optimized configurations from the sensitivity test 6 -Case 2 were used for final simulations.This figure is related to STAR Methods.
Figure S3.Overview of sensitivity test 1.(a) This test was to investigate the contribution of pure radial expansion (αxyz) [S1] of four cranial volumes (intracranial volume-ICV, orbital volume-OV, nasal cavity volume-NCV and upper intraoral volume-UIV) in modeling the postnatal craniofacial growth up to 48 months of age.Nodal constraints in all degrees of freedom were placed around the foramen magnum, as in prior studies [S1,S2], and bone formation at the sutures was modelled using the tissue differentiation algorithm described in STAR Methods.(b) 2D morphological comparison of inner skeleton between pre-aligned in silico models and average skulls at 6, 12, 24, 36 and 48 months of age using cross-sections from three planes.(c)-(f) Absolute differences between the predicted four cranial volumes and target values at each age.The results show the well-predicted calvarial growth (see sagittal plane in (b)) and volumetric changes (c), but less satisfactory predictions across the face and skull base, these include: (1) unexpected volume expansions displayed as bulges at facial component covers (OVC, NCVC and UIVC); (2) the differences in size and position of the spheno-occipital region between in silico models and the average in vivo individual became significant after 12 months (see sagittal plane in (b)); (3) a lack of prediction of skull base growth (i.e.temporal-occipital region; see coronal plane in (b)) after 12 months of age; (4) obvious mismatch of the entire midfacial complex after 6 months of age, including less growth of the nasal cavity and upper palatal region in the supero-inferior direction, and insufficient antero-posterior growth of the upper facial region (see sagittal and transverse planes in (b)).All insufficient predictions identified in this test were improved by the following tests, where sensitivity test 2 addresses problem (1), sensitivity test 3 addresses problem (2), sensitivity test 4 addresses problem (3), and sensitivity test 5&6 address problem (4).This figure is related to Figure S2 (sensitivity test 1) and STAR Methods.

Figure S4 .
Figure S4.Measurements of eight linear dimensions used as references for the values of nodal displacement constraints applied in sensitivity tests 2 and 3: (a) six average skulls with three planes defined for measurements, and two methods of calculating the values of displacement (D) at each load-step (age) derived from linear measurements.D1-D2 measured in the mid-sagittal plane between two aligned skulls (different ages) indicate relative anterior movement of the nasal complex.Similarly, D3-D4 and D5-D6 measured in the sagittal and transverse planes indicate relative anterior growth of the orbital rim.D7 was measured to quantify palatal movement.(b) D1-D8 measured based on Method 1: the maximum values of D1-D8 between the 3 and 48 months of age skull models were measured from the grid plot, and the average displacement for each month were calculated through linear interpolation of the total 45 months of intervals.(c)-(g) D1-D8 measured for method 2: the maximum values of D1-D8 between each of the two skull models at target ages (i.e. between 3 and 6 months) were measured from the corresponding grid plot, then the average value was calculated through the linear interpolation of the corresponding age interval (i.e. 3 months).This figure is related to Figure S5 and STAR Methods.

Figure S5 .
Figure S5.Overview of sensitivity test 2. (a) This test further expanded on sensitivity test 1 (Base case) by applying additional one-directional nodal displacement constraints to facial component covers (Case 1-3) to limit the undesirable volume expansion at covers (see 2D diagrams).Four local coordinate systems (in parallel with the global coordinate system) were created at the center of the external surface of each cover.Nodes on the external surface of each cover were selected and their coordinates were transformed to corresponding local coordinate systems for initialization (see 3D view).Displacement constraints in one specific direction were placed on selected nodes at the beginning of each load-step.The values of maximum nodal displacement applied to different covers in each load-step were characterized based on Method 1 (FigureS4) for Case 1, and Method 2 (FigureS4) for Case 2. In Case 3, the values for one-directional displacement constraints were optimized and defined through iterative tests at each load-step until no additional deformation of covers was observed (less than 10% changes in cover thickness).(b) 2D morphological comparison of inner and outer craniofacial skeleton between pre-aligned in silico models from all cases and average skulls at 6, 12, 24, 36 and 48 months of age shown as plots of cross-sections from three planes.(c)-(f) Absolute differences between predicted four cranial volumes from all cases and target values at each age.The results show that, by applying the nodal displacement constraints to facial component covers (OVC, NCVC and UIVC), there are no bulges observed at the covers due to the radial volume expansion (see (b)).Under the same expansion coefficients of four cranial volumes at each load-step from sensitivity test 1, the estimated values of displacement constraints applied in Case 1 and 2 led to significant changes in the predicted ICV, OV and UIV (see (c)-(e)), while Case 3 limited expansion of upper palate (see (f)).This figure is related to FigureS2and STAR Methods.

Figure S6 .
Figure S6.Overview of sensitivity test 3. (a) This test further expanded on sensitivity test 2 -Case 3 by altering nodal constraints on the skull base components (including sphenoid-S, basilar and later parts of occipital bone-O(b) and O(l))to optimize prediction in the spheno-occipital region.This was carried out through four case tests: (1) Case 1 -Nodal constraints in all degrees of freedom (DOF) were placed around the basilar portion of the occipital bone; (2) Case 2 -Nodal constraints in all DOF were placed around the presphenoid region, following a prior study[S3]; (3) Case 3 -Nodal constraints in one DOF were placed around the presphenoid region, while the one-directional displacement constraints (which allow a certain degree of displacement of constrained nodes in one direction) were placed around the foramen magnum.The foramen magnum displacement values (DO(l)) at each load-step were determined through iterative tests based on the measurement of D8 (see FigureS4), following the same method in sensitivity test 2 -Case 3; Case 4 -Nodal constraints in one DOF were placed around the basilar portion of the occipital bone with the same one-directional nodal displacement constraints placed around the foramen magnum.(b) 2D morphological comparison of inner skeleton between pre-aligned in silico models from all cases and average skulls at 6, 12, 24, 36 and 48 months of age using cross-sections from three planes.(c)-(f) Absolute differences between predicted four cranial volumes from all cases and target values at each age.The results show that, nodal constraints in all DOF placed on the basilar portion of occipital bone (Case 1) and presphenoid region (Case 2) led to excessive local deformation of the relevant bones, while rotation of the spheno-occipital region was observed in both case tests compared with the Base Case (see (b)).By combining the one-directional nodal displacement constraints on the foramen magnum with either nodal constraints in one direction on the presphenoid region (Case 3) or on the basilar portion of the occipital bone (Case 4), the excessive local deformation due to over constraining was eliminated (see (b)).The predicted cranial volumes in Case 4 are closer to the target volume than the others (see (c)-(f)).This figure is related to FigureS2and STAR Methods.
Figure S6.Overview of sensitivity test 3. (a) This test further expanded on sensitivity test 2 -Case 3 by altering nodal constraints on the skull base components (including sphenoid-S, basilar and later parts of occipital bone-O(b) and O(l))to optimize prediction in the spheno-occipital region.This was carried out through four case tests: (1) Case 1 -Nodal constraints in all degrees of freedom (DOF) were placed around the basilar portion of the occipital bone; (2) Case 2 -Nodal constraints in all DOF were placed around the presphenoid region, following a prior study[S3]; (3) Case 3 -Nodal constraints in one DOF were placed around the presphenoid region, while the one-directional displacement constraints (which allow a certain degree of displacement of constrained nodes in one direction) were placed around the foramen magnum.The foramen magnum displacement values (DO(l)) at each load-step were determined through iterative tests based on the measurement of D8 (see FigureS4), following the same method in sensitivity test 2 -Case 3; Case 4 -Nodal constraints in one DOF were placed around the basilar portion of the occipital bone with the same one-directional nodal displacement constraints placed around the foramen magnum.(b) 2D morphological comparison of inner skeleton between pre-aligned in silico models from all cases and average skulls at 6, 12, 24, 36 and 48 months of age using cross-sections from three planes.(c)-(f) Absolute differences between predicted four cranial volumes from all cases and target values at each age.The results show that, nodal constraints in all DOF placed on the basilar portion of occipital bone (Case 1) and presphenoid region (Case 2) led to excessive local deformation of the relevant bones, while rotation of the spheno-occipital region was observed in both case tests compared with the Base Case (see (b)).By combining the one-directional nodal displacement constraints on the foramen magnum with either nodal constraints in one direction on the presphenoid region (Case 3) or on the basilar portion of the occipital bone (Case 4), the excessive local deformation due to over constraining was eliminated (see (b)).The predicted cranial volumes in Case 4 are closer to the target volume than the others (see (c)-(f)).This figure is related to FigureS2and STAR Methods.

Figure S7 .
Figure S7.Overview of sensitivity test 4. (a) This test further expanded on sensitivity test 3 -Case 4 by applying additional (αxyz) or directional(αx, αy, αz)  expansion to skull base components, including the sphenoid-S, temporal bone-T, basilar and lateral parts of occipital bone-O(b) and O(l), to optimize the prediction of the entire skull base region.This was carried out through three case tests: (1) Case 1 -S, T, O(b) and O(l)) were assumed to expand radially, and the expansion coefficients for each bony component(αS,xyz ,αT,xyz ,αO(b),xyz and αO(l),xyz) were determined through iterative tests at each load-step, where the radial expansion coefficient was set to 0.1% at the beginning and then increased by 0.05% at each iteration, until the overall predicted shape matched the reference models at the skull base; (2) Case 2 -The sphenoid was assumed to expand directionally, more in the Y direction (antero-posterior direction) than in the X and Z directions (medio-lateral direction and superior-inferior direction).αS,y was set at 0.5% while αS,x and αS,z were both set at 0.1% at the beginning with the same increment rate (0.05%) at each iteration (other configurations were the same as Case 1); (3) Case 3 -The sphenoid was assumed to expand in Y and Z directions (other configurations were the same as Case 1).(b) 2D morphological comparison of inner skeleton between pre-aligned in silico models from all cases and average skulls at 6, 12, 24, 36 and 48 months of age using cross-sections from three planes.(c)-(f) Absolute differences between predicted four cranial volumes from all cases and target values at each age.(g)-(j) Volumes of bone components of all cases measured at each age.The results show that, by applying additional radial expansion to the skull base bones, Case 1 over predicts growth of the sphenoid and temporal bones in the medio-lateral direction, while Cases 2 and 3 show good predictions of the overall size of the sphenoid (see sagittal plane in (b)).When no medio-lateral expansion is applied to sphenoid, Case 3 shows better prediction of spheno-temporal regions (see coronal plane in (b)).This figure is related to FigureS2and STAR Methods.
Figure S7.Overview of sensitivity test 4. (a) This test further expanded on sensitivity test 3 -Case 4 by applying additional (αxyz) or directional(αx, αy, αz)  expansion to skull base components, including the sphenoid-S, temporal bone-T, basilar and lateral parts of occipital bone-O(b) and O(l), to optimize the prediction of the entire skull base region.This was carried out through three case tests: (1) Case 1 -S, T, O(b) and O(l)) were assumed to expand radially, and the expansion coefficients for each bony component(αS,xyz ,αT,xyz ,αO(b),xyz and αO(l),xyz) were determined through iterative tests at each load-step, where the radial expansion coefficient was set to 0.1% at the beginning and then increased by 0.05% at each iteration, until the overall predicted shape matched the reference models at the skull base; (2) Case 2 -The sphenoid was assumed to expand directionally, more in the Y direction (antero-posterior direction) than in the X and Z directions (medio-lateral direction and superior-inferior direction).αS,y was set at 0.5% while αS,x and αS,z were both set at 0.1% at the beginning with the same increment rate (0.05%) at each iteration (other configurations were the same as Case 1); (3) Case 3 -The sphenoid was assumed to expand in Y and Z directions (other configurations were the same as Case 1).(b) 2D morphological comparison of inner skeleton between pre-aligned in silico models from all cases and average skulls at 6, 12, 24, 36 and 48 months of age using cross-sections from three planes.(c)-(f) Absolute differences between predicted four cranial volumes from all cases and target values at each age.(g)-(j) Volumes of bone components of all cases measured at each age.The results show that, by applying additional radial expansion to the skull base bones, Case 1 over predicts growth of the sphenoid and temporal bones in the medio-lateral direction, while Cases 2 and 3 show good predictions of the overall size of the sphenoid (see sagittal plane in (b)).When no medio-lateral expansion is applied to sphenoid, Case 3 shows better prediction of spheno-temporal regions (see coronal plane in (b)).This figure is related to FigureS2and STAR Methods.

Figure S8 .
Figure S8.Overview of sensitivity test 5.(a) This test further expanded on sensitivity test 4 -Case 3 by applying additional (αxyz) or directional(αx, αy, αz)  expansion to facial components, to optimize the prediction of the midfacial region.This was carried out through two case tests: (1) Case 1 -Directional expansion applied to NCV and UIV; (2) Case 1 -Along with the directional facial volume expansion in Case 1, four facial bony components (upper portion of ethmoid bone-E(u), lateral wall of nasal cavity-LN, zygoma-Z and maxilla-M) were further expanded radially or directionally.Similar to sensitivity test 4 (FigureS7), the expansion coefficient of selected facial components was determined through iterative tests at each age (load-step), where radial or directional expansion coefficients were set at 0.1% and increased by 0.05% at each iteration until the predicted OV, NCV and UIV reached target values, and differences were less than 5% between predicted facial dimensions (upper facial height and biorbital breadth) and mean values from the normative dataset (n=217).(b) 2D morphological comparison of inner skeleton between pre-aligned in silico models from all cases and average skulls at 6, 12, 24, 36 and 48 months of age using cross-sections from three planes.(c)-(f) Absolute differences between four predicted cranial volumes from all cases and target values at each age.(g)-(j) Volume of bone components of all cases measured at each age.The results show that, the directional expansion of NCV and UIV in Case 1 led to better predictions of the nasal cavity region (see sagittal plane in (b)), but the size difference in the midfacial complex between predictions and the average models becomes prominent after 12 months of age (see coronal plane in (b)).In Case 2, radial and directional expansion of key facial components optimized the prediction of midfacial growth and the resulting facial shapes presented a better match to the average models up to 48 months of age (see sagittal and coronal planes in (b)).This figure is related to FigureS2and STAR Methods.

Figure S9 .
Figure S9.Overview of sensitivity test 6.(a) This test further expanded on sensitivity test 5 -Case 2 by altering the elastic modulus of NCV and UIV (referring to the volume of each defined space or cavity) to test the effect of material properties on predicted shapes.This was carried out through three case tests: (1) Case 1 -The elastic modulus (E) of NCV and UIV were set as 30MPa at 3 months of age (load-step 0) with an increment rate of 100 MPa/month (same as the material properties of suture components); (2) Case 2 -The E of NCV and UIV were set as 421MPa at 3 months of age (load-step 0) with E increased by 125 MPa/month (same as the material properties of bone components); (3) Case 3 -The E of NCV and UIV to 3000MPa at 3 months of age (load-step 0) with E increased by 125 MPa/month (an arbitrary extreme-large value of elastic modulus).Note that the expansion coefficients were adjusted to accommodate changes in material properties.(b) 2D morphological comparison of inner skeleton between pre-aligned in silico models from all cases and average skulls at 6, 12, 24, 36 and 48 months of age using cross-sections from four planes.(c)-(f) Absolute differences between four predicted cranial volumes from all cases and target values at each age.The results show that the difference in elastic modulus of NCV and UIV leads to non-significant effects on predicted facial shapes at specific ages.As expansion of these two volumes aims to expand the relevant regions to the correct shape and size rather than predicting the growth contributions of soft tissues such as the brain (intracranial volume), the configurations applied in Case 2 were used in final simulations.This figure is related to FigureS2and STAR Methods.

Figure S10 .
Figure S10.Scatter plots of the comparisons of the measurement of upper facial index (Upper facial height / Bizygomatic breadth x 100) between in silico models (in red) and in vivo average skulls (in blue) at 6, 12, 24, 36, and 48 months of age.Note that the corresponding data from the control group of 217 normal individuals (in grey) are also plotted in this scatter plot, indicating the overall changes and distributions in the normative dataset.This figure is related to Figure4.

Figure S11 .
Figure S11.Sensitivity tests of hydrostatic strain range.3D model highlights the overall hydrostatic strain across the in silico model after a single growth load-step from 3 to 4 months of age.The dashed box (in blue) shows the most applicable threshold of hydrostatic strain for this study.Note that regions shown in grey represent regions below or above the specified thresholds.This figure is related to STAR Methods.

.
Details of in silico baseline model including definitions of each component, abbreviations, and the assumptions or simplifications made during model development.This table is related to Figure 1.

Table S2 .
Details of 12 cranial measurements, 1 module and 3 indexes used in this study.This table is related to Figure2and STAR Methods.
a Defined by authors.

Table S3 .
(a) Four cranial volumes (intracranial volume-ICV, orbital volume-OV, nasal cavity volume-NCV and upper intraoral volume-UIV) predicted at each load-step (age) through regression analyses of the normative in vivo dataset (n=217).(b)Details of individuals selected as a baseline model at 3 months and as average skulls at 6, 12, 24, 36, and 48 months of age according to cranial volumes.This table is related to Figure2and STAR Methods.