A Computational Analysis of Bone Formation in the Cranial Vault in the Mouse

Bones of the cranial vault are formed by the differentiation of mesenchymal cells into osteoblasts on a surface that surrounds the brain, eventually forming mineralized bone. Signaling pathways causative for cell differentiation include the actions of extracellular proteins driven by information from genes. We assume that the interaction of cells and extracellular molecules, which are associated with cell differentiation, can be modeled using Turing’s reaction–diffusion model, a mathematical model for pattern formation controlled by two interacting molecules (activator and inhibitor). In this study, we hypothesize that regions of high concentration of an activator develop into primary centers of ossification, the earliest sites of cranial vault bone. In addition to the Turing model, we use another diffusion equation to model a morphogen (potentially the same as the morphogen associated with formation of ossification centers) associated with bone growth. These mathematical models were solved using the finite volume method. The computational domain and model parameters are determined using a large collection of experimental data showing skull bone formation in mouse at different embryonic days in mice carrying disease causing mutations and their unaffected littermates. The results show that the relative locations of the five ossification centers that form in our model occur at the same position as those identified in experimental data. As bone grows from these ossification centers, sutures form between the bones.


INTRODUCTION
Bones of the mammalian cranial vault are formed by the process of intramembranous ossification where condensations of multipotent mesenchymal cells differentiate directly into functioning osteoblasts to form bone (Tubbs et al., 2012;Percival and Richtsmeier, 2013). Figure 1 shows the estimated length and time scales associated with various processes of skull development in the mouse. The study of the formation of cranial vault bones is critical because it may help uncover fundamental mechanisms associated with birth defects such as craniosynostosis (involving premature closure of cranial vault sutures and skull dysmorphology) as well as answer basic questions in developmental biology and evolution.
Initially, undifferentiated mesenchymal cells migrate to future sites of bone formation situated on the brain and its meningeal layers (pia, arachnoid, dura mater). Over time, the proliferation and differentiation are regulated by growth factor signaling pathways and their downstream transcription factors in order for these cells to become committed to an array of different fates. In some cells, the intracellular signaling pathways cause differentiation of mesenchymal cells into osteoblasts, the cells that build bone Gordeladze et al., 2010;Tubbs et al., 2012). Differentiation of condensed groups of osteoblasts results in the formation of ossification centers that form in tissue membranes surrounding the brain. Next, osteoblasts begin to secrete a bone matrix, osteoid, which is then mineralized, eventually forming a bone of the cranial vault. For the mouse, the process of cell migration begins around embryonic day 9 (E9) and skull bone growth continues postnatally, involving length scales ranging from the nanometer to millimeter as schematically shown in Figure 1.
In order to understand the fundamental mechanisms of skull growth, both experimental and computational methods have been employed. Many studies have experimentally studied the roles of various proteins in cranial bone formation (Holleville et al., 2003;Wan and Cao, 2005) and growth of cranial bones (Martínez-Abadías et al., 2013;Motch Perrine et al., 2014;Percival and Richtsmeier, 2013) and compared the craniofacial bone growth patterns of normal mice and those carrying mutations that in humans cause disease. Figure 2 shows that cranial vault bone elements (frontal, parietal, and interparietal bones) appear from embryonic day 14.5 (E14.5) and continue to grow through postnatal day 0 (P0) and beyond. It also shows that the frontal bone forms first and the interparietal bone forms later, forming sutures between individual bones while they grow. Experimental studies are extremely valuable but can be costly and only so many variations can be explored. Therefore, it is valuable to also examine the possibility of using computational methods to understand fundamental mechanisms of morphogenesis.
Several computational studies have been conducted to model the process of skull bone formation. A mathematical model for reaction-diffusion controlled by two interacting chemical molecules, proposed by Turing (1952), has been employed in the study of biological pattern formation and development of biological systems. Kondo and Shirota (2009) analyzed the mechanism of skin pattern formation of animals using the Turing model and (Marcon and Sharpe, 2012) adopted the model to explain various biological development processes. Garzón-Alvarado et al. (2013) used the model to establish a computational framework for investigating bone formation in human cranial vault. The model, commonly referred to as the reaction-diffusion model, shows that through the regulatory loop of interacting molecules the concentration of the molecules forms an inhomogeneous special pattern in space. In this study, we adopt an approach similar to that of Garzón-Alvarado et al. (2013), to study growth of the skull in a mouse model of human disease and then propose an extension of the framework.
As depicted on Figure 1, we subdivide the process into two stages: (1) initiation (differentiation) of primary centers of ossification; and (2) bone growth. In the first stage, we focus on differentiation of osteoblast lineage cells (OLCs), which leads to the initial primary centers of ossification of the flat bones of the cranial vault. We assume that the interaction of extracellular molecules, which are associated with the differentiation process of OLCs along an osteogenic path, can be modeled using the reaction-diffusion model. Reaction-diffusion models can be further subdivided into activator-inhibitor and activator-substrate models according to how molecules interact with each other (Gierer and Meinhardt, 1972). The primary difference between the two models is the way one molecule inhibits the other molecule. In the activatorinhibitor model, one molecule (activator) enhances the other molecule (inhibitor) but the enhancement of inhibitor inhibits the action of activator so that the molecules are in phase. On the other hand, in the activator-substrate model, one molecule (activator) consumes the other molecule (substrate) to be enhanced and eventually is restricted by depletion of the other molecule (substrate) so that the molecules are not in phase. More details about these models can be seen in a work of Koch and Meinhardt (1994). Since many details remain to be discovered in order to fully understand the key molecular players in proliferation and differentiation of OLCs as they form intramembranous bones of the cranial vault, it is currently not known whether the activatorinhibitor or the activator-substrate approach more accurately models intramembranous bone formation. Since the activatorinhibitor model more closely models what has been observed experimentally and because the regulatory relation in the model is simpler, unlike Garzón-Alvarado et al. (2013), we employ the activator-inhibitor model.
In the second stage, we deal with the rapid proliferation of bone cells from the primary condensations and their outward migration Frontiers in Bioengineering and Biotechnology | Biomechanics (stage 1) to form bones of the cranial vault, particularly the frontal and parietal bones. Work from other laboratories collectively suggest a pattern of presumptive bone cells expanding outward from mesenchymal condensations, predominantly towards the apex of the head (Iseki et al., 1997;Rice et al., 2000;Ting et al., 2009). Under normal regulatory conditions, these expanding condensations represent the primary region of osteoblast differentiation, the source of differentiating cells, and locus of osteogenesis and bone growth. They also define the earliest shapes of developing bones. We modeled the action of the morphogen stimulating bone growth from these osteoblasts.
In this paper, we present a combined experimental and computational study of cranial bone growth. In Section 2, we present the computational framework for determining the location of primary ossification centers and simulating growth of individual bones of the cranial vault. Our computational predictions are guided and compared to a large archive of experimental data. Micro computed tomographic (µCT ) images of mouse heads of various embryonic ages were used for making a computational domain and testing and tuning the mathematical model by comparison between experimental and computational results. Section 3 gives results of the simulation and discussions about them are made in Section 4. Finally, in Section 5, we conclude and summarize our findings and provide suggestions for future work.

