A coupled reaction-diffusion-strain model predicts cranial vault formation in development and disease

Abstract

How cells utilize instructions provided by genes and integrate mechanical forces generated by tissue growth to produce morphology is a fundamental question of biology. Dermal bones of the vertebrate cranial vault are formed through the direct differentiation of mesenchymal cells on the neural surface into osteoblasts through intramembranous ossification. Here we join a self-organizing Turing mechanism, computational biomechanics, and experimental data to produce a 3D representative model of the growing cerebral surface, cranial vault bones, and sutures. We show how changes in single parameters regulating signaling during osteoblast differentiation and bone formation may explain cranial vault shape variation in craniofacial disorders. A key result is that toggling a parameter in our model results in closure of a cranial vault suture, an event that occurred during evolution of the cranial vault and that occurs in craniofacial disorders. Our approach provides an initial and important step toward integrating biomechanics into the genotype phenotype map to explain the production of variation in head morphology by developmental mechanisms.

Keywords: Brain; Computational morphogenesis; Craniosynostosis; Finite volume method; Intramembranous ossification; Mouse model; Skull growth and evolution.

Fig. 1

a Schematic diagram of the multi-scale reaction-diffusion-strain (RDE) model for cranial vault bone formation. b Schematic diagram showing the computational process. Structural analysis gives the strain field that is used to estimate the distribution of molecules and cells using the RDE model. Material properties of the domain are updated according to the distribution of cells, and then are used for structural analysis to estimate updated strain field.

Fig. 2

Process of establishing a computational domain that grows over time based on the data from experimental mice aged from embryonic day 13.5 (E13.5) to birth (P0). a µCT and MRM images of the mouse brain and cranial vault at E15.5, E17.5, and P0. b Simplified geometry of the brain and cranial vault, constructed from the µCT and MRM images. The brain shape at E13.5 is assumed as an ellipsoid. c Displacement vectors from points on the surface at the earlier time point to the surface at the later time point are computed using normal mapping method and represented with green arrows. d Schematic diagram of normal mapping. Displacement vectors (green arrows) are normal to the original surface at t₁ at each point on the surface (black asterisks) and end on the target surface at t₂ (red asterisks). The computed displacement vectors are used as boundary conditions for simulation to make the computational domain to grow over time.

Fig. 3

The RDE model predicts the location of primary centers of ossification and pattern of cranial vault bone growth. a Computational result of distribution of concentration of activator relative to inhibitor (a²/h) at E13.5 and E14.5. In superior view, anterior at top; in lateral view, anterior is left. b Computational prediction of distribution of differentiating osteoblasts and cranial vault bone formation by embryonic day. Superior view of skulls, anterior at top, posterior at bottom. Ossification centers for right and left frontal bones appear first (~ E14.5), followed by right and left parietal bones (~E15). Two more ossification centers representing the interparietal bone appear at E15.5. c Observed cranial vault bone formation and growth in embryonic mice (see Appendix A).

Fig. 4

Mechanical strain estimated from structural analysis reveals mechanism of spatio-temporal pattern of bone formation. Computational estimation of volumetric strain rate Ė_V (a) and accumulated volumetric strain E_V (b) on the domain (shown from above; anterior at top, posterior at bottom) arranged by embryonic day.

Fig. 5

Result of parametric study. Various phenotypes including typical crania, prematurely closing sutures, an increase or decrease in the number of bones, and a decrease in size of bones are predicted.

Fig. 6

Comparison of experimental observations and predictions by the RDE model about cranial dysgenesis due to molecular variants. a Experimental observation of vault bones (stained by alizarin red) of a wild-type (WT) and Gdf 6^−/− mouse. Lateral view, rostrum to left, eye at base of coronal suture (CS). CS is open in WT and closed in Gdf 6^−/− mouse (adopted from Clendenning and Mortlock (2012)). b Computational prediction of distribution of osteoblasts with the reference value (left, CS patent) and the reduced value of α_o (right, CS closed). c Experimental observation of alizarin-red stained cranial vault of WT and Axin2^−/− mice (superior view, rostrum at top). Inter-nasal suture between arrows at top; metopic suture (MS) between frontal bones at bottom (adopted from Yu et al. (2005)). d Computational prediction of distribution of osteoblasts, forming bone, and morphology of MS with the reference value (far left), reduced value of α_h, increased value of β_h, and reduced value of γ_h.

Fig. 7

Distribution of osteoblasts at E17.5 from simulation results using various computational models. a Result using a model with only assumptions A1 and A2, without the assumption pertaining to mechanical effects. Bones grow and form sutures through reaction-diffusion process. b Result using a model with assumptions A1, A2, and only the assumption about the mechanical effect on production of activator (A3a). Locations of primary ossification centers are specified by the mechanical effect on the production of activator. c Result using a model with all assumptions A1, A2, and A3. Relative speed of bone growth is achieved by the mechanical effect on cell differentiation (A3b).