In silico biology of bone modelling and remodelling: regeneration

L. Geris, J. Vander Sloten, H. Van Oosterwyck


Bone regeneration is the process whereby bone is able to (scarlessly) repair itself from trauma, such as fractures or implant placement. Despite extensive experimental research, many of the mechanisms involved still remain to be elucidated. Over the last decade, many mathematical models have been established to investigate the regeneration process in silico. The first models considered only the influence of the mechanical environment as a regulator of the healing process. These models were followed by the development of bioregulatory models where mechanics was neglected and regeneration was regulated only by biological stimuli such as growth factors. The most recent mathematical models couple the influences of both biological and mechanical stimuli. Examples are given to illustrate the added value of mathematical regeneration research, specifically in the in silico design of treatment strategies for non-unions. Drawbacks of the current continuum-type models, together with possible solutions in extending the models towards other time and length scales are discussed. Finally, the demands for dedicated and more quantitative experimental research are presented.

1. Introduction

Bone is a remarkable material as, under most circumstances, it is capable of truly regenerating itself, whereas soft tissue wound healing results in scar formation (Glowacki 1998). However, despite bone's natural healing capacity and the extensive amount of research that has been conducted in this area, delayed healing and non-unions often develop. For instance, 5–10 per cent of the over 6 million fractures annually occurring in the USA develop into delayed or non-unions (Einhorn 1995, 1998; Praemer et al. 1999), costing society large amounts of money. In 2000, in the European Union, the direct cost for the 3.8 million osteoporosis-related fractures was estimated at 32 billion euros (Reginster & Burlet 2006). This cost is sure to rise in the future, in light of the ageing population and the prediction that 40 per cent of all postmenopausal women will suffer one or more fractures during their remaining lifetime (Compston et al. 1998; Reginster & Burlet 2006). Therefore, prevention and effective treatment are highly desirable.

Potential causes of osteoporosis and the use of mathematical models in designing preventive treatment strategies are discussed elsewhere in this issue by Gerhard et al. (2009). The current paper focuses on the (in silico) biology of the bone regeneration process, taking place after bone traumata such as fractures or implant placement. It starts with a summary of the most important biological processes taking place during bone regeneration, and subsequently discusses the modelling efforts that have been undertaken in this research domain hitherto. Examples are given to illustrate the added value of mathematical models. Finally, drawbacks and possible solutions related to the current models are discussed, as well as the demands for dedicated and more quantitative experimental research.

2. Bone regeneration

(a) Biology of bone regeneration

Bone regeneration involves complex physiological processes, involving the participation of many cell types, regulated by a myriad of biochemical and mechanical factors. The principal mechanisms underlying the healing around implants are comparable to those of fracture healing (Schenk & Herrmann 1984; Gross 1988; Plenk & Zitter 1996; Overgaard 2000), although this view has triggered much debate since some authors have suggested that the presence of an implant induces a different healing mode (Schwartz et al. 1997; Steflik et al. 1998).

The regeneration response can be divided into primary and secondary healing (Overgaard 2000). Primary healing involves a direct healing without the appearance of inflammation or the formation of callus tissue. Primary healing seems to occur only under optimal conditions, i.e. mechanical stability and the absence of gaps (in fracture healing, the anatomical restoration is needed). When such conditions are present, the basic multicellular units will directly re-establish the Haversian canals in the bone. Secondary healing is the most common form of healing and takes place when the optimal conditions needed for primary repair are absent. The secondary healing of bone, whether it is in fracture healing or implant osseointegration, encompasses three overlapping phases: the inflammatory phase; the reparative phase; and the remodelling phase (Einhorn 1995, 1998; Overgaard 2000; Hadjiargyrou et al. 2002; Davies 2003; Gerstenfeld et al. 2003). Figure 1 gives a schematic overview of these different phases.

Figure 1

Schematic of different phases of fracture healing. (a) The inflammatory phase, (b) the soft callus stage of the reparative phase, (c) the hard callus stage of the reparative phase and (d) the remodelling phase (adapted from Bailón-Plaza & van der Meulen 2001).