Regulation of osteoblast differentiation
Osteoblasts differentiate from mesenchymal progenitors, going through distinct developmental stages, which are regulated by various developmental signals. Although the role of key developmental signals is well defined, most of this knowledge comes from studies of endochondral development of long bones and even there, little is known about how the signals execute the osteoblast-specific differentiation program (Long, 2011).
This study is based on the hypothesis that differentiation of mesenchymal cells takes place by chemical reactions triggered by the interaction of extracellular molecules (both local and systemic regulatory signals), with cell surface ligands, although it is important to point out that the specific proteins in the mouse system are not clearly determined. Thus, we quickly arrive at one challenge of this approach: determining quantitative values of the parameters that represent specific proteins. We acknowledge that parameters are the cause of great speculation throughout this study and maintain that developing this model as part of a collaborative project that seeks to precisely define the role of the critical molecular participants will improve our model over time. Conversely, it is clearly a challenge to obtain and measure these parameters experimentally, and we propose that computational studies presented here will help elucidate possible bounds of parameters and guide the design of future experiments to determine their values with increased certainty. To mitigate the speculative nature of parameterization, we conduct a stability analysis that helps to bound some values and parametric studies to estimate the effects of some of the parameters. Figure 3 is a schematic diagram showing the process of differentiation of mesenchymal cells to osteoblasts and the regulatory relationship between extracellular molecules, which influence cell differentiation. In our model, mesenchymal cells surrounding the brain express diffusible extracellular molecules, which play a key role in cell differentiation and the molecules are assumed to follow behaviors that can be described using the activator-inhibitor www.frontiersin.org FIGURE 3 | Schematic diagram of extracellular and cellular process associated with differentiation of mesenchymal cells to osteoblast cells. Undifferentiated mesenchymal cells surrounding the brain express diffusible extracellular molecules, which play a key role in cell differentiation ( and ). One of the molecules (activator) activates signaling pathways to initiate cell differentiation of mesenchymal cells into osteoblasts ( ). In an extracellular process, the activator simultaneously enhances itself ( ) and the other key molecule (inhibitor) ( ), while the inhibitor inhibits the activator ( ). These two proteins eventually establish a regulatory loop and diffuse in space with different speed ( ) to form an inhomogeneous spatial pattern of concentration. model (Figure 3, and ). The activator initiates chemical reactions that activate signaling pathways to initiate differentiation of mesenchymal cells into osteoblasts (Figure 3, ). In an extracellular process, the activator simultaneously enhances itself (Figure 3, ) and the inhibitor (Figure 3, ), while the inhibitor inhibits the activator (Figure 3, ). These two proteins eventually establish a regulatory loop and the diffusion of products of the inhibitor and activator (Figure 3, ) form an inhomogeneous spatial pattern of concentration that results in molecules condensing at specific locations in a domain. The features of an activator-inhibitor system can be mathematically modeled using a reaction-diffusion model, which is given by: where, C a and C h represent the concentration of activator and inhibitor, respectively. σ a and σ h are the constants quantifying the production of activator and inhibitor from mesenchymal cells. The parameters µ a and µ h quantify degradation or depletion of the proteins. The parameters ρ a and ρ h are constants associated with the non-linear interaction between activator and inhibitor. The interaction term shows that the activator enhances itself and the inhibitor [C a 2 in numerator in equations (1a) and (1b)] but is constrained by the inhibitor [C h in denominator in equation (1a)]. D a and D h represent the diffusion rate of each molecule and 2 is the Laplace operator i.e., 3 i=1 ∂ 2 ∂x i 2 describing the spatial diffusion of molecules. So equations (1a) and (1b)  In this model, parameters should satisfy a certain constraint in order to make an inhomogeneous spatial pattern from a very small perturbation on a homogeneous initial condition. If diffusion of a molecule is fast relative to the reaction between activator and inhibitor, a small perturbation cannot be amplified but the molecules will reach another homogeneous condition. In the future, actual values of parameters might be defined by way of laboratory experiments but experiments that quantify these types of parameters are currently limited. Therefore, parameters should be estimated with careful consideration of the biologically reasonable range. Additionally, Koch and Meinhardt (1994) suggested constraints that the parameters should satisfy for pattern formation using a linear stability analysis. Homogeneous, steady state initial concentration of each molecule can be achieved mathematically by setting time rate of change and spatial diffusion terms in equations (1a) and (1b) be to zero: where C a0 and C h0 signify the concentration of activator and inhibitor, respectively, at steady state. By adding a small perturbation to the homogeneous steady condition, the concentration of two molecules can be represented as equations (3a) and (3b) where, In equation (3b) perturbation δC a and δC h are assumed to change in time and space. The value ω can be a complex number and the imaginary part of it represents a frequency at which the perturbation changes in time. The change in space is characterized by wave number k, which is associated with ω. When the real part of the frequency ω is positive, the perturbation increases with time Frontiers in Bioengineering and Biotechnology | Biomechanics so that concentrations of molecules can form an inhomogeneous spatial pattern. By substituting equations (2), (3a), and (3b) into equations (1a) and (1b) and conducting a linear stability analysis, a condition, which parameters should satisfy for pattern formation (i.e., for making the real part of ω to be positive) can be obtained as below: Details can be found in Koch and Meinhardt (1994). In this study, we used a set of parameters, which meet the condition of equation (4) as listed in Table 1. Using values in Table 1 with equation (4) yields: 0.667 < 2.0 < 8.468. The effects of various model parameters on the results can be seen in Supplementary Material. Similar to Garzón-Alvarado et al. (2013), our model is based on the hypothesis that differentiation of mesenchymal cells into osteoblasts is triggered by an activator. Here, we acknowledge that there are several biological stages prior to the final differentiation of OLCs into osteoblasts (Figure 4), so our model focuses on cellular mechanisms occurring after the concentration of OLCs has formed and the cells have started along an osteoblast differentiation trajectory. These cells are more likely to form osteoblasts where concentration of activator is high. In addition to concentration of activator, spatial positions of cells are also found to affect the rate of cell differentiation. Tubbs et al. (2012) shows that cells on the inferior surface (that will form bones of the skull base) differentiate faster than cells on the superior surface (cranial vault bones). With these impacting factors on cell differentiation the rate of generation of osteoblast can be modeled by the equations (5a) and (5b).
C a n C a n + C aT n T a n T a n +t n f (x rel ) In the equations, C o represents the concentration of osteoblast and C a indicates the concentration of activator as before. C aT represents the threshold concentration of activator that allows mesenchymal cells to differentiate, which means only cells in the region where the concentration of activator exceeds the threshold value can differentiate into osteoblasts.

