Mechanistic PK-PD model of alendronate treatment of postmenopausal osteoporosis predicts bone site-specific response

Alendronate is the most widely used drug for postmenopausal osteoporosis (PMO). It inhibits bone resorption, affecting osteoclasts. Pharmacokinetics (PK) and pharmacodynamics (PD) of alendronate have been widely studied, but few mathematical models exist to simulate its effect. In this work, we have developed a PK model for alendronate, valid for short- and long-term treatments, and a mechanistic PK-PD model for the treatment of PMO to predict bone density gain (BDG) at the hip and lumbar spine. According to our results, at least three compartments are required in the PK model to predict the effect of alendronate in both the short and long terms. Clinical data of a 2-year treatment of alendronate, reproduced by our PK-PD model, demonstrate that bone response is site specific (hip: 7% BDG, lumbar spine: 4% BDG). We identified that this BDG is mainly due to an increase in tissue mineralization and a decrease in porosity. The difference in BDG between sites is linked to the different loading and dependence of the released alendronate on the bone-specific surface and porosity. Osteoclast population diminishes quickly within the first month of alendronate treatment. Osteoblast population lags behind but also falls due to coupling of resorption and formation. Two dosing regimens were studied (70 mg weekly and 10 mg daily), and both showed very similar BDG evolution, indicating that alendronate accumulates quickly in bone and saturates. The proposed PK-PD model could provide a valuable tool to analyze the effect of alendronate and to design patient-specific treatments, including drug combinations.


Introduction
In this document we provide a complete version of the mathematical model used in this work (the code in MATLAB is also accessible from this publication). With the intention of making a self-contained explanation, the mathematical development followed in the main document is enriched with 5 the features which were not explained there. More precisely, we explain here the construction of the PK model (we also provide the differential equations governing the rest of the PK models which are not included in the main text), the bone cell population model based on the original model by Martin et al. [22] with the modifications introduced by Martínez-Reina et al. [25] and the 10 inclusion of the alendronate PD model. The figures of the sensitivity analysis are also included at the end of this document.
This document is structured as follows: in section 2, the differential equations governing the PK models studied in this work are presented. Section 3 explains the bone cell population model with all its components, and there- 15 fore is divided in several subsections. In particular, in subsection 3.1 we present the general equations for competitive binding between ligands and receptors; the specific equations for the competitive binding of the complex RANK-RANKL-OPG are given in subsection 3.2, those for the com-plex Wnt-Scl-LRP5/6 are given in subsection 3.3 and the co-regulation of 20 RANKL levels via the antagonistic effect of PTH and nitric oxide in subsection 3.4. The inclusion of damage in the model is explained in subsection 3.5. The upregulation of RANKL expressed by osteocytes due to damage is given in subsection 3.6. The equations related to the regulatory effect of TGF-β are included in subsection 3.7. Subsection 3.8 describes the term of pro- 25 liferation of osteoblast precursors. The algorithm of bone mineralisation is explained in subsection 3.9. In subsection 3.10, the pharmacodynamic model of the alendronate is introduced. In subsection 3.11, the model constants are given in Table 1. Finally, in section 4, the figures of the sensitivity analysis conducted with the model parameters are presented. 30