Following bone trauma, the cortical bone, periosteum and the surrounding soft tissues are damaged, rupturing numerous blood vessels. This blood rapidly coagulates to form a clot enclosing the trauma site. In this clot, activated and degranulating platelets provide a source of growth factors, playing an important role in the wound-healing cascade. This is the start of the initial stage of bone regeneration, the inflammation phase. As a consequence of the vascular damage, the trauma site becomes hypoxic. Osteocytes at the trauma site become deprived of their nutrition and necrose, as do the damaged tissues in that area. This necrotic process triggers an immediate inflammatory response, bringing inflammatory cells, leucocytes and macrophages to the region. These are followed by the invasion of fibroblasts, mesenchymal stem cells and endothelial cells (Taguchi et al. 2005). Growth factors and cytokines, important regulators of the healing process, are produced by the cells present in the regeneration area as well as released into this area from the surrounding tissues (damaged bone ends, muscles, periosteum and marrow; Gerstenfeld et al. 2003; Malizos & Papatheodorou 2005).

The inflammation phase is followed by the reparative phase. In fracture healing, the reparative phase can be further divided into two subphases: the soft and the hard callus phases (Einhorn 1995, 1998; Hadjiargyrou et al. 2002; Gerstenfeld et al. 2003). The initially formed granulation tissue is gradually replaced by fibrous tissue, forming the soft callus. In this soft callus, mesenchymal stem cells will start differentiating, guided by cues from their microenvironment, such as the presence of biological factors (Einhorn 1998) and the perceived mechanical stimuli (Carter et al. 1998). A direct differentiation towards osteoblasts is observed near the cortex, away from the fracture site. These osteoblasts produce woven bone matrix, the hard callus, in a process called intramembranous ossification. Mesenchymal stem cells differentiate into chondrocytes (cartilage-forming cells) in the central fracture area, where the soft callus will gradually take on the appearance of cartilage, mechanically stabilizing the fracture zone. The chondrocytes mature towards hypertrophic chondrocytes, initiating the biochemical preparations in the cartilage matrix to undergo calcification. This concludes the soft callus phase of healing. In the next phase, the hard callus phase, blood vessels invade the calcified cartilage, bringing along osteoblasts. These osteoblasts, receiving enough oxygen and subjected to the proper mechanical stimuli, will produce a hard callus tissue consisting of mineralized woven bone matrix in a process called endochondral ossification. When the fracture ends are connected by a bony callus, clinical union is reached.

Already during the reparative phase, the final remodelling phase begins with osteoclastic resorption of unnecessary or poorly placed parts of the regenerated bone and the formation of lamellar bone. This remodelling takes place for a prolonged period of time, gradually reverting the blood supply to a normal state and restoring the bone at the regeneration site to its original shape and strength. More information on the biology and in silico modelling of the remodelling process can be found elsewhere in this issue (Gerhard et al. 2009).

(b) Challenges in fracture-healing research

Non-unions, hypertrophic or atrophic, and synovial pseudarthrosis generally are differentiated based on radiographic and histological appearances (Rodriguez-Merchan & Forriol 2004). Hypertrophic non-unions show an abnormal vascularity and callus formation, with a horseshoe or elephant-foot (abundant callus) configuration on radiographs. Atrophic non-unions typically show little callus formation around a fibrous tissue-filled fracture gap. Synovial pseudarthrosis is defined as a non-union in which a fluid-filled cavity with a synovial-like membrane is present at the fracture site. There is no universally accepted definition of non-union of a fracture. If a fracture fails to heal within the time usually required (i.e. known from clinical experience), it is called a delayed union. Generally, the failure of a fracture to heal in six to eight months (in humans) constitutes a non-union. The fracture gap of a non-union usually consists of fibrous tissue with varying amounts of cartilage.

Numerous adverse mechanical and biological factors influence the development of a non-union: excessive motion; a large interfragmentary gap; loss of blood supply; and severe periosteal and soft tissue trauma (McKibbin 1978; Whiteside 1978; Hulth 1989; Oni et al. 1989; Marsh et al. 1994; Cañadell & Forriol 1997; An et al. 1999; Park et al. 1999; Landry 2000). General factors such as old age, cachexia and malnutrition, steroids or anticoagulants, anti-inflammatory agents, burns and radiation may contribute to the establishment of non-unions, but are not the primary causes (Hulth 1989; Aaron 1996). Classical treatment strategies for non-unions are aimed at restoring optimal mechanical circumstances to induce healing. Excessive motion of the fracture fragments will be reduced by the use of plates, external or intramedullary fixators. However, the application of appropriate dynamic stimulation of the fracture has been shown to enhance the healing process (Goodship & Kenwright 1985; Goodship et al. 1998; Kenwright & Gardner 1998). The use of growth factors for enhancing fracture healing is a relatively novel treatment modality. Results of experimental and clinical studies have provided support for the clinical use of growth factors for treatment of acute and non-union fractures, as reviewed by Southwood et al. (2004). Despite the large amount of information existing in the literature, additional research is required to determine the exact mechanisms of non-union, and subsequently the optimal therapeutic strategies for each type of non-union.

