Cell Interaction and Mechanobiological Modeling of Bone Remodeling Process

Abstract

According to the structural and metabolic demands of the body, proportionate and accurate bone quantities are resorbed and formed, establishing what is known as bone remodeling process. This physiological process requires a highly coordinated regulation through a complex interconnected network involving several cells from diverse origins, in addition to various hormones, cytokines, growth factors and signaling pathways. One of the main factors initiating the remodeling process is the mechanotransduction mechanism, through which osteocytes translate the mechanical stimuli subjected to the bone into biochemical signals, generating thereby the activation of osteoclasts and osteoblasts that govern bone resorption and formation. This mechanically-induced behavior of bone tissue has been the target of computational modeling and numerical simulations, to address biomechanical questions and provide information that is not amenable to direct measurements. In this context, the current chapter aims to review the coupling and mechanotransduction mechanisms spearheading the remodeling process, in addition to the main mathematical models developed over recent years and their use in bone numerical simulations based on the finite element method.

1. Introduction

Bone is a living composite material, providing various mechanical and homeostatic functions and characterized by a constant adaptation of its structure according to the metabolic and physical demands of the body. This self-adaptation property is governed by the complementary activities of resorption and formation, establishing what is known as bone remodeling process. This process requires a highly coordinated regulation over time and space to consistently maintain the proper quality and quantity of bone forming the skeleton. This coordination mainly incorporates two types of cells: bone-resorbing osteoclasts and bone-forming osteoblasts, representing the two major players of the remodeling event and providing a delicate balance between the removed bone amount and the subsequent deposited amount, which is carried out by generating the appropriate number of osteoblasts, a phenomenon known as the coupling mechanism. The coordination between osteoblast and osteoclast activities also involves other cells of various origins, in addition to several hormones, cytokines and growth factors that tightly connect the osteoblast- and osteoclast-lineages through a complex communication network during the remodeling cycle.

The remodeling process maintains the calcium homeostasis and provides a crucial mechanism for old bone removal, as well as for damage repair and adaptation to physical stress, which helps preserving the mechanical integrity of the skeleton. The remodeling process takes place at an anatomically distinct sites called bone mass units (BMUs), with each BMU operating asynchronously and independently of other BMUs throughout the skeleton. Notably, BMUs in the cortical bone greatly differ from those in the trabecular bone, in terms of structure, but also of the resorption and formation activities. Bone remodeling is initiated in a canopy, defining what is called the bone remodeling compartment (BRC), where intercellular communication occurs among the component bone cells, from vascular and endothelial cells and probably from immune cells reaching the remodeling sites through the blood supply. To better understand bone behavior and biomechanics, several research have been conducted using clinical investigations, as well as numerical modeling.

Computational modeling and numerical simulation of represent an interesting tool to address biomechanical questions, particularly those targeting the biomechanical behavior of bone, allowing to provide information that is not amenable to direct measurements, such as bone strength and joint load. This kind of information is required for several clinical applications, including fracture prevention, implant design, and pathology analysis. In recent years, partition of unity methods, explicitly using finite element (FE) mesh, has become popular due to its easy applicability. The FE method is one of the most widely used numerical analysis techniques based on FE mesh. It provides approximate solutions to a wide range of engineering problems. Recently, the FE method has experienced a phenomenal expansion in the field of bone biomedical engineering, owing to its flexibility and diversity as an analysis tool, but also because it easily manages complex and evolving cellular domains, and can be generalized to multidimensions with little complication. Particularly, the FE method has been widely used to analyze and predict the mechanical behavior of bones, under physiologic and pathologic conditions, based on mathematical models that describe the cellular mechanisms governing the remodeling process. These mathematical models are including more and more factors and actors, to attempt affording a more realistic description of bone cell interactions.

In this context, the current paper provides a review of the processes of bone remodeling and cell interactions, as well as the recent findings about the mechanotransduction mechanisms, in addition to the mechanobiological models targeting bone dynamics and the main generated FE results in recent years.

2. Bone cells and coupling mechanisms

