Modeling osteoporosis to design and optimize pharmacological therapies comprising multiple drug types

Abstract
Editor's evaluation
eLife digest
Introduction
Results
Discussion
Methods
Appendix 1
Appendix 2
Appendix 3
Data availability
References
Article and author information
Metrics

Abstract

For the treatment of postmenopausal osteoporosis, several drug classes with different mechanisms of action are available. Since only a limited set of dosing regimens and drug combinations can be tested in clinical trials, it is currently unclear whether common medication strategies achieve optimal bone mineral density gains or are outperformed by alternative dosing schemes and combination therapies that have not been explored so far. Here, we develop a mathematical framework of drug interventions for postmenopausal osteoporosis that unifies fundamental mechanisms of bone remodeling and the mechanisms of action of four drug classes: bisphosphonates, parathyroid hormone analogs, sclerostin inhibitors, and receptor activator of NF-κB ligand inhibitors. Using data from several clinical trials, we calibrate and validate the model, demonstrating its predictive capacity for complex medication scenarios, including sequential and parallel drug combinations. Via simulations, we reveal that there is a large potential to improve gains in bone mineral density by exploiting synergistic interactions between different drug classes, without increasing the total amount of drug administered.

Editor's evaluation

The authors have developed a mathematical framework of drug interventions for postmenopausal osteoporosis using bisphosphonates, parathyroid hormone analogs, romosozumab, and denosumab. After calibrating and validating the model, authors demonstrated a predictive ability for complex clinical scenarios including sequential and parallel drug combinations. These data may be of great help in clinical practice.

https://doi.org/10.7554/eLife.76228.sa0

eLife digest

Our bones are constantly being renewed in a fine-tuned cycle of destruction and formation that helps keep them healthy and strong. However, this process can become imbalanced and lead to osteoporosis, where the bones are weakened and have a high risk of fracturing. This is particularly common post-menopause, with one in three women over the age of 50 experiencing a broken bone due to osteoporosis.

There are several drug types available for treating osteoporosis, which work in different ways to strengthen bones. These drugs can be taken individually or combined, meaning that a huge number of drug combinations and treatment strategies are theoretically possible. However, it is not practical to test the effectiveness of all of these options in human trials. This could mean that patients are not getting the maximum potential benefit from the drugs available.

Jörg et al. developed a mathematical model to predict how different osteoporosis drugs affect the process of bone renewal in the human body. The model could then simulate the effect of changing the order in which the therapies were taken, which showed that the sequence had a considerable impact on the efficacy of the treatment. This occurs because different drugs can interact with each other, leading to an improved outcome when they work in the right order.

These results suggest that people with osteoporosis may benefit from altered treatment schemes without changing the type or amount of medication taken. The model could suggest new treatment combinations that reduce the risk of bone fracture, potentially even developing personalised plans for individual patients based on routine clinical measurements in response to different drugs.

Introduction

Osteoporosis, a disease characterized by porous bone prone to fractures, affects hundreds of millions of people worldwide (Cooper and Ferrari, 2019; Hernlund et al., 2013). Most recent estimates place the global annual incidence of bone fragility fractures at 9 million in the year 2000 (Cooper and Ferrari, 2019); projections for the year 2050 suggest between 7 and 21 million annual hip fractures (Gullberg et al., 1997). Osteoporosis-associated bone fractures lead to disabilities, pain, and increased mortality (Cooper and Ferrari, 2019). In the United States, medical cost for osteoporosis, including inpatient, outpatient, and long-term care costs, has been estimated at US$17 billion in 2005 (Burge et al., 2007); in the European Union, the total cost of osteoporosis, including pharmacological interventions and loss of quality-adjusted life-years (QALYs), is projected to rise from about €100 billion in 2010 to €120 billion in 2025 (Odén et al., 2013).

Osteoporotic bone is the consequence of an imbalance of continuous bone resorption and bone formation, which—under close to homeostatic conditions—has the function to remove microfractures and renew the structural integrity of bone. Postmenopausal women are particularly at risk of osteoporosis: the rapid decline of systemic estrogen levels after menopause and other aging-related effects such as increased oxidative stress contribute to or drive the development of osteoporosis (Riggs et al., 1998; Manolagas, 2010). Moreover, osteoporosis can be a sequela of diseases affecting bone metabolism and remodeling such as primary hyperparathyroidism or gastrointestinal diseases (Painter et al., 2006). Osteoporosis can also be a side effect of treatments for other diseases; as a prime example, glucocorticoid administration is the most common cause of secondary osteoporosis (Weinstein, 2012). Over the last decades, an array of different osteoporosis treatments have emerged, from simple dietary supplementations such as calcium and vitamin D to specialized drugs targeting bone-forming and -resorbing cells and related signaling pathways (Tu et al., 2018). This entails a plethora of different medication options, including a large number of possible dosing schemes and combinations of drugs, administered in sequence or in parallel. Due to the huge number of such treatment schemes and the required time from study inception to completion, very few of them have been clinically tested so far when compared to the total number of available options.

Concomitant with the development of new osteoporosis drugs, mathematical and biophysical modeling approaches capturing bone-related physiology have advanced our quantitative understanding of the biological principles governing bone mineral metabolism, bone turnover, and development of osteoporosis. Pioneering work by Lemaire et al., 2004 describes the dynamics of bone-forming and -resorbing cell populations coupled through signaling pathways and could qualitatively reproduce the effects of senescence, glucocorticoid excess, and estrogen and vitamin D deficiency on bone turnover. Since then, compartment-based descriptions of the mineral metabolism, bone-forming and -resorbing cell populations, and related signaling factors have elucidated the role of essential regulatory mechanisms underlying mineral balance and bone turnover (Komarova et al., 2003; Lemaire et al., 2004; Pivonka et al., 2008; Pivonka et al., 2010; Peterson and Riggs, 2010; Zumsande et al., 2011; Schmidt et al., 2011; Graham et al., 2013; Tanaka et al., 2014; Komarova et al., 2015; Berkhout et al., 2015). Coarse-grained as well as detailed spatially extended descriptions of bone geometry have also addressed the effects of mechanical forces and the propagation of the multicellular units responsible for bone turnover (Ryser et al., 2009; Buenzli et al., 2011; Scheiner et al., 2013; Buenzli et al., 2014; Pivonka et al., 2013), as well as the influence of secondary diseases such as multiple myeloma (Ayati et al., 2010). Detailed models of bone remodeling and calcium homeostasis have become versatile and widely used tools in hypothesis testing, such as the seminal model by Peterson and Riggs, 2010, which includes submodels for various organs such as gut, kidney, and the parathyroid gland. Pharmacokinetic and pharmacodynamic (PK/PD) models of therapeutic interventions have mostly focused on capturing the mechanisms of action of a single or a few drugs and testing their dosing regimens (Marathe et al., 2008; Marathe et al., 2011; Ross et al., 2012; Scheiner et al., 2014; Eudy et al., 2015; Lisberg et al., 2017; Martínez-Reina and Pivonka, 2019; Zhang and Mager, 2019). Recent modeling efforts have also started addressing the effects of drug combinations on bone-forming and -resorbing cells, pointing out the need for corresponding model frameworks to include clinically relevant variables like bone mineral density (BMD) and bone turnover biomarkers (BTMs) (Lemaire and Cox, 2019), as well as combination therapies of physical exercise and drug treatment (Lavaill et al., 2020). An integrated mathematical framework for multiple drugs, which can also be used to quantitatively predict the effects of drug combinations in sequence and in parallel is not yet available.

Building on established mechanisms of bone turnover, we here present a quantitative model of bone turnover and postmenopausal osteoporosis treatment, unifying the description of multiple classes of drugs with different mechanisms of action, namely, bisphosphonates, parathyroid hormone (PTH) analogs, sclerostin antibodies, and receptor activator of NF- $κ$ B ligand (RANKL) antibodies. We calibrate the model using published population-level data from several clinical trials and assess its ability to predict the outcome of previously conducted clinical studies based on the medication scheme alone. We then use the model to demonstrate how medication schemes involving drug combinations can be optimized for a given medication load and discuss future model extensions.

Mechanisms of bone turnover and its regulation

Our model is based on a small set of key principles of bone turnover, which we briefly recapitulate here (Figure 1). As a composite tissue comprising hydroxyapatite, collagen, other proteins, and water (Boskey, 2013), bone is constantly turned over to renew its integrity and remove microdamage, at an average rate of about 4% per year in cortical bone and about 30% per year in trabecular bone (Manolagas, 2000).

Figure 1

Download asset Open asset

Schematic of the osteoporosis model describing the cell dynamics and signaling pathways within a ‘representative bone remodeling unit (BRU)’.

Regulatory interactions between different model components are indicated by colored boxes (see legend). TGFβ, transforming growth factor beta; BMP, bone morphogenetic protein; PDGF, platelet-derived growth factor; IGF, insulin-like growth factor; FGF, fibroblast growth factor.

Bone-resorbing and -forming cells

Bone resorption is performed by osteoclasts, multinucleated cells formed through the differentiation and fusion of their immediate precursors (pre-osteoclasts), which are derived from pluripotent hematopoietic stem cells via the myeloid lineage (Boyce and Xing, 2008). Osteoclasts attach to bone tissue and resorb it through the secretion of hydrogen ions and bone-degrading enzymes (Fuller and Chambers, 1995), which leads to the release of minerals and signaling factors stored in the bone matrix. New bone is formed by osteoblasts, a cell type derived from mesenchymal stem cells via several intermediate states that give rise to pre-osteoblasts and finally osteoblasts (Eriksen, 2010). Groups of osteoblasts organize into cell clusters (osteons) and collectively lay down an organic matrix (osteoid), which subsequently becomes mineralized over the course of months. Osteoblasts that are enclosed in the newly secreted bone matrix become osteocytes, nondividing cells with an average life span of up to several decades. Osteoclasts and osteoblasts organize into spatially defined local clusters termed ‘bone remodeling units’ (BRUs) (Figure 1), in which osteoblasts replenish the bone matrix previously resorbed by osteoclasts with a delay of several weeks. In cortical bone, the outer protective bone layer, BRUs migrate as a whole in ‘tunnels,’ whereas within the inner cancellous bone, BRUs propagate on the surfaces of the trabeculae, renewing the bone matrix in the process (Eriksen, 2010).

Signaling pathways

The differentiation and activity of osteoclasts and osteoblasts are regulated through several signaling pathways and hormones; recent reviews provide comprehensive descriptions of the various pathways (Siddiqui and Partridge, 2016). Osteoclast formation and activity are prominently regulated by RANKL and macrophage colony-stimulating factor (M-CSF) synthesized by bone marrow stromal cells. RANKL binds to receptor activator of NF- $κ$ B (RANK) on osteoclast precursors and promotes their differentiation into mature osteoclasts; osteoprotegerin (OPG) acts as a decoy receptor for RANKL and thus inhibits bone resorption (Boyce and Xing, 2008; Clarke, 2008). When laying down new bone, osteoblasts store signaling factors in the bone matrix, including transforming growth factor beta (TGFβ), bone morphogenetic protein (BMP), insulin growth factors (IGFs), platelet-derived growth factor (PDGF), and fibroblast growth factors (FGFs) (Solheim, 1998). Upon bone resorption, these factors are released and regulate cell fates and activity of osteoblasts and osteoclasts, thereby coupling bone resorption and formation (Houde et al., 2009; Eriksen, 2010). Osteocytes secrete sclerostin, a Wnt inhibitor interfering with extracellular binding of Wnt ligands (Li et al., 2005). Sclerostin inhibits bone formation and promotes resorption via downregulation of osteoblastogenesis and upregulation of osteoclastogenesis (Delgado-Calle et al., 2017; Maré et al., 2020). Since bone also acts as a mineral reservoir for the body, regulators of calcium homeostasis such as PTH and vitamin D also strongly affect the balance of bone formation and resorption alongside the intestinal absorption and renal reabsorption of calcium (Mundy and Guise, 1999).

Estrogen

The sex hormone estrogen inhibits bone resorption by inducing apoptosis of osteoclasts (Kameda et al., 1997) and lowering circulating sclerostin levels (Mödder et al., 2011). The rapid decline of estrogen levels after menopause is one known cause of postmenopausal osteoporosis (Riggs et al., 1998).

Results

Model overview

The primary purpose of our model is to provide an efficient representation of bone turnover on multiple time scales from weeks to decades that allows for the quantitative description of drug interventions. Of particular interest are the consequences of pharmacological therapies on long-term dynamics of the BMD in specific bone sites and biochemical markers of bone formation and resorption. To this end, we considered a minimal set of physiologically relevant dynamic components (Figure 1) that are sufficient to capture a large range of clinically observed population-level data on drug interventions. Thus, our model describes a ‘representative BRU’ that abstracts from the vast set of intricate regulatory mechanisms underlying calcium homeostasis or the complex bone geometry.