3. Mathematical models of bone regeneration

(a) From mechanics …

The first theoretical models to describe the mechanoregulation of skeletal tissue differentiation were introduced by Pauwels (1960). He recognized that physical factors cause stress and deformation of the mesenchymal cells and that these stimuli could determine the cell differentiation pathway. He hypothesized that elongation stimulates the formation of fibrous connective tissue. Hydrostatic pressure on the other hand is a specific stimulus for the formation of cartilaginous tissue. A specific stimulus for the formation of bone, however, is not present in the theory of Pauwels. According to Pauwels (1960), bony tissue proceeds on the basis of a rigid framework of fibrous tissue, cartilage or bone. Some years later, the concept of interfragmentary strain (IFS) was developed by Perren (1979) and Perren & Cordey (1980). They proposed that a tissue, with a certain failure strain, cannot be formed in a region that experiences strains higher than this level. This means that a fracture gap can only be filled with a tissue capable of sustaining the IFS without failure. The IFS concept provides a theoretical basis for evaluating fracture treatment strategies, but is not applicable for bone regeneration in general, since it disregards the structural and mechanical heterogeneity of the fracture callus.

Building on the theories of Pauwels (1960), Perren (1979) and Perren & Cordey (1980), many research groups have formulated a theoretical model relating tissue differentiation to mechanical loading. Some of the most widely used mechanoregulatory models are represented in figure 2. Carter (1987) and Carter et al. (1998) specifically discussed the importance of cyclic tissue loading and proposed a mechanical stimulus that takes into account the local stress or strain history. The stress acting on the regenerating tissue is described in terms of hydrostatic stress (leading to a change in volume) and distortional strain (a measure for the change in shape). Direct bone formation is permitted in regions experiencing low hydrostatic stresses and low distortional strains. Carter, however, did not provide values for which stimuli will favour bone formation. Claes & Heigele (1999) combined the result of their finite-element analyses of the mechanical stimuli on ossifying surfaces during fracture healing with a histological analysis of a real callus geometry. From this, they derived a quantitative mechanoregulatory model relating magnitudes of hydrostatic pressure and principal strain to the bone formation process. Intramembranous bone formation occurred at the regenerating bone surface for hydrostatic pressures smaller than 0.15 MPa and strains smaller than 5 per cent, while endochondral bone formation occurred when the compressive hydrostatic pressure exceeded this threshold (and strains remain smaller than 15%).

Figure 2

Schematic of the mechanoregulatory models proposed by (a) Pauwels (1960), (b) Prendergast et al. (1997), (c) Carter et al. (1998) and (d) Claes & Heigele (1999). Dark grey, endochondral ossification; light grey, intramembranous ossification; white, connective tissue or fibrocartilage.

The above-mentioned theories all considered tissues as solid elastic materials. Prendergast et al. (1997) proposed a mechanoregulatory model for tissue differentiation, based on the poroelastic (biphasic) behaviour of the tissues. Maximal distortional strain and relative fluid velocity constitute the stimulus that controls the differentiation process. High values of both solid strain and fluid velocity favour fibrous tissue formation, while intermediate values lead to cartilaginous tissue. Bone can be formed only if the values are sufficiently low. Huiskes et al. (1997) quantified the upper and lower limits of the mechanical stimulus for the different tissue phenotypes, which yielded a mechanoregulatory diagram for tissue differentiation. The model was extended to include mesenchymal cell migration (by means of a diffusion equation), as a first step to incorporate the underlying cellular processes (Lacroix & Prendergast 2002; Lacroix et al. 2002). Similar to the model developed by Prendergast et al. (1997), Kuiper et al. (2000) presented a model where shear strain and fluid shear stress regulate bone fracture healing. Callus growth was incorporated in the model controlled by the proliferation of granulation tissue and realistic healing patterns were obtained. Another approach was followed by Ament & Hofer (2000) who implemented a set of fuzzy logic rules guiding the tissue differentiation using strain energy density as mechanical stimulus. Table 1 summarizes the most important mathematical models of tissue differentiation, mentioning their most important characteristics. The numerical studies that apply these models without adding or changing them have been omitted from the table.

