Combined Effects of Exercise and Denosumab Treatment on Local Failure in Post-menopausal Osteoporosis–Insights from Bone Remodelling Simulations Accounting for Mineralisation and Damage

Denosumab has been shown to increase bone mineral density (BMD) and reduce the fracture risk in patients with post-menopausal osteoporosis (PMO). Increase in BMD is linked with an increase in bone matrix mineralisation due to suppression of bone remodelling. However, denosumab anti-resorptive action also leads to an increase in fatigue microdamage, which may ultimately lead to an increased fracture risk. A novel mechanobiological model of bone remodelling was developed to investigate how these counter-acting mechanisms are affected both by exercise and long-term denosumab treatment. This model incorporates Frost's mechanostat feedback, a bone mineralisation algorithm and an evolution law for microdamage accumulation. Mechanical disuse and microdamage were assumed to stimulate RANKL production, which modulates activation frequency of basic multicellular units in bone remodelling. This mechanical feedback mechanism controls removal of excess bone mass and microdamage. Furthermore, a novel measure of bone local failure due to instantaneous overloading was developed. Numerical simulations indicate that trabecular bone volume fraction and bone matrix damage are determined by the respective bone turnover and homeostatic loading conditions. PMO patients treated with the currently WHO-approved dose of denosumab (60 mg administrated every 6 months) exhibit increased BMD, increased bone ash fraction and damage. In untreated patients, BMD will significantly decrease, as will ash fraction; while damage will increase. The model predicted that, depending on the time elapsed between the onset of PMO and the beginning of treatment, BMD slowly converges to the same steady-state value, while damage is low in patients treated soon after the onset of the disease and high in patients having PMO for a longer period. The simulations show that late treatment PMO patients have a significantly higher risk of local failure compared to patients that are treated soon after the onset of the disease. Furthermore, overloading resulted in an increase of BMD, but also in a faster increase of damage, which may consequently promote the risk of fracture, specially in late treatment scenarios. In case of mechanical disuse, the model predicted reduced BMD gains due to denosumab, while no significant change in damage occurred, thus leading to an increased risk of local failure compared to habitual loading.

1. Introduction

Denosumab treatment of patients with post-menopausal osteoporosis (PMO) has been shown to increase bone mineral density (BMD) as assessed by dual-energy X-ray absorptiometry (DXA) and it was observed that denosumab reduced the risk of new radiographic vertebral fractures by 68%, with the risk of hip fractures and non-vertebral fractures decreasing by 40 and 20%, respectively (Cummings et al., 2009). The increase in BMD has been attributed to a decrease in bone turnover and the associated increase in the degree of bone mineralisation (Dempster et al., 2018), as the osteoclastic activity dissolves the bone matrix and prevents it from being mineralised further. Despite its high efficacy in the treatment of PMO, long-term (i.e., >4 years) treatment with denosumab has also been linked with certain risks. Among those, development of atypical femoral fractures (AFF) is of significant concern and has been associated with the accumulation of microcracks in the bone matrix due to suppression of bone remodelling (Aspenberg, 2014), which is detected by a decrease in the levels of bone turnover markers (Miller et al., 2008).

Action of anti-resorptive drugs has been extensively studied experimentally and linked to decreases in bone turnover and suppression of osteoclastic activity together with observed changes in BMD (Langdahl, 2020). Riggs and Parfitt (2005) suggested that there are potentially three mechanisms on how anti-resorptive drugs exert increases in BMD: (i) increases of bone mineralisation, (ii) accumulation of microcracks in the bone matrix, and (iii) positive bone balance due to net increase of osteoblasts compared to osteoclasts. The latter mechanism which theoretically could lead to overfilling of resorption cavities in trabecular bone has more recently been ruled out due to the fact that bone remodelling is a coupled process, and suppression of osteoclastic activity also leads to reduction of osteoblastic activity (Sims and Martin, 2015). Hence, understanding on how the first two mechanisms interact and regulate BMD is central to understanding action of anti-resorptive drugs. However, it is currently not known to which extent the two mechanisms of increases of bone mineralisation and accumulation of microcracks contribute to net BMD increases. The aim of the current paper is to add to our understanding of the relative contributions of these competing mechanism which can provide new insights into the efficacy and safety of long-term treatment of OP with denosumab (and other anti-resorptive drugs). In particular, this new knowledge would help design new treatment regimens with respect to drug duration, dose magnitude, and dosing intervals.

