Connecting the Dots: Linking Osteocyte Activity and Therapeutic Modulation of Sclerostin by Extending a Multiscale Systems Model.

The goal of this work was to extend a mathematical, multiscale systems model of bone function, remodeling, and health in order to explore hypotheses related to therapeutic modulation of sclerostin and quantitatively describe purported osteocyte activity within bone remodeling events. A pharmacokinetic model with first-order absorption and dual elimination pathways was used to describe the kinetics of romosozumab, a monoclonal antibody (mAb) against sclerostin. To describe total circulating sclerostin, an extended indirect response model of inhibition of offset was developed. These models were subsequently linked to the systems model, with sclerostin signaling changes in resorption and formation through established osteocyte-mediated mechanisms. The model proposes relative contributions of the osteocyte to the RANKL pool, a major player in feedback signaling, and is used to explore hypotheses surrounding attenuation of anabolic activity after multiple doses of sclerostin mAbs, a phenomenon whose mechanism is poorly understood.

Study Highlights

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC? ☑ The current systems pharmacology models that include osteocyte activity or sclerostin mAb intervention are not designed to predict quantitative clinical outcomes. • WHAT QUESTIONS DID THIS STUDY ADDRESS? ☑ Is it possible to leverage the clinical study data available with sclerostin mAbs to extend a systems model to predict responses to therapeutic modulation of sclerostin and describe osteocyte activity within bone remodeling events? • WHAT THIS STUDY ADDS TO OUR KNOWLEDGE ☑ The extended systems model can be used to examine hypotheses surrounding the mechanism for attenuation of anabolic activity after multiple doses of sclerostin mAbs. This has not been fully explored by laboratory experimentation. It is also used to investigate the relative contribution of osteocytes to feedback regulation within the bone. • HOW THIS MIGHT CHANGE CLINICAL PHARMACOLOGY AND THERAPEUTICS ☑ The extended model can be used to explore therapeutic target modulation in order to maximize and maintain increased BMD in osteoporosis patients.

Sclerostin has been identified as a target for osteoporosis treatment because preventing sclerostin inhibition of Wnt has been shown to both increase markers of bone formation and decrease markers of resorption, expanding net gain of bone calcification and increasing bone mineral density (BMD).1 This mechanism, which “decouples” bone formation and resorption, is differentiated from other osteoporosis treatment mechanisms that are either purely anabolic (both formation and resorption increase, e.g., intermittent parathyroid hormone (PTH)) or catabolic (both formation and resorption decrease, e.g., bisphosphonates, RANKL-inhibition). Furthermore, sclerostin is mainly expressed in the osteocyte, limiting off-target effects of inhibition in other tissues.

Questions remain about the mechanism of sclerostin inhibition and how this is linked to osteocyte activity and feedback regulation in bone remodeling.2 One clinical question is whether or not efficacy can be maintained after multiple doses of an anti-sclerostin monoclonal antibody (mAb). Identification of appropriate dosing regimens of sclerostin mAb and/or its combination with an antiresorptive to promote greater formation and prolonged maintenance of strong bone is a critical step in the advancement of this therapy. Potential for a model that is aimed at elucidating the mechanisms of sclerostin modulation includes exploration of dosing regimen and trial design considerations. Such inputs, although not meant at this stage to generate statistical probabilities, could generate hypotheses (learnings) for further experimental confirmation, e.g., through clinical investigation.

A multiscale bone model has been published3 that combines important aspects of three previous models of bone in order to combine quantitative aspects of bone physiology: all major organ systems involved in calcium handling,4 feedback control between osteoblasts (OB) and osteoclasts (OC) through the Receptor Activator of Nuclear Factor κ#/RANK-ligand/Osteoprogerin (RANK/RANKL/OPG) axis,5 and dynamics of intermittent PTH administration.6 In its current construct, the model lacks the osteocyte and sclerostin-related components necessary to predict effects of sclerostin mAb treatment on clinical outcomes like BMD. Other published models of sclerostin, osteocytes (OCY), and Wnt signaling are either qualitative in nature, with model variables lacking physiological meaning,7 or they are focused on mechanical strain analysis.8,9 The “strain” models depict quantitative changes in a single bone unit during loading, but they do not account for feedback signaling between bone cells, which largely contribute to signal transduction and remodeling. In contrast with other models, the multiscale bone model provides an evaluated platform to predict changes in BMD based on clinical markers of formation and resorption. It has already been used to predict changes in BMD after treatment with denosumab.10 The new model components were developed by incorporating knowledge of the Wnt/β-catenin signaling and its role in bone formation by leveraging data from recent clinical studies.1,11,12 The updated model promotes understanding of how OCY signals contribute to remodeling within the bone and how sclerostin mAbs can be used to harness these signals to maximize bone formation in patients with osteoporosis.

