*Article* **A General Mechano-Pharmaco-Biological Model for Bone Remodeling Including Cortisol Variation**

**Rabeb Ben Kahla 1,2, Abdelwahed Barkaoui 3, Moez Chafra <sup>1</sup> and João Manuel R. S. Tavares 4,\***


**Abstract:** The process of bone remodeling requires a strict coordination of bone resorption and formation in time and space in order to maintain consistent bone quality and quantity. Bone-resorbing osteoclasts and bone-forming osteoblasts are the two major players in the remodeling process. Their coordination is achieved by generating the appropriate number of osteoblasts since osteoblasticlineage cells govern the bone mass variation and regulate a corresponding number of osteoclasts. Furthermore, diverse hormones, cytokines and growth factors that strongly link osteoblasts to osteoclasts coordinated these two cell populations. The understanding of this complex remodeling process and predicting its evolution is crucial to manage bone strength under physiologic and pathologic conditions. Several mathematical models have been suggested to clarify this remodeling process, from the earliest purely phenomenological to the latest biomechanical and mechanobiological models. In this current article, a general mathematical model is proposed to fill the gaps identified in former bone remodeling models. The proposed model is the result of combining existing bone remodeling models to present an updated model, which also incorporates several important parameters affecting bone remodeling under various physiologic and pathologic conditions. Furthermore, the proposed model can be extended to include additional parameters in the future. These parameters are divided into four groups according to their origin, whether endogenous or exogenous, and the cell population they affect, whether osteoclasts or osteoblasts. The model also enables easy coupling of biological models to pharmacological and/or mechanical models in the future.

**Keywords:** biomechanics; mathematical model; cell dynamics; bone physiology; bone disorders

#### **1. Introduction**

Fragility fracture rates are growing exponentially, mainly due to population aging. The World Health Organization has recorded a substantial increase in population growth and aging, with a life expectancy rising from about 65 years old in 2005 up to 73 years old in 2019; while in Africa this latter age is around 64 years old, and it is around 78 years old in Europe and the Western Pacific. Recently, various countries worldwide have experienced a fragility fracture crisis with this increase in life expectancy. In 2017, the International Osteoporosis Foundation reported an estimated 2.7 million fragility fractures in six European countries, mainly Germany, Italy, France, Spain, Sweden and the United Kingdom, which resulted in an associated annual cost of 37.5 billion euros for their health care systems. On an individual level, these fragility fractures affected the independence and quality of life of thousands of people in each of these countries. By 2030, the number of annual fragility

 

**Citation:** Kahla, R.B.; Barkaoui, A.; Chafra, M.; Tavares, J.M.R.S. A General Mechano-Pharmaco-Biological Model for Bone Remodeling Including Cortisol Variation. *Mathematics* **2021**, *9*, 1401. https://doi.org/10.3390/ math9121401

Academic Editor: Mauro Malvè

Received: 9 April 2021 Accepted: 25 May 2021 Published: 17 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

fractures is expected to have increased by 23%, reaching 3.3 million, with projected costs of approximately 47.4 billion euros [1]. This makes bone fractures, mainly fragility fractures, a major public health problem.

The evolution of bone mineral density is related to bone metabolism and the different factors affecting it. In fact, throughout life, bone is constantly being renewed through the resorption of old and damaged bone and the formation of new bone. This renewal is what is known as bone remodeling, a process requiring a strict coordination in time and space to maintain bone quantity and quality. This coordination mainly involves boneresorbing osteoclasts and bone-forming osteoblasts, which are the two major players in the remodeling process. The delicate balance between the amount of resorbed bone and the subsequent deposited new bone requires a close coordination of the resorption and formation activities. This coordination controls the generation of the appropriate number of osteoblasts for remodeling, which is referred to as the coupling mechanism. Moreover, several hormones, cytokines and growth factors are involved in linking the osteoblast- and osteoclast-lineage cells (via a complex network) throughout the remodeling cycle.

The remodeling phenomenon has been a research and discussion subject for over a century. In 1892, Julius Wolff stated that bones will adapt to the degree of mechanical loading. This statement is currently known as Wolff's Law and represents the first explicit statement that directly links bone microstructure to mechanical loading. This law establishes that the tendencies of bone trabeculae are aligned with the principal directions of stress. Since then, the mathematical description of bone behavior has experienced a remarkable development. At first, most models were phenomenological or purely mechanical, lacking any biological foundation. Afterwards, models started gaining insight into the biological processes of bone remodeling and the relationship between them and the mechanical processes. Concepts such as bone multicellular unit activity, metabolic activity, mineralization processes as well as damage accumulation have now been introduced and are included in the new mathematical models. These new bone remodeling models have been classified into mechanological, biological and mechanobiological, according to the main mechanism of the process.

