Following the approach taken by Martin et al. (Martin et al., 2019), the bone remodelling process can be described as cell balance equations. The bone cell types (i.e. state variables) considered in the current model are: osteoblast precursor cells (Ob_p), active osteoblasts (Ob_a), osteoclast precursor cells (Oc_p), active osteoclasts (Oc_a) and osteocytes (Ot). The cell pools of uncommitted progenitor cells of both lineages (Ob_u, Oc_u) were assumed constant: d O b p d t = D Ob u ⋅ O b u ⋅ π a c t , O b u TGF−β − D Ob p ⋅ O b p ⋅ π r e p , O b p TGF−β (1) + P Ob p ⋅ O b p ⋅ π a c t , O b p Wnt d O b a d t = D Ob p ⋅ O b p ⋅ π r e p , O b p TGF−β − Δ Ob a ⋅ O b a (2) d O c p d t = D Oc u ⋅ O c u ⋅ π a c t , O c u RANKL − D Oc p ⋅ O c p ⋅ π a c t , O c p RANKL (3) d O c a d t = D Oc p ⋅ O c p ⋅ π a c t , O c p RANKL − A Oc a ⋅ O c a ⋅ π a c t , O c p TGF−β (4) d O t d t = η d f b m d t (5) where D Ob u , D Ob p , D Oc u and D Oc p are the differentiation rates of Ob_u, Ob_p, Oc_u and Oc_p, respectively; A Oc a is the apoptosis rate of Oc_a and Δ Ob a is the rate of clearance of active osteoblasts through apoptosis or differentiation into osteocytes. The variables π a c t , O b u TGF−β , π r e p , O b p TGF−β and π a c t , O c p TGF−β represent activator and repressor functions related to the binding of TGF-β to its receptor. Similarly, π a c t , O c u RANKL and π a c t , O c p RANKL are the activator functions related to the RANK-RANKL binding. Finally, P Ob p is the proliferation rate of Ob_p, a process which is mediated by the Wnt signalling pathway through the activator function π a c t , O b p Wnt and is described in detail in the Supplementary Material along with other features of the model. A schematic figure of the mechanistic PK-PD model is presented in Figure 1. Model parameters of the cell population model are given in the Supplementary Table S1.

Equation 5 establishes that the population of osteocytes varies proportional to the bone matrix fraction f _bm , given that the density of osteocytes is constant within the bone matrix, η, is assumed constant, as done in (Martin et al., 2019). Finally, the variation of bone matrix fraction is obtained through the balance between resorbed and formed tissue: d f b m d t = − k res ⋅ O c a + k form ⋅ O b a (6) where k _res and k _form are, respectively, the rates of bone resorption and osteoid formation (see Supplementary Table S1).

