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

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 Javier Mart´ınez-Reinaa,1,∗, Jos´e Luis Calvo-Gallegoa, Peter Pivonkab,2 aDepartmento de Ingenier´ıa Mec´anica y Fabricaci´on, Universidad de Sevilla, Seville 41092, Spain bSchool of Mechanical, Medical and Process Engineering, Queensland University of Technology, QLD 4000, Australia. Abstract Recent reviews by the clinical bone research community suggest caution with prescription of drug holidays for patients with postmenopausal osteoporosis (PMO) treated with denosumab for an extended period of time. Main reasons for this suggestion are based on the fact that discontinuation of denosumab treatment leads to a relapse of osteoclastic bone resorption and a loss of bone mineral density (BMD) to pre-treatment levels at only 12-28 months. The question remains what is the best treatment option for cases where it is required to discontinue and/or reduce the drug dose and what are the consequences on BMD and bone turnover markers (BTMs). The latter questions are difficult to be addressed using clinical trials alone given the large number of parameter combinations involved to answer this problem. In this paper, we apply a recently developed in silico mechanistic pharmacokinetic - pharmacodymamic (PK-PD) model of the effect of denosumab on bone remodelling in PMO. To address the above clinical relevant questions, we design a wide range of current and virtual treatment regimens to study ∗Corresponding author 1Departmento de Ingenier´ıa Mec´anica y Fabricaci´on, Universidad de Sevilla. Escuela T´ecnica Superior de Ingenier´ıa. Camino de los Descubrimientos s/n, 41092. Sevilla. [email protected] 2School of Chemistry, Physics and Mechanical Engineering, Queensland University of Technology, 2 George Street, Brisbane QLD 4000, Australia Accepted version in Journal of L A T EX Templates the effect of drug holiday duration and therapy resumption on the evolution of BTMs, BMD and mineral content. Our numerical simulation results indicate the symptomatic effect of denosumab, which is lost once treatment is stopped. This effect is most clearly seen on rapid loss of BMD to pretreatment levels 12 months after the last injection (8% and 3.6% per year in the lumbar spine and femoral neck, respectively). Also, we identify that independently of the duration of drug holiday (i.e. 12, 16 or 18 months) resuming treatment can restore BMD quite effectively. However, the latter result does not consider the possibility of potential fractures that can occur during the drug holiday. Finally, we identify a treatment case most promising for achieving maintenance of BMD and mineral content, while moderately increasing BTMs. The latter case uses no drug holiday, but reduces the most commonly prescribed denosumab dose (60 mg every 6 months) by half at same interval. Keywords: Postmenopausal osteoporosis; denosumab; RANK-RANKL-OPG pathway; PK-PD modelling; bone fragility; drug holidays 1. Introduction A recent review by Tsourdi et al. as part of development of a mission statement for The European Calcified Tissue Society (ECTS) indicates that optimal duration of osteoporosis (OP) treatment remains controversial [1]. With the large number of currently available drugs for osteoporosis treat-5 ment it is essential to understand how long these drugs can be used safely, with typical treatment periods ranging from 5 years up to 20 years [2]. Denosumab and bisphosphonates are two of those drugs. Denosumab is a fully monoclonal antibody that binds to RANKL with high affinity and specificity, inhibiting RANKL-RANK interaction [3]. Bisphosphonates constitute10 a group of drugs with a very high affinity for bone mineral to which they bind [4] and a mechanism of action consisting in the inhibition of osteoclastic activity through cytoskeletal disruption, changing intracellular protein traffic, blocking intracellular signal transduction pathways and induction of osteoclast apoptosis [5]. Discontinuation of bisphosphonates, one of the most15 commonly prescribed OP drug, is considered safe due to the fact that the antiresorptive effect after stopping treatment persists for months or years based on the fact that bisphosphonates are absorbed onto and into the bone 2 matrix with high affinity. However, the concept of treatment break does not apply to drugs other than bisphosphonates, since for drugs without skeletal20 retention the fracture risk is expected to increase after drug discontinuation. Despite the fact that the concept of a treatment break does not apply to drugs without skeletal retention, it was reported that many clinicians take patients off the drug after about 5 years of denosumab treatment (see [1] and references therein). Also, patients are often advised by dentists that25 medication should be stopped temporarily before a dental procedure to avoid the risk of osteonecrosis of the jaw (ONJ). In addition, physicians themselves may recommend that treatment be stopped after several years in patients whose BMD has increased sufficiently to move the patient above the threshold of osteoporosis and out of a high-risk category based on T−score.30 Recently, there has been concern that discontinuation of denosumab will lead to an increased risk of multiple vertebral fractures associated with rapid bone loss when treatment is stopped [6,7]. ECTS formed a working group to address this problem and perform a systematic review of existing literature on the effects of stopping denosumab and provide advice on management35 [1]. Data from phase 2 and 3 clinical trials underscore a rapid decrease of bone mineral density (BMD) and a steep increase in bone turnover markers (BTMs) after discontinuation of denosumab [8,9]. Clinical case series report multiple vertebral fractures after discontinuation of denosumab [6] and a renewed analysis of FREEDOM and FREEDOM extension Trial suggests,40 albeit does not prove, that the risk of multiple vertebral fractures may be increased when denosumab is stopped due to a rebound increase in bone resorption. In his editorial, McClung highlights the rebound in remodelling rates to values higher than pre-treatment levels, with rapid bone loss upon stopping45 denosumab which is well documented (see [10] and references therein). In the phase 2 study of Miller et al., follow-up after discontinuing denosumab treatment was available in 50 patients [9]. After two years of treatment with 210 mg denosumab every 6 months or 30 mg every 3 months, serum CTX levels increased to twice the placebo value at 6-12 months after stopping50 therapy and fell back to baseline two years after treatment was discontinued. The mechanism of the overshoot in bone resorption is not understood, but it is believed to be due to an expanded pool of osteoclast precursors that are simultaneously activated or due to a high RANK ligand/OPG ratio after denosumab is cleared from the circulation. BMD in the lumbar spine and55 total hip had increased, on average, by about 8 and 5 %, respectively, but 3 fell back to or near baseline within 12 months of stopping therapy. The important clinical question is whether this interval of high turnover and rapid bone loss results in an abnormally high fracture risk −a rebound to a level of risk higher than just a return to the pre-treatment fracture60 status. Based on the above clinical evidence it is suggested that patients considered at high fracture risk should either continue denosumab therapy for up to 10 years or be switched to an alternative treatment. For patients at low risk, a decision to discontinue denosumab could be made after 5 years,65 but bisphosphonates therapy should be considered to reduce or prevent the rebound increase in bone turnover. While there is not sufficient evidence that discontinuation of denosumab treatment is linked to increased fracture risk, it is imperative to investigate the effect of drug holiday on bone biomarkers including BTMs and BMD.70 However, the latter effects are difficult to be addressed using clinical trials alone, given the large number of parameter combinations, including dose level, duration of drug holiday and dose levels after resuming therapy, involved to answer this question. We further hypothesized that the safest option would be to continue denosumab treatment on a lower dose regimen.75 The rational behind this type of treatment is that any of the negative side effects of denosumab (e.g., ONJ, atypical fractures, delayed healing of dental implants etc.) have only been reported to occur on high dose treatment and that if one could reduce the dosing regimen one would be able to retain the positive effect of the drug on BMD without obtaining a resorption overshoot80 due to discontinuation of treatment. In order to investigate the efficacy of different denosumab treatment patterns, including discontinuation and reduction of dose on bone turnover (as expressed in bone cell numbers), bone mineral density (BMD) and mineral content (denoted as α), we utilized a recently developed mechanistic PK-PD85 model of denosumab treatment of PMO [11]. This model allows to quantify the effect of denosumab on bone turnover and BMD taking into account the mineralisation process. The latter mechanism has been shown of fundamental importance for interpretation of obtained BMD gains following antiresorptive treatment. The model was recalibrated to reproduce the clinical90 results of Miller et al. [9], by adjusting: 1) the accessibility of denosumab to bone tissue compartments, 2) the bone mineralisation rate and 3) the severity of PMO, measured by the abnormal rate of RANKL expression. Utilizing this model, we investigate a wide variety of different denosumab 4 treatment patterns (summarized in Table 1) including the treatments anal-95 ysed by Miller et al.[9] (T1 and T2), which were used to calibrate the model. We then designed a range of virtual denosumab treatment regimes (T3-T9) to investigate different aspects of drug holidays and the subsequent commencement of therapy. Physical activity has a strong influence on bone remodelling through dif-100 ferent mechanisms that could also be affected by the denosumab treatment. For this reason, the interaction of drug dosing and exercise, through variations in the applied external load, is investigated at the end of the Results section. The paper is organized as follows: Section 2 summarises the most signifi-105 cant features of the mechanistic PK-PD model, in subsection 2.1, and a brief description of the simulated cases in subsection 2.2. In Subection 3.1 of the Results the calibration of two parameters is performed so that the numerical results are as close as possible to the clinical data. Then, the numerical predictions for different denosumab treatments is analysed in Subsection 3.2. A110 comprehensive discussion on the results is provided in Section 4. Conclusions and outlook to future work is presented in Section 5. 2. Materials and Methods 2.1. Mechanistic PK-PD model for simulation of the effect of denosumab on bone remodelling115 A brief description of the model is provided next. For a more detailed explanation, the reader is referred to the Appendix and to [11] where it was recently published. Following the approach taken by Pivonka et al. [12–14], the bone remodelling process can be described by cell balance equations. The bone cell types considered as state variables in the current model are: (i) os-120 teoblast precursor cells (OBp), (ii) active osteoblasts (OBa) and (iii) active osteoclasts (OCa), while the cell pools of uncommitted osteoblasts (OBu) and osteoclast precursors (OCp) are assumed constant. A set of differential equations governs the temporal change of cell populations, regulated by the following biochemical factors: PTH, TGF-β, RANKL, RANK and OPG,125 whose concentrations are also variables of the model. Parathyroid hormone (PTH) is a hormone secreted by the parathyroid glands to regulate serum calcium concentration by its action on calcium absorption through the intestine, recirculation/filtration in the kidneys and retrieval from bone. Transforming 5 growth factor beta (TGF-β) is a cytokine involved in a great number of cel-130 lular functions including the upregulation of osteoclasts apoptosis and differentiation of uncommitted osteoblast progenitors into responding osteoblasts as well as the downregulation of differentiation of responding osteoblasts into active osteoblasts. Receptor activator of nuclear factor κB (RANK) and its ligand (RANKL) are cytokines that control the differentiation of osteoclasts135 precursors into active osteoclasts. Finally, osteoprotegerin (OPG) is a decoy receptor for RANKL and competes with RANK for binding to RANKL, so completing the signalling pathway RANKL-RANK-OPG. The main outcome of cell balance equations is to provide the concentrations OCaand OBato be used in the following bone volume balance140 equation: dfbm dt =−kres ·OCa+kform ·OBa; (1) where fbm is the bone volume fraction and the constants kres and kform are respectively the bone matrix volume resorption rate and osteoid volume formation rate. As done in previous works [12–14], the signalling RANK-RANKL-OPG145 pathway is considered in the model to control osteoclasts activation through RANK-RANKL binding, as well as the action of denosumab on bone remodelling via competitive binding reactions within the RANK-RANKL-OPG pathway [3,15,16]. So, higher concentrations of denosumab give rise to lower concentrations of RANKL-RANK complexes and, hence, to a lower150 osteoclasts activation rate. Denosumab concentration in plasma is estimated by means of the onecompartment model proposed by Marathe et al. [15]. The amount of drug Doseden received in a subcutaneous injection and given in ng per kg of body weigth is used in the model to calculate the blood serum concentration, CP,den, as a function of time (see Appendix). A body weight of 60 kg was assumed to calculate Doseden.3Once the denosumab concentration in blood serum is known, the concentration of RANKL is worked out through: RANKL = RANKLmax βRANKL +PRANKL βRANKL +˜ DRANKLRANKLmax · 3In such case, the typical dose of 60 mg results in Doseden = 106ng of denosumab per kg of body weight. 6 ·1 + Ka,[RANKL-OPG]OPG + Ka,[RANKL-RANK]RANK + ζKa,[RANKL-den]Cp,den−1; (2) where Ka,[RANKL-OPG],Ka,[RANKL-RANK], and Ka,[RANKL-den] are the equilibrium association 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. βRANKL is the production155 rate of RANKL on the surface of osteoblasts, ˜ DRANKL is the RANKL degradation rate, RANKLmax is the maximum concentration of RANKL and PRANKL models the increase in RANKL production either induced by PMO (PPMO RANKL) or mechanical disuse (Pmech RANKL), such that PRANKL =PPMO RANKL +Pmech RANKL. Model parameters for the competitive binding model are given in Table 2160 of the Appendix. In Eq. (2) all association binding constants Kawere obtained from in vitro experiments [17]. Furthermore, RANKL production and degradation rate constants were estimated from previous works [12,13,18]. Hence, the only remaining factor to estimate was ζ, the accessibility factor of denosumab165 from blood serum to bone tissue. This accessibility factor is defined such that ζ= 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 can be site specific. We have performed a parametric study to quantify appropriate170 values for ζ, which will eventually act as a calibration factor. It was shown in a previous in silico study [11] that bone density rises with denosumab treatments, with the mineralization of the old tissue playing a key role in that rise, since mineral remains longer within bone because of the hindering effect of denosumab on bone turnover. Other mechanisms175 must contribute to that increase, as bone mineral density continued to grow in those simulations after the mineral content reached its maximum value, but their contribution seemed less important. A mineralization algorithm developed in [19] was used in that in silico study to estimate the temporal evolution of tissue mineral content during the remodelling process. This algo-180 rithm distinguishes three phases of the mineralization process: (i) an initial phase, called mineralisation lag time, that lasts from 6 to 22 days [20,21] during which no deposition of mineral occurs; (ii) a primary phase, which is very quick (it takes a few days to reach the 70% of the maximum mineral content [22]), and (iii) a secondary phase, when mineral is added at a decreas-185 ing rate [23], as the tissue becomes saturated with mineral. This secondary 7 phase may last from 6 months [24] to several years [25]. The rate of mineral deposition during this phase, measured by parameter κ(see Appendix), has a strong influence on bone density gain produced by the treatment, an influence which will be analysed in the referred parametric study along with190 the influence of parameter ζ. The mineral content is usually given by the so-called ash fraction, which can be measured in cadaver specimens or in biopsies. This variable is defined as the ratio between the mass of mineral mm(or ash mass) and the dry mass (the sum of inorganic and organic mass): α=mm mm+mo (3) The reader is referred to the Appendix for a detailed calculation of α.195 Finally, stiffness affects the strain energy density (SED) bone is subjected to and SED, in turn, affects mechanoregulation via the anabolic and catabolic regulatory functions, Πmech act,OBpand πRANKL act,OCp(see Appendix and Scheiner at al [16] for details). Bone stiffness is correlated with bone volume fraction and mineral content (through ash fraction, α) using Hernandez expression [26],200 which assumes an isotropic behaviour and provides the following Young’s modulus: E(MPa) = 84370 f2.58 bm α2.74 (4) A constant uniaxial compressive stress σwith different values (see Table 2) was applied on trabecular bone. This value together with the Young’s modulus given in Eq. (4) allows calculating the strain energy density, Ψ =205 σ2 2E. This corresponds to a simplified load case. The implementation of the model to a whole vertebral body or femur would require considering the distribution of stresses, as well as the distribution of cells and biochemical factors, for example in a Finite Element subroutine, and is out of the scope of the present paper.210 In order to simulate disease progression in PMO with subsequent denosumab treatment we have followed the approach described in detail in [11– 13,27], by investigating bone cell numbers, apparent density and ash fraction as biomarkers. While the former are representative for non-specific bone resorption and formation markers, the latter two are bone specific and reflects215 the material properties. PMO was simulated by introducing PPMO RANKL(t), a disease-related increase in RANKL production over time [16]: PPMO RANKL(t) = PPMO,ini RANKL φRANKL PMO (t) (5) 8 where PPMO,ini RANKL is the RANKL production at the onset of the disease (t− tPMO,ini) and φRANKL PMO is a reduction factor defined as: φRANKL PMO (t) = ξ2 ξ2+ht−tPMO,ini τRANKL PMO i2(6) where the shape factor, ξ, and the characteristic time of RANKL production,220 τRANKL P MO , determine the shape of the Lorentz-type function given by Eq. (6). The values ξ= 65 and τRANKL P MO = 10 days were taken from [16]. One year of PMO was simulated prior to the beginning of the treatment, which consisted in subcutaneous injections of denosumab in different dosages and frequencies. 2.2. Description of the simulated cases225 The analysed treatment patterns are summarized in Tab. 1and motivated by the two experimental studies of Miller et al. [9] and Bone et al. [8] together with virtual treatment patterns that we have assigned for testing of denosumab efficacy (T3-T9). The terminology used throughout the paper is the following: D1denotes the initial dosage, F1denotes the (SC) injection230 frequency, LT1denotes the length of the first treatment regimen, LDH1 denotes the length of drug holiday, D2denotes the denosumab dose after drug holiday, and F2denotes the frequency of the newly commenced treatment. Note that T1 is the uninterrupted denosumab treatment regimen reported by Miller et al [9], that implements the most typical dose (60 mg every 6235 months) used in clinical practice. T2 is a regimen also reported in [9], with a lower initial dose, a drug holiday at year 2 and resuming the treatment with the typical dose. Among the virtual patterns, in T3 the treatment is indefinitely stopped after 4 years; in T4 through T8 the drug holiday is implemented after 4 years of treatment, with different lengths or dosage after240 recommencing treatment; finally, in T9 there is no drug holiday, but instead the dose is decreased after 4 years of treatment. The evolution of apparent density predicted by the model, ρ(t), was used to calculate the bone density gain achieved with the treatment, BDG, which is the main result given in clinical studies and is obtained by comparison245 with baseline, i.e. the density measured at the beginning of treatment, after 1 year of PMO in the present simulations: BDG(%) = ρ(t)−ρ(1 year) ρ(1 year)·100 (7) 9 6 12 24 36 48 60 72 84 96 108 120 months 0 5 10 15 20 BDG (%) 6 12 24 36 48 60 72 84 96 108 120 months 0.62 0.64 0.66 0.68 0.7 0.72 Figure 5: Effect of treatment pattern after the drug holidays. Evolution of BDG (left) and ash fraction (right) in the simulations of treatments T4, T7, T8 and T9 for PPMO,ini RANKL = 3000 pM, ζ= 0.006 and κ= 0.0045 over a 10 years period. density. Actually, exercise has a well-known anabolic effect and apparent density rises with the stress level, but BDG is measured in percentage, relative to the density at the beginning of the treatment, and this is much lower for σ= 4 MPa, especially in the lumbar spine. The influence of these small365 stress changes is not very pronounced on the mineral content, which is only slightly increased by the stress. 4. Discussion The PK-PD model developed in [11] was further calibrated to predict the clinical results obtained by Miller et al. [9] in patients with PMO and370 treated with denosumab, either uninterruptedly or with prescribed drug holidays. Parameter optimization was performed to investigate the effect of major factors regulating bone remodelling and drug action, among these we analysed the accessibility factor, ζ, the mineralisation rate, κ, and severity of PMO, PPMO,ini RANKL . The simulation results follow the global trend of the clinical375 results, though quantitative differences were found. These differences may be partly explained by two factors: the wide dispersion of the clinical results and some aspects not considered in the present study due to lack of data. Some of these aspects stand out from the rest: 16 6 12 24 36 48 60 72 84 96 108 120 months 0 5 10 15 20 25 BDG (%) 6 12 24 36 48 60 72 84 96 108 120 months 0.6 0.62 0.64 0.66 0.68 0.7 0.72 6 12 24 36 48 60 72 84 96 108 120 months 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 Figure 6: Combined effect of exercise and denosumab treatment. Evolution of BDG (top), ash fraction (middle) and apparent density (bottom) in treatment T8 for different stresses, PPMO,ini RANKL = 3000 pM, ζ= 0.006 and κ= 0.0045. 17  Patient’s body weight. Denosumab dose is given in mg of active com-380 pound per injection, but its serum concentration is controlled by the patient’s body weight. An average of 60 kg was assumed, but its variability could be large and was not measured in the data by Miller et al. [9]. Body weight should be recorded along with dose in future clinical studies.385  Severity of PMO disease progression. The increase in RANKL production caused by PMO (PPMO,ini RANKL ) could vary among patients indicating different disease severity, but it is hard to measure. It is also difficult to assess the time elapsed from the onset of disease to the beginning of the treatment. Both represent important sources of uncertainty for390 the numerical results as they have a strong influence on BDG. Particularly the onset of disease, because it determines the bone density at the beginning of the treatment with respect to which BDG is measured.  Patient’s activity level. As seen in Fig. 6, the stress level affects very importantly bone density and consequently BDG. An average value of395 5 MPa was arbitrarily chosen for the stress that trabecular bone is subjected to, but this could vary among patients, though it was not known from the clinical study.  Initial fbm. We have used fbm = 25% and fbm = 15% for hip and spine simulations, respectively, but these values could be patient specific and400 have a strong influence on the results.  Values of model parameters. Model parameters were taken from previous studies and represent average values of a population, but some of them could be patient specific and might have a strong influence on the numerical results. A sensitivity analysis of these constants on BDG405 and mineral content should be conducted in future studies to quantify that influence. All in all, the analysis of treatments T1 and T2 (see Figs. 2and 3) highlighted the fact that the predicted BDG is much higher in the porous trabecular bone (lumbar spine), because it represents a percentage increase410 with respect to the initial density, which is much lower than in the hip. However, the increase of bone density in absolute terms is primarily due to bone mineralisation, as shown in a previous study [11], and is bigger 18 in the hip, because the mineral content is higher (Fig. 3). The obtained ash fraction (around 0.73) is very high and could jeopardize the integrity of415 bone, as stated above. Stopping the treatment at month 24 has a negative effect on BDG, but it is beneficial for the mineral content, which is prevented from increasing dangerously. Nonetheless, soon after resuming the treatment, the initial trend is recovered with a steady increase in both, BDG and α. In the case of α, the value reached right before the interruption is clearly420 surpassed soon after resuming the treatment and the subsequent evolution almost mimics that of the treatment without interruption. In other words, the benefits of the drug holiday seem to be only temporary in this case and, for this reason, alternative strategies of the drug holiday need to be studied in order to achieve an optimum treatment.425 It can be noted that the error in the BDG estimation is higher for the spine and more notably for T2 (see Fig. 2). This can be due again to the variability of clinical results, which are higher for the spine and make them more difficult to adjust. Clinical studies of anti-catabolic drug efficacy generally only measure430 BMD (before and after treatment) and reduction of fractures as the final goal of the treatment. However, they do not measure the mineral content in the bone matrix, which is a key factor contributing to BMD and determining the mechanical properties of bone and its suitability to bear the loads under normal physiological activity. It is well known that anti-catabolic435 drugs reduce bone turnover, thus enabling the existing bone matrix to embed more mineral in the extrafibrillar space and so increasing bone tissue density and, hence, BMD. A highly mineralised bone matrix, in turn, gives rise to a more brittle material behaviour, increasing the stiffness and reducing the fracture toughness of bone [30]. This fact is confirmed by experimental440 studies that measured a higher amount of microstructural damage in interstitial bone, which has a higher mineral content [31–33]. Very importantly too, bone turnover suppression alters the tissue repair process [34] and the accumulation of unrepaired microstructural damage results in unimpeded crack progression and leads to the development of AFFs [35], which have445 been associated with anti-catabolic treatments [36]. Our numerical results also support that there could be a link between bone turnover and degree of mineralisation on the one hand and the potential development of AFFs on the other hand, due to long-term suppression of bone remodelling. More precisely, the predicted ash fraction was higher in the hip than in the lumbar450 spine, due to a lower bone turnover in the hip compared to lumbar spine. 19 This could explain why atypical fractures occur preferentially in the femur. From this point of view, an interruption of the anti-resorptive treatment (drug holidays) over a prolonged period of time may seem a good strategy to re-invoke bone remodelling, restore normal mineral content and repair455 microstructural damage, so reducing the risk of AFFs. In accordance with the clinical results, the in silico model showed that BMD decreases rapidly upon discontinuation of denosumab treatment and also that this BMD loss is effectively restored soon after the denosumab treatment is resumed. However, the clinical results do not provide any infor-460 mation about the changes in mineral content with and without treatment. Based on the experimental data of Bone et al. [8] and looking at changes in bone turnover based on histological data, we infer that the BMD gain obtained during denosumab treatment is mostly due to an increase of tissue mineralisation with porosity staying largely constant. On the other hand,465 introduction of a drug holiday allows for bone remodelling to be re-invoked at trabecular surfaces as demonstrated re-appearance of tetracycline labels [8]. Our numerical simulations reflect this behaviour, based on changes in porosity, mineral content and BMD. Hence, if the aim of drug discontinuation is to restore bone matrix properties to pre-treatment level, our simulations470 clearly show that this effect is quickly lost after resuming therapy, so only providing temporary benefits. We also showed that there is no benefit in extending the drug holidays over 12 months based on no significant changes in mineral content after that time. Based on this result it seems more adequate not to prescribe a drug475 holiday, but rather to reduce denosumab treatment at a lower dosage. This strategy was investigated in case T9 where we reduced the denosumab dosage to 30 mg at a 6 month injection interval. Using this strategy a compromise between minimizing BMD loss and stabilization of the mineral content to a moderate level could be achieved. The bone material properties for this latter480 treatment case (T9) are less brittle and have a higher fracture toughness. Consideration of the tissue mineral content is necessary to assess the likeliness of AFFs and fractures in general. Experimentally, this is difficult to achieve due to the fact that bone biopsies can only be performed on iliac crest samples and that one cannot access bone samples in vivo at the femoral neck485 and lumbar spine. One possibility may be to develop a database of diseased patients that have been on anti-resorptive therapies for a long time. This could provide an upper bound estimate for maximum degree of mineralisation at these respective bone sites. Also one would require information on baseline 20 mineral content measures in a healthy population as a reference measure.490 Thus, from a modelling perspective, assessing the likeliness of AFFs would be necessarily grounded on phenomenological considerations, for instance, by establishing a predefined bound for the degree of bone mineralisation, which in case of being exceeded would mean a high risk of fracture. However, this would be a strong simplification of the problem. An alternative approach to495 predict occurrence of (low energy) fractures would be to include a damage model within the bone remodelling algorithm, like those proposed in [19,37]. Utilizing such a model would allow to assess the combined effects of changes in porosity, mineral content and loading with or without drug treatments on damage accumulation in the bone matrix. The latter model parameter500 could be linked to the characteristic time of occurrence of AFFs. Note that in the current model this link was indirectly made using the bone mineral content parameter. The damage model referred to before [19,37] considers damage repairing through bone turnover. Therefore, the use of that model would allow to evaluate the incidence of a decreased bone turnover on the505 accumulation of microstructural damage and the occurrence of AFFs. It was also shown that physiological activity level, expressed as the mechanical stress acting on trabecular bone, had no significant influence on the mineral content. The latter quantity only slightly increased with increasing mechanical stress. However, including a damage law in our model would510 allow to create a relationship between bone matrix stress and damage accumulation in the bone matrix. We believe, that such a model formulation would provide mechanistic insight into the effects of physiological activity and occurrence of AFFs. In the task force report by ASBMR it was pointed out that bone regions subjected to tensile strains are associated with atypi-515 cal femoral fractures [38]. More recently, Martelli et al. showed that peakstrains and stresses in the femur strongly depend on the type of daily activity (e.g. walking, stair up, stair down etc.) and the number of repetitions of a certain activity [39]. In that work, walking (the most frequently performed daily activity) was identified as prime candidate which was shown to induce520 the highest tensile strain (2004 µϵ) in the lateral femoral shaft [39]. Linking the later information obtained from an organ scale musculoskeletal (MSK) model with a tissue scale model of bone remodelling including the mineralisation process and damage accumulation in the bone matrix would be a way to diagnose likelihood of AFF occurrence. Development of such a multiscale525 model is beyond the scope of the current paper and aimed for future research. The main limitation of the present study is related to the values assigned 21 to some of the model parameters. In particular, the parameters analysed here, i.e. denosumab accessibility factor, ζ, mineralisation rate constant, κ, and PMO severity, PPMO,ini RANKL have been estimated based on best fit with ex-530 perimental results. Other model parameters, not varied in current study, could also have a potential significant effect on numerical results and might help explain the variability of clinical results. A comprehensive parameter sensitivity analysis will be the scope of future studies. Nonetheless, due to the fact that many model parameters are directly related to bone site-535 specific remodelling, experimental data is difficult to obtain. Clinical studies do not provide the required detail on the remodelling process, but help identify patient-specific parameters such as patient’s body weight, physiological activity level and onset of the bone disease. However, clinical trial data do not provide individual patient data, but provide information on a population540 scale. Consequently, our modelling results refer to population scale. In the future, it would be beneficial to obtain individual subject data and then use the in silico model in a patient-specific way to predict BMD, mineral content, and site-specific bone turnover. As stated before, another limitation of this study is that the model is ap-545 plied at the global (or body) level, in a representative volume element (RVE) without any specificity. The applied mechanical loads were simplistic and the RVE was not affected by the surroundings. At the organ level the distribution of stresses would play a key role, which would affect the distribution of cells and biochemical regulatory factors. Such a study would be of much550 interest, particularly for the evaluation of fracture risk, but it would require the implementation of the model in a Finite Element subroutine and is left for future studies. 5. Summary and conclusions The numerical simulation results have provided the following novel in-555 sights into the effects of drug holiday and resuming therapy on bone health:  Duration of drug holiday. There is no use in extending the drug holidays over 1 year because the evolution of mineral content that is observed after resuming the treatment quickly catches up with those obtained in treatments resumed before.560  Resumed treatment effects. It seems adequate to reduce the dose after resuming the treatment in terms of reducing the risk of fracture due to 22 an excessive mineral content.  Site-specific bone effects. Drug treatment and drug holidays have different effects depending on bone site.565  Patient’s activity level. Simulating altered physiological activity (based on applied mechanical stress) only had a minor effect on BMD, BTMs, and mineral content. Increased mechanical stress leads to increased BMD and bone mineral content as a result of the mechanobiological feedback included in the model.570  Continuation of treatment at lower dose. The most promising virtual treatment turned out to be not including a drug holiday, but only reducing the denosumab dose by half. Switching from 60 mg every 6 months to 30 mg every 6 months provided the best compromise between maintaining BMD, reducing mineral content and allowing for some bone575 remodelling to take place. However, it must be noted that the benefit of this strategy should be confirmed by an evaluation of fracture risk that included damage and, more importantly, by clinical trials.  Existence of an optimal treatment. The results suggest the existence of an optimal treatment that could have a varying dose with time in580 order to keep bone density gain and mineral content within reasonable values. This optimal treatment could also be patient’s specific. We note that the current model does not include either damage formation in the bone matrix due to physiological activity or damage repairing through bone turnover. Damage (i.e., microcracks) has been linked to targeted re-585 modelling and could provide further insights into the link between mineralisation, BMD, physiological loading and denosumab treatment patterns with and without drug holidays. In future studies we aim to explore these effects in more detail, identifying safe dosing patterns in terms of BTMs, BMD, mineral content and damage accumulation/repair. We also aim to confirm590 the hypothesis that an excessive mineral content along with the accumulation of unrepaired microstructural damage could explain AFFs. Acknowledgement The authors gratefully acknowledge the research support from the Ministerio de Ciencia e Innovaci´on of the Spanish Government through the595 23 research project PID2019-106969R entitled Modelos de Remodelaci´on ´ Osea y su Aplicaci´on para Simular el Efecto de Tratamientos Antireabsortivos contra la Osteoporosis Postmenop´ausica. Appendix In this appendix a brief description of the mathematical model is included.600 For a more detailed explanation the reader is referred to [11] where it has been recently published. Model of bone cell interactions in bone remodelling The model describes bone cell interactions involved in bone remodelling, also taking into account bone mineralisation of the newly formed tissue. As in previous models, the RANK-RANKL-OPG pathway was considered, together with the action of several regulatory factors on bone cells, including TGF−βand mechanobiological feedback. In the approach proposed by Pivonka et al. [12–14] the bone remodelling process can be described as cell balance equations, using the following bone cell concentrations as state variables: (i) osteoblast precursor cells (OBp), (ii) active osteoblasts (OBa) and (iii) active osteoclasts (OCa). The cell pools of uncommitted osteoblasts (OBu) and osteoclast precursors (OCp) are assumed constant: dOBp dt =DOBu·πT GF −β act,OBu ·OBu+POBp·Πmech act,OBp·OBp+ −DOBp·πT GF −β rep,OBp·OBp; (9) dOBa dt =DOBp·πT GF −β rep,OBp·OBp−AOBa·OBa; (10) dOCa dt =DOCp·πRANKL act,OCp·OCp−AOCa·πT GF −β act,OCp·OCa; (11) DOBu,DOBp,DOCpare differentiation rates of OBuinto OBp,OBpinto OBa and OCpinto OCa, respectively. POBpis the proliferation rate of OBp.AOBa 605 and AOCaare apoptosis rates of OBaand OCarespectively. The variables πT GF −β act,OBu,πT GF −β rep,OBpand πT GF −β act,OCprepresent activator and repressor functions related to the binding of TGF-βto its receptor. Similarly, πRANKL act,OCpis the activator function related to the RANK-RANKL binding. We note that the cell balance equations are composed of a production term and a linear610 degradation term, which describes differentiation of one cell type into another 24 (or terminal cell fate, i.e. apoptosis). The linear degradation term is further regulated via sigmoidal activator/repressor functions which take values from zero to one. Model parameters of the cell population model are given in Table 2.615 Πmech act,OBpis the function regulating the anabolic part of the mechanobiological feedback. This variable quantifies the influence of mechanical stimulus on the proliferation of preosteoblasts, which eventually affects the population of active osteoblasts and the amount of tissue being formed: Πmech act,OBp=˜ Πmech act,OBp1 + λ−1Ψbm ˜ Ψbm −1 (12) where Ψbm is SED at the RVE and ˜ Ψbm is a reference SED that distinguishes620 normal load from overload. Following Scheiner et al. [16] the adopted values were: ˜ Πmech act,OBp= 0.5 and λ1= 0 if Ψbm ≤˜ Ψbm (normal load) or λ2= 1.25 if Ψbm >˜ Ψbm (overload). The catabolic part of the mechanobiological feedback is given by the RANKL production due to disuse:625 Pmech RANKL =λ21−Ψbm ˜ Ψbm (13) where λ2= 105pM/day if Ψbm ≤˜ Ψbm and λ2= 0 if Ψbm >˜ Ψbm [16]. 5.1. Denosumab action on RANK-RANKL-OPG: competitive binding Similar to previous mechanistic PK-PD models, we incorporate the action of denosumab on bone remodelling via competitive binding reactions within the RANK-RANK-OPG pathway [3,15,16]. As stated in section 2 denosumab competes with RANK (and OPG) for binding to RANKL. Consequently, higher concentrations of denosumab give rise to lower concentrations of RANKL-RANK complexes. The concentration of RANKL is given by Eq. (2) and, with that, the activator function of RANKL can be expressed as: πRANKL act,OCp=[RANKL-RANK] Kd,[RANKL-RANK] + [RANKL-RANK] (14) with Kd,[RANKL-RANK] as the corresponding equilibrium dissociation binding constant. Hence, lower concentrations of RANKL-RANK complexes leads to a lower OCaactivation rate (see Eq. (11)).630 25 [9] Miller PD, Bolognese MA, Lewiecki EM, McClung MR, Ding B, Austin M, et al. Effect of denosumab on bone density and turnover in postmenopausal women with low bone mass after long-term continued, dis-760 continued, and restarting of therapy: A randomized blinded phase 2 clinical trial. Bone 2008;43(2):222–9. URL: http://dx.doi.org/10. 1016/j.bone.2008.04.007. doi:10.1016/j.bone.2008.04.007. [10] McClung MR. Cancel the denosumab holiday. Osteoporosis International 2016;27(5):1677–1682. URL: https://www.ncbi.nlm.nih.gov/765 pubmed/26932443. doi:DOI10.1007/s00198-016-3553-3. [11] Mart´ınez-Reina J, Pivonka P. Effects of long-term treatment of denosumab on bone mineral density: insights from an in-silico model of bone mineralization. Bone 2019;125(125):87–95. URL: https://www.ncbi. nlm.nih.gov/pubmed/31055117. doi:10.1016/j.bone.2019.04.022.770 [12] Pivonka P, Zimak J, Smith D, Gardiner BS, Dunstan C, Sims NA, et al. Model structure and control of bone remodeling: A theoretical study. Bone 2008;43(2):249–63. URL: http://dx.doi.org/10.1016/j.bone. 2008.03.025. doi:10.1016/j.bone.2008.03.025. [13] Pivonka P, Zimak J, Smith D, Gardiner BS, Dunstan C, Sims NA,775 et al. Theoretical investigation of the role of the rank–rankl–opg system in bone remodeling. Journal of Theoretical Biology 2010;262(2):306–16. URL: http://www.sciencedirect.com/science/ article/pii/S002251930900441X. doi:https://doi.org/10.1016/j. jtbi.2009.09.021.780 [14] Scheiner S, Pivonka P, Hellmich C. Coupling systems biology with multiscale mechanics, for computer simulations of bone remodeling. Computer Methods in Applied Mechanics and Engineering 2013;254(1):181– 96. URL: http://www.sciencedirect.com/science/article/pii/ S0045782512003234. doi:https://doi.org/10.1016/j.cma.2012.10.785 015. [15] Marathe DD, Marathe A, Mager DE. Integrated model for denosumab and ibandronate pharmacodynamics in postmenopausal women. Biopharmaceutics and Drug Disposition 2011;32(8):471–81. URL: http: //dx.doi.org/10.1002/bdd.770. doi:10.1002/bdd.770.790 32 [16] Scheiner S, Pivonka P, Smith DW, Dunstan CR, C. H. Mathematical modeling of postmenopausal osteoporosis and its treatment by the anticatabolic drug denosumab. International Journal for Numerical Methods in Biomedical Engineering 2013;30(1):1–27. doi:10.1002/cnm.2584. [17] Kostenuik P, Nguyen H, McCabe J, Warmington K, Kurahara C, Sun795 N, et al. Denosumab, a fully human monoclonal antibody to rankl, inhibits bone resorption and increases bmd in knock-out mice that express chimeric (murine/human) rankl. Journal of Bone and Mineral Research 2009;2(24):182–95. [18] Pivonka P, Buenzli P. R. Scheiner S, Hellmich C, Dunstan C. The influ-800 ence of bone surface availability in bone remodelling - a mathematical model including coupled geometrical and biomechanical regulations of bone cells. Engineering Structures 2013;47(1):134–147. [19] Mart´ınez-Reina J, Garc´ıa-Aznar J, Dom´ınguez J, Doblar´e M. On the role of bone damage in calcium homeostasis. J Ther Biol 2008;254(3):704–12.805 [20] Eriksen E, Hodgson S, Eastell R, Cedel S, O’Fallon W, Riggs B. Cancellous bone remodeling in type I (postmenopausal) osteoporosis: quantitative assesment of rates of formation, resorption and bone loss at tissue and cellular levels. J Bone Miner Res 1990;5(4):311–9. [21] Need A, Horowitz M, Morris H, Moore R, Nordin C. Seasonal change810 in osteoid thickness and mineralization lag time in ambulant patients. J Bone Miner Res 2007;22(5):757–61. [22] Hernandez C, Beaupr´e G, Carter D. A model of mechanobiologic and metabolic influences on bone adaptation. J Rehabil Res Dev 2001;37(2):235–44.815 [23] Parfitt A. The physiologic and clinical significance of bone hisotomorphometric data. In: Recker R, editor. Bone Histomorphometry Techniques and Interpretation. Boca Raton, FL, USA: CRC Press; 1983, p. 143–223. [24] Parfitt A. Bone remodeling and bone loss: Understanding the patho-820 physiology of osteoporosis. Clin Obstet Gynecol 1987;30(4):789–811. 33 [25] Frost H. Bone Remodeling Dynamics. Springfield, IL: Charles C. Thomas Co.; 1963. [26] Hernandez C, Beaupr´e G, T.S. K, Carter D. The influence of bone volume fraction and ash fraction on bone strength and modulus. Bone825 2001;29(1):74–8. [27] Lemaire V, Tobin FL, Greller LD, Cho CR, Suva LJ. Modeling the interactions between osteoblast and osteoclast activities in bone remodeling. Journal of Theoretical Biology 2004;229(3):293–309. URL: http://www. sciencedirect.com/science/article/pii/S0022519304001407.830 doi:https://doi.org/10.1016/j.jtbi.2004.03.023. [28] Hildebrand T, Laib A, M¨uller R, Dequeker J, R¨uegsegger P. Direct three-dimensional morphometric analysis of human cancellous bone: microstructural data from spine, femur, iliac crest and calcaneus. Journal of Bone and Mineral Research 1999;7(14):1167–74. URL: https://www.835 ncbi.nlm.nih.gov/pubmed/10404017. doi:10.1359/jbmr.1999.14.7. 1167. [29] Currey J. Tensile yield in compact bone is determined by strain, postyield behaviour by mineral content. J Biomech 2004;37(4):549–56. [30] Bala Y, Farlay D, G. B. Bone mineralization: from tissue to crystal in840 normal and pathological contexts. Osteoporos Int 2013;24(8):2153–66. [31] Boyce T, Fyhrie D, Glotkowski M, Radin E, Schaffler M. Damage type and strain mode associations in human compact bone bending fatigue. J Orthopaed Res 1998;16(6):322–329. [32] O’Brien F, Taylor D, Dickson G, Lee T. Visualisation of three-845 dimensional microcracks in compact bone. J Anat 2000;197(3):413–420. [33] Qiu S, Rao D, Fyhrie D, Palnitkar S, Parfitt A. The morphological association between microcracks and osteocyte lacunae in human cortical bone. Bone 2005;37(1):10–5. [34] Mashiba T, Hirano T, Turner C, Forwood M, Johnston C, Burr D.850 Suppressed bone turnover by bisphosphonates increases microdamage accumulation and reduces some biomechanical properties in dog rib. J Bone Miner Res 2000;15(4):613–620. 34 [35] Ettinger B, Burr D, Ritchie R. Proposed pathogenesis for atypical femoral fractures: lessons from materials research. Bone 2013;55(2):495–855 500. [36] Saita Y, Ishijima M, Kaneko K. Atypical femoral fractures and bisphosphonate use: current evidence and clinical implications. Therapeutic Advances in Chronic Disease 2015;6(4):185–93. doi:10.1177/ 2040622315584114.860 [37] Mart´ınez-Reina J, Garc´ıa-Aznar J, Dom´ınguez J, Doblar´e M. A bone remodelling model including the directional activity of BMUs. Biomech Model Mechanobiol 2009;8(2):111–27. doi:10.1007/ s10237-008-0122-5. [38] Shane E, Burr D, Ebeling PR, Abrahamsen B, Adler RA, Brown TD,865 et al. Atypical subtrochanteric and diaphyseal femoral fractures: Report of a task force of the american society for bone and mineral research. Journal of Bone and Mineral Research 2010;25(11):2267– 94. URL: https://onlinelibrary.wiley.com/doi/abs/10.1002/ jbmr.253. doi:10.1002/jbmr.253.870 [39] Martelli S, Pivonka P, Ebeling PR. Femoral shaft strains during daily activities: Implications for atypical femoral fractures. Clinical Biomechanics 2014;29(8):869–76. URL: https://doi.org/10.1016/ j.clinbiomech.2014.08.001. doi:https://10.1016/j.clinbiomech. 2014.08.001.875 [40] Martin R. Porosity and specific surface of bone. Crit Rev Biomed Engl 1984;10(3):179–222. 35