The mechanological models show the mechanical stimulus and the final bone response (net resorption or formation) in a black-box manner [2–8]. These models are mainly based on the mechanostat theory formulated by Frost [9], which represents a sophisticated version of the statement of Wolff's law [10]. Additionally, they only focus on bone biomechanical properties, without taking the biochemical aspects and the cell dynamics into account. This does not provide a clear understanding of these internal processes, and there is no inclusion of bone tissue responses to drugs or hormonal dysfunctions. Later models were then developed as semi-mechanistic versions of this bone remodeling theory to provide more realistic outcomes [11–13].

Two research groups were mainly responsible for initiating the biological models. The first research group was Kamarova et al. [14], who combined the autocrine and paracrine factors into a single exponential parameter. Thus, bone remodeling could be studied using power-law approximations, which was an almost perfect fit, since osteoblast and osteoclast activities could be temporally describe throughout a single remodeling cycle. This Kamarova et al. model [14] has been further developed by several other authors, such as Ryser et al. [15], who predicted the spatiotemporal evolution of a single Bone Multicellular Unit (BMU), and Ayati et al. [16] who reformulated the approach to develop a diffusion bone remodeling model, which enabled the study of myeloma bone disease. The second research group was Lemaire et al. [17], who developed a model including the RANK-RANKL-OPG signaling pathway, combined with the influence of the transforming growth factor-β (TGFβ) and the parathyroid hormone (PTH) on bone cell activities, while taking into account the differentiation stages of the osteoblasts and osteoclasts. Owing to its success, this second approach has been developed further by other groups of authors, such as Maldonado et al. [18], who added the osteocyte function as a mechanotransducer, and Pivonka et al. [19], who used a set of activator/repressor functions for a more appropriate description of the receptor–ligand interactions that occur throughout the remodeling process. The first family of models is simplistic, neglecting the cell differentiation stages that has an important role in the osteoblast–osteoclast interplay, but has the advantage of requiring a relatively low number of parameters, facilitating the numerical and computational calculations. On the other hand, the second family of models reduces the number of simplifications, but requires a significantly high number of parameters, which may affect the results [20].

The mechanobiological models follow the models of either Lemaire et al. [17] or Komarova et al. [14]. Owing to the complexity of simultaneously modeling mechanical and biological components, progress in the implementation of this third category of mathematical models has been successful only in recent years. In this third category of mathematical bone remodeling models, the modeling process starts by calculating the mechanical parameters such as stress, strain, pressure and fluid velocity, which then regulate the biological process that includes cell activities and tissue formation [21–23]. However, a primary limitation of mechanobiological models is the need to solve a series of Partial Differential Equations (PDEs) and to develop a Finite Element (FE) formulation to implement within an iterative procedure [24–26].

Few models tackle bone remodeling as a mechano-chemo-biological process, going from the mechanical stimulus applied to bone up to the generated chemical reactions and followed by bone cell responses [27–33]. All these works were developed using a number of simple assumptions to model, the so-called mechanotransduction mechanism/process. In addition, coupling mechanical and chemical phenomenon together with mechanosensing [34], which is a lesser-known component of the remodeling process [35–37], requires various simplified assumptions.

Following these steps, with the aim to create a novel mechano-pharmaco-biological model, the current article provides a first step towards this goal. This work presents a pharmaco-biological approach that is coupled with a previously developed mechanical approach [38] and the hormonal effects of another previous study [39,40]. Both of these previous works were focused on modelling the effects of cyclic loading and sex hormones on the remodeling process throughout the lifecycle of a person. However, the current work focuses on the effects of sex hormones combined with that of cortisol, TGFβ and PTH on bone remodeling in the case of healthy persons, and in the case of persons with an endogenous hypercortisolism known as Cushing disease. Then, the effects of bisphosphonates and denosumab on the enhancement of the pathologic remodeling process were compared. Since osteoblastic-lineage cells have been found to govern the bone mass variation, the mechanical stimulus was included in the proposed model as an exogenous paracrine model acting on osteoblast concentrations, which, in turn, act on the osteoclast concentrations.

#### **2. Materials and Methods**

#### *2.1. Development of the Bone Remodeling Mathematical Model*