The RANK-RANKL-OPG signalling pathway controls the differentiation of uncommitted osteoclast progenitors and osteoclasts maturation, respectively through π a c t , O c u RANKL and π a c t , O c p RANKL (see Eqs. 3, 4. Thus, an imbalance in that pathway, such as that occurring after menopause, may result in the development of osteoporosis. The action of Dmab against the disease aims at restoring that balance. Following Martin et al. (Martin et al., 2019), the concentrations of OPG, RANK and RANKL are given by the following equations, which are explained in the Supplementary Material: O P G = P OPG D ̃ OPG + D ̃ OPG - RANKL R A N K L K OPG - RANKL (7) R A N K = N Oc p RANK O c p 1 + R A N K L K RANK - RANKL (8) R A N K L = P RANKL ⋅ D ̃ RANKL + D ̃ OPG - RANKL K OPG - RANKL ⋅ O P G + D ̃ RANK - RANKL K RANK - RANKL ⋅ R A N K + D ̃ Dmab - RANKL K RANKL - Dmab ⋅ [ D m a b ] BC − 1 (9) where the binding of Dmab to RANKL is considered in the last term of the bracket of Eq. 9; D ̃ X and D ̃ X - Y are the degradation rates of the factor X and the complex X-Y, respectively; K _X-Y is the dissociation constant of the complex X-Y and N Oc p RANK is the number of RANK receptors per osteoclast precursor. [Dmab]_BC is the concentration of Dmab in the bone compartment and will be analysed in Section 2.3. P _OPG is the production rate of OPG by active osteoblasts: P OPG = β OPG,Ob a π rep,Ob a PTH O b a 1 − O P G O P G max (10) where β OPG,Ob a is the OPG production rate, π rep,Ob a PTH is the repressor function that quantifies the effect of PTH on the production of OPG and [OPG_max] is the saturation concentration of OPG above which no further production takes place. To evaluate P _RANKL, the RANKL production rate of Eq. 9, we have assumed that RANKL is expressed by osteocytes and osteoblasts precursors, following experimental evidence (Nakashima et al., 2011; Xiong et al., 2015) and then: P RANKL = β RANKL,Ot O t 1 − [ R A N K L ] tot [ R A N K L ] max π act,RANKL dam + β RANKL,Ob p π act/rep,RANKL PTH,NO O b p 1 − [ R A N K L ] tot [ R A N K L ] max + P RANKL PMO (11) where P RANKL PMO is the RANKL production due to PMO; β _RANKL,Ot and β RANKL,Ob p are the RANKL production rate of osteocytes and osteoblast precursors, respectively; π act/rep,RANKL PTH,NO is a co-regulatory function that takes into account the up-regulation of RANKL transcription by the parathyroid hormone (PTH) and its inhibition by nitric oxide (NO) (Martin et al., 2019) and π act,RANKL dam is an activator function accounting for the upregulation of RANKL expression by osteocytes due to microstructural damage (Martínez-Reina et al., 2021b) (see the details of both regulatory functions in the Supplementary Material). Finally, [RANKL]_max is the saturation concentration of RANKL above which no further expression takes place and [RANKL] _tot is the total concentration of RANKL (bound and free) and is defined as follows: [ R A N K L ] tot = R A N K L ⋅ 1 + O P G K OPG - RANKL + R A N K K RANK - RANKL + [ D m a b ] BC K RANKL - Dmab (12)

Following (Martínez-Reina et al., 2021b), the RANKL production due to PMO, P RANKL PMO is modelled as a sigmoidal function of time: P RANKL PMO t = P RANKL , max PMO t − t onset γ t − t onset γ + δ PMO γ f o r t ≥ t onset 0 f o r t < t onset (13) where P RANKL , max PMO is the maximum (long-term) RANKL production due to PMO, t _onset is the time of onset of the disease, δ _PMO is a time constant establishing when the 50% of P RANKL , max PMO is reached and γ is the sigmoidicity factor.

Equations 11–13 can be substituted in Eq. 9 to work out the free RANKL concentration, i.e. [RANKL]. Then, the activator functions in Eqs. 3, 4 have the same structure: π act,X RANKL = [ R A N K L ] K act,X RANKL + [ RANKL ] with X  = O c u , O c p (14) and different constants: K a c t , O c u RANKL and K a c t , O c p RANKL . In previous works (Martínez-Reina and Pivonka, 2019; Martínez-Reina et al., 2021a; Martínez-Reina et al., 2021b) both constants K act,X RANKL were chosen equal in Eq. 3, so resulting in a constant pool of Oc_p in the steady state. One of the novelties of the current paper is to choose these constants different, so creating a variable pool of Oc_p. To analyse the effect of this change, the response of two models is compared: 1) the model proposed in this paper, termed as BCPM _ext (extended model) with different constants K act,X RANKL for Oc _u and Oc _p , leading to a variable pool of Oc_p and 2) a model, termed as BCPM _st (standard model) that coincides with the new model except that constants K _act,X are equal and thus produces a constant pool of Oc_p.

Several pharmacokinetic (PK) models of Dmab have been proposed including one- and two-compartment models (Dua et al., 2015). In recent papers (Martínez-Reina and Pivonka, 2019; Martínez-Reina et al., 2021a) we have used a one-compartment model with Michaelis-Menten kinetics in order to characterise the serum Dmab PK profiles, with the constants of that model fitted from a clinical study by Marathe et al. (Marathe et al., 2011). This model is framed in a dashed box in Figure 1. A first-order rate process (k _a ) governs the absorption of the drug (Dose) from the subcutaneous (SC) injection site into the central compartment ([Dmab] _CC ), being V _c /F the volume of the central compartment adjusted for bioavailability. The drug elimination from the central compartment is described by a combination of a linear first-order process (k _el1 ) and a non-linear saturation process (V _max, K _m ): d [ D m a b ] C C d t = D o s e V c / F k a e − k a t − k el1 [ D m a b ] C C + V m a x V c / F [ D m a b ] C C K m + [ D m a b ] C C (15)

