Royal Society Publishing

In silico biology of bone modelling and remodelling: adaptation

Friederike A. Gerhard, Duncan J. Webster, G. Harry van Lenthe, Ralph Müller


Modelling and remodelling are the processes by which bone adapts its shape and internal structure to external influences. However, the cellular mechanisms triggering osteoclastic resorption and osteoblastic formation are still unknown. In order to investigate current biological theories, in silico models can be applied. In the past, most of these models were based on the continuum assumption, but some questions related to bone adaptation can be addressed better by models incorporating the trabecular microstructure. In this paper, existing simulation models are reviewed and one of the microstructural models is extended to test the hypothesis that bone adaptation can be simulated without particular knowledge of the local strain distribution in the bone. Validation using an experimental murine loading model showed that this is possible. Furthermore, the experimental model revealed that bone formation cannot be attributed only to an increase in trabecular thickness but also to structural reorganization including the growth of new trabeculae. How these new trabeculae arise is still an unresolved issue and might be better addressed by incorporating other levels of hierarchy, especially the cellular level. The cellular level sheds light on the activity and interplay between the different cell types, leading to the effective change in the whole bone. For this reason, hierarchical multi-scale simulations might help in the future to better understand the biomathematical laws behind bone adaptation.

1. Introduction

Bone is a material able to adapt to its external mechanical environment. This insight was gained more than a century ago by a German anatomist called Wolff (1892) and has been since referred to as ‘Wolff's law’. The underlying mechanisms and triggering factors of bone reorganization, however, are not well understood. In the 1970s, new discoveries (White et al. 1972; Whedon et al. 1976a,b; Rambaut & Johnston 1979) were made when it was realized that astronauts returned from space flights with significant bone loss. Likewise, it was observed that people exposed to prolonged bed rest also suffered from considerable bone atrophy (Donaldson et al. 1970; Pavosky 1977; Stupakov et al. 1981). These observations led to the assumption that bone needs everyday physical activity and the gravitational forces on Earth to sustain a structure capable of coping with daily physical demands. This form of bone loss seems to be driven by mechanics and is typically classified as disuse osteoporosis, as opposed to the so-called post-menopausal osteoporosis (hormone driven). These days, this disease is considered one of the major health problems facing middle-aged women and elderly people of either sex (Riggs & Melton 1995). In 2000, the total direct costs in Europe were estimated at 32 billion euros, which are expected to increase to 77 billion euros in 2050 based on the expected demographic shift in Europe (Kanis & Johnell 2005).

To combat post-menopausal osteoporosis, a deeper understanding of the microstructural adaptation mechanisms (called bone modelling and remodelling) is sought. Bone modelling is necessary to prepare the bone for increased or changing loads and is defined as the process of bone formation and resorption on separate surfaces, i.e. uncoupled formation and resorption (Frost 1990a). Bone remodelling undertakes the replacement of aged and damaged bone tissue with new (and undamaged) bone to prevent the development of fatigue fractures. Here, bone resorption and formation need to take place at the same site, i.e. coupled formation and resorption (Frost 1990b). These mechanisms occur both in the cortical shell and in the trabecular compartment. However, remodelling is shown to be more sensitive in the trabecular compartment. For this reason, this paper concentrates on trabecular bone adaptation. Figure 1a,b highlights the difference between cortical and trabecular bone as well as introducing the cells implicated in bone adaptation, respectively. The cellular mechanisms will be described in more detail in §2.

Figure 1

(a) Bone can be divided into cortical shell and trabecular compartment. Bone modelling and remodelling take place in both, but more in the trabecular compartment. Thus, (b) the cellular processes of modelling (M) and remodelling (R1–R4) are shown.

