Mathematics121

International Journal for Numerical Methods in Biomedical Engineering
Wiley-Blackwell, John Wiley & Sons
Mathematical modeling of postmenopausal osteoporosis and its treatment by the anti-catabolic drug denosumab
S Scheiner, P Pivonka, […], and C Hellmich

Additional article information

Abstract
Denosumab, a fully human monoclonal antibody, has been approved for the treatment of postmenopausal osteoporosis. The therapeutic effect of denosumab rests on its ability to inhibit osteoclast differentiation. Here, we present a computational approach on the basis of coupling a pharmacokinetics model of denosumab with a pharmacodynamics model for quantifying the effect of denosumab on bone remodeling. The pharmacodynamics model comprises an integrated systems biology-continuum micromechanics approach, including a bone cell population model, considering the governing biochemical factors of bone remodeling (including the action of denosumab), and a multiscale micromechanics-based bone mechanics model, for implementing the mechanobiology of bone remodeling in our model. Numerical studies of postmenopausal osteoporosis show that denosumab suppresses osteoclast differentiation, thus strongly curtailing bone resorption. Simulation results also suggest that denosumab may trigger a short-term bone volume gain, which is, however, followed by constant or decreasing bone volume. This evolution is accompanied by a dramatic decrease of the bone turnover rate by more than one order of magnitude. The latter proposes dominant occurrence of secondary mineralization (which is not anymore impeded through cellular activity), leading to higher mineral concentration per bone volume. This explains the overall higher bone mineral density observed in denosumab-related clinical studies. Copyright © 2013 John Wiley & Sons, Ltd.

Keywords: bone remodeling, pharmacokinetics, systems biology, micromechanics
1 INTRODUCTION
Osteoporosis is a disease characterized by significant bone loss, leading to an increased risk of bone fracture. It is particularly common in postmenopausal women and the elderly. Often postmenopausal osteoporosis (PMO) is first apparent from a bone density scan, and for this reason it is sometimes called the “silent disease”, as there may be no noticeable symptoms associated with bone loss 1. Once detected, various ways to limit further bone loss or induce bone gain have been tried clinically. These include dietary changes and exercise, and different drug treatments that influence the turnover of bone, including administration of bisphosphonates and parathyroid hormone (PTH) 2–4.

Recently, Amgen 5 developed a monoclonal antibody named denosumab for the treatment of osteoporosis. Denosumab binds with high affinity and specificity to the molecule known as RANKL, the ligand of the receptor activator of nuclear factor kappa-B 6, which is essential for the development of active osteoclasts – osteoclasts are the cells responsible for bone resorption 7. Following the successful completion of more than 30 clinical trials involving a large number of patients 8–12, denosumab is now available for the treatment of osteoporosis and is marketed for this purpose under the name Prolia 13, whereas denosumab is also marketed under the name Xgeva for the purpose of treating of bone cancer metastases.

Anticipating the progression of PMO and the effects of therapeutic agents (such as denosumab) on the basis of experimental observations alone is difficult (if not impossible), owing to complex interactions between the cells governing the bone physiology in health and disease. Hence, experiment-based drug design is a lengthy and expensive process. Given the continuously increasing capability of computer simulations, supporting and complementing the design of drug administration regimes by computational modeling has proven promising 14–16. The gold standard of mathematical modeling of disease systems comprises coupling pharmacokinetic (PK) models (quantifying the “availability” of the drug at sites of interest) with pharmacodynamic (PD) models (quantifying the time-dependent effect of the drug on the overall physiological system). Most of previous PK/PD-type models relate, on a simplistic phenomenological basis, drug exposure to related levels of so-called biomarkers, which serve as indicators for drug effectiveness and safety 17–19. While state-of-the-art approaches, see for example the works of Marathe et al. 20, Pivonka et al. 21, or Peterson and Riggs 22, provide valuable qualitative insights on the effects of drug administration, two major shortcomings are evident in any of the state-of-the-art models we found in the literature:

The studied disease is usually driven through a constant change of certain parameters. However, as pointed out by Post et al. 19, it is essential to consider the experimentally observed time-dependence of the biological factors triggering progressive, chronic diseases such as PMO.
Bone remodeling is driven solely biochemically, without consideration of biomechanical feedback. It is, however, beyond any doubt 23 that bone remodeling is regulated both biochemically and biomechanically.
In this paper, we aim at resolving these shortcomings by developing and integrating computational models in an attempt to better understand the mechanism(s) of the action of denosumab on remodeling events in bone. We also take into consideration the clinically observed increase in bone density following treatment with denosumab. Remarkably, this observation is a counterintuitive clinical outcome as preventing bone resorption (which is the functional principle of anti-catabolic drugs such as denosumab) does not, by itself, imply any bone gain. In particular, we aim at explaining how the binding of RANKL by denosumab may be related to the clinically observed increase in bone density.

For these purposes, we first develop a PK model, on the basis of the experimental data of Bekker et al. 8, that predicts the blood serum concentration of denosumab following subcutaneous administration (Section 2). We then link this PK model to a recently developed integrated systems biology-continuum micromechanics model of bone remodeling 24. This model takes into account key features of the known physiology of bone remodeling, including the most significant biochemical and biomechanical regulatory mechanisms (Section 3.1). Here, the model of Scheiner et al. 24 is extended: (i) to enable adequate simulation of the effects of PMO on bone remodeling (Section 3.2), and (ii) to also provide the pharmacodynamical effects of the exogenous drug denosumab (Section 3.3). We utilize our model then for simulation of PMO without administration of denosumab (Section 4.1), and we study the effects of single administration and multiple administrations of denosumab (Sections 4.2 and 4.3), and interpret the results through comparison with experimental data from biomarker measurements (Section 5.1) and from bone mineral density measurements (Section 5.2). After discussing the influence of the mechanical loading on the disease progress (Section 5.3), we round the paper off with concluding remarks (Section 6).

2 PHARMACOKINETICS OF DENOSUMAB
Denosumab is injected subcutaneously and over time is absorbed into the blood serum. Through blood circulation, denosumab may then reach surfaces of metabolizing bone. Mathematical description of the effect of denosumab requires formulation of a PK model taking into account the delay between drug administration and effectivity.

To this end, in line with Marathe et al. 20, a two-compartment representation 25 of the physiological system is chosen, with compartment I relating to the subcutaneous tissue, in which denosumab is present at mass concentration mden,sub, and compartment II relating to the blood serum, in which denosumab is present at mass concentration mden,ser, see Figure 1. The transfer of denosumab from the subcutaneous tissue to the blood serum is governed by absorption rate kabs. We take into account that in the blood serum denosumab is subjected to degradation (also referred to as clearance) by both organs and enzymes. Degradation by organs is considered via degradation rate korg, whereas degradation by enzymes is considered via degradation rate kenz.

Figure 1
Figure 1
Pharmacokinetics model representation of denosumab: compartment I represents the subcutaneous tissue with denosumab being present at mass concentration mden,sub, whereas compartment II represents the blood serum with denosumb being present at mass concentration …
2.1 Governing equations for calculation of the denosumab concentrations
For development of a mathematical model describing the pharmacokinetics of denosumab, we assume that the absorption of denosumab from the subcutaneous tissue into the blood serum is linearly related to the actual mass (concentration) mden,sub in the subcutaneous tissue. Accordingly,

equation image 1
Solution of Equation (1) requires consideration of an appropriate initial condition. In a single-dose administration regime, with denosumab administered at t = 0, t being the time variable, Δmden,sub is the sudden increase of mden,sub in the subcutaneous tissue due to the administered dose Dden, and prior to denosumab administration the concentration of denosumab in the subcutaneous tissue equals zero. Based on this initial condition, Equation (1) can be solved analytically, see Equation (A.1) in Appendix A.

In the blood serum, denosumab degradation occurs proportional to its actual concentration mden,ser. However, experimental findings 26 suggest that above a certain limit concentration, hereafter denoted as An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu2.jpg, degradation is restricted. Following classical literature on enzyme kinetics 27, we include in our PK model that this degradation restriction can be attributed to a limited degradation capacity of the involved enzymes. In mathematical terms, the serum concentration of denosumab is hence governed by

equation image 2
Prior to administration of denosumab, at t = 0, the denosumab serum concentration equals zero, thus mden,ser(t = 0) = 0. Due to absorption of denosumab from the subcutaneous tissue mden,ser builds up over time, governed by Equation (2)1, and arrives at limit concentration An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu4.jpg at An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu5.jpg. For An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu6.jpg the evolution of mden,ser follows Equation (2)2 (because An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu7.jpg). After reaching a dose-dependent peak concentration, mden,ser decreases, because of the decreasing supply from the subcutaneous tissue, until it reaches again An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu8.jpg at An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu9.jpg. For An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu10.jpg the serum concentration is smaller than the limit concentration and hence again governed by Equation (2)1. On the basis of these conditions, Equation (2) can be solved for mden,ser, see Equations (A.2)–(A.4) in Appendix A.

Depending on the administration intervals and the delivered doses, a denosumab administration regime involving multiple administrations is characterized by a large number of possible initial and boundary conditions on the basis of which the Equations (1) and (2) have to be solved, see Appendix A.

2.2 Calibration of the absorption and degradation rates
Now, the PK model derived in Section 2.1 and Appendix A is calibrated through determination of parameters kabs, korg, and kenz. For this purpose, we use the clinical data of Bekker et al. 8, who measured the temporal evolutions of the serum concentration for six different administration doses An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu14.jpg mg/kg, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu15.jpg mg/kg, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu16.jpg mg/kg, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu17.jpg mg/kg, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu18.jpg mg/kg, and An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu19.jpg mg/kg. A two-membered evolution algorithm 28 is employed for determination of absorption and degradation rates, aiming at minimization of the relative error between computed and experimentally determined serum concentrations.

