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

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.

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.

Introduction
The part of the cranium that surrounds the brain is known as the neurocranium and can be divided into two parts. The cranial base is positioned beneath the brain and is composed of a series of bones that form first as a set of cartilaginous anlagen that are replaced by bone through a process known as endochondral ossification. The cranial vault covers and protects the superior and lateral aspects of the brain and consists of a series of flat bones that form through intramembranous ossification, where undifferentiated mesenchymal cells of either neural crest or mesoderm origin form condensations and differentiate directly into osteoblasts and begin to form bone. 1 Intramembranous ossification is based on epithelialmesenchymal interactions that are controlled by a cascade of reactions of locally produced and circulating factors (e.g., signaling molecules and their receptors, extracellular matrix proteins and receptors) 2-4 that play significant roles in the differentiation of cells into osteoblasts that make bone. 5,6 Though much is known about the details of signaling required for mesenchymal cells to differentiate and begin to produce collagen and then mineralize the collagen matrix, much less is known about how single bones acquire their characteristic shapes and how they cooperate to form an integrated skull. Research in tissue regeneration is revealing that biomechanical input, though not completely understood, plays an important role in these process. [7][8][9][10][11] Several computational studies have been conducted to model the process of bone formation. One common approach is the application of a mathematical model for reaction-diffusion controlled by two interacting chemical molecules, first proposed by Turing 12 , which has been employed in the study of biological pattern formation and the development of biological systems. [13][14][15][16][17] The Turing model, also commonly referred to as the reactiondiffusion model, shows that through the regulatory loop of interacting molecules, the concentration of molecules forms an inhomogeneous spatial pattern. We previously examined the reaction-diffusion model for pattern formation on simple 1D and 2D domains, 18 and used the approach to predict the locations of the primary centers of ossification in the cranial vault and the pattern of cranial bone growth. 19 However, the approach did not include effects of mechanical stimuli such as local stress and strain on cells that have been revealed to have a significant effect on the differentiation of bone cells and on bone formation by several groups. Claes and Heigele 20 reported that intramembranous bone formation occurs within a certain range of strain and stress, otherwise endochondral bone or connective tissue are formed. Guevara et al. 21 proposed a computational approach to find the effect of stress distribution during the development of long bone. Here, we propose a mechanobiological model that combines the influence of mechanical stimuli with a reactiondiffusion model to predict and explain the process of cranial vault bone formation and development in the mouse. We compute the estimated strain distribution across the layer of undifferentiated mesenchyme surrounding the brain that results from rapid expansion of the growing brain. 22 By examining the correlation between the estimated strain and the pattern of bone formation through an iterative comparison of computed and experimental cranial vaults, we present a mathematical model of the biological process of cranial vault bone formation.

Methods
We followed a similar approach that was used in our previous study 19 to describe the molecular activities in the process of the bone formation in the cranial vault. While our previous approach showed that development of cranial bone can be numerically simulated with the reaction-diffusion (R-D) model, details of which will be introduced later (Section 2.2), some parts of the model were based on the experimental observations rather than elucidating the mechanical mechanisms that might underlie the lack of ossification at the apex of the cranial vault dome and how local factors respond to subtle inputs as each bone takes shape. To explore the underlying mechanisms we hypothesize that mechanical stimuli initiated by growth of underlying brain plays a significant role in the process of bone formation. Here we test the hypothesis by conducting a structural analysis (Section 2.1), coupling a molecular reaction-diffusion model with the mechanical input resulting from the structural analysis (Section 2.2) and finally comparing computational results with experimental observations (Section 3).

Structural Analysis
We conducted a structural analysis using the finite volume method to determine how mechanical strain affects the development of cranial vault bone. This requires defining the surface geometry, a continuum model and the boundary conditions. 2.1.1. Geometry-A computational domain was constructed using an isosurface reconstructed from micro CT images of a mouse head at embryonic day 17.5 (E17.5) ( Fig.  1) (pixel size and slice thickness ranged from 0.014 to 0.025 mm). We segmented the five cranial vault bones (right and left frontal bones, right and left parietal bones and a single interparietal bone), made a global surface that covers the cranial vault and then made a thin volume around the surface to represent space occupied by mesenchymal cells. The domain has three surfaces: 1) the inner surface that acts as a surrogate for the surface of underlying brain; 2) the outer surface that represents the layer of mesenchymal cells; and 3) the bottom surface or edge where the inner surface and outer surface meet at the base that we are currently using to represent the chondrocranium, a cartilagenous structure formed prior to the cranial vault. 23 The mesh was generated using a commercial software ANSYS ICEM CFD (ANSYS, Inc., Canonsburg, PA, USA) and contains 51,177 tetrahedral elements.