Our model comprises the following dynamic components to describe the bone turnover through a representative BRU: cell densities of (i) pre-osteoclasts, (ii) osteoclasts, (iii) pre-osteoblasts, (iv) osteoblasts, (v) osteocytes, (vi) sclerostin concentration, (vii) total bone density, and (viii) bone mineral content (BMC). The BMD is given by the product of bone density and BMC. Osteoblasts and osteoclasts can undergo apoptosis and are derived from pre-osteoblasts and pre-osteoclasts, respectively, with differentiation rates that depend on regulatory factors such as estrogen and sclerostin (Figure 1). Pre-osteoblasts and pre-osteoclasts are formed at constant rates and undergo apoptosis. These progenitor populations provide a dynamic reservoir for rapid differentiation and activation of osteoblasts and osteoclasts, respectively, which can be temporarily depleted if stimulated by a drug intervention. Osteocytes are derived from osteoblasts and provide a source of sclerostin, which has a regulatory effect on osteoblasts, osteoclasts, and thus, bone density change. The gain and loss rates of bone density are proportional to the density of osteoblasts and osteoclasts, respectively. The BMC has a steady state whose level can be temporarily shifted through drug administration, effectively accounting for more complex underlying dynamics such as promotion of secondary mineralization. All rates of cell formation, differentiation, apoptosis, and bone formation and resorption generally depend on the concentration of sclerostin, estrogen, and a ‘resorption signal.’ These dependencies also implicitly account for regulation of bone remodeling via other routes, for example, the RANK–RANKL–OPG pathway. The effects of aging and the onset of menopause are represented through an age-dependent serum estrogen concentration, which has been determined from the literature (Sowers et al., 2008; Appendix 1). The resorption signal corresponds to the melange of signaling factors stored in the bone matrix. Therefore, its release is proportional to the rate of bone resorption. The serum concentration of BTMs such as the resorption marker C-terminal telopeptide (CTX), the formation markers procollagen type 1 amino-terminal propeptide (P1NP), and bone-specific alkaline phosphatase (BSAP) were identified with elementary functions of the bone resorption and formation rates in the model (Appendix 1).

We extended this core model of long-term bone turnover by a dynamic description of the mechanisms of action of several drug classes used in osteoporosis treatment: RANKL antibodies (denosumab), sclerostin antibodies (romosozumab), bisphosphonates (alendronate and others), and PTH analogs (teriparatide) (Appendix 2). We also included blosozumab, another sclerostin inhibitor, which was investigated in osteoporosis trials but not approved for osteoporosis treatment at the time the present work was conducted. PTH is known to exert anabolic or catabolic effects depending on whether administration is intermittent or continuous (Tam et al., 1982; Hock and Gera, 1992); PTH description in our model is restricted to the anabolic administration regimes relevant for osteoporosis treatment. A schematic overview of all model components, mechanisms, and regulatory interactions is given in Figure 1; a detailed formal description of the model and its extensions is provided in Appendix 1 and Appendix 2.

Capturing clinical study results with the model

The model and the corresponding medication modules rely on an array of physiological parameters (rates of cell formation, differentiation and death, concentration thresholds for signaling activity, medication efficacies and half-lives, etc.) many of which are not directly measureable. However, clinical measurements on physiological responses to medications with different mechanisms of action provide a wealth of indirect information about time scales of bone turnover and regulatory feedbacks. We calibrated the model using published clinical data from various seminal studies on both (i) long-term BMD age dependence and (ii) the response of BMD and BTMs to the administration of different drugs (see Appendix 3—table 2 for a comprehensive list of data sources). Although BMD constitutes the major target variable of our model, the dynamics of BTM concentrations carry important complementary information about the mode of action of the administered drugs (antiresorptive, anabolic, and combinations) that crucially informs the calibration procedure. To allow the model to capture the effects of medications as physiologically sensible modulations of the age-dependent bone mineral metabolism, we created hybrid datasets each of which comprised both aging-related BMD changes and the response to a treatment (see ‘Methods’ and Appendix 1—figure 1D).

We then determined a single set of model parameters through a simultaneous fit of the free 31 model parameters to capture a specified set of hybrid aging/treatment datasets containing different drug responses (Appendix 3). Without constraining the average rate of skeletal bone turnover, model calibration yielded an inferred value of about 6% per year on average, of the same order of values reported for cortical bone, which constitutes about 75% of the skeleton (Manolagas, 2000). The model was able to capture the BMD and BTM dynamics across all calibration datasets with remarkable accuracy (Appendix 3—figure 1), despite the model’s structural simplicity. To quantify the goodness of the fit, we computed the mean absolute percentage error (MAPE) between model simulations and clinical data; the MAPE for BMD was consistently below 1% for all calibration datasets (Appendix 3—table 3), indicating an excellent agreement between model and data. The qualitative behavior of BTMs (i.e., the direction of their excursions from baseline) was captured correctly in all calibration datasets, indicating an adequate description of the drugs’ mode of action in the model; relative deviations in the total magnitude of BTM excursions observed for some datasets were mostly due to slight offsets in the timing of peaks and troughs and low absolute values of the respective BTM concentrations, as highlighted by comparing different goodness measures (Appendix 3 and Appendix 3—table 3).

After obtaining the reference parameter set, we sought to validate the calibrated model by assessing its ability to predict the effects of drug dosing schemes that had not been used for calibration. Model validation included complex sequential and parallel drug combinations and therefore challenged the model to predict the effects of treatment schemes beyond those used in calibration (Appendix 3—table 2). To this end, the model received only drug dosing information used in the respective clinical trials but was not informed by BMD or BTM measurements, which instead it had to predict. With the single set of previously determined parameters, the model showed a remarkable capacity to quantitatively forecast the effects of a multitude of medication schemes, both during treatment and follow-up periods (Figure 2, Figure 2—figure supplement 1). Even in scenarios including sequential treatments with up to three different drug types and parallel treatments with two different drugs, respectively, the model was able to predict the complex progression of both BMD and biomarker levels with a high degree of accuracy (Figure 2). Across all validation datasets, MAPEs for BMD were consistently below 1.5% (Appendix 3—table 3), indicating an excellent predictive capacity of the model. In summary, this validation provided a strong corroboration of the model’s capacity to capture the physiological dynamics of bone turnover and the mechanisms of action of various drugs relevant to osteoporosis treatment using a single set of model parameters.

Figure 2 with 1 supplement see all

Download asset Open asset

With a single set of parameters, the calibrated model can quantitatively predict the effects of various drugs in different dosing regimens, alone and in combination.

(A) Comparison of simulated total hip bone mineral density (BMD, black curves) and clinical data (dots), including aging behavior (green dots) and treatment behavior (black dots) of various sequential drug treatments, including denosumab, romosozumab, alendronate, and teriparatide. Hybrid aging/treatment datasets were created combining data from Looker et al., 1998 (aging dataset, green dots in panel A; in total $N = 3251$ subjects 20 years and older), as well as Recknor et al., 2015 (blosozumab 180 mg Q2W: $N = 25$ ), McClung et al., 2018 (placebo/deno.: $N = 18$ , alendro./romo./deno.: $N = 21$ ), and Leder et al., 2015 (deno./teri.: $N = 27$ , teri. + deno./deno.: $N = 23$ ) (treatment datasets, black dots in panels A and B) as indicated, see ‘Methods.’ (B) Zoom into the treatment regions shown in panel (A) including BMD (black) and baseline changes of the bone resorption marker C-terminal telopeptide (CTX, red) and the bone formation marker procollagen type 1 amino-terminal propeptide (P1NP, blue). Colored bars above the plots indicate the medication scheme (see legend). Data points show population averages; average types and error bar types as reported in the respective original publication. In both panels, BMD is displayed as a fraction of its value at $t_{0} = 25$ years.

Testing alternative treatment schemes

Having established the predictive capacity of the model for the considered medications, we aimed to utilize the model to study and optimize hypothetic drug dosing regimens. As an example, we considered a sequential treatment with three drugs of different types: the bisphosphonate alendronate, the sclerostin inhibitor romosozumab, and the RANKL inhibitor denosumab. In a clinical trial reported by McClung et al., 2018, the sequence alendronate (70 mg per week for 1 year), followed by romosozumab (140 mg per month for 1 year), followed by denosumab (60 mg per 6 months for 1 year) had been studied (Figure 2). However, in principle there are six different sequences in which these drugs can be administered: ARD, ADR, DAR, DRA, RAD, and RDA (A: alendronate; R: romosozumab; D: denosumab). A priori, it is not obvious whether synergistic or antagonistic interactions between these drugs and the physiological state in which they leave the patient may lead to a differential short- and long-term evolution of BMD and biomarkers between different medication sequences. Probing all six sequences in a clinical trial would present a time- and resource-consuming endeavor and inevitably expose part of the study population to suboptimal treatment schemes. Instead, we probed these different treatment options using the present model (Figure 3A). To assess the predicted clinical success of different sequences, we compared two clinically relevant outcomes across different schemes: (i) the maximum achieved BMD increase (as compared to baseline at treatment start) irrespective of when it occurred (Figure 3B) and (ii) the residual long-term effects of treatment on BMD as monitored by the relative BMD 10 years after treatment end (Figure 3C) .

Figure 3

Download asset Open asset

The model predicts differential outcomes for different sequences of the same drugs at constant total medication load.

(A) Simulated progression of bone mineral density (BMD) and C-terminal telopeptide (CTX) and procollagen type 1 amino-terminal propeptide (P1NP) concentrations for different sequences (columns) of the three drugs denosumab (D), alendronate (A), and romosozumab (R) as indicated. Simulated treatment starts at age 67. The total amount of drug administered is identical among columns. Clinical results on the sequence ARD (column 5) were reported in McClung et al., 2018, see also Figure 2. (B) Maximum simulated BMD (relative to baseline at treatment start) achieved during the course of treatment for different drug sequences. (C) Simulated BMD 10 years after treatment end (relative to baseline at treatment start) for different drug sequences.

Indeed, we found that the outcomes of different medication sequences were markedly different despite the same total amount of drug administered (Figure 3A). Some sequences (such as ARD and RAD) reached a considerably higher maximum BMD during the course of the simulated treatment, which allowed us to rank treatments according to maximum BMD gain (Figure 3B). Notably, while some sequences were superior to others as measured by the maximum BMD increase during treatment, they performed markedly worse (as compared to, e.g., DRA and RDA) with regard to long-term BMD evolution as predicted by model simulations (Figure 3C). This behavior suggests that short-term BMD gains may be limited as a proxy for the clinical benefit of a treatment as a whole. Within our modeling scheme, the explanation for this behavior is found in differing ‘rebound’ effects after treatment end: simulated drug-mediated inhibition of osteoclastogenesis leads to a build-up of an undifferentiated osteoclast precursor pool. After treatment end, this precursor pool becomes licensed to differentiate and rapidly gives rise to a large active osteoclast population, leading to accelerated resorption of the bone matrix that had been built up during treatment. In this paradigm, specific drug sequences lead to an attenuation of this effect, for example, by enhancing osteoclast apoptosis during such a ‘rebound’ phase, thereby modulating bone turnover in the long run.

In summary, our model analysis suggests considerable potential in the improvement of dosing regimens and drug sequencing in osteoporosis treatment, especially combination therapies, to achieve an optimal effect for a given medication load. These improvements are possible because the mechanisms of action of one drug may act either favorably or adversely on the state of the bone mineral metabolism left behind by the preceding treatment with another drug.

Discussion

Treatment of osteoporosis is complex, expensive, and in many circumstances opinion-based. With bone physiology as our guiding principle, we have introduced a mathematical modeling framework that can quantitatively capture and predict the progression of osteoporosis in postmenopausal women with and without medical therapy. Our model is built on a small set of essential mechanisms of bone turnover. The effectivity of this approach suggests that—despite the complexity of the bone mineral metabolism—the dynamics relevant for osteoporosis medications can be condensed into only a few components. These components describe the biology of osteoblasts, osteoclasts, and osteocytes, as well as their precursor cell populations and a few essential regulatory feedbacks through hormones and signaling factors such as estrogen, sclerostin, and bone-matrix-derived factors.

The general nature of the model allowed us to capture the BMD and BTMs of a multitude of clinical treatment studies. Notably, the model can also predict the effects of a broad range of drug dosing regimens and complex drug combination therapies beyond those used for model development. This corroborates the model’s predictive capacity, supporting its use for the design of future clinical trials. However, it is important to note that some parameters (e.g., concentration thresholds for signaling factors) were inferred through the model calibration procedure from BMD and BTM dynamics alone. Hence, when used as a predictive tool, general quantitative limitations of the model have to be considered, especially when extrapolating into extreme dosing regimens, dosing frequencies, or age regions beyond the validated ones.

It is of clinical relevance that exemplary model predictions suggest a large potential for the development of optimized combination therapies involving different drug types and treatment schemes. These may range from a simple rearrangement of a sequence of drugs at given total drug doses (as shown in this article) to complex interwoven or cyclic administration schemes that exploit synergistic effects between different medication types. Notably, model simulations extrapolating the long-term BMD development after treatment end suggest that medication schemes eliciting a rapid BMD increase are not necessarily accompanied by a sustained elevation of the BMD. Instead, some initially successful treatment schemes may lead to a ‘rebound’ effect of accelerated bone loss after treatment end, a prediction that cautions against using short-term BMD increases as the sole proxy for treatment success. Such extrapolations into follow-up periods long after treatment end, which are mostly inaccessible to clinical studies, highlight the potential role of the model in considering long-term treatment success when optimizing treatment schemes.

In our research, we have focused on postmenopausal osteoporosis, the most widespread type of osteoporosis. However, the generic manner in which the model represents bone remodeling and the effects of medications renders it a general platform for the study of treatments that can be adapted to other types of primary and secondary osteoporosis. The modular nature of the model enables future extensions; besides additional medication types, these may include the effects of comorbidities that elicit osteoporosis or interact with it (such as primary and secondary hyperparathyroidism), medications that contribute to osteoporosis (such as glucocorticoid therapy), lifestyle-dependent factors such as smoking and alcohol consumption, the effects of dietary supplementation of osteoporosis treatment through calcium and vitamin D and effects of microgravity on bone, as experienced by astronauts on extended missions in space. Physical activity is another important contributor to bone remodeling, which we have not considered here. Detailed modeling approaches involving biomechanical feedback suggest synergistic effects between drug treatment of osteoporosis and physical activity (Lavaill et al., 2020). Such results call for a further exploration of integrated approaches to osteoporosis therapy combining pharmacological treatment and lifestyle adjustments.