View this table:
Table 1

Summary of mathematical models of tissue regeneration, indicating their major constituents. (SED, strain energy density; GF, growth factor; MSC, mesenchymal stem cell; FB, fibroblast; CC, chondrocyte; OB, osteoblast; CGGF, chondrogenic growth factor; OGGF, osteogenic growth factor; VGF, vascular growth factor; EC, endothelial cell.)

Over the past years, the mechanoregulatory models described above were applied to various experimental situations to assess their validity. Carter's model (Carter et al. 1988, 1998) was used in several computational studies to investigate oblique fracture healing (Blenman et al. 1989), pseudarthrosis formation (Loboa et al. 2001), asymmetric clinical fractures (Gardner et al. 2006) and distraction osteogenesis (Morgan et al. 2006). However, none of the studies predicted tissue differentiation adaptively over time. Lacroix & Prendergast (2002) successfully applied the model of Prendergast et al. (1997), which was calibrated originally with an experiment using a dynamic implant device in the femoral condyles of Labrador dogs (Søballe 1993), to normal fracture healing. Andreykiv et al. (2005) applied the same model to investigate glenoid component loosening in shoulder implants. This allowed them to draw conclusions on the importance of primary stability and glenoid component stiffness. Geris et al. (2003, 2004, 2008b) applied different mechanoregulatory models to study implant osseointegration in an in vivo repeated sampling bone chamber. Certain aspects of the numerical results were corroborated by the experiments, such as the amount of bone formed under specific implant loading conditions. However, other aspects of the numerical results, such as the numerically predicted presence of cartilage in the bone chamber, were contradicted by the experimental observations where no cartilage could be detected. Isaksson et al. (2006a,b) compared the predictive power of different mechanoregulatory models, applying them to an identical experimental set-up of fracture healing. All models were able to simulate the course of normal fracture healing correctly. However, under torsional loading (i.e. a different mechanical condition from the one applied in the experiments the models were derived from), only the model of Prendergast et al. (1997) was able to successfully predict fracture healing. Applying the latter mechanoregulatory model to the experimental set-up of tibial distraction osteotomy (Isaksson et al. 2007), spatial and temporal tissue distributions as seen experimentally, as well as variations in predicted tissue distributions due to alterations in distraction rate and frequency were obtained. A variation of the mechanoregulatory model of Prendergast et al. (1997) was used by Idelsohn et al. (2006) to investigate mandibular distraction osteogenesis. Kelly & Prendergast (2005) applied the same mechanoregulatory model to simulate tissue regeneration in osteochondral defects and were able to capture patterns of cellular differentiation observed experimentally in this healing process.

The use of mechanoregulatory models has increased our quantitative understanding of how bone regeneration is influenced by mechanical loading. They can be used to derive quantitative guidelines for in vivo loading regimes that can be tolerated by differentiating tissues (an idea that was already present in the work of Perren and co-workers) or that can even enhance regeneration (Isaksson et al. 2006b). A clear shortcoming of these mechanoregulatory models is their (almost) complete lack of biological mechanisms and how those mechanisms may be mediated by mechanical signals. Instead, the algorithms describe a relationship between mechanical signals and tissue differentiation (changes in tissue composition and corresponding mechanical properties) in a rather phenomenological way. This obviously limits their applicability to studies where the effect of mechanical loading parameters is investigated, instead of biological parameters (as, for example, in the case of cellular therapies or growth factor treatments). In addition, these algorithms do not allow us to formulate and investigate hypotheses on how mechanics may modulate biological processes.

(b) To biology …

On the other side of the modelling spectrum are the bioregulatory models that describe bone regeneration as a process directed only by biological factors such as growth factors and precursor cells.

A simple mathematical model for calvarial wound healing has been proposed by Adam (1999) and Arnold & Adam (1999). Their model described the spatial distribution of a generic growth stimulator concentration in the wound site. They investigated the influence of ranges of parameters related to this generic growth factor and the influence of the wound geometry on the occurrence of critical size defects (Arnold & Adam 1999; Adam 2002; Bellomo 2003).