Since the adaptation of trabecular bone due to changing loading conditions or post-menopausal osteoporosis is a process that takes weeks up to years, data acquisition often represents a critical factor in research studies. Computer simulations have therefore gained more and more importance during the last 20 years. A computer simulation tries to describe the complex interdependent biological processes involved using a set of mathematical equations and attempts to mirror the actual biological process in silico. The quality of such models is therefore determined by the relevance and accuracy of the underlying assumptions upon which the governing mathematical equations have been formulated. While many of the earlier research questions aimed at a deeper understanding of the macroscopic mechanical properties, the scientific community has started to focus more and more on how these mechanical properties are determined by the underlying cellular processes, a field referred to as mechanobiology. This evolution has improved our understanding of the cellular mechanisms and molecular pathways governing bone adaptation. However, we must be careful not to neglect the global, macroscopic impact of these biological processes when we are interested in determining bone competence. Thus, in order to gain a deeper understanding of the whole process of bone modelling and remodelling, there is a need for hierarchical simulation models investigating different levels of resolution and incorporating these levels in a multi-scale approach. This paper provides an overview of the most important underlying biological principles and reviews the state of the art of hierarchical computational approaches that are used to describe the adaptation of bone due to mechanical loading and disease. Although a number of in silico models have been proposed for bone regeneration (see Geris et al. (2009) in this issue), this paper will only focus on adaptation models. To illustrate the power of in silico models, we will extend an established microstructural bone loss algorithm to describe bone adaptation due to increased mechanical loading. Its purpose is to investigate the feasibility of modelling bone adaptation without specific knowledge of the local mechanical milieu within the bone. A comparison of these results with findings from actual in vivo loading data shows that overall changes in bone mass can be predicted very accurately without knowledge of the underlying mechanics. Nevertheless, differences in local architecture can be noted between experiment and prediction, demonstrating the importance of incorporating local regulation of mechanical adaptation.

2. Biology

The actual change in the trabecular microstructure during mechanical adaptation is realized at the cellular level where the executive cells involved in the bone remodelling process are osteoblasts and osteoclasts (figure 1b). Osteoblasts are responsible for bone formation while osteoclasts resorb bone. The biochemical mechanisms triggering the activity and recruitment of these cell types, however, are still not identified (Cowin 2007). Today, it is widely assumed that osteoblasts and osteoclasts are controlled by a sensing system of osteocytes (Klein-Nulend et al. 1995, 2005; Mullender & Huiskes 1997; Burger & Klein-Nulend 1999; Cowin 2007). Osteocytes are mechanosensitive cells (Klein-Nulend et al. 1995) lying within the so-called osteocyte lacunae deep within the mineralized matrix and are highly interconnected by cell processes, lying in small canaliculi (figure 1b). The canalicular network is extended throughout the whole bone matrix and thus possesses the characteristics to serve as a communication system (Mullender & Huiskes 1997). For this reason, research has turned more and more towards the role of osteocytes and how they communicate with osteoblasts and osteoclasts to investigate processes related to mechanical bone adaptation (Mullender et al. 2004; Cowin 2007). The transduction from an external mechanical force into a signal sensed by the osteocytes may take place by strain (Huiskes et al. 1987, 2000; Han et al. 2004), micro-damage (Prendergast & Taylor 1994) or fluid flow (Knothe Tate et al. 1998; Burger & Klein-Nulend 1999; Nicolella et al. 2006). For osteocytes to activate the executive cells, it is assumed that a cascade of biochemical processes incorporating an enzyme called nitric oxide synthase occurs (Turner et al. 1996; Armour et al. 2001; Samuels et al. 2001).

The interactions between osteocytes, osteoblasts and osteoclasts are not known in detail. What is known is that during the bone remodelling process, osteoclasts first resorb a certain amount of bone and osteoblasts follow to refill the resulting resorption cavity (figure 1b). This unit of osteoclasts and osteoblasts is often referred to as the basic multicellular unit (BMU). It has been proposed that the activity of the BMU is governed by feedback from mechanical load transfer (Huiskes et al. 2000). According to this theory, osteoblasts are recruited in response to an osteocytic signal associated with high strain energy density (SED) values, while osteoclasts are either recruited by a low SED signal (disuse) or they act in a spatially random fashion (in response to the random occurrence of microcracks resulting from daily loading; Vashishth et al. 2000). Other theories, however, assume high strains to produce microcracks and thus attract osteoclasts to remove the damaged bone (Prendergast & Taylor 1994; McNamara et al. 2006). This controversy points to the fact that the actual cell activation mechanism has not yet been identified.

In the presence of post-menopausal osteoporosis, bone cells are impaired in their response to mechanical stimulation (Mullender et al. 2005). This in turn leads to bone loss and increased fracture risk. The pathological factors characterizing osteoporosis have been a topic of interest of many researchers to date, but are not the basis of this review.

3. In silico modelling

(a) Organ-level approaches (continuum models)