Clearly, the goal of osteoporosis therapy is the reduction of fracture risk during and after therapy. While BMD has a prime role in the evaluation of osteoporosis therapies and can be measured rather easily using dual-energy x-ray absorptiometry (DXA), its relationship to fracture risk is complex. Fracture risk calculations used in clinical practice also involve demographic and lifestyle-related factors while mostly relying on BMD point measurements (Kanis et al., 2009). However, the quantitative associations between BMD, age, and fracture risk reported in many studies (Kanis et al., 2001; Berger et al., 2009; Austin et al., 2012; Krege et al., 2013; Black et al., 2018; Ensrud et al., 2022) can be used to create statistical models that may relate entire BMD time courses to a patient’s fracture risk. Combining such statistical models with the physiology-based model presented here would enable to optimize therapies directly for a minimized long-term fracture risk instead of maximized BMD gain. Thus, our model can serve as a quantitative starting point for the forecast of pharmacological therapies of osteoporosis but also highlights the role of mechanistic mathematical descriptions in understanding the biological principles of drug action.

Methods

Hybrid aging/treatment datasets

To create hybrid aging/treatment datasets, we merged a dataset comprising the BMD age dependence from Looker et al., 1998 with different clinical study datasets containing the BMD response to various medications (Appendix 3—table 2). The aging dataset from Looker et al., 1998 consisted of mean total femur BMD measurements in non-Hispanic white, non-Hispanic black, and Mexican American women, reported in 10-year age bins ranging from 20 to 80 years and older. We used bin averages as proxy BMD indicators for the center of the respective age window (Appendix 1—figure 1B). Rescaling the reported means for the three ethnic groups to their value for the earliest age bin revealed that relative changes in BMD were remarkably consistent among ethnic groups (Appendix 1—figure 1C) despite differing absolute baselines. Therefore, and since the model only addresses relative BMD changes, we resorted to the dataset with the largest underlying study population for calibration, which was the dataset comprising the non-Hispanic white female study population. Datasets on the response to medications from clinical trials on romosozumab, blosozumab, denosumab, alendronate, and teriparatide consisted of study population averages of total hip BMD and serum concentrations of one or more BTMs (CTX, P1NP, BSAP) during the treatment, and if available, during a follow-up period. Reported study population averages on the respective quantities were digitized directly from the data figures in the corresponding publications (Appendix 3—table 2).

To merge aging and treatment datasets, the BMD from treatment datasets was rescaled such that the BMD baseline at treatment start corresponds to the linearly interpolated age-dependent BMD at treatment start. The treatment start was placed at the average age of the study population upon study start (rounded to full years) as reported in the respective publication (Appendix 1—figure 1D). BTM measurements were normalized to baseline values.

Appendix 1

Model of long-term bone remodeling and osteoporosis

The model of bone remodeling underlying the present description of osteoporosis and its treatment was built with the aim to provide a minimal set of dynamic components necessary to quantitatively capture population averages of both aging-related changes in bone turnover and the response to osteoporosis medications with different mechanisms of action. The model is partitioned into a core model describing the patient physiology (Appendix 3—table 1) and separate extensions for different drug classes (Appendix 3—table 2).

The model is compartment-based and describes average cell densities, bone densities, and average concentrations of signaling factors within a ‘representative bone remodeling unit (BRU)’, that is, a fictitious BRU corresponding to an average over the considered bone type. All model variables are treated as nondimensional quantities, that is, the model only addresses relative changes in all variables, which simplifies the model structure and reduces the number of free parameters.

Model description

The model describes the dynamics of the cell densities of pre-osteoclasts ( $ρ_{C^{*}}$ ), pre-osteoblasts ( $ρ_{B^{*}}$ ), osteoclasts ( $ρ_{C}$ ), osteoblasts ( $ρ_{B}$ ), osteocytes ( $ρ_{Y}$ ), as well as functional sclerostin levels ( $s$ ), total bone density ( $ρ_{b}$ ), and BMC ( $c_{b}$ ). Functional estrogen levels ( $e$ ) are provided as an explicitly age-dependent function, described further below. Moreover, the model includes a ‘resorption signal’ ( $r$ ) corresponding to the composite concentration of bone matrix-derived signaling factors (TGFβ, BMPs, PDGF, IGFs, and FGFs) released upon bone resorption. Sclerostin, estrogen, and the resorption signal act as regulatory factors that modulate the rates of cell proliferation, differentiation, and apoptosis, as well as their bone-forming and -resorbing activity; for notational convenience, they are summarized in the vector $ϕ = (s, e, r)$ . Rates that depend on $ϕ$ are denoted with a tilde.

The dynamics of the cell density variables is given by

\begin{matrix} {\dot{ρ}}_{C^{*}} & = {\tilde{μ}}_{C^{*}} (ϕ) - f_{C^{*} \to C} (ϕ) - {\tilde{η}}_{C^{*}} (ϕ) ρ_{C^{*}}, \end{matrix}

\begin{matrix} {\dot{ρ}}_{B^{*}} & = {\tilde{μ}}_{B^{*}} (ϕ) - f_{B^{*} \to B} (ϕ) - {\tilde{η}}_{B^{*}} (ϕ) ρ_{B^{*}}, \end{matrix}

\begin{array}{ll} {\dot{ρ}}_{C} & = f_{C^{*} \to C} (ϕ) - {\tilde{η}}_{C} (ϕ) ρ_{C}, \end{array}

\begin{matrix} {\dot{ρ}}_{B} & = f_{B^{*} \to B} (ϕ) - {\tilde{η}}_{B} (ϕ) ρ_{B} - f_{B \to Y} (ϕ), \end{matrix}

\begin{aligned} {\dot{ρ}}_{Y} & = f_{B \to Y} (ϕ) - {\tilde{η}}_{Y} (ϕ) ρ_{Y}, \end{aligned}

where dots denote time derivatives, ${\tilde{μ}}_{x}$ denotes the formation rate of cell population $x$ , and ${\tilde{η}}_{x}$ denotes its apoptosis rate. Differentiation or conversion from one cell type to another is described by the absolute differentiation rates $f_{x \to y}$ of cell population $x$ giving rise to cells of type $y$ ; they have the generic form $f_{x \to y} (ϕ) = {\tilde{ω}}_{x} (ϕ) ρ_{x},$ with ${\tilde{ω}}_{x}$ denoting the differentiation rate. All rates are functions of the regulatory factors $ϕ$ as indicated.

Sclerostin is synthesized by osteocytes; the dynamics of Sclerostin levels is given by

\dot{s} = {\tilde{α}}_{s} (ϕ) ρ_{Y} - {\tilde{κ}}_{s} (ϕ) s,

where ${\tilde{α}}_{s}$ and ${\tilde{κ}}_{s}$ denote the synthesis and degradation rate, respectively.

The dynamics of the total bone density (BD) follows:

\begin{array}{ll} {\dot{ρ}}_{b} & = b^{+} - b^{-}, \\ b^{+} & = {\tilde{λ}}_{B} (ϕ) ρ_{B}, \\ b^{-} & = {\tilde{λ}}_{C} (ϕ) ρ_{C}, \end{array}

where $b^{+}$ and $b^{-}$ are the absolute rates of bone formation and resorption, respectively, with ${\tilde{λ}}_{B}$ and ${\tilde{λ}}_{C}$ being the formation and resorption rates per unit osteoblast and osteoclast density, respectively.

The BMC follows the dynamics

{\dot{c}}_{b} = \tilde{γ} (ϕ) [{\tilde{c}}_{0} (ϕ) - c_{b}],

where ${\tilde{c}}_{0}$ is the steady-state homeostatic BMC, and $\tilde{γ}$ is the equilibration rate.

Based on the physiological foundations summarized in the Section 'Mechanisms of bone turnover and its regulation' in the main text, we now specify the functional form of all rates. Upregulation and downregulation through a regulatory factor with concentration $x$ are described by a multiplicative or additive contribution $g^{+} (x / x_{0})$ or $g^{-} (x / x_{0})$ , respectively, where x₀ is the threshold concentration at which half-effect is reached (EC₅₀) and where $g^{\pm}$ are monotonic, saturating functions of the Hill-type (Keener and Sneyd, 2009),

g^{+} (u) = \frac{u}{1 + u}, g^{-} (u) = \frac{1}{1 + u},

Using these conventions, the dependencies of rates on regulatory factors are given by

\begin{aligned} {\tilde{μ}}_{C^{*}} (ϕ) = 1, {\tilde{μ}}_{B^{*}} (ϕ) = 1, \end{aligned}

\begin{aligned} {\tilde{ω}}_{C^{*}} (ϕ) = g^{-} (\frac{e}{e_{C^{*}}}) g^{+} (\frac{s}{s_{C^{*}}}) ω_{C^{*}}, {\tilde{ω}}_{B^{*}} (ϕ) = g^{-} (\frac{s}{s_{B^{*}}}) ω_{B^{*}}, \end{aligned}

\begin{matrix} {\tilde{ω}}_{B} (ϕ) = ω_{B}, . \end{matrix}

\begin{aligned} {\tilde{η}}_{C^{*}} (ϕ) = 0, {\tilde{η}}_{B^{*}} (ϕ) = 0, \end{aligned}

\begin{aligned} {\tilde{η}}_{C} (ϕ) = [1 + ν_{C} g^{+} (\frac{e}{e_{C}}) g^{+} (\frac{r}{r_{C}})] η_{C}, {\tilde{η}}_{B} (ϕ) = η_{B}, \end{aligned}

\begin{matrix} {\tilde{η}}_{Y} (ϕ) = η_{Y}, . \end{matrix}

\begin{aligned} {\tilde{α}}_{s} (ϕ) = g^{-} (\frac{e}{e_{s}}), & {\tilde{κ}}_{s} (ϕ) = κ_{s}, \end{aligned}

\begin{aligned} {\tilde{λ}}_{B} (ϕ) = λ_{B} g^{-} (\frac{s}{s_{Ω}}) [1 + ν_{Ω} g^{+} (\frac{r}{r_{Ω}})], {\tilde{λ}}_{C} (ϕ) = λ_{C}, \end{aligned}

\begin{aligned} \tilde{γ} (ϕ) = γ, {\tilde{c}}_{0} (ϕ) = c_{0}, \end{aligned}

where rates without a tilde denote model parameters. A full list of parameters introduced here and their description is provided in Appendix 3—table 4. This regulatory scheme comprises a multitude of interactions, many of which are simplified effective representations of indirect molecular mechanisms (e.g., through the RANK–RANKL–OPG pathway as described below). In the model building process, the effects of additional regulatory elements (not presented here) were probed and found to be nonessential within the scope of the present modeling aim or nonidentifiable with regard to related parameter values. Specific choices for regulatory interactions were partly based on insights from animal and culture studies and are motivated as follows:

Equation 9: Pre-osteoclasts and pre-osteoblasts are produced at constant rates, a simplifying assumption reflecting the fact that the main function of these populations in the present context is to provide a dynamic reservoir for the rapid supply with active osteoclasts and osteoblasts, respectively.
Equation 10: Estrogen suppresses pre-osteoclast to osteoclast differentiation; a consequence of suppression of RANKL expression (Streicher et al., 2017). Sclerostin upregulates pre-osteoclast to osteoclast differentiation; a consequence of upregulation of RANKL expression and downregulation of OPG expression (Wijenayaka et al., 2011). Sclerostin downregulates osteoblastogenesis; a consequence of the inhibition of osteoblast differentiation mediated by bone morphogenetic protein 2 (BMP2) and Wnt3a and possibly other pathways (Winkler et al., 2005; Thouverey and Caverzasio, 2015).
Equation 11: Sclerostin downregulates osteoblast to osteocyte conversion (Atkins et al., 2011).
Equation 13: Estrogen and the resorption signal induce osteoclast apoptosis; estrogen has been reported to induce osteoclast apoptosis both directly and mediated by TGFβ (Kameda et al., 1997; Hughes et al., 1996); TGFβ has been shown to upregulate Bim, a member of the (pro-apoptotic) Bcl2 family (Houde et al., 2009).
Equation 15: Estrogen reduces sclerostin production (Mödder et al., 2011).
Equation 16: Sclerostin downregulates bone formation (Li et al., 2008; Atkins et al., 2011). The resorption signal upregulates bone formation; a consequence of TGFβ1 enhancing bone collagen synthesis (Rydziel et al., 1997b); furthermore, TGFβ1, skeletal BMPs, and IGF1 have been reported to inhibit collagenase-3 expression in osteoblasts (Rydziel et al., 1997b; Gazzerro et al., 1999; Rydziel et al., 1997a).

Estrogen concentration is described through its explicit age dependence. Clinical data reported by Sowers et al., 2008 were used to construct a function capturing the key features of age-dependent decline of serum estradiol levels:

e (t) = {\begin{cases} 1 & t < t_{e} \\ \frac{1}{1 + (t - t_{e}) / τ_{e}} & t \geq t_{e} \end{cases} .

where $t$ denotes time. The parameters $t_{e}$ and $τ_{e}$ denote the age at the onset of estradiol decline and a characteristic time scale of the decline, respectively. The time scale $τ_{e}$ was determined using a fit of the function to the data reported in Sowers et al., 2008 (Appendix 1—figure 1A, Appendix 3—table 4).

The resorption signal $r$ corresponds to the concentration of bone matrix-derived signaling factors released upon bone resorption. Assuming a release rate proportional to the bone resorption rate $b^{-}$ and first-order degradation, we consider a highly simplified dynamics of the type $\dot{r} = b^{-} - κ r$ , where $κ$ is an effective average degradation rate of the components of the resorption signal. Given that the time scale of degradation, $κ^{- 1}$ , is much shorter (minutes to hours) than the time scale of osteoclast formation and death (weeks), the instantaneous concentration can be approximated to always follow its steady state, $r \approx b^{-} / κ$ , which is proportional to the osteoclast density, $r \propto b^{-} \propto ρ^{C}$ , via Equation 7. Since the resorption signal acts as a regulator of bone formation and is rescaled by individual concentration thresholds (see Equation 9–Equation 17), the proportionality constant can be absorbed in these thresholds, which enables us to set $r = ρ^{C}$ . Thus, the resorption signal concentration is approximated by the osteoclast density, so that no additional dynamic variable is required.

Description of BMD, bone turnover rate, and BTMs

To compare the model output to clinical data, we relate model variables to clinical observables frequently measured in clinical trials such as BMD and established biomarkers of bone turnover. The BMD follows from the model state as the product of total bone density and BMC:

BMD = ρ_{b} c_{b} .

In our model, levels of BTMs such as the bone formation markers P1NP and BSAP and the bone resorption marker CTX are related to the rates $b^{+}$ and $b^{-}$ of bone formation and resorption, respectively (see Equation 7). Here, we relate the BTMs P1NP, BSAP, and CTX to bone turnover rates by power laws with marker-specific exponents:

\begin{array}{ll} θ_{B S A P} & = (b^{+})^{q_{B S A P}}, \\ θ_{P 1 N P} & = (b^{+})^{q_{P 1 N P}}, \\ θ_{C T X} & = (b^{-})^{q_{C T X}} . \end{array}

The exponents q_x are obtained as fit parameters using clinical trial data, as described further below.

Appendix 1—figure 1

Download asset Open asset

Parameterization of the aging behavior and creation of hybrid aging/treatment datasets for model calibration and validation.

(A) Age dependence of estradiol serum levels. Clinical data (dots) modified from Sowers et al., 2008. The curve shows a fit of the function given by Equation 18 to determine the parameter τ_e (Appendix 3—table 4). (B) Bone mineral density (BMD) age dependence for different ethnic groups as indicated. Data modified from Looker et al., 1998; reported age bin averages have been used to represent the center of the age bin. (C) BMD age dependence shown in panel (B), where all curves have been normalized to their earliest value (t₀ = 25y). (D) Schematic of how hybrid aging/treatment datasets were generated by merging the same aging dataset with different treatment datasets; for details, see ‘Methods.’

Appendix 2

Model extensions for medications

We include a dynamic description of several drug classes through separate model extensions, which depend on the functional drug concentration. The pharmacodynamic description is drug-specific and represents the individual mechanism of action of the respective drug class. For the pharmacokinetic description of each drug, we resort to simple first-order kinetics with drug-specific half-lives, which reduces the amount of model parameters. More detailed pharmacokinetic descriptions involve drug absorption and transfer between different body compartments, depending on the route of administration (oral, intravenous, or subcutaneous). However, simulations of the calibrated model demonstrate that first-order kinetics yields an effective approximation of the pharmacokinetic features essential to capture a drug’s long-term effects on bone remodeling, as suggested by comparisons of simulated and measured BTM concentrations (Figure 2, Appendix 3—figure 1). A patient’s systemic concentration of a medication is represented by an effective (dimensionless) variable $ψ$ that indicates the relative concentration of the medication. Typically, $ψ$ is given in multiples of a threshold that parameterizes the effect of the drug (such as EC₅₀)—the precise interpretation of $ψ$ depends on the model extension that describes the pharmacodynamics of the drug; see ‘Pharmacodynamics for specific medications’.

Pharmacokinetics

The pharmacokinetics of a drug $x$ that is administered in intervals of weeks or months is described by two parameters: the efficacy $E_{x}$ and the half-life $T_{x}$ . Given repeated administrations with doses $c_{1}, \dots, c_{n}$ at times $t_{1}, \dots, t_{n}$ , the efficacy-weighted concentration variable of the drug $x$ therefore follows the exponential kinetics

ψ_{x} (t) = E_{x} \sum_{i = 1}^{n} c_{i} 2^{- (t - t_{i}) / T_{x}} Θ (t - t_{i}),

where $Θ$ is the Heaviside function, defined by

Θ (t) = {\begin{cases} 0 & t < 0 \\ 1 & t \geq 0 \end{cases} .

Drugs that are administered more frequently (e.g., daily or weekly) are more efficiently captured in a quasi-continuous scheme. The dynamics of BMD and BTM levels is much more inert than such fast administration/degradation dynamics and is well-described by their effective average action. In this quasi-continuous scheme, the drug is considered to be administered at a given average rate for a specified amount of time, so that its concentration evolves according to the dynamic equation

\frac{d ψ_{x}}{d t} = E_{x} \sum_{i = 1}^{n} c_{i} Θ (t - t_{i}) Θ (t_{i}^{*} - t) - \frac{\ln 2}{T_{x}} ψ_{x},

with initial condition ${ψ_{x} (t) |}_{t \to - \infty} = 0$ , where c_i are doses per unit time, and where t_i and $t_{i}^{*}$ are the start and end times of a treatment period, respectively. (Numerically, the quasi-continuous scheme has the advantage that model simulations do not have to resort to extremely small integration time steps to capture the details of short-term drug degradation, which considerably improves runtime.)

Different drugs $x_{1}, x_{2}, \dots$ of the same class (sclerostin inhibitors, bisphosphonates, PTH analogs, etc.) are considered through an effective concentration equivalent that is the sum of the efficacy-weighted doses of different drugs:

Ψ (t) = \sum_{i} ψ_{x_{i}} (t),

where the $ψ_{x_{i}}$ are given by Equation 21 in the case of discrete doses and Equation 22 in the case of quasi-continuous dosing.

Pharmacodynamics for specific medications

We now introduce separate model extensions that embody the essential mechanisms of action of different drug classes.

RANKL antibodies

Denosumab is a monoclonal antibody (mAb) that binds with affinity to receptor activator of NF- $κ$ B ligand (RANKL) and blocks its interactions with RANK (Kostenuik et al., 2009), which, in turn, decreases osteoclast formation (Boyce and Xing, 2008). Moreover, denosumab has been suggested to enable increased bone tissue mineralization, which leads not only to a halt of the BMD decline but also to a long-term increase in BMD (Scheiner et al., 2014). We include these mechanisms in our model through a modification of the reference pre-osteoclast differentiation rate $ω_{C^{*}}$ to include a downregulation by denosumab and the steady-state BMC c₀ to include an upregulation through denosumab in Equation 10 and Equation 17, respectively:

\begin{array}{ll} ω_{C^{*}} & \to [1 - β_{C^{*}}^{r A b} g^{+} (Ψ_{r A b})] ω_{C^{*}}, \\ c_{0} & \to c_{0} + β_{b}^{r A b} g^{+} (Ψ_{r A b}), \end{array}

where $Ψ_{rAb}$ is the RANKL antibody concentration in multiples of the half-maximal effective concentration (EC₅₀), determined through Equation 21 and Equation 23. For simplicity, EC₅₀ for the regulation of both differentiation and mineralization are taken to be identical, which is justified a posteriori by showing its effectivity in approximating the drug action. The scaling factors $β_{C^{*}}^{rAb}$ and $β_{b}^{rAb}$ parameterize the respective maximum effect strength and are subject to the constraints $β_{C^{*}}^{r A b} < 1$ and $c_{0} + β_{b}^{r A b} < 1$ to ensure positive rates and BMCs between 0 and 100%.

Sclerostin antibodies

Romosozumab and blosozumab are mAbs that bind to sclerostin and prevent its inhibitory effects on bone formation (Recknor et al., 2015; Lim and Bolster, 2017; McClung et al., 2018). (Note that blosozumab was not approved for osteoporosis treatment at the time this manuscript was written.) Accordingly, we represent the mechanism of action of sclerostin antibodies by adding a new variable $s^{*}$ corresponding to the level of antibody-bound sclerostin and adding the dynamics of antibody-binding and unbinding in Equation 6,

\begin{array}{ll} \dot{s} & \to \dot{s} - κ_{s} Ψ_{s A b} s + δ_{s} s^{*}, \\ {\dot{s}}^{*} & = κ_{s} Ψ_{s A b} s - (δ_{s} + κ_{s}) s^{*}, \end{array}

where $Ψ_{sAb}$ is the effective sclerostin antibody concentration equivalent determined through Equation 21 and Equation 23. Here, $κ_{s}$ denotes the sclerostin degradation rate and $δ_{s}$ is the sclerostin/antibody binding rate; this parameterization implies that $Ψ_{sAb}$ is given in multiples of the effective antibody levels needed to achieve a binding rate equal to the unperturbed degradation grade of sclerostin.

Bisphosphonates

Bisphosphonates (like alendronate, ibandronate, risedronate, and zoledronate) bind to hydroxyapatite on the bone surface, thereby preventing osteoclasts from bone resorption; they further inhibit osteoclast-mediated bone resorption by promoting osteoclast apoptosis (Sato et al., 1991; Rodan and Fleisch, 1996). To simplify the model extension, we effectively represent the mechanism of action of bisphosphonates through the upregulation of the osteoclast apoptosis rate in Equation 13:

η_{C} \to η_{C} + g^{+} (Ψ_{b p}) η_{C}^{b p},

where $Ψ_{bp}$ is the effective bisphosphonate concentration equivalent determined through Equation 22 and Equation 23, and $η_{C}^{b p}$ is the maximum additional apoptosis rate caused by the presence of bisphosphonates.

PTH analogs

Teriparatide and abaloparatide are recombinant human PTH analogs (Jiang et al., 2003; Hattersley et al., 2016). PTH is known to have an anabolic effect on bone if administered intermittently while exerting a catabolic effect if administered continuously (Tam et al., 1982; Hock and Gera, 1992). Here, we consider an effective representation of the action of PTH in the anabolic regime only; this leads to a highly simplified and efficient description of the effective action of PTH therapies on bone turnover. However, it implies that the scope of our model is restricted to anabolic administration schemes and cannot be expected to yield correct results if probed in inappropriate regimes. In the anabolic regime, teriparatide downregulates osteoblast apoptosis (Jilka et al., 1999). Moreover, bone turnover markers show a marked increase early after treatment start but decline while drug administration remains unaltered (Leder et al., 2014, see Appendix 3—figure 1); such an effect is achieved in our model by rapid upregulation of osteoclast/osteoblast differentiation. We therefore also include a regulatory effect on osteoclast differentiation (which indirectly affects osteoblast differentiation as well):

\begin{array}{ll} η_{B} & \to [1 - β_{B}^{p t h} g^{+} (Ψ_{p t h})] η_{B}, \\ ω_{C^{*}} & \to [1 + β_{C^{*}}^{p t h} g^{+} (Ψ_{p t h})] ω_{C^{*}}, \end{array}

where $Ψ_{pth}$ is the effective teriparatide concentration equivalent determined through Equation 22 and Equation 23, and the parameters $β_{B}^{pth}$ and $β_{C^{*}}^{pth}$ parameterize the maximum effect strength of osteoblast apoptosis and pre-osteoclast differentiation, respectively.

Appendix 3

Simulations and parameter fits

Simulation protocol for aging and treatment

A model simulation was implemented in Python using standard NumPy and SciPy packages (Oliphant, 2006; Virtanen et al., 2020) and solved using SciPy’s ‘solve_ivp’ function with BDF solver. To compare how the model predicts the bone turnover dynamics of a hybrid aging/treatment dataset (see ‘Methods’ and Appendix 1—figure 1D) for a given set of model parameters, simulations were structured as follows. Drug dosing information of the corresponding dataset was provided to the model through the set of administered doses and the administration times: For the case of discrete dosing, dosing information consisted of doses c_i and administration times t_i entering Equation 21 (used for the drugs blosozumab, romosozumab, and denosumab). For the case of quasi-continuous dosing, dosing information consisted of doses per unit time c_i and time windows $[t_{i}, t_{i}^{*}]$ entering Equation 22 (used for the drugs alendronate and teriparatide).

The model was initialized at $t = 0$ in its steady state for all dynamic variables (i.e., the state for which all time derivatives are zero); except for the bone density $ρ_{b}$ , which does not possess a unique steady state and which was set to unity to represent peak bone density. We then simulated the model until well after the treatment period. All aging-related effects in the model were mediated by explicitly time-dependent auxiliary functions, as explained in Appendix 1 . To compare simulation results with clinical data, all relevant model variables were rescaled such that the first recorded data point in the corresponding clinical dataset coincided with the corresponding time point in the simulation, so that relative changes from a reference time point could be compared.

Parameter fits

To systematically fit model parameters, we defined a cost function that takes into account multiple fit quantities depending on their availability in the datasets. For a hybrid dataset $α$ and clinical quantity $β$ (BMD and serum levels of CTX, P1NP, and BSAP), we defined the distance function

D^{α β} = \frac{1}{\sum_{i} w_{i}^{α β}} \sum_{i} w_{i}^{α β} z_{i}^{α β} {(x_{i}^{α β} - {\hat{x}}_{i}^{α β})}^{2}, .

where x_i denotes the clinical data point at time point $i$ , ${\hat{x}}_{i}$ denotes the respective simulated data point, w_i denotes the relative weight of the respective data point depending on its certainty, and z_i denotes the relative weight depending on the time interval represented by the respective data point. For BTMs, we used the weights $w_{i}^{α β} = 1 / (1 + e_{i}^{α β})$ , where $e_{i}^{α β}$ is the mean of the upper and lower error bars of the respective quantity $β$ ; for the BMD, we used unit weights ( $w_{i}^{α, BMD} = 1$ ). To account for the fact that time intervals between data points may vastly differ (e.g., between the coarsely sampled aging dataset and the densely sampled treatment datasets), we included an interval-dependent weighting factor z_i, such that each data point was weighted by the average distance to its neighboring data points: The time interval-related weighting factors $z_{i}^{α β}$ were defined as $z_{i} = (δ_{i - 1} + δ_{i}) / 2$ ( $1 \leq i \leq n$ ), where $δ_{i} = t_{i + 1} - t_{i}$ ( $1 \leq i \leq n - 1$ ), $δ_{0} = δ_{1}$ , $δ_{n} = δ_{n - 1}$ and t_i denotes the time point of measurement $i$ .

We defined the combined cost function over all considered datasets as

J = \sum_{α β} W^{β} D^{α β}, .

where $W^{β}$ is an additional weighting factor that determines the relative importance of the different fit quantities in the cost function (Appendix 3—table 1).

In a first step, all fit parameters were manually adjusted for the model to exhibit a roughly sensible aging behavior. In a second step, a selected subset of hybrid aging/treatment datasets (see ‘Methods’ and Appendix 1—figure 1D) were used to fit all free model parameters. To perform the fits, we used a Covariance Matrix Adaptation Evolution Strategy via the Python package ‘pycma’ (Hansen et al., 2019). Results of the parameter fits are shown in Appendix 3—figure 1. The full list of fit parameters, including their final fit values, is given in Appendix 3—table 4.

Appendix 3—table 1

Values of fit weights $W^{β}$ used in Equation 29.

Weight	Value
$W^{(BMD)}$	300
$W^{(CTX)}$	1
$W^{(P1NP)}$	1
$W^{(BSAP)}$	1

BMD, bone mineral density; CTX, C-terminal telopeptide; P1NP, procollagen type 1 amino-terminal propeptide; BSAP, bone-specific alkaline phosphatase

Goodness-of-fit measures

To assess the goodness of the parameter fit and model predictions, we considered complementary goodness measures. The MAPE between clinical and simulated results is defined by

G_{MAPE} = \frac{1}{n} \sum_{i = 1}^{n} \frac{| x_{i} - \hat{x} (t_{i}) |}{x_{i}},

where the sum runs over all clinically recorded time points for the respective quantity (baseline changes of BMD and BTM levels), x_i denotes the clinical data point, t_i the time it was taken, and $\hat{x} (t)$ denotes the model result, which is a continuous function of time.

As a complementary measure, we introduce a ‘windowed minimal absolute percentage error’ (WMAPE), which indicates the mean minimal distance between model results and the data within a time window around the data point that reflects the average time spacing between data points. Formally, the WMAPE is given by

G_{WMAPE} = \frac{1}{n} \sum_{i = 1}^{n} \frac{1}{x_{i}} \min_{t \in [t_{i} - τ, t_{i} + τ]} | x_{i} - \hat{x} (t) |,

Here, we choose the time window $τ$ as half the median distance between data points, $τ = {median}_{i} (t_{i} - t_{i - 1}) / 2$ .

Appendix 3—table 2

Data sources used to calibrate and validate the model.

Columns titled ‘Figure(s)’ indicate the plot panels in the respective publication that were digitized. BMD always refers to total hip bone mineral density.

Publication	Medication(s)	Dosings	Figure(s)				Table(s)
Publication	Medication(s)	Dosings	BMD	CTX	P1NP	BSAP	BMD
Black et al., 2006	Alendronate	5–10 mg Q1D	2	3	3	—	—
Bone et al., 2011	Denosumab	60 mg Q6M	3b	4b	4a	—	—
Cosman et al., 2016	Romosozumab	210 mg Q1M	3b	3e	3d	—	—
	Denosumab	60 mg Q6M	3b	3e	3d	—	—
Leder et al., 2014	Teriparatide	20 mcg Q1D	2d	4c,f	4b,e	—	—
Leder et al., 2015	Teriparatide	20 mcg Q1D	3	4	—	—	—
	Denosumab	60 mg Q6M	3	4	—	—	—
Lewiecki et al., 2019	Romosozumab	210 mg Q1M	3b	—	—	—	—
	Denosumab	60 mg Q6M	3b	—	—	—	—
Looker et al., 1998	[Age-dependent BMD]	—	—	—	—	—	7
McClung et al., 2006	Denosumab	6 mg Q3M, 14 mg Q6M, 210 mg Q6M	2b	2e	—	2f	—
McClung et al., 2017	Denosumab	6–14 mg Q3M, 14–210 mg Q6M	2b	—	—	—	—
McClung et al., 2018	Romosozumab	140 mg Q1M, 210 mg Q1M	3c	4b	4a	—	—
	Denosumab	60 mg Q6M	3c,d	4b,d	4a,c	—	—
	Alendronate	70 mg Q1W	3d	4d	4c	—	—
Recknor et al., 2015	Blosozumab	180 mg Q4W, 180 mg Q2W, 270 mg Q2W	3b	4d	4a	—	—
Saag et al., 2017	Alendronate	70 mg Q1W	3b	3d	3c	—	—

Q, every; M, month; D, day; W, week; CTX, C-terminal telopeptide; P1NP, procollagen type 1 amino-terminal propeptide; BSAP, bone-specific alkaline phosphatase.

Appendix 3—table 3

Goodness-of-fit measures for calibration and validation datasets.

Mean absolute percentage error (MAPE) and windowed minimal absolute percentage error (WMAPE) as defined in Equation 30 and Equation 31, respectively. The column ‘Shown in’ indicates the figure in this article that shows the respective simulation and data plot.

Medication(s)		MAPE				WMAPE				Shown in
Medication(s)	Data ref.	BMD (%)	CTX (%)	P1NP (%)	BSAP (%)	BMD (%)	CTX (%)	P1NP (%)	BSAP (%)	Shown in
Calibration datasets
Alendronate 5–10 mg Q1D	Black et al., 2006	0.9	7.0	21.1	—	0.6	4.7	13.5	—	Appendix 3—figure 1
Alendronate 70 mg Q1W	Saag et al., 2017	0.5	34.4	13.7	—	0.2	20.3	10.2	—	Appendix 3—figure 1
Blosozumab 180 mg Q4W	Recknor et al., 2015	0.7	16.5	23.8	—	0.4	10.8	18.7	—	Appendix 3—figure 1
Blosozumab 270 mg Q2W	Recknor et al., 2015	0.3	26.8	13.2	—	0.2	15.8	4.8	—	Appendix 3—figure 1
Denosumab 14 mg Q3M → denosumab 60 mg Q6M	McClung et al., 2017	0.3	—	—	—	0.1	—	—	—	Appendix 3—figure 1
Denosumab 14 mg Q6M	McClung et al., 2006	0.2	73.9	—	17.0	0.0	46.3	—	0.2	Appendix 3—figure 1
Placebo	Recknor et al., 2015	0.4	7.3	15.4	—	0.3	7.2	15.3	—	Appendix 3—figure 1
Teriparatide 20 mcg Q1D	Leder et al., 2014	0.4	17.1	8.5	—	0.2	9.1	1.6	—	Appendix 3—figure 1
Teriparatide 20 mcg Q1D→ denosumab 60 mg Q6M	Leder et al., 2015	0.4	65.5	—	—	0.2	8.6	—	—	Appendix 3—figure 1
Validation datasets
Alendronate 70 mg Q1W → romosozumab 140 mg Q1M → denosumab 60 mg Q6M	McClung et al., 2018	0.5	25.4	20.0	—	0.3	17.1	3.6	—	Figure 2
Blosozumab 180 mg Q2W	Recknor et al., 2015	0.3	21.3	20.3	—	0.1	13.2	12.6	—	Figure 2
Denosumab 60 mg Q6M → teriparatide 20 mcg Q1D	Leder et al., 2015	0.4	103.1	—	—	0.3	44.6	—	—	Figure 2
Placebo	McClung et al., 2018	0.6	6.2	10.5	—	0.5	6.2	10.5	—	Figure 2
Placebo→ denosumab 60 mg Q6M	McClung et al., 2018	0.6	13.2	10.1	—	0.4	5.1	7.4	—	Figure 2
Teriparatide 20 mcg Q1D + denosumab 60 mg Q6M → denosumab 60 mg Q6M	Leder et al., 2015	0.7	183.7	—	—	0.4	66.4	—	—	Figure 2
Alendronate 70 mg Q1W → romosozumab 140 mg Q1M → placebo	McClung et al., 2018	0.5	22.9	16.8	—	0.4	15.1	4.2	—	Figure 2—figure supplement 1
Placebo→ denosumab 60 mg Q6M	Cosman et al., 2016	0.5	81.9	9.6	—	0.0	44.4	4.0	—	Figure 2—figure supplement 1
Placebo→ denosumab 60 mg Q6M	Lewiecki et al., 2019	0.4	—	—	—	0.1	—	—	—	Figure 2—figure supplement 1
Placebo→ denosumab 60 mg Q6M	McClung et al., 2017	0.6	—	—	—	0.3	—	—	—	Figure 2—figure supplement 1
Romosozumab 210 mg Q1M → alendronate 70 mg Q1W	Saag et al., 2017	1.0	30.0	37.1	—	0.5	19.4	11.7	—	Figure 2—figure supplement 1
Romosozumab 210 mg Q1M → denosumab 60 mg Q6M	Cosman et al., 2016	1.2	92.5	22.4	—	0.6	50.7	6.9	—	Figure 2—figure supplement 1
Romosozumab 210 mg Q1M → denosumab 60 mg Q6M	Lewiecki et al., 2019	1.0	—	—	—	0.5	—	—	—	Figure 2—figure supplement 1
Romosozumab 210 mg Q1M → denosumab 60 mg Q6M	McClung et al., 2018	0.6	14.2	57.8	—	0.3	6.5	15.9	—	Figure 2—figure supplement 1
Romosozumab 210 mg Q1M → placebo	McClung et al., 2018	0.7	15.0	53.3	—	0.5	10.2	18.5	—	Figure 2—figure supplement 1

Q, every; M, month; D, day; W, week; BMD, bone mineral density; CTX, C-terminal telopeptide; P1NP, procollagen type 1 amino-terminal propeptide; BSAP, bone-specific alkaline phosphatase.

Appendix 3—table 4

Full list of parameters of the core model and the medication extensions.

Parameter	Description	Value	Unit	Origin	Model equation
Core model
$ω_{C^{*}}$	Reference pre-osteoclast to osteoclast differentiation rate	0.93	$d^{- 1}$	Calibration	Equation 10
$e_{C^{*}}$	Estrogen threshold for downregulation of pre-osteoclast to osteoclast differentiation	0.94	1	Calibration	Equation 10
$s_{C^{*}}$	Sclerostin threshold for upregulation of pre-osteoclast to osteoclast differentiation	$8.60 \times 10^{6}$	1	Calibration	Equation 10
$η_{C}$	Reference osteoclast apoptosis rate	0.02	$d^{- 1}$	Calibration	Equation 13
$e_{C}$	Estrogen threshold for upregulation of osteoclast apoptosis	0.99	1	Calibration	Equation 13
$r_{C}$	Resorption signal threshold for upregulation of osteoclast apoptosis	10.10	1	Calibration	Equation 13
$ν_{C}$	Max. rel. effect of regulatory factors on osteoclast apoptosis	$1.23 \times 10^{- 4}$	1	Calibration	Equation 13
$ω_{B^{*}}$	Reference pre-osteoblast to osteoblast differentiation rate	0.32	$d^{- 1}$	Calibration	Equation 10
$s_{B^{*}}$	Sclerostin threshold for downregulation of pre-osteoblast to osteoblast differentiation	$1.63 \times 10^{2}$	1	Calibration	Equation 10
$η_{B}$	Reference osteoblast apoptosis rate	$8.68 \times 10^{- 3}$	$d^{- 1}$	Calibration	Equation 13
$ω_{B}$	Reference osteoblast to osteocyte conversion rate	$6.24 \times 10^{- 4}$	$d^{- 1}$	Calibration	Equation 11
$η_{Y}$	osteocyte apoptosis rate	$1.10 \times 10^{- 4}$	$d^{- 1}$	Estimate	Equation 14
$κ_{s}$	Sclerostin degradation rate	0.05	$d^{- 1}$	Estimate; see Suen et al., 2015; Ominsky et al., 2015.	Equation 15
$e_{s}$	Estrogen threshold for downregulation of sclerostin secretion	9.60	1	Calibration	Equation 15
$λ_{C}$	Reference bone resorption rate per unit density osteoclast	$3.82 \times 10^{- 6}$	$d^{- 1}$	Calibration	Equation 16
$λ_{B}$	Reference bone formation rate per unit density osteoblast	$1.29 \times 10^{- 6}$	$d^{- 1}$	Calibration	Equation 16
$s_{Ω}$	Sclerostin threshold for downregulation of bone formation	$3.04 \times 10^{3}$	1	Calibration	Equation 16
$r_{Ω}$	Resorption signal threshold for upregulation of bone formation	$1.02 \times 10^{3}$	1	Calibration	Equation 16
$ν_{Ω}$	Max. rel. effect of the resorption signal on bone formation	$1.08 \times 10^{2}$	1	Calibration	Equation 16
$γ$	Equilibration rate of the bone mineral content	$6.65 \times 10^{- 3}$	$d^{- 1}$	Calibration	Equation 17
$c_{0}$	Reference bone mineral content	0.80	1	Estimate	Equation 17
$t_{e}$	Onset of estrogen decline	50.00	$y$	Estimate	Equation 18
$τ_{e}$	Time scale of estrogen decline	2.60	$y$	Indep. fit (Appendix 1—figure 1A)	Equation 18
Bone turnover markers
$q_{CTX}$	Exponent relating the bone resorption rate to CTX levels	1.16	1	Calibration	Equation 15
$q_{P1NP}$	Exponent relating the bone formation rate to P1NP levels	1.45	1	Calibration	Equation 15
$q_{BSAP}$	Exponent relating the bone formation rate to BSAP levels	0.92	1	Calibration	Equation 15
Medication extension: sclerostin antibodies
$E_{blosozumab}$	Efficacy: blosozumab	0.01	1	Calibration	Equation 21
$T_{blosozumab}$	Effective half-life: blosozumab	7.00	$d$	$T_{romosozumab}$	Equation 21
$E_{romosozumab}$	Efficacy: romosozumab	0.01	1	$E_{blosozumab}$	Equation 21
$T_{romosozumab}$	Effective half-life: romosozumab	7.00	$d$	Solling et al., 2018	Equation 21
$δ_{s}$	Sclerostin/antibody unbinding rate	0.05	$d^{- 1}$	$κ_{s}$	Equation 25
Medication extension: RANKL antibodies
$E_{denosumab}$	Efficacy: denosumab	$4.34 \times 10^{3}$	1	Calibration	Equation 21
$T_{denosumab}$	Effective half-life: denosumab	10.00	$d$	Bekker et al., 2004	Equation 21
$β_{C^{*}}^{rAb}$	Max. rel. effect of RANKL antibodies on pre-osteoclast to osteoclast differentiation	0.87	1	Calibration	Equation 24
$β_{b}^{rAb}$	Max. rel. effect of RANKL antibodies on mineralization	0.02	1	Calibration	Equation 24
Medication extension: bisphosphonates
$E_{alendronate}$	Efficacy: alendronate	$2.97 \times 10^{- 5}$	1	Calibration	Equation 22
$T_{alendronate}$	Effective half-life: alendronate	$1.53 \times 10^{2}$	$d$	Calibration	Equation 22
$η_{C}^{bp}$	Max. contribution of bisphosphonates to osteoclast apoptosis rate	1.00	$d^{- 1}$	Calibration	Equation 26
Medication extension: PTH analogs
$E_{teriparatide}$	Efficacy: teriparatide	0.27	1	Calibration	Equation 22
$T_{teriparatide}$	Effective half-life: teriparatide	0.04	$d$	Satterwhite et al., 2010	Equation 22
$β_{B}^{pth}$	Max. rel. effect of PTH analogs on osteoblast apoptosis	1.31	1	Calibration	Equation 27
$β_{C^{*}}^{pth}$	Max. rel. effect of PTH analogs on pre-osteoclast to osteoclast differentiation	4.28	1	Calibration	Equation 27

CTX, C-terminal telopeptide; P1NP, procollagen type 1 amino-terminal propeptide; BSAP, bone-specific alkaline phosphatase; PTH, parathyroid hormone.

Appendix 3—figure 1

Download asset Open asset

Calibration datasets comparing model predictions and clinical data from various studies.

All conventions identical to Figure 2. Drug administrations are provided in the bottom row. See Appendix 3—table 2 for a list of data sources and Appendix 3—table 3 for goodness-of-fit measures. Dosing: mg, milligrams; mcg, micrograms; Q $x$ M, dose administered every $x$ months; Q $x$ W, every $x$ weeks; Q $x$ D, every $x$ days; B, blosozumab; A, alendronate; D, denosumab; T, teriparatide.

Data availability

The current manuscript is a computational study, so no data have been generated for this manuscript. Modelling code including the code needed to produce simulation-related figures is part of the Source Code Files. Digitised data from previously published scientific articles are also part of the Source Code Files. All original source publications are specified in Appendix 3 Table 2 of the manuscript. Data were preprocessed as described in the 'Methods' section of the manuscript. Model parameters are also part of the Source Code Files. In addition, they are provided in Appendix 3 Table 4 of the manuscript.

References

1. Atkins GJ
2. Rowe PS
3. Lim HP
4. Welldon KJ
5. Ormsby R
6. Wijenayaka AR
7. Zelenchuk L
8. Evdokiou A
9. Findlay DM
(2011) Sclerostin is a locally acting regulator of late-osteoblast/preosteocyte differentiation and regulates mineralization through a MEPE-ASARM-dependent mechanism
Journal of Bone and Mineral Research 26:1425–1436.

https://doi.org/10.1002/jbmr.345
- PubMed
- Google Scholar
1. Austin M
2. Yang YC
3. Vittinghoff E
4. Adami S
5. Boonen S
6. Bauer DC
7. Bianchi G
8. Bolognese MA
9. Christiansen C
10. Eastell R
11. Grauer A
12. Hawkins F
13. Kendler DL
14. Oliveri B
15. McClung MR
16. Reid IR
17. Siris ES
18. Zanchetta J
19. Zerbini CAF
20. Libanati C
21. Cummings SR
22. FREEDOM Trial
(2012) Relationship between bone mineral density changes with denosumab treatment and risk reduction for vertebral and nonvertebral fractures
Journal of Bone and Mineral Research 27:687–693.

https://doi.org/10.1002/jbmr.1472
- PubMed
- Google Scholar
1. Ayati BP
2. Edwards CM
3. Webb GF
4. Wikswo JP
(2010) A mathematical model of bone remodeling dynamics for normal bone cell populations and myeloma bone disease
Biology Direct 5:1–17.

https://doi.org/10.1186/1745-6150-5-28
- PubMed
- Google Scholar
1. Bekker PJ
2. Holloway DL
3. Rasmussen AS
4. Murphy R
5. Martin SW
6. Leese PT
7. Holmes GB
8. Dunstan CR
9. DePaoli AM
(2004) A single-dose placebo-controlled study of AMG 162, A fully human monoclonal antibody to RANKL, in postmenopausal women
Journal of Bone and Mineral Research 19:1059–1066.

https://doi.org/10.1359/JBMR.040305
- PubMed
- Google Scholar
1. Berger C
2. Langsetmo L
3. Joseph L
4. Hanley DA
5. Davison KS
6. Josse RG
7. Prior JC
8. Kreiger N
9. Tenenhouse A
10. Goltzman D
11. CaMos Research Group
(2009) Association between change in BMD and fragility fracture in women and men
Journal of Bone and Mineral Research 24:361–370.

https://doi.org/10.1359/jbmr.081004
- PubMed
- Google Scholar
(2015) Application of a systems pharmacology-based placebo population model to analyze long-term data of postmenopausal osteoporosis
CPT 4:516–526.

https://doi.org/10.1002/psp4.12006
- PubMed
- Google Scholar
1. Black DM
2. Schwartz AV
3. Ensrud KE
4. Cauley JA
5. Levis S
6. Quandt SA
7. Satterfield S
8. Wallace RB
9. Bauer DC
10. Palermo L
11. Wehren LE
12. Lombardi A
13. Santora AC
14. Cummings SR
(2006) Effects of continuing or stopping alendronate after 5 years of treatment
JAMA 296:2927.

https://doi.org/10.1001/jama.296.24.2927
- PubMed
- Google Scholar
1. Black DM
2. Cauley JA
3. Wagman R
4. Ensrud K
5. Fink HA
6. Hillier TA
7. Lui LY
8. Cummings SR
9. Schousboe JT
10. Napoli N
(2018) The ability of a single BMD and fracture history assessment to predict fracture over 25 years in postmenopausal women: The study of osteoporotic fractures
Journal of Bone and Mineral Research 33:389–395.

https://doi.org/10.1002/jbmr.3194
- PubMed
- Google Scholar
1. Bone HG
2. Bolognese MA
3. Yuen CK
4. Kendler DL
5. Miller PD
6. Yang YC
7. Grazette L
8. San Martin J
9. Gallagher JC
(2011) Effects of denosumab treatment and discontinuation on bone mineral density and bone turnover markers in postmenopausal women with low bone mass
The Journal of Clinical Endocrinology and Metabolism 96:972–980.

https://doi.org/10.1210/jc.2010-1502
- PubMed
- Google Scholar
1. Boskey AL
(2013) Bone composition: relationship to bone fragility and antiosteoporotic drug effects
BoneKEy Reports 2:447.

https://doi.org/10.1038/bonekey.2013.181
- PubMed
- Google Scholar
1. Boyce BF
2. Xing L
(2008) Functions of RANKL/RANK/OPG in bone modeling and remodeling
Archives of Biochemistry and Biophysics 473:139–146.

https://doi.org/10.1016/j.abb.2008.03.018
- PubMed
- Google Scholar
(2011) Spatio-temporal structure of cell distribution in cortical bone multicellular units: A mathematical model
Bone 48:918–926.

https://doi.org/10.1016/j.bone.2010.12.009
- PubMed
- Google Scholar
(2014) Bone refilling in cortical basic multicellular units: insights into tetracycline double labelling from a computational model
Biomechanics and Modeling in Mechanobiology 13:185–203.

https://doi.org/10.1007/s10237-013-0495-y
- PubMed
- Google Scholar
1. Burge R
2. Dawson-Hughes B
3. Solomon DH
4. Wong JB
5. King A
6. Tosteson A
(2007) Incidence and economic burden of osteoporosis-related fractures in the United States, 2005-2025
Journal of Bone and Mineral Research 22:465–475.

https://doi.org/10.1359/jbmr.061113
- PubMed
- Google Scholar
1. Clarke B
(2008) Normal bone anatomy and physiology
Clinical Journal of the American Society of Nephrology 3:S131–S139.

https://doi.org/10.2215/CJN.04151206
- PubMed
- Google Scholar
Report
1. Cooper C
2. Ferrari S
(2019) IOF Compendium of Osteoporosis
International Osteoporosis Foundation.

https://www.osteoporosis.foundation/sites/iofbonehealth/files/2020-01/IOF-Compendium-of-Osteoporosis-web-V02.pdf
- Google Scholar
1. Cosman F
2. Crittenden DB
3. Adachi JD
4. Binkley N
5. Czerwinski E
6. Ferrari S
7. Hofbauer LC
8. Lau E
9. Lewiecki EM
10. Miyauchi A
11. Zerbini CAF
12. Milmont CE
13. Chen L
14. Maddox J
15. Meisner PD
16. Libanati C
17. Grauer A
(2016) Romosozumab treatment in postmenopausal women with osteoporosis
The New England Journal of Medicine 375:1532–1543.

https://doi.org/10.1056/NEJMoa1607948
- PubMed
- Google Scholar
(2017) Role and mechanism of action of sclerostin in bone
Bone 96:29–37.

https://doi.org/10.1016/j.bone.2016.10.007
- PubMed
- Google Scholar
1. Ensrud KE
2. Lui LY
3. Crandall CJ
4. Orwoll ES
5. Langsetmo L
6. Schousboe JT
7. Fink HA
8. Lane NE
9. Kado DM
10. Cauley JA
11. Stefanick ML
12. Cawthon PM
(2022) Repeat bone mineral density screening measurement and fracture prediction in older men: a prospective cohort study
The Journal of Clinical Endocrinology and Metabolism 19:dgac324.

https://doi.org/10.1210/clinem/dgac324
- PubMed
- Google Scholar
1. Eriksen EF
(2010) Cellular mechanisms of bone remodeling
Reviews in Endocrine & Metabolic Disorders 11:219–227.

https://doi.org/10.1007/s11154-010-9153-1
- PubMed
- Google Scholar
(2015) Connecting the dots: linking osteocyte activity and therapeutic modulation of sclerostin by extending a multiscale systems model
CPT 4:527–536.

https://doi.org/10.1002/psp4.12013
- PubMed
- Google Scholar
1. Fuller K
2. Chambers TJ
(1995) Localisation of mRNA for collagenase in osteocytic, bone surface and chondrocytic cells but not osteoclasts
Journal of Cell Science 108:2221–2230.

https://doi.org/10.1242/jcs.108.6.2221
- PubMed
- Google Scholar
(1999) Skeletal bone morphogenetic proteins suppress the expression of collagenase-3 by rat osteoblasts
Endocrinology 140:562–567.

https://doi.org/10.1210/endo.140.2.6493
- PubMed
- Google Scholar
(2013) The role of osteocytes in targeted bone remodeling: a mathematical model
PLOS ONE 8:e63884.

https://doi.org/10.1371/journal.pone.0063884
- PubMed
- Google Scholar
(1997) World-wide projections for hip fracture
Osteoporosis International 7:407–413.

https://doi.org/10.1007/pl00004148
- PubMed
- Google Scholar
Software
(2019) CMA-ES/Pycma on Github, version r3.2.2
Zenodo.

https://doi.org/10.5281/zenodo.2559634
1. Hattersley G
2. Dean T
3. Corbin BA
4. Bahar H
5. Gardella TJ
(2016) Binding selectivity of abaloparatide for PTH-Type-1-receptor conformations and effects on downstream signaling
Endocrinology 157:141–149.

https://doi.org/10.1210/en.2015-1726
- PubMed
- Google Scholar
1. Hernlund E
2. Svedbom A
3. Ivergård M
4. Compston J
5. Cooper C
6. Stenmark J
7. McCloskey EV
8. Jönsson B
9. Kanis JA
(2013) Osteoporosis in the European union: medical management, epidemiology and economic burden
Archives of Osteoporosis 8:1–115.

https://doi.org/10.1007/s11657-013-0136-1
- Google Scholar
1. Hock JM
2. Gera I
(1992) Effects of continuous and intermittent administration and inhibition of resorption on the anabolic response of bone to parathyroid hormone
Journal of Bone and Mineral Research 7:65–72.

https://doi.org/10.1002/jbmr.5650070110
- PubMed
- Google Scholar
1. Houde N
2. Chamoux E
3. Bisson M
4. Roux S
(2009) Transforming growth factor-beta1 (TGF-beta1) induces human osteoclast apoptosis by up-regulating Bim
The Journal of Biological Chemistry 284:23397–23404.

https://doi.org/10.1074/jbc.M109.019372
- PubMed
- Google Scholar
1. Hughes DE
2. Dai A
3. Tiffee JC
4. Li HH
5. Mundy GR
6. Boyce BF
(1996) Estrogen promotes apoptosis of murine osteoclasts mediated by TGF-beta
Nature Medicine 2:1132–1136.

https://doi.org/10.1038/nm1096-1132
- PubMed
- Google Scholar
1. Jiang Y
2. Zhao JJ
3. Mitlak BH
4. Wang O
5. Genant HK
6. Eriksen EF
(2003) Recombinant human parathyroid hormone (1-34) [teriparatide] improves both cortical and cancellous bone structure
Journal of Bone and Mineral Research 18:1932–1941.

https://doi.org/10.1359/jbmr.2003.18.11.1932
- PubMed
- Google Scholar
(1999) Increased bone formation by prevention of osteoblast apoptosis with parathyroid hormone
The Journal of Clinical Investigation 104:439–446.

https://doi.org/10.1172/JCI6610
- PubMed
- Google Scholar
1. Kameda T
2. Mano H
3. Yuasa T
4. Mori Y
5. Miyazawa K
6. Shiokawa M
7. Nakamaru Y
8. Hiroi E
9. Hiura K
10. Kameda A
11. Yang NN
12. Hakeda Y
13. Kumegawa M
(1997) Estrogen inhibits bone resorption by directly inducing apoptosis of the bone-resorbing osteoclasts
The Journal of Experimental Medicine 186:489–495.

https://doi.org/10.1084/jem.186.4.489
- PubMed
- Google Scholar
1. Kanis JA
2. Johnell O
3. Oden A
4. Dawson A
5. De Laet C
6. Jonsson B
(2001) Ten year probabilities of osteoporotic fractures according to BMD and diagnostic thresholds
Osteoporosis International 12:989–995.

https://doi.org/10.1007/s001980170006
- PubMed
- Google Scholar
1. Kanis JA
2. Oden A
3. Johansson H
4. Borgström F
5. Ström O
6. McCloskey E
(2009) FRAX and its applications to clinical practice
Bone 44:734–743.

https://doi.org/10.1016/j.bone.2009.01.373
- PubMed
- Google Scholar
Book
1. Keener J
2. Sneyd J
(2009) Mathematical Physiology (2nd edition)
New York, NY: Springer-Verlag.

https://doi.org/10.1007/978-0-387-75847-3
- Google Scholar
1. Komarova SV
2. Smith RJ
3. Dixon SJ
4. Sims SM
5. Wahl LM
(2003) Mathematical model predicts a critical role for osteoclast autocrine regulation in the control of bone remodeling
Bone 33:206–215.

https://doi.org/10.1016/s8756-3282(03)00157-1
- PubMed
- Google Scholar
1. Komarova SV
2. Safranek L
3. Gopalakrishnan J
4. Ou MJY
5. McKee MD
6. Murshed M
7. Rauch F
8. Zuhr E
(2015) Mathematical model for bone mineralization
Frontiers in Cell and Developmental Biology 3:51.

https://doi.org/10.3389/fcell.2015.00051
- PubMed
- Google Scholar
1. Kostenuik PJ
2. Nguyen HQ
3. McCabe J
4. Warmington KS
5. Kurahara C
6. Sun N
7. Chen C
8. Li L
9. Cattley RC
10. Van G
11. Scully S
12. Elliott R
13. Grisanti M
14. Morony S
15. Tan HL
16. Asuncion F
17. Li X
18. Ominsky MS
19. Stolina M
20. Dwyer D
21. Dougall WC
22. Hawkins N
23. Boyle WJ
24. Simonet WS
25. Sullivan JK
(2009) Denosumab, a fully human monoclonal antibody to RANKL, inhibits bone resorption and increases BMD in knock-in mice that express chimeric (murine/human) RANKL
Journal of Bone and Mineral Research 24:182–195.

https://doi.org/10.1359/jbmr.081112
- PubMed
- Google Scholar
1. Krege JH
2. Wan X
3. Lentle BC
4. Berger C
5. Langsetmo L
6. Adachi JD
7. Prior JC
8. Tenenhouse A
9. Brown JP
10. Kreiger N
11. Olszynski WP
12. Josse RG
13. Goltzman D
14. CaMos Research Group
(2013) Fracture risk prediction: importance of age, BMD and spine fracture status
BoneKEy Reports 2:404.

https://doi.org/10.1038/bonekey.2013.138
- PubMed
- Google Scholar
(2020) Study of the combined effects of PTH treatment and mechanical loading in postmenopausal osteoporosis using a new mechanistic PK-PD model
Biomechanics and Modeling in Mechanobiology 19:1765–1780.

https://doi.org/10.1007/s10237-020-01307-6
- PubMed
- Google Scholar
1. Leder BZ
2. Tsai JN
3. Uihlein AV
4. Burnett-Bowie SAM
5. Zhu Y
6. Foley K
7. Lee H
8. Neer RM
(2014) Two years of Denosumab and teriparatide administration in postmenopausal women with osteoporosis (The DATA Extension Study): A randomized controlled trial
The Journal of Clinical Endocrinology and Metabolism 99:1694–1700.

https://doi.org/10.1210/jc.2013-4440
- PubMed
- Google Scholar
1. Leder BZ
2. Tsai JN
3. Uihlein AV
4. Wallace PM
5. Lee H
6. Neer RM
7. Burnett-Bowie SAM
(2015) Denosumab and teriparatide transitions in postmenopausal osteoporosis (the DATA-Switch study): extension of a randomised controlled trial
Lancet 386:1147–1155.

https://doi.org/10.1016/S0140-6736(15)61120-5
- PubMed
- Google Scholar
1. Lemaire V
2. Tobin FL
3. Greller LD
4. Cho CR
5. Suva LJ
(2004) Modeling the interactions between osteoblast and osteoclast activities in bone remodeling
Journal of Theoretical Biology 229:293–309.

https://doi.org/10.1016/j.jtbi.2004.03.023
- PubMed
- Google Scholar
1. Lemaire V
2. Cox DR
(2019) Dynamics of bone cell interactions and differential responses to PTH and antibody-based therapies
Bulletin of Mathematical Biology 81:3575–3622.

https://doi.org/10.1007/s11538-018-0533-0
- PubMed
- Google Scholar
(2019) One year of romosozumab followed by two years of denosumab maintains fracture risk reductions: results of the FRAME extension study
Journal of Bone and Mineral Research 34:419–428.

https://doi.org/10.1002/jbmr.3622
- PubMed
- Google Scholar
1. Li X
2. Zhang Y
3. Kang H
4. Liu W
5. Liu P
6. Zhang J
7. Harris SE
8. Wu D
(2005) Sclerostin binds to LRP5/6 and antagonizes canonical Wnt signaling
The Journal of Biological Chemistry 280:19883–19887.

https://doi.org/10.1074/jbc.M413274200
- PubMed
- Google Scholar
1. Li X
2. Ominsky MS
3. Niu QT
4. Sun N
5. Daugherty B
6. D’Agostin D
7. Kurahara C
8. Gao Y
9. Cao J
10. Gong J
11. Asuncion F
12. Barrero M
13. Warmington K
14. Dwyer D
15. Stolina M
16. Morony S
17. Sarosi I
18. Kostenuik PJ
19. Lacey DL
20. Simonet WS
21. Ke HZ
22. Paszty C
(2008) Targeted deletion of the sclerostin gene in mice results in increased bone formation and bone strength
Journal of Bone and Mineral Research 23:860–869.

https://doi.org/10.1359/jbmr.080216
- PubMed
- Google Scholar
1. Lim SY
2. Bolster MB
(2017) Profile of romosozumab and its potential in the management of osteoporosis
Drug Design, Development and Therapy 11:1221–1231.

https://doi.org/10.2147/DDDT.S127568
- PubMed
- Google Scholar
1. Lisberg A
2. Ellis R
3. Nicholson K
4. Moku P
5. Swarup A
6. Dhurjati P
7. Nohe A
(2017) Mathematical modeling of the effects of CK2.3 on mineralization in osteoporotic bone
CPT 6:208–215.

https://doi.org/10.1002/psp4.12154
- PubMed
- Google Scholar
1. Looker AC
2. Wahner HW
3. Dunn WL
4. Calvo MS
5. Harris TB
6. Heyse SP
7. Johnston CC
8. Lindsay R
(1998) Updated data on proximal femur bone mineral levels of US adults
Osteoporosis International 8:468–489.

https://doi.org/10.1007/s001980050093
- PubMed
- Google Scholar
1. Manolagas SC
(2000) Birth and death of bone cells: basic regulatory mechanisms and implications for the pathogenesis and treatment of osteoporosis
Endocrine Reviews 21:115–137.

https://doi.org/10.1210/edrv.21.2.0395
- PubMed
- Google Scholar
1. Manolagas SC
(2010) From estrogen-centric to aging and oxidative stress: a revised perspective of the pathogenesis of osteoporosis
Endocrine Reviews 31:266–300.

https://doi.org/10.1210/er.2009-0024
- PubMed
- Google Scholar
(2008) Integrated cellular bone homeostasis model for denosumab pharmacodynamics in multiple myeloma patients
The Journal of Pharmacology and Experimental Therapeutics 326:555–562.

https://doi.org/10.1124/jpet.108.137703
- PubMed
- Google Scholar
(2011) Integrated model for denosumab and ibandronate pharmacodynamics in postmenopausal women
Biopharmaceutics & Drug Disposition 32:471–481.

https://doi.org/10.1002/bdd.770
- PubMed
- Google Scholar
(2020) The role of sclerostin in bone and ectopic calcification
International Journal of Molecular Sciences 21:3199.

https://doi.org/10.3390/ijms21093199
- PubMed
- Google Scholar
1. Martínez-Reina J
2. Pivonka P
(2019) Effects of long-term treatment of denosumab on bone mineral density: insights from an in-silico model of bone mineralization
Bone 125:87–95.

https://doi.org/10.1016/j.bone.2019.04.022
- PubMed
- Google Scholar
1. McClung MR
2. Lewiecki EM
3. Cohen SB
4. Bolognese MA
5. Woodson GC
6. Moffett AH
7. Peacock M
8. Miller PD
9. Lederman SN
10. Chesnut CH
11. Lain D
12. Kivitz AJ
13. Holloway DL
14. Zhang C
15. Peterson MC
16. Bekker PJ
(2006) Denosumab in postmenopausal women with low Bone mineral density
New England Journal of Medicine 354:821–831.

https://doi.org/10.1056/NEJMoa044459
- PubMed
- Google Scholar
1. McClung MR
2. Wagman RB
3. Miller PD
4. Wang A
5. Lewiecki EM
(2017) Observations following discontinuation of long-term denosumab therapy
Osteoporosis International 28:1723–1732.

https://doi.org/10.1007/s00198-017-3919-1
- PubMed
- Google Scholar
1. McClung MR
2. Brown JP
3. Diez-Perez A
4. Resch H
5. Caminis J
6. Meisner P
7. Bolognese MA
8. Goemaere S
9. Bone HG
10. Zanchetta JR
11. Maddox J
12. Bray S
13. Grauer A
(2018) Effects of 24 months of treatment with romosozumab followed by 12 months of denosumab or placebo in postmenopausal women With low bone mineral density: a randomized, double-blind, phase 2, parallel group study
Journal of Bone and Mineral Research 33:1397–1406.

https://doi.org/10.1002/jbmr.3452
- PubMed
- Google Scholar
1. Mödder UIL
2. Clowes JA
3. Hoey K
4. Peterson JM
5. McCready L
6. Oursler MJ
7. Riggs BL
8. Khosla S
(2011) Regulation of circulating sclerostin levels by sex steroids in women and in men
Journal of Bone and Mineral Research 26:27–34.

https://doi.org/10.1002/jbmr.128
- PubMed
- Google Scholar
1. Mundy GR
2. Guise TA
(1999) Hormonal control of calcium homeostasis
Clinical Chemistry 45:1347–1352.

https://doi.org/10.1093/clinchem/45.8.1347
- PubMed
- Google Scholar
(2013) Assessing the impact of osteoporosis on the burden of hip fractures
Calcified Tissue International 92:42–49.

https://doi.org/10.1007/s00223-012-9666-6
- PubMed
- Google Scholar
Book
1. Oliphant TE
(2006)
A Guide to NumPy

Trelgol Publishing USA.
- Google Scholar
1. Ominsky MS
2. Vlasseros F
3. Jolette J
4. Smith SY
5. Stouch B
6. Doellgast G
7. Gong J
8. Gao Y
9. Cao J
10. Graham K
11. Tipton B
12. Cai J
13. Deshpande R
14. Zhou L
15. Hale MD
16. Lightwood DJ
17. Henry AJ
18. Popplewell AG
19. Moore AR
20. Robinson MK
21. Lacey DL
22. Simonet WS
23. Paszty C
(2015) Two doses of sclerostin antibody in cynomolgus monkeys increases bone formation
Bone Mineral Density, and Bone Strength. J. Bone Miner. Res 25:948–959.

https://doi.org/10.1002/jbmr.14
- PubMed
- Google Scholar
(2006) Secondary osteoporosis: a review of the recent evidence
Endocrine Practice 12:436–445.

https://doi.org/10.4158/EP.12.4.436
- PubMed
- Google Scholar
1. Peterson MC
2. Riggs MM
(2010) A physiologically based mathematical model of integrated calcium homeostasis and bone remodeling
Bone 46:49–63.

https://doi.org/10.1016/j.bone.2009.08.053
- PubMed
- Google Scholar
1. Pivonka P
2. Zimak J
3. Smith DW
4. Gardiner BS
5. Dunstan CR
6. Sims NA
7. John Martin T
8. Mundy GR
(2008) Model structure and control of bone remodeling: a theoretical study
Bone 43:249–263.

https://doi.org/10.1016/j.bone.2008.03.025
- PubMed
- Google Scholar
1. Pivonka P
2. Zimak J
3. Smith DW
4. Gardiner BS
5. Dunstan CR
6. Sims NA
7. Martin TJ
8. Mundy GR
(2010) Theoretical investigation of the role of the RANK-RANKL-OPG system in bone remodeling
Journal of Theoretical Biology 262:306–316.

https://doi.org/10.1016/j.jtbi.2009.09.021
- PubMed
- Google Scholar
(2013) The influence of bone surface availability in bone remodelling—A mathematical model including coupled geometrical and biomechanical regulations of bone cells
Engineering Structures 47:134–147.

https://doi.org/10.1016/j.engstruct.2012.09.006
- Google Scholar
1. Recknor CP
2. Recker RR
3. Benson CT
4. Robins DA
5. Chiang AY
6. Alam J
7. Hu L
8. Matsumoto T
9. Sowa H
10. Sloan JH
11. Konrad RJ
12. Mitlak BH
13. Sipos AA
(2015) The effect of discontinuing treatment with blosozumab: follow-up results of a phase 2 randomized clinical trial in postmenopausal women with low bone mineral density
Journal of Bone and Mineral Research 30:1717–1725.

https://doi.org/10.1002/jbmr.2489
- PubMed
- Google Scholar
(1998) 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 13:763–773.

https://doi.org/10.1359/jbmr.1998.13.5.763
- PubMed
- Google Scholar
1. Rodan GA
2. Fleisch HA
(1996) Bisphosphonates: mechanisms of action
The Journal of Clinical Investigation 97:2692–2696.

https://doi.org/10.1172/JCI118722
- PubMed
- Google Scholar
1. Ross DS
2. Battista C
3. Cabal A
4. Mehta K
(2012) Dynamics of bone cell signaling and PTH treatments of osteoporosis
Discrete and Continuous Dynamical Systems - Series B 17:2185–2200.

https://doi.org/10.3934/dcdsb.2012.17.2185
- Google Scholar
(1997a) Insulin-like growth factor I inhibits the transcription of collagenase 3 in osteoblast cultures
Journal of Cellular Biochemistry 67:176–183.

https://doi.org/10.1002/(SICI)1097-4644(19971101)67:2<176::AID-JCB3>3.0.CO;2-U
- Google Scholar
(1997b) Transforming growth factor β1 inhibits collagenase 3 expression by transcriptional and post-transcriptional mechanisms in osteoblast cultures
Journal of Cellular Physiology 170:145–152.

https://doi.org/10.1002/(SICI)1097-4652(199702)170:2<145::AID-JCP6>3.0.CO;2-O
- Google Scholar
(2009) Mathematical modeling of spatio-temporal dynamics of a single bone multicellular unit
Journal of Bone and Mineral Research 24:860–870.

https://doi.org/10.1359/jbmr.081229
- PubMed
- Google Scholar
1. Saag KG
2. Petersen J
3. Brandi ML
4. Karaplis AC
5. Lorentzon M
6. Thomas T
7. Maddox J
8. Fan M
9. Meisner PD
10. Grauer A
(2017) Romosozumab or alendronate for fracture prevention in women with osteoporosis
New England Journal of Medicine 377:1417–1427.

https://doi.org/10.1056/NEJMoa1708322
- PubMed
- Google Scholar
1. Sato M
2. Grasser W
3. Endo N
4. Akins R
5. Simmons H
6. Thompson DD
7. Golub E
8. Rodan GA
(1991) Bisphosphonate action: alendronate localization in rat bone and effects on osteoclast ultrastructure
Journal of Clinical Investigation 88:2095–2105.

https://doi.org/10.1172/JCI115539
- PubMed
- Google Scholar
1. Satterwhite J
2. Heathman M
3. Miller PD
4. Marín F
5. Glass EV
6. Dobnig H
(2010) Pharmacokinetics of teriparatide (rhPTH[1-34]) and calcium pharmacodynamics in postmenopausal women with osteoporosis
Calcified Tissue International 87:485–492.

https://doi.org/10.1007/s00223-010-9424-6
- PubMed
- Google Scholar
(2013) Coupling systems biology with multiscale mechanics, for computer simulations of bone remodeling
Computer Methods in Applied Mechanics and Engineering 254:181–196.

https://doi.org/10.1016/j.cma.2012.10.015
- Google Scholar
(2014) Mathematical modeling of postmenopausal osteoporosis and its treatment by the anti-catabolic drug denosumab
International Journal for Numerical Methods in Biomedical Engineering 30:1–27.

https://doi.org/10.1002/cnm.2584
- PubMed
- Google Scholar
(2011) Coping with time scales in disease systems analysis: application to bone remodeling
Journal of Pharmacokinetics and Pharmacodynamics 38:873–900.

https://doi.org/10.1007/s10928-011-9224-2
- PubMed
- Google Scholar
1. Siddiqui JA
2. Partridge NC
(2016) Physiological bone remodeling: systemic regulation and growth factor involvement
Physiology 31:233–245.

https://doi.org/10.1152/physiol.00061.2014
- PubMed
- Google Scholar
1. Solheim E
(1998) Growth factors in bone
International Orthopaedics 22:410–416.

https://doi.org/10.1007/s002640050290
- PubMed
- Google Scholar
(2018) The clinical potential of romosozumab for the prevention of fractures in postmenopausal women with osteoporosis
Therapeutic Advances in Musculoskeletal Disease 10:105–115.

https://doi.org/10.1177/1759720X18775936
- PubMed
- Google Scholar
1. Sowers MR
2. Zheng H
3. McConnell D
4. Nan B
5. Harlow SD
6. Randolph JF Jr
(2008) Estradiol rates of change in relation to the final menstrual period in a population-based cohort of women
The Journal of Clinical Endocrinology and Metabolism 93:3847–3852.

https://doi.org/10.1210/jc.2008-1056
- PubMed
- Google Scholar
1. Streicher C
2. Heyny A
3. Andrukhova O
4. Haigl B
5. Slavic S
6. Schüler C
7. Kollmann K
8. Kantner I
9. Sexl V
10. Kleiter M
11. Hofbauer LC
12. Kostenuik PJ
13. Erben RG
(2017) Estrogen regulates bone turnover by targeting RANKL expression in bone lining cells
Scientific Reports 7:6460.

https://doi.org/10.1038/s41598-017-06614-0
- PubMed
- Google Scholar
1. Suen PK
2. Zhu TY
3. Chow DHK
4. Huang L
5. Zheng LZ
6. Qin L
(2015) Sclerostin antibody treatment increases bone formation, bone mass, and bone strength of intact bones in adult male rats
Scientific Reports 5:15632.

https://doi.org/10.1038/srep15632
- PubMed
- Google Scholar
(1982) Parathyroid hormone stimulates the bone apposition rate independently of its resorptive action: differential effects of intermittent and continuous administration
Endocrinology 110:506–512.

https://doi.org/10.1210/endo-110-2-506
- PubMed
- Google Scholar
1. Tanaka S
2. Adachi T
3. Kuroda T
4. Nakamura T
5. Shiraki M
6. Sugimoto T
7. Takeuchi Y
8. Saito M
9. Bilezikian JP
(2014) New simulation model for bone formation markers in osteoporosis patients treated with once-weekly teriparatide
Bone Research 2:14043.

https://doi.org/10.1038/boneres.2014.43
- PubMed
- Google Scholar
1. Thouverey C
2. Caverzasio J
(2015) Sclerostin inhibits osteoblast differentiation without affecting BMP2/SMAD1/5 or Wnt3a/β-catenin signaling but through activation of platelet-derived growth factor receptor signaling in vitro
BoneKEy Reports 4:757.

https://doi.org/10.1038/bonekey.2015.126
- PubMed
- Google Scholar
1. Tu KN
2. Lie JD
3. Wan CKV
4. Cameron M
5. Austel AG
6. Nguyen JK
7. Van K
8. Hyun D
(2018)
Osteoporosis: a review of treatment options

Pharm. Ther 43:92–104.
- Google Scholar
1. Virtanen P
2. Gommers R
3. Oliphant TE
4. Haberland M
5. Reddy T
6. Cournapeau D
7. Burovski E
8. Peterson P
9. Weckesser W
10. Bright J
11. van der Walt SJ
12. Brett M
13. Wilson J
14. Millman KJ
15. Mayorov N
16. Nelson ARJ
17. Jones E
18. Kern R
19. Larson E
20. Carey CJ
21. Polat İ
22. Feng Y
23. Moore EW
24. VanderPlas J
25. Laxalde D
26. Perktold J
27. Cimrman R
28. Henriksen I
29. Quintero EA
30. Harris CR
31. Archibald AM
32. Ribeiro AH
33. Pedregosa F
34. van Mulbregt P
35. SciPy 1.0 Contributors
(2020) Author Correction: SciPy 1.0: fundamental algorithms for scientific computing in Python
Nature Methods 17:352.

https://doi.org/10.1038/s41592-020-0772-5
- PubMed
- Google Scholar
1. Weinstein RS
(2012) Glucocorticoid-induced osteoporosis and osteonecrosis
Endocrinology and Metabolism Clinics of North America 41:595–611.

https://doi.org/10.1016/j.ecl.2012.04.004
- PubMed
- Google Scholar
1. Wijenayaka AR
2. Kogawa M
3. Lim HP
4. Bonewald LF
5. Findlay DM
6. Atkins GJ
7. Carvalho DP
(2011) Sclerostin stimulates osteocyte support of osteoclast activity by a RANKL-dependent pathway
PLOS ONE 6:.

https://doi.org/10.1371/journal.pone.0025900
- Google Scholar
1. Winkler DG
2. Sutherland MSK
3. Ojala E
4. Turcott E
5. Geoghegan JC
6. Shpektor D
7. Skonier JE
8. Yu C
9. Latham JA
(2005) Sclerostin inhibition of Wnt-3a-induced C3H10T1/2 cell differentiation is indirect and mediated by bone morphogenetic proteins
The Journal of Biological Chemistry 280:2498–2502.

https://doi.org/10.1074/jbc.M400524200
- PubMed
- Google Scholar
1. Zhang L
2. Mager DE
(2019) Systems modeling of bortezomib and dexamethasone combinatorial effects on bone homeostasis in multiple myeloma patients
Journal of Pharmaceutical Sciences 108:732–740.

https://doi.org/10.1016/j.xphs.2018.11.024
- Google Scholar
(2011) General analysis of mathematical models for bone remodeling
Bone 48:910–917.

https://doi.org/10.1016/j.bone.2010.12.010
- PubMed
- Google Scholar

Article and author information

Author details

David J Jörg

Biomedical Modeling and Simulation Group, Global Research and Development, Fresenius Medical Care Germany, Bad Homburg, Germany

Contribution
Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review and editing

For correspondence
david.joerg@fmc-ag.com

Competing interests
DJJ is an employee of Fresenius Medical Care Germany. DJJ is an inventor on a patent application named "Virtual Osteoporosis Clinic" (WO 2021/231374 A1). The author has no other competing interests to declare

"This ORCID iD identifies the author of this article:" 0000-0001-5960-0260
Doris H Fuertinger

Biomedical Modeling and Simulation Group, Global Research and Development, Fresenius Medical Care Germany, Bad Homburg, Germany

Contribution
Conceptualization, Writing – review and editing

Competing interests
DHF is an employee of Fresenius Medical Care Germany. DHF is an inventor on a patent application named "Virtual Osteoporosis Clinic" (WO 2021/231374 A1). The author has no other competing interests to declare

"This ORCID iD identifies the author of this article:" 0000-0001-6321-2367
Alhaji Cherif

Renal Research Institute, New York, United States

Contribution
Methodology, Writing – review and editing

Competing interests
AC is an employee of the Renal Research Institute, a wholly owned subsidiary of Fresenius Medical Care. AC is an inventor on a patent application named "Virtual Osteoporosis Clinic" (WO 2021/231374 A1). The author has no other competing interests to declare
David A Bushinsky

Department of Medicine, University of Rochester School of Medicine and Dentistry, Rochester, United States

Contribution
Biological and clinical expertise, Writing – review and editing

Competing interests
DAB reports a grant from NIH, personal fees, stock and stock options from Tricida; personal fees from and stock in Amgen; personal fees from Relypsa/Vifor/Fresenius; personal fees from Sanifit outside the submitted work. The author has no other competing interests to declare
Ariella Mermelstein

Renal Research Institute, New York, United States

Contribution
Data curation, Writing – review and editing

Competing interests
AM is an employee of the Renal Research Institute, a wholly owned subsidiary of Fresenius Medical Care. The author has no other competing interests to declare
Jochen G Raimann

Renal Research Institute, New York, United States

Contribution
Data curation, Writing – review and editing

Competing interests
JGR is an employee of the Renal Research Institute, a wholly owned subsidiary of Fresenius Medical Care. The author has no other competing interests to declare
Peter Kotanko
1. Renal Research Institute, New York, United States
2. Icahn School of Medicine at Mount Sinai, New York, United States
Contribution
Biological and clinical expertise, Conceptualization, Investigation, Methodology, Visualization, Writing – review and editing

Competing interests
PK is an employee of the Renal Research Institute, a wholly owned subsidiary of Fresenius Medical Care. PK holds stock in Fresenius Medical Care. PK is an inventor on a patent application named "Virtual Osteoporosis Clinic" (WO 2021/231374 A1). The author has no other competing interests to declare

"This ORCID iD identifies the author of this article:" 0000-0003-4373-412X

Funding

No external funding was received for this work.

Acknowledgements

We thank Friederike E Thomasius for insightful discussions and critical comments on the manuscript.

Copyright

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.