Continuum
Model-For now, we assume that mesenchymal cells behave as linear elastic, isotropic, homogeneous media. We acknowledge that soft tissue behaves as viscoelastic media in nature but our assumption is used as a first guess to model the relation between biomechanical strain and the pattern of bone formation. Considering the cells as continuum elements, the balance of momentum with no body forces which governs the mechanical behavior of the material is given by: ∂ 2 (ρu)/∂t 2 − ∇ · σ = 0. A constitutive relation between stress and strain for elastic material is described by: σ = 2με+λtr(ε)I. The strain tensor, ε, is defined in terms of displacement u assuming small strain theory by . The equations were solved by the finite volume method using an elastic solid mechanical module in the open source code OpenFOAM. 24 .

Boundary and Initial
Condition-Boundary conditions for three surfaces of the domain for the structural analysis are shown in Fig. 2. A higher pressure of 1 kPa is applied as the boundary condition on the inner surface, representing pressure due to the growth of underlying brain. The value of pressure produced by a growing embryonic brain is unknown in biology but here we adopt a value of intracranial pressure due to brain and cerebrospinous fluid volume as measured in adult mice presented by Oshio et al. 25 . A lower pressure (0 Pa) is adopted for use on the outer surface. The bottom surface is treated as fixed at the moment of analysis.
Material properties were adopted from Pillarisetti et al. 26 who measured mechanical characteristics of mouse embryonic stem cells. Values of 1000 kg/m 3 for density (ρ), 20 kPa for Young's modulus (E) and 0.4 for Poissons ratio (ν) were used for the material properties of the soft tissue. Fig. 3. Figure 3(b) shows the distribution of hydrostatic strain ε hyd which can be estimated by (ε xx + ε yy + ε zz )/3. We focused on the hydrostatic strain based on the study of Carter et al. 27 suggesting that hydrostatic stimuli, rather than shear stimuli, may play an important role in tissue differentiation. When the contour of the strain field ( Fig. 3(b)) and isosurface of the murine cranial vault segmented from micro CT images of an E17.5 mouse ( Fig. 3(a)) are superimposed ( Fig. 3(c)) an obvious correspondence between the distribution of strain and the pattern of bone formation is observed. Less bone forms in areas of high strain. From this observation, we formulate a hypothesis; Hypothesis 1 (H1): Mechanical stimuli produced by growth of the underlying brain play a significant role in the process of bone formation by inhibiting cell differentiation local to areas of high hydrostatic strain.

Results: Strain Field-Results of solid mechanical analysis are shown in
This hypothesis agrees with the study of Claes and Heigele 20 indicating that magnitudes of local stress and strain determine the type of differentiation of mesenchymal cells. Specifically, intramembranous ossification, which is responsible for formation of cranial vault bones, occurs within a certain range of strain and an excessive magnitude of strain restricts bone formation. 6 shows an important role of bone morphogenetic proteins (BMPs) in bone formation and the feedback mechanism between BMPs and their antagonists such as noggin. From this observation we formulate our second hypothesis; Hypothesis 2 (H2): The differentiation of mesenchymal cells is activated by a key molecule (the activator, e.g. BMP) whose behavior is regulated by another molecule (the inhibitor, e.g. noggin) so that the regulatory loop of the two molecules can be modeled using an activator-inhibitor model, a type of Turing's reaction-diffusion model. 12 as used in our previous study. 19 In this study we modified the previously used activatorinhibitor model 28 by adding a strain effect to the expression of key molecules. The coupled reaction-diffusion-strain (R-D-ε) model for pre-patterning is described mathematically as:

Pre-patterning model-The experimental evidence observed by Wan and Cao
where C a and C h indicate concentration of activator and inhibitor, respectively and α, β, γ and D are constant parameters which control various behaviors of the activator and inhibitor system. That is, Equations (1a) and (1b) indicate that the rate of change of concentration of each molecule (∂C a /∂t and ∂C h /∂t) is controlled by its production (α a and α h ε hyd 2 ), degradation (β a C a and β h C h ), interaction between the molecules ( and ) and diffusion in space (D a ∇ 2 C a and D h ∇ 2 C h ). Based on the hypothesis (H1) that high strain inhibits intramembranous bone formation, a term in Eq. (1b) related to the production rate of the inhibitor molecule is given as proportional to the square of the hydrostatic strain (α h ε hyd 2 ) so that an increase in the inhibitor molecule will occur local to regions where the magnitude of strain (either positive or negative) is high. The proportional relation between the order of strain and the production of inhibitor was determined by comparisons of computational results with experimental data. A quadratic term (α h ε hyd 2 ) turned out to be the most appropriate as an α h ε hyd p term does not show a significant effect when p < 2 (e.g., linear term α h ε hyd ) and the term overestimates the effect when p > 2 (e.g., cubic term α h ε hyd 3 ).
It is convenient to nondimensionalize the equations to reduce the number of parameters  where α = (β a α h )/(γ a α a ), β = β h /β a , γ = (α a γ h )/(β a γ a ), D = D h /D a and ∇ *2 = ∂ 2 /∂x *2 . Nondimensional parameters were chosen to satisfy constraints for pattern formation using a linear stability analysis. 29 The values of the parameters for the pattern formation of activator and inhibitor (Eq. (2a) and (2b)) are provided in Table 1.  Table 2.

