Assessment of Strategies for Safe Drug Discontinuation and Transition of Denosumab Treatment in PMO—Insights From a Mechanistic PK/PD Model of Bone Turnover

Denosumab (Dmab) treatment against postmenopausal osteoporosis (PMO) has proven very efficient in increasing bone mineral density (BMD) and reducing the risk of bone fractures. However, concerns have been recently raised regarding safety when drug treatment is discontinued. Mechanistic pharmacokinetic-pharmacodynamic (PK-PD) models are the most sophisticated tools to develop patient specific drug treatments of PMO to restore bone mass. However, only a few PK-PD models have addressed the effect of Dmab drug holidays on changes in BMD. We showed that using a standard bone cell population model (BCPM) of bone remodelling it is not possible to account for the spike in osteoclast numbers observed after Dmab discontinuation. We show that inclusion of a variable osteoclast precursor pool in BCPMs is essential to predict the experimentally observed rapid rise in osteoclast numbers and the associated increases in bone resorption. This new model also showed that Dmab withdrawal leads to a rapid increase of damage in the bone matrix, which in turn decreases the local safety factor for fatigue failure. Our simulation results show that changes in BMD strongly depend on Dmab concentration in the central compartment. Consequently, bone weight (BW) might play an important factor in calculating effective Dmab doses. The currently clinically prescribed constant Dmab dose of 60 mg injected every 6 months is less effective in increasing BMD for patients with high BW (2.5% for 80 kg in contrast to 8% for 60 kg after 6 years of treatment). However, bone loss observed 24 months after Dmab withdrawal is less pronounced in patients with high BW (3.5% for 80kg and 8.5% for 60 kg). Finally, we studied how to safely discontinue Dmab treatment by exploring several transitional and combined drug treatment strategies. Our simulation results indicate that using transitional reduced Dmab doses are not effective in reducing rapid bone loss. However, we identify that use of a bisphosphonate (BP) is highly effective in avoiding rapid bone loss and increase in bone tissue damage compared to abrupt withdrawal of Dmab. Furthermore, the final values of BMD and damage were not sensitive to the time of administration of the BP.


Introduction
In this document we provide those details of the mathematical model not provided in the main document. More precisely, we explain here those features of the original model by Martin et al. [1] which remained unchanged; while those that have been modified or are more relevant to the study are 5 presented in the main document.
This document is structured as follows: In section 2 we present the general equations for competitive binding between ligands and receptors; the specific equations for the competitive binding of the complex Wnt-Scl-LRP5/6 are given in section 3 and the co-regulation of RANKL levels via the antagonistic 10 effect of PTH and nitric oxide in section 4. The equations for the competitive binding of the complex RANK-RANKL-OPG-Dmab were given in the main document. The upregulation of RANKL expressed by osteocytes due to damage is given in section 5. The equations related to the regulatory effect of TGF-β are given in section 6. Section 7 describes the term of proliferation 15 of osteoblast precursors. Finally, the model constants are given in Table 1.

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 20 that compete to bind to LRP5/6 to control the proliferation of osteoblast precursors and also the case of RANK, OPG and Dmab 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 25 A and B to form, respectively, the complexes A-R and B-R. Let us consider 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.