As aforementioned, the proposed approach is based on combining both of the biological models of Komarova et al. [14] and Pivonka et al. [19]. Figure 1 summarizes the cell dynamics according to the model of Komarova et al. [14], and Figure 2 summarizes that of Pivonka et al. [19].

The proposed approach retains the structure of the Komarova et al. model [14] and includes the model parameters according to the Pivonka et al. approach [19]. Then, the effects of sex hormones (estradiol in women and testosterone in men), cortisol and antiresorptive drugs (bisphosphonates and denosumab) are taken into consideration. The evolution of each of these parameters is described according to the concept of Hill functions.

**Figure 1.** Diagram of the cell dynamics according to the model of Komarova et al. [14]. Here, g11 reflects the action of TGFβ and g22 the action of IGF, both of them representing autocrine factors secreted by one cell to influence the other cells of its own lineage. g12 reflects the actions of TGFβ and IGF, acting as paracrine factors secreted by osteoclasts to modulate osteoblasts. g21 reflects the actions of RANKL and OPG, acting as paracrine factors secreted by osteoblasts to modulate osteoclasts.

**Figure 2.** Diagram of the cell dynamics according to the model of Pivonka et al. [19], which considers that osteoblasts modulate osteoclast differentiation through g21. In this diagram, g12 only reflects the action of TGFβ stored in the bone matrix and released by osteoclasts during the bone resorption phase, and g21 reflects the actions of RANKL and OPG expressed by osteoblasts to modulate osteoclast differentiation. The model also involves the effect of PTH, as an external factor to the BMU.

#### 2.1.1. Development of the Pharmaco-Biological Model

Several hormones, cytokines and growth factors influence cell behavior within an active BMU. The development of the proposed biological model was based on including each factor according to the concept of Hill functions [41].

#### Transforming growth factor beta (TGFβ)

TGFβ plays a critical role in bone remodeling. It stimulates the synthesis of matrix protein, dramatically affects osteoblasts and osteoclasts, and is abundant in bone. During bone formation, osteoblasts accumulate TGFβ in a latent form in the bone matrix. During bone resorption, osteoclasts release and activate TGFβ that, in turn, induces the activation of preosteoblasts recruited for the following formation phase and their migration to prior resorptive sites [42–44]. Thus, the preosteoblast concentration increases due to the differentiation of osteoblast progenitors into preosteoblasts, promoted by TGFβ [45,46], and decreases due to the differentiation of preosteoblasts into active osteoblasts, suppressed by TGFβ [45–49]. Additionally, the active osteoclast concentration decreases owing to active osteoclast apoptosis (Table 1). Mathematically, these effects can be expressed through functions: *Ta*(*OC*) expressing the stimulation of osteoclast apoptosis by TGFβ; *Ta*(*OB*) expressing the stimulation of osteoblast progenitor differentiation to preosteoblasts by

TGFβ; and *Tr*(*OB*) conveying the repression of preosteoblast differentiation into active osteoblasts by *TGFβ*:

$$\text{Ta(OC)} = \frac{TGF\beta}{K\_{\text{Ta1}} + TGF\beta} \tag{1}$$

$$\text{Tra}(OB) = \frac{TGF\beta}{K\_{\text{Ta2}} + TGF\beta} \tag{2}$$

$$Tr(OB) = \frac{K\_{Tr2}}{K\_{Tr2} + TGF\beta} \tag{3}$$

where *TGFβ* denotes *TGFβ* concentration, *KTa*<sup>1</sup> the activation coefficient of active osteoclast apoptosis mediated by *TGFβ*, *KTa*<sup>2</sup> the activation coefficient of osteoblast progenitor differentiation mediated by *TGFβ* and *KTr*<sup>2</sup> the repression coefficient of preosteoblast differentiation mediated by *TGFβ*.

**Table 1.** TGFβ actions and their mathematical formulation.

Parathyroid hormone (PTH)

PTH is one of the main endocrine regulators of extracellular phosphate and calcium levels. It up-regulates the RANKL expression by osteoblasts and down-regulates the OPG expression by the same cells, which leads to osteoclast maturation and activity. Moreover, PTH stimulates the production of osteocytic soluble RANK (sRANK) that also promotes osteoclast activity. Meanwhile, PTH suppresses the production of osteocytic sclerostin (SOST) that suppresses osteoblast activity. Although PTH is likely to promote bone resorption and formation, its precise mechanism remains unclear, because of the limited in vivo data [50,51].