In Eq. 15 Dose is given in ng per kg of body weight and then [Dmab] _CC is calculated in ng/ml and subsequently converted into pmol/l, through the molecular weight of Dmab M _Dmab = 149 kDa (Amgen). Supplementary Table S1 summarises the PK model parameters of the one-compartment model adjusted by Marathe et al. (Marathe et al., 2011).

In previous models (Martínez-Reina and Pivonka, 2019; Martínez-Reina et al., 2021a) we have assumed that a fraction of the Dmab present in the central compartment was available in the bone compartment to compete with RANK to bind to RANKL. However, this availability actually implied a reversible exchange of Dmab between both compartments. Given the affinity of Dmab for RANKL, which is expressed by osteoblasts precursors within the bone compartment, a flux from the central compartment to the bone compartment seems more plausible than a reversible exchange. To this end, we have added the bone compartment to the PK model (see Figure 1). The term in square brackets in Eq. 15 represents the elimination from the central compartment in the model of Marathe et al. We have assumed that only a fraction (1 − ζ) of this term is actually eliminated via urine and the rest, ζ, is the flux of Dmab into the bone compartment. In turn, the latter fraction can be considered as P D m a b B C , the production term of Dmab in the competitive binding reactions between RANKL, OPG and Dmab: P Dmab B C = ζ k el1 [ D m a b ] C C + V m a x V c / F [ D m a b ] C C K m + [ D m a b ] C C (16)

This production rate can be replaced in the expression that gives the concentration of ligands in competitive binding reactions (see Supplementary Material), i.e.: [ D m a b ] B C = P Dmab B C D ̃ Dmab B C + D ̃ RANKL - Dmab K RANKL - Dmab ⋅ [ R A N K L ] (17)

This expression can be used in Eq. 9 to give the concentration of free RANKL. All the parameters are known from previous works (Martínez-Reina and Pivonka, 2019), except for D ̃ Dmab B C and ζ that needs to be adjusted in this model as explained later on.

Targeted bone remodelling theories hypothesise that one of the major functions of bone remodelling is to remove microcracks from bone matrix, so avoiding an excessive accumulation of the latter, which could result in macroscopic failure (Parfitt, 2002). The accumulation of microcracks in a particular volume of material is addressed here using a Continuum Damage Mechanics approach (Lemaitre and Chaboche, 1990). This 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 Eq. 18. This 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: C = 1 − d C 0 (18) where C and C ₀ are, respectively, the stiffness tensors of damaged and undamaged bone (Lemaitre and Chaboche, 1990). In the isotropic damage theory, Eq. 18 can be rewritten in terms of the respective Young’s moduli, E and E ₀, as E = (1 − d) E ₀ (Pattin et al., 1996; Zioupos and Currey, 1998) (see Eq. 38).

A balance of microdamage is considered through the accumulation due to fatigue loading and the removal due to 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: d ̇ = d ̇ A − d ̇ R (19) where d ̇ A is the rate of damage accumulation by fatigue loading and d ̇ R is the rate of damage removal by bone remodelling. The latter is assessed by assuming that damage is uniformly distributed throughout the representative volume element (RVE). So, the amount of repaired damage is proportional to the damage present in that volume and to the volume of tissue being resorbed, V ̇ r , through the fraction that this volume represents within the bone matrix volume: d ̇ R = d V ̇ r V b m = d k res ⋅ O c a f b m (20)