Bone growth
After the formation of primary ossification centers, essentially condensations of osteoblasts differentiated from mesenchymal cells, bone grows from these centers by a combination of the production of osteoid and its mineralization by osteoblasts and continued differentiation of rapidly proliferating mesenchymal cells that migrate outwards from the condensation (Iseki et al., 1997;Rice et al., 2000;Yoshida et al., 2008;Ting et al., 2009). The process can be shown schematically in Figure 6. Osteoblasts in the ossification centers begin to secret various morphogens. Some of these morphogens (e.g., Dlx5) are produced from osteoblasts and diffuse into the neighboring space allowing adjacent mesenchymal cells to differentiate into osteoblasts. Experimental data reveal in mice that five primary ossification centers, one each for the bones of the cranial vault (right and left frontal, right and left parietal and interparietal bones) that appear first on the more rostral surface of the embryonic mouse brain as shown in Figure 2. The squamous portion of the occipital bone develops later. In normal conditions, the bones of the skull vault are separated by sutures, fibrous joints that accommodate the expanding brain and allow the skull to undergo reshaping during birth (Ting et al., 2009). Although the mechanisms of suture formation are only partly known, most research favors molecular mechanisms that balance the proliferation and differentiation of osteogenic cells in the developing sutures (Lee et al., 1999;Yousfi et al., 2001Yousfi et al., , 2002Chen et al., 2003;Bialek et al., 2004) and/or that establish tissue boundaries (Merrill et al., 2006;Ting et al., 2009) in suture formation and premature closure of sutures. Similar to Garzón-Alvarado et al. (2013) in this study, suture formation between bone elements is modeled by assuming that differentiation of mesenchymal cells adjacent to bone is restricted by the morphogen secreted from the adjacent bone elements, although we realize that the signaling could come from, or be transmitted through dura mater surrounding the developing brain (Opperman, 2000;Richtsmeier and Flaherty, 2013). The differentiation process in the bone growth stage is modeled by following equation in this study.
In this equation, C o indicates the concentration of osteoblast and subscription i indicates any quantities of i-th bone where i can represent the number of different regions of skull bone (in our case, 5 regions: right and left frontal, right and left parietal and interparietal bone). C os is saturation concentration of osteoblast that sets the upper limit of generation of osteoblast. C m represents the concentration of the morphogen, which is secreted from osteoblasts and allows differentiation of adjacent mesenchymal cells. C mT is the threshold value of the morphogen that triggers differentiation of mesenchymal cells. C ml in equation (6) represents the limitation concentration of the morphogen, in which cell differentiation of i-th bone is restricted when it senses morphogen above this value secreted from j-th bone (i = j). λ is a constant quantifying the amount of osteoblast differentiated by action of the morphogen. n is a constant characterizing how sharply the rate of cell differentiation increases after the concentration of the morphogen reaches the threshold value, or decreases after concentration of osteoblast and morphogen from other regions of bone elements reach a saturation and limitation value, respectively. So equation (6)  In our model, the production and diffusion of the morphogen associated with bone growth are modeled by another diffusion equation: In this equation, C m represents the concentration of the morphogen. C o indicates the concentration of osteoblasts from which the morphogen is expressed. C mS is the saturation concentration of the morphogen, in which, production of the morphogen is FIGURE 6 | Schematic diagram of bone growth from primary ossification centers. Osteoblasts in the ossification centers begin to secret a morphogen, which diffuses into the neighboring space and allows adjacent mesenchymal cells to differentiate into osteoblasts. Suture formation between bone elements is modeled by assuming that differentiation of mesenchymal cells adjacent to bone is restricted by the morphogen secreted from the adjacent bone elements. Morphogen limitation C ml 1.0 × 10 −3 ng/mm 3 restricted after this value. D m represents diffusion rate of the morphogen and 2 is the Laplace operator describing spatial diffusion of the morphogen. α is a constant quantifying the amount of the morphogen expressed by osteoblasts. n is a constant characterizing how sharply the expression of the morphogen decreases after the concentration of the morphogen reaches the saturation value. So equation (6) shows that the time rate of change of the morphogen concentration [equation (7)-] is affected by concentration of already differentiated osteoblasts [equation (7)-], which secrete the morphogen, restriction when it reaches saturation [equation (7)-] and diffusion into space occupied by undifferentiated mesenchymal cells [equation (7)-]. The model parameters used with equations (6) and (7) are listed in Table 3.