In the current model, two actions of PTH are considered: its action on stimulating RANKL production (Equation (4)) and its action on inhibiting OPG production (Equation (5)):

$$\Pr(OC) = \frac{PTH}{K\_{P41} + PTH} \tag{4}$$

$$Pr(OC) = \frac{K\_{Pr1}}{K\_{Pr1} + PHH} \tag{5}$$

where *PTH* denotes the *PTH* concentration, *KPa*<sup>1</sup> the activation coefficient of *PTH* action on RANKL production rate, increasing the preosteoclast differentiation rate, and *KPr*<sup>1</sup> the repression coefficient of *PTH* action on OPG production rate, which also increases the preosteoclast differentiation rate, Table 2.

**Table 2.** PTH actions and their mathematical formulation.

#### Estradiol (Es)

Although the mechanism of estrogen transcriptional activity is not fully understood, it has been suggested that estrogen regulates bone matrix formation by acting on key factors involved in osteoblast differentiation and maturation [52]. The main effect of estrogens is the suppression of bone turnover, probably via osteocytes. Yet, they also inhibit bone resorption, through direct effects on osteoclasts, although the estrogen effects on osteoblasts/osteocytes are also likely to take part. A gap between the resorption and the formation activities has been linked to estrogen deficiency, probably because of the loss of estrogen effects on reducing the osteoblast apoptosis rate, decreasing NF-κB osteoblastic activity, lowering oxidative stress and perhaps other as yet undefined effects [53].

Estrogens may be classified into three main classes: estriol, estrone and estradiol. The latter represents the most potent human endogenous estrogen, with a high affinity for estrogen receptors. Thus, the proposed model is based on the effect of estradiol concentration (Table 3) on cells and bone remodeling through the following two repressive functions:

$$\text{Es(OC)} = \frac{K\_{\text{Es1}}}{K\_{\text{Es1}} + \text{Es(t)}} \tag{6}$$

$$\text{Es(OB)} = \frac{K\_{\text{Es2}}}{K\_{\text{Es2}} + Es(t)}\tag{7}$$

where *Es*(*t*) denotes the estradiol concentration function, *KEs*<sup>1</sup> the repression coefficient of estradiol action on osteoclast differentiation and *KEs*<sup>2</sup> the repression coefficient of estradiol action on osteoblast apoptosis.

**Table 3.** Estradiol actions and their mathematical formulation.

#### Testosterone (Ts)

Androgen receptors have been identified in cultured human fetal osteoblasts using a nuclear-binding assay. Subsequently, mRNA and the androgen receptor protein have been identified in osteoblasts. Almost all studies have shown that androgens enhance the expression of androgen receptors in osteoblasts. Testosterone and 5α-dihydrotestosterone have been found to stimulate the proliferation of cultured preosteoblasts in distinctive species, and the collected evidence generally suggests that androgens stimulate osteoblast differentiation [54].

Moreover, androgen deficiency is most likely to be associated with osteoclast proliferation after orchiectomy. Androgens exert their bone defensive effects indirectly via osteoblasts, and orchidectomy generates preosteoblast proliferation, which increases RANKL secretion, thereby stimulating osteoclast proliferation and activation and resulting in bone loss. In vitro studies have shown that dihydrotestosterone interacts with androgen receptors on osteoclasts and inhibits bone resorption in human osteoclasts [54].

Since testosterone is the main androgen produced by Leydig cells and represents the most well-known male sex hormone, the effect of its concentrations (Table 4) on bone cells and remodeling is included in the model, through the repression function and the activation function, as follows:

$$\text{Ts(OC)} = \frac{K\_{\text{Ts1}}}{K\_{\text{Ts1}} + \text{Ts(t)}} \tag{8}$$

$$\text{Ts(OB)} = \frac{\text{Ts(t)}}{K\_{\text{Ts2}} + \text{Ts(t)}} \tag{9}$$

where *Ts*(*t*) denotes the testosterone concentration function, *KTs*<sup>1</sup> the repression coefficient of testosterone action on osteoclast activity and *KTs*<sup>2</sup> the activation coefficient of testosterone action on preosteoblasts.

**Table 4.** Testosterone actions and their mathematical formulation.

#### Cortisol (Co)