Increases in BMD due to anti-resorptive treatments are explained by the conceptual model of bone mineralisation. This model links the rate of bone remodelling (i.e., bone turnover) with the degree of bone tissue mineralisation (BTM) (Bala et al., 2013; Boivin and Meunier, 2002). Bone mineralisation has two phases: a fast primary phase, which takes place over several days to weeks and achieves a degree of mineralisation of approximately 70% and a slow secondary phase, which can take from months to years and may achieve degrees of mineralisation of up to 95%. On the other hand, mineral is removed from the bone matrix by osteoclastic action, which dissolves it, returning it to the bloodstream. In this manner, the model predicts that bone sites undergoing high turnover are characterised by a lower BTM (and BMD) based on the fact that continuous remodelling prevents excessive secondary mineralisation. On the contrary, at sites of low turnover there is sufficient time for secondary mineralisation to occur and the mineral content is quite high, as occurs for example in interstitial bone (O'Brien et al., 2000).

The mechanism of targeted remodelling was first discovered in the late 90s looking at histological sections of normal, healthy cortical bone which showed microcracks around osteons which have been associated with exposure to dynamic habitual loading (e.g., walking) (Burr et al., 1985; Parfitt, 2002). In engineering, these microcracks are commonly referred to as fatigue microcracks, which accumulate in a material exposed to dynamic loading, eventually coalescing into macrocracks and leading to fracture (i.e., structural failure). Targeted bone remodelling is the process by which these microcracks are removed from the bone matrix in order to avoid development of fatigue failure or occurrence of stress fractures. While originally observed mostly in cortical bone, fatigue damage may also occur in trabecular bone (Dendorfer et al., 2006; Rapillard et al., 2006). However, in trabecular bone microcracks cannot easily accumulate due to the fact that cancellous bone exhibits a higher bone turnover compared to cortical bone. Fatigue loading induces formation of microcracks in areas of cortical bone which are subsequently resorbed (Verborgt et al., 2000). Furthermore, it was shown that both mechanical disuse and fatigue loading increase osteocyte apoptosis in specific bone regions (Verborgt et al., 2000; Aguirre et al., 2006). Indeed, tail-suspension stimulates osteocyte apoptosis, which is followed by bone resorption targeted to areas containing the apoptotic osteocytes in mice (Aguirre et al., 2006).

Bone remodelling has been also shown to play a key role in calcium and phosphorus homeostasis (Peterson and Riggs, 2010). The osteoclastic action that returns bone mineral into the bloodstream is enhanced by the secretion of parathyroid hormone (PTH), which controls calcium homeostasis at several levels and is increased if calcium deficiency is detected in the serum. In such case, calcium must be retrieved from bone matrix, where it is stored. Martínez-Reina et al. (2008) hypothesised that calcium retrieval is potentially more effective if it takes place at highly mineralised bone sites. This fact could be linked to targeted bone remodelling, as highly mineralised bone also accumulates a larger amount of damage (O'Brien et al., 2000; Qiu et al., 2005). Thus, the target of bone remodelling would not only be repairing damage, but also returning calcium to the bloodstream as efficiently as possible.

The mechanobiological link between microdamage and remodelling was established via discovery of increased remodelling around apoptotic osteocytes in rat ulnar fatigue-loading experiments (Verborgt et al., 2000). In the latter case, inhibition of osteocyte apoptosis prevents the intra-cortical resorption that occurs in response to microcracks (Cardoso et al., 2009), suggesting that osteocyte apoptosis controls osteoclast recruitment to the damaged area. Consistent with this idea, Tatsumi et al. (2007) have demonstrated that stimulation of osteocyte apoptosis, in and of itself, is sufficient to stimulate bone resorption that is associated with an increase in RANKL production in bone, but the cellular source of the RANKL was not determined.

There are limited data on bone quality for patients treated with denosumab compared to other anti-resorptive drugs. The pharmacodynamics of denosumab is different to that of other drugs, but it is their anti-resorptive action on bone remodelling that results in similar effects on bone quality. A recent study on the use of bisphosphonate (BP) treatment showed that the anti-resorptive therapy did not result in a detectable mechanical benefit in the trabecular bone specimens (from hip fracture patients) examined (Jin et al., 2017). Instead, BP use was associated with substantially reduced bone strength. This low strength may be due to the greater accumulation of microcracks and a lack of any discernible improvement in bone volume or microarchitecture. That study suggested that the clinical impact of BP-induced microcrack accumulation may be significant. Major current limitations on detecting significant effects of anti-resorptive treatments on microcrack density in bone is due to the fact that bone biopsies are only taken at non-load bearing bone sites, e.g., iliac crest, hence, masking the effect of mechanical loading and anti-resorptive therapy. However, it is very plausible from a material science point of view to hypothesise that a higher mineralised bone matrix (as observed with anti-resorptive treatment) will accumulate more microcracks during dynamic loading.

Based on the above described mechanisms observed for action of denosumab treatment of PMO, we have developed a comprehensive model of bone remodelling incorporating the effect of bone mineralisation, microdamage, and mechanobiological feedback. The current model of bone remodelling is an extension of our previous mechanistic PK-PD model of the effects of denosumab treatment on PMO (Martínez-Reina and Pivonka, 2019; Martínez-Reina et al., 2021) with respect to accounting for the accumulation of microdamage in the bone matrix. This model incorporates the relation between bone turnover and BMD, together with the bone mineralisation process and the formation of microdamage in the bone matrix. The evolution law for damage accumulation is formulated within the framework of continuum damage mechanics (Lemaitre and Chaboche, 1990). The bone remodelling model accounts for bone cell interactions via the RANK-RANKL-OPG pathway, the action of TGF–β and mechanobiological feedback (Martínez-Reina and Pivonka, 2019; Martínez-Reina et al., 2021; Pivonka et al., 2008, 2010; Scheiner et al., 2013). Mechanical overuse is simulated via increase of osteoblast precursors proliferation (Pivonka et al., 2012; Scheiner et al., 2013, 2014), while mechanical disuse is simulated via RANKL production by osteoblast precursor cells. Furthermore, bone matrix damage was linked to increased RANKL production by (apoptotic) osteocytes. The mineralisation model takes into account the balance of mineral within bone tissue and is based on the work of Martínez-Reina et al. (2008). As in previous studies, the PK model of denosumab is a one compartment model including a drug saturation term for high doses (see Marathe et al., 2011).

Utilising this model, we investigate a variety of treatment scenarios with emphasis on combined effects of mechanical loading (including overuse and disuse) together with denosumab treatment in PMO.

2. Mechanistic PK-PD Model of Bone Adaptation Including Damage

2.1. Model of Bone Cell Interactions in Bone Adaptation

A brief description of the mathematical model describing bone cell interactions is provided. As in previous models, the RANK-RANKL-OPG pathway, together with the action of several regulatory factors on bone cells, including TGF–β, PTH, and mechanobiological feedback is given (for details on original models see Pivonka et al., 2008, 2010; Scheiner et al., 2013; Pivonka et al., 2012; Martínez-Reina and Pivonka, 2019). The new model has been designed following the structure of the original model, adding the population of osteocytes, as done in Martin et al. (2019), slightly modifying the mechanoregulation feedback and adding new features relevant to the formulation of damage, the last two modifications being dealt with in subsections 2.3 and 2.4, respectively.

Following the approach taken by Pivonka et al., the bone adaptation process can be described as cell balance equations. The bone cell types (i.e., state variables) considered in the current model are: (i) osteoblast precursor cells (Ob_p), (ii) active osteoblasts (Ob_a), (iii) active osteoclasts (Oc_a), and (iv) osteocytes (Ot). The cell pools of uncommitted osteoblasts (Ob_u) and osteoclast precursors (Oc_p) are assumed constant:

$\begin{array}{ll}\begin{array}{l}\frac{dO{b}_p}{dt}=D_{O{b}_u}·{\Pi }_{act}^{TGF-\beta }·O{b}_u+P_{O{b}_p}·{\Pi }_{act}^{{\psi }_{bm}}·O{b}_p\\ -D_{O{b}_p}·{\Pi }_{rep}^{TGF-\beta }·O{b}_p;\end{array}& (1)\end{array}$

$\begin{array}{ll}\frac{dO{b}_a}{dt}=D_{O{b}_p}·{\Pi }_{rep}^{TGF-\beta }·O{b}_p-A_{O{b}_a}·O{b}_a;& (2)\end{array}$

$\begin{array}{ll}\frac{dO{c}_a}{dt}=D_{O{c}_p}·{\Pi }_{act}^{RANKL}·O{c}_p-A_{O{c}_a}·{\Pi }_{act}^{TGF-\beta }·O{c}_a;& (3)\end{array}$

$\begin{array}{ll}\frac{dOt}{dt}=\eta \frac{d{f}_{bm}}{dt}& (4)\end{array}$

where D_Obu, D_Obp, and D_Ocp are the differentiation rates of Ob_u, Ob_p, and Oc_p respectively. The second term in the right-hand side of Equation (1) corresponds to proliferation of osteoblast precursors and P_Obp gives the maximum proliferation rate. A_Oba and A_Oca are the apoptosis rates of Ob_a and Oc_a respectively. The variables ${\Pi }_{act}^{TGF-\beta }$, and ${\Pi }_{rep}^{TGF-\beta }$ represent activator and repressor functions related to the binding of TGF–β to its receptor. Similarly, ${\Pi }_{act}^{RANKL}$ is the activator function related to the RANK-RANKL binding. ${\Pi }_{act}^{{\psi }_{bm}}$ is a function of the mechanical stimulus that regulates the anabolic part of the mechanobiological feedback in the proliferation term and will be addressed in section 2.3. Finally, η is the concentration of osteocytes in bone matrix, which is assumed constant as in Martin et al. (2019), thus leading to proportional variations of osteocytes population and fraction of bone matrix volume per total volume, f_bm (Equation 4). The variation of f_bm over time is precisely one of the main outcomes to be derived from the set of cell population equations and is defined through the balance between resorbed and formed tissue:

$\begin{array}{ll}\frac{d{f}_{bm}}{dt}=-k_{res}·O{c}_a+k_{form}·O{b}_a;& (5)\end{array}$

where k_res and k_form are, respectively, the bone matrix volume resorption rate and osteoid volume formation rate. This distinction is important with regard to the mineralisation algorithm, since the bone matrix resorbed by active osteoclasts is mineralised, while the osteoid deposited by active osteoblasts contains no mineral. Cell balance equations (Equations 1–4) are composed of a production term and a degradation one, which describes differentiation of one cell type into another (or terminal cell fate, i.e., apoptosis). A schematic figure of the mechanistic PK-PD model is presented in Figure 1. Model parameters of the cell population model are given in Table A1.

FIGURE 1

Figure 1. Schematic representation of the mechanistic PK-PD model: bone cell differentiation stages along with biochemical and biomechanical interactions are presented. Subcutaneous injection of denosumab leads to distribution of the drug into the central compartment where it interacts with the RANK-RANKL-OPG pathway (red arrow between Ob_p and Oc_p). The latter interactions are accounted for via competitive binding reactions. The mineralisation of osteoid is shown in orange.

2.2. Denosumab Action on RANK-RANKL-OPG: Competitive Binding

The action of denosumab on bone adaptation is included via competitive binding reactions within the RANK-RANK-OPG pathway (Marathe et al., 2008, 2011; Scheiner et al., 2014). In these models, action of denosumab is taken into account via the RANKL activator function ${\Pi }_{act}^{RANKL}$ (Equation 3, first term on the right). Denosumab competes with RANK (and OPG) for binding to RANKL. Thus, higher concentrations of denosumab give rise to lower concentrations of RANKL-RANK complexes and, hence, lower values of ${\Pi }_{act}^{RANKL}$. Adapting the approach of Scheiner et al. (2014),

$\begin{array}{l}RANKL=RANK{L}_{eff}\frac{{\beta }_{RANKL}+P_{RANKL}}{{\beta }_{RANKL}+{\tilde{D}}_{RANKL}RANK{L}_{eff}}·\\ ·[1+\frac{OPG}{K_{d,\left[RANKL-OPG\right]}}+\frac{RANK}{K_{d,\left[RANKL-RANK\right]}}\\ +\frac{\zeta C_{p,den}}{K_{d,\left[RANKL-den\right]}}{]}^{-1};& (6)\end{array}$

where K_{d, [RANKL−OPG]}, K_{d, [RANKL−RANK]}, and K_{d, [RANKL−den]} are the equilibrium dissociation binding constants for binding of OPG, RANK and denosumab to RANKL. OPG, RANK, and RANKL are the concentrations of respective regulatory factors in the bone tissue compartment, while C_p,den is the concentration of denosumab in the central compartment (see Equation 36 in Appendix) and ζ is the accessibility factor of denosumab from the central compartment to the bone tissue compartment¹. In the original model (Lemaire et al., 2004) all concentrations were formulated with respect to a pseudo central compartment and, consequently, no distinction between site-specific bone tissue compartments was needed. However, formulation of mechanobiological PK-PD models requires specification of a particular bone site which is exposed to physiological mechanical loading. ${\tilde{D}}_{RANKL}$ is the RANKL degradation rate, P_RANKL provides the RANKL production rate induced by PMO and mechanical underloading. β_RANKL is the production rate of endogenous RANKL on the surface of osteoblasts precursors and osteocytes. We have assumed that RANKL is expressed by those cells, following experimental evidence (Nakashima et al., 2011; Xiong and O'Brien, 2012) and a previous model (Martin et al., 2019). So, RANKL_eff is the total effective carrying capacity of those cells that controls the maximum expression of RANKL:

$\begin{array}{ll}RANK{L}_{eff}=R_{RANKL}·O{b}_p·{\Pi }_{act}^{PTH}+R_{RANKL}·Ot·{\Pi }_{act}^{dam}& (7)\end{array}$

where R_RANKL is the carrying capacity of the individual cells of both types, that we have assumed equal. We have also assumed that the expression of RANKL on the surface of osteoblast precursors is upregulated by PTH, following previous models (Pivonka et al., 2008, 2010), and we have introduced in the present model an upregulation factor of RANKL expression by osteocytes due to damage, through a sigmoidal function:

$\begin{array}{ll}{\Pi }_{act}^{dam}=\frac{d^{\xi }}{d^{\xi }+{\delta }_{50}^{\xi }}& (8)\end{array}$

where d denotes the damage variable, described more in detail in section 2.4. The shape factor, ξ = 3, and the value of damage leading to a 50% of the maximum response, δ₅₀ = 0.1, were chosen in such a way that both terms in Equation 7 are typically of the same order of magnitude². Verborgt et al. (2000) showed that osteocyte apoptosis occurs after fatigue induced bone matrix damage. Moreover, they found that osteocyte apoptosis was highly localised to sites of microdamage that are subsequently remodelled. Osteocytes in the vicinity of a microcrack would express both Bax (a proapoptotic gene product) and Bcl-2 (an antiapoptotic gene product), with the peak of Bax expression observed immediately at the microcrack locus and the peak of Bcl-2 expression at some distance (1–2 mm) from microcracks (Verborgt et al., 2002). Seemingly, distant osteocytes would protect themselves from matrix injury induced cell death, thereby exercising an additional level of control in the regulation of osteocyte apoptosis and bone remodelling. This expression of apoptotic signals would be related to the expression of come and eat me signals (Jin and El-Deiry, 2005) to attract macrophages to the site of apoptotic osteocytes. Kurata et al. (2006) showed later that focal wounding of osteocyte-like cells (MLO-Y4) in vitro triggered release of RANKL and macrophage-colony stimulating factor (M-CSF), although whether these key signaling molecules come from dying cells or the non-apoptotic surviving cells was not examined. Here we have assumed by using Equation (7) that apoptotic osteocytes near microcracks would express RANKL so causing their surrounding bone matrix to be resorbed.

In Equation (6), the RANKL production rate P_RANKL is given by two terms that define the contribution of mechanical underloading, $P_{RANKL}^{mech}$, and a disease-related increase in RANKL production over time, $P_{RANKL}^{PMO}$:

$\begin{array}{ll}{P}_{RANKL}=P_{RANKL}^{mech}+P_{RANKL}^{PMO}& (9)\end{array}$

The first term is explained in section 2.3, while the term $P_{RANKL}^{PMO}$ is a consequence of the onset of menopause though increasing gradually over time, through the following sigmoidal function:

$\begin{array}{ll}{P}_{RANKL}^{PMO}(t)=P_{RANKL}^{PMO,max}\frac{{(t-t_{onset})}^2}{{(t-t_{onset})}^2+{\delta }_{PMO}^2}\text{\hspace{1em}for}t>t_{onset}& (10)\end{array}$

where $P_{RANKL}^{PMO,max}$ is the maximum (long-term) RANKL production rate due to PMO, t_onset is the time of onset of the disease and δ_PMO is a time constant that establishes when the 50% of $P_{RANKL}^{PMO,max}$ is reached.

Finally, the activator function of RANKL in Equation (3) can be expressed as:

$\begin{array}{ll}{\Pi }_{act}^{RANKL}=\frac{RANKL·RANK}{K_{d,\left[RANKL-RANK\right]}+RANKL·RANK};& (11)\end{array}$

The concentrations of denosumab, RANK, OPG, PTH, and TGF−β along with the model parameters of the binding reactions are provided in the Appendix.

2.3. Mechanoregulation

The model includes the mechanical feedback regulation of bone through the Mechanostat Theory proposed by Frost (2003) (see Figure 2). This theory postulates the existence of 4 zones or “windows” in Frost's terminology: (1) disuse window, where net bone loss is observed for a low level of “Minimally Effective Strains” (MES) or other stimuli; (2) adapted window, where no net effect of BMUs on bone mass is seen for intermediate values of MES; (3) mild overload window, where net bone formation occurs for high MES, and (4) pathologic overload window, leading to fracture, for very high values of MES. This last window is not directly considered in the mechanical regulation through the definition of ${\Pi }_{act}^{{\psi }_{bm}}$ (see Figure 3), but indirectly through the accumulation of microstructural damage.

FIGURE 2

Figure 2. The Mechanostat Theory (adapted from Frost, 2003).

FIGURE 3

Figure 3. Function of proliferation of osteoblasts precursors that establishes the relation between the anabolic factor ${\Pi }_{act}^{{\psi }_{bm}}$ and the SED.

Mechanical disuse is assumed to enhance the production of RANKL on osteoblasts precursors, through the term $P_{RANKL}^{mech}$, which is modulated by the strain energy density (SED) of bone matrix, designated as ψ_bm:

$\begin{array}{ll}{P}_{RANKL}^{mech}=\{\begin{array}{ll}{P}_{RANKL}^{mech,max}\left(1-\frac{{\psi }_{bm}}{{\psi }_r}\right)& \text{for\hspace{1em}}{\psi }_{bm}<{\psi }_r\\ 0& \text{for\hspace{1em}}{\psi }_{bm}\ge {\psi }_r\end{array}& (12)\end{array}$

where ψ_r is the SED below which underuse increases RANKL production and $P_{RANKL}^{mech,max}$ is the maximum RANKL production rate due to underuse. RANKL production is upregulated by PTH and downregulated by nitric oxide (NO), which is produced by osteocytes and, in turn, upregulated by mechanical stimulus. However, this process is only indirectly considered through Equation (12), which assumes a maximum RANKL production rate for total disuse.

Overload is assumed to promote bone formation by proliferation of osteoblasts precursors through the activator function ${\Pi }_{act}^{{\psi }_{bm}}$, which is given by the piecewise linear function of SED defined in Figure 3. The less steep piece of the function would correspond to the disuse and adapted windows of the Mechanostat Theory, where bone formation is not particularly promoted. In the case of the disuse window, this would be added to the increased RANKL production ($P_{RANKL}^{mech}$). Obviously, the steepest piece of function ${\Pi }_{act}^{{\psi }_{bm}}$ would correspond to the overload window.

The SED, ψ_bm, was used as a measure of the mechanical stimulus sensed by bone cells to drive bone adaptation, as traditionally done in the literature (Beaupré et al., 1990; Huiskes et al., 1987). ψ_bm was used here as an alternative to the strains MES, used in the Mechanostat Theory (Frost, 2003). In a uniaxial stress state they are related by:

$\begin{array}{ll}{\psi }_{bm}=\frac{1}{2}E·ME{S}^2& (13)\end{array}$

being E the Young's modulus. The parameter ψ_r in Equation (12) as well as ψ_bm1 and ψ_bm2, used in Figure 3 to define function ${\Pi }_{act}^{{\psi }_{bm}}$, were defined using Equation (13), respectively with MES_r = 1,000 με, MES₁ = 800 με, and MES₂ = 1,600 με. These values and the values of ${\Pi }_{act}^{{\psi }_{bm}}$ corresponding to ψ_bm1 and ψ_bm2 were adjusted to reproduce the Mechanostat Theory along with the Principle of Cellular Accomodation (see Appendix).

2.4. Damage of Bone Matrix

In this section we address how to estimate microstructural damage, that we assumed to drive bone remodelling through Equation (7) and to affect mechanical properties as will be discussed shortly. It has long been hypothesised that one of the major functions of bone remodelling is to remove microcracks from the bone matrix and so to avoid accumulation of the latter, which could result in macroscopic failure. One way to describe the accumulation of microcracks in a particular volume of material is via use of Continuum Damage Mechanics (Lemaitre and Chaboche, 1990). The latter theory introduces a damage variable, d, which is linked to the density of microcracks in a volume of material and to the loss of stiffness through Equation (14). This is variable is such that d ∈ [0, 1], with d = 0 corresponding to an undamaged state and d = 1 to a local fracture or failure situation:

$\begin{array}{ll}\text{C}=(1-d){\text{C}}_0& (14)\end{array}$

where C and C₀ are, respectively, the stiffness tensors of damaged and undamaged bone (Lemaitre and Chaboche, 1990)³.

Microdamage accumulates in the bone matrix due to fatigue loading and is repaired by bone remodelling, as osteoclasts resorb the damaged tissue while the osteoid deposited by osteoblasts is initially intact. The evolution law for damage can be expressed as:

$\begin{array}{ll}\dot{d}={\dot{d}}_A-{\dot{d}}_R& (15)\end{array}$

where $\dot{d}$_A is the rate of damage accumulation by fatigue loading and $\dot{d}$_R is the rate of damage removal by bone remodelling. The latter is assessed by assuming that damage is uniformly distributed throughout the representative volume element (RVE). Thus, the amount of damage repaired by remodelling is proportional to the damage present in that volume and to the volume of tissue being resorbed (see Equation 5), through the fraction that this volume represents within the bone matrix volume:

$\begin{array}{ll}{\dot{d}}_R=d\frac{k_{res}·O{c}_a}{f_{bm}}& (16)\end{array}$

Damage accumulation is evaluated following the procedure described in Martínez-Reina et al. (2008) and Martínez-Reina et al. (2009), which, in turn, are based on the works by Pattin et al. (1996) and García-Aznar et al. (2005). Experimental fatigue tests provide the evolution of d with the strain or the stress level and the number of cycles (Pattin et al., 1996), as well as fatigue life, N_f, which is typically given by expressions such as:

$\begin{array}{ll}\begin{array}{ll}{N}_f=\frac{K_i}{{\epsilon }^{{\delta }_i}}& i=c(\text{compresion}),t(\text{tension})\end{array}& (17)\end{array}$

where K_i and δ_i stand for constants that are different in tension and compression and ε is the uniaxial strain expressed in με = μm/m. García-Aznar et al. (2005) correlated Equation (17) with the experimental results obtained by Pattin et al. (1996) to get: $K_c=9.333·10^{40}$ and δ_c = 10.3 in compression, and $K_t=1.445·10^{53}$ and δ_t = 14.1 in tension. The loss of stiffness E/E₀ was also experimentally measured by those authors as a function of the applied constant strain and the number of cycles. Again, García-Aznar et al. (2005), fitted the experimental curves obtained by Pattin and co-workers with the following expressions:

$\begin{array}{ll}{d}_c=-\frac{1}{{\gamma }_c}\left[ln(1-C_c{\epsilon }^{{\delta }_c}N)\right]& (\mathrm{18a})\end{array}$

$\begin{array}{ll}{d}_t=1-{\left[\frac{1}{C_{t2}}ln(e^{C_{t2}}-C_{t1}{\epsilon }^{{\delta }_t}N)\right]}^{\frac{1}{{\gamma }_t}}& (\mathrm{18b})\end{array}$

where N is the number of cycles and C_t2 = −20 was fitted from the experimental curves along with:

$\begin{array}{l}{\gamma }_c=-5.238(\epsilon -6100)10^{-3}+7;\\ C_c=\frac{1-e^{-{\gamma }_c}}{K_c}\text{in compression}\\ {\gamma }_t=-0.018(\epsilon -4100)+12;C_{t1}=\frac{e^{C_{t2}}-1}{K_t}\text{in tension}& (19)\end{array}$

In the damage model proposed by Martínez-Reina et al. (2009), cracks were assumed to grow normal to the maximum strain direction and only under tensile strains. This allows to apply the model to a general strain state, by replacing ε with the maximum principal strain, ε_max. The tests performed by Pattin et al. (1996) and fitted with Equation (18) were conducted under constant strain. However, they can be applied to a general loading history using the procedure described in Martínez-Reina et al. (2009) and explained next.

Let us assume that, at a given moment, damage is equal to d and a maximum principal strain ε_max is applied N cycles in the next step. Let us calculate the increment of damage accumulated by fatigue, Δd_A, with those cycles. Likely, the current damage was not produced by a constant strain ε_max, but we can assume that it was so without loss of generality. Then, we can use Equation (18b) to work out the number of cycles Ñ that would have been needed to reach the current damage d with the current strain ε_max.

$\begin{array}{ll}d=1-{\left[\frac{1}{C_{t2}}ln(e^{C_{t2}}-C_{t1}{\epsilon }_{max}^{{\delta }_t}Ñ)\right]}^{\frac{1}{{\gamma }_t}}\Rightarrow Ñ& (20)\end{array}$

The increment of damage Δd_A would have been reached with the additional N cycles applied at the present step and can be assessed from:

$\begin{array}{ll}d+\Delta d_A=1-{\left[\frac{1}{C_{t2}}ln\left(e^{C_{t2}}-C_{t1}{\epsilon }_{max}^{{\delta }_t}(Ñ+N)\right)\right]}^{\frac{1}{{\gamma }_t}}& (21)\end{array}$

This procedure allows working out the increment of damage, Δd_A, but requires that Equation (15) be rewritten in incremental form and integrated using an explicit integration scheme, as done in Martínez-Reina and Pivonka (2019) and Martínez-Reina et al. (2021).

2.5. Degradation of Fatigue Properties With the Mineral Content

The mineral phase contributes to increase the stiffness of bone, but also makes it more brittle (Currey, 2004). As far as we know, no experimental study has provided a correlation between mineral content and bone fatigue properties, though some studies have confirmed that interstitial bone, with the highest mineral content, is where microcracks can be more easily found (Boyce et al., 1998; O'Brien et al., 2000; Qiu et al., 2005). For this reason, we have followed the damage model previously proposed by Martínez-Reina et al. (2009) in which fatigue properties are degraded as the mineral content rises. According to this idea the following assumptions are made in the model (see Martínez-Reina et al., 2009 for more details):

1. The shape of d − N curves, expressed by the Equations (18), is maintained regardless of the mineral content.

2. Only the fatigue life is affected by the mineral content, by redefining K_t in Equation (17), while keeping constant the exponent δ_t. This modifies the d − N law, as C_t1 depends on K_t. In this way, increasing K_t results in a longer fatigue life and a slower damage accumulation rate.

3. A life of 10⁷ cycles was assigned to the fatigue limit. This fatigue limit is usually assumed to occur for a given fraction of the ultimate tensile strain, ε_u/β, where the parameter β may depend on the type of material (Juvinall, 1967). So, K_t was obtained from Equation (17) as:

$\begin{array}{ll}{K}_t([Ca])=10^7{\left(\frac{{\epsilon }_u([Ca])}{\beta }\right)}^{{\delta }_t}& (22)\end{array}$

Here, the most typical value β = 2 was chosen following Martínez-Reina et al. (2009).

4. As Currey showed (Currey, 2004), ε_u depends on the calcium concentration of bone matrix, [Ca]. The following regression ε_u = ε_u([Ca]) was fitted in Martínez-Reina et al. (2009) from the experimental results presented by Currey (2004):

$\begin{array}{ll}log{\epsilon }_u=31.452-11.341log\left[Ca\right]& (23)\end{array}$

where ε_u is expressed in με and the concentration [Ca] is expressed in mg of calcium per g of bone matrix and is related to the ash fraction, α, the ratio between the ash mass and the dry mass (see Appendix for details). More precisely, the relation [Ca] = 398.8 · α was assumed, based on the molecular weigths of hydroxyapatite and type I collagen. Equations (20), (22), and (23) allowed to define C_t1 as a function of α to be used in (18b) and the related equations.

The value fitted by García-Aznar et al. (2005), $K_t=1.445·10^{53}$, corresponds to a normal value of ash fraction, α = 0.72. The importance of the degradation of fatigue properties with the mineral content on the risk of local failure will be evaluated in a simulation in which two cases will be compared: (1) a model implementing Equations (22) and (2) a model using $K_t=1.445·10^{53}$ constant.

2.6. Bone Apparent Density and Stiffness

Bone apparent density changes as a consequence of the variation of porosity, accounted by Equation (5), and mineralisation. The latter process controls tissue density, ρ_t, given by:

$\begin{array}{ll}{\rho }_t=\frac{m}{V_{bm}}={\rho }_m{v}_m+{\rho }_o{v}_o+{\rho }_w{v}_w& (24)\end{array}$

where m and V_bm are, respectively, the mass and volume occupied by bone matrix, while v_m, v_o, and v_w stand for the specific volumes of the three phases that compose bone matrix (namely, mineral, organic, and water) and ρ_i stand for the corresponding densities. While v_o can be assumed constant, v_m and v_w vary throughout the mineralisation process (see Appendix for details). Bone apparent density is then given by porosity (or alternatively bone volume fraction) and bone tissue density:

$\begin{array}{ll}\rho =\frac{m}{V_{bm}}\frac{V_{bm}}{V_{RVE}}={\rho }_t{f}_{bm}& (25)\end{array}$

Finally, bone stiffness is needed to assess SED. In this study, we have assumed that bone tissue is an isotropic material with a Poisson's ratio ν = 0.3 and a Young's modulus given in MPa by the following correlations:

$\begin{array}{ll}E(\rho ,d)=\{\begin{array}{ll}2014{\rho }^{2.5}(1-d)& \text{\hspace{1em}if }\rho <1.2g/c{m}^3\\ 1763{\rho }^{3.2}(1-d)& \text{\hspace{1em}if }\rho \ge 1.2g/c{m}^3\end{array}& (26)\end{array}$

These expressions are based on the correlations experimentally obtained by Jacobs (1994), which are multiplied by the factor (1 − d), to consider microstructural damage as usually done in Continuum Damage Mechanics (Lemaitre and Chaboche, 1990) (recall Equation 14).

2.7. Estimation of the Risk of Local Failure

The likelihood of osteoporotic patients suffering a fracture depends on many factors, such as: bone apparent density, amount of microstructural damage, mineral content and brittleness of bone matrix, trabecular microarchitecture, magnitude and orientation of the load, among many others. For this reason, evaluating the risk of fracture or the strength of a certain bone is a complex task. Several works have addressed this problem at a structural level, estimating the bone strength and/or failure patterns using a FE modelling approach (Hambli, 2013a,b; Harrison et al., 2013; Fan et al., 2016; Hambli et al., 2016), but it is out of the scope of the present work.

Having recognised this limitation, it would be interesting, however, to have a tool to compare the risk of failure in different scenarios, at least at a local, i.e., material point of view. To this end, we have defined a variable to estimate the risk of failure by taking into account only bone apparent density, damage, mineral content and magnitude of the load from the factors referred to above. In order to define that variable we have used the equations presented in subsection 2.4.

More precisely, we have defined the variable NDtF(σ, t), which gives the number of days needed to reach local failure (d = 1) at a given instant, t, if an overload consisting in a uniaxial stress σ is applied from that moment on. We use the variables f_bm(t), d(t), Oc_a(t), and α(t) corresponding to that instant and assume that all of them, except damage, remain constant until failure, which is a strong simplification, especially if NDtF yields a high value. Next, we define a time variable τ commencing at the instant t and describing the evolution that damage would follow until failure if the conditions at time t were kept constant, i.e., d(t + τ). NDtF is the number of days that fulfills the following equation d(t + NDtF) = 1.⁴ The algorithm used to assess NDtF(σ, t) consists of the following steps:

1. Start from f_bm(t), d(t), Oc_a(t), and α(t).

2. Apparent density, ρ(t), is obtained from f_bm(t) and α(t) (see Appendix).

3. For a given instant t + τ, d(t + τ) is known (commencing with τ = 0 and d(t)). Thus, evaluate the Young's modulus from Equation (26), so to obtain E(ρ(t), d(t + τ)).

4. Assuming a uniaxial stress state, calculate the maximum principal strain:

$\begin{array}{ll}{\epsilon }_{max}=\{\begin{array}{ll}\frac{\sigma }{E\left(\rho (t),d(t+\tau )\right)}& \text{if}\sigma \ge 0\\ -\frac{\nu \sigma }{E\left(\rho (t),d(t+\tau )\right)}& \text{if}\sigma <0\end{array}& (27)\end{array}$

5. With d(t + τ) and ε_max, use Equation (20) to work out Ñ.

6. Equation (15), (16), and (21) are used to update damage when N additional cycles (corresponding to Δτ = 1 day) are applied:

$\begin{array}{ll}\begin{array}{l}d(t+\tau +\Delta \tau )=1-{\left[\frac{1}{C_{t2}}ln\left(e^{C_{t2}}-C_{t1}{\epsilon }_{max}^{{\delta }_t}(Ñ+N)\right)\right]}^{\frac{1}{{\gamma }_t}}\\ -d(t+\tau )\frac{k_{res}·O{c}_a(t)·\Delta \tau }{f_{bm}(t)}\end{array}& (28)\end{array}$

7. If d(t + τ + Δτ) ≥ 0.999, then assign NDtF(σ, t) = τ + Δτ and exit the algorithm. Otherwise, go to step 3 and commence a new iteration.

It must be emphasised that NDtF is related to the risk of local failure (within the RVE) and that the risk of fracture at the organ level is influenced by other factors, such as the microstructure and the distribution of loads, mechanical properties, density, damage and mineral content throughout the bone among others. This makes it necessary to take into account the whole organ, for example with a FE model, in order to assess the true risk of fracture of the specific organ.

2.8. Homeostatic Initial Conditions

All the simulations started from an equilibrium point corresponding to a homeostasis situation under physiological (healthy) conditions. The governing equations constitute a set of differential-algebraic equations whose equilibrium points could be assessed by setting the derivatives equal to zero in Equations (1)-(5) and solving the resulting set of algebraic equations. However, the equations that model the mineralisation process are recursive and cannot be solved either in closed form or numerically. Alternatively, the equilibrium point was obtained by running a simulation for a given stress state and letting the system to reach homeostasis in the long-term. For a given stress, initial conditions were arbitrarily set in a first iteration. The values at the end of an iteration were used as the initial conditions in the next one and the process was repeated until convergence (two successive iterations producing differences less than 0.1% in every variable). We have assumed here a uniaxial stress state, which completely defines the homeostatic condition for a given value of the applied stress.

3. Results

The procedure outlined in subsection 2.8 was used to obtain the equilibrium (homeostatic) points for a range of tensile uniaxial stress. Figure 4 shows the homeostatic values of f_bm, strain and damage d for each value of the stress. Analogous graphs can be obtained for the rest of variables of the model and for compressive stress.

FIGURE 4

Figure 4. Homeostatic values of f_bm, strain and damage d for different values of uniaxial tensile stress.

A tensile stress of 12 MPa was chosen and the homeostatic state for that stress was obtained as explained above. This homeostatic state was used to define the set of initial conditions for all the simulations.

PMO was simulated by a gradual increase in the RANKL production rate (see Equation 10) where the parameter $P_{RANKL}^{PMO,max}$ controls the degree of incidence of the disease, starting from 600 pM day⁻¹ to simulate a moderate degree of incidence. The WHO-approved denosumab dosage for the treatment of osteoporosis is 60 mg administered every 6 months (60Q6) via subcutaneous injection, being higher doses restricted for the treatment of other diseases such as multiple myeloma or giant cell tumor. Figure 4 shows the evolution of bone apparent density when the treatment is commenced 3 years after the onset of the disease⁵ and the results are compared with the evolution in a non-treatment group. Patient's body weight (BW) influences the concentration of denosumab present in the central compartment (Marathe et al., 2011) and strongly affects the results. For this reason, a BW = 60 kg has been taken as the reference value.

In order to evaluate the risk of local failure that anti-resorptive treatments may cause in the long-term, the temporal evolution of microstructural damage, ash fraction and NDtF are also shown in Figure 5. It can be seen that damage and ash fraction increase simultaneously at the beginning of the treatment. This makes the risk of local failure to rise as well⁶. After a given time, as the mineral content gets stabilised, damage is slowed down and the risk of failure falls below the values obtained for the non-treatment group.

FIGURE 5

Figure 5. Evolution of density (top left), damage (bottom left), mineral content (bottom right), and NDtF corresponding to an overload σ = 24 MPa (top right) for the WHO-approved treatment 60Q6 and comparison with the non-treatment group.

The time elapsed from the onset of disease to the beginning of the treatment was varied from 1 to 10 years in a set of simulations (see Figure 6) aimed at highlighting the importance of starting the treatment as early as possible. It can be seen how late treatments can notably increase the risk of failure.

FIGURE 6

Figure 6. Influence of the time elapsed between the onset of disease and the beginning of the treatment. Evolution of density (left, solid lines), damage (left, dashed lines), and NDtF-σ = 24 MPa (right) for treatment 60Q6.

Figure 7 compares the effect of changes in the mechanical loading (i.e., physiological exercise), hereafter denoted as overloading, through the applied stress. Six cases were analysed, all except one (case 6, “no treatment”) including the 60Q6 treatment. In all the cases with treatment, this starts 3 years after the onset of the disease, except for case 5, and the change of load coincides with the beginning of the treatment. The cases are: (1) the nominal case with no change in mechanical load; (2) a stepwise increase of 1.2 MPa (10% of the homeostatic load) in the applied tensile stress; (3) a stepwise decrease of 1.2 MPa; (4) a stepwise increase of 3.6 MPa; (5) a stepwise increase of 3.6 MPa in a late treatment, commencing 10 years after the onset of the disease; (6) no treatment⁷. The results show that bone density rises with the increase of applied stress, though damage can also rise. In this regard, damage can be reduced with a decrease of the applied stress, but this does not imply an improvement of bone quality, since bone density and consequently stiffness may drop significantly, as in case 3, which exhibits a very low NDtF. In general, the simulations predicted that an increment of stress coincident with the treatment reduces the risk of local failure, though an excessive increment of stress could be dangerous in a late treatment, in which the initial bone condition can be very deteriorated by the prolonged disease.

FIGURE 7

Figure 7. Influence of an increment in the applied stress coinciding with the beginning of the treatment. Evolution of density (left, solid lines), damage (left, dashed lines), and NDtF-σ = 24 MPa (right) for treatment 60Q6. Four cases of changes in the applied stress are compared with the nominal case (no change) and the non-treatment group.

The importance of accumulation of unrepaired damage and degradation of fatigue properties with the mineral content on the risk of failure is compared next (see Figure 8). To this end, two cases were compared: (1) the model implementing the degradation of fatigue properties through Equation (22), as in the previous simulations, and (2) a model using $K_t=1.445·10^{53}$ constant. In the latter, only the accumulation of unrepaired damage plays a role; while in the former, both factors are important. This will be done for a late treatment 60Q6 beginning 10 years after the onset of the disease and coinciding with a stepwise increase of stress of two different magnitudes: (a) 3.6 MPa (30% of the nominal value) and (b) 4.8 MPa (40%). It can be seen that the embrittlement of bone matrix due to an excessive mineral content plays a key role in the failure experienced with a 30% overload, as the simulation without the degradation of properties does not predict that failure. However, in the case of a 40% overload, the accumulation of unrepaired damage seems to play a more important role, as the simulation without the degradation of properties also predicts that failure.

FIGURE 8

Figure 8. Influence of the degradation of fatigue properties with the mineral content. Evolution of density (left, solid lines), damage (left, dashed lines), and NDtF-σ = 24 MPa (right) for treatment 60Q6. Two scenarios are considered: (a) a stepwise increase of 30% of the nominal stress coinciding with a late treatment and (b) a stepwise increase of 40%; and two models with and without considering the degradation of fatigue properties with the mineral content. Note: the red and blue solid lines in the left figure coincide, respectively, with the cyan and green solid lines.

4. Discussion

4.1. General Comments

Hernandez et al. (2001b) developed a computational model of bone remodelling to compare the contributions of focal bone balance and mineralisation on BMD by simulating alendronate treatment using a bone balance method (decreased remodelling space, increased focal bone balance, uniform bone mineralisation) and a mineralisation method (decreased remodelling space, neutral focal bone balance, varying bone mineralisation). Their results suggested that the mineralisation method is more descriptive of long-term alendronate treatment indicating that adequate modelling of the mineralisation process is essential to explain observed BMD changes caused by alendronate. While the authors suggest that this model may be used to identify improved dosing regimens and to predict which osteoporosis treatments are more effective, we note that the latter model did not include mechanical loading and microcrack formation explicitly. Hence, the combined effects of physical exercise and treatment with anti-resorptive agents cannot be adequately investigated.

A previous study by Nyman et al. addressed a similar question as we posed in the current paper, i.e., how does mechanical loading and long-term anti-resorptive treatment affect BMD and risk of fracture (Nyman et al., 2004). They utilised a bone remodelling model incorporating micro-damage formation due to mechanical loading together with simulating effects of bisphosphonate treatment. Both disuse and fatigue microdamage were assumed to stimulate the activation frequency of basic multicellular units (BMUs) such that bone remodelling served to remove excess bone mass and microdamage. Bisphosphonate effects were simulated as suppression of BMU activation frequency either without a change in resorption by the BMU or with an independent decrease in resorption while the bone formation process was unaffected (i.e., formation initially exceeded resorption). We note that the initial increase of bone formation with respect to resorption is only temporarily and a very small contributor to increase in BMD in the long term (Martínez-Reina and Pivonka, 2019; Martínez-Reina et al., 2021). Based on the fact that the work of Nyman et al. did not include the mineralisation process they had to include the overfilling of resorption cavities as a mechanism of action of bisphosphonates in order to obtain BMD increases. Their model predicted a plateau in the bone mass gain that typically occurs in clinical studies of bisphosphonate treatment. However, the mechanism of overloading and the mineralisation process were not incorporated into that model.

The above limitations are overcome in the present study in which a new model has been developed based upon a previous one (Martínez-Reina et al., 2021) to include the level of microstructural damage as a new parameter regulating mechanobiological feedback. The also new proliferation function (see Figure 3) was calibrated to reproduce the Mechanostat Theory and the Principle of Cellular Accomodation. Once the new model was calibrated, it was used to simulate the response of bone to PMO, modelled as a gradual increase in RANKL expression, and subsequent treatments with denosumab. As discussed in a previous work (Martínez-Reina and Pivonka, 2019), bone density gain in denosumab treatments is mainly explained by the bone mineralisation process, which makes the mineral content to reach abnormally high values once bone turnover is blocked by the drug. But such a high mineral content jeopardises bone integrity as it makes bone matrix more brittle. Moreover, microstructural damage begins to accumulate and remains unrepaired due to the suppression or decrease of bone resorption. Figure 5 (bottom) showed how mineral content and damage rise simultaneously after the treatment commences, which is probably due to the concurrence of both factors: increased brittleness and unrepaired damage. In fact, damage is greater in the treatment group than in the non-treatment one for a long time, despite the greater BMD reached with the treatment, which increases stiffness and reduces strains. In the long-term, this reduction of strains predominates over the two negative factors commented before and damage falls below the values of the non-treatment group.

The risk of suffering an atypical femoral fracture (AFF) might also be enhanced by other factors not considered in this work. For example, at the micro (local) level, it could be affected by alterations of bone microarchitecture or by the higher concentrations of advanced glycation end-products within the extracellular collagen matrix, which can also raise bone brittleness (Vashishth et al., 2001). At the macro (organ) level, AFF might be affected by a redistribution of loads, mechanical properties, density, damage and mineral content. For example, the increased thickness of the femoral cortex observed in patients who have suffered an AFF (Larsen and Schmal, 2018) has been associated to an augmented strength and stiffness of the bone, which makes the skeletal structure more brittle at the same time (Donnelly et al., 2012). The consideration of these factors requires the use of FE models in combination with PK-PD models, as done by Hambli et al. (2016). In the present work, the risk of failure was only estimated at the local level through NDtF, but this limitation could be overcome by applying the present algorithm to such FE models.

4.2. Disease and Treatment 60Q6

The simulation of the non-treatment case predicted that the onset of the disease triggers the risk of local failure after a couple of years (see black line in Figure 5 top right). This occurs because damage is initially reduced by an increasing bone turnover rate (bone resorption is able to repair damage quite efficiently) and this decrease in damage compensates for the decrease in bone density, making NDtF to be approximately constant. In the mid-term, density has fallen so much that bone stiffness is significantly diminished. Strain increases and consequently damage rises, with bone resorption unable to repair the large amount of damage. The concurrence of damage increase and density decrease makes NDtF fall abruptly (higher risk of failure). In the long-term, though there is still an increasing risk of failure, NDtF reaches a stable value.

Anti-resorptive treatments lead to a decrease of bone turnover rate entailing an increase of mineral content and a subsequent increase of damage (see red line in Figure 5 bottom left), which is initially even higher than in the non-treatment group. However, in the long-term, bone density gain counteracts this effect by increasing stiffness and reducing strains, making damage slow down and, eventually, reducing the risk of failure.

4.3. Time Elapsed From the Onset of Disease to the Beginning of the Treatment

One factor that could limit the benefits of the treatment is the delay in its commencement. As stated before, non-treated patients could suffer a rapid bone loss right after the onset of disease, which would get stabilised after 10–12 years. However, damage did not stop rising in our simulations (though it did quite slowly in the long-term) and this made the risk of failure increase a little. This result highlights the importance of starting the treatment as soon as possible, before damage goes up excessively. Figure 6 analysed the effect of the time elapsed from the onset of disease to the beginning of the treatment. It could be seen how a late treatment can have severe consequences, as bone density gain is lower than in early treatments and, more importantly, because damage may rise even above the values reached in non-treated patients. This damage increase is quite remarkable in late treatments and is probably due to the deterioration of fatigue properties, the decrease of damage repairing and an insufficient increase of bone density and stiffness. NDtF shows how the combination of all these effects increased the risk of failure in late treatments. Although all of the NDtF curves seemed to converge to the same evolution, the period right after the start of the treatment was significantly more dangerous in late treatments.

4.4. Combination of Treatment and Exercise

The influence of exercise in combination with the denosumab treatment was analysed to conclude that, in general, an increment of the applied stress coincident with the commencing of the treatment is beneficial since it produces a greater bone density gain while keeping damage under control. The simulations showed that an increment of exercise is followed by a decrease of the risk of failure. On the contrary, a decrease of mechanical loading reduces the effectiveness of the treatment and this could seriously compromise bone integrity [see case 3 (dark blue line), in Figure 7]. In this case, bone formation is not promoted by the low stress and cannot compensate the bone density fall at the end of every treatment cycle, when bone turnover is slightly reactivated.

Case 5 (cyan in Figure 7) involved a strong increase of the stress in a late treatment and was analysed to illustrate three important ideas: (1) how dangerous late treatments can be; (2) the enhanced risk of failure that a significant increase of stress might produce and (3) how apparently small values of d ~ 0.05 are indeed high and can easily lead to failure in a few number of cycles. Late treatments were shown in Figure 6 to enhance damage accumulation and produce a limited bone density gain. This effect, in combination with a strong increment of stress, that also contributes to increase damage, led to a fatigue failure (d = 1 and NDtF = 0) soon after the beginning of the treatment.

It must be noted that all the cases analysed in Figure 7 implied overloads resulting in normal/moderate constant stresses. No high overloads or traumatic events were considered, which could cause significantly higher stresses, eventually leading to high-energy fractures. These high overloads were indeed considered in the calculation of NDtF, which provided the remaining fatigue life if the stress was increased up to 24 MPa, i.e., in case of an extra overload. In this regard, it is interesting to note that case 3 (underload by 10%) presented the highest risk of failure up to day 5,000 and yet case 5 was the only one to undergo failure. This is explained by the different performance of bone in both situations. In case 3, the risk of failure was determined by a low stiffness and seemed to be not so sensitive to overloads. Meanwhile, in case 5 it was determined by the high amount of unrepaired damage, which made bone more prone to failure in the event of a sudden overload, at least in the local scale, at the RVE. Certainly, the situation could be different at the organ level, as the reduction of stiffness of case 3 could redistribute the loads within the bone, what could initiate fracture elsewhere in the same bone.

The choice of 24 MPa as the overload to calculate NDtF was arbitrary, but allowed an easier comparison of the presented results and a more comprehensive analysis of the different factors that affect the risk of failure. The results were qualitatively similar for overloads around 24 MPa. For higher overloads, NDtF was very small regardless of the case, while lower overloads did not lead to a significant increase of the risk of failure in any case.

4.5. Importance of Brittleness vs. Unrepaired Damage in the Risk of Failure

Decrease in bone turnover due to anti-catabolic treatments has been associated to the development of AFF (Saita et al., 2015), due to the alteration of the tissue repair process (Mashiba et al., 2000). The accumulation of unrepaired microstructural damage results in unimpeded crack progression and may eventually lead to that type of fracture (Ettinger et al., 2013). However, another factor could also contribute to that increased occurrence of AFF, as hypothesised in a recent work (Martínez-Reina et al., 2021). This factor is bone mineral content, which increases as a consequence of the suppression of bone turnover, i.e., while mineral is not prevented from being accumulated within bone matrix and is not returned to blood serum by bone resorption. The mineral phase makes bone matrix more brittle, increasing the stiffness but reducing the fracture toughness of bone (Bala et al., 2013). This fact was confirmed by experimental studies that measured a higher amount of microstructural damage in interstitial bone, which has a higher BMC (Boyce et al., 1998; O'Brien et al., 2000; Qiu et al., 2005) and it was modelled in subsection 2.5 through the degradation of bone fatigue properties with the mineral content. In the above-mentioned work (Martínez-Reina et al., 2021), the risk of failure was not assessed and it was only associated to a high BMC. The present model was intended to evaluate the risk of failure by adding the damage level to the previous model as a new variable and considering damage accumulation by fatigue and damage repair by resorption. The new hypothesis is that it must be the concurrence of both factors (accumulation of unrepaired damage and increased brittleness) what would explain the high risk of AFF in antiresorptive treatments. The comparison made in Figure 6 aimed at discerning the relevance of those factors on the risk of AFF. From those simulations it could be concluded that both would play a role in the occurrence of AFF. If the stress is increased by 30%, the degradation of fatigue properties plays a key role in the occurrence of failure, as it only occurs if those properties are degraded. However, if the stress is increased by 40%, the accumulation of unrepaired damage is more important, as it is so fast that failure would occur regardless of whether the fatigue properties are degraded or not.

4.6. Limitations of the Study and Final Comments

The main limitation of this study is that the model was applied at the RVE level and the equations were solved only at this local level, without any interaction with the surroundings. Apart from diffusive terms of cell populations and biochemical factors, the distribution of stresses would play a key role in the behaviour at the organ level. Local changes in porosity, mineral content and damage would produce local changes in stiffness that might redistribute loads and affect the organ globally. To consider this, the model could be implemented in a Finite Element (FE) code. This would allow to evaluate that redistribution, which is particularly important in damage propagation and consequently in the assessment of the risk of fracture at the organ level. Moreover, the FE model would allow to simulate stress states more general than the simplistic uniaxial state modelled here. Nonetheless, this is out of the scope of the present paper and is left for future studies.

The risk of failure was also evaluated in a simplistic way: by estimating the remaining fatigue life of the RVE at a given time point in case of a constant overload and if all the variables except damage remained constant until failure. This estimated risk of failure should be considered only for the purpose of qualitative comparisons, as done here, and making a clear distinction between the risk of local failure, assessed at the RVE, and the risk of fracture, assessed at the organ level. As stated above, damage propagation and redistribution of stresses at the organ level would have a fundamental influence on the true risk of fracture of a given bone, but other factors such as microstructure or the multiaxiality of loads, not considered here, would also have it.

In our mineralisation algorithm, calcium (and phosphorus) availability is unlimited for deposition in bone matrix. This is only a simplification, as its availability depends on calcium (and phosphorus) homeostasis at the body level. Peterson and Riggs (2012) developed a model that accounted for calcium balance in blood serum. This process is controlled mainly by three mechanisms: (1) absorption in the intestine, (2) filtration/recirculation in the kidneys, and (3) deposition or retrieval from bone matrix through bone resorption. That balance is affected by numerous factors that could reduce mineral availability for bone deposition, among which there could be dietary restrictions and pathologies such as hypoparathyroidism or renal dysfunction. The incorporation of these mechanisms into the present model is currently under development.

Clinical results of bone density gain in denosumab treatments for PMO (Miller et al., 2008) exhibit a great variability that could be explained by the high number of factors that affect bone response. The present model is based on a previous one (Martínez-Reina et al., 2021) that was validated with the mentioned clinical results. Though not shown here, the algorithmic novelties incorporated into the present model did not produce significant changes in the model predictions and, particularly, the predicted bone density gain remained within the commented variability. Among the factors that may explain this variability we have investigated two of them: time elapsed from the onset of the disease and level of physical activity, though other factors such as treatment dose and frequency, incidence of the disease, patient's age and body weight, etc. could also have a great impact on the evolution of bone density gain and should be analysed in future studies.

Another aspect that this variability suggests is the need for a patient specific treatment, provided that this variability comes from patients' features that can be identified and measured. Obviously, this specific treatment should not conflict with WHO recommendations and therefore, they could not include dosage and frequency as variables, because these are fixed to 60Q6 for the treatment of PMO; however, other variables such as the activity level certainly could. For example, a gradual increase in activity level could be prescribed to old women who commence the treatment late. The parameters of the training regimen could be optimised as a function of patient's age, body weight, recent DEXA scans, etc. For this reason, it is of paramount importance to analyse the influence of more factors in order to design patient specific anti-resorptive treatments, which is the global aim of the present study.

Data Availability Statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

Author Contributions

JM-R and PP conceived, planned the computational simulations, and wrote the manuscript. JC-G and JM-R carried out the computational simulations. All authors discussed the results, reviewed, and commented on the manuscript.

Funding

Funding was provided by the European Regional Development Fund, the Consejería de Economía, Conocimiento, Empresas y Universidad de la Junta de Andalucía y la Universidad de Sevilla, dentro del Programa Operativo FEDER 2014-2020. It was awarded to the research project P18-RT-3611, entitled Efecto combinado del ejercicio físico y el denosumab en el tratamiento de la osteoporosis. Diseño de un tratamiento farmacológico específico de paciente for which this paper has been prepared.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

PP gratefully acknowledges support from the Australian Research Council (ARC ITTC, IC190100020).

Footnotes

1. ^ζ = 1 represents unrestricted access to denosumab, whereas ζ < 1 reflects restricted access, for example due to bone marrow being present or low blood perfusion. Hence, the denosumab concentration is bone site specific.

2. ^Typical values of variable d are around 0.01. This makes ${\Pi }_{act}^{dam}~10^{-3}$, while in Equation 7 ${\Pi }_{act}^{PTH}~10^{-2}$ (see Pivonka et al., 2008). With these values, the contributions of both terms in Equation (7) are similar given that the population of osteocytes (Ot) is typically one order of magnitude greater than that of osteoblasts precursor cells (Ob_p).

3. ^In the isotropic damage theory, Equation (14) can be similarly written in terms of the respective Young's moduli, E and E₀, as E = (1 − d)E₀ (Pattin et al., 1996; Zioupos, 1998). Subsequently, we will use the latter formulation.

4. ^In practical terms this condition was replaced by d(t + NDtF) = 0.999 to avoid the negative stiffness that could arise from d > 1 (see Equation 26).

5. ^We will identify the onset of the disease with the instant when the increase in RANKL production starts, though the disease is likely not evident at that point.

6. ^It must be recalled that NDtF is inversely related to the risk of failure: the fewer days are needed to reach d = 1, the higher is the risk of failure.

7. ^These changes of load must not be confused with the overload considered in the assessment of NDtF. For example, in case 3 the load does change from the homeostatic value to 90% of that value at the moment the treatment begins. Then, NDtF(σ, t) is calculated for the hypothetical case that the load were increased to σ.

8. ^The WHO-approved dose of 60 mg administered to a 60 kg patient would result in $Dos{e}_{den}=10^6$ ng of denosumab per kg of body weight.

References

Aguirre, J., Plotkin, L., Stewart, S., Weinstein, R., Parfitt, A., Manolagas, S., et al. (2006). Osteocyte apoptosis is induced by weightlessness in mice and precedes osteoclast recruitment and bone loss. J. Bone Miner. Res. 21, 605–615. doi: 10.1359/jbmr.060107

PubMed Abstract | CrossRef Full Text | Google Scholar

Aspenberg, P. (2014). Denosumab and atypical femoral fractures. Acta Orthop. 85:1. doi: 10.3109/17453674.2013.859423

PubMed Abstract | CrossRef Full Text | Google Scholar

Bala, Y., Farlay, D., and G., B. (2013). Bone mineralization: from tissue to crystal in normal and pathological contexts. Osteoporos. Int. 24, 2153–2166. doi: 10.1007/s00198-012-2228-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Beaupré, G., Orr, T., and Carter, D. (1990). An approach for time-dependent bone modeling and remodeling theoretical development. J. Biomech. 8, 651–661. doi: 10.1002/jor.1100080506

PubMed Abstract | CrossRef Full Text | Google Scholar

Boivin, G., and Meunier, P. J. (2002). The degree of mineralization of bone tissue measured by computerized quantitative contact microradiography. Calcif. Tissue Int. 70, 503–511. doi: 10.1007/s00223-001-2048-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Boyce, T., Fyhrie, D., Glotkowski, M., Radin, E., and Schaffler, M. (1998). Damage type and strain mode associations in human compact bone bending fatigue. J. Orthop. Res. 16, 322–329. doi: 10.1002/jor.1100160308

PubMed Abstract | CrossRef Full Text | Google Scholar

Burr, D., Martin, R., Schaffler, M., and Radin, E. (1985). Bone remodeling in response to in vivo fatigue microdamage. J. Biomech. 18, 189–200. doi: 10.1016/0021-9290(85)90204-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Cardoso, L., Herman, B., Verborgt, O., Laudier, D., Majeska, R., and Schaffler, M. (2009). Osteocyte apoptosis controls activation of intracortical resorption in response to bone fatigue. J. Bone Miner. Res. 24, 597–605. doi: 10.1359/jbmr.081210

PubMed Abstract | CrossRef Full Text | Google Scholar

Cummings, S. R., Martin, J. S., McClung, M. R., Siris, E. S., Eastell, R., and Reid, I. R. (2009). Denosumab for prevention of fractures in postmenopausal women with osteoporosis. N. Engl. J. Med. 361, 756–765. doi: 10.1056/NEJMoa0809493

PubMed Abstract | CrossRef Full Text | Google Scholar

Currey, J. (2004). Tensile yield in compact bone is determined by strain, post-yield behaviour by mineral content. J. Biomech. 37, 549–556. doi: 10.1016/j.jbiomech.2003.08.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Dempster, D., Brown, J., Fahrleitner-Pammer, A., Kendler, D., Rizzo, S., Valter, I., et al. (2018). Effects of long-term denosumab on bone histomorphometry and mineralisation in women with postmenopausal osteoporosis. J. Clin. Endocrinol. Metab. 103, 2498–2509. doi: 10.1210/jc.2017-02669

PubMed Abstract | CrossRef Full Text | Google Scholar

Dendorfer, S., Maier, H., and Hammer, J. (2006). Fatigue damage in cancellous bone: an experimental approach from continuum to micro scale. J. Mech. Behav. Biomed. Mater. 2, 113–119. doi: 10.1016/j.jmbbm.2008.03.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Donnelly, E., Meredith, D., Nguyen, J., Gladnick, B., Rebolledo, B., Shaffer, A. D., et al. (2012). Reduced cortical bone compositional heterogeneity with bisphosphonate treatment in postmenopausal women with intertrochanteric and subtrochanteric fractures. J. Bone Miner. Res. 27, 672–678. doi: 10.1002/jbmr.560

PubMed Abstract | CrossRef Full Text | Google Scholar

Ettinger, B., Burr, D., and Ritchie, R. (2013). Proposed pathogenesis for atypical femoral fractures: lessons from materials research. Bone 55, 495–500. doi: 10.1016/j.bone.2013.02.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, R., Gong, H., Zhang, X., Liu, J., Jia, Z., and Zhu, D. (2016). Modeling the mechanical consequences of age-related trabecular bone loss by XFEM simulation. Comput. Math. Method 2016:3495152. doi: 10.1155/2016/3495152

PubMed Abstract | CrossRef Full Text | Google Scholar

Frost, H. (2003). Bone's mechanostat: a 2003 update. Anat. Rec. Part A 275, 1081–1101. doi: 10.1002/ar.a.10119

CrossRef Full Text | Google Scholar

García-Aznar, J., Rueberg, T., and Doblaré, M. (2005). A bone remodelling model coupling microdamage growth and repair by 3D BMU-activity. Biomech. Model. Mechanobiol. 4, 147–167. doi: 10.1007/s10237-005-0067-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Hambli, R. (2013a). A quasi-brittle continuum damage finite element model of the human proximal femur based on element deletion. Med. Biol. Eng. Comput. 51, 219–231. doi: 10.1007/s11517-012-0986-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Hambli, R. (2013b). Micro-CT finite element model and experimental validation of trabecular bone damage and fracture. Bone 56, 363–374. doi: 10.1016/j.bone.2013.06.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Hambli, R., Boughattas, M., Daniel, J., and Kourta, A. (2016). Prediction of denosumab effects on bone remodeling: a combined pharmacokinetics and finite element modeling. J. Mech. Behav. Biomed. Mater. 60, 492–504. doi: 10.1016/j.jmbbm.2016.03.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Harrison, N., McDonell, P., Mullins, L., Wilson, N., O'Mahoney, D., and McHugh, P. (2013). Failure modelling of trabecular bone using a non-linear combined damage and fracture voxel finite element approach. Biomech. Model. Mechanbiol. 12, 225–241. doi: 10.1007/s10237-012-0394-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Hernandez, C., Beaupré, G., and Carter, D. (2001a). A model of mechanobiologic and metabolic influences on bone adaptation. J. Rehabil. Res. Dev. 37, 235–244.

PubMed Abstract | Google Scholar

Hernandez, C., Beaupré, G., Marcus, R., and Carter, D. (2001b). A theoretical analysis of the contributions of remodeling space, mineralisation, and bone balance to changes in bone mineral density during alendronate treatment. Bone 29, 511–516. doi: 10.1016/S8756-3282(01)00613-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Huiskes, R., Weinans, H., Grootenboer, H., Dalstra, M., Fudala, B., and Sloof, T. (1987). Adaptive bone-remodelling theory applied to prosthetic design analysis. J. Biomech. 20, 1135–1150. doi: 10.1016/0021-9290(87)90030-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Jacobs, C. (1994). Numerical Simulation of BoneAdaptation to Mechanical Loading (Ph.D. thesis). Stanford, CA: Stanford University.

Google Scholar

Jin, A., Cobb, J., Hansen, U., Bhattacharya, R., Reinhard, C., Vo, N., et al. (2017). The effect of long-term bisphosphonate therapy on trabecular bone strength and microcrack density. Bone Joint Res. 6, 602–609. doi: 10.1302/2046-3758.610.BJR-2016-0321.R1

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, Z., and El-Deiry, W. (2005). Overview of cell death signaling pathways. Cancer Biol. Ther. 4, 139–163. doi: 10.4161/cbt.4.2.1508

PubMed Abstract | CrossRef Full Text | Google Scholar

Juvinall, R. (1967). Engineering Considerations of Stress, Strain and Strength. New York, NY: McGraw-Hill.

Google Scholar

Kurata, K., Heino, T., Higaki, H., and Vaananen, H. (2006). Bone marrow cell differentiation induced by mechanically damaged osteocytes in 3D gel-embedded culture. J. Bone Miner. Res. 21, 616–625. doi: 10.1359/jbmr.060106

PubMed Abstract | CrossRef Full Text | Google Scholar

Langdahl, B. (2020). Treatment of postmenopausal osteoporosis with bone-forming and antiresorptive treatments: combined and sequential approaches. Bone 139:115516. doi: 10.1016/j.bone.2020.115516

PubMed Abstract | CrossRef Full Text | Google Scholar

Larsen, M., and Schmal, H. (2018). The enigma of atypical femoral fractures. A summary of current knowledge. EFORT Open Rev. 3, 494–500. doi: 10.1302/2058-5241.3.170070

PubMed Abstract | CrossRef Full Text | Google Scholar

Lemaire, V., Tobin, F. L., Greller, L. D., Cho, C. R., and Suva, L. J. (2004). Modeling the interactions between osteoblast and osteoclast activities in bone remodeling. J. Theor. Biol. 229, 293–309. doi: 10.1016/j.jtbi.2004.03.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Lemaitre, J., and Chaboche, J. (1990). Mechanics of Solid Materials. Cambridge: Cambridge University Press.

Google Scholar

Marathe, A., Peterson, M. C., and Mager, D. E. (2008). Integrated cellular bone homeostasis model for denosumab pharmacodynamics in multiple myeloma patients. J. Pharmacol. Exp. Ther. 326, 555–562. doi: 10.1124/jpet.108.137703

PubMed Abstract | CrossRef Full Text | Google Scholar

Marathe, D. D., Marathe, A., and Mager, D. E. (2011). Integrated model for denosumab and ibandronate pharmacodynamics in postmenopausal women. Biopharm. Drug Dispos. 32, 471–481. doi: 10.1002/bdd.770

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, M., Sansalone, V., Cooper, D., Forwood, M., and Pivonka, P. (2019). Mechanobiological osteocyte feedback drives mechanostat regulation of bone in a multiscale computational model. Biomech. Model. Mechanobiol. 18, 1475–1496. doi: 10.1007/s10237-019-01158-w

PubMed Abstract | CrossRef Full Text | Google Scholar

Martin, R. (1984). Porosity and specific surface of bone. Crit. Rev. Biomed. Engl. 10, 179–222.

Google Scholar

Martínez-Reina, J., Calvo-Gallego, J., and Pivonka, P. (2021). Are drug holidays a safe option in treatment of osteoporosis? – Insights from an in silico mechanistic PK–PD model of denosumab treatment of postmenopausal osteoporosis. J. Mech. Behav. Biomed. Mater. 113:104140. doi: 10.1016/j.jmbbm.2020.104140

PubMed Abstract | CrossRef Full Text | Google Scholar

Martínez-Reina, J., García-Aznar, J., Domínguez, J., and Doblaré, M. (2008). On the role of bone damage in calcium homeostasis. J. Ther. Biol. 254, 704–712. doi: 10.1016/j.jtbi.2008.06.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Martínez-Reina, J., Garía-Aznar, J., Domínguez, J., and Doblaré, M. (2009). A bone remodelling model including the directional activity of BMUs. Biomech. Model. Mechanobiol. 8, 111–127. doi: 10.1007/s10237-008-0122-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Martínez-Reina, J., and Pivonka, P. (2019). Effects of long-term treatment of denosumab on bone mineral density: insights from an in-silico model of bone mineralization. Bone 125, 87–95. doi: 10.1016/j.bone.2019.04.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Mashiba, T., Hirano, T., Turner, C., Forwood, M., Johnston, C., and Burr, D. (2000). Suppressed bone turnover by bisphosphonates increases microdamage accumulation and reduces some biomechanical properties in dog rib. J. Bone. Miner. Res. 15, 613–620. doi: 10.1359/jbmr.2000.15.4.613

PubMed Abstract | CrossRef Full Text | Google Scholar

Miller, P. D., Bolognese, M. A., Lewiecki, E. M., McClung, M. R., Ding, B., Austin, M., et al. (2008). Effect of denosumab on bone density and turnover in postmenopausal women with low bone mass after long-term continued, discontinued, and restarting of therapy: a randomized blinded phase 2 clinical trial. Bone 43, 222–229. doi: 10.1016/j.bone.2008.04.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakashima, T., Hayashi, M., Fukunaga, T., Kurata, K., Oh-hora, M., Feng, J., et al. (2011). Evidence for osteocyte regulation of bone homeostasis through RANKL expression. Nat. Med. 17, 1231–1234. doi: 10.1038/nm.2452

PubMed Abstract | CrossRef Full Text | Google Scholar

Nyman, J., Yeh, O., Hazelwood, S., and Martin, R. (2004). A theoretical analysis of long-term bisphosphonate effects on trabecular bone volume and microdamage. Bone 35, 296–305. doi: 10.1016/j.bone.2004.03.007

PubMed Abstract | CrossRef Full Text | Google Scholar

O'Brien, F., Taylor, D., Dickson, G., and Lee, T. (2000). Visualisation of three-dimensional microcracks in compact bone. J. Anat. 197, 413–420. doi: 10.1046/j.1469-7580.2000.19730413.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Parfitt, A. (2002). Targeted and nontargeted bone remodeling: relationship to basic multicellular unit origination and progression. Bone 30, 5–7. doi: 10.1016/S8756-3282(01)00642-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Pattin, C., Caler, W., and Carter, D. (1996). Cyclic mechanical property degradation during fatigue loading of cortical bone. J. Biomech. 29, 69–79. doi: 10.1016/0021-9290(94)00156-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Peterson, M. C., and Riggs, M. M. (2010). A physiologically based mathematical model of integrated calcium homeostasis and bone remodeling. Bone 46, 49–63. doi: 10.1016/j.bone.2009.08.053

PubMed Abstract | CrossRef Full Text | Google Scholar

Peterson, M. C., and Riggs, M. M. (2012). Predicting nonlinear changes in bone mineral density over time using a multiscale systems pharmacology model. CPT Pharmacometr. Syst. Pharmacol. 1, 1–8. doi: 10.1038/psp.2012.15

PubMed Abstract | CrossRef Full Text | Google Scholar

Pivonka, P., Buenzli, P., Scheiner, S., Hellmich, C., and Dunstan, C. (2012). The influence of bone surface availability in bone remodellin a mathematical model including coupled geometrical and biomechanical regulations of bone cells. Eng. Struct. 47, 134–147. doi: 10.1016/j.engstruct.2012.09.006

CrossRef Full Text | Google Scholar

Pivonka, P., Zimak, J., Smith, D. W., Gardiner, B. S., Dunstan, C. R., Sims, N. A., et al. (2008). Model structure and control of bone remodeling: a theoretical study. Bone 43, 249–263. doi: 10.1016/j.bone.2008.03.025

PubMed Abstract | CrossRef Full Text | Google Scholar

Pivonka, P., Zimak, J., Smith, D. W., Gardiner, B. S., Dunstan, C. R., Sims, N. A., et al. (2010). Theoretical investigation of the role of the rank rankl opg system in bone remodeling. J. Theor. Biol. 262, 306–316. doi: 10.1016/j.jtbi.2009.09.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Qiu, S., Rao, D., Fyhrie, D., Palnitkar, S., and Parfitt, A. (2005). The morphological association between microcracks and osteocyte lacunae in human cortical bone. Bone 37, 10–15. doi: 10.1016/j.bone.2005.01.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Rapillard, L., Charlebois, M., and Zysset, P. (2006). Compressive fatigue behavior of human vertebral trabecular bone. J. Biomech. 39, 2133–2139. doi: 10.1016/j.jbiomech.2005.04.033

PubMed Abstract | CrossRef Full Text | Google Scholar

Riggs, B. L., and Parfitt, A. M. (2005). Drugs used to treat osteoporosis: the critical need for a uniform nomenclature based on their action on bone remodeling. J. Bone Miner. Res. 20, 177–184. doi: 10.1359/JBMR.041114

PubMed Abstract | CrossRef Full Text | Google Scholar

Saita, Y., Ishijima, M., and Kaneko, K. (2015). Atypical femoral fractures and bisphosphonate use: current evidence and clinical implications. Ther. Adv. Chron. Dis. 6, 185–193. doi: 10.1177/2040622315584114

PubMed Abstract | CrossRef Full Text | Google Scholar

Scheiner, S., Pivonka, P., and Hellmich, C. (2013). Coupling systems biology with multiscale mechanics, for computer simulations of bone remodeling. Comput. Methods Appl. Mech. Eng. 254, 181–196. doi: 10.1016/j.cma.2012.10.015

CrossRef Full Text | Google Scholar

Scheiner, S., Pivonka, P., Smith, D. W., Dunstan, C. R., and Hellmich, C. (2014). Mathematical modeling of postmenopausal osteoporosis and its treatment by the anti-catabolic drug denosumab. Int. J. Numer. Methods Biomed. Eng. 30, 1–27. doi: 10.1002/cnm.2584

PubMed Abstract | CrossRef Full Text | Google Scholar

Sims, N., and Martin, T. (2015). Coupling signals between the osteoclast and osteoblast: how are messages transmitted between these temporary visitors to the bone surface? Front. Endocrinol. 6:41. doi: 10.3389/fendo.2015.00041

PubMed Abstract | CrossRef Full Text | Google Scholar

Tatsumi, S., Ishii, K., Amizuka, N., Li, M., Kobayashi, T., Kohno, K., et al. (2007). Targeted ablation of osteocytes induces osteoporosis with defective mechanotransduction. Cell Metab. 5, 464–475. doi: 10.1016/j.cmet.2007.05.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Turner, C. (1999). Toward a mathematical description of bone biology: the principle of cellular accommodation. Calcif. Tissue Int. 65, 466–471. doi: 10.1007/s002239900734

PubMed Abstract | CrossRef Full Text | Google Scholar

Vashishth, D., Gibson, G., Khoury, J., Schaffler, M., Kimura, J., and Fyhrie, D. (2001). Influence of nonenzymatic glycation on biomechanical properties of cortical bone. Bone 28, 195–201. doi: 10.1016/S8756-3282(00)00434-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Verborgt, O., Gibson, G., and Schaffler, M. (2000). Loss of osteocyte integrity in association with microdamage and bone remodeling after fatigue in vivo. J. Bone Miner. Res. 15, 60–67. doi: 10.1359/jbmr.2000.15.1.60

PubMed Abstract | CrossRef Full Text | Google Scholar

Verborgt, O., Tatton, N., Majeska, R., and Schaffler, M. (2002). Spatial distribution of Bax and Bcl-2 in osteocytes after bone fatigue: complementary roles in bone remodeling regulation? J. Bone Miner. Res. 17, 907–914. doi: 10.1359/jbmr.2002.17.5.907

PubMed Abstract | CrossRef Full Text | Google Scholar

Xiong, J., and O'Brien, C. (2012). Osteocyte rankl: new insights into the control of bone eemodeling. J. Bone Miner. Res. 27, 499–505. doi: 10.1002/jbmr.1547

PubMed Abstract | CrossRef Full Text | Google Scholar

Zioupos, P., and znd Currey, J. (1998). Changes in the stiffness, strength, and toughness of human cortical bone with age. Bone 22, 57–66. doi: 10.1016/S8756-3282(97)00228-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Appendix

Algorithm of Bone Mineralisation

Bone tissue is made up of bone matrix and pores. Thus, the representative volume element, V_RVE, can be divided into the bone matrix volume, V_bm, and the volume of vascular pores, V_vas. In turn, bone matrix volume is divided into inorganic (mineral), organic (mainly collagen), and water phases, respectively designated as V_m, V_o, and V_w:

$\begin{array}{ll}{V}_{RVE}=V_{bm}+V_{vas}=V_m+V_o+V_w+V_{vas}& (29)\end{array}$

The volume fractions of extravascular bone matrix and vascular pores are respectively given by f_bm = V_bm/V_RVE and f_vas = 1 − f_bm = V_vas/V_RVE. The mineral content is usually measured by the so-called ash fraction, the ratio between mass of mineral m_m (or ash mass) and dry mass (the sum of inorganic and organic mass):

$\begin{array}{ll}\alpha =\frac{m_m}{m_m+m_o}=\frac{{\rho }_m{V}_m}{{\rho }_m{V}_m+{\rho }_o{V}_o}& (30)\end{array}$

where density of hydroxyapatite is taken for ${\rho }_m=3.2g/c{m}^3$ (Currey, 2004). Organic phase is mainly composed of type I collagen, but other non-collagenous proteins are also present. Thus, ${\rho }_o=1.2g/c{m}^3$ was adjusted to provide a tissue density ${\rho }_t=2.1g/c{m}^3$ for the completely mineralised tissue α = 0.73 (García-Aznar et al., 2005).

Specific volumes are defined by: v_o = V_o/V_bm, v_m = V_m/V_bm, and v_w = V_w/V_bm, which implies that v_o + v_w + v_m = 1. Then, Equation (30) can be given in terms of the specific volumes:

$\begin{array}{ll}\alpha =\frac{{\rho }_m{v}_m}{{\rho }_m{v}_m+{\rho }_o{v}_o}& (31)\end{array}$

Bone tissue density is given by:

$\begin{array}{ll}{\rho }_t=\frac{m}{V_{bm}}={\rho }_m{v}_m+{\rho }_o{v}_o+{\rho }_w{v}_w& (32)\end{array}$

Mineral accumulates by displacing water present in bone matrix (Hernandez et al., 2001a). Therefore, the volume ratio of organic phase is assumed constant during the mineralisation process, v_o = 3/7 (Martin, 1984); while the variations of mineral and water volume ratios hold Δv_m = −Δv_w. We have followed Hernandez et al. (2001a) to assess the increase of v_m with time, by distinguishing the mineralisation lag time; the primary phase, with a linear increase, and the secondary phase, with an exponentially decreasing rate:

$\begin{array}{ll}{v}_m(t)=\{\begin{array}{ll}0& \text{if\hspace{1em}}t\le t_{mlt}\\ v_{mprim}\frac{t-t_{mlt}}{t_{prim}}& \text{if\hspace{1em}}{t}_{mlt}<t\le t_{prim}+t_{mlt}\\ v_{mmax}-(v_{mmax}-v_{mprim})\\ e^{-\kappa ·(t-t_{prim}-t_{mlt})}& \text{if\hspace{1em}}{t}_{prim}+t_{mlt}<t\end{array}& (33)\end{array}$

where t_mlt and t_prim are, respectively, the length of the mineralisation lag time and the primary phase; v_{m prim} is the mineral specific volume at the end of the primary phase, corresponding to α = 0.45 (Hernandez et al., 2001a); v_{m max} is the mineral specific volume corresponding to the maximum calcium content, 300 mg/g (Currey, 2004); and κ is a parameter measuring the rate of mineral deposition during the secondary phase (see Table A1). Note that, with the assumptions made above, Equations (31) and (32) establish a biunivocal relation between α and ρ_t. Also, according to Equation (25), apparent density is univocally determined by f_bm and α.

TABLE A1

Table A1. Values taken for the constants of the PK-PD model.

The amount of mineral contained in a RVE depends on the age of the tissue through Equation (33), but the RVE can be made up of tissue patches formed in the recent history, viz. of different ages. Moreover, the tissue within the RVE can be resorbed, which puts the mineral back into the blood flow. The amounts of tissue of different ages contained in the RVE are estimated using the algorithm depicted in Figure 9 (Martínez-Reina et al., 2009). ${\overline{V}}_{form}(t,\tau )$ provides the bone volume formed τ days ago and still present (not yet resorbed) at time t. Knowing the distribution of tissue patches of different ages at day t (left column) and the volume formed (V_form(t) = k_formOb_a(t)) and resorbed that day (V_res(t) = k_resOc_a(t)), the distribution at day t + 1 (right column) can be estimated:

$\begin{array}{ll}{\overline{V}}_{form}(t+1,i+1)={\overline{V}}_{form}(t,i)-V_{res}(t)\frac{{\overline{V}}_{form}(t,i)}{V_{bm}(t)}& (34)\end{array}$

Finally, the mineral content of each patch is summed to estimate the average mineral content of the RVE at day t + 1.

$\begin{array}{ll}{v}_m(t+1)=\frac{{\displaystyle {\sum }_{i=0}^{t_R}}{\overline{V}}_{form}(t+1,i)·v_m(i)}{V_{bm}(t+1)}& (35)\end{array}$

where the mineral contents of the patches, v_m(i), are calculated through Equation (33). t_R represents the residence time, e.g., the typical time the patch tissue remains within the bone before being resorbed. This residence time can be very large but the queue can be truncated at a shorter time to reduce the computational cost. See (Martínez-Reina and Pivonka, 2019) for more details.

FIGURE 9

Figure 9. FIFO (first in - first out) queue algorithm used to update the distribution of tissue patches of different ages within the RVE.

It must be noted here that the recursive character of this mineralisation algorithm makes it necessary to integrate the set of differential equations that governs the model using an explicit integration scheme.

Adjustment of Mechanoregulation

The points that delimit the piecewise linear function ${\Pi }_{act}^{{\psi }_{bm}}$ in Figure 3 were adjusted to accomplish the Mechanostat Theory and the Principle of Cellular Accommodation. The Mechanostat Theory establishes and “adapted window” between 800 and 1,200 με (Frost, 2003), as were the homeostatic strains obtained in Figure 4.

The Principle of Cellular Accommodation establishes that bone cells react strongly to variations in their environment, but eventually “accommodate” to steady state signals (Turner, 1999), in such a way that those variations in the stress would produce a transient response to adapt to the new environment. To check if the proposed model accomplished with that principle, a set of simulations were run starting from the homeostatic situation and implementing variations of the stress that were kept constant to let bone reach a new steady state (see Figure 10). Increments of stress lead to an increase in bone density and stiffness (reversely in the case of decrements of stress), while strains lie out of the “adapted window” defined by the Mechanostat Theory. Bone tends to return to that strain range and accommodates to the new environment when the “adapted window” is reached. The limits established in Figure 3 to define ${\Pi }_{act}^{{\psi }_{bm}}$ were adjusted to achieve the behaviour of Figure 10, with an equilibrium strain around 1,200 με, as stated before, and a slope for the adaptation response small enough not to produce numerical instabilities and large enough to achieve a new equilibrium in a reasonable time (70% of the total change reached after 8–10 years).

FIGURE 10

Figure 10. Variation of density and strain under uniaxial loading simulating underuse and overuse states.

One Compartment PK Model of Denosumab

Among the pharmacokinetic (PK) models of denosumab that can be found in the literature, we have followed the one proposed by Marathe et al. (2011). This is a one-compartment model with Michaelis-Menten kinetics that estimates serum denosumab profiles. This model includes a first-order rate process (constant k_a) governing the absorption of the drug (input parameter Dose_den) from the subcutaneous injection site into the central compartment (variable C_{P, den} and constant V_c). Drug elimination from the central compartment is described by a combination of a linear first-order process (constant k_el) and a non-linear saturation process (constants V_max, K_m):

$\begin{array}{ll}\frac{d{C}_{p,den}}{dt}=k_a\frac{Dos{e}_{den}}{V_c/F}·e^{-k_{a}t}-\frac{C_{p,den}}{K_m+C_{p,den}}\frac{V_{max}}{V_c/F}-k_{el}·C_{p,den}& (36)\end{array}$

where, V_c is the volume of the central compartment and the factor F is the bioavailability, which is equal to 1 when the drug is administered intravenously, as assumed here. In Equation (36) Dose_den is given in ng of denosumab per kg of body weight⁸. C_{p, den} is the concentration of denosumab in the central compartment, calculated (in ng/ml) as a function of time by solving the differential equation (36) and subsequently converted into pmol/l, through the molecular weight of denosumab M_den = 149 kDa (Amgen). The initial condition for Equation (36) is set to zero, indicating the absence of drug. The prolonged absorption phase and the absence of intravenous data precludes the need for including distribution of the drug to a non-specific tissue compartment and thus reduces the number of parameters in this model.

Biochemical Regulatory Factors

TGF–β is stored in the bone matrix and released during resorption by osteoclasts. Its concentration is calculated following Pivonka et al. (2008):

$\begin{array}{ll}TGF-\beta =\frac{{\alpha }_{TGF-\beta }{k}_{res}O{c}_a}{{\tilde{D}}_{TGF-\beta }}& (37)\end{array}$

where α_TGF−β is the concentration of TGF–β in bone matrix and ${\tilde{D}}_{TGF-\beta }$ is the TGF–β degradation rate. The concentration of TGF–β is used to define the activator/repressor functions in Equations (1)–(3):

$\begin{array}{ll}{\Pi }_{act}^{TGF-\beta }=\frac{TGF-\beta }{K_{act}^{TGF-\beta }+TGF-\beta }& (38)\end{array}$

$\begin{array}{ll}{\Pi }_{rep}^{TGF-\beta }=\frac{K_{rep}^{TGF-\beta }}{K_{rep}^{TGF-\beta }+TGF-\beta }& (39)\end{array}$

with $K_{act}^{TGF-\beta }$ and $K_{rep}^{TGF-\beta }$ the activation and repression constants, respectively. RANK is expressed by osteoclasts precursors:

$\begin{array}{ll}RANK=R_{RANK}·O{c}_p& (40)\end{array}$

where R_RANK is the carrying capacity. Concentration of OPG is downregulated by PTH and is calculated following Pivonka et al. (2012):

$\begin{array}{ll}OPG=\frac{{\beta }_{OPG}O{b}_a{\Pi }_{rep}^{PTH}}{{\tilde{D}}_{OPG}+\frac{{\beta }_{OPG}O{b}_a{\Pi }_{rep}^{PTH}}{OP{G}_{max}}}& (41)\end{array}$

where β_OPG and ${\tilde{D}}_{OPG}$ are, respectively, the production and degradation rates of OPG and the activator/repressor functions that govern PTH regulation on the RANKL-RANK-OPG signaling pathway (Equations 41, 7) are given by:

$\begin{array}{ll}{\Pi }_{act}^{PTH}=\frac{PTH}{K_{act}^{PTH}+PTH}& (42)\end{array}$

$\begin{array}{ll}{\Pi }_{rep}^{PTH}=\frac{K_{rep}^{PTH}}{K_{rep}^{PTH}+PTH}& (43)\end{array}$

where the concentration PTH = 2.91 pM has been assumed constant in this case (Pivonka et al., 2008) and $K_{act}^{PTH}$ and $K_{rep}^{PTH}$ are the corresponding activation and repression constants.

The values of the constants of the PK-PD model are given in Table A1. A detailed discussion of these values can be consulted in (Martínez-Reina and Pivonka, 2019; Marathe et al., 2011; Scheiner et al., 2014).