To gain better insight into the physiological interplay between the mechanical environment and the biological response of bone at a macroscopic level, early theoretical models aimed to provide a more detailed mathematical description of the law established by Wolff (1892). In these approaches, bone was considered as a continuum material and was described by its apparent density only. Cortical and trabecular components were accounted for by assigning two different ranges of apparent density (Fyhrie & Carter 1986), thus avoiding the need to describe their heterogeneous microstructures. Density-based simulations of trabecular bone adaptation have been addressed in many continuum approaches (Frost 1964; Cowin 1984; Fyhrie & Carter 1986; Reeve 1986; Huiskes et al. 1987; Brown et al. 1990; Martin 1991; Thomsen et al. 1996; Adachi et al. 1998; Langton et al. 1998; Hernandez et al. 2000) and have been reviewed before (Cowin 1993; Jacobs 2000; Hart 2001). To name only some of the approaches, several groups (Carter & Hayes 1976; Cowin 1984; Fyhrie & Carter 1986; Huiskes et al. 1987; Beaupre et al. 1990; Weinans et al. 1992; Van Rietbergen et al. 1993) proposed strain as the driving force for remodelling. A theoretical model using micro-damage as a driving force was published by Prendergast & Taylor (1994). Thomsen et al. (1994) focused on the simulation of osteoporotic bone loss on a continuum level by choosing a stochastic simulation approach.

The quantification of apparent bone density was without doubt a major advancement towards a better understanding of bone characteristics. However, with a greater understanding of the structural behaviour of bone, it became clear that the trabecular microstructure could not be ignored. With the arrival of new microstructural imaging techniques (Feldkamp et al. 1989), new methods of investigation became accessible and made it possible to visualize the organization of the trabecular network. This enabled researchers to better investigate the effects of post-menopausal osteoporosis on the trabecular architecture or even the realignment of individual trabeculae as a consequence of changes in loading. These enhanced imaging capabilities also allowed new ways to validate existing theories, including tissue-level properties.

(b) Tissue-level approaches

One of the first publications taking the intricate trabecular microstructure into account was published by Weinans et al. (1992). They presented a two-dimensional model based on finite-element (FE) analysis, which focused on the relationship between applied loads and the structural changes in single trabeculae. Trabeculae were modelled as a set of pixel elements and their morphological changes were simulated by removal or addition of the elements on the trabecular surface according to the local stresses and strains predicted by FE analysis. The theory behind the implemented algorithm was based on earlier work by Huiskes et al. (1987). However, the resulting two-dimensional artificial structures suffered from a so-called chessboard effect. This effect was expressed by pixel-sized holes within the two-dimensional representation of the bone structure. It was shown later that this behaviour was an effect of the mathematical implementation, and a solution to this problem was presented by Mullender et al. (1994), who allowed the number of osteocytes to be independent of the number of structural elements. Here, the authors assumed a certain sensing radius of the osteocytic signal, which was modelled as a decaying exponential function. From a biological point of view, this meant that several osteocytes could contribute to the stimulus an osteoblast receives. Evaluation of the model took place by visual inspection and comparison against the Weinans model.

Several studies based on the Mullender model followed (Smith et al. 1997; Zidi & Ramtani 1999). Mullender & Huiskes (1997) demonstrated in a follow-up study that the mechanosensory system indeed is more likely to consist of osteocytes than of surface lining cells (figure 1b).

Adachi et al. (1997) presented a competing two-dimensional FE model for trabecular surface remodelling, which assumed that an affinity for uniform stress conditions was the driving force behind bone remodelling. This algorithm was validated by applying it to a whole two-dimensional proximal femur using an artificially created trabecular microstructure.

Although these models appear to provide reasonable results when assessed qualitatively, one major limitation of all of the above models is that they lack quantitative validation using experimental datasets. A two-dimensional model validated against real experimental data was published for the first time by Tabor & Rokita (2002). They extended the stochastic continuum approach for osteoporosis that had been presented by Thomsen et al. (1994) and modified the original algorithm, so that the state of bone was represented as the state of all pixels. Thus, it became possible to account for the structure of real trabeculae. Comparing the results of this approach with histomorphometric data resulted in a good agreement. A limitation of this model, however, was its two-dimensional design. Because trabecular bone is a three-dimensional structure, only three-dimensional simulation approaches are capable of predicting microstructural changes reliably. Thus, the evolution of microstructural approaches advanced more towards three-dimensional simulations. Two of the early three-dimensional approaches were presented by Tayyar et al. (1999) and by Guo & Kim (2002). Both approaches focused on the simulation of bone loss over time by using artificial grids as input structures for successive bone removal. One of the first microstructural models using three-dimensional structures generated via micro-computed computer tomography (μCT) was presented by van der Linden et al. (2001). In this study, the authors investigated whether bone loss occurs predominantly by a disconnection of trabeculae, by an accrual of loose fragments or by a formation deficit of osteoblasts. The simulated three-dimensional bone structures were applied to a quantitative bone morphometric analysis. The resulting bone morphometric parameters were then successfully compared with experimental results. From this study, the authors found that the osteoblastic formation deficit was responsible for the majority of total bone loss, while bone loss owing to a disconnection of trabeculae and loose fragments seemed to play a minor role. This result disagrees with a competing study (Tayyar et al. 1999) on artificial structures, indicating that bone loss occurs predominantly due to a perforation and subsequent removal of single trabeculae. In a further study, van der Linden et al. (2003) used the developed model to predict the effects of an antiresorptive drug taken after menopause and found that late treatment can result in almost the same bone mass as early treatment after menopause; however, early treatment is shown to better conserve bone strength and stiffness.