Cortisol has well-established implications for diverse body systems. Specifically, it exerts direct negative effects on bone mineral density (BMD) by affecting bone cell growth, disrupting the bone remodeling process, impairing calcium intestinal absorption and renal reabsorption, as well as by inhibiting activity of sex steroids [55]. Indeed, an imbalance in the serum calcium homeostasis increases osteoclast-resorptive activity and eventually reduces BMD. Even a short bout of elevated cortisol levels may result in a reduced bone formation rate and a lower BMD. Prolonged cortisol oversecretion is consistently associated with a high prevalence of osteoporosis and may be linked to a decrease in BMD with age and to an increase in fragility fracture risk in elderly people [56]. The negative association between cortisol and BMD and the positive association between cortisol and fracture risk have been reported in several studies conducted on healthy older adults [57,58]. Moreover, signs of glucocorticoid-induced osteoporosis (GIOP), including a BMD reduction in the spine, altered ultrasound bone characteristics, as well as a higher number of morphometric

fractures than healthy individuals, were found in patients with adrenal incidentaloma, and diagnosed as having subclinical hypercortisolism [59].

Indeed, it is unclear to date whether physiological cortisol levels also contribute to bone diseases or not. When performing a four-year analysis of single-point serum cortisol levels, no correlation was found between BMD, bone markers and bone loss. However, when analyzing integrated serum cortisol levels over 24 h, a correlation was found with BMD at the femur and the spine. These findings point out that physiological cortisol concentrations affect bone density. However, analyzing its effects may be difficult owing to diurnal variations of serum cortisol [60].

Therefore, the proposed model considers that the cortisol concentration (Table 5) increases osteoclast activity through the activation function; *Co*(*OC*) inhibits osteoblast formation through the repression function; and *Co*(*OB*) inhibits estradiol actions and testosterone actions according to:

$$\text{Co}(\text{OC}) = \frac{\text{Co}(t)}{\text{K}\_{\text{Co}1} + \text{Co}(t)} \tag{10}$$

$$\text{Co}(OB) = \frac{K\_{\text{Co2}}}{K\_{\text{Co2}} + \text{Co}(t)} \tag{11}$$

$$CoE(OC) = -Es(OC)\tag{12}$$

$$\text{CoE}(OB) = -\text{Es}(OB) \tag{13}$$

$$\text{CoT(OC)} = -\text{Ts(OC)}\tag{14}$$

$$\text{CoT}(OB) = -\text{Ts}(OB) \tag{15}$$

where *Co*(*t*) denotes the cortisol concentration function, *KCo*<sup>1</sup> the activation coefficient of cortisol action on preosteoclasts and *KCo*<sup>2</sup> the repression coefficient of cortisol action on preosteoblasts.

#### Bisphosphonates (Bp)

Bisphosphonates affect osteoclasts, but not their precursors. In fact, bisphosphonates are internalized in osteoclasts, probably by endocytosis, and inhibit the synthesis of a key enzyme in the mevalonate (MVA) pathway. This alters the intracellular proteins, accumulates cytotoxic intermediates, and alters the function of osteoclasts, which presumably increase their apoptosis rate [61]. Thus, bone resorption is inhibited (Table 6). The inhibition of osteoclast activity is given by the following repressive function:

$$Bp(OC) = k\_{Bp} \times Bp \tag{16}$$

where *Bp* denotes the administered bisphosphonate concentration and *kBp* the Hill function parameter for drug regulation.

#### Denosumab (De)

Denosumab acts in the same way as OPG does, which is the natural antagonist receptor for RANKL. Denosumab binds to RANKL, thereby deterring the binding of RANKL to its receptor, RANK, expressed on both osteoclast and preosteoclast membranes. Subsequently, the RANK-RANKL signaling pathway is not activated, which impairs preosteoclast differentiation and osteoclast function, and possibly increases the osteoclast apoptosis rate. All of these effects inhibit bone resorption [61]. Thus, the following repression function describes the negative effect of denosumab on osteoclast differentiation and activity (Table 7):

$$De(OC) = \frac{C\_D^{scr}}{0.6667 \times d\_D \times 3 \times 10^{-16}}\tag{17}$$

where *Cser <sup>D</sup>* denotes the serum concentration of denosumab and *dD* its administered dose.


**Table 5.** Cortisol actions and their mathematical formulation.


**Table 7.** Denosumab actions and their mathematical formulation.

#### 2.1.2. Mechanical Model

The mechanical approach developed by Bonfoh et al. [62] was adopted on a macroscopic scale. The local dimension sites of the BMU, the total number of osteocytes and their sensitivity as well as the interactions between bone cells were not considered when formulating the expression of the mechanical stimulus. However, when assuming an elastic isotropic behavior for bone, the mechanical stimulus can be expressed by:

$$
\Delta\psi(\stackrel{\rightarrow}{\dot{\mathfrak{x}}}) = w \frac{\begin{pmatrix} \stackrel{\rightarrow}{\dot{\mathfrak{x}}}(i) \\ \stackrel{\rightarrow}{\rho} \end{pmatrix}}{}-W^\* \tag{18}
$$

where *W*∗ denotes the balance stimulus [63], *w* <sup>→</sup> *x* (*i*) the strain energy density and *ρ* the apparent density of bone. The mechanical signal *w* detected by an osteocyte *i* at its location → *x i* is described by:

$$w\left(\stackrel{\rightarrow}{x}^{(i)}\right) = \frac{1}{2}\sigma\left(\stackrel{\rightarrow}{x}^{(i)}\right) : \mathfrak{e}\left(\stackrel{\rightarrow}{x}^{(i)}\right) \tag{19}$$

where *σ* and *ε* denote the stress and strain tensors, respectively.

Bone is naturally damaged under the effect of the daily stresses it is subjected to, which leads to its fatigue and thereby to its aging and the occurrence and propagation of microcracks. To describe the evolution of bone mechanical properties over time, a fatigue damage formulation can be used [64–66], and the damage resulting from cyclic loading can be modeled using the life cycle approach suggested by Chaboche [67]. The damage can reach a maximum value of 1 (one), referring to material failure, and can be expressed in terms of the number of failure cycles, *Nf* , given by [68–70]:

$$N\_f = \mathbb{C}\Delta\varepsilon^{-\vartheta} \tag{20}$$

where Δ*ε* denotes the amplitude of the applied microstrains, which refer to the equivalent strain defined by *εeq* = 2 <sup>3</sup> *εijεij*, and *C* and *ϑ* are constants obtained from experimental data in the literature. These two variables, *C* and *ϑ*, depend on the nature of the solicitation [71]:

$$N\_{f, \ c} = 1.479 \cdot 10^{-21} \Delta \varepsilon^{-10.3} \text{ for compressor loads} \tag{21}$$

$$N\_{f, \ t} = 3.630 \cdot 10^{-32} \Delta \varepsilon^{-14.1} \text{ for tensile loads} \tag{22}$$

The damage at failure (*d* = 1) represents the accumulated damage at a given cycle expressed as:

$$d^{n+1} = d^n + \delta d \tag{23}$$

where *δd* denotes the incremental damage to the cycle (*n* + 1).

Nonlinear cumulative damage is generally characterized using the expression *δd* = 1 *Nf β* . However, here the following nonlinear simplistic evolution is used:

$$
\delta d = \frac{1}{N\_f} \tag{24}
$$

Afterwards, cumulative damage is expressed using accumulated stress cycles:

$$d^{u+1} = d^u + \delta d = \frac{N+1}{N\_f} \tag{25}$$

where *N* denotes the number of loading cycles.

When isotropic properties are assigned to the adopted material, the incremental damage can be expressed using the compressive and/or the tensile fatigue cycles. Hence, total damage *<sup>δ</sup>di*, which is induced by an osteocyte *<sup>i</sup>* at its location <sup>→</sup> *x* (*i*) , refers to the sum of both compressive *δdi*, *<sup>c</sup>* and tensile *δdi*, *<sup>t</sup>* damages. The latter depends on the microstrain amplitude Δ*ε<sup>i</sup>* detected by each osteocyte *i* [69]:

$$
\delta d\_{i,\ c} = \frac{1}{N\_{f,\ c}} = \frac{1}{1.479 \cdot 10^{-21} \Delta \varepsilon^{-10.5}}\tag{26}
$$

$$
\delta d\_{i, \, t} = \frac{1}{N\_{f, \, t}} = \frac{1}{3.630 \cdot 10^{-32} \Delta \varepsilon^{-14.1}} \tag{27}
$$

$$
\delta d\_{\bar{i}} = \delta d\_{\bar{i}\_{\prime} \cdot \mathcal{c}} + \delta d\_{\bar{i}\_{\prime} \cdot \mathcal{t}} \tag{28}
$$

Using the example given by Baste et al. [72], an isotropic formulation for the damage affecting bone properties was developed by multiplying blank modules by (1 − *di*). Therefore, the elastic moduli are expressed as:

$$E = E\_i^0 (1 - d\_i)^2 \tag{29}$$

Since the damage is only considered in the longitudinal directions, the values of the shear moduli remain constant.

#### *2.2. Overview of the Whole Model*