Damage accumulation is evaluated following the procedure described in (Martínez-Reina et al., 2008; Martínez-Reina et al., 2009). This procedure makes use of the results of the experimental fatigue tests performed by Pattin et al. (Pattin et al., 1996), who provided the evolution of damage with the strain level and the number of cycles. This evolution was mathematically modelled by García-Aznar et al. (García-Aznar et al., 2005) to yield the following differential equation under tensile stresses: d ̇ a = N ̇ C 1 C 2 γ f 1 − d 1 − γ f ε m a x δ f exp − C 2 1 − d γ f (21) where N ̇ is the number of cycles applied per unit time and ɛ _max is the maximum principal strain expressed in μɛ. ¹ The rest of parameters and constants of the model are: C 1 = e C 2 − 1 K f C a ; δ f = 14.1 ; γ f = − 0.018 ε m a x − 4100 + 12 ; C 2 = − 20 ; (22) where K _f ([Ca]) is a function of the mineral content which will be defined next. The experimental tests perfomed by Pattin et al. (Pattin et al., 1996) included an estimation of fatigue life, N _f , which was related to the deformation by the following expression: N f = K f ε m a x δ f (23) where K _f was assumed constant and equal to 1.445 ⋅ 10⁵³ in tension. Martínez-Reina et al. (Martínez-Reina et al., 2008) introduced a correction in K _f to consider the degradation of the fatigue properties with the increase in mineral content. A life N _f = 10⁷ cycles was assigned to the fatigue limit, which is usually assumed to occur for a given fraction of the ultimate tensile strain, ɛ _u /β, where the parameter β depends on the type of material (Juvinall, 1967) and β = 2 was assumed for bone (Martínez-Reina et al., 2008) with good results. So, K _f was obtained from Eq. 23 as: K f C a = 1 0 7 ε u C a β δ f (24) where the ultimate tensile strain depends on the calcium concentration of bone matrix, [Ca], as Currey (Currey, 2004) showed. The following regression was fitted in (Martínez-Reina et al., 2008) from the experimental results presented by Currey (Currey, 2004): l o g ε u = 31.452 − 11.341 l o g C a (25) where ɛ _u is expressed in μɛ and the concentration [Ca] is expressed in mg of calcium per g of bone matrix. This concentration is directly related to the ash fraction, α, which will be defined in the next section. More precisely, the relation [Ca] = 398.8 α was assumed, based on the molecular weigths of hydroxyapatite and type I collagen (Martínez-Reina et al., 2008).

The mineralisation model used in this work is based on that presented in (Martínez-Reina et al., 2008) and implemented in a model which is similar to the present one (Martínez-Reina and Pivonka, 2019). That model was a mixture of differential and recursive equations, difficult to implement in a system of ODEs. For this reason, it has been simplified to yield an explicit set of differential equations, explained next.

Bone is made up of a solid bone matrix and pores filled with marrow. A certain representative volume element, V _RVE , can be divided into the bone matrix volume, V _bm , and the volume of pores, V _p . The bone matrix volume is divided into inorganic (mineral), organic (mainly collagen) and water phases, designated as V _m , V _o and V _w , respectively: V RV E = V b m + V p = V m + V o + V w + V p (26)

The composition of bone matrix is defined in terms of the volume fractions of the three phases as: v i = V i V b m ∀ i = m , o , w (27)

Thus, the following condition holds v m + v o + v w = 1 (28)

The mineral content is usually measured by the so-called ash fraction, the ratio between mass of mineral (or ash mass) and dry mass (the sum of inorganic and organic mass): α = m m m m + m o = ρ m V m ρ m V m + ρ o V o = ρ m v m ρ m v m + ρ o v o (29) where Eq. 27 have been used and ρ _i are the corresponding densities of the three phases, being ρ _m = 3.2 g/cm³ (Currey, 2004) and ρ _o = 1.41 g/cm³ (Hernandez et al., 2001a). The tissue density is then given by: ρ t = m w + m m + m o V b m = ρ w v w + ρ m v m + ρ o v o = 1 + ρ o − 1 v o + ρ m − 1 v m (30) where ρ _w = 1 g/cm³ and Eqs. 27, 28 have been used to derive the right-hand side of Eq. 30.

Osteoid, the tissue laid by osteoblasts, contains only the organic phase and no mineral. Mineral accumulates in bone matrix afterwards, during the mineralisation process, which consists of three phases: 1) an initial phase, called mineralisation lag time, that lasts from 6 to 22 days (Eriksen et al., 1990; Need et al., 2007) during which no deposition of mineral occurs; 2) a primary phase, which is very quick (it takes a few days to reach the 70% of the maximum mineral content (Hernandez et al., 2001b)), and 3) a secondary phase, when mineral is added at a decreasing rate (Parfitt, 1983), as the tissue becomes saturated with mineral. Mineral accumulates in bone matrix by displacing water (Hernandez et al., 2001a). Thus, the volume fraction of organic phase is approximately constant during the mineralisation process and fixed here at v _o = 3/7 (Martin, 1984); while the variations of mineral and water volume fractions would hold Δv _m = −Δv _w . So, the mineralisation process is accounted for through the temporal variation of v _m . As v m = V m V b m (recall Eq. 27), the temporal derivative of this expression gives: v ̇ m = V ̇ m V b m − v m V ̇ b m V b m = v ̇ m V b m = c o n s t a n t − v m V ̇ b m V b m (31)