Computational domain
In our previous study, the reaction-diffusion model of activator and inhibitor (BMP2 and Noggin in the paper) was solved in a one-dimensional and two-dimensional domain (Lee et al., 2014).
In this study, we used three-dimensional reconstructions of micro CT images of mouse heads at embryonic day 17.5 (E17.5) to conduct a computational simulation on a more realistic domain. At E17.5, some parts of mouse skull are already formed as shown in Figure 7. We considered only bones of the cranial vault and simplified them to make a surface on which reaction-diffusion of extracellular molecules takes place. Then we made a shellshaped domain around the surface and generated a mesh using commercial software ANSYS ICEM. The length of the vault is approximately 7 mm and the thickness of the shell is about 0.3 mm. In total, 84,200 elements are used for the calculation. In future applications, we will move these simulations to surfaces of earlier embryonic ages before bone is formed to more realistically model the domain prior to the initiation of cranial vault bone formation.

Initial condition
From a computational point of view, initial conditions are a significant aspect to consider in a simulation. The model in this study also depends on initial concentration of activator and inhibitor molecules. Different initial concentrations of activator and inhibitor result in significant changes in the final bone growth patterns. The effect of various initial conditions can be seen in Supplementary Material. Experimental research reveals that mesenchymal condensations that will form the primary centers of ossification initiate at the supra orbital region just above the globe of the eye and develop into right and left frontal bones (Tubbs et al., 2012). Based on these observations, we initiated a small perturbation of activator (0.5% of steady state condition) only at two points corresponding to the location of the right and left frontal bones on the domain. Consequently, the entire surrounding domain has steady state concentration of activator and inhibitor [C a0 = 1.20 ng/mm 3 , C h0 = 1.44 ng/mm 3 in equation (3a)] except the two points as shown in Figure 8. Zero-gradient of concentration of molecules was applied as boundary condition at domain boundaries, which means there is no flux of molecules by diffusion across the boundaries. www.frontiersin.org