MATERIALS AND METHODS

Data

A phase I study reporting time–concentration profiles of romosozumab after a single dose1 was used to build the pharmacokinetic (PK) model. Total sclerostin concentrations measured in two phase I studies over a range of single and multiple doses of blosozumab12 were used to estimate parameters in the pharmacodynamics (PD) model. C-terminal telopeptide (CTx) and procollagen type 1 N propeptide (P1NP) data from these three studies and two additional studies11,13 were used to build the sclerostin-related components and BMD changes into the model. A sixth was used as a qualification dataset.14

The bone formation marker, serum P1NP, as a marker of bone formation, has replaced serum bone-specific-alkaline-phosphatase (BSAP) in newer clinical trials. Since the existing model only described changes in BSAP, a regression model was developed to calculate P1NP as a function of BSAP. Data collected from 21 different studies (Supporting Information refs. 4-22) identified by a PubMed search conducted on or around March 6, 2014, using keywords “osteoporosis” or “postmenopausal,” and “alkaline phosphatase,” or “ALP,” “amino-terminal propeptide,” or “P1NP,” and “clinical trials,” and “humans.” Graphical presented data from these publications were digitized using GraphClick (v. 3.0 Arizona Software).

PK/PD model

A PK model was used to describe nonlinear kinetics of sclerostin mAbs and this was linked to a PD model describing changes in levels of circulating sclerostin. Typical values for PK parameters were estimated using a maximum likelihood (MLE) approach, based on the reported treatment arm-level data. Unexplained residual variation was described with a proportional residual variance model, with residual random effect assumed to be normally distributed with mean zero and variance, σ.2 Four subcutaneous dose levels of 1, 3, 5, and 10 mg/kg and two intravenous (i.v.) dose levels of 1 mg/kg and 5 mg/kg were used to fit a two-compartment model with first-order absorption and parallel linear and nonlinear clearance pathways. An indirect response PD model was used to describe changes in total circulating sclerostin. The PK/PD parameters were estimated using nonlinear-mixed effects models (NONMEM software, v. 7.2, ICON Development Solutions, Hanover, MD).

Translational relationship between formation markers

Regression models were used to explore the relationship between time-matched BSAP and P1NP, both normalized to baseline. Model development considered linear and nonlinear relationships between P1NP and BSAP. Inter-arm variance for BSAP was weighted by the inverse of the sample size for each arm15: where ε represents the random effect and n = the number of subjects contributing to the data point for the ith arm.

A random effect was applied to the slope term in linear and nonlinear models to account for longitudinal differences between arms. Model selection was based on successful model minimization and the Akaike Information Criteria (AIC).16

Multiscale model expansion to include osteocyte function and sclerostin modulation

Code to represent the multiscale bone model has been developed in R.17 Simulations were performed using the DLSODA differential equation solver provided through ODEPACK within C++ interfaced to R (“Rcpp” package18). Parameters were coarsely estimated by tuning individual parameters involved in indirect changes in OB and OC over physiological ranges, solvable by the ODEsolver. Model performance was evaluated by visual inspection, compared to the clinical data. After model structure and initial parameter estimates were in place, the parameters were optimized, individually or in groups, using the R package “minqa.”19 This is a derivative-free optimization algorithm by quadratic approximation, used to minimize an ordinary least-squares objective function. BMD-related parameters were tuned and optimized last. The order in which parameters were optimized and the data used to optimize each parameter is summarized in Table 3.

Table 3

Parameters describing changes in turnover makers and BMD associated with sclerostin mAbs