The first term corresponds to the variation of the mineral content due to mineralisation or resorption and it would be equal to the variation of mineral content if V _bm were constant. The second term is due to the variation of porosity. To understand this term it must be noted that the mineral content decreases when osteoid is deposited, since osteoid does not contain any mineral and it only contributes to increase the bone matrix volume, so reducing the concentration of mineral. Both terms will be analysed separately: v ̇ m V b m = c t e = v ̇ m m i n e r a l i s a t i o n − v ̇ m r e s o r p t i o n (32)

In a previous model (Martínez-Reina et al., 2008; Martínez-Reina and Pivonka, 2019) the variation of v _m due to mineralisation was defined in a piecewise manner, considering the three phases separately: with no variation of v _m during the mineralisation lag time, a linear increase during the primary phase and an exponential increase during the secondary phase. Here, this procedure will be simplified by assuming that it is governed by a saturation model: v ̇ m m i n e r a l i s a t i o n = K v m m a x − v m (33) leading to an exponential solution which approximates rather well the global response of the previous model for a value of the constant K = 0.007, with a fast initial mineralisation rate (primary phase) that slows down in the mid-long-term (secondary phase). The maximum mineral content v m max = 0.516 was fixed such that, together with the aforementioned v _o = 3/7, it yields the maximum tissue density ρ t max = 2.31 g / c m 3 (Hernandez et al., 2001a) through Eq. 30. The rate constant K will be assumed fixed in this work, though it may depend on the amount of calcium and phosphorus available in the serum, which, in turn, may depend on diverse physiological factors (Peterson and Riggs, 2010).

The variation of mineral content due to resorption is similar to the damage repair term in the damage model (see Eq. 20). It must be taken into account that mineral is dissolved by osteoclasts and so removed from the bone matrix as damage was previously assumed. Thus, the amount of mineral removed by resorption is proportional to the mineral content of the tissue being resorbed, v _m , and to the proportion of tissue being resorbed within the bone matrix, V ̇ r V b m : v ̇ m r e s o r p t i o n = v m V ̇ r V b m = v m V ̇ r / V RV E V b m / V RV E = v m V ̇ r / V RV E f b m = v m f b m k res O c a (34)

Using Eqs. 32-34 Eq. 31 can be rewritten as: v ̇ m = K v m m a x − v m − v m f b m k res O c a − v m f ̇ b m f b m (35) where it has been used that V ̇ b m V b m = f ̇ b m f b m . Taking into account the balance between formation and resorption (Eq. 6): v ̇ m = K v m m a x − v m − v m f b m k form O b a (36)

Once v _m is updated, the ash fraction can be derived from Eq. 29 and the tissue density from Eq. 30. Then, the apparent density is given by: ρ = m w + m m + m o V RV E = m w + m m + m o V b m V b m V RV E = ρ t f b m (37)

Bone was assumed to be an isotropic material with a Poisson’s ratio ν = 0.3 and a Young’s modulus given in MPa by the following expresions: E ρ , d = 2014 ρ 2.5 1 − d i f ρ < 1.2 g / c m 3 1763 ρ 3.2 1 − d i f ρ ≥ 1.2 g / c m 3 (38)

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

The main effect of the treatment is the bone density gain (BDG), which is measured with respect to baseline, i.e. at the beginning of treatment, when bone apparent density is ρ ₀: B D G % = ρ t − ρ 0 ρ 0 ⋅ 100 (39)

The monthly time derivative of apparent density (MTDAD) is defined in Eq. 40 and used to measure the intensity of bone loss and, more precisely, the spike after Dmab discontinuation, as bone loss is accelerated a few months after the last Dmab injection, as compared to the disease state before treatment: M T D A D % m o n t h − 1 = ρ ̇ ρ t = 0 ⋅ 100 (40)

Finally, the fracture risk is evaluated by assessing the critical load σ ^crit (t), i.e. the overload that would lead to local fracture (d = 1) in a few days for a given bone. This overload would be given by the present model as a function of the porosity, mineral content, bone turnover rate and current damage level. So, for a given instant t, the overload is calculated by integrating Eq. 19 and using the following algorithm:

1. Extract the results f _bm (t), d(t), Oc _a (t) and α(t) for the considered instant t.

To compare the deterioration of bone mechanical properties with the disease and after discontinuation of the treatment, the critical load at each time point was compared to the critical load at the onset of disease giving the following safety factor: S F t = σ crit t σ crit t = 0 (42)