FIGURE 8 | Initial condition.
A small perturbation of activator on two points on the frontal side of domain was given as an initial condition.
Elsewhere the domain has homogeneous concentration of activator and inhibitor. This condition is based on the biological observation showing that mesenchymal condensations initiate at the supra orbital region just above the globe of the eye and develop into right and left frontal bones.

PRIMARY OSSIFICATION CENTERS
The reaction-diffusion model [equations (1a) and (1b)] and the cell differentiation model [equations (5a) and (5b)] were solved using the finite volume method and changes in concentrations of activator, inhibitor, and osteoblasts in the domain were examined over time. We assume mesenchymal cells are present near E0. This assumption is not the case in real development because there is neither brain formed nor mesenchymal cells surrounding the brain, but helps to establish a modeling process. Simulation using a brain shape at earlier time is planned in a subsequent research. Figure 9 shows the concentrations of activator and inhibitor at time 0 and after E15.4. At time E0 activator and inhibitor are distributed uniformly through the domain with only 0.5% of disturbances of activator on the two points at frontal side. These small perturbations at time E0 increase over time by regulatory reaction between activator and inhibitor to make a specific pattern of concentration in the domain so that six highly concentrated regions appear as shown in result at E15.4. Concentrations of activator and inhibitor are in phase and reveal the key feature of activatorinhibitor models: the activator is enhanced by itself and activates the inhibitor, but the inhibitor acts only to inhibit the activator. Based on our hypothesis that differentiation of mesenchymal cells to osteoblasts will occur local to points of elevated concentration of activator, our cell differentiation model [equations (5a) and (5b)] predicts five primary centers of ossification as shown in Figure 10. Figure 10 shows the concentrations of osteoblasts at E0 and E15.4. We consider the regions where the concentration of osteoblast exceeds a threshold value (1.0 ng/mm 3 ) to develop into a bony element via the formation of primary ossification centers and these regions are shown as opaque while other regions are translucent in Figure 10. Compared to the distribution of the activator (Figure 9), the primary ossification centers appear where the concentration of the activator is high except the top of the domain. The top of the domain has high concentration of the activator but due to spatial effect on cell differentiation, which we assume in the model [governed by equations (5a) and (5b)] ossification does not initiate in this region. Biologically, this may be an important observation as this is the location where the anterior fontanelle (soft spot) forms in developing heads. Five primary ossification centers form and are located on the right and left side of the anterior region, on the right and left side of the middle region, and in the center of the rear side of the domain, respectively. These computationally generated centers of ossification correspond with the biological location of the right and left frontal bones, right and left parietal bones, and interparietal bone, respectively.

