The influence of bone surface availability in bone remodelling—A mathematical model including coupled geometrical and biomechanical regulations of bone cells
Institute for Mechanics of Materials and Structures, Vienna University of Technology, A-1040 Vienna, Austria
Mechanical & Mechatronic Engineering, University of Sydney, Sydney NSW 2006, Australia
July 2, 2019
Bone is a biomaterial undergoing continuous renewal. The renewal process is known as bone remodelling and is operated by bone-resorbing cells (osteoclasts) and bone-forming cells (osteoblasts). An important function of bone remodelling is the repair of microcracks accumulating in the bone matrix due to mechanical loading. Cell-cell communication between cells of the osteoclastic lineage and cells of the osteoblastic lineage is thought to couple resorption and formation so as to preserve bone integrity and achieve homeostatic bone renewal. Both biochemical and biomechanical regulatory mechanisms have been identified in this coupling. Many bone pathologies are associated with an alteration of bone cell interactions and a consequent disruption of bone homeostasis. In osteoporosis, for example, this disruption leads to long-term bone loss and fragility, and can ultimately result in fractures.
Here we focus on an additional and poorly understood potential regulatory mechanism of bone cells, that involves the morphology of the microstructure of bone. Bone cells can only remove and replace bone at a bone surface. However, the microscopic availability of bone surface depends in turn on the ever-changing bone microstructure. The importance of this geometrical dependence is unknown and difficult to quantify experimentally. Therefore, we develop a sophisticated mathematical model of bone cell interactions that takes into account biochemical, biomechanical and geometrical regulations. We then investigate numerically the influence of bone surface availability in bone remodelling within a representative bone tissue sample. Biochemical regulations included in the model involve signalling molecules of the receptor–activator nuclear factor B pathway (rank–rankl–opg), macrophage colony-stimulating factor (mcsf), transforming growth factor (tgf) and parathyroid hormone (pth). For the biomechanical regulation of bone cells, a multiscale homogenisation scheme is used to determine the microscopic strains generated at the level of the extravascular matrix hosting the osteocytes by macroscopic loading. The interdependence between the bone cells’ activity, which modifies the bone microstructure, and changes in the microscopic bone surface availability, which in turn influences bone cell development and activity, is implemented using a remarkable experimental relationship between bone specific surface and bone porosity. Our model suggests that geometrical regulation of the activation of new remodelling events could have a significant effect on bone porosity and bone stiffness. On the other hand, geometrical regulation of late stages of osteoblast and osteoclast differentiation seems less significant. We conclude that the development of osteoporosis is probably accelerated by this geometrical regulation in cortical bone, but probably slowed down in trabecular bone.
Key words: mechanical feedback, geometrical feedback, specific surface, bone remodelling, bone stiffness, osteoporosis
Bone is a biomaterial that has a variety of physiological functions. In addition to load bearing and support for locomotion, bone protects internal organs and participates in calcium and phosphorous homeostasis. From an engineering perspective the structural function of bone is most importantly characterised by its stiffness and strength. Daily activities (such as walking and running) subject bone to periodical loads which, over extended periods of time (weeks, months and years), can lead to fatigue damage and the formation of microcracks. If these microcracks are not removed in due time, their evolution may result in a macroscopic structural failure, i.e., a fragility fracture. To prevent the occurrence of fatigue fractures, nature has equipped bone tissues with a cellular mechanism of self-repair , referred to by biologists as ‘bone remodelling’ [2, 3]. Bone remodelling is a coordinated process of bone resorption by cells called osteoclasts, and bone formation by cells called osteoblasts. Osteoclasts and osteoblasts usually operate together in self-contained groups processing the renewal of a localised portion of the bone tissue. These groups are called bone multicellular units (bmus) and constitute a single ‘remodelling event’. There are about such bmus in a normal adult skeleton [4, 3, 2]. Cell population and cell activity in a bmu are tightly controlled to establish local bone homeostasis (i.e., balanced bone resorption and bone formation). In bone pathologies, this cellular control is perturbed and homeostatic bone renewal is disrupted. In osteoporosis, bone is progressively lost, which results in reduced bone stiffness and strength.
Over the last decades, bone biologists have identified a large number of biochemical regulatory factors influencing bone remodelling. The formation of osteoclasts has been shown to rely crucially on macrophage colony-stimulating factor (mcsf) and on the receptor-activator nuclear factor B (rank) cell signalling pathway, which involves the receptor rank, the ligand rankl and osteoprotegerin (opg) [5, 6]. Rankl activates the rank receptor on precursor osteoclasts, which triggers their development and sustains their activity. The soluble molecule opg is a decoy receptor of rankl which can prevent rankl from binding to rank. Another important molecule mediating the communication between osteoblasts and osteoclasts is transforming growth factor (tgf). Tgf is stored in high concentrations in the bone matrix. It is released into the bone microenvironment, where it exerts its action on several bone cells, during bone matrix resorption by active osteoclasts . The existence of a mechanical regulation of bone remodelling has long been suspected. It is now well established that mechanical feedback is a key regulatory mechanism to maintain bone mass [7, 8, 9, 10, 11, 12]. The commonly accepted view is that osteocytes act as mechanosensors that transduce local mechanical signals into biochemical responses. These biochemical responses are thought to regulate the initiation of bone remodelling processes and to modulate the coupling between bone resorption and formation (see e.g.  and Refs cited therein).
The existence of biochemical and biomechanical regulations of bone cells is well-established and has been extensively studied. However, the notion that the morphology of the microstructure of bone may induce an additional regulation of bone cells of purely geometrical nature is not often mentioned in the recent literature. This may be due to the experimental difficulty of assessing the importance of a geometrical regulation. Biochemical and biomechanical regulations can experimentally be partially or fully repressed by selective gene knock-outs or monoclonal antibodies targeting key components in the bone cell signalling pathways. By contrast, one cannot simply “switch off” a geometrical regulation of bone remodelling when this self-repair process modifies the microstructure (and so the geometry) of the material.
Bone tissue is diverse and exhibits a broad variety of microstructures. However, two distinctive types of bone tissue are usually identified: cortical bone and trabecular bone  (see images in Figure 1). Cortical bone has typical porosities of – while trabecular bone has typical porosities of – [3, 2]. Mathematical models for the estimation of mechanical properties of bone tissue have shown that bone stiffness is predominantly determined by the porosity ,111The total porosity of bone is made of a vascular porosity, which contains marrow components, blood vessels, bone cells and their precurors, and the lacunae-canaliculi porosity, which contains osteocytes and their processes. The lacunae-canaliculi porosity is only a small fraction of the total porosity (see e.g. [13, Table 1]) and no remodelling occurs at these surfaces. Therefore, the lacunae-canaliculi porosity will not be consider in this work and we will refer to the vascular porosity simply as the bone porosity. Similarly, in the present context, we are not interested in the intercrystalline and intermolecular porosities, which we simply regard as part of the ‘solid bone matrix’. the interaction of the different material phases and pore shape, while other microstructural characteristics such as the exact pore distribution play a secondary role [14, 15, 16]. For biochemical processes, pore morphology can be expected to play a significant role. Indeed, pore morphology determines the so-called specific surface (i.e., the amount of bone surface available in a representative volume element), which is an essential geometrical factor for the bone cells. Bone cells require a bone surface to fulfill their functions, whether to initiate a bone remodelling process or to operate resorption and formation. Osteoclasts require attachment to a particular area of the bone surface before resorbing. Osteoblasts are observed to only secrete osteoid (a collagen-rich substance which later mineralises and becomes new bone matrix) on existing bone surfaces. Finally, mechanical signals sensed by osteocytes embedded in the bone matrix are passed on to bone cells in the vascular cavity through the bone surface. Effects similar to chemical exchange reactions between pore walls and solutes in fluid-saturated porous materials can be expected to occur in this context.
The issue of quantifiying the role of bone surface availability in bone remodelling was raised by bone biologists already some time ago [17, 13, 2]. In Ref. , Martin provides a first attempt to investigate theoretically the effect of a geometrical regulation of bone remodelling in osteoporosis (see Figure 1). Osteoporosis is associated with increased porosity in both cortical and trabecular bone . In Martin’s own words: “In [cortical] bone, increased porosity provides more surface area on which cells can work, thereby increasing the capacity for further porosity changes. In [trabecular] bone, increased porosity decreases the amount of surface available to the cells, thereby decreasing the capacity for further remodelling.”
While the proposed mechanism of geometrical feedback on bone remodelling seems plausible, it is difficult to test its validity experimentally and to determine its importance quantitatively. Some researchers have employed the concept of geometrical feedback for simulations of bone remodelling . However, to our knowledge, there is no systematic study in the literature of the effects of a possible geometrical regulation at several stages of the remodelling sequence. Also, the interplay between geometrical feedback and mechanical feedback in bone remodelling has not been investigated. A mechanical feedback has the potential to stabilise bone loss or gain [7, 22] and may either compete with or enhance the effect of the geometrical feedback seen in Figure 1 depending on the type of bone.
The aim of this paper is to address the above questions using a state-of-the-art computational model of bone remodelling. The following questions related to geometrical feedback are investigated: (i) At which stage of the bone remodelling sequence (i.e., activation, resorption, formation) does geometrical feedback have the strongest effect?; (ii) How do geometrical and mechanical feedbacks interact?; (iii) What is the impact of geometrical feedback in osteoporosis in terms of bone porosity and bone stiffness?
To address these questions we extend a previously developed mathematical model of bone remodelling [23, 24, 22]. This multiscale model takes into account both biochemical and biomechanical regulations of bone remodelling. Biochemical regulatory factors include the rank–rankl–opg pathways together with the action of tgf on bone cells. Biomechanical regulation of bone formation and bone resorption is mediated by the microscopic strain energy density (sed) of the bone matrix. This strain energy density is calculated from a micromechanical homogenisation scheme. The introduction of geometrical feedback due to microscopic bone surface availability is elucidated through a phenomenological relationship between the specific surface and the vascular porosity obtained from various types of bone .
2 Mathematical model of bone remodelling
Few mathematical models of bone remodelling include explicitly biochemical interactions of bone cells that couple bone resorption and bone formation. Lemaire et al.  have proposed a bone cell population model that incorporates some of the most important known bone biology. This cell population model was further developed by Pivonka et al. to investigate the effect of rankl expression on osteoblasts of varying maturity . These models have been shown to correctly capture important physiological behaviours of bone remodelling both in bone homeostasis and in bone pathologies [24, 26, 27, 28]. Recently, we have proposed an extension of the model of Ref.  to include a biomechanical regulatory mechanism to trigger bone cell responses, leading to a fully coupled model of biochemical and biomechanical regulations . This model uses a novel multiscale approach based on micromechanical homogenisation of bone stiffness. It allows to consistently calculate bone matrix strains at the osteocyte level. Osteocytes convert this micromechanical signal into biochemical signals to bone cells in the vascular space, effectively providing a biomechanical feedback regulation of bone remodelling.
In the following, we present an extension of the model of Ref.  that incorporates the influence of microscopic bone surface availability at various stages of the bone remodelling sequence.
2.1 Fundamental biochemical and biomechanical regulation of bone remodelling
The biochemical regulatory mechanisms considered in the model have been described in detail in Refs. [23, 24]. Three stages of osteoblast development are considered: uncommitted osteoblast precursors (obs), pre-osteoblasts (obs) and active osteoblasts (obs). Similarly, three stages of osteoclast development are considered: uncommitted osteoclast precursors (ocs), pre-osteoclasts (ocs) and active osteoclasts (ocs). Figure 2 schematically depicts these bone cell types along with their biochemical, biomechanical and geometrical regulations. These regulations are summarised sequentially in the following.
The generation of osteoblasts is assumed to be regulated by transforming growth factor (tgf) (a factor released from the bone matrix during osteoclastic bone resorption) while the generation of osteoclasts is assumed to be regulated by rankl and opg (molecules in the receptor-activator nuclear factor B (rank) system, that are expressed by osteoblasts). Tgf promotes the differentiation of uncommitted osteoblast progenitors (ob) into pre-osteoblasts (ob), but it inhibits the differentiation of pre-osteoblasts into active osteoblasts (ob). Furthermore, tgf promotes osteoclast apoptosis (programmed cell death). Rankl is a protein expressed on the surface of pre-osteoblasts. The binding of rankl to the receptor rank found on pre-osteoclasts promotes the differentiation of pre-osteoclasts into active osteoclasts. However, this binding may be prevented by the presence of opg, a decoy receptor of rankl produced in soluble form by active osteoblasts. Furthermore, the circulating parathyroid hormone pth induces rankl expression and downregulates opg expression by osteoblasts to produce a systemic increase in available rankl, resulting in an increase in the formation and activity of osteoclasts. The differentiation of uncommitted osteoclast progenitors (ocs) into pre-osteoclasts (ocs) requires signalling by both rankl and macrophage colony-stimulating factor (mcsf). All these signalling actions are accounted for using mass action kinetics for the chemical bindings between ligands and their receptors . The fraction of occupied receptors on a cell (i.e., bound to a ligand) is assumed to determine the strength of the signal received by the cell, and so is assumed to determine the strength of the cell’s response.
The biomechanical regulatory mechanisms considered in the model are described in detail in Ref. . Mechanical disuse is known to increase bone resorption by increasing the rankl/opg ratio, which increases the differentiation of pre-osteoclasts into active osteoclasts. In the model, this is implemented using a mechanically controlled rankl-production term, see Eq. (10). Mechanical overuse is believed to increase bone formation by stimulating wnt signalling to pre-osteoblasts, which increases their proliferation and ultimately leads to an increased population of active osteoblasts [22, 30]. In the model, this is implemented using a mechanically controlled proliferation of pre-osteoblasts, see Eq. (11).
2.2 Geometrical and morphological characteristics of bone: porosity and specific surface
To represent the microstructure of bone at the tissue level, the most important geometrical and morphological parameters are the vascular porosity () and the specific surface (). In cortical bone, vascular porosity corresponds to the so-called ‘Haversian canal’ system. In trabecular bone, vascular porosity corresponds to the marrow space around the trabecular struts . The vascular compartment contains all the bone cells considered in our model. Vascular porosity is defined as the volume fraction of vascular pores, i.e., the volume of vascular pores () per tissue volume ():222The tissue volume is assumed to be of the order of 1–3 \milli\metre. This volume corresponds to a representative volume element (RVE) large enough to comprise several remodelling events (bmus) , yet small enough to represent spatial heterogeneity of bone tissue (in particular, small enough to distinguish cortical bone and trabecular bone) .
The bone matrix volume fraction is defined in the same way as the volume of bone matrix () per tissue volume, i.e.:
We recall that all porosities at observation scales below the vascular porosity, such as the lacunar porosity, and the canaliculi connecting the lacunae, are not involved in the remodelling process, i.e. the intricate biochemical processes take place within the vascular pore space (see also footnote 1). From Eqs (1)–(2), it follows that . Typically, cortical bone exhibits a range of porosities – while trabecular bone exhibits a range of porosities –.
The specific surface (or surface density) of a porous material is defined as the interstitial surface area of the pores () per tissue volume, i.e.:
with dimensions . Generally speaking, the specific surface is an important quantity for a variety of phenomena in porous media. For example, it determines the adsorption capacity of industrial adsorbents and plays an important role in determining the effectiveness of catalysts and ion exchange columns and filters. It is also related to the fluid conductivity or permeability of porous media (see e.g. Ref. ). Experimentally, is commonly estimated by adsorption methods, quantitative stereology,333Stereology is the study of three-dimensional properties of objects observed in two-dimensional sections. fluid flow and micro-computed tomography. As mentioned above, in bone the specific surface is important as it determines the available working area for osteoblasts and osteoclats. The specific surface can also be expected to have an influence on the transmission of specific signalling by osteocytes in the bone matrix to osteoblasts and osteoclasts developing in the vascular pores.
The microstructure of a material determines both the porosity and the specific surface. Depending on the microstructure, different materials exhibit different values for these quantities. Bone tissue spans a wide range of porosities, each of which is characteristic of a particular micro-architecture, and so of a particular value of the specific surface. Based on a large number of experimental data, Martin has provided a remarkable phenomenological relationship between bone specific surface () and vascular porosity () [13, Eq. (68)]:
where the polynomial coefficients are estimated as , and (in ). In Figure 3, the relation (4) is plotted together with experimental data obtained from various types of human bone (femur, iliac crest, vertebra, rib) both in health and disease. The data and the curve show two important characteristics: (i) all specimens approximately follow the same curve independently of bone type and with no significant difference between diseased and healthy bone. This remarkable universality establishes the curve as an intrinsic property of bone; (ii) the specific surface exhibits a maximum at a porosity of about 0.37, intermediate between cortical and trabecular bone (bone at such porosity is denoted as ‘porous-compact’ bone).
In this work, we use Eq. (4) to include geometrical regulation in the model as follows. The evolution of the bone cell populations predicts the evolution of the vascular porosity due to resorption by osteoclasts and formation by osteoblasts. Eq. (4) then enables us to estimate changes in the specific surface associated to the changes in porosity. This change in microscopic bone surface availability is in turn assumed to influence the evolution of the cell populations.
2.3 Bone cell governing equations
In the model, osteoblasts and osteoclasts of different developmental stages are considered. The populations of osteoblasts and osteoclasts are therefore heterogeneous. The composition of these populations is determined by following individually the populations of each of the developmental stages of osteoblasts and osteoclasts mentioned above. The evolution of these cell (sub)populations is transcribed mathematically as so-called ‘rate equations’ [23, 22]. The populations of uncommitted progenitor osteoblasts (obs) and uncommitted progenitor osteoclasts (ocs) are assumed constant and so are not state variables. These uncommitted cells represent a pool of progenitor (or stem) cells that is assumed to be maintained by self-renewal unlimitedly. The possibility for geometrical regulation is included at each stage of osteoblast and osteoclast development through functions of the specific surface, namely , , and (see Figure 2). In the following, we denote the bone cell densities within a tissue sample (number of cells per unit volume) by their symbol ob, ob, ob, oc, oc, oc. The concentrations of the biochemical signalling molecules within a tissue sample (number of molecules per unit volume) is also denoted by their symbol tgf, rankl, etc.444To align with common practice, we shall use the terminology ‘density’ for cells and ‘concentration’ for signalling molecules, even if the units are chosen the same. Based on the above descriptions of the biochemical, biomechanical, and geometrical regulatory mechanisms, the governing equations of the bone cell densities in the model are expressed as:
Below, we discuss the various quantities occurring in these equations (see Refs [23, 22] for more details). The parameters , , , denote differentiation rate parameters for uncommitted osteoblast progenitors, pre-osteoblasts, uncommitted osteoclast progenitors and pre-osteoclasts, respectively; is a proliferation rate parameter for pre-osteoblasts; and are apoptosis rate parameters for active osteoblasts and active osteoclasts. Biochemical regulation of cell development is achieved through so-called “activator” and “repressor” functions of the biochemical signalling molecules tgf, mcsf and rankl as follows. The functions , , and are activator and repressor functions regulating osteoblast differentiation and osteoclast apoptosis based on the concentration of tgf. The function is an activator function regulating differentiation of ocs into ocs based on the concentration of macrophage colony stimulation factor (mcsf). The functions and are activator functions regulating osteoclast differentiation based on the concentration of rankl. For simplicity, we assume that . The equations governing the evolution of the biochemical signalling molecules tgf, mcsf, rankl, rank, opg and pth and the form of the activator and repressor functions are presented in Appendix A, see Eqs (15)–(22). Tables 1–3 in Appendix A list the dynamic quantities, the biochemical parameters and the biomechanical quantities of the model.
2.3.1 Mechanical feedback regulation
Biomechanical studies suggest that the strain energy density is an important quantity that determines bone adaptation to various mechanical loads. This quantity is commonly chosen in the literature due its scalar nature as a measure of mechanical stimulus sensed by the bone cells to drive bone adaptation [3, 32].
The model of biomechanical regulation developed in Ref.  similarly uses the strain energy density as the signal conveying mechanical information to the bone cells. However, as discussed in Ref. , it is paramount to estimate consistently the bone matrix strain energy density at the micro-scale where osteocytes sense this mechanical signal and transduce it into a biochemical response. In Ref. , we used a homogenisation procedure based on Eshelby’s classical matrix inclusion problem and the Mori-Tanaka scheme to estimate the microscopic strains generated at the level of the extravascular matrix hosting the osteocytes, by the macroscopic loading of the tissue. This homogenisation procedure leads to an intricate dependence of the microscopic strain energy density of the bone matrix, , upon the macroscopic stress tensor and the vascular porosity :
In this paper, we will exemplify the effect of biomechanical regulation for the situation of a constant uniaxial compressive loading only (i.e., with ). For this situation, the dependence of upon is plotted in Figure 4. The reader is referred to Refs [22, 14, 15, 33] for the general mathematical specification of and a detailed presentation of the homogenisation procedure.
As mentioned in Section 2.1, the biomechanical regulation of bone remodelling is believed to operate through different pathways for resorption and formation. Following Ref. , the biomechanical regulation of bone resorption is realised in the model via modulation of the rank–rankl–opg signalling pathway by the microscopic strain energy density . The production rate of rankl on pre-osteoblasts, , is assumed to be enhanced during mechanical disuse:
where is a parameter quantifying the strength of the biomechanical transduction and is the steady-state value of the strain energy density. The steady state is assumed to be an initial homeostatic state of bone remodelling, with no bone gain or loss. The increase in rankl production rate during mechanical disuse increases both and in Eqs (7) and (8), and consequently leads to increased osteoclast generation (see Appendix A, Eqs. (18), (22)).
Following Ref. , the biomechanical regulation of bone formation is realised in the model via modulation of pre-osteoblast proliferation by the microscopic strain energy density . In Eq. (5), pre-osteoblasts are generated both by differentiation from obs and by self-expansion through proliferation. The modulation of proliferation by the strain energy density is expressed by an ‘activator’ function , defined as:
where is a constant parameter quantifying the strength of the biomechanical transduction and, is the minimum value of the strain energy density for which .
2.3.2 Geometrical feedback regulation
The four regulatory functions , , and in Eqs (5)–(8) include a geometrical feedback at various stages of osteoblast and osteoclast development. This enables us to distinguish two types of geometrical action. Indeed, the modulation of the bone cell developmental stages by the specific surface can be interpreted in the model in terms of modulation of the initiation of new remodelling events (bmu creation) and modulation of resorption and formation within existing bmus, as explained in the following:
Initiation of new remodelling events. The initiation of a new remodelling event is a localised process that creates a new bmu. While the exact biochemical mechanisms that lead to the creation of a new bmu are poorly understood, it is believed that first steps in this process are the recruitment of pre-osteoclasts and pre-osteoblasts at the bone surface . This recruitment is thought to be controlled by osteocytes sensing the local mechanical state of the bone matrix and communicating with progenitor cells in the marrow through the bone surface. The complex dependence of bone surface availability on the recruitment of pre-osteoblasts and pre-osteoclasts is modelled by the geometrical regulation of cell differentiation exerted by the functions and . Therefore, the geometrical regulation by and models the influence of bone surface availability on the initiation of new remodelling events, and so on the number of bmus in the representative volume element, which in turn determines the bone turnover rate . This type of geometrical regulation of bone remodelling is similar to the geometrical regulation of the ‘activation frequency’ of bmus used by Hazelwood et al. .
Modulation of resorption and formation within existing bmus. Active osteoclasts can only resorb bone from the bone surface. Similarly, active osteoblasts are only observed to deposit new bone at the bone surface. The maturation of pre-osteoblasts and pre-osteoclasts into active cells thus depends on bone surface availability. This dependence is modelled by the geometrical regulation of cell differentiation exerted by the functions and . Therefore, these functions dictate how many active osteoclasts and active osteoblasts can form in bmus that are already remodelling bone (i.e., which contain already pre-osteoblasts and pre-osteoclasts), and so how much bone is resorbed and formed in existing bmus, which in turn determines bone balance. Note that while bone surface availability is necessary for cell activation, it is not sufficient. For instance, osteoclasts can only become active if, in addition, their receptor rank is activated by the ligand rankl.
Due to the different ways of action of bone surface on cell differentiation explained above, the regulatory functions , , and may take different forms and can be complicated functions of the specific surface . We take here a phenomenological approach and assume that each of these functions can be represented by a power law of :
In Eq. (12), denotes the specific surface in the bone remodelling steady state, which is assumed to be homeostatic (no bone gain or loss, but a steady bone turnover). The normalisation of by its steady-state value ensures that in the steady state, , and so that the steady state of the model is consistent with previous models of bone remodelling without geometrical regulation [23, 22]. The benefit of using a power-law function of in Eq. (12) is that geometrical feedback can be “switched off” by choosing in Eq. (12). In this situation, not only the steady state, but also the dynamical behaviour of the model of Ref.  (including mechanical feedback) is retrieved. Mechanical feedback is “switched off” by choosing , in Eqs (10) and (11). Note that when both geometrical and mechanical feedbacks are “switched off”, the bone cell population model of Refs [23, 24] is retrieved, except for a remaining pre-osteoblast proliferation term in Eq. (5) that was not accounted for previously.555Compared to Refs [23, 24], the differentiation rate of obs to obs is reduced accordingly to ensure that the model converges to the same steady state.
To reveal in which ways the morphology of the microstructure of bone may influence bone remodelling, we will investigate the effects of several combinations of , , and in Section 3 and determine combinations that lead to physiologically meaningful results.
2.4 Changes in porosity and bone matrix fraction due to cell activity
The activity of osteoclasts and osteoblasts leads to the removal and deposition of new bone. This activity modifies the volume fraction of bone matrix in the tissue. Osteoblasts deposit osteoid, a collagen-rich substance which later mineralises into new bone. Primary mineralisation of osteoid is relatively fast: 70% of the maximum mineral density is reached within a few days in humans . Given the much larger time spans involved in the remodelling processes, it is fair to model the osteoblasts “instantaneously” depositing “fully” mineralised new bone matrix, as was assumed in Refs [23, 24, 22]. We further assume that the resorption rate of bone matrix by an individual active osteoclast (in volume per unit time) is constant, and that the rate of new bone matrix deposition by an individual osteoblast (in volume per unit time) is constant. The evolution of the vascular porosity and bone matrix volume fraction are thus given by
2.5 Comparison with the model of geometrical regulation of bone remodelling by Martin
It is informative at this point to compare our formulation of geometrical regulation of bone remodelling with that proposed by Martin in Ref. . The evolution of the vascular porosity proposed by Martin is, in our notations (see [13, Eq. (67)]):666In Ref. , Martin uses the symbol to denote only the fraction of 3D specific surface that is capable of remodelling. This excludes surfaces of the lacunae-canaliculi system . In this paper, we do not consider surfaces of the lacunae-canaliculi system, and so in Ref.  is equal to our definition of .
where and are the active osteoblast and active osteoclast surface densities at an active bone site (number of cells per unit surface), and and are the fractions of the total available bone surface in which there is osteoblastic and osteoclastic activity. The quantity is the formation rate of bone matrix and corresponds to ob in our model. The quantity is the resorption rate of bone matrix and corresponds to in our model. To model an osteoporotic pathological condition, Martin assumes a constant bone imbalance between resorption and formation in remodelling bmus. This imbalance is modelled by setting /year in Eq. (14), meaning that in average, a 2 -thick layer of the bone surface is resorbed each year. The results presented in Figure 1 were obtained by Martin in Ref.  by integrating Eq. (14) with this constant imbalance and with either (no geometrical feedback) or using the function presented in Eq. (4) to include the effect of geometrical feedback.
A limitation of the formulation by Martin is that the dependence of the rate of change of porosity upon does not account for microscopic effects of bone surface availability on the recruitment and development of bone cells. Indeed, in Martin’s formulation, the complexity of osteoblast and osteoclast recruitment and development at a remodelling bone surface is embodied by the fractions and , but these fractions are assumed to be independent of . Martin and others have interpreted Eq. (14) as the influence of geometrical regulation on the activation frequency of bmus [13, 2]. Indeed, the densities of active osteoblasts and osteoclasts in Martin’s model are proportional to the specific surface. Comparing Eqs (13) and (14), one sees that and . Thus, an increased specific surface represents increased numbers of remodelling obs and ocs, and so an increased number of bmus.
The strength of our model is to consider several stages in the recruitment and development of osteoblasts and osteoclasts explicitly. This enables us to include a more detailed geometrical feedback acting directly on these different stages of the remodelling sequence. In our model, even steady-state values of obs and ocs are complex functions of the specific surface. This complexity represent the implicit dependence of and upon not accounted for by Martin. Additionally, no mechanical feedback has been considered by Martin. In Ref. , Scheiner et al. have discussed the problem of neglecting mechanical feedback in models of bone remodelling applied to catabolic bone diseases (such as osteoporosis), which leads to a unlimited bone loss which is physiologically not observed. This behaviour can be seen in Figure 1 for cortical bone where the vascular porosity continues to increase (both with and without geometrical regulation). It is known from bone mineral density measurements in osteoporotic patients that the increase in bone porosity in osteoporosis slows down with time. The porosity eventually reaches an upper bound whose value depends on the patient. We note here that the long-term value of the vascular porosity reached with mechanical feedback in our model depends in particular on the initial value of (see Section 3).
3 Numerical Results and discussion
Having incorporated biochemical, biomechanical and geometrical regulation of bone remodelling we are now in a position to investigate the effects of various regulatory parameters on changes in bone cell numbers, vascular porosity and bone stiffness. To compare our model with the model suggested by Martin in Ref. , we simulate an underlying osteoporotic condition as a perturbation from the original (homeostatic) steady state situation.
Osteoporosis is a bone disease that leads to an increase in porosity in both cortical and trabecular bone . This increase in porosity generates a progressive reduction in bone stiffness and a higher fracture risk. To simulate an osteoporotic condition in our model, we perturb the homeostatic steady state (a state with no bone loss or gain and a constant bone turnover rate) by increasing the rankl/opg ratio. This can be achieved in our model by prescribing an excess of pth concentration (see Refs. [25, 23]).777The physiological effect of pth administration for bone remodelling is complex and in particular, depends on the time course of the administration . Continuous pth administration (infusion) leads to bone loss with increased turnover. However, intermittent pth administration (daily injections) leads to bone gain [34, 35]. Only the ‘continuous’ action of pth is represented in our model [25, 23]. We have shown previously that increasing pth leads to an increase in the rankl/opg ratio, and so to bone loss with higher turnover rate compared to the homeostatic bone remodelling state [24, 28]. Increasing pth concentration thus consists in an adequate representation of osteoporosis in the present model. We also note that the effect of pth in bone remodelling is known to interact synergistically with mechanical loading [36, 37]. Consistently with the interpretation of increasing pth as a representation of osteoporosis in our model, we do not consider this synergistic interaction here. The osteoporotic condition is assumed to develop instantly from an initial homeostatic state at time (corresponding to a middle biological age). To obtain a steady increase in of 0.01/year without geometrical regulation, as has been assumed by Martin (see Figure 1), a continuous pth administration rate has been applied at all times . We denote by the time elapsed since the onset of the osteoporotic condition. The evolution of the system is followed for a period of time of 20 years from the onset of osteoporosis, i.e., . The initial values for vascular porosity of cortical and trabecular bone have been chosen as and . For the simulations including biomechanical regulation, a uniaxial compressive stress with has been assumed to exert on the representative volume element.
3.1 Simulation of osteoporosis: Evolution of bone porosity and bone stiffness properties
Figure 5 shows the simulated evolution of bone cell densities, vascular porosity () and selected components of the bone stiffness matrix (i.e., axial stiffness and shear modulus ) in osteoporosis considering no geometrical feedback and no mechanical feedback. The evolution of bone cell densities exhibits a short transient for a period of days after the onset of osteoporosis (Figure 5a). After this initial transient, ocs and obs reach a new steady state, in which resorption and formation are imbalanced and lead to osteoporotic bone loss. Note that a short time interval has been chosen in Figure 5a to demonstrate the short transient cellular response, while a larger time interval of 20 years is chosen to follow the evolution of the vascular porosity, axial stiffness and shear modulus. The progressive increase in vascular porosity is shown in Figure 5b together with the associated reduction in bone stiffness. After 20 years of osteoporosis, the normal stiffness in the longitudinal direction of cortical bone, , has dropped by about 20%, whereas the shear modulus of cortical bone, , has dropped by about 40%. From these results it is clear that the evolution of vascular porosity drives the changes in bone stiffness differently for the different components of the stiffness tensor . Generally speaking, an increase in is always associated with a reduction of bone stiffness properties. For conciseness, in the following we will only present the effects of geometrical and biomechanical regulation on the vascular porosity alone.
3.2 Geometrical regulation—effect on individual stages of osteoblast and osteoclast developments
We first present how geometrical regulation influences the osteoporotic increase in bone porosity without accounting for mechanical feedback, i.e., by setting and in Eqs (10), (11). Based on the four regulatory functions , , , and accounting for geometrical regulation, we first performed simulations ‘switching on’ or ‘off’ each of these regulatory functions to investigate all different possible combinations. However, it turns out that among all possible simulations there are only three patterns which reflect the bone systems behaviour. These patterns can best be studied by looking first at the cases in which only a single regulatory function is ‘switched on’ while all others are ‘switched off’ (i.e., identically set to one) (Figure 6). We will present a selection of combinations of geometrical regulation on several stages of osteoblast and osteoclast development in Figures 7 and 8.
Figure 6a shows the influence of on the increase of vascular porosity in osteoporosis compared to the case of no geometrical and no biomechanical regulations in both cortical and trabecular bone. The time evolution of the porosity in Figure 6a clearly shows that while reduces bone loss in cortical bone it accelerates bone loss in trabecular bone. The strength of this regulation depends on the exponent of the normalised geometrical regulatory function in Eq. (12). Exponents in the range exhibit a relatively strong regulatory effect on the evolution of the disease, while exponents in the range only exhibit a moderate effect. Interestingly, the mechanism of action of is opposite to the geometrical regulation obtained by Martin  (see Figure 1). In cortical bone, the osteoporotic increase in induces an increase in (see Figure 4), and so an increase of . This increases in turn the generation of osteoblasts, which has an anabolic (i.e., bone forming) effect and stabilises the osteoporotic loss of bone [25, 23]. By contrast, in trabecular bone the osteoporotic increase in induces a decrease in (see Figure 4), and so a decrease of . This decreases in turn the generation of osteoblasts, which accelerates the osteoporotic loss of bone.
While this scenario can be understood from a theoretical point of view, it has to be emphasised that it is probably physiologically unrealistic. As argued in Section 2.3.2, the geometrical regulation of ob differentiation is associated with the creation of a new remodelling event, i.e., of a new bmu. However, both ob differentiation and oc differentiation are activated in such an event. Here, the consideration of a geometric regulation on ob differentiation alone represents the initiation of a formation event and is not associated with a joint initiation of a resorption event. To represent the geometrical regulation of bmu creation, both and should act together (see Figures 7 and 8).
Figure 6b shows the influence of on the increase of in osteoporosis. In this figure, is only weakly affected by the geometrical feedback in both cortical and trabecular bone. The effect of bone surface availability on the differentiation of pre-osteoblasts into active osteoblasts modelled by does not seem to affect the bone remodelling balance significantly in the first 10 years. In fact, neither the population of active osteoblasts nor the population of active osteoclasts is affected significantly, while the population of pre-osteoblasts is decreased in cortical bone and increased in trabecular bone (not shown). The increased differentiation rate of pre-osteoblasts into active osteoblasts in cortical bone (due to the increase in the available bone surface) depletes the population of pre-osteoblasts. The population of active osteoblasts is thus derived from a smaller pool of pre-osteoblasts that are differentiating faster, and therefore stays relatively constant. The opposite effect occurs in trabecular bone, i.e., a larger pool of pre-osteoblasts is created but they are differentiating into active osteoblasts more slowly.
Figure 6c shows the influence of on the increase of in osteoporosis. The behaviour of obtained for this type of geometrical regulation is in general agreement with the geometrical regulation obtained by Martin  (see Figure 1), i.e., geometrical regulation leads to increased bone loss in cortical bone and decreased bone loss in trabecular bone. In cortical bone, the increase in due to the osteoporotic condition leads to an increase in osteoclastogenesis, which by biochemical coupling also increases (although to a lesser extent) the population of osteoblasts. This results in a high turnover rate with a catabolic bias, i.e., a high rate of bone loss. In trabecular bone, the situation is reversed since decreases with the evolution of the osteoporotic condition. A state with lower turnover rate and a reduction in resorption is reached. The low turnover rate explains the reduced amplitude of the response in trabecular bone compared to that seen in cortical bone.
The strength of the geometrical regulation is determined by the value of and is seen to strongly affect bone balance. In particular for the case of strong modulation, i.e., , a rapid acceleration of cortical bone loss occurs over the first 10 years of osteoporosis. At nearly 10 years, a transition is taking place leading to a reduction in the rate of bone loss. This behaviour is due to the fact that at this time, the bone porosity reaches the critical value at which the specific surface is maximum (see Figure 3). In the first 10 years of osteoporosis, is in the ascending branch of and so increases, while after reaching the maximum specific surface, decreases, which leads to a reduction in the rate of bone loss. Over 20 years, the original cortical bone has been resorbed enough to reach trabecular porosities. Such a strong effect of bone remodelling on bone porosity is normally not physiologically seen. However, it is possible that in osteoporosis, a locally strong geometrical feedback may help initiate or accentuate the observed ‘trabecularisation’ of bone, by which cortical bone is progressively lost and transformed into trabecular bone at the endocortical wall, leading to cortical wall thinning and expansion of the medullary cavity [38, 39, 19, 40]. Indeed, in bone, the endocortical wall exhibits the highest specific surface and is known to be highly remodelling. It can be expected that geometrical regulation plays a particularly significant role in this region of bone.
Figure 6d shows the influence of on the increase of in osteoporosis. This influence is qualitatively similar to that of shown in Figure 6a, i.e., a reduction in bone loss in cortical bone and an increase in bone loss in trabecular bone. However, the geometrical regulation modelled by is less pronounced than the geometrical regulation modelled by .
3.3 Geometrical regulation—combined effect on several stages of osteoblast and osteoclast developments
As mentioned previously, a geometrical regulation of the creation of new remodelling events (i.e. new bmus) should involve a regulation of both the recruitment of osteoclasts and the recruitment of osteoblasts (see Section 2.3.2). In Figure 7a, we show the combined influence of both and on the increase of in osteoporosis. The exponents and measuring the strength of the geometrical regulation are assumed identical. By comparing with the individual influences of and in Figure 7a and 7c, one sees that the geometrical regulation of osteoblast recruitment overrides that of osteoclast recruitment. As a consequence, the overall behaviour is opposite to that obtained by Martin  for the geometrical regulation of the creation of new bmus (see Figure 1).
By contrast, in Figure 7b, the geometrical regulation of the last stage of osteoclast differentiation by and the geometrical regulation by of the last stage of osteoblast differentiation (modelling a regulation of the activation of cells within existing bmus) seem to compensate each other and to result in an evolution of the porosity that is almost unaffected by geometrical feedback. The geometrical regulations assumed in Figures 7a and 7b represent different natures of the influence of bone surface availability on bone remodelling (see Section 2). Therefore, our simulations suggest that for the simulated evolution of bone porosity in osteoporosis, the influence of surface availability is significantly stronger on the creation of new remodelling events (new bmus) (Figure 7a) than on the activation of bone cells within already active remodelling sites (within existing bmus).
In Figure 8, we show that a similar influence of geometrical regulation of bmu creation as that obtained by Martin  can be retrieved in our model by modifying the relative strengths of the regulatory functions and via the exponents and . The vascular porosity can exhibit a wide range of behaviours, interpolating between the situation of Figure 6a (, ) and the situation of Figure 6c (, ). We conclude that geometrical feedback has the potential to significantly influence the evolution of bone diseases. However, complex biochemical coupling between osteoblasts and osteoclasts makes it difficult to predict the relative strength of the influence of geometrical regulation on osteoblast and osteoclast developments.
Finally, we note that because the curve exhibits a maximum at (see Figure 3), geometrical feedback always induces opposite behaviours for the evolution of vascular porosity in cortical bone (for which ) and in trabecular bone (for which ).
3.4 Coupled geometrical and mechanical regulations
Figure 9 shows the effect of adding a mechanical regulation of bone remodelling for the evolution of the vascular porosity. The geometrical regulation considered in Figure 9 (see dotted line) is assumed to represent the influence of bone surface availability for bmu creation as in Figure 8, solid line. Two strenghts of the biomechanical transduction are illustrated (i.e., and in Eq. (11)). Figure 9 suggests that mechanical feedback has the potential to progressively override both the evolution of osteoporosis and the influence of geometrical regulation. Indeed, mechanical feedback is seen to stabilise bone balance, irrespective of whether geometrical feedback stabilises or destabilises bone balance in osteoporosis. Such a stabilisation of bone loss is clinically observed in osteoporotic patients.
Similarly to the simulations by Scheiner et al. , mechanical feedback counteracts bone loss in osteoporosis both in cortical and trabecular bone, on a time scale that depends on the strength of mechanical regulation (i.e., on the parameter ). Interestingly, in our simulations with geometrical feedback, the strength of the mechanical regulation also has an influence on the steady-state value of the porosity attained. This is to be contrasted to the situation without geometrical feedback in which no such influence was observed . In fact, without geometrical feedback, the vascular porosity enters the right hand side of the governing equations (5)–(8), (13) only implicitly via the strain energy density , see Eq. (9). Consequently, the steady-state value of is uniquely determined by the steady-state value of the strain energy density, , and the macroscopic loading , which are themselves independent of the strength of biomechanical regulation . With geometrical feedback, additional dependences upon the vascular porosity are added in the governing equations for the bone cells, Eqs (5)–(8). Consequently, the steady-state value of now depends in addition on the biochemical and cellular state of the system, which depends in turn on the strength of the biomechanical regulation of the bone cells.
From the above discussion, it is clear that the coupling between the geometrical and mechanical feedbacks is effectively mediated by the biochemical and cellular state of the system. In osteoporotic patients, both a geometrical and a mechanical feedback are present at the same time, as well as biochemical and hormonal dysregulations underlying the establishment of osteoporosis. The interdependence between the biochemical and hormonal dysregulations in osteoporosis, geometrical feedback and mechanical feedback, is not trivial to elucidate. Physiologically, it is expected that all these influences play a role for the evolution of bone vascular porosity. Our mathematical model is a first attempt to integrate these influences and help understand their contribution for clinically observed changes of bone porosities in osteoporotic patients.
In this paper we developed a computational model of bone remodelling that takes into account biochemical, biomechanical and geometrical regulations of bone cells. The biochemical regulation of the bone cells is based on the model developed in Refs [23, 24]. The biomechanical regulation of the bone cells is based on the model developed in Ref. , in which a continuum micromechanical approach is used to consistently link bone cell responses with mechanical properties of bone. The new contribution of this paper is the inclusion of a geometrical regulation of bone cells in the model. Geometrical feedbacks were included at several developmental stages of osteoblasts and osteoclasts to represent the influence of microscopic bone surface availability for various bone microstructures. We investigated the influence of a geometrical regulation of bone cells in bone remodelling both without and with consideration of biomechanical regulation for a simulated osteoporotic condition. From the numerical simulations, we identified the following actions of geometrical and mechanical regulation on bone remodelling:
Geometrical regulation of bone remodelling may play an important role for the initiation of new bmus as described by the combined effect of the geometrical regulatory functions and acting on ob differentiation and oc differentiation; in particular geometrical regulation of ocs via seems to be most important to retrieve similar evolutions of bone porosities in osteoporosis as obtained by Martin ;
Geometrical regulation of bmu creation affects cortical and trabecular bone in opposite ways in osteoporosis: while bone resorption is enhanced in cortical bone due to the increase in specific surface with increasing porosity, bone resorption is slowed down in trabecular bone due to the decrease in specific surface with increasing porosity;
Geometrical regulation of the activation of osteoblasts and osteoclasts in existing bmus seems to play a secondary role for the evolution of osteoporosis. While the specific surface can influence the differentiation of bone precursor cells into active resorbing/forming cells, no significant influence on bone porosity was observed in our simulations.
Our simulations suggest that geometrical regulation may play a role in osteoporosis for the initiation and/or accentuation of the observed ‘trabecularisation’ of bone at the endocortical wall. At the endocortical wall, the specific surface and bone turnover are high and so the effects of a geometrical feedback can be expected to be significant;
Our simulations of coupled geometrical and mechanical regulations suggest that the stabilisation of bone loss observed clinically in osteoporotic patients is probably accelerated by geometrical feedback in trabecular bone, but is probably slowed down by geometrical feedback in cortical bone.
Both mechanical and geometrical feedbacks are important to account for in our model of bone remodelling. Mechanical feedback enables the local porosity of bone tissue to stabilise to a well-defined value within . Geometrical feedback enables this value to be not only determined by the external loading of the tissue, but also by the biochemical and cellular state of the system, as would be expected physiologically.
Financial support by the Australian Research Council (ARC), in the framework of the project Multiscale modelling of transport through deformable porous materials (project number DP-0988427) and by the European Research Council (ERC) in the framework of the project Multiscale poro-micromechanics of bone materials, with links to biology and medecine (project number FP7-257023) are gratefully acknowledged.
Appendix A Model description
In this appendix, we complete the mathematical description of the model and give the values of the parameters. The nomenclature used in the paper is split into three tables. Table 1 lists dynamics quantities involved in the governing equations of the bone cells and bone porosity, Eqs (5)–(8), (13). Table 2 lists all the parameters relevant to the biochemical regulation of the model. Table 3 lists quantities involved in the biomechanical regulation of the model.
Concentrations of the biochemical signalling molecules
As in Ref. [23, 24, 22, 26], the concentration of the biochemical signalling molecules are governed by rate equations based on mass action kinetics. Ligand–receptor binding reactions occur on a time scale much faster than the characteristic times of cellular response (such as differentiation, apoptosis). The rate equations for the biochemical signalling molecules can therefore be taken in their steady state (see Refs. [23, 26] for details). This gives:
In Eq. (18), is the maximum concentration of rankl (also referred to as effective carrying capacity), which is regulated by the parathyroid hormone pth:
where is the maximum number of membrane-bound rankl molecules that can be expressed on a single pre-osteoblast (see Refs. [23, 24] for more details). In this work, the concentration of mcsf, and so the quantity , are assumed constant (see below and Table 2). The additional production rate of pth in Eq. (19), , is used to increase the normal systemic levels of pth to simulate an osteoporotic condition (see comments in Section 3).
Activator and repressor functions
The regulation of the bone cell behaviours (such as differentiation, apoptosis, expression rate of a ligand) by the biochemical signalling molecules is modelled in Eqs (5)–(8) by so-called “activator” functions and “repressor” functions , where denotes the signalling molecule (ligand) and the signalled cell. These activator and repressor functions represent the strength of the response of the cell to the signal mediated by the ligand and are assumed to be given by:
where is the dissociation binding constant between the ligand and its receptor on the cell. The quantity represents the fraction of the receptors on the cell that are bound to a ligand (see Refs. [25, 23] for more details). For example, the activator functions and are defined as:
where is the dissociation binding constant between rankl and the rank receptor on ocs and ocs, and rankl is the free (unbound) rankl concentration given by Eq. (18).
The geometrical and biomechanical regulations of the bone cells is normalised by the steady-state values. These values depend in particular on the initial bone porosity, and so are calculated prior to evolving the system.
|oc||pM||density of pre-osteoclasts|
|oc||pM||density of active osteoclats|
|ob||pM||density of pre-osteoblasts|
|ob||pM||density of active osteoblasts|
|tgf||pM||concentration of transforming growth factor|
|rank||pM||concentration of receptor-activator nuclear factor B|
|rankl||pM||concentration of receptor-activator nuclear factor B ligand|
|opg||pM||concentration of osteoprotegerin|
|pth||pM||concentration of parathyroid hormone|
|–||volume fraction of vascular pores|
|MPa||microscopic strain energy density of the bone matrix|
|oc||pM||density of uncommitted osteoclast progenitors|
|ob||pM||density of uncommitted osteoblast progenitors|
|200||daily volume of bone matrix resorbed per osteoclast|
|40||daily volume of bone matrix formed per osteoblast|
|differentiation rate parameter|
|differentiation rate parameter|
|oc apoptosis rate parameter|
|differentiation rate parameter|
|differentiation rate parameter|
|ob proliferation rate parameter|
|ob apoptosis rate|
|0.5||value of the activator function of mcsf for differentiation|
|5.68 pM||dissociation binding constant for rankl binding on oc and oc|
|pM||dissociation binding constant for tgf binding on oc|
|pM||dissociation binding constant for tgf binding on ob|
|pM||dissociation binding constant for tgf binding on ob|
|150 pM||dissociation binding constant for pth binding on ob (in )|
|0.222 pM||dissociation binding constant for pth binding on ob (in )|
|0.034/pM||association binding constant for rankl and rank|
|0.001/pM||association binding constant for rankl and opg|
|production rate of rankl|
|production rate of systemic pth|
|continous administration rate of pth to model osteoporosis|
|production rate of opg per ob|
|maximum number of rankl per ob|
|number of rank receptors per oc|
|opg density at which endogeneous production stops|
|density of tgf stored in the bone matrix|
|degradation rate of tgf|
|degradation rate of rankl|
|degradation rate of opg|
|degradation rate of pth|
|macroscopic stress tensor|
|microscopic strain tensor of the bone matrix|
|microscopic strain energy density of the bone matrix|
|strength of biomechanical transduction of bone resorption|