Trabecular bone fracture healing simulation with finite element analysis and fuzzy logic

Trabecular bone fractures heal through intramembraneous ossification. This process differs from diaphyseal fracture healing in that the trabecular marrow provides a rich vascular supply to the healing bone, there is very little callus formation, woven bone forms directly without a cartilage intermediary, and the woven bone is remodelled to form trabecular bone. Previous studies have used numerical methods to simulate diaphyseal fracture healing or bone remodelling, however not trabecular fracture healing, which involves both tissue differentiation and trabecular formation. The objective of this study was to determine if intramembraneous bone formation and remodelling during trabecular bone fracture healing could be simulated using the same mechanobiological principles as those proposed for diaphyseal fracture healing. Using finite element analysis and the fuzzy logic for diaphyseal healing, the model simulated formation of woven bone in the fracture gap and subsequent remodelling of the bone to form trabecular bone. We also demonstrated that the trabecular structure is dependent on the applied loading conditions. A single model that can simulate bone healing and remodelling may prove to be a useful tool in predicting musculoskeletal tissue differentiation in different vascular and mechanical environments.     

A 3D computational simulation of fracture callus formation: influence of the stiffness of the external fixator

The stiffness of the external fixation highly influences the fracture healing pattern. In this work we study this aspect by means of a finite element model of a simple transverse mid-diaphyseal fracture of an ovine metatarsus fixed with a bilateral external fixator. In order to simulate the regenerative process, a previously developed mechanobiological model of bone fracture healing was implemented in three dimensions. This model is able to simulate tissue differentiation, bone regeneration, and callus growth. A physiological load of 500 N was applied and three different stiffnesses of the external fixator were simulated (2300, 1725, and 1150 N/mm). The interfragmentary strain and load sharing mechanism between bone and the external fixator were compared to those recorded in previous experimental works. The effects of the stiffness on the callus shape and tissue distributions in the fracture site were also analyzed. We predicted that a lower stiffness of the fixator delays fracture healing and causes a larger callus, in correspondence to well-documented clinical observations.     

The influence of mechanical stimulus on the pattern of tissue differentiation in a long bone fracture--an FEM study

2D, coronal plane, finite elements models (FEMs) were developed from orthogonal radiographs of a diaphyseal tibial fracture and its reparative tissue at four different time points during healing. Each callus was separated into regions of common tissue histology by computerised radiographic analysis. Starting point values of tissue material properties from the literature were refined by the model to simulate exactly the mechanical behaviour of the subject's callus and bone during loading. This was achieved by matching measured inter-fragmentary displacements with calculated inter-fragmentary forces. Stress and strain distributions in the callus and bone were calculated from peak inter-fragmentary displacements measured during natural walking activity, and were correlated with the subsequently observed pattern of tissue differentiation and maturation of the callus. The growth and stiffening of the external callus progressively reduced the inter-fragmentary gap strain. Partial maturation of the gap tissue was apparent only one week before fixator removal. Principal stresses in the callus were compared with 'yield stresses' in corresponding tissue from the literature. This indicated the presence of stress concentrations medial and lateral to the fracture gap, which probably caused tissue damage during normal activity levels. Tissue damage may also have precipitated partial structural failure of the callus, both of which were believed to have delayed healing during the middle third of the fixation period. Had the fixation device provided greater inter-fragmentary support during early healing, this may have prevented callus failure and the consequent delay in healing. A further benefit of this would have been the reduction of the initially high intra-gap tissue strains to a magnitude more conducive to earlier maturation of the bridging tissue that united the bone.     

Comparison of biophysical stimuli for mechano-regulation of tissue differentiation during fracture healing