BONE GROWTH
The bone growth model [equations (6) and (7)] was solved using the results of primary ossification centers as an initial condition and the change of concentration of differentiated osteoblast was examined. Figure 11 shows only regions where concentration of osteoblasts exceeds the threshold value (1.0 ng/mm 3 ), which can develop into bone eventually. It shows the regions of high concentration of osteoblasts expand from the primary ossification centers from E15.4 to E22.9. The region of condensed osteoblasts on the posterior end of the domain grows slow compared to the other four regions because this primary ossification center is the last to form and is smaller in size. Spaces between the ossifying regions become reduced as the bones grow but they do not meet each other because of the repulsive effect that occurs between bones in the model. This is an important aspect of our model, which may model suture formation.

DISCUSSION
The number and locations of ossification centers found in the simulation results (Figure 10) agree well with experimental data showing five ossification centers representing the two frontal, the two parietal, and the single interparietal bones of the cranial vault (Figures 2 and 7B). Based on the fact that these results were computed by not specifying any locations, which would eventually become ossification centers, but by only solving the diffusionreaction model with proper initial perturbations and model parameters, the computational framework using Turing's model can be seen as one reasonable approach for exploring processes required for the formation of primary ossification centers and for cranial vault growth.
The results of simulation of bone growth (Figure 11) show that the region of osteoblasts expands from the ossification centers. Our bone growth model explains this process as differentiation of adjacent undifferentiated mesenchymal cells near primary centers of ossification into osteoblasts by the action of a morphogen that is expressed from pre-existing osteoblasts and diffused in space with FIGURE 11 | Change of region of high concentration of osteoblast over time. The regions originally marked by the differentiation of osteoblasts expand from the primary centers of ossification over time. Because these osteoblasts differentiate into osteocytes and eventually become trapped within mineralized bone, it can be said to show pattern of bone growth. The results agree well with experimental observation showing two frontal bones, two parietal bones, and one interparietal bone. Sutures form between bones as bones grow according to repulsive effect between bones in our model.