Bailón-Plaza & van der Meulen (2001) have developed a more extensive mathematical model describing the reparative phase of the secondary fracture healing process. The mathematical model studies the effects of growth factors during both intramembranous and endochondral ossification. The model constitutes seven variables, namely the concentration of mesenchymal stem cells, chondrocytes and osteoblasts, the concentration of an osteogenic and a chondrogenic growth factor and the density of a combined fibrous/cartilaginous extracellular matrix and the bone matrix. The cell densities change due to cell migration, proliferation, differentiation, endochondral replacement and cell removal, while the changes in the matrix density are the result of synthesis and resorption. Most of these processes are mediated by the osteogenic and chondrogenic growth factors. The concentration of growth factors, which are produced by the cells, changes due to diffusion synthesis and decay.

Based on the implant osseointegration experiments (Søballe 1993) that were used by Prendergast et al. (1997) to develop their model, Ambard & Swider (2006) proposed a bioregulatory model of bone implant healing, where the migration of osteoblasts towards the implant surface is mediated both by the presence of an osteogenic growth factor and the local matrix density.

Geris et al. (2008a) started from the model of Bailón-Plaza & van der Meulen (2001) and introduced key aspects of healing such as angiogenesis and directed cell migration. A schematic overview of this model is presented in figure 3. The model describes the spatial and temporal evolution of the concentration of mesenchymal stem cells, fibroblasts, chondrocytes, osteoblasts and endothelial cells. These cells are responsible for the production of the corresponding matrix types (fibrous tissues, cartilage, bone and vascular tissue). Similar to the model of Bailón-Plaza & van der Meulen (2001), many processes are mediated by the presence of a generic representative of the chondrogenic, osteogenic or vascular growth factor family. Additionally, in the model of Geris et al. (2008a), angiogenesis plays a major role in controlling the regeneration process. In the description of endochondral ossification, mature chondrocytes start producing the vascular growth factor that attracts endothelial cells, which in turn allow osteoblasts to exist in that area, leading to cartilage degradation and concurrent bone formation.

Figure 3

Schematic of the bioregulatory model developed by Geris et al. (2008b).

The authors also evaluated the models' ability to predict certain pathological cases of fracture healing and took a first step towards the in silico design of therapeutic strategies (Geris et al. 2008c). An example of an atrophic non-union based on Reed et al. (2002) and several treatment strategies are shown in figure 4. Non-union was achieved experimentally by removing the periosteum and bone marrow at both sides of the fracture gap, which was numerically simulated by removing the boundary condition for mesenchymal stem cells in these areas. Figure 4a shows the domain and the boundary conditions for the simulation of atrophic non-unions. A good correspondence was observed between experimental observations (Reed et al. 2002) and numerical simulations regarding the distribution of the various tissues in the callus (fibrous tissue, cartilage, bone and vasculature) at different time points (results not shown). Six weeks after fracture induction, the callus is filled with a well-vascularized fibrous tissue while cartilage and bone are absent. The proposed treatment strategies shown in figure 4 are the administration of mesenchymal stem cells at fracture induction (figure 4b), the increase of the proliferative capacity of mesenchymal stem cells at fracture induction (corresponding to the experiments by Makino et al. (2005), where transforming growth factor-β3 was applied at fracture induction; figure 4c) and the administration of mesenchymal stem cells three weeks after fracture induction (figure 4d). Both cases (b) and (c) result in union after five to six weeks, although in case (b) there is still a substantial amount of fibrous tissue present. Clinically, this indicates that by the administration of cell and/or growth factors very early in the healing process, the normal regeneration process could be restored in fractures with substantial soft tissue damage. The model further predicts that cellular therapy is much less effective when cells are administered three weeks after fracture induction (leading to only partial regeneration after 45 days), instead of this being done immediately (figure 4d). These results show the added value of bioregulatory models for the exploration of cellular or growth factor therapies. At the same time, it is clear that there is still a long way to go to validate (and further refine) these models before they can be applied (pre-)clinically.

Figure 4

Simulation of potential treatment strategies of atrophic non-unions. (a(i)) Geometrical domain derived from a rodent fracture callus at post-fracture week 3 (Harrison et al. 2003) and (ii) indication of the boundary conditions used for the numerical simulations. (b–d(i)) Tissue fractions in the fracture callus over time (solid curve, bone; dotted curve, cartilage; dashed curve, fibrous tissue) and (ii) the distribution of the bone matrix in the callus at various time points for different treatment strategies (black, bone; grey, soft tissue). (b) Administration of mesenchymal stem cells at fracture induction, (c) increase in the proliferative capacity of mesenchymal stem cells at fracture induction and (d) administration of mesenchymal stem cells three weeks after fracture induction.