The new PK-PD model proposed for Dmab (see Figure 1) was adjusted first to obtain D ̃ Dmab B C and ζ. This was done by comparing the in-silico results of BCPM _ext with the clinical results obtained by Bone et al. (Bone et al., 2011). In the latter, the patients were subjected to a 2 years treatment of Dmab consisting in the usual dose approved by the WHO, 60 mg injections every 6 months (60Q6). The average time elapsed since menopause for those patients was around 10 years. These conditions were simulated with BCPM _ext for a combination of parameters. More precisely, in a 20 × 20 grid of equally spaced values defined in the range: D ̃ Dmab B C ∈ [ 0,10 ] d a y − 1 and ζ ∈ [0, 1]. The comparison was performed at two different sites: the lumbar spine, for which a representative porosity of 85% (f _bm = 15%) was chosen and a uniaxial compression state was assumed (σ = 0.15 MPa), and the hip, for which f _bm = 25% was chosen and a uniaxial tension state was simulated (σ = 0.5 MPa). The latter stress state would approximate the area of the femoral neck subjected to tensile stresses due to bending. The value of the stresses were adjusted to produce a homeostasis state for the corresponding value of f _bm , as done in Figure 4 of (Martínez-Reina et al., 2021b). The following error was computed from the comparison: E = ∑ i ∑ j B D G i j model − B D G i j clin 2 (43) where the superscripts model and clin refer to the in-silico and clinical results, respectively; the subscript i stands for lumbar spine and hip and the subscript j refers to the time points for which the results by Bone et al. (Bone et al., 2011) were available and could be compared (see Figure 2). This error was computed for all the points of the aforementioned grid and the subregion where the error was lower was refined to search for a global minimum of E, obtained for the pair of values ( D ̃ Dmab B C , ζ ) = ( 24.5 d a y − 1 , 0.85 ) , which we adopted for the rest of simulations. Figure 2 compares the in-silico and clinical BDG results corresponding to that pair of values. First of all, it can be seen that the predictions of BCPM _st (blue) are far from the clinical results as the parameters D ̃ Dmab B C and ζ were fitted for BCPM _ext . Thus, focusing on the results obtained with BCPM _ext it can be seen that BDG is greater in the lumbar spine than in the femoral neck and that bone loss is very acute after the discontinuation of the treatment. This rebound of bone loss is similar at both bone sites, though it appears to stabilise in the lumbar spine, as opposed to the femoral neck.

Next we investigate the cause of the BDG rebound behaviour, which is due to the underlying changes of bone cell numbers in the BCPM. Figure 3 shows the temporal evolution of osteoclasts and osteoblasts populations for the case of the hip during Dmab treatment, normalised by the cell numbers at the beginning of the treatment. The models BCPM _st (dashed) and BCPM _ext (solid) are compared.

Focusing on the results of BCPM _ext , one can see that the predicted population of osteoclasts falls with Dmab injections, as the drug suppresses the effect of RANKL, but as the concentration of Dmab decreases, a rebound in osteoclastic numbers is observed. When the treatment is discontinued, this rebound is prolonged and pre-treatment values are reached.

The osteoblasts population is also reduced and a rebound is observed as well. The changes in osteoblast numbers follow those of osteoclastic cells though slower and deferred in time. Consequently, the rebound is noticeable some months after stopping the Dmab treatment.

On the other hand, the behaviour predicted by the BCPM _st is significantly different: no rebound in osteoclastic cell numbers can be observed. On the contrary, a decrease after discontinuation of the treatment can be seen, and there are no evident changes in osteoblastic cell numbers.

Figure 4 compares further the predictions of the two models with respect to the evolution of apparent density (top left), MTDAD (top right), damage level (bottom left) and the safety factor (bottom right), which are presented for the case f _bm = 25% subjected to uniaxial tension and under a treatment 60Q6 beginning 10 years after menopause and discontinued after 5 years of treatment. It can be seen that menopause produces an immediate and pronounced bone loss (apparent density falls and MTDAD is negative) after menopause, when RANKL levels increase. Bone mass gain is evident after the Dmab injection, when the increase of RANKL levels is counteracted by the drug. Later, a series of cycles occur during the treatment and bone mass gain slows down as Dmab is cleared away from blood to lead again to bone loss at the end of each treatment cycle.