Pharmacokinetics
In this work, we investigated 5 different PK models in order to find the most suitable one, which minimizes error between experimental data and model results. The schemes of the models can be seen in figures 1 and 2 (the model in the latter figure is separated as it was the final selected model). In 35 the following, CC stands for the central compartment (plasma), NCT for the non-calcified tissues, BC for the bone compartment and IC for the inactive compartment, with the gut compartment added in the case of oral doses. The arrows in figures 1 and 2 indicate the flow of alendronate, k i are the absorption rate constants and k el,i the elimination rate constants, with i = 40 CC, NCT, BC, IC or urine. The studied models are: 1. One-compartment model with 2 elimination mechanisms, proposed by Chae et al. [5] (see Fig. 1a). This model has an important drawback. It is known from the literature that alendronate is only excreted by urine [38], but the model consists of a CC and 2 exits, one to urine 45 and another to a non-specified organ (in the original model, the outflow rates were termed k urine and k non−urine ). The latter was renamed here as k BC in order to have an analogous terminology in all models.
To overcome this drawback the following two-compartment model was proposed. 50 2. Two-compartment model. The only difference with the one-compartment model is that one of the exits from the CC flows into the BC and there is a return flow from BC to CC (see Fig. 1b).
3. Three-compartment model "in series". This model adds to the previous model a compartment connected to the BC, as seen in Fig. 1c.
The model was called "in series" to distinguish it from the threecompartment model "in parallel". This model considers the different availability of the alendronate deposited in bone. To this end, it distinguishes the drug that has been deposited near the bone matrixmarrow interface from the alendronate that was buried deeper into the 60 bone matrix. The former, termed as active alendronate (and contained in the BC), is more accessible for osteoclasts to resorb it, since resorption occurs mainly on the bone matrix -marrow interface. Thus, it can affect the osteoclastic activity through endocytosis. On the contrary, the latter is assumed to be inaccessible to resorption and was termed where V c is the volume of distribution of the central compartment.
The differential equations which govern the temporal evolution of the alendronate in the one-compartment model with two elimination mechanisms 90 (see Fig. 1a) are: The differential equations which govern the temporal evolution of the alendronate in the two-compartment model (see Fig. 1b) are: The differential equations which govern the temporal evolution of the alendronate in the three-compartment model "in series" (see Fig. 1c) are: The differential equations governing the temporal evolution of alendronate in the three-compartment PK model "in parallel" (see Fig. 2) are: The differential equations which govern the temporal evolution of the alendronate in the four-compartment model (see Fig. 1d in the main document) are:

Bone cell population model
A previously published mathematical BCPM describing the bone cell interactions was used [22]. This model considers catabolic (RANK-RANKL-OPG) and anabolic (Wnt-Scl-LRP5/6) signalling pathways, together with the action of parathyroid hormone (PTH), nitric oxide (NO), transforming growth 105 factor beta (TGF-β) and mechanobiological feedback on bone cells. The effect of bone mineralisation was added following Martínez-Reina and Pivonka [29] and Martínez-Reina et al. [24]. The accumulation and repair of microstructural damage was also taken into account as in Martínez-Reina et al. [27].

110
The bone cell types, whose concentrations are the state variables of the 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 osteoblasts (Ob u ) and osteoclasts (Oc u ) are assumed constant as in [22].
where D Obu , D Obp , D Ocu and D Ocp are the differentiation rates of Ob u , Ob p , Oc u and Oc p , respectively; A Oca is the apoptosis rate of Oc a and ∆ Oba is the rate of clearance of active osteoblasts through apoptosis or differentiation into osteocytes. The variables π TGF−β act,Obu , π TGF−β rep,Obp and π TGF−β act,Ocp represent activator and repressor functions related to the binding of TGF-β to its re-120 ceptor. Similarly, π RANKL act,Ocu and π RANKL act,Ocp are the activator functions related to the RANK-RANKL binding. Finally, P Obp is the proliferation rate of Ob p , a process which is mediated by the Wnt signalling pathway through the activator function π W nt act,Obp . These functions, as well as the remaining equations needed to complete the model, are described in the following subsections.

125
The values of the model constants are given in Table 1.
Finally, Eq. (29) establishes that the population of osteocytes varies as the bone matrix fraction f bm , if the density of osteocytes trapped within bone matrix, η, is assumed constant, as done in Martin et al. [22]. Bone matrix fraction is defined as the volume of bone matrix, V b , per total volume of the 130 bone sample (i.e. the representative volume element, V RV E ), expressed as a percentage, i.e.: Its evolution is obtained through the balance between resorbed and formed tissue: where k res and k f orm are, respectively, the rates of bone resorption and osteoid 135 formation.