Order	Parameter	Definition	Value (units)	Data used to estimate
1	γROB	Exponent on sclerostin effect on differentiation rate of OB	0.0703 (unitless)	P1NP
2	EMAXSCLER	Max sclerostin effect on OB	5.22 (hours)	P1NP
2	γSCLEROB	Hill coefficient for sclerostin effect on OB	1.15 (unitless)	P1NP
2	EC50SCLER	Half-max effect of sclerostin on OB	667 (hr−1)	P1NP
3	kout_T	Depletion rate on translation compartment	0.00692 (hr−1)	P1NP
Fixed*	FracOCY	Fraction of OB becoming OCY	0.50 (unitless)	–
Fixed*	ρ	Exponent on osteoblast effect on RANKL	0.0100 (unitless)	–
4	δ	Exponent on osteocyte effect on RANKL; scales PTH effect on RANKL	0.0592 (unitless)	P1NP, CTx
4	γOCY	Exponent on sclerostin effect on OCY apoptosis rate	0.405 (unitless)	P1NP, CTx
5	γOPG	Exponent on sclerostin effect on OPG	1.61 (unitless)	CTx
Fixed*	kout_BMDdel	Delay parameter on anabolic effects on BMD	0.001 (unitless)	Lumbar spine BMD
6	γOB LS	Exponent on anabolic contribution to overall lumbar spine BMD	0.711 (unitless)	Lumbar spine BMD
6	γOB FN	Exponent on anabolic contribution to overall femoral neck BMD	0.0934 (unitless)	Femoral neck BMD
6	γOB TH	Exponent on anabolic contribution to overall total hip BMD	0.686 (unitless)	Total hip BMD

*These parameters were shown to be relatively insensitive within the tested range so these values were fixed prior to optimization.

Qualification of expanded multiscale model

To evaluate model structure, simulated data were plotted against the model-building dataset.1,11–13 After the final model structure was in place and the univariate parameter optimization step was performed, a local sensitivity analysis was conducted on all of the 23 estimated sclerostin-related parameters (Tables 2, 3; Supplementary Figure 1). This analysis determined the influence of each parameter estimate on clinical endpoints, since simultaneous optimization in this context was not possible. Monte Carlo simulations over a range of parameter estimates spanning ±0.8 * the value of the final parameter estimate were plotted against time-matched changes in P1NP and CTX from a new clinical study,14 which had not been used for the initial parameter estimation. This same analysis was also performed for each of the lumbar spine (LS) and total hip (TH) BMD parameters (femoral neck (FN) BMD was not reported in the new study).

Table 2

Six candidate models describing the relationship between time-matched BSAP and P1NP

Model	Model equation	Random effects	AIC
Linear		On A	5,035.57
Exponential		On A	4,904.99
Emax		On	5,036.15
Sigmoid Emax		On	4,707.93
Sigmoid Emax + Intercept		On	4,710.38
Quadratic		On A	4,721.53

Data from 21 clinical trials were evaluated. Random effects were added on the “slope” parameters in each model to account for variability between arms.

Simulations were performed to verify that changes made in the model did not negatively impact the denosumab therapy simulations from previous works.10 The simulated denosumab dataset contained nine treatment arms with doses taken from several clinical trials of denosumab: placebo, 6, 60, 140, 100, and 210 mg at dosing intervals of 3 and 6 months.20–24 Simulated P1NP and CTx concentrations for each of these arms were overlaid with time-matched observations and validated by visual inspection (Supplementary Figure 2). Dose-matched simulations were performed to explore the impact of dosing interval on clinical outcomes. The final model was evaluated according to the criteria outlined by Agoram25 for large systems pharmacology models.

RESULTS

The PK of rosozumab and blosozumab were used to drive the PD response in sclerostin (Figure 1). Identifiability of PK and PD parameter estimates for each drug was dictated by data availability. Currently, public-source data included only PK data for romosozumab, but not blosozumab, and only sclerostin data for blosozumab, but not romosozumab. Therefore, the PK parameters were estimated from the antibody concentration–time profiles of romosozumab1 and the PD parameters were estimated from blosozumab12 (Table 1).

Figure 1