Osteoblasts are cuboidal cells located on the newly synthesized bone interface and have two types of functions: a bone building function, through which they produce matrix proteins, and endocrine functions, through which they release a wide range of regulatory factors that influence energy metabolism, male fertility and cognition [1, 2]. Osteoblasts originate from mesenchymal stromal cells (MSCs) according to two distinct embryonic populations: the first one is when osteoblast-lineage cells derive from the neural ectoderm, the mesenchymal progenitors directly differentiate into preosteoblasts and the subsequent mature osteoblasts form the calvarian bones and clavicles through intramembranous ossification, whereas the second one is when the MSCs differentiate into perichondral cells and chondrocytes, a subsequent chondrocyte hypertrophy drives the perichondral cell differentiation into preosteoblasts, and mature osteoblasts form bone in the extremities and in the axial skeleton through endochondral ossification [3, 4, 5].

Osteoclasts are giant multinucleated cells that hydrolyze and solubilize both of the inorganic and organic components of bone, owing to their polarized secretion of acid and proteolytic enzymes. Thus, osteoclasts provide bone with a unique characteristic of being the only tissue in the body able to undergo a self-destruction, fulfilling thereby a crucial physiologic process for bone homeostasis [6]. Osteoclasts originate from the fusion of monocyte–macrophage [7] precursors that derive from hematopoietic stem cells (HSCs) in the marrow [8]. The hematopoietic precursors required for osteoclast formation are provided through the capillary blood supply closely associated with and penetrating the BRC [9], as well as from nearby marrow precursors [10]. Interestingly, the programming of osteoclast formation involve the actions of osteoblast-lineage cells, endothelial cells, and BRC microenvironment [9].

Osteocytes are the most abundant and the dominant mechanosensory cells of bone and are easily identified in a bone section, owing to the matrix in which they are completely embedded during skeletal maturation of previous remodeling cycles [1, 11]. These long-living cells are characterized by a stellar/spider shape, due to which they form and extensive interconnected network of long branched cellular processes within a fluid-filled canalicular system [12]. These processes contact each other, and probably other cell populations through gap junctions, providing thereby a cell–cell interaction by an intercellular exchange of small signaling molecules. This network plays a crucial role in regulating bone material turnover, by coordinating bone response to mechanical signals, as well as to endocrine and paracrine biological signals [13, 14]. Osteocytes have both local and systemic effects through responding to the strains generated by mechanical stimuli and translating the load into biochemical signals via a mechanotransduction mechanism [11, 15], which is one of the main factors initiating the remodeling process. Osteocytes emerge from the differentiation of a subset of osteoblasts that are trapped in the newly synthetized osteoid matrix before its mineralization [16], by undergoing a four-stage differentiation process to become mature osteocytes: (i) type I preosteocytes, named osteoblastic osteocytes, (ii) type II preosteocytes, named osteoid-osteocytes, (iii) type III preosteocytes, named young osteocytes, and (iv) old osteocytes [17].

Reversal cells were found to belong to osteoblast-lineage cells that seem to particularly be preosteoblasts [18] that progressively mature into bone-forming osteoblasts and play a crucial role in the resorption-to-formation coupling mechanism [19, 20]. The intermediate position of the reversal cells between osteoblasts and osteoclasts suggests their obvious contribution to the osteoblast–osteoclast interplay [21]. Interestingly, the early reversal cells located proximal to osteoclasts are less mature than the late ones located proximal to osteoblasts [18] and appear to be morphologically, ultrastructurally and immunohistochemically quite different from the late ones located next to osteoid surfaces. These differences reflect their diverse cellular interactions, varied functions, and distinct differentiation states [18, 19, 20].

Bone lining cells are an osteoblast subpopulation and their interaction with osteoclasts attached to bone, due to close contact, is tightly associated with the initiation of osteoclastogenesis [22, 23]. Besides, the lining cells establish physical homotypic connections with osteocytes through gap functions, which suggests that these lining cells form a functional membrane separating bone from interstitial fluids [24]. Particularly during the remodeling cycle, bone lining cells persist over the remodeling sites to isolate osteoblasts and osteoclasts from the bone marrow [25]. Indeed, at the very beginning of a remodeling cycle, the lining cells separate and raise from the underlying bone surface, forming a canopy above the site to be resorbed. This initiating phenomenon might result from osteocyte signaling to surface cells through their canaliculae when they recognize the need for replacing a specific bone area [26]. Subsequent signals could come from osteocyte apoptosis or from the lining cells themselves, generating the release of paracrine factors and chemokines, which attract the precursors of both osteoblasts and osteoclasts, as well as other required vascular elements [27].

3. Coupling mechanism