Competitive binding
Many biological processes are controlled by binding of biochemical factors which act as receptor and ligand. In some of them two or more ligands compete to bind to the receptor. This is the case of Wnt and sclerostin 140 that compete to bind to LRP5/6 to control the proliferation of osteoblast precursors and also the case of RANK and OPG which compete to bind to RANKL to control the differentiation of osteoclast precursors into mature active osteoclasts.
Let us consider separately the binding of a given receptor R to its ligands for each species X=A,B,R a production term P X and a degradation term D X , along with a degradation term for the complex D X−Y . Let K r X−Y and K f X−Y be the reverse and forward binding reaction constants, respectively.
The law of mass action provides the following set of differential equations: where [X] and [X − Y] represent, respectively, the concentration of species X and complex X-Y. The degradation terms are assumed proportional to the concentration of the species, i.e. D X =D X [X], withD X being the degradation rate.
Following Pivonka et al. [37] we assume that the binding reactions are much faster than the cell responses they produce and hence a quasi-steady state can be assumed, implying that the time derivatives of Eqs. (34) are null. This condition in Eqs. (34a) and (34b) yield: where: The stationarity condition of Eqs. (34c)-(34e) yields: If the degradation (D X ,D X−Y ), production (P X ) and dissociation con- The total concentration of a receptor is the sum of the concentrations of the receptor which is found free and bound to ligands, i.e.: where Eqs. (35) have been used. The stationarity condition is equivalent to 170 establish that the production rate of a species must equal the degradation rate, including all its forms, free and bound. For instance, in the case of the receptor that can bind to different ligands, this condition reads: which can be obtained from Eqs. (34a), (34b) and (34e) by imposing that the stationarity condition is met, In the 175 case of a ligand, that only binds to the receptor, that condition reads: The degradation rates are usually assumed as constants but the production rates are modelled in a more complex way. For instance, the production rate of ligand L can be split into a term corresponding to endogenous production, P L,b , and a term accounting for external dosage, P L,d : The endogenous production is sometimes modelled by the following equation: where Y is the concentration of the cell type Y producing L with a production rate β L,Y , regulated by the species X through the activator or repressor function π X act/rep,Y . The parenthesis establishes a saturation condition in such 185 a way that ligand is not produced if its concentration reaches the maximum or saturation value, [L] max . As described by Pivonka et al. [35] the activation of a certain biological process regulated by the formation of the complex L-R is given by the ratio between the receptors R occupied by ligands L and the total number of 190 receptors: Similarly, the repressor action of the binding is given by the complementary to one of the latter: In case of a single ligand the latter expressions yield the first-order Hill activator and repressor functions:

Competitive RANK-RANKL-OPG binding
The RANK-RANKL-OPG signalling pathway controls the differentiation of uncommitted osteoclast progenitors and osteoclasts maturation, respectively through π RAN KL act,Ocu and π RAN KL act,Ocp (see Eqs.10 and 11 in the main document). Thus, an imbalance in that pathway, such as that occurring after 200 menopause, may result in the development of osteoporosis. Following Martin et al. [22], the concentrations of OPG, RANK and RANKL are given by the following equations: whereD X andD 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 205 X-Y and N RANK Ocp is the number of RANK receptors per osteoclast precursor. P OPG is the production rate of OPG by active osteoblasts: where β OPG,Oba is the OPG production rate, π PTH rep,Oba 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 210 takes place. To evaluate P RANKL , the RANKL production rate of Eq. (49), we have assumed that RANKL is expressed by osteocytes and osteoblasts precursors, following experimental evidence [30,44] and then: where P PMO RANKL is the RANKL production due to PMO, which expression is given in the main document; β RANKL,Ot and β RANKL,Obp are the RANKL pro-215 duction rate of osteocytes and osteoblast precursors, respectively; π PTH,NO act/rep,RANKL 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) [22] and π dam act,RANKL is an activator function accounting for the upregulation of RANKL expression by osteocytes due to microstructural 220 damage [27] (see the details of both regulatory functions in sections 3.4 and 3.6 respectively). 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: with X = Oc u , Oc p (53) and different constants: K RANKL act,Ocu and K RANKL act,Ocp , as in Martínez-Reina et al. [25]. In previous works [29,26,27] both constants K RANKL act,X were chosen equal, so resulting in a constant pool of Oc p in the steady state.