Most long-bone fractures heal through indirect or secondary fracture healing, a complex process in which endochondral ossification is an essential part and bone is regenerated by tissue differentiation. This process is sensitive to the mechanical environment, and several authors have proposed mechano-regulation algorithms to describe it using strain, pore pressure and/or interstitial fluid velocity as biofeedback variables. The aim of this study was to compare various mechano-regulation algorithms' abilities to describe normal fracture healing in one computational model. Additionally, we hypothesized that tissue differentiation during normal fracture healing could be equally well regulated by the individual mechanical stimuli, e.g. deviatoric strain, pore pressure or fluid velocity. A biphasic finite element model of an ovine tibia with a 3mm fracture gap and callus was used to simulate the course of tissue differentiation during normal fracture healing. The load applied was regulated in a biofeedback loop, where the load magnitude was determined by the interfragmentary movement in the fracture gap. All the previously published mechano-regulation algorithms studied, simulated the course of normal fracture healing correctly. They predicted (1) intramembranous bone formation along the periosteum and callus tip, (2) endochondral ossification within the external callus and cortical gap, and (3) creeping substitution of bone towards the gap from the initial lateral osseous bridge. Some differences between the effects of the algorithms were seen, but they were not significant. None of the volumetric components, i.e. pore pressure or fluid velocity, alone were able to correctly predict spatial or temporal tissue distribution during fracture healing. However, simulation as a function of only deviatoric strain accurately predicted the course of normal fracture healing. This suggests that the deviatoric component may be the most significant mechanical parameter to guide tissue differentiation during indirect fracture healing.     

The spatio-temporal arrangement of different tissues during bone healing as a result of simple mechanobiological rules

During secondary bone healing, different tissue types are formed within the fracture callus depending on the local mechanical and biological environment. Our aim was to understand the temporal succession of these tissue patterns for a normal bone healing progression by means of a basic mechanobiological model. The experimental data stemmed from an extensive, previously published animal experiment on sheep with a 3?mm tibial osteotomy. Using recent experimental data, the development of the hard callus was modelled as a porous material with increasing stiffness and decreasing porosity. A basic phenomenological model was employed with a small number of simulation parameters, which allowed comprehensive parameter studies. The model distinguished between the formation of new bone via endochondral and intramembranous ossification. To evaluate the outcome of the computer simulations, the tissue images of the simulations were compared with experimentally derived tissue images for a normal healing progression in sheep. Parameter studies of the threshold values for the regulation of tissue formation were performed, and the source of the biological stimulation (comprising e.g. stem cells) was varied. It was found that the formation of the hard callus could be reproduced in silico for a wide range of threshold values. However, the bridging of the fracture gap by cartilage on the periosteal side was observed only (i) for a rather specific choice of the threshold values for tissue differentiation and (ii) when assuming a strong source of biological stimulation at the periosteum.     

The mechanical heterogeneity of the hard callus influences local tissue strains during bone healing: a finite element study based on sheep experiments

During secondary fracture healing, various tissue types including new bone are formed. The local mechanical strains play an important role in tissue proliferation and differentiation. To further our mechanobiological understanding of fracture healing, a precise assessment of local strains is mandatory. Until now, static analyses using Finite Elements (FE) have assumed homogenous material properties. With the recent quantification of both the spatial tissue patterns (Vetter et al., 2010) and the development of elastic modulus of newly formed bone during healing (Manjubala et al., 2009), it is now possible to incorporate this heterogeneity. Therefore, the aim of this study is to investigate the effect of this heterogeneity on the strain patterns at six successive healing stages. The input data of the present work stemmed from a comprehensive cross-sectional study of sheep with a tibial osteotomy (Epari et al., 2006). In our FE model, each element containing bone was described by a bulk elastic modulus, which depended on both the local area fraction and the local elastic modulus of the bone material. The obtained strains were compared with the results of hypothetical FE models assuming homogeneous material properties. The differences in the spatial distributions of the strains between the heterogeneous and homogeneous FE models were interpreted using a current mechanobiological theory (Isakson et al., 2006). This interpretation showed that considering the heterogeneity of the hard callus is most important at the intermediate stages of healing, when cartilage transforms to bone via endochondral ossification.     

Fracture healing as a post-natal developmental process: molecular,spatial, and temporal aspects of its regulation

Fracture healing is a specialized post-natal repair process that recapitulates aspects of embryological skeletal development. While many of the molecular mechanisms that control cellular differentiation and growth during embryogenesis recur during fracture healing, these processes take place in a post-natal environment that is unique and distinct from those which exist during embryogenesis. This Prospect Article will highlight a number of central biological processes that are believed to be crucial in the embryonic differentiation and growth of skeletal tissues and review the functional role of these processes during fracture healing. Specific aspects of fracture healing that will be considered in relation to embryological development are: (1) the anatomic structure of the fracture callus as it evolves during healing; (2) the origins of stem cells and morphogenetic signals that facilitate the repair process; (3) the role of the biomechanical environment in controlling cellular differentiation during repair; (4) the role of three key groups of soluble factors, pro-inflammatory cytokines, the TGF-beta superfamily, and angiogenic factors, during repair; and (5) the relationship of the genetic components that control bone mass and remodeling to the mechanisms that control skeletal tissue repair in response to fracture.     