Müller (2005) focused on the same issue and introduced an algorithm for simulated bone atrophy (SIBA). In his model, bone loss was, for the first time, accounted for in a realistic time frame, spanning several decades of remodelling. The process of osteoporotic bone loss was realized by means of a global iterative Gaussian filtration (osteoclastic resorption) and a subsequent thresholding procedure (osteoblastic refilling), employing several parameters to describe cellular activities in the different phases of hormonal development in women, i.e. pre-, peri- and post-menopausal phases. To validate this approach, the simulated structures were analysed and the resulting morphometric indices (bone volume density (BV/TV), specific bone surface (BS/BV), trabecular thickness (Tb.Th), trabecular spacing, trabecular number (Tb.N), degree of anisotropy and structure model index) from one simulated age group were compared with experimental results gained from the corresponding age group of seven human biopsies. Müller found no significant changes between the parameters from the simulation and experimental groups. Figure 2a shows the simulated course of BV/TV over time for seven individual specimens, resembling the actual course in ageing humans. Additionally, the simulated bone structures provided a realistic appearance (figure 2b). The conclusion of this study was that applying the thinning procedure without incorporation of strain knowledge may reflect a case where the osteocytes no longer sense strain- or flow-related impulses.

Figure 2

(a) BV/TV in the course of age-related bone loss as determined by SIBA, starting at the individual ages of each of the seven donors. The onset of menopause was assumed to be at 50 years. (b) Series of simulated images of an iliac crest bone biopsy of a 47-year-old woman. The simulation illustrates the changes in bone architecture over time, where with increasing age the plate-like structure is partly resolved and the rods forming the network are more and more thinned until eventual perforation.

Most of the three-dimensional approaches simulating bone loss were developed to provide functionality in the case of post-menopausal osteoporosis or age-related changes; however, they do not simulate the remodelling process of healthy, normal bone in response to increased (or decreased) loading conditions. To simulate the behaviour of normal bone in the presence of changing loads, these approaches may be extended. There also exist approaches that were originally developed to simulate normal bone behaviour in response to external loads. In the following paragraphs, two of the most promising three-dimensional approaches for bone remodelling will be described along with an extension of the model presented by Müller (2005).

Adachi et al. (2001) published a three-dimensional microstructural extension of their two-dimensional model (Adachi et al. 1997) for trabecular surface adaptation to new loading directions, which is based on the local uniform stress hypothesis. This considers the non-uniformity in the local stress distribution on the trabecular surfaces as the driving force of bone adaptation, directing structural changes towards a structure in which the stresses are distributed more homogeneously. In the first step, an FE analysis is performed in order to determine the local stress σc at any particular location xc. The non-uniformity of the stress distribution is then expressed by a functionEmbedded Imagewhere σd denotes the stress in the neighbourhood of σc and can be written asEmbedded Imagewhere S is the trabecular surface and w(l) is the weight function, considering the distance of a neighbouring point to the current one, e.g.Embedded Imageas illustrated in figure 3a. The total amount of surface movement is determined by substituting the value Γ(xc) into a function M as shown in figure 3b. The whole process is iterated until a remodelling equilibrium is attained. The authors validate their algorithm on cancellous bone (figure 3c) by comparing the morphometric indices of the artificially simulated structures with real ex vivo three-dimensional bone structures. This approach has not yet been used to investigate questions related to osteoporosis, although the potential might be there.

Figure 3

(a) The driving force of remodelling is defined as the relative difference between stress σc at the location xc and σd determined by the weighted sum of the neighbouring stresses σr. (b) Remodelling rate equation Embedded Image as a function of the driving force of remodelling Γ. (c) Changes in the three-dimensional architecture of a trabecular bone cube due to surface remodelling under compressive loading, starting from an initial voxel FE model based on a μCT digital image obtained from the canine distal femur. (Figure courtesy of Taiji Adachi.)