35
Following Pivonka et al. [2] 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. (3) are null. This condition in Eqs. (3a) and (3b) yield: where: The stationarity condition of Eqs. (3c)-(3e) yields: If the degradation (D X ,D X−Y ), production (P X ) and dissociation constants, (K X−Y ) are known, (6) constitutes a non-linear system of three equations with three unknowns, namely 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. (4) have been used. The stationarity condition is equivalent to 50 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. (3a), (3b) and (3e) by imposing that the stationarity condition is met, . In the case of 55 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 a way that ligand is not produced if its concentration reaches the maximum or saturation value, [L] max . As described by Pivonka et al.
[3] 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 70 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 binding Wnt-Scl-LRP5/6
The previous equations can be used to describe the 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 pro-80 tein LRP5/6, so triggering intracellular activation of β-catenin. Sclerostin, produced by osteocytes, 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. (7) 85 reads for LRP5/6: The production of sclerostin is given by an equation like (9), which now 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 90 is given by the receptor-ligand binding equation (3a), which reads here: The external dosage of sclerostin, P Scl,d , is set to zero and the endogenous 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 95 is downregulated by the mechanical stimulus through the repressor function π Ψ bm rep,Scl (see Eq. (36) later on). Replacing (19) and (18) into (17) yields: where the constants coefficients of the polynomial are: Some studies have shown an increase in serum sclerostin after menopause [4,5], while sclerostin expression (local mRNA levels) was found to decrease 110 in animal models of menopause [5]. Following [1] we have assumed an exponential decay of the degradation rate of sclerostin to acknowledge this 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. π PTH,NO act/rep,RANKL = λ s π PTH act,RANKL + π NO rep,RANKL + λ c π PTH act,RANKL · π NO rep,RANKL where the activator function accounting for the effect of PTH is: being K PTH act and K PTH rep constants and the concentration of PTH given by: which comes from Eqs. (8) to (11) when there is no ligand for the species, the 125 external dosage is null (P PTH,d = 0), the endogenous production rate is not regulated (π X act/rep,Y · Y = 1 see Eq. (11)) and the saturation value [PTH] max 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. (8) in the absence of ligands. The external 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 135 activator function that accounts for the production of NO by osteocytes. This function and the repressor function affecting the production of sclerostin by osteocytes (see Eq. (19)) 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 [6]. 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 Beaupré et al. [7]:

Upregulation of RANKL expressed by osteocytes due to microstructural damage
As proposed in [8] 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 155 osteoclasts. Its concentration is calculated following Pivonka et al. [2]: 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. (1)-(4) of the main document. These functions control the upregulation of the differentiation of Ob u into 160 Ob p , the upregulation of osteoclast apoptosis and the downregulation of the differentiation of Ob p into Ob a : with K TGF−β act and K TGF−β rep the activation and repression constants, respectively.

Proliferation of osteoblast precursors 165
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 170 Wnt-Scl-LRP5/6 signalling pathway through π Wnt act,Obp (see Eq. (42)). As discussed in Buenzli et al. [9], 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. [9] P Obp was defined considering a saturation 175 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. (44) is fulfilled.

Model constants
The model constants, except for those related to the damage and mineralisation algorithms, which were given in the main document, are provided in the following Units γ act 7 - Table 1: Values taken for the constants of the PK-PD model. All the constants were taken from the model developed by Martin et al. [1], except those marked with an asterisk.
All the constants in Table 1 were taken from the model developed by Martin et al. [1], except those marked with an asterisk. These changes were 185 motivated by the fact that Martin et al. did not considered damage and the variable mineral content of bone matrix. Hence, those constants needed to be readjusted to reproduce the behaviour of the previous model. So, β Scl,Ot and β NO,Ot were readjusted to achieve the same sclerostin and NO levels. The constants of the block "Upregulation of RANKL via damage" are new 190 since damage were not considered in the previous model. Their values were fitted as in [8], so that the contribution of damage to RANKL production is similar in normal situations (homeostasis) to the production of RANKL by osteoblast precursors (see Eq. (11) in the main document). The constants of the block "Dmab PK-PD constants" are specific of the Dmab PK model 195 and are taken from previous works [8,10]. The degradation rate of the complex RANKL-Dmab,D Dmab−RANKL , was chosen equal to the degradation rates of the other complex involving RANKL, while ζ andD Dmab BC were adjusted as explained in section 3.1 of the main document. The constants of the block "PMO related constants" were readjusted as done in [1], i.e. to 200 reproduce the results of the longitudinal experimental study of Nordin et al. [11] on the evolution of the bone mineral density (BMD) in the forearm of post-menopausal women. Martin et al. [1] did not account for variations in mineral content and assumed the bone matrix fraction to evolve in the same way as the BMD. Now that this limitation has been overcome by including 205 the mineral content, it was necessary to readjust the constants.