Mineral and matrix contributions to rigidity in fracture healing

The purpose of this study was to investigate the relationships among selected properties of fracture callus: bending rigidity, tissue density, mineral density, matrix density and mineral-to-matrix ratio. The experimental model was an osteotomized canine radius in which the development of the fracture callus was modified by electrical stimulation with various levels of direct current. This resulted in a range of values for the selected properties of the callus, determined post mortem at 7 weeks after osteotomy. We found that the rigidity (R) of the bone-callus combination obeyed relationships of the form R = axb, where x is the tissue density, mineral density, matrix density or the mineral-to-matrix ratio of the repair tissue. These are analogous to power-law relationships found in studies of compact and cancellous bone. The results suggest that fracture callus at 7 weeks after osteotomy in canine radius behaves more like immature compact bone than cancellous bone in its mechanical and physicochemical properties. The present study demonstrates the feasibility of developing non-invasive in vivo densitometric methods to monitor fracture healing, since models may be developed that can predict mechanical properties from densitometric data. Further studies are needed to develop a refined model based on experimental data on the mechanical and physicochemical properties and microstructure of fracture callus at different stages of healing.     

A mechano-regulation model for tissue differentiation during fracture healing: analysis of gap size and loading

Bone has a capability to repair itself when it is fractured. Repair involves the generation of intermediate tissues, such as fibrous connective tissue, cartilage and woven bone, before final bone healing can occur. The intermediate tissues serve to stabilise the mechanical environment and provide a scaffold for differentiation of new tissues. The repair process is fundamentally affected by mechanical loading and by the geometric configuration of the fracture fragments. Biomechanical analyses of fracture healing have previously computed the stress distribution within the callus and identified the components of the stress tensor favouring or inhibiting differentiation of particular tissue phenotypes. In this paper, a biphasic poroelastic finite element model of a fracture callus is used to simulate the time-course of tissue differentiation during fracture healing. The simulation begins with granulation tissue (post-inflammation phase) and finishes with bone resorption. The biomechanical regulatory model assumes that tissue differentiation is controlled by a combination of shear strain and fluid flow acting within the tissue. High shear strain and fluid flows are assumed to deform the precursor cells stimulating formation of fibrous connective tissue, lower levels stimulate formation of cartilage, and lower again allows ossification. This mechano-regulatory scheme was tested by simulating healing in fractures with different gap sizes and loading magnitudes. The appearance and disappearance of the various tissues found in a callus was similar to histological observation. The effect of gap size and loading magnitude on the rate of reduction of the interfragmentary strain was sufficiently close to confirm the hypothesis that tissue differentiation phenomena could be governed by the proposed mechano-regulation model.     

Simulation of the nutrient supply in fracture healing

The healing process for bone fractures is sensitive to mechanical stability and blood supply at the fracture site. Most currently available mechanobiological algorithms of bone healing are based solely on mechanical stimuli, while the explicit analysis of revascularization and its influences on the healing process have not been thoroughly investigated in the literature. In this paper, revascularization was described by two separate processes: angiogenesis and nutrition supply. The mathematical models for angiogenesis and nutrition supply have been proposed and integrated into an existing fuzzy algorithm of fracture healing. The computational algorithm of fracture healing, consisting of stress analysis, analyses of angiogenesis and nutrient supply, and tissue differentiation, has been tested on and compared with animal experimental results published previously. The simulation results showed that, for a small and medium-sized fracture gap, the nutrient supply is sufficient for bone healing, for a large fracture gap, non-union may be induced either by deficient nutrient supply or inadequate mechanical conditions. The comparisons with experimental results demonstrated that the improved computational algorithm is able to simulate a broad spectrum of fracture healing cases and to predict and explain delayed unions and non-union induced by large gap sizes and different mechanical conditions. The new algorithm will allow the simulation of more realistic clinical fracture healing cases with various fracture gaps and geometries and may be helpful to optimise implants and methods for fracture fixation.     

Regulation of fracture repair by growth factors.