The second three-dimensional FE model for bone adaptation was proposed by Ruimerman et al. (2005a). This model had been originally designed for growth and mechanical adaptation in trabecular bone. The authors follow the theory presented by Huiskes et al. (2000). According to this theory, osteoblasts form new bones relative to the mechanical signals they receive from osteocytes where the osteocytes emit signals according to the SED sensed in their environment. Osteoclasts are recruited by osteocyte apoptosis or micro-cracks, which occur in a spatially random manner. The coupling factor between osteoblast formation and osteoclast resorption is of a mechanical origin. Embedding these biological assumptions into mathematical equations is as follows. The total amount of bone formation at a particular location in the bone surface is given byEmbedded Imagewhere the amount of bone resorbed by osteoclasts is a functionEmbedded Imageand the amount of bone formed by osteoblasts isEmbedded Imagewhere τ is the formation rate and ktr is the formation threshold. The osteocyte stimulus function P(x, t) is described byEmbedded Imagewhere f(x, xi) is the decay function representing the decrease in signal strength according to distance; μi is a mechanosensitivity factor; and R(xi, t) is the SED rate calculated by FE analysis.

Implementation of these equations resulted in a computer simulation model, which was validated using a 3×3×3 mm3 cube of artificially generated trabecular bone with respect to the following aspects. In a first study, Ruimerman et al. (2005a) successfully reproduced trabecular-like structures from a three-dimensional structured grid, which is illustrated in figure 4. Resulting morphometric indices resulted in reasonable dimensions; they also showed trabecular realignment after a change in the loading direction. In a second simulation series, alternative loading conditions were successfully applied to check whether the trabecular structure would adapt to a new equilibrium. Drastic reduction in loads resulted in disuse osteoporosis-like patterns with reduced bone mass and Tb.Th. In a third series, post-menopausal osteoporosis was simulated by increasing the osteoclast activation frequency. This entailed an increase in the osteoclast resorption frequency due to the mathematical formulation of the theory, and an increase in osteoblast recruitment due to the underlying mechanics. Still, the simulation led to a new equilibrium with reduced Tb.Th and connectivity density and an increase in average SED. In a second study, Ruimerman et al. (2005b) compared alternatives with SED as the osteocytic stimulus, i.e. maximal principal strain and volumetric strain. The authors found by comparison with results from the literature that SED is the most likely signal sensed by osteocytes. Furthermore, they found that oestrogen deficiency after menopause can only be successfully simulated by assuming that oestrogen inhibits osteoclast but not both osteoclast and osteoblast activities under normal conditions. The novelty of this theory was that it explained the effects of mechanical forces on trabecular-bone growth, maintenance and adaptation by relating mechanical stimuli in the bone matrix to assumed expressions in the cells involved in bone metabolism (Ruimerman et al. 2005a). The authors base their model on a number of assumptions that need still to be validated; however, the fact that they produce realistic results, speaks in favour of the model. One limitation might be that the volume of bone under investigation was relatively small, and thus it remains to be tested whether this model represents whole bone behaviour as well.

Figure 4

(a) Starting from a porous initial configuration representing bone in the post-mineralized foetal stage. During the simulation, (b–d) the structure developed in approximately 8 years into (e) a mature one. (f) From this point on the structure is maintained, while no more large architectural changes occur. Reprinted from Ruimerman et al. (2005a) with permission from Elsevier.

(c) Extending and validating a model for bone adaptation

As in silico models provide an opportunity to test different hypotheses, this section gives an example of how simulations can be used to further our understanding of the biological behaviour of bone in the presence of external loading.

In silico approaches should, in general, be able to reflect as many physiological and pathological situations as possible. Here, SIBA as presented by Müller (2005) was adapted in order to test the hypothesis that bone adaptation could be modelled by knowing the applied load only, without resolving the microstructural stresses and strains. This hypothesis was investigated by comparing the results of the simulation with real cross-sectional animal data gained from an experimental loading model, which was established to investigate the behaviour of bone in response to superphysiological loads. More details with respect to this model can be found elsewhere (Webster et al. 2008). Briefly, the fifth caudal vertebrae (C5) of 4 groups of 10 C57BL/6 (B6), 12-week-old, female mice were loaded in vivo with amplitudes of 0 (control), 2, 4 and 8 N, respectively, using a mechanical loading device (3000 cycles, 10 Hz, three times a week for four weeks). This treatment resulted in a substantial increase in BV/TV and other structural parameters for the different loading groups. The solid line in figure 5a illustrates the increase in BV/TV of the experiment.