(c) And back again

Some of the most recent models in the literature combine both mechanical and biological factors. Bailón-Plaza & van der Meulen (2003) adapted their bioregulatory model and made a number of the model parameters, such as proliferation, endochondral differentiation and matrix production, dependent on the local mechanical environment. After calibration to well-defined in vivo experiments, the beneficial and adverse effects of moderate and excessive (delayed) mechanical stimulation were investigated. A fuzzy logic model was proposed by Simon et al. (2003) and Shefelbine et al. (2005) to simulate fracture healing, based on local mechanical factors, as described by Claes & Heigele (1999), and the local vascularity. To model the progress of angiogenesis in the fracture callus, a number of basic mechanoregulatory fuzzy rules were applied.

Gomez-Benito et al. (2005) and Garcia-Aznar et al. (2007) presented a model describing bone regeneration in terms of the evolution of the densities of different cell types in a growing callus. Several cellular events (migration, proliferation, differentiation and death) as well as tissue synthesis, damage, calcification and remodelling over time are represented in the model. They further aimed to analyse the main components that constitute the matrix of the different tissues, such as collagen, proteoglycans, mineral and water, and derived mechanical properties of the tissues from that composition. They applied the second invariant of the deviatoric strain tensor as the stimulus steering the differentiation process, thereby neglecting the influence of biological (growth) factors. They used their model to investigate the influence of gap size and interfragmentary movement on the callus size.

Andreykiv et al. (2007) started from the mechanoregulatory model proposed by Prendergast et al. (1997) and included diffusion, proliferation and differentiation of cell populations (mesenchymal stem cells, fibroblasts, chondrocytes and osteoblasts). Model parameters were selected based on the comparison with an in vivo experiment. The model is able to show asymmetric distribution of tissues in the callus under high bending moments. A similar model was proposed by Isaksson et al. (2008a). It was applied among others to investigate the influence of periosteal stripping and impaired endochondral ossification potential on the fracture healing process under normal loading conditions.

A final model was presented by Geris et al. (2008d). Starting from the bioregulatory model (Geris et al. 2008b), a number of biological parameters were made dependent on the local mechanical environment. Figure 5 shows the implementation scheme for such a model. As an example, the proliferation rates of osteoblasts and endothelial cells were made dependent on the local fluid flow according to the graph in figure 5 (bottom right). Figure 6 shows the results for the experiments and simulations of various implant loading amplitudes in a bone chamber set-up. This repeated sampling bone chamber is implanted in rabbit tibiae and has a cylindrical implant positioned in the centre. Cells and tissues can enter the chamber through perforations in the chamber wall. The dual chamber design enables us to investigate the effect of various implant loading conditions on the peri-implant tissue dynamics, while excluding site- and animal-specific influences (Duyck et al. 2004; figure 6a). Numerical results were obtained using both the model of Prendergast et al. (1997) and the mechanobioregulatory model of Geris et al. (2008d). Figure 6b shows the distribution of various tissues at various time points. The main difference between the simulation results from the two models is the prediction of cartilage formation in the bone chamber by the model of Prendergast et al. (1997), whereas the mechanobioregulatory model does not predict any cartilage formation at all as observed experimentally. Looking at quantitative results (figure 6c,d), no good correspondence between either of the models and the experiments is obtained. For the 30–90 μm loading cases, the model of Prendergast indicated higher bone content (expressed as a bone area fraction) for higher implant displacements, similar to the experiments. Contrarily, for this same range of implant displacement, the mechanobioregulatory model gives a better prediction of the amount of bone to implant contact.

Figure 5

Overview of the implementation of the coupled mechanobioregulatory model (Geris et al. 2008d).

Figure 6

Bone ingrowth in a repeated sampling in vivo bone chamber, a comparison of experimental and simulation results. (a) Schematic of the bone chamber and its components (top) (Duyck et al. 2004). Indication of the domain used in the numerical simulations, as well as the boundary conditions (bottom). (b) Simulation results using the model of (i) Prendergast et al. (1997) and (ii) Geris et al. (2008d) showing the predicted tissues in the bone chamber after nine weeks for various implant displacement magnitudes. (c) Comparison of the experimental and simulation results (dark grey, Duyck et al. 2006; grey, Geris et al. 2008b; light grey, Prendergast et al. 1997) in terms of (i) the amount of bone formed inside the chamber and (ii) the amount of bone to implant contact after six weeks.