The spatial and temporal arrangement of bone cells within the BMU is crucial to the remodeling process, ensuring its distinct and sequential phase coordination, under the control of several hormones, cytokines, and growth factors. The communication between bone cells is led through at least three modes: (i) a communication through a direct cell–cell contact, (ii) a communication through gap junction formation, and (iii) a communication through diffusible paracrine factors.

Indeed, osteoblasts regulate osteoclast formation, differentiation, and maturation. In turn, osteoclasts exert positive and negative regulatory effects on osteoblast activity. Moreover, osteocytes, reversal cells, lining cells and bone matrix are all actors in the osteoblast–osteoclast crosstalk, among others. The interplay between osteoblasts and osteoclasts is also mediated by bidirectional signaling pathways, including Eph/ephrin pathways, Semaphorins/Plexins pathways, as well as RANKL/W9/RANK pathway. Figure 1 summarizes the communication network between bone cells.

Figure 1.

Intercellular communication pathways within the BMU whereby the remodeling process takes place. (1) mechanical stimulus affecting OCY (osteocytes). (2) stimulatory effects from OCY to act OB (active osteoblasts). (3) bone matrix signals to act OB. (4) signaling within osteoblast-lineage cells. (5) Stimulatory and inhibitory signaling from act OB to OCP (osteoclast precursors), OC (osteoclasts), and act OC (activated osteoclasts). (6) bidirectional signaling between act OB and act OC. (7) Stimulatory and inhibitory signals from OCP to OBP (osteoblast precursors). (8) stimulatory signals from late RvCs (reversal cells) to act OB. (9) Stimulatory signals from early RvCs to act OC.

4. Mechanotransduction mechanisms

As mentioned above, osteocytes are the dominant mechanosensory cells in bone, owing to their location in the matrix and their complex dendritic network. However, the lining cells were also suggested to be important mechanosensory cells in the adult skeleton, which thereby needs more investigation to clarify their exact role in bone biology. Unlike in vivo loading, where the ideal intensity, frequency, and timing to increase bone mass [28, 29, 30] are well characterized, in vitro experiments to identify these parameters and replicate in vivo results remain a major challenge. The fact that in vivo osteocyte gene expression changes under mechanical loading and unloading shows that the load affects osteocyte function [31, 32, 33, 34].

Early in vitro experiments were based on hydrostatic pressure and substrate stretching, whereas current ones are based on fluid flow shear stress (FFSS), owing to the higher sensitivity of primary osteocytes to shear stress than to substrate stretching [35, 36].

Mechanosensation ability allows osteocytes to orchestrate osteoblast activation and osteoclast partial suppression under increased loading, and osteoblast partial suppression and osteoclast activation under reduced loading [37]. However, the way external forces are transmitted at the cellular and molecular levels is still unclear. Mechanical stressors include hydrostatic pressure, FFSS, and direct cellular deformation [38]. These mechanical stresses are driven by microstrain of bone matrix generated by loading and gravitational forces. Cell responses are also influenced by the specific components of these stressors, such as amplitude, frequency and rate.

Although the calcified bone matrix is a mechanically rigid material, mechanical loading induces poro-elastic interactions and microstrains of the matrix, of up to 0.2% [39, 40]. These microstrains drive the interstitial fluid flow (IFF) within the lacunocanalicular spaces [41, 42]. Indeed, real-time measurement of load-induced solute transport has been demonstrated, which suggested a peak shear stress of 5 Pa on osteocyte processes [43]. Loading of long bones also increases the pressure at the intramedullary cavity and induces IFF at the endosteal surface, as well as within the lacunocanalicular network [41]. Moreover, intramedullary pressurization-derived IFF can induce fluid shear stress-related responses not only in osteocytes but also in osteoblasts and osteoclasts on the endosteal surface, which suggests that osteoblast and osteoclasts are also mechanosensitive [44, 45]. Several research works have shown that osteocytes are connected to the canalicular wall through transverse tethering elements and transmembrane molecules that provide physical connections between the extracellular matrix, the intracellular protein complexes and cytoskeletal structures [46, 47, 48, 49].