Figure 5

(a) Experimental (CT, solid line) and simulated (SIBA, dashed line) increase in BV/TV due to cyclic loading of murine trabecular C5 vertebrae for different loads applied and simulated, respectively. (b) Regional analysis by dividing the trabecular compartment into 15 overlapping subregions. (c) One subregion as computed from the simulations of mechanical loading with (ii) 4 N and (iii) 8 N loads. The original (i) 0 N sample was used as the starting structure.

With respect to the simulations, the Müller algorithm for bone loss (Müller 2005) was used in such a way that binary input structures were blurred by applying the Gaussian filter and, in contrast to the original approach, a low subsequent threshold was chosen to simulate bone formation instead of bone loss. In order to account for directionality of the applied mechanical forces (in the axial direction of the vertebrae), the three-dimensional Gaussian filter kernel was compressed in the loading direction. The choice of the threshold and the size and extent of the Gaussian filter were obtained by optimally matching the simulated change in global mean BV/TV to the experimental change in global mean BV/TV for all loading groups, as is illustrated in figure 5a. To compare the simulated datasets with in vivo loading datasets and to investigate how bone formation varied within the bone, the trabecular compartments of the caudal mouse vertebrae were divided into 15 subregions (figure 5b). Figure 5c shows simulated changes in the fifth subregion for loads of 4 and 8 N when using the 0 N group as input datasets. Quantitative results of the regional variations between the in vivo and the in silico datasets are presented in figure 6. As a first result, the experimental and simulated curves of BV (figure 6a) and BV/TV (figure 6b) through the regions seem to match very well. Thus, the agreement between the in vivo and in silico datasets strengthened the hypothesis that knowledge about the global applied load may be sufficient to simulate bone adaptation due to mechanical loading.

Figure 6

Regional analysis between experimental (CT, solid line) and simulated (SIBA, dashed line) structures of (a) BV, (b) BV/TV, (c) Tb.Th and (d) Tb.N.

However, the simulated structures overestimated and underestimated Tb.Th (figure 6c) and Tb.N (figure 6d), respectively. An explanation of this second result may be given through the nature of the applied algorithm, which was designed for a uniform trabecular thickening and not for a structural reorganization including the creation of new trabeculae. Thus, the experimental increase in bone mass coming from structural reorganization was simply compensated by an increase in Tb.Th.

The fact that differences in local architecture could be noted between experiment and prediction manifests a current limitation of the model and demonstrates the importance of incorporating local regulation mechanisms of mechanical bone adaptation. Nevertheless, the simulations revealed that bone adaptation in response to mechanical loading is more complicated than uniform trabecular thickening. Instead, structural reorganization including the formation of new trabeculae seems to take place. Little is known, however, about how new trabeculae are formed under physiological conditions. To tackle this question, it would be necessary to perform longitudinal studies where in vivo bone adaptation can be tracked in individual mice over time.

(d) Multi-cell-level approaches

If we follow a hierarchical multi-scale approach, the next level of simulation would include cell-level and single-cell approaches. This section and §3e therefore give a short overview of current state of the art of in silico simulation models on cellular and single-cell levels.

At the tissue level, the localized activities of osteoclasts, osteoblasts and osteocytes are often implicitly merged to a scalar value, denoting bone formation or removal at a particular location. However, for a deeper understanding of how these cell types interact with each other to eventually renew or adapt the bone structure at the surface of a particular trabecula, it is important to depict these processes at a higher scale of resolution.

One of the first studies on a microscopic, but still multi-cellular, level was presented by Smit & Burger (2000). They built an FE model of a resorption cavity in a trabecula by means of a sphere residing in a simplified solid cylinder. Loading this FE model along the trabecular orientation resulted in elevated strain levels at the bottom of the resorption cavity. At the cavity border, however, low strains appeared in the direction of loading, while high strains occurred perpendicular to loading. From a biological point of view and under the assumption that osteoclasts remove bone at areas of low strain, this would mean that lateral osteoclastic excavation is inhibited, while remodelling in the longitudinal direction will proceed further. Although the authors assumed an idealized geometric structure and simplified homogeneous material properties, these findings support the hypothesis that BMU coupling is regulated by strains.