In order to visualize the effects of all the included parameters more clearly, Figure 3 summarizes the developed model, and it depicts the level at which each parameter acts and whether this action is a stimulation or an inhibition.

**Figure 3.** Summary of the proposed bone remodeling model. The shown numbers indicate the involved parameter; the positive and negative signs indicate stimulation or inhibition, respectively, exerted by each factor; the blue color indicates the biological parameters of the model, the green the pharmacological parameters and the purple the mechanical loading.

The pharmaco-biological parameters mentioned above are grouped into autocrine parameters, *A<sup>i</sup>* , and paracrine parameters that are subdivided into endogenous factors, *Pi EN* , produced by the human body, and exogenous ones, *<sup>P</sup><sup>i</sup> EX* , introduced into the human body. Hence, the differential equations of cell dynamics can be written as follows:

$$\Rightarrow \begin{cases} \begin{array}{c} \frac{d\text{OC}}{dt} = a\_{\text{OC}} \text{OC}^{A^{\text{OC}}} \text{OB}^{P\_{\text{FN}}^{\text{OC}}} + P\_{\text{E}\_{\text{N}}}^{\text{OC}} - \beta\_{\text{OC}} \text{OC} \\\ \frac{d\text{OB}}{dt} = a\_{\text{OB}} \text{OB}^{A^{\text{OB}}} \text{OC}^{P\_{\text{E}\_{\text{N}}}^{\text{OB}}} - \beta\_{\text{OB}} \text{OB} \end{array} \tag{30}$$

where *OC* denotes the osteoclast concentration, *OB* the osteoblast concentration, *α*<sup>i</sup> the cell production activities and *β*<sup>i</sup> the cell elimination activities, with:

$$P\_{E\_N}^{\text{OC}} \gets \text{Ta(OC)} + \text{Pr(OC)} + \text{Pa(OC)} + \text{Es(OC)} + \text{Ts(OC)} + \text{Co(OC)} + \text{CoE(OC)} + \text{CoT(OC)}$$

$$P\_{E\_N}^{\text{OB}} \gets \text{Ta(OB)} + \text{Tr(OB)} + \text{Es(OB)} + \text{Ts(OB)} + \text{Co(OB)} + \text{CoE(OB)} + \text{CoT(OB)}$$

$$P\_{E\_N}^{\text{OC}} \gets \text{Bp(OC)} + \text{De}(\text{OC})$$

#### *2.3. Mechanobiological Coupling*

Nowadays, osteoclast activity is acknowledged as being modulated by the osteoblastic cell lineage. By adapting the model of Pastrama et al. [31] to that of Bonfoh et al. [62], the proposed approach considers that the mechanical stimulus, *ψ*, that the bone is subjected to, acts as an exogenous paracrine factor on the variation of osteoblast concentration, according to:

$$\Rightarrow \begin{cases} \begin{array}{c} \frac{dOC}{dt} = \mathfrak{a}\_{\mathsf{OC}} \mathsf{OC}^{A^{\mathsf{OC}}} \mathcal{O} \mathcal{B}^{P^{\mathcal{C}}\_{\mathcal{N}} + P^{\mathcal{O}}\_{\mathcal{E}\_{\mathcal{N}}}}\_{\mathcal{E}\_{\mathcal{N}}} - \mathcal{f}\_{\mathsf{OC}} \mathcal{O} \mathcal{C} \\\ \frac{dOB}{dt} = \mathfrak{a}\_{\mathsf{OB}} \mathcal{O} \mathcal{B}^{A^{\mathcal{O}}} + P^{\mathcal{O}}\_{\mathcal{E}\_{\mathcal{N}}} \mathcal{O} \mathcal{C}^{P^{\mathcal{O}}\_{\mathcal{E}\_{\mathcal{N}}}} - \mathcal{f}\_{\mathsf{OB}} \mathcal{O} \mathcal{B} \end{array} \tag{31}$$

with:

$$P\_{E\_X}^{OB} \leftarrow a + b e^{-\gamma \Delta \psi} \tag{32}$$