This global behaviour is seen both in BCPM _ext and BCPM _st , but there is a clear difference between them, which is more evident in MTDAD: bone loss predicted after discontinuation by BCPM _ext (black line) is more pronounced and prolonged than that predicted by BCPM _st (red line). The risk of treatment discontinuation is quite evident looking at the evolution of damage, which presents a peak some months after the last Dmab injection, a peak which is only predicted by BCPM _ext .

The safety factor is a variable that measures how far is the tissue to local failure and is mainly influenced by the damage level and the magnitude of strain. The latter is, in turn, influenced by the stiffness and thus by the apparent density. In summary, the safety factor is controlled by the apparent density and the damage level. However, in view of the results, it appears that the former has more influence, for such values of damage and SF shows a similar trend to apparent density though it decreases by 35% in 20 years while apparent density decreases only by 16% in the case shown in Figure 4.

Currently, Dmab is administered at a constant dose of 60 mg every 6 months. Consequently, dose is not adjusted for the patient’s body weight (BW) though it influences the concentration of drug present in the central compartment. Marathe et al. (Marathe et al., 2011) showed a strong influence of BW on Dmab efficacy. For this reason, a BW = 60 kg has been taken as the reference in the previous simulations. In the following, we analyse the effect of BW on the response of bone to the usual Dmab treatment 60Q6 (see Figure 5). It can be seen that the rebound in osteoclastic cell numbers seems more pronounced for lighter patients (see both apparent density and MTDAD plot, Figure 5). Furthermore, the increment in damage after discontinuation follows a similar trend. Nonetheless, bone density gain and consequently the safety factor, are clearly lower for heavier patients, for whom the dose of 60 mg appears to be insufficient.

An alternative for heavier patients could be an adjustment of the dosage, either with a greater dose or with a more frequent injection. The rationale for this option is based on the clinical results obtained by Marathe et al. (Marathe et al., 2011) on the temporal evolution of Dmab concentration in the central compartment. These authors observed a sharp drop around 5 months after the injection of 1 mg/kg of BW. This could justify the more frequent dosage in heavier patients, so to limit the bone loss rebound seen at the end of each injection cycle. Figure 6 compares the treatment 60Q5 for a patient of BW = 80 kg, which yields a quite similar response to that of the WHO-approved 60Q6 for BW = 60 kg.

Given that an abrupt withdrawal of Dmab treatment leads to rapid bone loss and increases bone fracture risk, an important question remains to be answered: “How can one safely discontinue Dmab treatment”? In this section we perform in-silico simulations to address this question.

One strategy would be to define a time period where Dmab is reduced to zero dose and then create a transition from the original dose (60 mg) to 0 mg by incrementally reducing the dose. Here, we simulate this transitional treatment by maintaining the frequency at 6 months and progressively decreasing the Dmab dose in the last 5 injections to 50, 40, 30, 20 and 10 mg. Figure 7 compares the 60Q6 treatment (black line) with the decreasing dose treatment (green line). It can be seen that there is no significant difference between the finally attained values for apparent bone density and damage. However, there is some difference on how quick bone is lost with the decreasing dose treatment showing a slower bone loss compared to an abrupt withdrawal. Similar, the damage variable increases more rapidly with the abrupt withdrawal compared with the decreasing dose strategy.

Several clinical studies have investigated an alternative Dmab discontinuation strategy which uses a bisphosphonate (BP) prior to or at the time for Dmab withdrawal (Kendler et al., 2020; Burckhardt et al., 2021; Tsourdi et al., 2021). BP have been reported to induce apoptosis of active osteoclasts (Halasy-Nagy et al., 2001; Takagi et al., 2021). Here we simulate the effect of BP treatment by increasing the apoptosis rate of active osteoclasts, A Oc a , by 10%.

Based on these clinical studies, two other alternative treatments were considered and compared in Figure 7: 1) the intake of BP 6 months after completing the 60Q6 constant dose Dmab treatment, as recommended by (Tsourdi et al., 2021) (Dmab + BP: Switch, red line); 2) the intake of BP 48 months after the Dmab treatment has commenced, or equivalently, 12 months after the Dmab dose was decreased following the pattern described above (i.e., 50, 40, 30, 20 and 10 mg for the last 5 injections) (Dmab + BP: Decr. dose, blue line).

It can be seen that the use of BP avoids the rapid bone loss and the increase of damage compared to the abrupt withdrawal of Dmab. Same values of apparent bone density and damage are obtained independent of the time of BP administration. The incremental reduction of Dmab does not appear to have a significant benefit compared to the instantaneous Dmab discontinuation.