Competitive binding Wnt-Scl-LRP5/6
Wnt signaling is an anabolic pathway promoting the proliferation of osteoblast precursors and hence bone formation. Extracellular Wnt binds to Frizzled and the lipoprotein receptor-related protein LRP5/6, so triggering intracellular activation of β-catenin. Sclerostin, produced by osteocytes, 235 modulates the signaling pathway by its interaction with LRP5/6 receptors. This prevents the formation of the Wnt-Frizzled-LRP5/6 complex and therefore hinders preosteoblasts proliferation. Competitive Wnt-Scl-LRP5/6 binding is modelled as follows. First, Eq. (38) reads for LRP5/6: The production of sclerostin is given by an equation like (40), which now 240 reads: whereD Scl andD Scl−LRP5/6 are the degradation rates of sclerostin and the sclerostin-LRP5/6 complex, respectively. The concentration of this complex is given by the receptor-ligand binding equation (34a), which reads here: The external dosage of sclerostin, P Scl,d , is set to zero and the endogenous 245 production of sclerostin by osteocytes is: where β Scl,Ot and [Scl] max are, respectively, the sclerostin production rate and its maximum concentration. The production of sclerostin by osteocytes is downregulated by the mechanical stimulus through the repressor function π Ψ bm rep,Scl (see Eq. (74) later on). Replacing (57) and (56) into (55) yields: Following Martin et al. [22] we assumed that the total number of LRP5/6 receptors per osteoblast precursor (N LRP5/6 OBp ) is constant and thus: Solving for [LRP5/6] in Eq. (54) and replacing it into (58) yields the following second-order polynomial of the free sclerostin, [Scl]: where the constants coefficients of the polynomial are: We can use (59) in the previous expression together with the previously calculated cell populations and [Wnt], which is assumed constant, to work out the three coefficients. Only one solution of (60) is positive as shown by Martin et al. [22] and this solution [Scl] is then used in (54) to calculate [LRP5/6]. Finally, using Eqs. (43) and (59), the activator function in the 260 Ob p proliferation term can be calculated as: Some studies have shown an increase in serum sclerostin after menopause [1,16], while sclerostin expression (local mRNA levels) was found to decrease in animal models of menopause [16]. Following [22] we have assumed an exponential decay of the degradation rate of sclerostin to acknowledge this 265 discrepancy between the serum levels and the local expression of sclerostin: where t onset is the time of onset of menopause,D 0 Scl is the pre-menopause value of the degradation rate of sclerostin and τ P M O is a time constant.

Co-regulation of RANKL via PTH and NO concentration
RANKL transcription is upregulated by parathyroid hormone (PTH) and 270 downregulated by nitric oxide (NO). In the model developed by Martin et al. [22] this antagonistic influence was merged into a co-regulatory function capturing both effects.
where the activator function accounting for the effect of PTH is: and the repressor effect on OPG (see Eq.(50)) is accounted for through the 275 function: being K PTH act and K PTH rep constants and the concentration of PTH given by: which comes from Eqs. (39) to (42) when there is no ligand for the species, the external dosage is null (P PTH,d = 0), the endogenous production rate is not regulated (π X act/rep,Y · Y = 1 see Eq. (42)) and the saturation value [PTH] max 280 is large enough to assume the parenthesis equal to 1. β PTH andD PTH are the endogenous production and degradation rate of PTH, respectively. On the other hand, the factor corresponding to nitric oxide is: with K NO rep a constant and the concentration of NO given by: which also comes from Eq. (39) in the absence of ligands. The external 285 dosage of nitric oxide P NO,d is set to zero in this study, β NO,Ot ,D NO and [NO] max are, respectively, the endogenous production and degradation rate of nitric oxide and its maximum content. The factor π Ψ bm act,NO is the mechanical feedback activator function that accounts for the production of NO by osteocytes. This function and the repressor function affecting the production of 290 sclerostin by osteocytes (see Eq. (57)) are defined by the following sigmoidal functions: where ρ ∼ and α ∼ are, respectively, the minimum and maximum anticipated response, γ ∼ is the sigmoidicity, influencing the steepness of the response, and δ ∼ is the value of the stimulus producing the half-maximal response [34].

295
Finally, the strain energy density (SED) at the bone matrix level is given by the SED at the continuum level, Ψ = 1 2 σ : ε, and the bone matrix fraction through the following expression proposed by Beaupre et al. [2]:

Damage
Targeted bone remodelling theories hypothesise that one of the major 300 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 [32]. The accumulation of microcracks in a particular volume of material is addressed here using a Continuum Damage Mechanics approach [20]. This theory introduces a damage variable, d, which is linked 305 to the density of microcracks in a volume of material and to the loss of stiffness through Eq. (76). 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: where C and C 0 are, respectively, the stiffness tensors of damaged and undamaged bone [20]. In the isotropic damage theory, Eq. (76) can be rewritten in terms of the respective Young's moduli, E and E 0 , as E = (1 − d) E 0 [46, 33] (see Eq. (104)). A balance of microdamage is considered through the accumulation due to fatigue loading and the removal due to bone remodelling, as osteoclasts resorb 315 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 320 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:ḋ Damage accumulation is evaluated following the procedure described in 325 [24,28]. This procedure makes use of the results of the experimental fatigue tests performed by Pattin et al. [33], 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. [11] to yield the following differential equation under tensile stresses: whereṄ is the number of cycles applied per unit time and ε max is the maximum principal strain expressed in µε. 1 The rest of parameters and constants of the model are: where K f ([Ca]) is a function of the mineral content which will be defined next. The experimental tests perfomed by Pattin et al. [33] included an estimation of fatigue life, N f , which was related to the deformation by the following expression: where K f was assumed constant and equal to 1.445·10 53 in tension. Martínez-
[24] 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 7 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 [17] and β = 2 was assumed for bone [24] 340 with good results. So, K f was obtained from Eq. (81) as: where the ultimate tensile strain depends on the calcium concentration of bone matrix, [Ca], as Currey [7] showed. The following regression was fitted in [24] from the experimental results presented by Currey [7]: where ε u is expressed in µε and the concentration [Ca] is expressed in mg of 345 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 [24].
3.6. Upregulation of RANKL expressed by osteocytes due to microstructural damage As proposed in [27] we have assumed that RANKL expression by osteocytes is upregulated by the presence of microstructural damage in the bone matrix through the factor π dam act,RANKL , which is defined as a sigmoidal function of damage, d: where K dam is a constant, ρ dam is the minimum value of the factor π dam act,RANKL , corresponding to d = 0, while α dam is its maximum value, corresponding to d = 1.

Regulatory role of TGF-β
TGF-β is stored in the bone matrix and released during resorption by 360 osteoclasts. Its concentration is calculated following Pivonka et al. [37]: 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 Eqs.

Proliferation of osteoblast precursors
Let us recall the differential equation of osteoblast precursors.
We can rewrite this equation as: where the terms in the right-hand side correspond, respectively, to the differentiation of Ob u into Ob p , the differentiation of Ob p into Ob a and the proliferation of Ob p at a rate P Obp which is determined by P Obp and the 375 Wnt-Scl-LRP5/6 signalling pathway through π Wnt act,Obp (see Eq. (88)). As discussed in Buenzli et al. [3], a necessary condition for the Ob p population to stay bounded and to converge to a meaningful steady-state (with finite, positive cell densities) is that: Following Buenzli et al. [3] P Obp was defined considering a saturation 380 factor: where P 0 Obp is a constant and Ob sat p is the maximum concentration of osteoblast precursors above which no proliferation occurs. This saturation and the choice of P 0 Obp (see Table 1) ensures that Eq. (90) is fulfilled.

Algorithm of bone mineralisation 385
The mineralisation model used in this work is based on that presented in [24] and implemented in a model which is similar to the present one [29].
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.

390
Bone is made up of a solid bone matrix and pores filled with marrow. A certain representative volume element, V RV E , can be divided into the bone matrix volume, V bm , and the volume of vascular pores, V vas . 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: The composition of bone matrix is defined in terms of the volume fractions of the three phases as: Thus, the following condition holds 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 400 inorganic and organic mass): where Eq. (93) have been used and ρ i are the corresponding densities of the three phases, being ρ m = 3.2 g/cm 3 [7] and ρ o = 1.41 g/cm 3 [14]. The tissue density is then given by: where ρ w = 1 g/cm 3 and Eqs. (93) and (94) have been used to derive the right-hand side of Eq. (96). 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: (i) an initial phase, called mineralisation lag time, that lasts from 6 to 22 days [8,31] during 410 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 [13]), and (iii) 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 [14]. Thus, the volume 415 fraction of organic phase is approximately constant during the mineralisation process and fixed here at v o = 3/7 [23]; 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 = Vm V bm (recall Eq. (93)), the temporal derivative of this expression gives: 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 min-425 eral and it only contributes to increase the bone matrix volume, so reducing the concentration of mineral. Both terms will be analysed separately: In a previous model [24,29] 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 430 the primary phase and an exponential increase during the secondary phase. Here, as in Martínez-Reina et al. [25], this procedure will be simplified by assuming that it is governed by a saturation model: 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 435 fast initial mineralisation rate (primary phase) that slows down in the midlong-term (secondary phase). The maximum mineral content v max m = 0.516 was fixed such that, together with the aforementioned v o = 3/7, it yields the maximum tissue density ρ max t = 2.31 g/cm 3 [14] through Eq. (96). The rate constant K will be assumed fixed in this work, though it may depend on the 440 amount of calcium and phosphorus available in the serum, which, in turn, may depend on diverse physiological factors [34].
The variation of mineral content due to resorption is similar to the damage repair term in the damage model (see Eq. (78)). It must be taken into account that mineral is dissolved by osteoclasts and so removed from the 445 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 bm : Using Eqs. (98)-(100) Eq. (97) can be rewritten as: where it has been used thatV bm V bm =ḟ bm f bm . Taking into account the balance between formation and resorption (Eq.14 in the main document): Once v m is updated, the ash fraction can be derived from Eq. (95) and the tissue density from Eq. (96). Then, the apparent density is given by: Bone was assumed to be an isotropic material with a Poisson's ratio 455 ν = 0.3 and a Young's modulus given in MPa by the following expresions: These are based on the correlations experimentally obtained by Jacobs [15], which were multiplied by the factor (1 − d) to consider microstructural damage as usually done in Continuum Damage Mechanics [20].