Several mechanosensors have been identified in bone, including cilia, integrins, calcium channels and G-protein coupled receptors (GPCRs). Integrins, among which αvβ3 is highly expressed in osteocytes, are composed of an α and β dimer and FFSS generates conformational changes in the β-subunit and the cascade signaling activation [46, 47]. The primary cilium is a non-motile structure needed for both mechano- and chemosensation, and represents another cellular moiety required to perceive FFSS. According to physical laws, FFSS occurs around the cell processes and not on the cell body where the primary cilium is located. The latter was proposed to perceive hydrostatic pressure applied on the cell body and not FFSS [50]. TAZ/YAP was also found to be an important signaling pathway for mechanosensation. Furthermore, glycocalyces on the dendritic process surfaces, but not on the cell body, were found to play a crucial role in mechanotransduction, whereas a different mechanosensing mechanism is active on the cell body [51]. Moreover, osteocytes sense load through cilia, which are single flagellar-like structures found on every cell [52, 53] and have specific functions. Particularly, cilia in bone cells induce the PGE2 release [53]. It was also reported that polycystin 1 (PC-1) in osteocytes is important for the anabolic response of bone to load [54]. Indeed, applying 2000 microstrain to bone sample at the macroscopic scale induced over 30,000 microstrain surrounding the osteocyte lacunae [55]. Besides, osteocyte processes are extremely responsive to mechanical loadings of piconewton-level, which is not the case for their cell body and processes with no local attachments [56].

In osteocytes, among other cells, biophysical stressors are transmitted to the cells by coupling the extracellular matrix to the actin cytoskeleton through focal adhesions. The actin cytoskeleton transmits mechanical forces from a focal adhesion site to a mechanosensing site within the cell and to the neighboring ones. Focal adhesion kinases (FAK) are major components of focal adhesions and are required for the mechanotransduction mechanism by osteocytes. Besides, a structural cytoskeletal protein named spectrin is required for osteoblast-to-osteocyte differentiation and was recently identified as a mechanosensitive element within the osteocyte [11, 57]. Other potential mechanosensors are ephrins, Connexin 43 (Cx43) hemichannels and ion channels, as well as gap junctions. A response of bone to mechanical loading and unloading also requires the action of an intact axis of parathyroid hormone (PTH)-related peptide (PTHrP) and its receptor (PTHR), or PTH-PTHrP-PTHR axis [58].

Mechanotransduction is the mechanism of transducing the mechanical signal sensed by osteocytes via the mechanosensation mechanism into biological cues. Ca++, ATP, NO, PGE2 and Wnts are the best described mechanically-induced pathways. Deleting one of these molecules inhibits the anabolic response of bone to mechanical stimulation.

Ca++ is an exclusively intracellular signal and the opening of stretch-activated calcium channel (TRPV6) results in a quick increase in intracellular Ca++, subsequently followed by cellular response to mechanical cues. Ca++ is also required for ATP response and ATP concentration rapidly increases upon mechanical forces. However, the exact mechanism by which Ca++ regulates ATP release is still not completely elucidated [49]. ATP can be released from osteocytes in response to whether mechanical stimulation or extracellular calcium [59, 60]. An ATP-gated ion channel, known as P2X7 nucleotide receptor, is expressed in many cell types and significantly influences the mechanosensation mechanisms, and P2X7 receptor is suggested to be essential for PGE2 release in response to mechanical strain [11].

In turn, NO and PGE2 are also expressed in and affect osteoblasts and osteoclasts. In bone, NO suppresses bone resorption and promotes its formation. Indeed, both osteoblasts and osteocytes release NO in response to mechanical strain or FFSS [61], but osteoblasts are less sensitive to FFSS than osteocytes [35, 36]. Under mechanical forces, osteocytes synthetize and release PGE2, which has been known to be one of the earliest responses to loading [62, 63, 64, 65]. FFSS promotes gap junction-mediated intercellular communication and stimulates Cx43 expression, which in turn forms hemichannels, allowing thereby to release PGE2 [62]. The latter acts in an autocrine fashion and activates EP2-EP4 receptors expressed on osteocytes, while acting in a paracrine fashion to modulate osteoblast and osteoclast activities. The activation of EP2-EP4 is associated with increased intracellular cAMP and activated protein kinase A (PKA), which regulates the expression of many downstream effectors, including RANKL, Dmp1, and Sost. In vivo, new bone formation is induced by PGE2 and anabolic loading effects are blocked by indomethacin [66]. PGE2 seem to be released in response to shear stress through hemichannels unopposed halves of gap-junction channels [62]. These hemichannels in osteocytes exert multiple functions, including the protection of cell viability and the release of signaling factors [67, 68].