www.frontiersin.org
time. This is in agreement with previous research that demonstrates that the cells that add to the primary ossification centers by differentiating into osteoblasts come from the condensations rather than being recruited from other mesenchymal populations surrounding the brain (Yoshida et al., 2008;Ting et al., 2009). Once differentiated, osteoblasts function to produce osteoid along collagen bundles and then mineralize that matrix eventually becoming trapped within mineralized lacunae and differentiating into osteocytes. Here, we have not modeled the entrapment of osteocytes into a mineralized matrix, but the expansion of regions of high concentration of osteoblasts is considered as bone growth. Results shown at E20.4 (Figure 11) are representative of data at P0 because gestation period of mouse, although variable, is approximately 20 days. The simulation result shows smaller volume of bones compared to experimental data (Figure 2), which means a slower growth rate. In addition, the relative size and shapes of the individual bones are somewhat different from experimental observations. These differences can be overcome by subtle adjustment of model parameters and, importantly, the adjustments may add to our knowledge of the processes involved in the formation of ossification centers and their growth. The effects of different model parameters on the simulation are presented in Supplementary Material. Our hypothesis that morphogens from one bone element inhibits the growth of contiguous bones [equation (6)], may not be true, but can be tested experimentally. That our results show suture formation similar to what is observed experimentally suggests that suture formation involves some kind of mechanism of repulsive factors between bones, and this does not counter a hypothesis of boundary formation between cellular compartments that serve as signaling interfaces as suggested by others (Merrill et al., 2006;Ting et al., 2009). Here, we suggest that the factors driving the formation of sutures may include chemical substances or mechanical stimuli by the growing brain or by opposing bones as the gaps between them narrow. Since the processes that control normal suture development and the mechanisms underlying abnormal premature suture closure are not well understood, further research into the effects of mechanical stress on bone growth is warranted.
Although the simulation captures many features of the developmental process, it has some limitations as well. Key molecules play important roles in cell differentiation and their identity, as well as the real values of these parameters remain to be determined. We expect this to be achieved in a follow-up study that will elucidate the cellular-level changes that occur in cranial development providing the basis for joining molecular cues and cell behavior with 3D shape changes that occur during ontogeny. The number of cells in initial ossification centers, rate of OLC differentiation and proliferation, intracranial pressure gradients from growth induced skull-soft tissue interaction, and rate of suture closure can be parameterized and modified in the model. The results can be continually quantitatively compared to our extensive image archive of developing cranial soft and hard tissues.
In this study, the computational domain is fixed in both size and shape while in biological systems the surface on which bone development occurs expands and changes shape as the underlying brain grows. Indeed, the size and shape of the domain does affect the number and location of ossification centers because diffusion is strongly related with geometry of the domain. In the future, we plan to let the shape of the computational domain change (or "grow") over time. This would offer a more realistic model, but we realize that our diffusional parameters will need to change over time as well. Some parts of our model are based on the experimental observation rather than clear mechanical mechanisms of how sutures form, why ossification does not occur at the top of the domain and how subtle, complicated shape of each bone forms. We expect that mechanical stimuli such as pressure and stress may affect these phenomena and they will be a focus of a future study. In a forthcoming study, our proposed framework will be improved by including effects of brain growth and mechanical stimuli (e.g., stress) between the brain and growing bone and pressure distribution in head.

CONCLUSION
Growth of the cranial vault is coordinated through tissue-tissue interactions between the brain, the developing meninges, the bone primordia, and the cranial sutures (Han et al., 2007), but the specifics are not well defined. Here, the processes associated with formation of primary centers of ossification and bone growth were mathematically modeled and solved using the finite volume method. The results show that five primary ossification centers form at positions like those identified in experimental data (Figure 10). Our results reveal bone growing from the primary ossification centers forming sutures between bones (Figure 11). Our study shows that the development of the cranial vault can be numerically simulated using the established computational framework. We expect that changes in model parameters, when examined in parallel with laboratory experimentation, will help clarify some of the key players and mechanisms of skull development, both normal and pathological.