Another application of the mechanobioregulatory model (Geris et al. 2008b) is shown in figure 7 for fracture healing. The absence of bone formation was predicted under severe overloading conditions (figure 7b; overload equal to five times physiological load), mainly due to the adverse effect of high loads on the angiogenic process, as specified in the model (figure 5; excessive fluid flow reduces endothelial cell proliferation). Figure 7b,c shows the results of two different treatment strategies, both involving the stabilization of the fracture at post-fracture day 21. Whereas clinically this would be achieved by the application of a fixator, the axial force acting on the fracture was reduced to half of the physiological load in the simulations. Upon removal of the overloading situation, the normal fracture healing process was restored (figure 7b). Additional administration of osteogenic growth factors (within the entire callus area) in combination with the stabilization did not alter the predicted outcome (figure 7c). A somewhat different tissue distribution pattern was observed, yet leading to a similar end result, namely an ossified callus bridging the fracture gap several weeks after stabilization.

Figure 7

Recapitulation of normal healing after reduction of callus overload. (a(i)) Geometrical domain derived from a rodent fracture callus at post-fracture week 3 (Harrison et al. 2003) and (ii) indication of the boundary conditions used for the numerical simulations. (b–d(i)) Schematic of the simulated conditions and (ii) resulting predicted distribution of the bone matrix in the callus at various time points for different treatment strategies. (b) Overload, (c) overload followed by load reduction at post-fracture week 3 and (d) overload followed by load reduction and simultaneous osteogenic growth factor administration at post-fracture week 3.

Similar to mechanoregulatory models, this mechanobioregulatory model allows us to numerically test various treatment strategies in which the mechanical conditions, such as the fixator stiffness and the time of application of the fixator, are modulated. In addition, it is able to evaluate the effectiveness of ‘combined’ treatments, and, for example, investigate whether there could be an additive or synergistic effect of adding growth factors. It also permits us to formulate hypotheses with respect to the mechanisms by which mechanics may modulate biological processes (such as angiogenesis). Obviously, this enhanced model potential has a certain price in terms of (significantly increased) model complexity.

4. Prospects

The models currently used for the simulation of bone regeneration are all of the continuum type, describing only a single time and length scale. In the case of bioregulatory models, many of the mathematical expressions used to describe phenomena such as migration and proliferation are derived from actual processes taking place at a lower level. One example is the diffusion term describing random cell migration, which is derived from a probabilistic micro-model of a random walk process (Murray 2002). An explicit incorporation of other levels might lead to additional insights into the dynamics of the fracture healing process. In the simulation of tumour development and treatment, multi-scale models are frequently used to explicitly model the cell cycle and angiogenesis (Anderson & Chaplain 1998; Alarcon et al. 2004) with the specific aim of the in silico design and testing of new therapies (Byrne et al. 2006; McDougall et al. 2006; Anderson & Quaranta 2008).

For the coupled mechanobioregulatory models discussed in this review, biological processes such as cell proliferation and differentiation are directed by tissue-level stresses and strains. These tissues are currently still assumed to be linear elastic or poroelastic at best (table 1), which gives, at best, a rough approximation of the actual tissue behaviour. The values used for the stiffness of the tissues featuring in the mathematical models for bone regeneration (typically, granulation tissue, fibrous tissue, cartilage and bone) have been shown to influence the simulation predictions. The stiffness of the soft tissues has been shown to influence not so much the order in which the typical fracture healing stages are predicted but rather the time that is needed to reach complete healing (Isaksson et al. 2008b). With the focus of mathematical modelling in bone regeneration shifting towards the in silico design of treatment strategies, this is a limitation that needs to be addressed. Experimentally measuring the material properties of the various tissues formed during the regeneration process is a challenging task due to the high degree of tissue heterogeneity (Moukoko et al. 2007). Another issue is the currently still unanswered question as to how mechanical signals are transduced from the tissue level down to the cell or even the intracellular level. The calculated tissue-level mechanical stimuli are used to adapt the biological parameters, whereas in reality a substantial difference is measured between tissue-level and cellular-level stimuli. You et al. (2001) proposed a model to explain (and calculate) the observed strain amplification phenomenon in the osteocytes' processes (see also Gerhard et al. 2009). These osteocytes are considered the main mechanosensors in bone (Mullender & Huiskes 1997; Burger & Klein-Nulend 1999). However, mechanosensing by cells in a regenerating tissue remains to be unravelled as yet.