(a) The PK model describes circulating sclerostin mAb concentrations that drive changes in circulating total sclerostin protein (SCLER) in the PD model. C represents drug concentration; At is the peripheral compartment; sc = subcutaneous; iv = intravenous; Vc, Vp, CL, Kin, Kout, Vmax, Km are model parameters defined in Table 1. (b) Parameter estimates for the PK model were generated by fitting the model to romosozumab data, and parameter estimates for the PD model were generated by fitting the model to blosozumab data.

Table 1

Estimated PK/PD parameters for antibodies against sclerostin

Parameter	Value (95% CI)
Absorption rate constant ka	0.187 (0.142, 0.233) day−1
Linear clearance, CL	0.254 (0.228, 0.281) L/day
Maximum elimination rate constant, Vmax	5.87 (2.49, 9.26) L/day−1
Michaelis-Menten constant, km	0.423 (0−1.64) nM for blosozumab and 9.93 (0.777−1.03) nM for romosozumab
Volume of the central compartment, Vc	2.9 (2.31, 3.48) L
Volume of peripheral compartment, Vp	3.29 (2.43, 4.16) L
Intercompartmental clearance, Q	0.467 (0.326, 0.609) L/day
Bioavailability, F	0.904 unitless
Synthesis rate constant, kin	3.68 (0−15.3) nM/day
Degradation rate constant, kout	26.0 (0–95.9) day−1
Internalization rate constant, k0	0.195 (0.0349–0.356) day−1

CIs were calculated as the estimate ± 1.96* asymptotic standard error of the estimate; symmetric CIs were truncated at 0 for rate constants km, kin, and kout.

For the PD estimates, the PK parameters, including the V parameter describing binding kinetics, were assumed to be the same for both mAbs and the PD parameters were subsequently estimated. In this sequence, the K parameter in the nonlinear clearance component of the PK model was estimated based on romosozumab concentration and was estimated separately for blosozumab using the sclerostin response. Therefore, all the PK parameters derived for romosozumab were fixed and used to simulate mAb concentration in the PD model, with the exception of k, which was reestimated for blosozumab. The rationale for this is discussed further in the Supporting Information.

Various linear and nonlinear relationships were explored to describe the BSAP-P1NP translation (Table 2). The lowest calculated AIC value resulted from the sigmoid E model (AIC = 4707.93), and this model structure was carried forward. An intercept parameter of int = 20.4181 was also added to this model within the multiscale bone model, so that at the initial conditions, percentage of baseline BSAP = P1NP = 100%.

The estimated parameters for this model were E = 2,050, 95% CI: (295, 3,800), γ = 1.8, 95% CI: (1.52, 2.08), and EC50 = 467, 95% CI: (128,805). ΩA = 2.27, 95% CI: (0.376–4.16), %CV = 38.2, σ = 0.515, 95% CI: (0.0756–0.955), %CV = 18.0, for n = 50 subjects.

Developing model structure to describe changes in turnover markers and BMD

Based on supporting literature, six points of intersection were identified within the multiscale bone model where changes in sclerostin have a known effect. These are (i) the depletion rate of pre-osteoblasts (ROB); (ii) the formation rate of OB; (iii) the rate of OCY apoptosis; (iv) the level of OCY effect applied to RANKL; (v) the accumulation of OPG; and (vi) the differential effects of sclerostin on regional changes in BMD. See Figure 2 to see schematically how these pieces fit into the model structure. Because identifiability of new model parameters is problematic if all are estimated simultaneously, P1NP-associated parameters were estimated first, followed by those parameters affecting CTx. This follows our understanding of the mechanism of Wnt pathway upregulation; β-catenin accumulates and causes changes in transcription factors regulating differentiation of pre-OB to OB. This effect is upstream, or takes place ahead, of modification of bone resorption through this pathway, which is largely controlled by signaling through RANK/RANKL/OPG in OB26 and OCY.27

Figure 2

Schematic of the bone-remodeling systems model. Intersection points of sclerostin signaling effects within the model are identified with numbers corresponding to description in the text. New model compartments are indicated with white text and shading and corresponding equation numbers.