Bone growth model-Bones
grow from a primary ossification center by the differentiation of mesenchymal cells around the periphery due to a signal input from a morphogen that diffuses and is received by adjacent cells to promote further differentiation. Komori 30 observed the details of this process by finding that runt-related transcription factor 2 (Runx2) is upregulated in osteoblasts and induces the differentiation of mesenschymal cells. From this observation, we formulate our third hypothesis;

Hypothesis 3 (H3):
Bones grow from the primary ossification centers by the action of molecules produced by osteoblasts (e.g., Runx2) that contribute to the formation of extracellular matrix proteins of forming bone and mineralization of the matrix.
differentiated osteoblasts according to the diffusion of the morphogen (H3) and the distribution of strain (H1) (Eq. (4) and (5)): where C o and C m represent concentration of the osteoblasts and the morphogen produced by the osteoblasts, respectively, while subscript i indicates any quantities of the i-th bone region (e.g. i = 1; right frontal, i = 2; left frontal, i = 3; right parietal, i = 4; left parietal and i = 5; interparietal bone region). The variable C o is computed for a short initial time (t < E15) using the initial cell differentiation model (Eq. (3)), and then is split up into the C o-i based on its location. Equation (4) indicates that the mesenchymal cells, exposed to higher concentrations of the osteoblast-produced morphogen whose expression is mediated by Runx2 relative to a threshold value (C mT ), differentiate into osteoblasts. The parameter C oS indicates the saturated concentration of the osteoblast. Adopting the approach of Garzón-Alvarado et al. 16 to model suture formation, the differentiation of mesenchymal cells at the periphery of the i-th bone is modeled to be restricted when it senses higher concentration of the morphogen secreted from the contiguous edge of the j-th bone (i ≠ j) relative to a limitation value (C ml ). The parameter κ is a constant quantifying the amount of osteoblast differentiated by action of the morphogen. The effect of strain is applied to the bone growth model (Eq. (4)) following the procedure for the initial cell differentiation model (Eq. (3)) in the way that mesenchymal cells differentiate when they experience lower strain than a certain limit (ε l ).
Equation (5) models the behavior of the morphogen that is expressed by the osteoblasts and diffuses into space, a behavior that makes the regions of osteoblasts expand over time. The parameter C mS indicates saturated concentration of the morphogen. The parameter ω is a constant quantifying the amount of the morphogen expressed by osteoblasts and D m is the diffusion rate of the morphogen. The values of the parameters for the bone growth model are given in Table 3.
The reaction-diffusion-strain (R-D-ε) model for cranial vault bone formation (Eq. (2)-(4)) is developed from our previously published reaction-diffusion (R-D) model 19 by being modified with mechanical considerations. Differences between the two models are summarized in Table 4 and computational results predicted by the two models will be compared in Section 3 to investigate the effects of mechanical input.
(www.openfoam.org). Since initial conditions may differ from individual to individual, we applied different sets of randomly distributed concentrations of activator (uniform random with 0.5%) to see the influence of our initial perturbation. The initial concentration of inhibitor is uniform over the domain. A zero gradient of concentration of molecules was adopted as the boundary condition for all surfaces, which means molecules are not able to move through the surface. Equation (1)-(3) were solved first until t = E15 to find the locations of the first differentiated osteoblasts (C o > 0.001kg/m 3 ). The C o computed at t = E15 using the initial cell differentiation model (Eq. (3)) is split up into the C o-i based on its location and the C o-i is used as the initial condition to solve Eq. (4) and (5). Figure 4 shows the distribution of activator at initial time (E0) and E15 estimated using two computational models. Looking at the results predicted by our previous reaction-diffusion (R-D) model 19 that does not consider biomechanical input (Eq. (2) with the boxed term replaced by a constant (α const = 1)), variation in the initial perturbation of the activator at the initial time ( Fig. 4(a) and (d)) leads to different spatial patterns of activator at E15 (Fig. 4(b) and (e)). The results show that the distribution of activator predicted by the R-D model depends strongly on the initial condition of the activator. Alternatively, the reactiondiffusion-strain (R-D-ε) model (Eq. (2)) shows that different initial perturbations establish a consistent spatial pattern of activator after some time when reaction-diffusion of interacting molecules and effects of mechanical stimuli are combined. Several (more than 10) cases with different initial perturbations with the model were tested to show a consistent final spatial pattern. The results of only two different cases are shown in Fig. 4(c) and (f). Comparing these results from two models, we conclude that a reaction-diffusion model is sensitive to the location of the initial perturbations of molecules. In nature, the shape of the domain surrounding a brain and the initial perturbations of molecules may vary slightly between individuals within a species but the locations and the pattern of bone formation are similar across individuals. The results from the coupled reaction-diffusion-strain model suggest that local strain can be a factor that guides or contributes to the consistent distribution of key molecules.  Fig. 1 (a)). Figure 6 shows the pattern of cranial vault bone growth predicted by the bone growth model (Eq. (4) and (5)). The results from the model (Fig. 6 (b)) are compared with the result from our previous model 19 (Fig. 6 (a)) that does not consider strain effects (Eq. (4) and (5) with the boxed term removed). The regions of high concentration of differentiated osteoblasts (> 0.001 kg/m 3 ) as represented in the computational results expand from the primary centers of ossification over time revealing areas absent of osteoblasts, that we consider the sutures between bones. The results predicted by the reaction-diffusion model that does not consider mechanical stimuli (Fig. 6 (a)) show that five bones, including the interparietal bone, grow with similar speed and meet at the top of the head at a faster rate than when growth is predicted using the coupled reaction-diffusion-strain model (Fig. 6 (b)). The results from the coupled model show that the growth of the interparietal bone is restricted due to high strain located apical to the bone (Fig. 3), modeling what is observed in experimental animals (Fig.  6 (c)) such that the interparietal bone is relatively smaller than the other bones (i.e., frontal and parietal bones) being constrained superoinferiorly. Growth of the other bones (frontals, parietals) towards the apex of the head is delayed, which is also observed experimentally. The computational results show some key features that accurately reflect our experimental observations. The left and right frontal bones grow faster in the circumferential direction than apically towards the center of the domain. Accelerated growth of the left and right frontal bones toward the parietal bones results in the formation of a suture that mimics the placement and morphology of the coronal suture that forms between the frontal and parietal bones in experimental observations. We present these features of the computational results as evidence of the key role of the distribution of mechanical strain in determining the locations of primary centers of ossification, the spatial-temporal pattern of bone growth, and the formation of the coronal suture (Fig. 3).