The canonical Wnt-signaling plays a significant role in bone homeostasis and mainly targets osteocytes among bone cells. Wnt activity increases under mechanical loading and decreases during unloading. Most of these effects are, indeed, mediated by sclerostin [49]. Mechanical transduction is also associated with sex hormones that are estrogens and androgens, and a close relationship has been determined between them and skeletal mechanobiology [69, 70, 71]. Insulin-like growth factor 1 (IGF-1) was also found to form another signaling pathway required for proper mechanotransduction [72].

5. Mathematical models

5.1 Main biological approaches of bone cell dynamics

Komarova et al. [73] developed a first mathematical model for the interaction between osteoblasts and osteoclasts, where the temporal osteoblast and osteoclast population dynamics and the associated changes in bone mass at a single BMU, were constructed. The originality of this work was the incorporation of autocrine and paracrine interactions among osteoblasts and osteoclasts, allowing to investigate the cooperative roles of both of these regulation mechanisms in bone remodeling control. The cell population dynamics were described using the following system of differential equations:

$$\frac{d{x}_1}{\mathit{dt}}={\alpha }_1{x}_1^{g_{11}}{x}_2^{g_{21}}-{\beta }_1{x}_1$$E1

$$\frac{d{x}_2}{\mathit{dt}}={\alpha }_2{x}_1^{g_{12}}{x}_2^{g_{22}}-{\beta }_2{x}_2$$E2

where $1$ and $2$ denote the osteoclasts and osteoblasts, respectively, $x_i$ the cell number, ${\alpha }_i$ the cell production activity, ${\beta }_i$ the cell removal activity, and $g_{\mathit{ij}}$ the net effectiveness of osteoclast- or osteoblast-derived autocrine or paracrine factors.

The changes in bone mass was determined using the following equation:

$$\frac{\mathit{dz}}{\mathit{dt}}=-k_1{y}_1+k_2{y}_2$$E3

where $z$ denotes the total bone mass, $k_i$ the normalized resorption and formation activities, and $y_i$ the active cell numbers, which is calculated according to the following conditions, with ${\overline{x}}_i$ being the cell number at steady state.

The findings revealed that the remodeling dynamic behavior mode mainly depends on osteoclast autocrine regulation parameter and the model suggests that preosteoblast availability may be a limiting factor in bone formation under certain conditions. This study revealed that modeling the simultaneous processes of osteoblast and osteoclast regulations and interactions, even in a simplistic form, results in a highly complex nonlinear behavior, and that the intrinsic properties of the osteoblast–osteoclast system can generate complex remodeling modes observed in vivo. However, only two cell types were taken into account, local autocrine and paracrine factors were supposed to only regulate osteoblast and osteoclast formation, and the parameters describing the autocrine and paracrine regulation effectiveness included actions of several factors.

Then, Lemaire et al. [74] developed a theoretical framework able to explain the experimental observations of bone biology. In their paper, a mathematical model of bone remodeling cellular control was proposed to particularly examine the biochemical control network failures leading to bone diseases, such as osteoporosis. The model consists of a synthetic system including the cellular and biomechanical feedback mechanisms that spearhead bone turnover regulation, taking into account the PTH action in the remodeling process. The originality of this model was the incorporation of the RANK-RANKL-OPG pathway, representing an essential regulation mechanism of osteoclast formation.

The reaction scheme of PTH binding with its receptor was formulated as follows, without taking osteoblastic interactions into account:

$$\begin{array}{c}\begin{array}{c}{p}_P\\ \downarrow \end{array}\\ P\\ \begin{array}{c}\downarrow \\ d_P\end{array}\end{array}+P_r\begin{array}{c}{k}_5\\ \rightleftharpoons \\ k_6\end{array}{P}_r\diamond P$$E4

where $P$ denotes PTH, $P_r$ the PTH receptor, $p_P$ and $d_P$ are the PTH production and dissociation fluxes, respectively, and $P_r\diamond P$ is the complex formed by PTH and its receptor.

The reaction schemes of the bindings of RANKL with RANK and of OPG with RANKL were giving by (Eq. 9) and (Eq. 10), respectively:

$$\begin{array}{c}\begin{array}{c}{p}_L\\ \downarrow \end{array}\\ L\\ \begin{array}{c}\downarrow \\ d_L\end{array}\end{array}+K\begin{array}{c}{k}_3\\ \rightleftharpoons \\ k_4\end{array}K\diamond L$$E5