where *a* = 1.6, *b* = −0.49 and *γ* = 16.67 *g*/*J*.

As mentioned above, the mechanical stimulus was included as an exogenous paracrine model acting on the osteoblast concentration since osteoblastic-lineage cells govern the bone mass variation. The osteoblast ability to influence osteoclast formation in a paracrine manner has been clearly demonstrated over the years. In fact, osteoblasts modulate osteoclast formation and activity by synthesizing a number of cytokines and growth factors. This is achieved through direct contact between the two cells, established by exchanging small water-soluble molecules through gap junctions [73].

In summary, the mechanical loading that the bone is subjected to activates osteoblasts and increases their concentration according to its intensity. The increase in osteoblast concentration generates higher RANKL and OPG production rates, which modulates the appropriate osteoclast concentration for the subsequent bone resorption activity. Bone resorption releases the matrix-embedded factors, TGFβs, that act on both cell populations to regulate their differentiation, activity and apoptosis. Hormones, including PTH, cortisol, estradiol and testosterone, are also important factors that are external to the BMU, but which regulate the cell dynamics during the remodeling event (Figure 4). The estradiol parameter is considered when the study is focused on the remodeling process in women, and the testosterone parameter is considered when the study is focused on the male remodeling process. Furthermore, when the cortisol concentration is higher than under physiologic conditions, RANKL production by osteoblastic cells increases, inducing an increase in the osteoclast differentiation rate, and thereby leading to a higher resorption rate. This affects the bone response to the mechanical loading by decreasing the bone mass density and altering its mechanical response.

**Figure 4.** Diagram summarizing the proposed mechano-pharmaco-biological model.

#### *2.4. Bone Mass Evolution*

Bone mass evolution throughout time is giving by:

$$\frac{dm}{dt} = (k\_{\rm OB} \times OB) - (k\_{\rm OC} \times OC) \tag{33}$$

where *kOB* and *kOC* denote the normalized formation and resorption activities, respectively, and *OB* and *OC* denote the osteoblast and osteoclast concentrations, respectively.

#### **3. Discussion and Conclusions**

The focus of this article was to provide a pharmaco-biological bone remodeling model that could be easily coupled with mechanical models and was extendable to be able to include various parameters, and consequently allowing the simulation of bone physiologic metabolism for pathologic disorders. At present, the model is designed to simulate the influence of an endogenous hypercortisolism, which is caused by an excessive secretion of cortisol, on bone response to mechanical stimulus and fatigue damage to which it is subjected. The current model combines two of the most current biological models used to predict the evolution of bone mechanical properties.

The most important parameters considered were the autocrine, endogenous paracrine or exogenous paracrine parameters; in short, these are the different ways any parameter may act on bone cell dynamics during the remodeling event. Previous models either neglected various parameters [14,16,74], or the number of parameterswere fixed [17,19,75,76]. However, the current model is easily extendable and can include various other parameters that can provide support in the remodeling process, since these parameters primarily act on osteoblasts or on osteoclasts, whether they are endogenous or exogenous. The proposed approach is able to track the changes in bone remodeling that are specific to each parameter. Consequently, it gives a better overview of the remodeling process by regrouping several parameters at once, instead of simulating one or a limited number of parameters each time. Furthermore, it is able to choose to neglect any unneeded parameter, according to the goal of the study in question.

The current model has some limitations since it assumes an isotropic homogenous material, which idealizes the bone behavior and can affect the outcome. Additionally, the model does not consider any differentiation stages of the bone cells in the mathematical formulations of the cell dynamics. Still, the aim was to provide a mathematical model that applied to any metabolic bone disorder, any drug administered to treat such disorder, and at the same time to track the evolution of hormones and growth factors incorporated in the model, in order to be able to adjust the drug dosage specifically for each patient.

In summary, the current model was developed based on the role of osteoblasts and osteoclasts in renewing bone, as the main players in the remodeling process. All the other factors involved were considered according to their effects on each of the cell lineages involved. Yet, it is well known that a number of these factors, such as the hormonal factors, are provided to osteoblasts and osteoclasts through blood microvessels across bone tissue. This draws the attention to the general problem of the pressure that the microvessel networks exert on bone, and their role in the bone remodeling response to the mechanical stimuli applied to it. Therefore, more accurate results may be found when taking the effect of microvessel networks into account in the mathematical and numerical bone modeling equations.

In forthcoming studies, the results of the current model will be analyzed, and the work will be extended to implement the finite element method, and to visualize the effects of hypercortisolism on a virtual 3D bone model. The proposed model can be confirmed and validated by conducting an experimental study. This will reveal its accuracy in simulating and predicting bone strength under cyclic loading, considering the physiologic conditions and the disorder related to hypercortisolism as described.

**Author Contributions:** All the authors contributed to this work: R.B.K. and A.B. for the mathematical model formulation, search for information about the model, and writing; M.C. and J.M.R.S.T. in the analysis, correction, discussion and writing. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

