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.

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.

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 fatigueloading 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 antiresorptive 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., 2008Pivonka et al., , 2010Scheiner et al., 2013). Mechanical overuse is simulated via increase of osteoblast precursors proliferation (Pivonka et al., 2012;Scheiner et al., 2013Scheiner et al., , 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.

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., 2008Pivonka et al., , 2010Scheiner 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: (1) where D Ob u , D Ob p , and D Oc p 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 Ob p gives the maximum proliferation rate. A Ob a and A Oc a are the apoptosis rates of Ob a and Oc a respectively. The variables TGF−β act , and TGF−β rep represent activator and repressor functions related to the binding of TGF-β to its receptor. Similarly, RANKL act is the activator function related to the RANK-RANKL binding. ψ bm act 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: 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.

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(Marathe et al., , 2011Scheiner et al., 2014). In these models, action of denosumab is taken into account via the RANKL activator function RANKL act (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 RANKL act . Adapting the approach of Scheiner et al. (2014), 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 1 . 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.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: 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(Pivonka et al., , 2010, and we have introduced in the present model an upregulation factor of RANKL expression by osteocytes due to damage, through a sigmoidal function: 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, δ 50 = 0.1, were chosen in such a way that both terms in Equation 7 are typically of the same order of magnitude 2 . Verborgt et al. (2000) showed 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 dam act ∼ 10 −3 , while in Equation 7 PTH act ∼ 10 −2 (see Pivonka et al., 2008). With these values, the Frontiers in Bioengineering and Biotechnology | www.frontiersin.org 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 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 ). 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 mech RANKL , and a disease-related increase in RANKL production over time, P PMO RANKL : The first term is explained in section 2.3, while the term P PMO RANKL is a consequence of the onset of menopause though increasing gradually over time, through the following sigmoidal function: where P PMO,max RANKL 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 PMO,max RANKL is reached.
Finally, the activator function of RANKL in Equation (3) can be expressed as: The concentrations of denosumab, RANK, OPG, PTH, and TGF−β along with the model parameters of the binding reactions are provided in the Appendix.

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)  act (see Figure 3), but indirectly through the accumulation of microstructural damage.
Mechanical disuse is assumed to enhance the production of RANKL on osteoblasts precursors, through the term P mech RANKL , which is modulated by the strain energy density (SED) of bone matrix, designated as ψ bm : where ψ r is the SED below which underuse increases RANKL production and P mech,max RANKL 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 ψ bm act , 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 mech RANKL ). Obviously, the steepest piece of function ψ bm act 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: being E the Young's modulus. The parameter ψ r in Equation (12) as well as ψ bm 1 and ψ bm 2 , used in Figure 3 to define function ψ bm act , were defined using Equation (13), respectively with MES r = 1,000 µε, MES 1 = 800 µε, and MES 2 = 1,600 µε. These values and the values of ψ bm act corresponding to ψ bm 1 and ψ bm 2 were adjusted to reproduce the Mechanostat Theory along with the Principle of Cellular Accomodation (see Appendix).

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: where C and C 0 are, respectively, the stiffness tensors of damaged and undamaged bone (Lemaitre and Chaboche, 1990) 3 . 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:ḋ whereḋ A is the rate of damage accumulation by fatigue loading andḋ 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:ḋ 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: 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 0 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: where N is the number of cycles and C t2 = −20 was fitted from the experimental curves along with: 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 .
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: This procedure allows working out the increment of damage, d A , but requires that Equation (15)

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 7 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: Here, the most typical value β = 2 was chosen following  (2004): 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), , 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.

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

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. 4 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: 5. With d(t + τ ) and ε max , use Equation (20) to work outÑ. 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).
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.

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.

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. 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 PMO,max RANKL controls the degree of incidence of the disease, starting from 600 pM day −1 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 5 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 6 . 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.
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 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, 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. 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 7 . 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.
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 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 σ . 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. 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.

General Comments
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.

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 midterm, 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.

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

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

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.

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.

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

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 : 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): where density of hydroxyapatite is taken for ρ m = 3.2 g/cm 3 (Currey, 2004). Organic phase is mainly composed of type I collagen, but other non-collagenous proteins are also present. Thus, ρ o = 1.2 g/cm 3 was adjusted to provide a tissue density ρ t = 2.1 g/cm 3 for the completely mineralised tissue α = 0.73 (García-Aznar et al., 2005). Specific volumes are defined by: 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: Bone tissue density is given by: 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: 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 α.
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). V form (t, τ ) 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 form Ob a (t)) and resorbed that day (V res (t) = k res Oc a (t)), the distribution at day t +1 (right column) can be estimated: Finally, the mineral content of each patch is summed to estimate the average mineral content of the RVE at day t + 1.
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.
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 ψ bm act 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 ψ bm act 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).

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 ): dC p,den dt = k a Dose den V c /F · e −k a t − C p,den K m + C p,den V max V c /F − k el · C p,den 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 8 . 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 8 The WHO-approved dose of 60 mg administered to a 60 kg patient would result in Dose den = 10 6 ng of denosumab per kg of body weight.
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): where α TGF−β is the concentration of TGF-β in bone matrix andD TGF−β is the TGF-β degradation rate. The concentration of TGF-β is used to define the activator/repressor functions in Equations (1) the activation and repression constants, respectively. RANK is expressed by osteoclasts precursors: where R RANK is the carrying capacity. Concentration of OPG is downregulated by PTH and is calculated following Pivonka et al. where the concentration PTH = 2.91 pM has been assumed constant in this case (Pivonka et al., 2008) and K PTH act and K PTH rep 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).