$$\begin{array}{c}\begin{array}{c}{p}_O\\ \downarrow \end{array}\\ O\\ \begin{array}{c}\downarrow \\ d_O\end{array}\end{array}+\begin{array}{c}\begin{array}{c}{p}_L\\ \downarrow \end{array}\\ L\\ \begin{array}{c}\downarrow \\ d_L\end{array}\end{array}\begin{array}{c}{k}_1\\ \rightleftharpoons \\ k_2\end{array}O\diamond L$$E6

where $K$ denotes the RANK receptor, $L$, the RANKL cytokine, $O$ the OPG protein, $K\diamond L$ the RANK-RANKL complex, and $O\diamond L$ the OPG-RANKL complex.

The model was able to simulate the coupling mechanism between osteoblasts and osteoclasts, the catabolic effect related to PTH continuous administration, the RANKL catabolic action and the OPG anti-catabolic action, in addition to metabolic bone diseases, such as glucocorticoid excess, senescence, vitamin D deficiency, and estrogen deficiency. The model also confirmed that bone formation therapies yielded better results that anti-resorptive therapies in restoring bone loss, and that combining anabolic and anti-resorptive therapies may provide better benefits than monotherapy.

Later, Pivonka et al. [75] developed a model to investigate and incorporate an optimal model structure for RANKL and OPG expression on osteoblast lineage at different maturation stages. Afterwards, the investigation dealt with optimal changes in differentiation rates able to provide effective functional control within an active BMU. The cell population model proposed in this study was mainly based on that of Lemaire et al. [74], but incorporating a rate equation describing changes in bone volume, a rate equation describing TGF-β concentration in terms of the resorbed bone volume, RANKL and OPG expressions on osteoblast-lineage cells at different maturation stages, as well as activator/repressor functions based on enzyme kinetics. The model does not refer to a single BMU. It includes spatial averages of cell numbers over a finite bone volume that contains many BMUs. But, it may be contrasted in the case of studying a single BMU, since temporal and spatial sequences define the type of the present bone cells.

The cell population dynamics were described by the following cell balance equations:

$$\frac{d{\mathit{OB}}_p}{\mathit{dt}}=D_{{\mathit{OB}}_u}·{\pi }_{\mathit{act},{\mathit{OB}}_u}^{\mathit{TGF}-\beta }-D_{{\mathit{OB}}_p}·{\mathit{OB}}_p·{\pi }_{\mathit{act},{\mathit{OB}}_p}^{\mathit{TGF}-\beta }$$E7

$$\frac{d{\mathit{OB}}_a}{\mathit{dt}}=D_{{\mathit{OB}}_p}·{\mathit{OB}}_p·{\pi }_{\mathit{rep},{\mathit{OB}}_p}^{\mathit{TGF}-\beta }-A_{{\mathit{OB}}_a}·{\mathit{OB}}_a$$E8

$$\frac{d{\mathit{OC}}_a}{\mathit{dt}}=D_{{\mathit{OC}}_p}·{\mathit{OC}}_p·{\pi }_{\mathit{act},{\mathit{OC}}_p}^{\mathit{\text{RANKL}}}-A_{{\mathit{OC}}_a}·{\mathit{OC}}_a·{\pi }_{\mathit{act},{\mathit{OC}}_p}^{\mathit{TGF}-\beta }$$E9

where ${\mathit{OB}}_u$ denotes the uncommitted osteoblast progenitors, ${\mathit{OB}}_p$ the preosteoblast cells, ${\mathit{OC}}_p$ the preosteoclast cells, ${\mathit{OB}}_a$, the active osteoblasts, ${\mathit{OC}}_a$ the active osteoblasts, $D_i$ the cell differentiation rate, $A_i$ the cell apoptosis rate, ${\pi }_{\mathit{act},{\mathit{OB}}_u}^{\mathit{TGF}-\beta }$, ${\pi }_{\mathit{rep},{\mathit{OB}}_p}^{\mathit{TGF}-\beta }$, and ${\pi }_{\mathit{act},{\mathit{OC}}_p}^{\mathit{TGF}-\beta }$ the activator/repressor functions related to the binding of TGF-β to its receptors on osteoblasts and osteoclasts, and ${\pi }_{\mathit{act},{\mathit{OC}}_p}^{\mathit{\text{RANKL}}}$ the activator function related to the binding of RANKL to its RANK on preosteoclasts. The above cell balance equations represent the changes in each cell population owing to the addition and removal of the respective cell lineage. Several activator and repressor function regulate the differentiation and apoptosis rates. For instance, the binding of TGF-β on its receptors expressed on uncommitted osteoblast progenitors promotes their differentiation, whereas its binding on its receptors expressed on preosteoblasts inhibits their differentiation.