To apply a sclerostin signal to the OB, two methods were tried. First, a “precursor pool” compartment was applied upstream of the ROB. This physiologically represents cells earlier in the OB lineage being recruited for differentiation which cause an increase in formation activity, as the flux of substrate moving through the OB differentiation pathway increases.28 The precursor-pool approach resulted in a very rapid depletion of the ROB pool and a peak P1NP concentration–time profile that was right-shifted relative to the true response. An alternative way to signal the flux of substrate through the OB pathway was to use a direct approach whereby a sclerostin effect was applied to the depletion rate of ROB and a separate effect was applied directly to the OB accumulation rate. Together, these increased OB activity, but caused an attenuation of OB activity over time as the ROB became depleted, which is consistent with the pattern of anabolic turnover markers seen after multiple doses of sclerostin mAb.12 A “translation” compartment was also added in the direct approach to delay the effect of sclerostin on the increase in OB activity only slightly, relative to the delay that was elicited by the precursor-pool approach. Together, these additions represent the first two points of interaction in the modified model.

OCY are the major source of sclerostin,2 and the cell type responsible for sclerostin modulation in response to mechanical stimulation.29 Studies have also shown that OCY are a major contributor to RANKL trafficking and regulation of osteoblast- and osteoclast-ogenesis at the level of the RANK/RANKL/OPG axis.30 Membrane-bound RANKL is provided by dendritic processes on the osteocyte to signal to osteoclast precursors and upregulate osteoclastogenesis. Conversely, upregulating the Wnt pathway can suppress osteoclast-mediated bone resorption through reduced expression of RANKL.26,31 For the third and fourth points of intersection, a sclerostin effect was applied to the depletion rate (or rate of apoptosis) of OCY and an osteocyte effect was applied to the RANKL production rate. The combined result is that when sclerostin mAb is dosed, OCY activity is purportedly decreased, concurrent with diminished SOST expression29 and RANKL production falls as well, gradually inhibiting osteoclastogenesis.

For the fifth point of intersection, a direct sclerostin effect, was applied to the production rate of OPG to simulate the effect of Wnt-pathway regulation of OPG secretion. In the literature, this level of regulation was demonstrated by knocking out osteocytic β-catenin in mice. These mice became osteoporotic due to increased numbers of OC and increased bone resorption27 and had diminished expression of osteocytic OPG, demonstrating that signaling through β-catenin in the osteocytes is required for regulation of resorption activity. Modeling a direct sclerostin signal on OPG production rate had the desired effect of rapidly decreasing the simulated CTx levels, immediately following a dose of sclerostin mAb.

To model sclerostin antibody effects on FN, TH, and LS BMD, the same model structure as was used to describe denosumab effects on BMD10 was adapted, by reestimating the power term reflecting the impact of OB activity on BMD. The predictions resulting from this model yielded a peak increase in BMD that occurred earlier than it did in clinical reports. In order to correct for this, a delay compartment was used to adjust the rate of anabolic bone formation under the premise that it takes more time for OB to become embedded in the bone matrix, calcify, and contribute to an increase in overall BMD than it does for catabolic (OC-driven) impact on BMD.

Model equations and parameter estimates describing turnover markers and BMD

For all of the intersection points added to the model, the sclerostin effect was applied as a power model in the form of: where γ = a parameter modulating the sclerostin effect.

In general, the steady state assumption used in the original publication3 was used and compartmental equations were solved at initial conditions to minimize identifiability problems for estimated rate constants. The OB effect took the form:

OBfast and OBslow represent the two pathways of active OB elimination described in the original model publication.3

The osteocyte compartment was represented as: where Frac represents the rate of conversion of OB to OCY and

The osteocyte effect took the same form as the sclerostin effect: where δ = a parameter modulating the osteocyte effect.

The final BMD model for sclerostin has separate compartments for TH, FN, and LS BMD, where the kout and γ were estimated separately for each region and each region's set of compartmental equations took the form:

Because simultaneous parameter optimization was not possible, considerations with regards to the parameters describing the turnover markers and BMD included the order in which they were implemented and optimized, and the data used to fit each parameter (Table 3).

Qualification