Using multi-scale models, it is possible to realistically incorporate the events taking place at the different length scales (global and local) on the condition that the additions to the model are based on validated biophysical relationships (see review by Hunter & Borg 2003).

This obviously leads to the final issue of this review: the sparseness of the experimental data needed to verify model assumptions and predictions at all levels. For almost all of the models discussed in §3, experimental results were used which were obtained from studies designed with specific questions in mind, other than the corroboration of the mathematical model results. Leucht et al. (2007) present a study where an implant osseointegration set-up is developed. Experimental results, measured in terms of cellular and biological response, are compared with numerical calculations of the mechanical situation inside the regeneration zone, allowing for a direct coupling between the mechanics and the fate decisions of the cells in the regeneration zone.

In order to populate the parameter sets of the developed mathematical models, not just more data but also more quantitative data are needed. As the models contain variables related to tissues, cells and molecules, in principle, quantitative data are needed at these different levels. Either in vitro or in vivo experiments could be designed that enable quantification. In vitro set-ups offer the advantage of a highly controllable and quantifiable environment. However, care should be taken when translating the conclusions to the actual in vivo environment under investigation, as cells and tissues are isolated from their natural environment. Either cell or tissue (organ) culture models could be developed. Cell behaviour, such as stem cell proliferation, commitment, differentiation and organization, seems to be mediated by cell morphology (McBeath et al. 2004; Nelson et al. 2005), which in turn may differ between two- and three-dimensional cell culture models. As cells are likely to be anchored to a three-dimensional matrix during fracture healing, three-dimensional culture models, such as natural polymers (fibrin and collagen), seem more appropriate to study biological interactions relevant for fracture healing (Gabbay et al. 2006; Weinand et al. 2006; Kasper et al. 2007). In contrast to skeletal development, where organ culture models are frequently used to study cellular and molecular interactions (see Minina et al. (2001) as an example of a limb culture system), a very limited number of tissue or organ culture studies can be found related to fracture healing (Einhorn et al. 1995; Igarashi & Yamaguchi 1999).

The use of in vivo experimental models of fracture healing obviously has the advantage of resembling reality, but quantitative measurements will be much more challenging. As the mathematical models predict molecular, cell and tissue dynamics as a function of time and space, there is a need for temporal and spatial experimental data. The use of imaging techniques presents an elegant way to non-invasively monitor and quantify the dynamics of fracture healing. At the tissue level, in vivo microquantitative X-ray computed tomography has been used to quantify the bone tissue distribution during fracture healing, especially in small animals such as rats and mice (Kalpakcioglu et al. 2008). At the cellular and molecular levels, molecular imaging techniques can provide data that seem essential for the development (such as parameter identification) and validation of the mathematical models presented here. As reviewed by Mayer-Kuckuk & Boskey (2006), imaging techniques such as positron emission tomography (PET), single photon emission computed tomography (SPECT), magnetic resonance imaging (MRI), bioluminescence and fluorescence imaging have been used in orthopaedic research (among others fracture healing) to monitor and quantify real-time biology, such as gene expression, cell migration and cell death. Not only are these techniques quantitative, but some of them such as PET, SPECT and MRI are also tomographic, so that they can provide spatial information.

5. Conclusion

In this review three different classes of mathematical models of bone regeneration were presented: mechanoregulatory; bioregulatory; and mechanobioregulatory models. The appropriateness of a certain (class of) model(s) or the required model complexity strongly depends on the details of the research study, such as research goal and questions, application, selected outcome variables, etc. Therefore, a general conclusion on these aspects cannot be drawn.

In order to advance bone regeneration research in general, a close coupling of both experimentally and mathematically driven research is extremely important. This coupling is needed during both the model development phase (‘model fitting’, parameter identification) and the model validation phase (evaluate model predictions for experimental situations other than the ones used during model development). Only in this way can truly predictive models be obtained that can be used to in silico test research hypotheses and design treatment strategies.


L.G. is a postdoctoral research fellow of the Research Foundation Flanders (FWO Vlaanderen). The authors declare that they have no potential conflict of interest and they had and will have no financial relationships with companies whose products are relevant to the subject of this study.


  • One contribution of 15 to a Theme Issue ‘The virtual physiological human: tools and applications I’.


View Abstract