McNamara et al. (2006) extended this study by examining the strain distribution around a resorption lacuna by means of two-dimensional FE models of real μCT-scanned trabeculae. Their results supported the findings from Smit & Burger (2000), and further disclosed that the magnitudes of the peak strains in trabeculae are sufficiently high to initiate micro-damage and propagate to structures local to the resorption site. In another study, McNamara & Prendergast (2005) investigated whether an increased depth of resorption cavity, as occurs in osteoporosis, leads to a strain concentration capable of creating more micro-damage, such that osteoclasts ‘chase’ newly formed damage, leading to perforation. In this study, the mechanosensors are assumed to be the surface cells, i.e. osteoblasts and bone-lining cells. The results favour the concept of a critical cavity depth; larger cavities inevitably would lead to the loss of whole trabeculae. However, one limitation of this model might be that the sensing cells are located only on the trabecular surface and therefore will not be able to sense micro-damage occurring in areas deeply embedded in the bone matrix. In a subsequent study (McNamara & Prendergast 2007), it was hypothesized that bone remodelling may be regulated by signals due to both strain and micro-damage. The authors concluded from their results that an algorithm of strain-adaptive remodelling below a critical damage threshold and damage-adaptive remodelling above the critical amount presents the best computational algorithm for simulation of bone remodelling and repair. With respect to the sensing cells, it was found that with osteocyte mechanosensation, damaged surface bone tissue was not completely removed nor was refilling of resorption cavities predicted. However, in a separate study that was mentioned in their paper, only the osteocyte sensor model, but not the bone lining cell model, was able to communicate a signal to the bone surface, which originated from damaged tissue residing in the core of trabeculae.

Mulvihill et al. (2007) enhanced the above-proposed two-dimensional algorithms using a three-dimensional FE model of a cylindrical trabecular strut and found that, without restriction of osteoblast attachment at the resorption cavity, unusual bone formation behaviour was observed. At the border of the cavity, strain magnitudes were sufficient to form bone but not high enough to create micro-damage. These findings are in conflict with the suggestions made by Smit & Burger (2000).

Another in silico approach investigating the remodelling process at a multi-cellular level for both trabecular and cortical bone was presented by van Oers et al. (2008). They propose a unifying theory for osteonal and hemi-osteonal remodelling in a two-dimensional computer simulation model where osteoclasts are modelled explicitly using the Potts model for cell simulation (Glazier & Graner 1993). Here, the authors follow the theory proposed by Huiskes et al. (2000) and Ruimerman et al. (2005a) with enhanced attention to the osteoclast resorption function. It is postulated that osteoclast activity is inhibited by osteocytes instead of presuming random resorption along the bone surface. Osteoblasts only form bone at sites where the osteocyte signal exceeds a certain threshold for a particular period of time as occurring at the borders of a trabecular resorption cavity. In their study, it is found that stress concentrations around resorption cavities prevent osteoclasts from perforating trabeculae, and osteoblasts are recruited to refill resorption cavities.

(e) Single-cell-level approaches

In an effort to gain a better understanding of how mechanoregulation in bone remodelling works, it could be of special interest to investigate more questions on the level of a single cell. In particular, osteocytes are widely assumed to play the role of mechanosensors in the bone matrix, based on the strain magnitude they perceive and thus they could play a pivotal role in bone adaptation. Since the shape of an osteocyte lacuna influences the mechanical strain perceived by the osteocyte, a study was performed by McCreadie et al. (2004) to quantify the size and shape of osteocyte lacunae in trabecular bone near a typical fracture site. In their study, osteocytes are for the first time quantified in three dimensions, which was made possible by confocal microscopy images. In addition to the quantification of size and shape of the osteocytes, they proposed an ellipsoid as an appropriate geometric structure for an osteocyte lacuna. However, this study still neglected the interconnecting canaliculi. Bonivtch et al. (2007) used the proposed geometric specification to perform a three-dimensional FE analysis on this structure, in order to investigate the microstructural response to the imposed macroscopic apparent strains. They developed a parametric model consisting of the osteocyte lacuna, a region of perilacunar tissue, canaliculi and the surrounding bone tissue. It was found that the maximum perilacunar tissue strain increased in a model without canaliculi by a factor of 1.5 and by a factor of 3.1 in a model incorporating the presence of canaliculi. The maximum strain always occurred at sites where the canaliculi entered the lacuna. Furthermore, it was found that the predicted maximum perilacunar strain increased with a decrease in perilacunar tissue modulus and decreased with an increase in the tissue modulus, which entailed an increase in canalicular diameter. The authors concluded from this result that osteocytes may be able to change their local environment, i.e. the perilacunar tissue modulus. Whether the canalicular diameter indeed changes with an increase in tissue modulus remains to be verified experimentally.