Fractured bones heal by a cascade of cellular events in which mesenchymal cells respond to unknown regulators by proliferating, differentiating, and synthesizing extracellular matrix. Current concepts suggest that growth factors may regulate different steps in this cascade (10). Recent studies suggest regulatory roles for PDGF, aFGF, bFGF, and TGF-beta in the initiation and the development of the fracture callus. Fracture healing begins immediately following injury, when growth factors, including TGF-beta 1 and PDGF, are released into the fracture hematoma by platelets and inflammatory cells. TGF-beta 1 and FGF are synthesized by osteoblasts and chondrocytes throughout the healing process. TGF-beta 1 and PDGF appear to have an influence on the initiation of fracture repair and the formation of cartilage and intramembranous bone in the initiation of callus formation. Acidic FGF is synthesized by chondrocytes, chondrocyte precursors, and macrophages. It appears to stimulate the proliferation of immature chondrocytes or precursors, and indirectly regulates chondrocyte maturation and the expression of the cartilage matrix. Presumably, growth factors in the callus at later times regulate additional steps in repair of the bone after fracture. These studies suggest that growth factors are central regulators of cellular proliferation, differentiation, and extracellular matrix synthesis during fracture repair. Abnormal growth factor expression has been implicated as causing impaired or abnormal healing in other tissues, suggesting that altered growth factor expression also may be responsible for abnormal or delayed fracture repair. As a complete understanding of fracture-healing regulation evolves, we expect new insights into the etiology of abnormal or delayed fracture healing, and possibly new therapies for these difficult clinical problems.     

Computational simulation of fracture healing: influence of interfragmentary movement on the callus growth

Bone fractures heal through a complex process involving several cellular events. This healing process can serve to study factors that control tissue growth and differentiation from mesenchymal stem cells. The mechanical environment at the fracture site is one of the factors influencing the healing process and controls size and differentiation patterns in the newly formed tissue. Mathematical models can be useful to unravel the complex relation between mechanical environment and tissue formation. In this study, we present a mathematical model that predicts tissue growth and differentiation patterns from local mechanical signals. Our aim was to investigate whether mechanical stimuli, through their influence on stem cell proliferation and chondrocyte hypertrophy, predict characteristic features of callus size and geometry. We found that the model predicted several geometric features of fracture calluses. For instance, callus size was predicted to increase with increasing movement. Also, increases in size were predicted to occur through increase in callus diameter but not callus length. These features agree with experimental observations. In addition, spatial and temporal tissue differentiation patterns were in qualitative agreement with well-known experimental results. We therefore conclude that local mechanical signals can probably explain the shape and size of fracture calluses.     

A numerical framework for mechano-regulated tendon healing—Simulation of early regeneration of the Achilles tendon

Mechano-regulation during tendon healing, i.e. the relationship between mechanical stimuli and cellular response, has received more attention recently. However, the basic mechanobiological mechanisms governing tendon healing after a rupture are still not well-understood. Literature has reported spatial and temporal variations in the healing of ruptured tendon tissue. In this study, we explored a computational modeling approach to describe tendon healing. In particular, a novel 3D mechano-regulatory framework was developed to investigate spatio-temporal evolution of collagen content and orientation, and temporal evolution of tendon stiffness during early tendon healing. Based on an extensive literature search, two possible relationships were proposed to connect levels of mechanical stimuli to collagen production. Since literature remains unclear on strain-dependent collagen production at high levels of strain, the two investigated production laws explored the presence or absence of collagen production upon non-physiologically high levels of strain (>15%). Implementation in a finite element framework, pointed to large spatial variations in strain magnitudes within the callus tissue, which resulted in predictions of distinct spatial distributions of collagen over time. The simulations showed that the magnitude of strain was highest in the tendon core along the central axis, and decreased towards the outer periphery. Consequently, decreased levels of collagen production for high levels of tensile strain were shown to accurately predict the experimentally observed delayed collagen production in the tendon core. In addition, our healing framework predicted evolution of collagen orientation towards alignment with the tendon axis and the overall predicted tendon stiffness agreed well with experimental data. In this study, we explored the capability of a numerical model to describe spatial and temporal variations in tendon healing and we identified that understanding mechano-regulated collagen production can play a key role in explaining heterogeneities observed during tendon healing.     

Midkine-Deficiency Delays Chondrogenesis during the Early Phase of Fracture Healing in Mice