In the model qualification step, simulated P1NP and CTX were consistent with changes in the turnover markers from three clinical trials with romosozumab and two clinical trials with blosozumab (Supplementary Figure 3). TH, FN, and LS BMD profiles from these same studies were also well described by the model (Supplementary Figure 4). A qualification dataset from a phase II clinical trial with blosozumab in postmenopausal women was simulated using the final model and changes in turnover markers and LS and TH BMD were plotted (Figure 3). The predicted mean change from baseline for the 180 mg dosed Q2W arm at 4 and 52 weeks were 229% (observed median; 95% CI: 231; 199–262%); for P1NP and 73.5% (83.9; 73.7; 94.2); for CTX and 152% (90.2; 81.4–99.0%); for P1NP, and 79.2 (84.6; 69.7–99.6%); for CTX, respectively. The predicted mean change from baseline at 52 weeks (180 mg dosed Q2W) for LS BMD was 16.5% (observed; 95% CI: 14.9; 12.6–17.1%) and 6.8% (4.5; 3.2–5.8%) for TH BMD. Despite the small dataset, for most data points the model predicted the change in the clinical endpoint that fell within the 95% CI reported in the literature.

Figure 3

Simulated P1NP (a), CTx (b), lumbar spine BMD (c), and total hip BMD (d) (blue line) overlaying data from a clinical trial with blosozumab (red points). This qualification dataset was not used in constructing the model (n = 29, 31, 30, 30 for arms PBO, 180 mg Q2W, 180 mg Q4W, and 270 mg Q2W, respectively).

In order to demonstrate how the model can be used to investigate the role of dosing protocol on clinical outcomes, dose-matched administrations of sclerostin mAb were simulated at several dosing intervals (Figure 4). Simulations of larger dosing intervals promote greater increases in P1NP (Figure 4a), due to precursors also achieving higher levels (Figure 4b). Maximum simulated resorption activity, however, is also increased with a larger dosing interval (Figure 4c), resulting in lower increases in total hip BMD when compared to smaller dosing intervals (Figure 4d).

Figure 4

Simulations of dose-equivalent sclerostin mAb administered at different dosing intervals demonstrate that changing the dosing interval may alter clinical outcomes. P1NP (a), osteoblast precursors (b), CTx (c), total hip BMD (d). Colored marks indicate the final dose in each arm.

DISCUSSION

Despite having a limited availability of clinical data, the updated model was able to capture the central tendency of the changes in turnover markers and BMD after single and multiple doses of sclerostin mAb, over a large dose range. The model, however, greatly underpredicted the peak percentage change of P1NP in the highest dosing groups with blosozumab (Supplementary Figure 1, arms 16 and 23). This was a consequence of the blosozumab data showing that a single 750 mg dose did not greatly increase sclerostin change from baseline over a single dose of 225 mg. The 270 mg dose elicited much greater changes than the 750 mg dose with the same dosing interval (see ref.12, fig 4). Because the model relies on the change in sclerostin as the signal for increase in P1NP, this apparent discrepancy can be seen in the model predictions at the 750 mg dose level, and to some degree in the 540 mg dosing arms.

There are many regulating elements of cellular feedback that are not explicitly defined in the model. This may explain why the model failed to recapitulate some of the dynamics of the CTX-time profile, which was also highly variable between arms. Recent data has shown that OCY and OB can also control bone remodeling directly through bone morphogenetic protein (BMP), which can affect endogenous sclerostin, independent of Wnt signaling,32 but sclerostin is also a potent antagonist of some BMPs.33 This feedback mechanism needs to be fully elucidated by experimental data before implementation into the model. OCY undergoing apoptosis also secrete interleukin (IL)−1 and IL-6, which promote osteoclastogenesis34 and provide another layer of feedback control that is not yet represented in the model. OB are also a source of IL-6, as well as macrophage-colony-stimulating factor (M-CSF), which also contributes to osteoclastogenesis.34 Supporting data are needed to fully characterize effects of sclerostin mAb on changes in CTx, as even the PBO arms show relatively large fluctuations in CTx over several months.1,11,12