Pharmacodynamic model of alendronate 460
This BCPM was integrated into the PK model by including it in the bone compartment (BC, see Fig. 2). All the absorption rate constants k i and the elimination rate constants k el,i fitted for the PK model (see Fig. 2) were assumed constant, except for k el,BC , which is variable and depends on the osteoclast-mediated release of alendronate from bone matrix. Note that the 465 latter variable depends on bone turnover [45,9], whereas in the PK model of section 2 k el,BC was constant. Now, in the PK-PD model k el,BC is given by: We have assumed that it is the concentration of alendronate in the bone compartment, BC, that controls the action of the drug on the bone remodelling process: where V Bone is the total volume of bone tissue in the body, 2.23 · 10 −3 m 3 , adapted from [43] for a 60kg adult female. The amount of alendronate accumulated within the bone increases with each weekly dose and therefore, if the effect of the drug were proportional to that amount, such effect would be increasingly pronounced. But this is not observed in the clinical trials [39], where the increase of BMD is high during the first months after the beginning of the treatment and slows down in the subsequent months, both in the lumbar spine and the hip. This behaviour could be explained by several facts highlighted in the literature: Alendronate is preferentially deposited in areas undergoing active re-480 sorption [38].
80% of normal bone turnover occurs in trabecular bone and 20% in cortical bone [10].
A considerably higher amount of alendronate is deposited in trabecular bone compared to cortical bone [21,19].

485
The volume of cortical bone is much higher than that of trabecular bone [10,43].
Therefore, a higher proportion of alendronate would be initially deposited in trabecular bone, which explains the fast increase of BMD in the lumbar spine and the hip observed in the clinical results. However, due to the fact 490 that the volume of trabecular bone is much lower than that of cortical bone, in a long treatment it is expected that the proportion of alendronate that reaches the former will be reduced over time in favour of the latter, due to the saturation of the tissue. This would explain the long-term stabilisation of bone mass gain observed clinically.