The evolution, over time, of bone volume, $\mathit{BV}$ was later formulated as follows, assuming that bone resorption and formation rates are proportional to the active cell numbers:

$$\frac{\mathit{dBV}}{\mathit{dt}}=-k_{\mathit{res}}{\stackrel{\sim }{\mathit{OC}}}_a+k_{\mathit{\text{form}}}{\stackrel{\sim }{\mathit{OB}}}_a$$E10

where $\mathit{BV}$ here denotes the percentage of normalized bone volume, $k_{\mathit{res}}$ the relative bone resorption rate, and $k_{\mathit{\text{form}}}$ the relative bone formation rate, with ${\stackrel{\sim }{\mathit{OC}}}_a={\mathit{OC}}_a\left(t\right)-{\mathit{OC}}_a\left(t_0\right)$ and ${\stackrel{\sim }{\mathit{OB}}}_a={\mathit{OB}}_a\left(t\right)-{\mathit{OB}}_a\left(t_0\right)$, where ${\mathit{OC}}_a\left(t_0\right)$ and ${O\mathit{C}}_a\left(t_0\right)$ denote the numbers of active osteoclasts and osteoblasts at the initial state, $t_0$. This formulation allows to link the evolution of cell numbers to the changes in bone volume.

The outcomes of this study suggested that RANKL expression profile provides BMUs with a best functional responsiveness, and that TGF-β is included in the up-regulation of osteoblast progenitor differentiation rate, in the down-regulation of preosteoblast differentiation rate, and in the up-regulation of active osteoclast apoptosis rate, which partially explains the particular suitability of TGF-β physiological actions in bone.

5.2 Main mechanical models of bone remodeling

Miller et al. [76] used an orthotropic material model to provide a 2D representation of the effective properties of the trabecular bone, with the aim of explaining its structure in the proximal femur. The proposed model explained the directionality of the trabecular bone and provided a quite well prediction of the directional material properties, which supported the consideration of anisotropy in adaptation algorithms. However, the study was based on a 2D plane stress assumption, which cannot reflect the 3D reality, the investigated problem was simplified, and the number of elastic constants per element was reduced.

Bonfoh et al. [77] proposed a bone remodeling model based on an external mechanical stimulus and tested its predictions for simple loading scenarios. Their approach allowed to describe osteoblast and osteoclast interactions when bone is subjected to external loading. The latter was expressed in terms of strain energy density $\omega$ detected by an osteocyte $i$ at its location ${\overrightarrow{x}}^{\left(i\right)}$:

$$\omega \left({\overrightarrow{x}}^{\left(i\right)}\right)=\frac{1}{2}\underset{¯}{\underset{¯}{\sigma }}\left({\overrightarrow{x}}^{\left(i\right)}\right):\underset{¯}{\underset{¯}{\epsilon }}\left({\overrightarrow{x}}^{\left(i\right)}\right)$$E11

where $\underset{¯}{\underset{¯}{\sigma }}$ and $\underset{¯}{\underset{¯}{\epsilon }}$ are the tensors of stress and strain, respectively.

The model predictions provided a coherent remodeling description. Applied to the dental implant osseointegration simulation, the model forecasted plausible results. However, the comparison of the obtained results with previous ones from literature showed significant differences regarding the remodeling area and the stress/strain fields.

Geraldes et al. [78] proposed an orthotropic strain-driven adaptation algorithm to assess the distribution of the volumetric material properties at the femur and the directionality of its internal structures within a continuum. The proposed algorithm included multiple load cases (Figure 2) and the maximum strain components across all frames from the daily physical activities were selected to generate a strain field envelope involving the maximum driving stimuli for the material properties and orientations.

Figure 2.

Key steps in updating the orthotropic material properties and directionality for the multiple load case adaptation process [78].