Dosing of denosumab occurs relative to the patient body mass. For conversion of the administered dose into a corresponding increase of mden,sub, we consider an average subject body weight of 70 kg 29, and a blood serum volume of 3 liters 30, resulting in An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu20.jpg. Because, for patient convenience, denosumab administration regimes preferably involve a small number of high-dose administrations rather than a large number of low-dose administrations, we consider only doses An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu21.jpg, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu22.jpg, and An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu23.jpg for model calibration, to focus on practically relevant doses. Furthermore, the experimental data of Bekker et al. 8 show that above a serum concentration of ≈ 750 ng/ml, denosumab degradation kinetics are restricted, thus An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu24.jpg ng/ml. On this basis, the implemented evolution algorithm yields kabs = 1.8498 month  − 1, korg = 0.9500 month  − 1, and kenz = 0.9932 month  − 1, see Figure 2 for a comparison of clinical data to PK model results.

Figure 2
Figure 2
Calibration of the denosumab pharmacokinetics model against the clinical data of Bekker et al. 8: (a) comparison between temporal evolutions of model-predicted and measured serum concentrations, and (b) model-predicted versus measured serum concentrations; …
3 THEORETICAL BASIS FOR MODELING THE EFFECT OF DENOSUMAB ON REMODELING OF BONE SUBJECTED TO POSTMENOPAUSAL OSTEOPOROSIS
3.1 A bio-chemo-mechanically coupled mathematical model of bone remodeling
The main novelty presented in this paper, namely a mathematical model enabling simulation of the effect of denosumab on bone remodeling within a representative volume element of bone subjected to PMO, is centered around a recently developed model integrating bone cell population kinetics and multiscale bone mechanics 24. For conciseness, we here present the model components as follows: While model extensions that are original contributions of this paper are presented at appropriate length, any model components borrowed from previous publications are only outlined briefly.