The growth and differentiation factor midkine (Mdk) plays an important role in bone development and remodeling. Mdk-deficient mice display a high bone mass phenotype when aged 12 and 18 months. Furthermore, Mdk has been identified as a negative regulator of mechanically induced bone formation and it induces pro-chondrogenic, pro-angiogenic and pro-inflammatory effects. Together with the finding that Mdk is expressed in chondrocytes during fracture healing, we hypothesized that Mdk could play a complex role in endochondral ossification during the bone healing process. Femoral osteotomies stabilized using an external fixator were created in wildtype and Mdk-deficient mice. Fracture healing was evaluated 4, 10, 21 and 28 days after surgery using 3-point-bending, micro-computed tomography, histology and immunohistology. We demonstrated that Mdk-deficient mice displayed delayed chondrogenesis during the early phase of fracture healing as well as significantly decreased flexural rigidity and moment of inertia of the fracture callus 21 days after fracture. Mdk-deficiency diminished beta-catenin expression in chondrocytes and delayed presence of macrophages during early fracture healing. We also investigated the impact of Mdk knockdown using siRNA on ATDC5 chondroprogenitor cells in vitro. Knockdown of Mdk expression resulted in a decrease of beta-catenin and chondrogenic differentiation-related matrix proteins, suggesting that delayed chondrogenesis during fracture healing in Mdk-deficient mice may be due to a cell-autonomous mechanism involving reduced beta-catenin signaling. Our results demonstrated that Mdk plays a crucial role in the early inflammation phase and during the development of cartilaginous callus in the fracture healing process.     

Influence of fracture gap size on the pattern of long bone healing: a computational study

Following fractures, bones restore their original structural integrity through a complex process in which several cellular events are involved. Among other factors, this process is highly influenced by the mechanical environment of the fracture site. In this study, we present a mathematical model to simulate the effect of mechanical stimuli on most of the cellular processes that occur during fracture healing, namely proliferation, migration and differentiation. On the basis of these three processes, the model then simulates the evolution of geometry, distributions of cell types and elastic properties inside a healing fracture. The three processes were implemented in a Finite Element code as a combination of three coupled analysis stages: a biphasic, a diffusion and a thermoelastic step. We tested the mechano-biological regulatory model thus created by simulating the healing patterns of fractures with different gap sizes and different mechanical stimuli. The callus geometry, tissue differentiation patterns and fracture stiffness predicted by the model were similar to experimental observations for every analysed situation.     

In-vivo imaging of the fracture healing in medaka revealed two types of osteoclasts before and after the callus formation by osteoblasts

The fracture healing research, which has been performed in mammalian models not only for clinical application but also for bone metabolism, revealed that generally osteoblasts are induced to enter the fracture site before the induction of osteoclasts for bone remodeling. However, it remains unknown how and where osteoclasts and osteoblasts are induced, because it is difficult to observe osteoclasts and osteoblasts in a living animal. To answer these questions, we developed a new fracture healing model by using medaka. We fractured one side of lepidotrichia in a caudal fin ray without injuring the other soft tissues including blood vessels. Using the transgenic medaka in which osteoclasts and osteoblasts were visualized by GFP and DsRed, respectively, we found that two different types of functional osteoclasts were induced before and after osteoblast callus formation. The early-induced osteoclasts resorbed the bone fragments and the late-induced osteoclasts remodeled the callus. Both types of osteoclasts were induced near the surface on the blood vessels, while osteoblasts migrated from adjacent fin ray. Transmission electron microscopy revealed that no significant ruffled border and clear zone were observed in early-induced osteoclasts, whereas the late-induced osteoclasts had clear zones but did not have the typical ruffled border. In the remodeling of the callus, the expression of cox2 mRNA was up-regulated at the fracture site around vessels, and the inhibition of Cox2 impaired the induction of the late-induced osteoclasts, resulting in abnormal fracture healing. Finally, our developed medaka fracture healing model brings a new insight into the molecular mechanism for controlling cellular behaviors during the fracture healing.     

Determining the most important cellular characteristics for fracture healing using design of experiments methods