It was hypothesized that diminishing the pre-OB pool with sclerostin inhibition contributes to attenuation in bone-formation activity after multiple doses. A consequence of this from a treatment perspective is that if a sclerostin mAb is dosed with a small dosing interval, and the mesenchymal-derived pre-OB do not have sufficient time to be replenished, the anabolic effect of the mAb will be diminished with repeated doses. Indeed, simulations of sclerostin mAb administration at equivalent dose levels, but at different dosing intervals, demonstrate that higher P1NP levels can be achieved by dosing in longer intervals (Figure 4a), because levels of pre-osteoblasts have time to recover between doses (Figure 4b). However, achieving the greatest magnitude of anabolic activity does not necessarily translate into larger gains in BMD, due to feedback signaling to OC. Simulations also show that at longer dosing intervals peak levels of resorption activity are also higher than at smaller dosing intervals, blunting the net gain in total BMD that can be achieved within the dosing interval (Figure 4c,d). This finding suggests that osteocyte-controlled feedback involved in sclerostin modulation has a considerable impact on turnover markers and that dosing protocol of sclerostin mAbs impacts these markers as well as changes in BMD. This model therefore offers a potential platform for simulating different regimens in order to explore the effect of dose strength, interval, and/or combination with other treatment combinations. The multiscale functionality enables consideration of feedback effects, for example, and how these mechanisms relate to BMD changes over a prolonged treatment period.

The model underpredicts BMD changes at the lumbar spine and total hip after multiple doses in the low dosing arm (Supplementary Figure 2A,B, arms 32, 33). This is likely due to the simulated peak P1NP responses, which reflect OB activity, at these later timepoints falling below baseline values (see Supplementary Figure 1A). Although it is difficult to assess what the true peak P1NP values are under this dosing regimen, it appears that P1NP response is attenuated after multiple doses in this arm, for reasons discussed previously. It is unclear what is responsible for maintaining the increase in BMD seen in the clinic at these later timepoints. It is possible that the decreased rates of resorption play more of a role at these later timepoints and account for a greater amount of the increase in BMD than the model currently predicts.

Another interesting finding of this work is the interaction of sclerostin modification with the RANK/RANKL/OPG axis and resulting dynamic changes in resorption/OC activity. The model suggests the immediate and significant decrease in CTx results from osteocytic Wnt signaling, which is known to increase OPG directly.33 This effect is coupled with a more gradual decline in CTx due to a drop in osteocytic contribution to RANKL. Finally, a mechanism in the model was needed to delay the anabolic effects of sclerostin inhibition on overall BMD. We conjecture that it takes more time for OB to embed themselves in the matrix and form new bone than it does for a decrease in resorption to produce an increase in BMD. This finding may help efforts to characterize the effects of other anabolic therapies like PTH and teriparatide on regional changes in BMD within the modeling framework.

Although simultaneous parameter fitting and full parameter identifiability could not be accomplished in this context, the updated model fulfilled the Agoram25 criteria: (i) Fit for purpose: only compartments and parameters that were necessary to improve predictions of specific measurable clinical endpoints were added to the model. (ii) Justification of model structure: every point of intersection for a sclerostin effect has an experimental or clinical basis in the published literature. (iii) Evaluation of component submodels: models were evaluated by goodness-of-fit and individual parameter sensitivity analyses were conducted. (iv) Qualification of the emergent properties of the system: the final model was cross-validated with denosumab therapy simulations to ensure the changes in turnover markers were consistent with what was reported in previous work.10 The model also predicts a significant increase in circulating PTH, after sclerostin mAb administration, which is consistent with clinical reports.12,14

CONCLUSION

A previously published multiscale model of bone function, control, and health was extended to include OCY, sclerostin, and effects of the upregulated Wnt signaling pathway elicited by sclerostin mAb administration. The model was used to simulate P1NP, CTX, and BMD profiles that resembled those in clinical studies with sclerostin mAb treatment. The utility of the model to explore biological implications of Wnt pathway modification and the role of the OCY in bone remodeling was demonstrated. Model simulations point to differential effects of osteocyte-driven feedback-driven increases in resorption activity and suggest possible benefits of shorter dosing intervals in future clinical trials with sclerostin mAbs. The model also suggests osteoporosis patients may benefit from combination therapy to mitigate these feedback effects. The extended model offers a broader in silico model-based platform to explore therapeutic target modification towards the goal of more efficiently addressing clinical development considerations and improving long-term outcomes.