A COMPUTATIONAL ANALYSIS OF BONE FORMATION IN THE CRANIAL VAULT USING A COUPLED REACTION-DIFFUSION-STRAIN MODEL

Abstract

Bones of the murine cranial vault are formed by differentiation of mesenchymal cells into osteoblasts, a process that is primarily understood to be controlled by a cascade of reactions between extracellular molecules and cells. We assume that the process can be modeled using Turing's reaction-diffusion equations, a mathematical model describing the pattern formation controlled by two interacting molecules (activator and inhibitor). In addition to the processes modeled by reaction-diffusion equations, we hypothesize that mechanical stimuli of the cells due to growth of the underlying brain contribute significantly to the process of cell differentiation in cranial vault development. Structural analysis of the surface of the brain was conducted to explore the effects of the mechanical strain on bone formation. We propose a mechanobiological model for the formation of cranial vault bones by coupling the reaction-diffusion model with structural mechanics. The mathematical formulation was solved using the finite volume method. The computational domain and model parameters are determined using a large collection of experimental data that provide precise three dimensional (3D) measures of murine cranial geometry and cranial vault bone formation for specific embryonic time points. The results of this study suggest that mechanical strain contributes information to specific aspects of bone formation. Our mechanobiological model predicts some key features of cranial vault bone formation that were verified by experimental observations including the relative location of ossification centers of individual vault bones, the pattern of cranial vault bone growth over time, and the position of cranial vault sutures.

Keywords: computational morphogenesis; developmental biology; finite volume method; intramembranous ossification; mechanobiology; skull growth.

Figures

Fig. 1

Computational Domain: A computational domain is constructed using 3D reconstruction of micro CT images of skull of mouse at embryonic day 17.5 (E17.5). Bones of the cranial vault (two frontal, two parietal and a single inter-parietal bone) were extracted to establish a global surface and a 3D space around the surface is constructed and meshed (shown in blue). (a) Superior view, anterior at top, posterior at bottom. (b) Left lateral view, facial skeleton to the left.

Fig. 2

Schematic showing boundary conditions for solid mechanical analysis. Difference between higher pressure of 1 kPa on the inner surface and lower pressure of 0 Pa on the outer surface represents pressure coming from growth of the underlying brain. The bottom surface, which represents the superior edge of the chondrocranium that is formed earlier than cranial vault bones, is treated as fixed at the moment of analysis. (a) isometric view (b) superior view and (c) section view of the domain.

Fig. 3

Distribution of hydrostatic strain (εhyd) and comparison with experimental data. (a) isosurface of mineralized cranial vault bone segmented from micro CT image of a mouse at E17.5. Superior view, anterior at top, posterior at bottom. (b) Field of hydrostatic strain. High tensile strain appears on the top of the domain. (c) Overlap of isosurface of cranial vault bone segmented from micro CT image and strain field. The locations of bones correspond with regions of low strain.

Fig. 4

Distribution of activator. (a) and (d) - two different sets of randomly distributed activator with 1% of perturbation at initial time. (b) and (e) - distribution of activator at E15 predicted by reaction-diffusion (R-D) model without consideration of strain effect for two different initial conditions. Different initial perturbations lead to different patterns of activator distribution. (c) and (f) - distribution of activator at E15 predicted by coupled reaction-diffusion-strain (R-D-ε) model for two different initial condition. Different initial perturbations lead to a consistent pattern of activator distribution.

Fig. 5

Primary centers of ossification (a) superior view (b) lateral view. The region of high concentration of osteoblasts (> 0.001 kg/m3) is shown opaque but elsewhere are transparent. Five locations of high concentration of osteoblast agree well with experimental observation (two frontal, two parietal and a single interparietal bone).

Fig. 6

Bone growth. Change of region of high osteoblast concentration over time predicted by (a) reaction-diffusion model without consideration of mechanical strain and (b) coupled reaction-diffusion-strain model. The region of high concentration of osteoblasts (>0.001 kg/m3) is represented in the computational results. The regions of high concentration of osteoblast expand from primary centers of ossification over time resulting in sutures between bones. Computational results are compared with (c) experimental observations of the pattern of cranial vault bone growth in mouse as visualized by micro CT.