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

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 osteoblastic-lineage 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.


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 priate 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][22][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][25][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][28][29][30][31][32][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][36][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.

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.  [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.  [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.
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.

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,43,44]. Thus, the preosteoblast concentration increases due to the  [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.   [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.
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.

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,43,44]. Thus, the preosteoblast concentration increases due to the  [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.

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][43][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][46][47][48][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β: where TGFβ denotes TGFβ concentration, K Ta1 the activation coefficient of active osteoclast apoptosis mediated by TGFβ, K Ta2 the activation coefficient of osteoblast progenitor differentiation mediated by TGFβ and K Tr2 the repression coefficient of preosteoblast differentiation mediated by TGFβ. Table 1. TGFβ actions and their mathematical formulation.
Diagram 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,46,47,48,49]. Additionally, the active osteoclast concentration decreases owing to active osteoclast apoptosis (Table 1). Mathematically, these effects can be expressed through functions: ( ) expressing the stimulation of osteoclast apoptosis by TGFβ; ( ) expressing the stimulation of osteoblast progenitor differentiation to preosteoblasts by TGFβ; and ( ) conveying the repression of preosteoblast differentiation into active osteoblasts by TGFβ: where denotes TGFβ concentration, the activation coefficient of active osteoclast apoptosis mediated by TGFβ, the activation coefficient of osteoblast progenitor differentiation mediated by TGFβ and the repression coefficient of preosteoblast differentiation mediated by TGFβ. 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 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)): where PTH denotes the PTH concentration, K Pa1 the activation coefficient of PTH action on RANKL production rate, increasing the preosteoclast differentiation rate, and K Pr1 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.

Diagram
Mathematics 2021, 9, x FOR PEER REVIEW 6 of 19 where denotes the PTH concentration, the activation coefficient of PTH action on RANKL production rate, increasing the preosteoclast differentiation rate, and the repression coefficient of PTH action on OPG production rate, which also increases the preosteoclast differentiation rate, Table 2. 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: where ( ) denotes the estradiol concentration function, the repression coefficient of estradiol action on osteoclast differentiation and the repression coefficient of estradiol action on osteoblast apoptosis.
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: where Es(t) denotes the estradiol concentration function, K Es1 the repression coefficient of estradiol action on osteoclast differentiation and K Es2 the repression coefficient of estradiol action on osteoblast apoptosis.  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:  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].

Concentration evolution
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: Formulation K Es1 +Es(t) ; Es(OB) = K Es2 K Es2 +Es(t)

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: Ts(OB) = Ts(t) K Ts2 + Ts(t) (9) where Ts(t) denotes the testosterone concentration function, K Ts1 the repression coefficient of testosterone action on osteoclast activity and K Ts2 the activation coefficient of testosterone action on preosteoblasts. Table 4. Testosterone actions and their mathematical formulation.

Diagram
Mathematics 2021, 9, x FOR PEER REVIEW 8 of 19 where ( ) denotes the testosterone concentration function, the repression coefficient of testosterone action on osteoclast activity and the activation coefficient of testosterone action on preosteoblasts. 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,

Concentration evolution
Mathematics 2021, 9, x FOR PEER REVIEW 8 of 19 where ( ) denotes the testosterone concentration function, the repression coefficient of testosterone action on osteoclast activity and the activation coefficient of testosterone action on preosteoblasts. 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,

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 Mathematics 2021, 9, 1401 9 of 18 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].

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: where Bp denotes the administered bisphosphonate concentration and k Bp 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): where C ser D denotes the serum concentration of denosumab and d D its administered dose.

Actions
Stimulation of preosteoclast differentiation into active osteoclasts (Co(OC) ). Inhibition of osteoblast activity (Co(OB) ). Inhibition of estradiol actions on osteoblasts and osteoclasts (CoE ). Inhibition of testosterone actions on osteoblasts and osteoclasts (CoT ).

Diagram
Mathematics 2021, 9, x FOR PEER REVIEW 10 of 19 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: where denotes the administered bisphosphonate concentration and the Hill function parameter for drug regulation.

Concentration evolution
Mathematics 2021, 9, x FOR PEER REVIEW 10 of 19 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: where denotes the administered bisphosphonate concentration and the Hill function parameter for drug regulation.

Diagram
Mathematics 2021, 9, x FOR PEER REVIEW 11 of 19 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

Actions
Inhibition of preosteoclast differentiation into active osteoclasts, and their activity (De(OC) ).

Diagram
Mathematics 2021, 9, x FOR PEER REVIEW 12 of 19 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: where * denotes the balance stimulus [63], ⃗ ( ) the strain energy density and the apparent density of bone. The mechanical signal detected by an osteocyte at its location ⃗ is described by: 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,65,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, , given by [68,69,70]: where Δ denotes the amplitude of the applied microstrains, which refer to the equivalent

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: where W * denotes the balance stimulus [63], w → 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: 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][65][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, N f , given by [68][69][70]: where ∆ε denotes the amplitude of the applied microstrains, which refer to the equivalent strain defined by ε eq = 2 3 ε 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 , t = 3.630·10 −32 ∆ε −14.1 for tensile loads (22) The damage at failure (d = 1) represents the accumulated damage at a given cycle expressed as: where δd denotes the incremental damage to the cycle (n + 1). Nonlinear cumulative damage is generally characterized using the expression However, here the following nonlinear simplistic evolution is used: Afterwards, cumulative damage is expressed using accumulated stress cycles: 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 δd i , which is induced by an osteocyte i at its location , refers to the sum of both compressive δd i, c and tensile δd i, t damages. The latter depends on the microstrain amplitude ∆ε i detected by each osteocyte i [69]: 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 − d i ). Therefore, the elastic moduli are expressed as: Since the damage is only considered in the longitudinal directions, the values of the shear moduli remain constant.

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.

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: with: where = 1.6, = −0.49 and = 16.67 ⁄ . 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 The pharmaco-biological parameters mentioned above are grouped into autocrine parameters, A i , and paracrine parameters that are subdivided into endogenous factors, P i E N , produced by the human body, and exogenous ones, P i E X , introduced into the human body. Hence, the differential equations of cell dynamics can be written as follows: where OC denotes the osteoclast concentration, OB the osteoblast concentration, α i the cell production activities and β i the cell elimination activities, with:

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: 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. 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.

Bone Mass Evolution
Bone mass evolution throughout time is giving by: where and denote the normalized formation and resorption activities, respectively, and and denote the osteoblast and osteoclast concentrations, respectively.

Discussion and Conclusions
The focus of this article was to provide a pharmaco-biological bone remodeling

Bone Mass Evolution
Bone mass evolution throughout time is giving by: where k OB and k OC denote the normalized formation and resorption activities, respectively, and OB and OC denote the osteoblast and osteoclast concentrations, respectively.

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 parameters were 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.

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