The study particularly highlighted the importance of stair climbing in the evolution of the properties of the fracture-prone femoral neck region, with implications for prevention strategies of non-pharmacological fracture based on exercise. However, the model was based on geometrical definition of certain muscles through straight lines, which may result in non-physiological lines of action and moment arms. Besides, the used femur geometry is not personalized for the studied case, but extracted from the muscle standardized femur, which may contribute to overassessment of the hip contact forces. The proposed multiple loading scenario did not take into account any impact activity and the performed analysis did not consider time- or frequency-dependent adaptation.

6. Finite element simulations

The primary results (Figure 3) of the study of Miller et al. [76] are the two local Young’s moduli $E_1$ and $E_2$, as well as the predicted density distribution. The model converged after 25 iterations. The results show the predicted properties of the trabecular elements, with the cortical shell shown as an outline.

The results generated by the model of Bonfoh et al. [77] showed that the magnitude of the signal received by pre-osteoblasts and pre-osteoclasts depends on the load intensity the osteocyte concentration that varies according to the age, sex, type of the considered bone, etc. Bone density is also related to the type of the considered bone and mainly to its initial apparent density. For the considered load case, the tops of the threads of the simulated implant are the most loaded areas where the stimulus is sufficient to lead the cell activity into an opposition. However, the strain energy density for the other zones remains deficient. Therefore, the concerned areas are in resorption or in a steady state.

Figure 3.

(Top) Predicted material properties after model converged. The arrow directions show the material axes of each element, which are the predicted trabecular directions. The lengths of the arrows represent the effective elastic moduli E1 and E2 of each element. The legend arrow represents a Young’s modulus of 2.5 GPa. (Bottom) Predicted density distribution. Density was estimated from the average elastic modulus of each element. All recognized major features of the proximal femur are visible. Maximal density is 0.88 g/cm3, which corresponds to a volume fraction of approximately 48%. The division of the model into constant cortical elements and variable trabecular elements prevents unrealistic densities in regions that are close to the cortices. [76].

Among the obtained results of Geraldes et al. [78], Figure 4 shows the predicted density (right) and dominant material orientations (middle) for coronal section of the proximal femur, compared to μCT slices of the same regions (left).

Figure 4.

Predicted density (right, in g/cm³) and dominant material orientations (middle) for a coronal (top) and transverse (bottom) section of the converged proximal femur undergoing multiple load cases. Legends highlighting the most interesting features identified by Singh et al. [79] (top, left) and Tobin [80] (bottom, left) were superimposed onto a μCT slice of the same region. The material orientations associated with E₁ are shown in red and E₃ in blue [78].

7. Conclusion

The BMU, previously seen as a changing structure where homogenous cell populations enter, act, and leave, is now revealed to be a tiny microenvironment where heterogeneous cell types mingle and interact with each other to provide a delicate preservation of bone quality and quantity, with highly coordinated activities between the different actors through a complex communication network. Besides, the integrity of bone homeostasis is significantly based on a proper detection and transduction of the mechanical stimulus to biochemical signals. This is mainly performed by osteocytes that were found to be the major mechanosensory and most abundant cells in an adult skeleton. However, the precise mechanisms by which osteocytes perceive and transduce mechanical forces are still unclear. What have emerged is the complexity and multiplicity of the signaling systems activated by subjecting bone to mechanical strains. Owing to the unique in vivo environment of an osteocyte, it is difficult to establish in vitro models that faithfully replicate the required processes. Still, recent technological advances have demonstrated an impressive progress in understanding osteocyte biology and functions, and further elucidation on the mechanobiological mechanisms of osteocytes holds promises of biological and medical implications.

Mathematical models have recently been widely largely used in biomechanical engineering to better understand bone dynamics and mechanical behavior, and to prevent thereby bone pathologies and fractures. Mathematically describing bone cellular and biomechanical behavior is not an easy task to perform, since bone is a living heterogeneous composite material and attempting to model it strongly depends on numerous mechanical and biochemical factors. It should be noted that accurate predictions are tightly linked to realistic characterization of the material behavior. Many of these models are currently combined with FE method to numerically investigate specific scenarios. Although the FE method presents several limitations, the accuracy of numerical results may be improved due to emerging advancing techniques and FE analysis still provides an important tool for assessing bone properties and understanding its behavior. It also offers a personalized analysis of bone stiffness based on specific measurements. The development of the FE method could successfully prevent and treat age- and disease-related bone fractures. The accurate assessment of bone structure and mechanical properties, with a more realistic prediction of applied stresses, should significantly improve the outcome of FE modeling and make this method effective in treating patients and developing new implant designs.