A recent study (Anderson & Knothe Tate 2008) has examined quantitatively, based on earlier findings (Tami et al. 2002; Steck & Knothe Tate 2005), the effect of geometric and dimensional idealization on prediction of the mechanical signals imparted by fluid drag to cell surfaces. The study revealed that the physiologic geometry per se has a profound influence on the nano-scale flow regimes in bone, and that this influence imparts locally to cells through load-induced fluid flow.

McGarry & Prendergast (2004) went one level further and tried to simulate an adherent general eukaryotic cell by an FE modelling approach. Applying forces in the piconewton range to the model membrane nodes suggested a key role for the cytoskeleton in determining cellular stiffness. According to the authors, such a model could be useful when computing cellular structural behaviour in response to various in vitro mechanical stimuli, such as fluid flow or substrate strain to simulate mechanobiological processes. In any case, the simulation of mechanobiological and nanostructural coherences of cells and their properties, although not yet accounting for real bone adaptation but more focusing on micro-mechanics, has great potential to improve our understanding of the cellular interrelationships at a nanostructural level.

4. Conclusion and outlook

This paper gave an overview of current simulation models for macro-, meso- and nano-scales, and presented an extension of one of the existing simulation algorithms for trabecular bone remodelling. During the last few decades, more and more insights have been gained into the processes involved in bone adaptation, but also more and more questions have arisen regarding the underlying biochemical and mechanobiological pathways. With the availability of different imaging techniques to visualize and quantify different aspects of bone remodelling, we can create and validate simulation tools to reflect the behaviour of the human body at organ level, micro-level and single-cell level. The models presented here, which simulate the mechanical behaviour of bone at the organ level, represent bone by using its apparent material properties and apparent mechanical strength. This plays a significant role in the understanding of the stability of our bones within the body. By neglecting the trabecular microstructure, however, these models can no longer be used when questions related to the microstructural arrangement of trabeculae arise. Knowledge and mathematical theories established at organ level were therefore developed to account for the trabecular microstructure.

A current limitation of the proposed microstructural models is that they are only validated qualitatively via visual inspection and quantitatively by comparing the simulation with cross-sectional experimental studies, due to a lack of in vivo imaging techniques. Furthermore, some of these approaches are based on artificially generated input structures, and thus it would need to be proven that the governing principles used in these artificial grids can be transferred to actual bone micro-architectures. The application of these models to real bone and their validation after several simulation steps is the next step towards a better understanding of the underlying biological processes involved in bone adaptation.

The extension of the bone loss simulating algorithm also includes bone adaptive mechanisms that incorporate actual bone architectures from in vivo loading studies, and shows that overall changes in bone mass can be predicted with knowledge of the applied load only. However, differences in local architecture could be noted between experiment and prediction, thereby demonstrating the importance of incorporating local regulation mechanisms of mechanical adaptation, i.e. the formation of new trabeculae. Some of the presented simulation approaches account implicitly for reorganization by the growth of new trabeculae, i.e. by tunnelling of an existing structure. It is, however, not clear whether new trabeculae really arise owing to tunnelling of existing structures or rather by bridging between existing trabeculae. The only way to check the exact mechanism is in vivo imaging of the bone microstructure. These imaging techniques are now available and allow tracking of microstructural changes over time. Thus, they would facilitate the vigorous validation and comparison of competing theories. In that sense, in vivo imaging provides new forms of evaluation, i.e. the validation of a model with respect to its ability to predict trabecular changes over time both spatially and temporally.

Nevertheless, caution must still be exercised when attempting to validate, particularly when such models use parameters to represent certain mechanisms at the nano-scale, which cannot be accurately quantified by experimental means. Furthermore, the issue of validation becomes more complicated in more complex models employing a large number of parameters: as the number of parameters increases so does the number of experiments required to verify the model's applicability to scientific questions.

Despite an inherent degree of uncertainty, the nano- to macro-scale models discussed here have great value from a speculative point of view, as during the development process, they have uncovered many biological questions and can thus guide research accordingly. Further research will then, in turn, provide invaluable feedback for the development of more sophisticated and more reliable hierarchical models, forming an iterative process that will progressively further our understanding of bone adaptation. This process could eventually lead to the development of novel strategies for the management of bone diseases such as osteoporosis.


The authors acknowledge funding from the Swiss National Science Foundation (PP 104317/1).


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


View Abstract