3.1.1 General remarks The previously developed and here utilized integrated approach to computational simulation of bone remodeling 24 is based on the notion of a representative volume element (RVE). Within an RVE of bone microstructure (with a characteristic length of 1 to 3 mm), bone remodeling, together with the underlying biochemistry and the mechanical environment, occurs (quasi-)homogeneously. Hence, spatial effects (such as transport processes and local variations of biochemical factors or mechanical properties) can be neglected.
3.1.2 Bone cell population kinetics The biochemically and biomechanically guided evolutions of osteoblasts (the cells responsible for bone formation) and osteoclasts (the cells responsible for bone resorption) are mathematically modeled by means of a so-called bone cell population model (BCPM), the foundation of which was laid by Lemaire et al. 31, explicitly considering various distinct developmental stages of osteoblasts and osteoclasts. Namely, the model takes into account the populations of uncommitted osteoblast progenitor cells (abbreviated to OBu), osteoblast precursor cells ( OBp), active osteoblasts ( OBa), osteoclast precursor cells ( OCp), and active osteoclasts ( OCa). The populations of these cell types are expressed in terms of respective molar concentrations Ci (describing the “amount” of the species per unit volume and being thus equivalent to corresponding average cell numbers). Progression of osteoblasts and osteoclasts along the chosen developmental stages ( OBu → OBp → OBa and OCp → OCa) is implemented considering the (regulatory) mechanisms described in the following paragraphs, see Figure 3 and 21,24,32, as well as references therein, for details on the mathematical implementation.
Figure 3
Figure 3
Graphical sketch, adapted from Pivonka et al. 33, showing all mechanisms considered in the model presented in this paper, with novel contributions colored green and red: guidance of cell developments (OBu → OBp → OB …
The concentration of the OBps, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu31.jpg, increases because of differentiation of OBus (which are assumed to exhibit a constant concentration), governed by maximum differentiation rate An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu32.jpg and activated through binding of transforming growth factor β (TGF- β) to its receptors 34,35, via activator function An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu33.jpg; it further increases because of proliferation of OBps, governed by maximum proliferation rate An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu34.jpg and activated through increased mechanical loading 36–39, by means of activator function An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu35.jpg, see Equation (7) and 24 for details; and it decreases because of differentiation of OBps to OBas, governed by maximum differentiation rate An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu36.jpg and inhibited by binding of TGF- β to its receptors 34,35, via repressor function An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu37.jpg:

equation image 3
The concentration of the OBas, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu39.jpg increases because of differentiation of OBps, governed by maximum differentiation rate An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu40.jpg and inhibited through binding of TGF- β to its receptors, via repressor function An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu41.jpg; and it decreases because of apoptosis governed by constant apoptosis rate An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu42.jpg:

equation image 4
The concentration of the OCas, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu44.jpg increases because of differentiation of OCps (which are assumed to exhibit a constant concentration), governed by maximum differentiation rate An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu45.jpg and activated through binding of receptor activator nuclear factor kappa β (RANK) to its ligand, RANKL, via activator function An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu46.jpg; and it decreases because of apoptosis, governed by maximum apoptosis rate An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu47.jpg and activated through binding of TGF- β to its receptors 40, via activator function An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu48.jpg:

equation image 5
Note that An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu50.jpg, further specified in Section 3.3, also considers the influence of the parathyroid hormone (PTH) and osteoprotegerin (OPG) on binding of RANK to RANKL 6, as well as the reduction of RANKL production because of increasing mechanical loading 41–44, see Equation (8) and 24,32 for details.

Activator functions An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu51.jpg, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu52.jpg, and An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu53.jpg, as well as repressor function An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu54.jpg are defined according to the concept of Hill functions 45, governed by concentrations of the respective substances, namely TGF- β and the RANK-RANKL complex. The complete formulations of these functions, as well as respective derivations, are given in full detail in 24,32.

3.1.3 Relation of bone cell populations to bone composition
The considered RVE of cortical bone is composed of extravascular bone matrix and vascular pore space. 1 Importantly, we consider the quasi-instantaneous character of primary mineralization 46–48 by having the modeled osteoblasts deposit directly mineralized solid bone matrix, thereby omitting consideration of the initially laid down unmineralized osteoid. For the mineralized bone matrix, we consider a constant, organ-dependent mineralization state, see Appendix B and 24,49 for further details. The “amount” of the aforementioned components within the studied RVE, namely extravascular bone matrix and vascular pore space, is quantified by means of volume fractions: the volume fraction of extravascular bone matrix is defined as its volume within the RVE divided by the total volume of the RVE, fbm = Vbm / VRVE; the volume fraction of vascular pore space, fvas, is defined analogously, fvas = Vvas / VRVE. A numerical value of fbm = 1 thus expresses that the RVE consists solely of extravascular bone matrix, in turn implying fvas = 0 (because fbm + fvas ≡ 1). Dealing with cortical bone only, the initial bone composition in all simulations presented in Section 4 is defined through An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu55.jpg. The change of volume fractions is governed by bone resorption and bone formation. To relate the bone cell concentrations An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu56.jpg and An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu57.jpg, obtained from the BCPM, to corresponding changes of the bone constituent volume fractions, we define resorption rate kres, quantifying how much bone is resorbed per concentration of active osteoclasts, and formation rate kform, quantifying how much bone is formed per concentration of active osteoblasts:

equation image 6
Accordingly, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu59.jpg represents equilibrated bone remodeling (with constant bone constituent volume fractions fbm and fvas), An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu60.jpg represents a catabolic bone remodeling regime (where fvas increases), while An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu61.jpg represents an anabolic bone remodeling regime (where fbm increases).

3.1.4 Mechanical regulation of bone remodeling In line with Scheiner et al. 24, we consider the strain energy density that is experienced by the extravascular bone matrix, Ψbm, as adequate mechanoregulatory quantity; in other words, the current magnitude of Ψbm dictates whether and to which extent addition or removal of bone is triggered by the mechanical loading. Estimation of Ψbm is accomplished by means of a continuum micromechanics-based 50 multiscale bone mechanics model 49: Based on the morphology of the RVE, on the intrinsic, experimentally accessible stiffnesses of the bone constituents (extravascular bone matrix and pore space), and on the corresponding volume fractions accessible through evaluation of the bone cell concentrations obtained from the BCPM, see Equation (6), the microscopic strain tensor of the extravascular bone matrix, Ψbm, is estimated. Ψbm follows then through An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu62.jpg, with An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu63.jpg denoting the microscopic stiffness tensor of the extravascular bone matrix (see Appendix B for the numerical values of the components of An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu64.jpg), and “ :” denoting a mathematical operation called second-order tensor contraction. Note that An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu65.jpg relates Ψbm to the microscopic stress tensor σbm through a linear constitutive relation, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu66.jpg.
Importantly, continuum micromechanics allows to “downscale” the stresses Σcort, acting upon the bone on the macroscale, to the related strain tensor Ψbm on the microscale. For the sake of simplicity, we prescribe, in all simulations presented in Section 4, a constant, uniaxial compression loading scenario, with Σ33 = − 30 MPa. 2 In our model, both mechanically governed activation of OBp proliferation and mechanically governed production of RANKL are regulated by Ψbm 24, establishing a feedback-type relation between bone biology and bone mechanics:

equation image 7
equation image 8
In Equation (7), An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu69.jpg denotes the value of An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu70.jpg related to “normal” mechanical loading, which, in turn, induces a corresponding SED An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu71.jpg, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu72.jpg; note that An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu73.jpg governs the osteoblast precursor proliferation in Equation (3). The sensitivity of OBp proliferation to increased mechanical loading can be adjusted by the so-called anabolic strength parameter λ. This parameter has been found to reasonably replicate mechanically stimulated bone gain if set to λ = 1.25 if An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu74.jpg, and λ = 0 if An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu75.jpg 24. In Equation (8), An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu76.jpg denotes the production of RANKL because of decreased mechanical loading (compared with normal loading conditions). In our model, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu77.jpg co-governs the concentration of RANKL, see Equation (12) and below, this way modulating the differentiation of OCp to OCa. Inhibition parameter κ allows for adjustment of the sensitivity of RANKL production to a reduced mechanical loading. Setting κ = 105 pM/day if An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu78.jpg and κ = 0 if An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu79.jpg has proven to give rise to reasonably simulated bone loss if subjected to disuse scenarios 24.

3.2 Consideration of the effect of postmenopausal osteoporosis on bone remodeling
Experimental studies on the pathophysiology of PMO reveal that several mechanisms might occur simultaneously, together inducing the disease pattern known as PMO. Among several other potential pathogenic mechanisms, estrogen deficiency 3 has been widely accepted as the main cause of PMO 51,52, resulting in increased osteoclast and osteoblast concentrations (and thus in an increased bone turnover) 53. Relative to the cell concentrations in healthy (premenopausal) bone, the increase of An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu80.jpg is higher than the increase of An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu81.jpg that leads to (potentially) significant bone loss after onset of PMO. This initial phase of PMO is accompanied by an increased ratio of the RANKL concentration over the OPG concentration, compared with normal bone physiology 54. After some time (some months up to a few years), the rate of bone loss decreases, and postmenopausal bone is subjected to long-lasting moderate bone loss thereafter.

Driving PMO biochemically can be achieved in different ways; see 21 for a related study. Estrogen and its role in bone remodeling is not considered explicitly in our model, so PMO cannot be simulated by feeding experimentally observed estrogen levels into the model. Instead, PMO is initiated by introducing disease-related increased production of RANKL, leading to increased osteoclast differentiation. In order to account for the moderate bone loss in the second phase of PMO, it is assumed that the excess production of RANKL reduces over time. In detail, we prescribe the following RANKL production regime for An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu82.jpg, with tPMO,ini being the point in time when PMO is initiated:

equation image 9
with An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu84.jpg as PMO-initiating excess production rate of RANKL, and with An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu85.jpg as a reduction factor of An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu86.jpg, which we define as

equation image 10
Parameter ξ and the characteristic time of the RANKL production decrease, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu88.jpg, together determine the shape of the Lorentz-type function given by Equation (10). In our model, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu89.jpg co-governs the concentration of RANKL, see Equation (12) and below, this way modulating the differentiation of OCps to OCas. Furthermore, in vivo data suggests that the decreasing mechanoresponsiveness of bone due to increased osteocyte apoptosis, another effect of estrogen deficiency, also plays a major role for the progression of PMO 52,56. In order to incorporate this effect into our model we assume that, after onset of PMO at t = tPMO,ini, the parameters governing the sensitivity of bone remodeling to a changing mechanical loading, κ and λ, see Equations (7) and (8), are decreased via reduction factor An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu90.jpg,

equation image 11
with τPMO as the characteristic time of the PMO-related mechanoresponsiveness decrease of cortical bone. Hence, contrary to conventional PK/PD approaches, in our model the disease-governing parameters are not held constant, but are varied over time, following the pathophysiological trends reported in the literature 19.

The set of Equations (9)–(11) can be considered as a (semi-)phenomenological approach to the sum of processes that in reality contribute to initiation and maintenance of PMO. While such an approach is common practice for simulation of complex physiological systems, the model parameters introduced in this section cannot be measured or directly deduced from experimental data. Thus, for the time being, we are content with adjusting these parameters such that simulation results agree well with corresponding experimental data, as demonstrated in Section 4.1.

3.3 Incorporation of the action of denosumab
In our model, denosumab contributes to regulation of osteoclast differentiation via activator function An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu92.jpg. Denosumab competes with RANK (and OPG) as for binding to RANKL, meaning the higher the concentration of denosumab the lower the concentration of RANKL-RANK complexes, hence the lower the numerical value of An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu93.jpg and the lower the fraction of the osteoclast precursor population that differentiates into active osteoclasts. Adapting the approach of Pivonka et al. 32, the concentration of RANKL follows

equation image 12
where Ka,[RANKL-OPG], Ka,[RANKL-RANK], and Ka,[RANKL-d] are the equilibrium association binding constants for binding of OPG, RANK, and denosumab to RANKL, COPG, CRANK, and Cden,ser are the molar concentrations of OPG, RANK, and denosumab, βRANKL is the intrinsic RANKL production rate, PRANKL is the RANKL dosage term, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu95.jpg is the constant degradation rate, and An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu96.jpg is the maximum concentration of RANKL. The concentration of RANKL-RANK complexes is considered through

equation image 13
While the BCPM and respective concentrations of regulatory factors are formulated for a specific RVE of bone tissue, the concentration of denosumab is given in the blood serum. In Equation (12), the varying accessibility of denosumab to different RVEs of bone tissue is taken into account through accessibility factor ζ: ζ = 1 represents unrestricted access to denosumab, whereas ζ < 1 reflects access restrictions (for example because of low or uneven blood circulation). In the following, we only consider the case of uniform and sufficient blood circulation, ζ = 1. Finally, activator function An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu98.jpg is governed by 32

equation image 14
with Kd,[RANKL-d] as the corresponding equilibrium dissociation binding constant. In the present model, the external production rate of RANKL, PRANKL, is made up of two components: An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu100.jpg (stemming from decreased mechanical loading), see Equation (8), and An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu101.jpg (for simulation of postmenopausal osteoporosis, see Equation (9), An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu102.jpg.

3.4 Governing model parameters
The interaction of RANKL and denosumab is considered via the denosumab serum concentration mden,ser, see Section 2, and the association binding constant Ka,[RANKL-d] which quantifies the affinity of denosumab-binding to human RANKL. Choice of the numerical value of Ka,[RANKL-d] is based on the in vitro experiments performed by Kostenuik et al. 57. Solution equilibrium binding analysis revealed an equilibrium dissociation binding constant of Kd,[RANKL-d] ≈ 3 pM. The equilibrium association binding constant is defined as the inverse of the corresponding dissociation constant, Ka,[RANKL-d] = (Kd,[RANKL-d]) − 1 ≈ 0.33 pM  − 1. Furthermore, the denosumab concentrations provided by the PK model are given in ng/mg, whereas the unit of the equilibrium association binding constant is pM  − 1 = (10 − 12 mol/l) − 1. Hence, PK model-derived mass concentrations mden,ser have to be converted to the corresponding molar concentration Cden,ser. This conversion follows from fundamental principles of physical chemistry 58: Given the molecular weight of denosumab, Mden = 147 kDa 59, the molar concentration corresponding to a mass concentration of mden,ser = 1 ng/ml is Cden,ser = 6.8027 pM.

Further model calibration has been thoroughly dealt with in 21,24,32; the numerical values of the parameters that are explicitly mentioned in this paper are summarized in Appendix B.

4 NUMERICAL STUDIES
4.1 Simulation of postmenopausal osteoporosis without administration of denosumab
In this section, the model introduced in Section 3.2 is calibrated such that porosity evolutions observed clinically during PMO progression 60 can be replicated, yielding An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu103.jpg pM/d, ξ = 65, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu104.jpg d, and An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu105.jpg d, compare Equations (9) and (10). Considering the such defined PMO-related RANKL-production and the decline of the mechanoresponsiveness of cortical bone, see Figure 4(a), yields corresponding evolutions of osteoclast and osteoblast concentrations, see Figure 4(b). The corresponding increase of the vascular porosity over time agrees well with physiologically observed porosity evolutions in osteoporotic bone 60, see Figure 4(c).

Figure 4
Figure 4
Simulation of postmenopausal osteoporosis (PMO) with disease initiation at tPMO,ini = 0: (a) prescribed temporal evolutions of the reduction factor of the disease-related RANKL production rate, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu106.jpg, and of the mechanoresponsiveness reduction …
Alternatively, to elucidate the bone remodeling kinetics leading to such porosity evolution, the simulation results can be plotted in terms of a comparison of resorbed bone per time (bone resorption rate times concentration of active osteoclasts, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu112.jpg) versus formed bone per time (bone formation rate times concentration of active osteoblasts, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu113.jpg), see Figure 4(d) – such graphical presentation is also referred to as “phase diagram”. A balanced bone turnover implies An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu114.jpg, and related data points lie on the diagonal line in Figure 4(d). A catabolic bone remodeling regime, indicated as light grey-shaded area in Figure 4(c), is characterized by dominating bone resorption, thus An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu115.jpg, and corresponding data points lie below the diagonal line.

On the other hand, an anabolic regime, indicated as light dark-shaded area in Figure 4(d), is characterized by dominating bone formation, thus An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu116.jpg, and corresponding data points lie above the diagonal line.

Onset of PMO, at An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu117.jpg, leads to drastically increasing bone resorption, because of a significant, quasi-instantaneous increase of osteoclast differentiation, and thus implying an increased osteoclast concentration, see Figure 4(b). The maximum ratio of An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu118.jpg over An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu119.jpg (and thus the maximum rate of bone loss) is reached after 5.05 days, indicated by the diamond-shaped markers in Figure 4. Subsequently, osteoclast differentiation decreases because of decreased disease-related RANKL-production, see Figure 4(a), thus the path of the phase diagram moves towards the diagonal line, see Figure 4(d). Furthermore, TGF- β is released during bone resorption, which leads to up-regulation of osteoblast differentiation, indicated by the delayed increase of the osteoblast concentration after onset of PMO in Figure 4(b). Osteoblasts produce RANKL, thus an increased osteoblast concentration implies also increased RANKL production, which, in turn, leads to increased osteoclast differentiation and in further consequence to deceleration of the decrease of the osteoclast concentration because of decreasing PMO-related production of RANKL, see the circle-shaped markers in Figure 4(b). This deceleration provides the explanation for the kink observed in the phase diagram at 76 days, see Figure 4(d). This kink is followed by decreasing bone turnover, meaning that both bone resorption and bone formation slowly diminish. Thus, while still in the catabolic domain, bone remodeling converges to balanced bone resorption and formation, also indicated by the continuously flattening slope of the solid graph in Figure 4(c).

It should be noted that micromechanics-based stiffness homogenization does not only serve as an integral part for implementing mechanoregulation of bone remodeling in our model (as described in Section 3) but also allows for mathematical modeling-based tracking of the bone stiffness over time (not plotted in this paper), in terms of orthotropic Young's moduli E1, E2, and E3, shear moduli G12, G13, and G23, and Poisson's ratios ν12, ν13, and ν23 as stiffness-defining material properties, see 24,61 for details. For example, after 10 years of progressing PMO, the vascular porosity is increased by 17.04%, see Figure 4(c), whereas Young's moduli are decreased by up to 33.55%, shear moduli are decreased by up to 34.64%, and Poisson's ratios are decreased by up to 16.43%. Hence, the effect of PMO on the mechanical integrity of bone can be significantly more severe than the underlying change of the bone composition – in order to simulate this effect, a sound multiscale model of bone mechanics, such as the one utilized in this paper, is essential.

4.2 Single-dose administration of denosumab
Now, we consider a regime of denosumab administration involving only one single administration, to thoroughly elucidate general characteristics of the model response. We prescribe that denosumab is injected 6 months 4 after initiation of PMO (at An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu120.jpg), and we investigate three different administration doses, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu121.jpg mg/kg, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu122.jpg mg/kg, and An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu123.jpg mg/kg. Based on these input data, the model is evaluated in terms of bone cell concentrations, as well as in terms of corresponding phase diagrams and evolutions of the bone constituent volume fractions.

For t < 6 months (prior to administration of denosumab), all graphs (representing doses An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu124.jpg, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu125.jpg, and An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu126.jpg, as well as the zero dose) obviously coincide, see the paths between the square-shaped and the circle-shaped markers in Figures 5(a)–(d). For An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu127.jpgmonths, the presence of denosumab entails inhibition of the differentiation of osteoclast precursor cells to active osteoclasts, consequently leading to fast decrease of An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu128.jpg, see the paths in Figures 5(b)–(d) after passing the circle-shaped markers. At the same time, a temporary increase of An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu129.jpg is observed, leading to a short but steep increase of fbm following denosumab administration, see Figure 5(e). This increase, indicated by the phase diagram paths crossing the balanced turnover-representing diagonal lines in Figures 5(b)–(d), stems from how the action of TGF- β is considered in our model. On the one hand, TGF- β inhibits differentiation of osteoblast precursor cells to active osteoblasts. Administration of denosumab leads to downregulation of the concentration of active osteoclasts, entailing a reduced release of TGF- β, and in further consequence an increase, due to reduced differentiation inhibition, of the concentration of active osteoblasts. Moreover, the presence of TGF- β is required for maintaining differentiation of osteoblast progenitor cells to osteoblast precursor cells. This means that reduction of the concentration of TGF- β leads to downregulation of the aforementioned differentiation process, and thus, with a certain time delay, to downregulation of the concentration of active osteoblasts, completing the “anabolic loops” observed in Figures 5(b)–(d).

Figure 5
Figure 5
Simulation results obtained for administration regimes involving one administration of denosumab: (a)–(d) phase diagrams (in logarithmic scales) for the simulated doses, (a) Dden = 0, (b) An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu134.jpgmg/kg, (c) An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu135.jpgmg/kg, and (d) An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu136.jpgmg/kg, with the …
Subsequent to the temporary dominance of An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu130.jpg over An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu131.jpg, the action of denosumab comes to full effect, and the bone turnover is decreased to a very low level (as compared to the initial state), An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu132.jpg and An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu133.jpg, see the region around the leftwards-pointing triangles in Figures 5(b)–(d). Notably, the higher the dose of denosumab the lower the minimum bone turnover, and the longer the duration of low bone turnover, compare with the time characteristics of the phase diagrams depicted in Figures 5(b)–(d), and the corresponding evolutions of fbm, see Figure 5(e). After clearance of denosumab, the paths of the dose-dependent bone turnover converge to the grey-colored path in Figure 5(a), representing PMO without administration of denosumab. In terms of the simulated evolutions of fbm, this means that after the temporary bone gain “peak”, related to denosumab administration, subsequent evolutions of fbm exhibit the same slope as the evolution of fbm related to PMO only (with Dden = 0), see Figure 5(e).

The temporal evolutions of fbm, depicted in Figure 5(e), reveal a prominent characteristic of the model. While the long-term evolution of the volume fraction of extravascular bone matrix is adequate – denosumab leads, as long as the internal degradation mechanisms described in Section 2 have not led to clearance of denosumab, to a delay of further porosity increase – the short-term response to drug administration does not quite resemble the experimentally and clinically observed trends. The previously described temporary dominance of osteoblast activity over osteoclast activity leads, according to the model, to a short period of very fast bone gain (Δfbm = 0.0125 within Δt ≈ 20 days), see Figure 5(e). After this short period of significant bone gain, max (dfbm / dt) = 8.12 × 10 − 2 % day − 1, the model predicts constant volume fractions, relating to the decrease of bone turnover to a very low level as shown in Figure 5(a), after which the effect of denosumab wears off, and the volume fraction of bone matrix decreases further (consequently the volume fraction of porosity increases), because of uninhibited progress of PMO. While increasing the dose of denosumab leads to lengthening of the volume fraction “plateau”, no dosing effect on the overall bone gain can be observed.

4.3 Multidose administration of denosumab
In order to assess the predictive capabilities of our model, Section 4.3 is devoted to simulation of denosumab administration regimes involving multiple administrations, as commonly investigated in experimental (in vivo) studies and in clinical practice. We present the results of simulations considering three different administered doses, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu140.jpg mg/kg, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu141.jpg mg/kg, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu142.jpg mg/kg, and four different intervals between administrations, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu143.jpg month, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu144.jpg months, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu145.jpg months, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu146.jpg months.

The simulation results are presented in terms of bone matrix volume fraction evolutions, see Figures 6(a) and (c), and corresponding phase diagrams, see Figures 6(b) and (d), to also shed light on the underlying bone turnover kinetics. Obviously, decreasing the administered dose leads to decreasing the anti-catabolic effect of denosumab because of shortened suppression of PMO-driven bone loss, see Figure 6(a). Intuitively, the same should be true for increasing the administration interval – this is for example observed when increasing the administration interval from 9 to 12 months. However, Figure 6(c) also shows that increasing the administration interval from 1 to 6 months results, for the administered dose of Dden = 1.0 mg/kg, in a slight increase in bone gain. On the one hand, this effect is caused by the considered mechanoregulatory mechanisms (see Section 3.1), or in other words by the attempt of bone to re-establish its original composition (for the investigated cortical bone related to fbm = 0.95). On the other hand, administration regimes including multiple administrations of denosumab also exhibit the short period of fast bone gain discussed in Section 4.2, which can, if administered doses and administration intervals are specifically coordinated, as well create an anabolic bone remodeling response – the paths of the phase diagram related to such administration regimes remain, on average, in the anabolic domain for a longer stretch of time.

Figure 6
Figure 6
Simulation results obtained for administration regimes involving multiple administrations of denosumab: temporal evolutions of volume fractions of the extravascular bone matrix, fbm, during progression of PMO (a) for different administered doses (An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu147.jpg mg/kg, …
For intervals An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu154.jpg months and An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu155.jpg months, the long-term influence is negligible due to the opposite effect (a short period of very fast bone loss) observed after clearance of denosumab. For interval An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu156.jpgmonth, the concentration of denosumab constantly remains at a sufficiently high level to maintain suppression of bone remodeling activities so the bone volume fraction remains constant. Increasing the interval of denosumab administration from 1 to 6 months, with the administered dose fixed at Dden = 1.0 mg/kg, results in a higher bone matrix volume fraction: An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu157.jpg, whereas An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu158.jpg. This trend is mainly caused by bone cell kinetics, a certain combination of bone remodeling parameters (characterizing for instance a specific patient), and denosumab administration parameters entail. Hence, our model allows to optimize administration intervals and doses to maximize bone gain after initiation of the administration regime. Clearly, careful tuning of administration interval and dose is the key to obtain such optimized administration regime, which underpins the potential benefit of supporting the design of drug administration regimes by means of computational modeling.

While we have performed a large number of simulations (not shown here) considering different combinations of administration doses and administration intervals, we have never observed any significant bone volume gain. We will further discuss the relevance of this observation with respect to experimental and clinical studies in Section 5.2.

5 DISCUSSION
5.1 Simulation results versus biomarker measurements
In physiological conditions, the direct effect of drug administration on bone remodeling at a particular bone site is hardly measurable. To nevertheless monitor the success of anti-catabolic drugs in the treatment of PMO, so-called bone turnover markers, also referred to as biomarkers, are commonly used in clinical practice 53,62, measured in the urine or in the serum. The evolution of bone resorption is tracked on the basis of bone resorption markers, such as collagen type I cross-linked N-telopeptide (NTx), or collagen type I cross-linked C-telopeptide (CTx). Analogously, bone formation can be monitored by means of bone formation markers, such as osteocalcin (OC or BGLAP), bone-specific alkaline phosphatase (BALP or BSAP), or procollagen type I N-terminal propeptide (PINP). While bone turnover markers provide valuable information on the overall effect of drug administration on bone remodeling, inferring that bone remodeling is altered to the quantitatively same extent on specific sites seems inaccurate (the reason for this is discussed in the next paragraph). We will hence compare the model predictions with biomarker measurements in qualitative fashion only.

Lewiecki 26 summarizes the biochemical effect of denosumab for different single-dose administrations. Compared with the results observed for placebo administration, the “reduction of NTx levels was observed to be dose-dependent, rapid, profound (up to 84%), sustained (for up to six months), and reversible with discontinuation”. The model-predicted evolutions of the concentration of active osteoclasts, serving as basis for computing the results depicted in Figure 5, fit reasonably well to these experimental observations. The only substantial difference is that after administration of denosumab the model-predicted concentration of active osteoclasts is reduced by more than 99% (in contrast to the reported reduction of NTx levels by up to 84% 26). This difference can be explained by considering the involved observation scales. The computational results are based on an RVE in which the denosumab serum concentration is assumed to be uniform, thus denosumab is fully active throughout the RVE. In contrast, blood circulation and thus the distribution of denosumab is not uniform throughout the body, and clearly regions of low supply with denosumab do not contribute to the reduction of NTx as regions where the supply is higher. Taking this into consideration, the agreement between model predictions and experimental results is deemed satisfactory. 5

A qualitatively and quantitatively similar trend is observed for bone formation markers during a multiple denosumab administration regime. For example Eastell et al. 63 observed a reduction of the PINP level of up to ≈ 90% and a reduction of the BALP level of up to ≈ 50%. Considering the importance of the actual observation scale of the results, as proposed in the previous paragraph, the model-predicted osteoblast concentration is assumed to be (qualitatively) consistent with experimental observations.

5.2 Simulation results versus bone mineral density measurements
Clearly, the most tangible effect of PMO is the related decrease of mechanical properties (such as stiffness and strength), leading to an increasing susceptibility to bone fractures. In clinical practice, the fracture risk is usually assessed by measuring the bone mineral density (BMD) by means of dual-energy X-ray absorptiometry on sites suspected or known to be subjected to PMO and comparing this value with the corresponding BMD in healthy bone, for example, via the so-called T-score 64. Often, the success of anti-resorptive drugs, such as denosumab, in the treatment of osteoporosis is assessed by measuring the BMD evolution after drug administration. These measurements are then interpreted as follows: the higher the increase of the BMD the more potent the investigated drug. For example, after administration of denosumab the BMD has been reported to increase by up to 2 % in cortical bone and up to 9 % in trabecular bone (over a time span of typically 1 to 2 years), see for example 11,9,65–68. Mineral concentrations are not included into the present version of our mathematical model, so that a direct comparison of BMD increase with simulation results is out of reach for the moment. However, it is very instructive to compare the increase in BMD measured on denosumab-treated patients with the decreasing or constant bone volume predicted by our model in Figures 5 and ​and66.

By definition, denosumab (and in general anti-resorptive drugs) are intended to decrease the activity of osteoclasts. Our simulation results (see Section 4), as well as measurements of bone resorption markers 26 confirm that denosumab indeed fulfills this intention. Simultaneously, because of communication between osteoclasts and osteoblasts, denosumab also entails significant decrease of the activity of osteoblasts – confirmation of this computationally observed behavior follows from measurements of bone formation markers 63. Evidence is thus overwhelming that the increase of bone mineral density after administration of denosumab simply cannot arise from bone gain in the sense that osteoblasts lay down “new” osteoid. The only mechanism remaining that is able to cause the clinically observed BMD increase after administration of denosumab is increased bone tissue mineralization, triggered through denosumab. The following mechanism seems very probable. Denosumab suppresses the action of RANKL, which is itself promoting the dissolution actions by osteoclasts and osteocytes. In particular, the latter are responsible for keeping the tissue mineralization degree at some limit value 69,70, that is, inhibiting purely chemically driven “anorganification” of bone tissue. Suppressed osteocytic dissolution action, however, would allow for progression of the so-called secondary mineralization process, and also increase the overall BMD, as evidenced in clinical studies 71–78. Hence, our model gives valuable insight into the biological and chemical processes leading to denosumab-triggered BMD increase. For the sake of completeness, we remark that secondary mineralization-induced increase of bone tissue stiffness with respect to Equation (A.9) is not explicitly accounted for in the present approach, adding that, from a mechanical point of view, a rigorous multiscale approach translating mineralization degree to tissue stiffness is at hand 79,70.

5.3 Influence of the mechanical loading on the disease progress
In the computational studies presented in this paper, we have considered a constant (that is time-invariant) uniaxial mechanical loading. It is perfectly clear that (i) such mechanical load case is a strong simplification of the “real” loading acting onto human bone (which is presumably always of three-dimensional nature, including nonzero shear components of the stress and strain tensors); and that (ii) mechanical loading varies (partly significantly) between different locations distributed across a human skeleton. In the context of the scope of this paper, the influence of a different load case is negligible as long as it is still time-invariant, and the resulting strain energy density is of comparable magnitude. This is because in our model, mechanical regulation is stimulated by changes of the strain energy density, irrespective of the exact shape of the related strain tensor. A temporarily changing mechanical loading, on the other hand, could potentially have a great influence on the disease progress (in terms of the resulting bone matrix volume fraction evolution); imposing such mechanical loading “histories” is, however, beyond the scope of this paper, which is restricted to mathematical modeling of PMO and its pharmaceutical intervention by means of the anti-resorptive drug denosumab. For a parametric study on the influence of the type of mechanical loading on the related bone remodeling behavior, see reference 24.

Another important model aspect, related to the mechanical loading, is the sensitivity of the bone remodeling response to a changing (local) mechanical loading. In the simulations presented in this paper, the macroscopic mechanical loading is held constant at all times (as discussed in the previous paragraph). However, the disease-related increase of the vascular porosity (after the onset of PMO, see Figures 5 and ​and6)6) implies an increase of the strains experienced by the extravascular bone matrix, and thus an increase of the related strain energy density, which is the quantity that has been chosen as mechanoregulatory stimulus. An increase of the strain energy density, beyond the value that occurs if the studied piece of cortical bone is healthy, leads to an increase of OBp proliferation, with the anabolic strength parameter λ as governing property, compare Equations (3) and (7). In the following, because of the apparent importance of parameter λ, we elucidate the sensitivity of the model predictions to variations of λ. In particular, we perform the simulations based on which Figure 6(a) has been created (that is a denosumab administration regime involving multiple administrations, with an administration interval of 9 months) for two alternative values of λ: λlow = 0.125 (representing bone with low mechanoresponsiveness) and λhigh = 12.5 (representing bone with high mechanoresponsiveness). In Figure 7, the volume fraction of the extravascular bone matrix after 10 years of PMO progress is plotted, confirming the striking sensitivity of the model predictions to parameter λ.

Figure 7
Figure 7
Study of the sensitivity of the model-predicted bone matrix volume fraction fbm after 10 years of postmenopausal osteoporosis progress with respect to a varying anabolic strength parameter λ, computed for an administration regime involving …
This parametric study also shows that the sensitivity of the model prediction to λ is indirectly proportional to the administered dose of denosumab: the higher the administered dose, the lower the sensitivity of the bone remodeling response to variations of λ. Thus, careful choice of the model parameters accounting for the mechanoresponsiveness of bone remodeling is crucial for the mathematical model presented here in order to replicate adequate predictions.

6 CONCLUSIONS
In this paper, an integrated model has been presented, allowing for simulation of the effects of postmenopausal osteoporosis on bone remodeling. For this purpose, a PK model of denosumab has been developed and coupled to a bone cell population model of bone remodeling, which itself communicates with a continuum micromechanics-based homogenization scheme of bone stiffness, in order to soundly take into account mechanobiological regulatory mechanisms of bone remodeling and corresponding changes of the bone stiffness.

The numerical studies presented in Sections 4 and 5 have shown that the new model proposed here allows for (qualitatively and quantitatively) reasonable simulation of osteoporotic scenarios, and of the pharmaceutical intervention through administering the anti-catabolic drug denosumab. In particular, the simulation results have given rise to the following specific conclusions:

In order to simulate PMO, the underlying mathematical model should be fed by temporally varying disease-triggering factors – in our model, PMO is driven, on the one hand, by disease-related RANKL production, see Equations (9) and (10), and, on the other hand, by disease-related reduction of the mechanoresponsiveness, see Equation (11). This way, not only the postmenopausal evolution of bone volume agrees well with corresponding clinical data, see Figure 4(c), but also the related bone turnover satisfyingly resembles the clinically and experimentally observed dynamics – a very high bone turnover in the initial phase of PMO is followed by a slow, long-lasting decrease converging to the normal, premenopausal bone turnover, see Figure 4(d).
Administration of denosumab entails, because of inhibition of osteoclast differentiation, down-regulation of both bone resorption and, with some time delay, bone formation to an ineffective level. One administration of denosumab has the following effects on the bone turnover: Initially, a short but significant phase of bone gain is observed, related to the time lag between drug-related inhibition of bone resorption and bone formation. Once bone formation is sufficiently down-regulated (approximately 5 to 10 days after denosumab administration, depending on the administered dose), the bone turnover reduces to a negligible level. Due to denosumab clearance, the bone turnover finally approaches the disease-related evolution, thus altogether denosumab administration causes a delayed progress of PMO. As expected, simulation results confirm that increasing the administered dose and/or decreasing the administration interval improves the anti-catabolic effect of denosumab, compare Figures 5(b)–(e).
The model predictions obtained through simulating denosumab administration regimes involving multiple administrations, see Section 4.3, agree, in qualitative terms, well with corresponding biomarker measurements.
One central component of the integrated model is the continuous estimation of the macroscopic bone stiffness by means of a continuum micromechanics-based homogenization scheme, appropriately taking into account the anisotropic mechanical behavior of bone. The such computed evolution of bone stiffness could serve as basis for assessing the risk of bone failure.
These findings highlight that the model proposed here could be a promising tool for computer simulation-based assessment of the effects of denosumab administration on osteoporotic bone. However, the results presented and discussed in Sections 4 and 5 also motivate to invest in explicit modeling of the secondary mineralization process and its linkage to RANKL and osteocytic action.

In addition to their important role for the mineralization process (as discussed previously), osteocytes are known to act as key players for sensing changes of the mechanical loading and for transducing these changes into corresponding biochemical events, leading to increase or decrease of the bone matrix volume fraction 80,81,23,82. In the systems biology approach employed in this paper, osteocytes are not considered explicitly, which clearly constitutes a limitation of the model. Introduction of the population of osteocytes into our model, as demonstrated in alternative modeling approaches, see for example references 83–86, is thus deemed as promising model extension, in particular to bring the considered mechanoregulatory mechanisms closer to experimental evidence 81,82,87,88. Furthermore, in future model extensions, we will aim at reduction of the number of model components to the necessary minimum; this will go along with reduction of the currently considerable number of model parameters, being another limitation of the model. These issues are topics of ongoing research.

Acknowledgments
Financial support by the Australian Research Council, in the framework of the project Multi-scale modeling of transport through deformable porous materials (project number DP-0988427), by the Australian Academy of Science, in the framework of Scientific Visits to Europe with the COST Action program and by the European Research Council, in the framework of the project Multiscale poro-micromechanics of bone materials, with links to biology and medicine (project number FP7-257023), is gratefully acknowledged.

APPENDIX A: ANALYTICAL SOLUTIONS OF THE PHARMACOKINETICS MODEL
A.1 Single-dose administration regimes
The analytical solution of differential equation (1), under consideration of initial condition mden,sub(t = 0) = Δmden,sub, see Section 2.1, reads

equation image A.1
Note that Δmden,sub is related to the administered dose Dden through consideration of an average patient body mass and serum volume, as discussed in Section 2.2. On the other hand, the solution of differential equations (2) is based on conditions mden,ser(t = 0) = 0, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu163.jpg, and An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu164.jpg. Solving Equations (2) for mden,ser yields

equation image A.2
equation image A.3
and

equation image A.4
Numerical evaluation of Equations (A.3) and (A.4) requires knowledge of time instants An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu168.jpg and An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu169.jpg. If the numerical value of An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu170.jpg is known, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu171.jpg and An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu172.jpg are accessible implicitly via Equations (A.2) and (A.3): Equation (A.2), evaluated for An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu173.jpg, allows for determination of An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu174.jpg, whereas Equation (A.3), evaluated for An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu175.jpg, allows for determination of An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu176.jpg.

A.2 Multidose administration regimes We now want to derive analytical expressions describing the denosumab serum concentration over time for denosumab administration regimes involving multiple injections.
For this purpose, administration intervals and doses are assumed to be constant, and the corresponding initial conditions are periodic. 6 Because of this regularity, only two different scenarios have to be considered: at the end of the previous administration interval the denosumab serum concentration mden,ser is (i) lower than the critical concentration An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu177.jpg, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu178.jpg, and (ii) it is higher (or equal) An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu179.jpg, see Figure 8.

Figure 8
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-f8.jpg
Two scenarios investigated for denosumab administration regimes involving multiple injections: in scenario (i), mden,ser starts from a level below An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu180.jpg in each administration interval ς, while in scenario (ii), mden,ser remains higher than An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu181.jpg at all times; within one administration interval scenario (i) comprises three domains ①, ②, and ③ based on whether mden,ser is lower or higher than An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu182.jpg (the transition between the domains is indicated by the white-faced circle-shaped marker, representing An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu183.jpg and black-faced circle-shaped marker, representing An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu184.jpg), following Equations (A.5)–(A.7), whereas scenario (ii) follows Equation (A.8).

First, we investigate scenario (i). Denosumab is administered at the beginning of each administration interval ς, at t = tad,ς, with tad,1 = 0. The administered dose and thus Δmden,sub are constant, and at the end of each interval the denosumab serum concentration has dropped to a level below the critical concentration, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu185.jpg. Consequently, each administration interval is characterized by domains of different enzyme kinetics, analogous to the single-administration regime elaborated in Sections 2.2 and A, which are delimited by the time instants at which the critical concentration is reached, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu186.jpg and An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu187.jpg, see also Figure A.1. The following notation is chosen for the different enzyme degradation regimes: An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu01.jpg denotes the denosumab serum concentration if An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu188.jpg, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu02.jpg denotes the denosumab serum concentration if An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu189.jpg, and An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu03.jpg denotes the denosumab serum concentration if An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu190.jpg. While An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu04.jpg, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu05.jpg, and An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu06.jpg are defined by Eqs.(A.2)–(A.4), An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu07.jpg, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu08.jpg, and An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu09.jpg follow as

equation image A.5
equation image A.6
and

equation image A.7
In scenario (ii) the denosumab serum concentration remains, after exceeding An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu194.jpg at An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu195.jpg in the first administration interval, above An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu196.jpg at all times. Thus, distinction between different domains is, unlike scenario (i), not necessary, and only one analytic expression has to be derived. If An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu197.jpg, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu198.jpg reads as

equation image A.8
where function An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu200.jpg is defined as An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu201.jpg if An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu202.jpg and An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu203.jpg otherwise.

APPENDIX B: MODEL PARAMETERS
In Table B.1, the model parameters of the integrated systems biology/micromechanical model, introduced in Section 3, that are explicitly mentioned in this paper, are summarized.

Table B.1
Parameters of the integrated systems biology-micromechanical model, see also 21,32.

parameter numerical value unit
kres 2 (pMday)  − 1
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu204.jpg 7 × 10 − 2 d  − 1
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu205.jpg 1.6570 × 10 − 1 d  − 1
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu206.jpg 2.1 × 100 d  − 1
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu207.jpg 2.1107 × 10 − 1 d  − 1
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu208.jpg 5.6487 × 10 − 4 d  − 1
Kd,[RANKL-d] 5.6797 × 100 pM
βRANKL 1.6842 × 102 pMd  − 1
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu209.jpg 1.0132 × 101 d  − 1
Ka,[RANKL-RANK] 3.4118 × 10 − 2 pM  − 1
Ka,[RANKL-OPG] 1 × 10 − 3 pM  − 1
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu210.jpg 0.5 –
λ 1.25 –
κ 105 pM/day
The bone formation rate kform is calibrated such that for the bone resorption rate, kres, defined in Table 1, steady-state cell concentrations, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu211.jpg, imply a balanced bone turnover; that is, fvas = const. and fbm = const.: An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu212.jpg. Calibration of the proliferation rate An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu213.jpg is described in detail in 24.

Furthermore, the stiffness tensor of the extravascular bone matrix, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu214.jpg, is defined in the line of Fritsch and Hellmich 79; based on the ultrasonics tests by Ashman et al. 89, conducted on human femurs, An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu215.jpg reads in compressed notation 90

equation image A.9
APPENDIX C: NOMENCLATURE
Abbreviations:
BALP bone-specific alkaline phosphatase
BCPM bone cell population model
BMD bone mineral density
CTx collagen type I cross-linked C-telopeptide
DXA dual-energy X-ray absorptiometry
NTx collagen type I cross-linked N-telopeptide
OBa active osteoblasts
OBp osteoblast precursor cells
OBu uncommitted osteoblast progenitor cells
OC osteocalcin
OCa active osteoclasts
OCp osteoclast precursor cells
OPG osteoprotegerin
PD pharmacodynamics
PINP procollagen type I N-terminal propeptide
PK pharmacokinetics
PMO postmenopausal osteoporosis
PTH parathyroid hormone
RANK receptor of nuclear factor kappa-B
RANKL ligand of RANK
RVE representative volume element
SED strain energy density
TGF- β transforming growth factor beta
Latin symbols:
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu217.jpg constant apoptosis rate of active osteoblasts
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu218.jpg maximum apoptosis rate of active osteoclasts
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu219.jpg stiffness tensor of the extravascular bone matrix
Cden molar concentration of denosumab
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu220.jpg molar concentration of active osteoblasts
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu221.jpg molar concentration of active osteoclasts
COPG molar concentration of OPG
CRANK molar concentration of RANK
CRANKL molar concentration of RANKL
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu222.jpg maximum molar concentration of RANKL
Dden dose of denosumab administration
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu223.jpg constant RANKL degradation rate
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu224.jpg maximum differentiation rate of osteoblast precursor cells
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu225.jpg maximum differentiation rate of uncommitted osteoblast progenitor cells
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu226.jpg maximum differentiation rate of osteoclast precursor cells
Ei Young's moduli of bone on the macroscopic observation scale, i = 1,2,3
fbm volume fraction of extravascular bone matrix
Gij shear moduli of bone on the macroscopic observation scale, ij = 12,23,13
kabs absorption rate
kenz degradation rate related to enzymes
kform bone formation rate
korg degradation rate related to organs
kres bone resorption rate
Ka,[RANKL-i] equilibrium association binding constant for the binding of RANKL to substance i, i = denosumab, OPG, RANK
Kd,[RANKL-i] equilibrium dissociation binding constant for the dissociation of substance i from RANKL, i = denosumab, RANK
mden,ser denosumab mass concentration in the blood serum
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu227.jpg limit denosumab mass concentration in the blood serum
Mden molecular weight of denosumab
Δmden,sub increase of denosumab mass concentration in the blood serum due to administration
mden,sub denosumab mass concentration in the subcutaneous tissue
PRANKL production rate of RANKL
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu228.jpg PMO-related production rate of RANKL
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu229.jpg PMO-related production rate of RANKL directly after disease initiation
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu230.jpg mechanical disuse-related production rate of RANKL
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu231.jpg maximum proliferation rate of osteoblast precursor cells
t time variable
tPMO,ini point in time when PMO is initiated
Vbm volume of the extravascular bone matrix within a representative volume element
VRVE volume of the representative volume element
Vvas volume of the vascular pore space within a representative volume element
Greek symbols:
βRANKL intrinsic RANKL production rate
Ψbm microscopic strain tensor experienced by the extravascular bone matrix
ζ denosumab accessibility factor
κ RANKL production inhibition parameter
λ anabolic strength parameter
νij Poisson's ratios of bone on the macroscopic observation scale, ij = 12,13,23
ξ parameter introduced for definition of An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu232.jpg
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu233.jpg activator function for the differentiation of uncommitted osteoblast progenitor cells due to the action of TGF- β
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu234.jpg activator function for the differentiation of osteoblast precursor cells due to the action of TGF- β
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu235.jpg activator function for the differentiation of osteoclast precursor cells due to the action of RANKL
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu236.jpg activator function for the proliferation of osteoblast precursor cells due to mechanical loading
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu237.jpg activator function for the proliferation of osteoblast precursor cells due to normal mechanical loading
σbm microscopic stress tensor experienced by the extravascular bone matrix
ς administration interval index in a multidose administration regime
Σ macroscopic stress tensor
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu238.jpg characteristic time of the RANKL production decrease
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu239.jpg mechanoresponsiveness reduction factor
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu240.jpg reduction factor An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu241.jpg
Ψbm strain energy density experienced on the observation scale of extravascular bone matrix
An external file that holds a picture, illustration, etc.
Object name is cnm0030-0001-mu242.jpg strain energy density experienced on the observation scale of extravascular bone matrix related to normal mechanical loading
Footnotes
Note that, in terms of evaluation of the model presented in Sections 3 and 4, we restrict ourselves to cortical bone. Extension of the model to trabecular bone is straightforward, merely requiring adjustment of some model parameters.

However, the micromechanical model employed here is able to straightforwardly consider any desired load case; see reference 24 and Section 5.3 for a discussion on the influence of the mechanical loading on the model predictions.

Estrogen inhibits osteoclastogenesis; that is, maturation of osteoclastic cells 55,54. Menopause entails a significant drop of the estrogen level, consequently leading to omission of this inhibition, followed by increased osteoclast activity and thus loss of bone.

The delay between onset of PMO and the administration of denosumab takes into account that the onset of PMO is not noticeable in real life, thus its treatment would be initiated after PMO has already been active for some time. Increasing or decreasing this delay implies quantitatively different simulation results; qualitative features of respective simulation results conform, however, with the ones described here.

In our model, the “uneven” distribution of denosumab across the human body could be taken into account by prescribing a corresponding distribution of the accessibility factor ζ, as introduced in Equation (12). At present, such accessibility distribution is, however, not available from experimental data.

For administration intervals and/or doses varying over time, the governing differential equations have to be solved in piecewise fashion (administration interval after administration interval). In such case aid of numerical solution methods is advisable. In clinical practice, such approach is, however, not standard, due to limited patient-friendliness.

Article information
Int J Numer Method Biomed Eng. 2014 Jan; 30(1): 1–27.
Published online 2013 Aug 30. doi: 10.1002/cnm.2584
PMCID: PMC4291103
S Scheiner,1,*† P Pivonka,2 D W Smith,3 C R Dunstan,4 and C Hellmich1
1Institute for Mechanics of Materials and Structures, Vienna University of Technology, Austria
2Australian Institute for Musculoskeletal Science, The University of Melbourne, Australia
3Faculty of Engineering, Computing and Mathematics, The University of Western Australia, Australia
4School of Aerospace, Mechanical and Mechatronic Engineering, The University of Sydney, Australia
*Correspondence to: S. Scheiner, Institute for Mechanics of Materials and Structures, Vienna University of Technology, Karlsplatz 13/202, A-1040 Vienna, Austria.
†E-mail: ta.ca.neiwut@reniehcs.nafets
Received 2013 Feb 25; Revised 2013 Jul 3; Accepted 2013 Jul 9.
Copyright © 2013 John Wiley & Sons, Ltd.
REFERENCES
Gerend M, Erchull M, Aiken L, Maner J. Reasons and risk: factors underlying women's perceptions of susceptibility to osteoporosis. Maturitas. 2006;55(3):227–237. [PubMed]
Delmas P. Treatment of postmenopausal osteoporosis. Lancet. 2002;359(9322):2018–2026. [PubMed]
Rubin M, Bilezikian J. New anabolic therapies in osteoporosis. Current Opinion in Rheumatology. 2002;14(4):433–440. [PubMed]
Black D, Greenspan S, Ensrud K, Palermo L, McGowan J, Lang T, Garnero P, Bouxsein M, Bilezikian J, Rosen C. The effects of parathyroid hormone and alendronate alone or in combination in postmenopausal osteoporosis. New England Journal of Medicine. 2003;349(13):1207–1215. [PubMed]
Amgen Inc. 2011. . Available from: http://www.amgen.com [Accessed 29 July 2013]
Boyce B, Xing L. Functions of RANKL/RANK/OPG in bone modeling and remodeling. Archives of Biochemistry and Biophysics. 2008;473(2):139–146. [PMC free article] [PubMed]
Martin R, Burr D, Sharkey N. Skeletal Tissue Mechanics. New York, USA: Springer Verlag; 1998.
Bekker P, Holloway D, Rasmussen A, Murphy R, Martin S, Leese P, Holmes G, Dunstan C, DePaoli A. A single-dose placebo-controlled study of AMG 162, a fully human monclonal antibody to RANKL, in postmenopausal women. Journal of Bone and Mineral Research. 2004;19(7):1059–1066. [PubMed]
McClung M, Lewiecki E, Cohen S, Bolognese M, Woodson G, Moffett A, Peacock M, Miller P, Lederman S, Chesnut C, et al. Denosumab in postmenopausal women with low bone mineral density. New England Journal of Medicine. 2006;354(8):821–831. [PubMed]
Body JJ, Facon T, Coleman R, Lipton A, Geurs F, Fan M, Holloway D, Peterson M, Bekker P. A study of the biological receptor activator of nuclear factor- κB ligand inhibitor, denosumab, in patients with multiple myeloma or bone metastases from breast cancer. Clinical Cancer Research. 2006;12(4):1221–1228. [PubMed]
Lewiecki E, Miller P, McClung M, Cohen S, Bolognese M, Liu Y, Wang A, Siddhanti S, Fitzpatrick L. Two-year treatment with denosumab (AMG162) in a randomized phase 2 study of postmenopausal women with low BMD. Journal of Bone and Mineral Research. 2007;22(12):1832–1841. [PubMed]
Ellis G, Bone H, Chlebowski R, Paul D, Spadafora S, Smith J, Fan M, Jun S. Randomized trial of denosumab in patients receiving adjuvant aromatase inhibitors for nonmetastatic breast cancer. Journal of Clinical Oncology. 2008;26(30):4875–4882. [PubMed]
U.S. Food and Drug Administration. 2011. . Available from: http://www.fda.gov [Accessed 29 July 2013]
Holford N, Kimko H, Monteleone J, Peck C. Simulation of clinical trials. Annual Review of Pharmacology and Toxicology. 2000;40:209–234. [PubMed]
Mould D, Sweeney K. The pharmacokinetics and pharmacodynamics of monoclonal antibodies – mechanistic modeling applied to drug development. Current Opinion in Drug Discovery and Development. 2007;10(1):84–96. [PubMed]
Post T, Cremers S, Kerbusch T, Danhof M. Bone physiology, disease and treatment – towards disease system analysis in osteoporosis. Clinical Pharmacokinetics. 2010;49(2):89–118. [PubMed]
Sheiner L. Learning versus confirming in clinical drug development. Clinical Pharmacology and Therapeutics. 1997;61(3):275–291. [PubMed]
Sharma A, Jusko W. Characterization of four basic models of indirect pharmacodynamic responses. Journal of Pharmacokinetics and Biopharmaceutics. 1996;24(6):611–635. [PubMed]
Post T, Freijer J, DeJongh J, Danhof M. Disease system analysis: basic disease progression models in degenerative disease. Pharmaceutical Research. 2005;22(7):1038–1049. [PubMed]
Marathe A, Peterson M, Mager D. Integrated cellular bone homeostasis model for denosumab pharmacodynamics in multiple myeloma patients. Journal of Pharmacology and Experimental Therapeutics. 2008;326(2):555–562. [PMC free article] [PubMed]
Pivonka P, Zimak J, Smith D, Gardiner B, Dunstan C, Sims N, Martin T, Mundy G. Theoretical investigation of the role of the RANK-RANKL-OPG system in bone remodeling. Journal of Theoretical Biology. 2010;262(2):306–316. [PubMed]
Peterson M, Riggs M. Predicting nonlinear changes in bone mineral density over time using a multscale systems pharmacology model. Pharmacometrics & Systems Pharmacology. 2012;1:1–8. [PMC free article] [PubMed]
Robling A, Castillo A, Turner C. Biomechanical and molecular regulation of bone remodeling. Annual Review of Biomedical Engineering. 2006;8:455–498. [PubMed]
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:181–196.
Shargel L, Wu-Pong S, Yu A. Applied Biopharmaceutics & Pharmacokinetics. 5th ed. New York, USA: McGraw-Hill; 2005.
Lewiecki E. Treatment of osteoporosis with denosumab. Maturitas. 2010;66(2):182–186. [PubMed]
Murray J. Mathematical Biology – I: An Introduction. 3rd ed. Berlin, Germany: Springer Science and Business Media, Inc; 2002.
Schwefel HP. Numerische Optimierung von Computer-Modellen mittels der Evolutionsstrategie [Numerical Optimization of Computer Models by means of the Evolution Strategy] Basel, Switzerland and Stuttgart, Germany: Birkhäuser Verlag; 1977. . In German.
McDowell M, Fryar C, Ogden C, Flegal K. National Health Statistics Reports. 2008. Anthropometric reference data for children and adults: United states, 2003–2006 , No. 10, Centers for Disease Control and Prevention,
Cirillo M, Laurenzi M, Tervisan M, Stamler J. Hematocrit, blood pressure, and hypertension: the Gubbio population study. Hypertension. 1992;20(3):319–326. [PubMed]
Lemaire V, Tobin F, Greller L, Cho C, Suva L. Modeling of the interactions between osteoblast and osteoclast activities in bone remodeling. Journal of Theoretical Biology. 2004;229(3):293–309. [PubMed]
Pivonka P, Zimak J, Smith D, Gardiner B, Dunstan C, Sims N, Martin T, Mundy G. Model structure and control of bone remodeling: a theoretical study. Bone. 2008;43(2):249–263. [PubMed]
Pivonka P, Buenzli P, Scheiner S, Hellmich C, Dunstan C. The influence of bone surface availability in bone remodelling – a mathematical model including coupled geometrical and biomechanical regulations of bone cells. Engineering Structures. 2013;47:134–147.
Janssens K, ten Dijke P, Janssens S, Van Hul W. Transforming growth factor- β1 to the bone. Endocrine Reviews. 2005;26(6):743–774. [PubMed]
Erlebacher A, Filvaroff E, Ye JQ, Derynck R. Osteoblastic responses to TGF- β during bone remodeling. Molecular Biology of the Cell. 1998;9(7):1903–1918. [PMC free article] [PubMed]
Jones D, Nolte H, Scholübbers JG, Turner E, Veltel D. Biochemical signal transduction of mechanical strain in osteoblast-like cells. Biomaterials. 1991;12(2):101–110. [PubMed]
Owan I, Burr D, Turner C, Qiu J, Tu Y, Onyia J, Duncan R. Mechanotransduction in bone: osteoblasts are more responsive to fluid forces than mechanical strain. American Journal of Physiology – Cell Physiology. 1997;273(3):C810–C815. [PubMed]
Kaspar D, Seidl W, Neidlinger-Wilke C, Beck A, Claes L, Ignatius A. Proliferation of human-derived osteoblast-like cells depends on the cycle number and frequency of uniaxial strain. Journal of Biomechanics. 2002;35(7):873–880. [PubMed]
Weyts F, Bosmans B, Niesing R, Van JL, Weinans H. Mechanical control of human osteoblast apoptosis and proliferation in relation to differentiation. Calcified Tissue International. 2003;72(4):505–512. [PubMed]
Fuller K, Lean J, Bayley K, Wani M, Chambers T. A role for TGF- β1 in osteoclast differentiation and survival. Journal of Cell Science. 2000;113(13):2445–2453. [PubMed]
Pitsillides A, Rawlinson S, Suswillo R, Bourrin S, Zaman G, Lanyon L. Mechanical strain-induced NO production by bone cells: a possible role in adaptive bone (re)modeling? FASEB Journal. 1995;9(15):1614–1622. [PubMed]
Mullender M, El Haj A, Yang Y, van Duin M, Burger E, Klein-Nulend J. Mechanotransduction of bone cells in vitro: mechanobiology of bone tissue. Medical and Biological Engineering and Computing. 2004;42(1):14–21. [PubMed]
Fan X, Roy E, Zhu L, Murphy T, Ackert-Bicknell C, Hart C, Rosen C, Nanes M, Rubin J. Nitric oxide regulates receptor activator of nuclear factor- κB ligand and osteoprotegerin expression in bone marrow stromal cells. Endocrinology. 2004;145(2):751–759. [PubMed]
Liu C, Zhao Y, Cheung WY, Gandhi R, Wang L, You L. Effects of cyclic hydraulic pressure on osteocytes. Bone. 2010;46(5):1449–1456. [PMC free article] [PubMed]
Alon U. An Introduction to Systems Biology: Design Principles of Biological Circuits. London, UK: Chapman & Hall; 2007.
Wergedal J, Baylink D. Electron microprobe measurements of bone mineralization rate in vivo. American Journal of Physiology. 1974;226(2):345–352. [PubMed]
Grynpas M. Age and disease-related changes in the mineral of bone. Calcified Tissue International. 1993;53(S1):57–64. [PubMed]
Busa B, Miller L, Rubin C, Qin YX, Judex S. Rapid establishment of chemical and mechanical properties during lamellar bone formation. Calcified Tissue International. 2005;77(6):386–394. [PubMed]
Hellmich C, Kober C, Erdmann B. Micromechanics-based conversion of CT data into anisotropic elasticity tensors, applied to FE simulations of a mandible. Annals of Biomedical Engineering. 2008;36(1):108–122. [PubMed]
Zaoui A. Continuum micromechanics: survey. Journal of Engineering Mechanics (ASCE) 2002;128(8):808–816.
Riggs B, Khosla S, MeltonIII L. A unitary model for involutional osteoporosis: estrogen deficiency causes both type I and type II osteoporosis in postmenopausal women and contributes to bone loss in aging men. Journal of Bone and Mineral Research. 1998;13(5):763–773. [PubMed]
Manolagas S. Birth and death of bone cells: basic regulatory mechanisms and implications for the pathogenesis and treatment of osteoporosis. Endocrine Reviews. 2000;21(2):115–137. [PubMed]
Garnero P, Sornay-Rendu E, Chapuy MC, Delmas P. Increased bone turnover in late postmenopausal women is a major determinant of osteoporosis. Journal of Bone and Mineral Research. 1996;11(3):337–349. [PubMed]
Hofbauer L, Schoppet M. Clinical implications of the osteoprotegerin/RANKL/RANK system for bone and vascular diseases. Journal of the American Medical Association. 2004;292(4):490–495. [PubMed]
Riggs B. The mechanisms of estrogen regulation of bone resorption. The Journal of Clinical Investigation. 2000;106(10):1203–1204. [PMC free article] [PubMed]
Tomkinson A, Gevers E, Wit J, Reeve J, Noble B. The role of estrogen in the control of rat osteocyte apoptosis. Journal of Bone and Mineral Research. 1998;13(8):1243–1250. [PubMed]
Kostenuik P, Nguyen H, McCabe J, Warmington K, Kurahara C, Sun N, Chen C, Li L, Cattley R, Van G, 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;24(2):182–195. [PubMed]
Atkins P. Physical Chemistry. 6th ed. Oxford, UK: Oxford University Press; 1998.
Green W. Denosumab (Prolia) injection – a new approach to the treatment of women with postmenopausal osteoporosis. Pharmacy and Therapeutics. 2010;35(10):553–559.
Bonnet N, Ferrari S. Exercise and the skeleton: how it works and what it really does. IBMS BoneKey. 2010;7(7):235–248.
Lekhnitskii S. Theory of Elasticity of an Anisotropic Elastic Body. San Francisco, USA: Holden-Day, Inc; 1963.
Looker A, Bauer D, Chesnut III C, Gundberg C, Hochberg M, Klee G, Kleerekoper M, Watts N, Bell N. Clinical use of biochemical markers of bone remodeling: current status and future directions. Osteoporosis International. 2000;11(6):467–480. [PubMed]
Eastell R, Christiansen C, Grauer A, Kutilek S, Libanati C, McClung M, Reid I, Resch H, Siris E, Uebelhart D, Wang A, Weryha G, Cummings S. Effects of denosumab on bone turnover markers in postmenopausal osteoporosis. Journal of Bone and Mineral Research. 2011;26(3):530–537. [PubMed]
Lu Y, Genant HK, Sheperd J, Zhao S, Mathur A, Fuerst T, Cummings SR. Classification of osteoporosis based on bone mineral densities. Journal of Bone and Mineral Research. 2001;16(5):901–910. [PubMed]
Bone H, Bolognese M, Yuen C, Kendler D, Wang H, Liu Y, Martin J. Effects of denosumab on bone mineral density and bone turnover in postmenopausal women. Journal of Clinical Endocrinology and Metabolism. 2008;93(6):2149–2157. [PubMed]
Kendler D, Roux C, Benhamou C, Brown J, Lillelstol M, Siddhanti S, Man HS, Martin J, Bone H. Effects of denosumab on bone mineral density and bone turnover in postmenopausal women transitioning from alendronate therapy. Journal of Bone and Mineral Research. 2010;25(1):72–81. [PubMed]
Miller P. Denosumab: anti-RANKL antibody. Current Osteoporosis Reports. 2009;7(1):18–22. [PubMed]
Miller P, Wagman R, Peacock M, Lewiecki E, Bolognese M, Weinstein R, Ding B, Martin J, McClung M. Effect of denosumab on bone mineral density and biochemical markers of bone turnover: six-year results of a phase 2 clinical trial. Journal of Clinical Endocrinology and Metabolism. 2011;96(2):394–402. [PubMed]
Bala Y, Farlay D, Delmas P, Meunier P, Boivin G. Time sequence of secondary mineralization and microhardness in cortical and cancellous bone from ewes. Bone. 2010;46(4):1204–1212. [PubMed]
Vuong J, Hellmich C. Bone fibrillogenesis and mineralization: quantitative analysis and implications for tissue elasticity. Journal of Theoretical Biology. 2011;287:115–130. [PubMed]
Chavassieux P, Arlot M, Reda C, Wei L, Yates A, Meunier P. Histomorphometric assessment of the long-term effects of alendronate on bone quality and remodeling in patients with osteoporosis. The Journal of Clinical Investigation. 1997;100(6):1475–1480. [PMC free article] [PubMed]
Stepan J, Alenfeld F, Boivin G, Feyen J, Lakatos P. Mechanisms of action of antiresorptive therapies of postmenopausal osteoporosis. Endocrine Regulations. 2003;37(4):225–238. [PubMed]
Seeman E. Reduced bone formation and increased bone resorption: rational targets for the treatment of osteoporosis. Osteoporosis International. 2003;14(S2-S8) [PubMed]
Boivin G, Meunier P. The mineralization of bone tissue: a forgotten dimension in osteoporosis research. Osteoporosis International. 2003;14:S19–S24. [PubMed]
Bone H, Hosking D, Devogelaer JP, Tucci J, Emkey R, Tonino R, Rodriguez-Portales J, Downs R, Gupta J, Santora A, Liberman U. Ten years' experience with alendronate for osteoporosis in postmenopausal women. The New England Journal of Medicine. 2004;350(12):1189–1199. [PubMed]
Lewiecki E, Silverman S. Redefining osteoporosis treatment: who to treat and how long to treat. Arquivos Brasileiros de Endocrinologia & Metabologia. 2006;50(4):694–704. [PubMed]
Fuchs R, Faillace M, Allen M, Phipps R, Miller L, Burr D. Biphosphonates do not alter the rate of secondary mineralization. Bone. 2011;49(4):701–705. [PubMed]
Muschitz C, Kocijan R, Fahrleitner-Pammer A, Lung S, Resch H. Antiresorptives overlapping ongoing teriparatide treatment result in additional increases in bone mineral density. Journal of Bone and Mineral Research. 2013;28(1):196–205. [PubMed]
Fritsch A, Hellmich C. ‘Universal’ microstructural patterns in cortical and trabecular, extracellular and extravascular bone materials: micromechanics-based prediction of anisotropic elasticity. Journal of Theoretical Biology. 2007;244(4):597–620. [PubMed]
Bonewald L, Johnson M. Osteocytes, mechanosensing and Wnt signaling. Bone. 2008;42(4):606–615. [PMC free article] [PubMed]
Bonewald L. The amazing osteocyte. Journal of Bone and Mineral Research. 2011;26(2):229–238. [PMC free article] [PubMed]
Jacobs C, Temiyasathit S, Castillo A. Osteocyte mechanobiology and pericellular mechanics. Annual Review of Biomedical Engineering. 2010;12:369–400. [PubMed]
Moroz A, Crane M, Smith G, Wimpenny D. Phenomenological model of bone remodeling cycle containing osteocyte regulation loop. Biosystems. 2006;84(3):183–190. [PubMed]
Rieger R, Hambli R, Jennane R. Modeling of biological doses and mechanical effects on bone transduction. Journal of Theoretical Biology. 2011;274(1):36–42. [PubMed]
Hambli R, Rieger R. Physiologically based mathematical model of transduction of mechanobiological signals by osteocytes. Biomechanics and Modeling in Mechanobiology. 2012;11(1-2):83–93. [PubMed]
Graham J, Ayati B, Holstein S, Martin J. The role of osteocytes in targeted bone remodeling: a mathematical model. PLoS ONE. 2013;8(5):1–10. ):e63884, [PMC free article] [PubMed]
Klein-Nulend J, Van der Plas A, Semeins C, Ajubi N, Frangos J, Nijweide P, Burger E. Sensitivity of osteocytes to biomechanical stress in vitro. FASEB Journal. 1995;9(5):441–445. [PubMed]
Klein-Nulend J, Bacabac R, Bakker A. Mechanical loading and how it affects bone cells: the role of the osteocyte cytoskeleton in maintaining our skeleton. European Cells and Materials. 2012;24:278–291. [PubMed]
Ashman R, Cowin S, VanBuskirk W, Rice J. A continuous wave technique for the measurement of the elastic properties of cortical bone. Journal of Biomechanics. 1984;17(5):349–361. [PubMed]
Cowin S, Mehrabadi M. The structure of the linear anisotropic elastic symmetries. Journal of the Mechanics and Physics of Solids. 1992;40(7):1459–1471.

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s