Computational models are employed as tools to investigate possible mechanoregulation pathways for tissue differentiation and bone healing. However, current models do not account for the uncertainty in input parameters, and often include assumptions about parameter values that are not yet established. The objective of this study was to determine the most important cellular characteristics of a mechanoregulatory model describing both cell phenotype-specific and mechanobiological processes that are active during bone healing using a statistical approach. The computational model included an adaptive two-dimensional finite element model of a fractured long bone. Three different outcome criteria were quantified: (1) ability to predict sequential healing events, (2) amount of bone formation at early, mid and late stages of healing and (3) the total time until complete healing. For the statistical analysis, first a resolution IV fractional factorial design (L₆₄) was used to identify the most significant factors. Thereafter, a three-level Taguchi orthogonal array (L₂₇) was employed to study the curvature (non-linearity) of the 10 identified most important parameters. The results show that the ability of the model to predict the sequences of normal fracture healing was predominantly influenced by the rate of matrix production of bone, followed by cartilage degradation (replacement). The amount of bone formation at early stages was solely dependent on matrix production of bone and the proliferation rate of osteoblasts. However, the amount of bone formation at mid and late phases had the rate of matrix production of cartilage as the most influential parameter. The time to complete healing was primarily dependent on the rate of cartilage degradation during endochondral ossification, followed by the rate of cartilage formation. The analyses of the curvature revealed a linear response for parameters related to bone, where higher rates of formation were more beneficial to healing. In contrast, parameters related to fibrous tissue and cartilage showed optimum levels. Some fibrous connective tissue- and cartilage formation was beneficial to bone healing, but too much of either tissue delayed bone formation. The identified significant parameters and processes are further confirmed by in vivo animal experiments in the literature. This study illustrates the potential of design of experiments methods for evaluating computational mechanobiological model parameters and suggests that further experiments should preferably focus at establishing values of parameters related to cartilage formation and degradation.     

Predicting the formation of different tissue types during Achilles tendon healing using mechanoregulated and oxygen-regulated frameworks

During Achilles tendon healing in rodents, besides the expected tendon tissue, also cartilage-, bone- and fat-like tissue features have been observed during the first twenty weeks of healing. Several studies have hypothesized that mechanical loading may play a key role in the formation of different tissue types during healing. We recently developed a computational mechanobiological framework to predict tendon tissue production, organization and mechanical properties during tendon healing. In the current study, we aimed to explore possible mechanobiological related mechanisms underlying formation of other tissue types than tendon tissue during tendon healing. To achieve this, we further developed our recent framework to predict formation of different tissue types, based on mechanobiological models established in other fields, which have earlier not been applied to study tendon healing. We explored a wide range of biophysical stimuli, i.e., principal strain, hydrostatic stress, pore pressure, octahedral shear strain, fluid flow, angiogenesis and oxygen concentration, that may promote the formation of different tissue types. The numerical framework predicted spatiotemporal formation of tendon-, cartilage-, bone- and to a lesser degree fat-like tissue throughout the first twenty weeks of healing, similar to recent experimental reports. Specific features of experimental data were captured by different biophysical stimuli. Our modeling approach showed that mechanobiology may play a role in governing the formation of different tissue types that have been experimentally observed during tendon healing. This study provides a numerical tool that can contribute to a better understanding of tendon mechanobiology during healing. Developing these tools can ultimately lead to development of better rehabilitation regimens that stimulate tendon healing and prevent unwanted formation of cartilage-, fat- and bone-like tissues.

     

Mechanobiological modeling of endochondral ossification: an experimental and computational analysis

Long bone formation starts early during embryonic development through a process known as endochondral ossification. This is a highly regulated mechanism that involves several mechanical and biochemical factors. Because long bone development is an extremely complex process, it is unclear how biochemical regulation is affected when dynamic loads are applied, and also how the combination of mechanical and biochemical factors affect the shape acquired by the bone during early development. In this study, we develop a mechanobiological model combining: (1) a reaction–diffusion system to describe the biochemical process and (2) a poroelastic model to determine the stresses and fluid flow due to loading. We simulate endochondral ossification and the change in long bone shapes during embryonic stages. The mathematical model is based on a multiscale framework, which consisted in computing the evolution of the negative feedback loop between Ihh/PTHrP and the diffusion of VEGF molecule (on the order of days) and dynamic loading (on the order of seconds). We compare our morphological predictions with the femurs of embryonic mice. The results obtained from the model demonstrate that pattern formation of Ihh, PTHrP and VEGF predict the development of the main structures within long bones such as the primary ossification center, the bone collar, the growth fronts and the cartilaginous epiphysis. Additionally, our results suggest high load pressures and frequencies alter biochemical diffusion and cartilage formation. Our model incorporates the biochemical and mechanical stimuli and their interaction that influence endochondral ossification during embryonic growth. The mechanobiochemical framework allows us to probe the effects of molecular events and mechanical loading on development of bone.