**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\mathbf{x}\_1}{dt} = \alpha\_1 \mathbf{x}\_1^{\mathbf{g}\_{11}} \mathbf{x}\_2^{\mathbf{g}\_{21}} - \beta\_1 \mathbf{x}\_1 \tag{1}$$

$$\frac{d\mathbf{x}\_2}{dt} = a\_2 \mathbf{x}\_1^{\mathbf{g}\_{12}} \mathbf{x}\_2^{\mathbf{g}\_{22}} - \beta\_2 \mathbf{x}\_2 \tag{2}$$

where 1 and 2 denote the osteoclasts and osteoblasts, respectively, *xi* the cell number, *α<sup>i</sup>* the cell production activity, *β<sup>i</sup>* the cell removal activity, and *gij* the net effectiveness of osteoclast- or osteoblast-derived autocrine or paracrine factors.

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

$$\frac{dz}{dt} = -k\_1 \mathbf{y}\_1 + k\_2 \mathbf{y}\_2 \tag{3}$$

where *z* denotes the total bone mass, *ki* the normalized resorption and formation activities, and *yi* the active cell numbers, which is calculated according to the following conditions, with *xi* 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}{ll} P\_P \\ \downarrow \\ P & + P\_r \rightleftharpoons P\_r \circ P \\ \downarrow \\ d\_P & \end{array} \tag{4}$$

where *P* denotes PTH, *Pr* the PTH receptor, *pP* and *dP* are the PTH production and dissociation fluxes, respectively, and *Pr* ⋄ *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}{ll} p\_L \\ \downarrow \\ L & +K \rightleftharpoons K \bullet L \\ d\_L \\ d\_L \\ \end{array} \tag{5}$$
 
$$\begin{array}{ll} p\_O & p\_L \\ \downarrow & & \downarrow \\ O & + & L \end{array} \rightleftharpoons O \circ L \\ d\_O & d\_L \end{array} \tag{6}$$

where *K* denotes the RANK receptor, *L*, the RANKL cytokine, *O* the OPG protein, *K* ⋄ *L* the RANK-RANKL complex, and *O* ⋄ *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.

*Cell Interaction and Mechanobiological Modeling of Bone Remodeling Process DOI: http://dx.doi.org/10.5772/intechopen.95045*

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

$$\frac{dOB\_p}{dt} = D\_{OB\_u} \cdot \pi\_{act, OB\_u}^{TGF-\beta} - D\_{OB\_p} \cdot OB\_p \cdot \pi\_{act, OB\_p}^{TGF-\beta} \tag{7}$$

$$\frac{d\text{OB}\_{a}}{dt} = D\_{\text{OB}\_{p}} \cdot \text{OB}\_{p} \cdot \pi\_{\text{rep,OB}\_{p}}^{\text{TGF}-\beta} - A\_{\text{OB}\_{a}} \cdot \text{OB}\_{a} \tag{8}$$

$$\frac{d\text{OC}\_{a}}{dt} = D\_{\text{OC}\_{p}} \cdot \text{OC}\_{p} \cdot \pi\_{\text{act,OC}\_{p}}^{\text{RANKL}} - A\_{\text{OC}\_{a}} \cdot \text{OC}\_{a} \cdot \pi\_{\text{act,OC}\_{p}}^{\text{TGF}-\beta} \tag{9}$$

where *OBu* denotes the uncommitted osteoblast progenitors, *OBp* the preosteoblast cells, *OCp* the preosteoclast cells, *OBa*, the active osteoblasts, *OCa* the active osteoblasts, *Di* the cell differentiation rate, *Ai* the cell apoptosis rate, *πTGF*�*<sup>β</sup> act*,*OBu* , *πTGF*�*<sup>β</sup> rep*,*OBp* , and *πTGF*�*<sup>β</sup> act*,*OCp* the activator/repressor functions related to the binding of TGF-β to its receptors on osteoblasts and osteoclasts, and *πRANKL act*,*OCp* 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, *BV* was later formulated as follows, assuming that bone resorption and formation rates are proportional to the active cell numbers:

$$\frac{dBV}{dt} = -k\_{res}\bar{O}\mathbf{C}\_a + k\_{form}\bar{O}\mathbf{B}\_a \tag{10}$$

where *BV* here denotes the percentage of normalized bone volume, *kres* the relative bone resorption rate, and *kform* the relative bone formation rate, with *OC*<sup>~</sup> *<sup>a</sup>* <sup>¼</sup> *OCa*ðÞ�*<sup>t</sup> OCa*ð Þ *<sup>t</sup>*<sup>0</sup> and *OB*<sup>~</sup> *<sup>a</sup>* <sup>¼</sup> *OBa*ðÞ�*<sup>t</sup> OBa*ð Þ *<sup>t</sup>*<sup>0</sup> , where *OCa*ð Þ *<sup>t</sup>*<sup>0</sup> and *<sup>O</sup>Ca*ð Þ *<sup>t</sup>*<sup>0</sup> 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 upregulation 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 *ω* detected by an osteocyte *i* at its location *x* !ð Þ*<sup>i</sup>*

$$\rho(\overrightarrow{\boldsymbol{\mathfrak{x}}}^{(i)}) = \frac{1}{2}\underline{\boldsymbol{\mathfrak{e}}}(\overrightarrow{\boldsymbol{\mathfrak{x}}}^{(i)}) : \underline{\boldsymbol{\mathfrak{e}}}(\overrightarrow{\boldsymbol{\mathfrak{x}}}^{(i)}) \tag{11}$$

:

where *σ* and *ε* 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

**Figure 2.**

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

*Cell Interaction and Mechanobiological Modeling of Bone Remodeling Process DOI: http://dx.doi.org/10.5772/intechopen.95045*

envelope involving the maximum driving stimuli for the material properties and orientations.

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*<sup>1</sup> 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.

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 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].*

#### *Biomechanics and Functional Tissue Engineering*

**Figure 4.**

*Predicted density (right, in g/cm3 ) 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 E1 are shown in red and E3 in blue [78].*