Discussion
One major difference between the computational predictions and experimental observations is the large area of undifferentiated mesenchymal cells between the parietal and interparietal bones in the experimental results that is not present in our computational results. Limitations of the current model may account for the observed difference. First, the initial morphology of the domain is based on data from a mouse aged E17.5, a time when ossification centers have formed and bones are in the process of growing. This morphology may differ enough from the morphology of the head at the initiation of ossification center formation to affect certain details of bone growth in our model. Second, the size and shape of the computational domain are fixed over time in this computational model, while the growing mouse head changes in size and shape over time (Fig. 6 (c)). Third, the model of bone growth is based on a baseline strain field in a homogeneous domain of mesenchymal tissue. As bones form in the pattern shown in Fig. 6, the strain distribution would change dramatically, since the new bone would have different material properties from those of mesenchymal cells. We expect the differences between predictions and observations to be reduced in future work by using a domain based on younger embryos, imposing the effect of growth of the domain and calculating strain field in each new time step based on the current distribution of mineralized bone and non mineralized tissue. Another improvement of our model is expected when we implement comparisons of our computational results with experimental observations using statistical approaches for the comparison of biological shapes 31 or the quantitative evaluation of overlapping percentages of bone volume.

Conclusions
We have developed a computational framework to simulate the formation and growth of cranial vault bones in the embryonic mouse. We conducted a structural analysis to explore the mechanical effects on bone formation in the cranial vault. The results of the structural analysis show that the distribution of mechanical strain over the cranial surface has a strong effect on the pattern of bone formation such that bone does not form in areas of high strain. Based on the hypothesis that high strain suppresses differentiation of mesenchymal cells, the previously published reaction-diffusion model 19 Table 1 Model parameters in coupled reaction-diffusion-strain model for pattern formation of activator and inhibitor (Eq. (2a) and (2b)).

Non-dimensional parameters Value
Production α  Table 2 Model parameters in coupled reaction-diffusion-strain model for primary centers of ossification (Eq. (3) Table 3 Model parameters in coupled reaction-diffusion-strain model for bone growth (Eq. (4) and (5)).

Parameters Value Units
Growth Constant κ