495
In view of the clinical results, we hypothesized that the BC can be divided in two parts: one part termed active, which is the closest to the bone-marrow interface and where the retrieval of alendronate is immediate through bone resorption; the second part, termed inactive, corresponds to the innermost regions of bone, where alendronate is buried and, to some extent, inaccessible 500 to bone resorption [40,38]. An updated compartmental model is shown in Fig. 3. The active subcompartment would predominate in trabecular bone, where most of the tissue is superficial, and its proportion would decrease with porosity, as the tissue becomes more compact.
We need a variable to measure this compacity, i.e. the proportion of 505 superficial tissue. The specific surface, S v , is defined as the area of bone matrix-marrow interface, S i , per total volume of the bone sample, V T . The following correlation between S v and vascular porosity, f vas , was given by [23] and used by [36]: where the vascular porosity is f vas = 1− f bm 100 and S v is expressed in mm 2 /mm 3 , 510 being this variable equal to 0 both for f vas = 0 and f vas = 1, when no free surface exists. In order to measure the compacity of the tissue, we would need to express that specific surface per bone matrix volume instead of per total volume: This quotient is not defined for f bm = 0, but Eq.(108) is not needed for 515 very low values of f bm , as we can assume that all the tissue is superficial and accessible to bone resorption in a very porous bone. More precisely, we have established f bm0 = 5% as the minimum bone matrix below which all the tissue is considered superficial. We have normalized the quotient in Eq.(108) and used the following expression for the active alendronate, i.e. the amount 520 of alendronate contained in the active subcompartment: where S v0 is the specific area corresponding to the reference value f bm0 . The concentration of active alendronate can be assessed as in Eq.(106), i.e.
[Ale BC,act ] = Ale BC,act V bone . The function Sv f bm · f bm0 S v0 , that controls the division between the active and inactive parts of the BC, is plotted against f bm in Fig.   525 4. If the tissue has no pores (f bm = 100%), all the BC is inactive as there is no surface where osteoclasts can resorb bone and, thus, there is no release of alendronate from bone matrix. As porosity increases (f bm decreases), the ratio between free surface and bone volume rises, as does the exposure of the drug (the active subcompartment becomes predominant). In contrast, 530 if more tissue is formed than resorbed, f bm increases and the inactive part becomes predominant, as the alendronate that was on the surface is buried into the bone matrix and is thus less accessible to osteoclasts.
So, Eqs. (14) and (15) are rewritten as follows, by changing the second terms on the right-hand side of both: turnover changes locally for any reason it will change accordingly in the whole skeleton. The adopted value f aver bm = 43.7% was estimated by assuming an average f bm = 93% for cortical bone [4], f bm = 14% for trabecular one [42] and that 80% of bone tissue volume is cortical and the rest is trabecular [43].
The active subcompartment is predominant in trabecular bone whose 550 volume is considerably lower than that of cortical bone. Thus, it is plausible to assume that this subcompartment will become saturated and relatively soon, whereas the total alendronate in the BC will not, as indicated in the previous subsection, because there is a large volume of cortical bone in the body able to admit high doses of the drug. It could also become saturated 555 in cases of very long treatments and then the excess of alendronate would be eliminated via urine, but the authors found no information on urine excretion in long treatments. The active alendronate has been saturated also as a function of S v /f bm , being its maximum concentration: with f a constant to be fitted. Thus, as alendronate enters the BC (through 560 Eq.(111)), it is divided between the active and inactive subcompartments following Eq.(109) 2 until [Ale BC,act ] = [Ale BC,act ] max , when the active subcompartment becomes saturated and no more drug is allowed in it. Alendronate has the following two effects on bone cells. On the one hand, it causes the disappearance of clear zones and ruffled borders disrupting 565 the cytoskeleton of osteoclasts by inhibiting farnesyl pyrophosphate (FPP) synthase, which leads to these structural changes and loss of function [41,12]. In other words, it limits the resorbing capacity of osteoclasts, which is measured in the model by k res (see Eq.(31)). Consequently, this parameter is reduced as follows: where k res,nom is the nominal value of k res , Π Ale rep is a constant quantifying the effect of alendronate on the resorbing capacity of osteoclasts and k el,BC · [Ale BC,act ] is the amount of alendronate per unit volume released from bone matrix through resorption, i.e. the concentration of alendronate that affects the surrounding osteoclasts. If Eq.(105) is replaced in (113): from which k res can be worked out: The second effect caused by alendronate is the inhibition of the FPP synthase in the mevalonate pathway and the reduction of protein prenylation, an essential post-translational lipid modification required for the function of numerous proteins, thereby inducing apoptosis in osteoclasts [45,41,12,40].

580
To account for this effect the apoptosis rate was increased in the model as follows: where Π Ale act quantifies the effect of alendronate on the apoptosis of osteoclasts. If Eq.(105) is replaced in (116):

Model constants
The model constants, except for those related to the damage and mineralisation algorithms, which were given in their respective sections, and for those adjusted in this work, which were given in the main document in the result section, are provided in the following In this section, the figures corresponding to the variation of each constant by 10% are presented.

+10%
-10% E Figure 12: Sensitivity analysis for a ±10% variation of Π Ale act . A: Alendronate plasma concentration vs time; B: Alendronate short-term urinary excretion vs time; C: Alendronate long-term urinary excretion vs time; D: BDG vs time for hip, i.e. f bm = 15% and E: BDG vs time for lumbar vertebra, i.e. f bm = 25%, for a once-weekly 70 mg dose. In images A, B and C the three curves are practically overlapping.

+10%
-10% E Figure 13: Sensitivity analysis for a ±10% variation of Π Ale rep . A: Alendronate plasma concentration vs time; B: Alendronate short-term urinary excretion vs time; C: Alendronate long-term urinary excretion vs time; D: BDG vs time for hip, i.e. f bm = 15% and E: BDG vs time for lumbar vertebra, i.e. f bm = 25%, for a once-weekly 70 mg dose. In images A, B and C the three curves are practically overlapping.