Application of a Systems Pharmacology-Based Placebo Population Model to Analyze Long-Term Data of Postmenopausal Osteoporosis.

Osteoporosis is a progressive bone disease characterized by decreased bone mass resulting in increased fracture risk. The objective of this investigation was to test whether a recently developed disease systems analysis model for osteoporosis could describe disease progression in a placebo-treated population from the Early Postmenopausal Intervention Cohort (EPIC) study. First, we qualified the model using a subset from the placebo arm of the EPIC study of 222 women who had similar demographic characteristics as the 149 women from the placebo arm of the original population. Second, we applied the model to all 470 women. Bone mineral density (BMD) dynamics were changed to an indirect response model to describe lumbar spine and total hip BMD in this second population. This updated disease systems analysis placebo model describes the dynamics of all biomarkers in the corresponding datasets to a very good approximation; a good description of an individual placebo response will be valuable for evaluating treatments for osteoporosis.

RCT Entities:   Population Interventions Outcomes

Study Highlights

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC? ☑ A complete mechanism-based model describing osteoblast and osteoclast activity has been reduced in order to apply it in a population approach. It was shown that the model reduction did not jeopardize the dynamical properties of the model. The reduced model was successfully applied to describe responses in treatment with various doses of tibolone and/or calcium in postmenopausal women. • WHAT QUESTIONS DID THIS STUDY ADDRESS? ☑ Can the current disease systems analysis model be used to describe the disease progression in another external placebo-treated population? • WHAT THIS STUDY ADDS TO OUR KNOWLEDGE ☑ The developed placebo model could describe the data of this new population to a very good approximation, showing the strength of a population systems pharmacology approach. • HOW THIS MIGHT CHANGE CLINICAL PHARMACOLOGY AND THERAPEUTICS ☑ The effect of a drug involves understanding of the progression of the disease and the effect of placebo. The developed model will therefore allow for a better understanding of disease progression of osteoporosis in treated and untreated patients.

Osteoporosis is defined as a progressive systematic skeletal disorder that is characterized by the loss of bone tissue, disruption of bone architecture, and bone fragility that leads to an increased risk of fractures.1 Removal and formation of bone occurs in a continuous remodeling cycle, which is a carefully regulated process involving many local and systemic factors.2 Both men and women lose bone mass when they age. However, in postmenopausal women, mainly due to the decline in estrogen levels during menopausal transition: 1) the activity of osteoclasts (cells responsible for resorbing mineralized bone) increases compared to the activity of osteoblasts (cells responsible for synthesis of new bone matrix); 2) osteoblast activity also increases, but with a delay. The combination of the two processes result in higher bone turnover. The rapid (menopause) and slow (aging) components of disease progression underlie elevated fracture rates in vertebrae and hip joints, increased mortality, and significant healthcare costs.3,4 Various treatments exist that aim to either promote bone formation or inhibit bone resorption, yet postmenopausal osteoporosis still represents an ongoing clinical challenge with decreases in quality of life and increase in healthcare costs. Additionally, given the rapidly aging population (both men and women) the burden of this disease is expected to increase.

Disease systems analysis models are used to study treatment effects on disease progression within the context of interacting networks.5 By doing so, these models allow for the comparison of effectiveness and safety of different treatment options. In order to estimate the actual treatment effect it is essential to account for the time course of the drug and/or placebo effects and disease progression.5,6 Although it is sometimes difficult to differentiate the placebo response and the underlying disease progression, this distinction is important for designing clinical trials and interpreting results. Furthermore, one of the premises of mechanism-based disease systems analysis models is that the model structure and the values of the system-specific parameters would ultimately be independent of the type of treatment. Therefore, in principle, they provide the opportunity to predict disease progression in another comparable population than the population that was used for model development and parameter fitting. In this study we will use the term “disease progression” for disease progression of the placebo arm of the Early Postmenopausal Intervention Cohort (EPIC) study in which all subjects took (at least) 500 mg calcium per day.

The objective of this study was to test whether a current disease systems analysis model for postmenopausal osteoporosis can adequately describe disease progression in a placebo-treated population of postmenopausal women. To this end, we focused on the contribution of the postmenopausal transition in combination with the placebo effect (i.e., calcium effect) on the osteoblast and osteoclast activity. For the current model, the placebo effect captures the response to daily calcium intake and the underlying disease progression. Variation in placebo response within and among clinical trials can substantially affect conclusions about the efficacy of new medications. Developing a robust population model to describe the placebo arm is of interest to enable quantification of drug effects, and eventually to guide the design of clinical trials for osteoporosis treatment. This outcome of the placebo arm is the result of individual placebo response (including the eventual effect of calcium intake) and progression of the disease.

For osteoporosis, different models have been proposed. The group of Peterson and Riggs have published a physiology-based mathematical representation of integrated calcium homeostasis and bone remodeling,7 which they later extended to show nonlinear changes in lumbar spine (LS) bone mineral density (BMD) upon denosumab treatment, a monoclonal antibody that decreases bone resorption,8 and to describe the impact of progressive loss of kidney function.9 This group has also described estrogen-related changes in bone and calcium balance through menopause transition10,11 with an extension to fracture risks.12 In their current form these extensive physiology-based models are challenging to apply in a population approach, which captures the existing variability in addition to the structural dynamics. In order to do so, Danhof and colleagues have shown that a full mechanism-based model of interacting osteoblast and osteoclasts can be reduced to a simpler model, which can describe the dynamics of the full model13 to a very good approximation (for details on the simplification, see Ref.14). The reduced model was then applied to clinical data from postmenopausal women receiving calcium and tibolone.15 This reduced systems disease model ensures parameter identifiability and could describe the dynamics of biomarkers that respond at widely different timescales to drug treatment. Even though the Peterson and Riggs model describes cellular through organ level aspects involved in bone and calcium homeostasis, it also contains a high number of differential equations and parameters to be estimated. Due to this model complexity, we will here focus on the disease systems analysis model developed by Post et al.15

In this study, we applied an osteoporosis disease systems analysis model15 to the placebo arm of a clinical study from the EPIC study.16 It should be noted that the EPIC study was intended as an osteoporosis prevention study.17 As such, patients enrolled were thought to be at risk of osteoporotic bone loss, but they did not meet the usual criteria for low bone BMD that would define osteoporosis (see Methods). The range of “years since menopause” (YSM) in the EPIC study population (0.5–27 years) is larger than in the population of the tibolone study (1–4 years). A characteristic of the disease systems analysis model is that the timescale of the disease process is incorporated rather than the timescale relative to the start of the study. Therefore, this systems pharmacology approach allowed us to 1) qualify the disease systems analysis model using an external dataset; 2) develop a robust placebo model; 3) challenge the model with data over a longer time after the onset of the disease (e.g., the start of the menopause); 4) relate various timescales (short- and long-term biomarkers); and 5) directly address and describe the covariate YSM instead of it being a covariate at baseline. The latter was possible because YSM and not study time was used as the timescale in the model.

METHODS

Subject population and study design

Data for the current analysis were obtained from the EPIC study, a clinical trial of oral alendronate in 1,609 postmenopausal women who were randomly assigned in a double-blind manner to receive alendronate, placebo, or open-label estrogen-progestin in order to evaluate the potential to prevent osteoporosis. In this study we only used data from the placebo arm (n = 470). Briefly, all participants were between 45 and 59 year of age and at least 6 months past menopause at baseline, were in good general health, and had no clinical or laboratory evidence of confounding systemic disease.1 Four study centers (two in the United States and two in Europe) were involved in this trial. To ensure that most women who entered the study did not yet have osteoporosis, only 10% of the women enrolled at each center were allowed to have an LS-BMD below 0.8 g/cm2, as measured by dual-energy x-ray absorptiometry. All women adhered to therapy (had taken at least 80% of the prescribed number of tablets, confirmed by tablet count). Dietary calcium intake was estimated at baseline and annually during the study on the basis of a food-frequency questionnaire.

Measurement of BMD and biochemical markers of bone turnover

BMD of LS and total hip (TH, defined as the femoral neck plus trochanter and intertrochanteric area) was measured by dual-energy x-ray absorptiometry (model 2000, Hologic, Waltham, MA) twice at baseline and annually thereafter. BMD is reported in g/cm2. Blood and morning second-void urine samples were collected after an overnight fast at baseline and every 6 months thereafter. Bone resorption was estimated by using urine N-telopeptide cross-links of type I collagen (NTX) (Osteomark, Ostex, Seattle, WA). NTX is reported as nmol bone collagen equivalents (bce) and corrected for creatinine excretion (nmol bce/mmol cr). In addition, the serum level of bone-specific alkaline phosphatase (BSAP) was measured at baseline and at months 6, 12, 24, 36, 42, and 48 in a random sample of 205 women to estimate bone formation (Ostase, Hybritech, San Diego, CA). BSAP is reported in ng/mL.

Mechanism-based disease systems analysis model

In an earlier publication it was demonstrated that the completely mechanistic bone cell interaction model proposed by Lemaire et al.,13 which involves responding but not yet active osteoblasts (R), active osteoblasts (B), and active osteoclasts (C), can be mathematically reduced while maintaining its dynamic properties.14 The reduced model is shown in Figure 1 and only involves active B and C, which is given by the following system of equations: where y and z are defined as y=B/B0 and z=C/C0 or, in other words: the activity of the osteoblast (B) and osteoclast (C) relative to their baseline values, B0 and C0, respectively. kB is the apoptosis rate of active osteoblasts. Cs is the value of C for which approximately half of the transforming growth factor-β (TGF-β) receptors are occupied.13 The function f(t) models disease progression. We used the function as proposed by Post et al.:15 where kestrogen is the first-order rate constant of disease progression. In Eq. 1 PCa is the placebo effect that starts at tstart and the effect of placebo is assumed to wear off over time with rate constant kCa: where Pmax is a measure for the calcium-induced inhibition of the RANK- RANKL-OPG through parathyroid hormone (PTH),7,13,18 which we take to be unity (100%).

Figure 1

Schematic representation of the mechanism-based disease systems analysis model. Active osteoblast and osteoclast cells and the indicated interactions form the mechanism-based core (shown in gray) of this model, which are linked to the biomarkers, NTX, BSAP, LS-BMD, and TH-BMD as shown in the dotted area. PTH stands for parathyroid hormone, TGF-β for transforming growth factor-β, OPG for osteoprotegerin, RANK for receptor activator of NF-κB, and RANKL for receptor activator of NF-κB ligand. RANKL binds to RANK and promotes osteoclast differentiation, while OPG inhibits this differentiation by binding RANKL. Figure and legend were adapted from Ref. 15.

The mechanistic core model (Eq. 1) is linked to the corresponding biomarkers (Figure 1) using the following transducer functions:

These transducer functions involve 1) a baseline parameter (NTX0, BSAP0), which links the relative cell activity to the baseline value of the marker, and 2) a positive transduction parameter (ρ), which links relative changes in cell activity to those in the corresponding bone turnover markers.

In comparison to the biomarkers NTX and BSAP, changes in LS-BMD and TH-BMD are relatively slow and have dynamics and timescales of their own, which were described by the zeroth-order process (15):

At baseline (healthy bone status), y = 1 and z = 1 and hence S = 1. In light of Eq. 5 this means that if y and z stay at baseline, then LS-BMD and TH-BMD do not change with time. However, if either or both y and z change the overall change in BMD depends on S.

Body composition is known to induce changes in bone morphology.19 Body mass index (BMI) was incorporated as a fraction of the baseline of LS-BMD (LS-BMD0) and baseline TH-BMD (TH-BMD0) using the median BMI of 25.2 kg/m2.

In order to initialize the model at a healthy normal state (z, y, and S = 1), individual timescales were normalized using time-since-onset-of-menopause as the characteristic time frame (Figure 2). Due to this normalization, individual subject's disease trajectories are harmonized to an identical start point and allowed to set the onset of menopause for the entire population at t = 0. When combining the individual disease trajectories in a population approach, the overall population's disease trajectory can then be defined and the need for a disease status covariate on baseline is not required as its influence on the variability is already described.

Figure 2

Normalization of individual timescales based on years since menopause (YSM). (a) Typical imaginary disease progression curve for osteoporosis (black solid line). Disease progression is plotted as relative fraction of healthy status (=1) vs. years since menopause. Five subjects (no real data) are shown that have different YSM at the start of the study. For these subjects we assume that they complete the study duration of 4 years as indicated by the colored lines. At the start of the study, these five subjects are thus at different locations on the disease progression curve. The shaded gray area represents the exclusion criteria for YSM. The range of YSM for EPIC 1 (1≤YSM≤5) and for EPIC 2 (all YSM) is indicated by the black arrows. (b) Normalization based on study time for the five subjects (same colors as in a). This normalization harmonizes the individual subjects disease trajectory to an identical study starting point.

Data analysis

Description of the data analysis, including the model file, can be found in the Supplementary Information.

RESULTS

Comparison of baseline demographic characteristics

Before applying the model, the baseline demographic characteristics for the placebo arms of the tibolone (originator) and EPIC studies were compared (Table 1). Data from the EPIC study were split into two different datasets: EPIC 1, which only includes women with YSM between 1 and 5 years (n = 222), and EPIC 2, which contains the entire range of YSM (n = 470). EPIC 1 was created to have demographic characteristics as close as possible to the tibolone study and was used for qualification of the disease systems analysis model. EPIC 2 was used to challenge the model for a longer time after the onset of the menopause. This led to mean YSM for EPIC 1 and EPIC 2 that were, respectively, comparable and higher when compared to the tibolone study. Furthermore, mean age was somewhat lower in EPIC 1 but similar for EPIC 2 in comparison to the tibolone study. There was no meaningful difference between the two clinical studies in terms of BMI. Additionally, for EPIC 1 and EPIC 2 we observed a lower mean value for LS-BMD and TH-BMD. The biomarker for resorption, NTX, was higher in the two EPIC studies than in the tibolone study. NTX, a specific collagen degradation product, is produced by osteoclastic bone resorption. When bone resorption is accelerated, which is known to happen in response to estrogen deficiency,20 NTX levels in urine increase as well. BSAP, on the other hand, was measured using two different assays and is reported in different units. To account for these different units a scaling parameter was introduced (assuming a linear relation between U/L and ng/mL, see below).

Table 1

Comparison of demographic characteristics of the placebo arms from the tibolone and EPIC study. The EPIC study was split up in two groups (EPIC 1 and EPIC 2) based on the years since menopause (YSM)

		Placebo arm of EPIC study
Characteristic	Placebo arm of tibolone study	EPIC 1 1<YSM<5	EPIC 2 All YSM
Subjects, n	149	222	470
Mean age at baseline ± SD, years	53.0 ± 3.3	50.3 ± 3.3a	53.3 ± 3.7c
Mean BMI at baseline ± SD, kg/m2	25.3 ± 3.8	25.2 ± 3.5a	25.2 ± 3.6a
Mean time since menopause at baseline ± SD, years	2.3 ± 0.9	2.5 ± 1.4c	5.7 ± 5.4b
Mean LS-BMD at baseline ± SD, g/cm2	1.04 ± 0.16	0.93 ± 0.12b	0.94 ± 0.12b
Mean TH-BMD at baseline ± SD, g/cm2	0.89 ± 0.13	0.84 ± 0.11d	0.85 ± 0.12d
Mean BSAP at baseline ± SD (U/Le and ng/mlf)	106.7 ± 29.6	11.2 ± 4.4b	11.1 ± 4.4b
Mean NTX at baseline ± SD, nmol bce/mmol cr	50.1 ± 23.0	67.1 ± 38.1b	88.0 ± 45.0b

Not significant when either of the two EPIC datasets is compared to tibolone study, using two-sample Student's t-test assuming equal variance.

P < 0.001 when either of the two EPIC datasets is compared to tibolone study using Welch's t-test.

P < 0.05 when compared to tibolone study using Welch's t-test.

P < 0.001 when either of the two EPIC datasets is compared to tibolone study using two-sample Student's t-test assuming equal variance.

Determination of BSAP activity was based on selective inhibition of the three common isoenzymes of alkaline phosphatase (bone, liver, and a third group of isoenzymes from intestinal mucosa, placenta, and neoplastic tissue that are sensitive to l-phenylalanine) and reported in U/L.22

Determination of BSAP activity was based on a solid phase, monoclonal antibody immunoenzymetric assay and is reported in ng/mL.16

Qualification of the disease systems analysis model using EPIC 1

Previous studies have indicated that most women have a rapid phase of bone loss following menopause, followed by a protracted period of slower bone loss (age-related) that continues into old age.18,19 The model takes the timescale of the disease process into account (i.e., YSM, see Methods and Figure 2).

Due to the design of the tibolone study (see Supplementary Information), only the effect on disease progression between years 1 and 4 after menopause could be analyzed. Therefore, to qualify the disease systems analysis model we first analyzed the results using data from EPIC 1 (Table 1). In order to apply the model (see Methods) to this dataset two modifications were applied: first, a scaling parameter was introduced to account for the different units for BSAP in the two clinical studies. The baseline value for BSAP was fixed at the value found in the tibolone study (in units U/L) and we only estimated the value of the scaling parameter. This scaling parameter (λBSAP,0) affects the baseline value of BSAP:

Second, the placebo effect during treatment with calcium was implemented using a separate calcium elimination rate constant (kCa) for the onset and offset of the placebo function:

This adaptation in placebo function allows for a placebo effect that wears off over time much more slowly than the original model, something we also observed during data analysis of the placebo response for the EPIC study.

Since it was our aim to qualify the model, we fixed the system-specific parameters and estimated the parameters related to the placebo function, transducer functions for bone turnover markers, and BMD as well as the interindividual (IIV) and residual variability observed for the bone turnover markers and BMD (see Table 2). All parameter values were in accordance with the original values or with the available data and could be estimated with good precision. Goodness of fit plots indicated that there were no underlying trends of model misspecifications. Overall, we found that the disease systems analysis model was able to satisfactorily predict the course of the changes for all four measurements in this external dataset.

Table 2

Comparison of the population parameter estimates for the placebo arm of the tibolone model and the models applied to of EPIC 1 and EPIC 2

Parameter (unit)	Description	Value (%CV) tibolone study	Value (%CV) EPIC 1 1≤YSM≤5	Value (%CV) EPIC 2 All YSM
System-related parameters
zs (fraction)	Constant in σ(z)	fixed at 0.659
kB (day−1)	Elimination rate constant of osteoblast	fixed at 0.0109
kestrogen (day−1)	Estrogen elimination rate constant	fixed at 0.00763
DA (day−1)	Osteoclast apoptosis rate constant	Could not be estimated, fixed at 1
b (%)	Status of the disease process at baseline	Could not be estimated, fixed at 1
Placebo-related parameters
kCa (day−1)	Calcium elimination rate constant	0.00237 (14.43)	NA	NA
kCa,onset (day−1)	Calcium elimination rate constant for onset	NA	0.00128 (7.2)	0.0009 (12.4)
KCa,offset (day−1)	Calcium elimination rate constant for offset	NA	0.000374 (10.7)	0.000226 (21.9)
Transducer function bone turnover markers
BSAP0 (U/L)	BSAP baseline value	fixed at 97.4a
λBSAP0 (−)	BSAP baseline scaling parameter	NA	−0.894 (0.4)	−0.896 (0.3)
NTX0 (nmol bce/mmol cr)	NTX baseline value	35.9 (2.58)	53.4 (2.6)	49.5 (5.7)
ρBSAP (−)	BSAP transduction parameter	fixed at 97.4a
ρNTX (−)	NTX transduction parameter	0.366 (11.61)	0.564 (7.6)	0.56 (15.6)
Transducer function bone mineral density
kLS (mg/day)	Zero-order turnover rate constant for LS	0.11 (34.45)	fixed at 0.11	NA
kin,ls (mg/day)	Zero-order production rate constant for LS-BMD	NA	NA	1.13 (22.7)
LS-BMD0 (g/cm2)	LS-BMD baseline value	0.97 (0.8)	0.98 (0.8)	0.99 (0.7)
BMI-LS-BMD0 fraction (−)	BMI fraction of LS-BMD baseline	0.00792 (14.90)	0.00892 (24.8)	0.0111 (14.5)
kTH (mg/day)	Zero-order turnover rate constant for TH	0.0821 (34.71)	fixed at 0.0821	NA
kin,th (mg/day)	Zero-order production rate constant for TH-BMD	NA	NA	0.295 (14.7)
TH-BMD0 (g/cm2)	TH-BMD baseline value	0.87 (0.8)	0.88 (0.8)	0.88 (0.7)
BMI-TH-BMD0 fraction (−)	BMI fraction of TH-BMD baseline	0.0133 (8.95)	0.0118 (22.1)	0.0154 (11.4)
ρBMD (−)	BMD transduction parameter	0.784 (33.80)	fixed at 0.784	NA
DAOB (−)	Coefficients for stimulation by relative osteoblast activity	NA	NA	0.121 (6.0)
DAOC(−)	Coefficients for stimulation by relative osteoclast activity	NA	NA	0.0456 (10.6)
Interindividual variability
IIV NTX0 (%)	IVV NTX baseline	29 (7.1)	39 (5.2)	40 (3.5)
IIV BSAP0 (%)	IVV BSAP baseline	25 (5.8)	33 (6.6)	32 (5.6)
IIV corr NTX0-BSAP0 (−)	IVV correlation NTX-BSAP baseline	0.48 (10.0)	0.50 (15.7)	0.50 (12.1)
IIV BMDLS,0 (%)	IVV LS-BMD baseline	11 (5.5)	12 (4.7)	12 (3.4)
IIV BMDTH,0 (%)	IVV TH-BMD baseline	11 (5.7)	12 (4.5)	12 (3.4)
IIV corr BMDLS,0-TH,0 (−)	IVV correlation LS-TH-BMD baseline	0.62 (7.0)	0.59 (2.3)	0.60 (1.4)
Residual variability
εBSAP (SD)	Residual variability BSAP	0.164 (2.1)	0.174 (4.4)	0.184 (4.6)
εBSAP,extremes (SD)	Residual variability BSAP extremes	0.632 (13.6)	0.536 (22.9)	0.521 (22.8)
εNTX (SD)	Residual variability NTX	0.312 (1.8)	0.307 (2.5)	0.314 (2.0)
εNTX, extremes (SD)	Residual variability NTX extremes	0.984 (8.1)	0.570 (14.4)	0.511 (19.2)
εLS (SD)	Residual variability LS-BMD	0.019 (2.2)	0.024 (3.6)	0.022 (2.7)
εTH (SD)	Residual variability LS-BMD	0.015 (3.0)	0.021 (3.5)	0.020 (2.3)

NA: Not applicable.

Fixed at the tibolone value, see main text for explanation.

Could not be estimated.

Application of disease systems analysis model to all data in the placebo arm of the EPIC study

The change in LS-BMD and TH-BMD is described by a zeroth-order process (Eq. 5), which can be viewed as a simplification of classical turnover models in which the loss term is omitted. For EPIC 1 (and the original tibolone study) this zeroth-order process turned out to be a good approximation. However, when analyzing all 470 women in the placebo arm of EPIC 2 a large deviation between model predictions and observed values was found for LS-BMD and TH-BMD, but not for the resorption (NTX) and formation (BSAP) markers. Likely this is due to the fact that a longer disease trajectory is being looked at. Therefore, we used an indirect response model for the BMD dynamics proposed by others.8,23 Changes in BMD are described with osteoblasts and osteoclasts stimulating the production and degradation processes as: where kinLS(TH) is a zeroth-order production rate constant for LS-BMD (TH-BMD), koutLS(TH) is the first-order degradation rate constant for LS-BMD (TH-BMD), and DAOB and DAOC are the coefficients for stimulation by the relative osteoblast (y) and osteoclast (z) activity. The model, including the indirect response equations for BMD, was then applied to EPIC 2. The resulting parameters values are shown in Table 2. The two models yielded similar values for baseline parameter values and comparable values for the residual variability. Parameter estimates for the indirect response equations for BMD were consistent with physiological values reported in the literature.23,24

To test the appropriateness of the current model, a visual predictive check (VPC) and numerical predictive check (NPC, see Supplementary Information) was performed. In Figure 3 the relative changes from baseline are shown on the study timescale for the degradation marker NTX, bone formation marker BSAP, LS-BMD, and TH-BMD. The results presented in Figure 3 show that the model adequately describes all biomarker data as the 5th, 50th, and 95th percentiles of the real data (red lines) overlap with the 5th, 50th, and 95th percentiles of the simulated data (black lines) and lay within the respective prediction intervals (blue and red areas). Furthermore, the 5th, 50th, and 95th percentiles of the model prediction (black lines) follow that of the real data (red lines). Overall, we conclude that the model is able to describe the dynamics of the biomarker in the entire placebo arm to a very good approximation.

Figure 3

Visual predictive check plots of the degradation marker NTX, the bone formation marker BSAP, LS-BMD, and TH-BMD on the study timescale for the model applied to all women in the placebo arm (EPIC 2). The blue dots represent the percentage change from baseline of the available observations. The 5th, 50th, and 95th percentiles of the real data in the bins are presented by the red dashed, red solid, and red dashed line, respectively. The 5th, 50th, and 95th percentiles of the simulated data (n = 500) in the bins are presented by the black dashed, black solid, and black dashed line, respectively. The confidence interval for the simulated data 5th, 50th, and 95th percentiles for each of the bins is presented by the blue, red, and blue area, respectively. Note that BSAP was measured in a random sample of 205 women and not in the entire population.

Sensitivity of the system-specific parameters

As a second level of model qualification, we tested the sensitivity of the system-specific parameters for both models (EPIC 1 and EPIC 2). To find how these parameters vary, the model was optimized with all other parameter values fixed at the values as shown in Table 2, and each system-specific parameter was estimated one at a time. Parameters values for zs and kB were lower compared to the tibolone study and could be estimated with acceptable precision. Kestrogen could not be identified for EPIC 1 and only with a very high coefficient of variation for EPIC 2 (see Table 3).

Table 3

Change in single system-specific parameter values when the disease systems analysis model was optimized with all other parameters fixed according to the values as shown in Table 2

Parameter (unit)	Value (%CV) tibolone study	Value (%CV) EPIC 1 1≤YSM≤5	Fold change difference	Value (%CV) EPIC 2 All YSM	Fold change difference
zs (−)	0.659 (24.0)	0.497 (27.6)	0.75	0.366 (19.9)	0.55
kB (day−1)	0.0109 (11.9)	0.0032 (36.5)	0.29	0.0019 (22.4)	0.17
kestrogen (day−1)	0.0076 (21.3)	(−)**	NA	0.018 (201)	2.4
Da (day−1)	1*(−)	1*(−)	NA	1*(−)	NA
b (%)	1*(−)	1*(−)	NA	1*(−)	NA

*Fixed at 1.

**Could not be estimated.

NA, not applicable.

Simulations of BMD dynamics with zeroth-order and indirect response model

Simulations with both types of BMD equations were performed in order to compare the effect on the dynamics. The parameter values as shown in Table 2 were used for these simulations. For both equations the change of BMD is determined by the balance between the relative changes in y and z (i.e., S = z/y). In Figure 4a it is shown how these relative quantities evolve as a result of the placebo treatment. Both models evolve towards a similar maximal value of S. The underlying dynamics for both models is different, however: the zeroth-order model has reached this maximal value already after 2 years, whereas the indirect response model reached this value only after 4 years. Note that S represents the ratio of relative osteoclast activity over relative osteoblast activity. Therefore, there is net bone loss if this ratio is bigger than 1 and net bone formation when it is smaller than 1. This can also be seen in Figure 4b, where the change of LS-BMD (gray lines) and TH-BMD (black lines) over time is plotted and decreases as a result of disease progression. The difference in the two BMD equations can be seen from this plot: the zeroth-order response function (solid lines) has a concave shape, whereas the indirect response model (dashed lines) has a convex shape. This implies that right after onset of disease progression (e.g., beginning of the menopause) the change in BMD is largest. In addition, this change is bigger for LS-BMD compared to TH-BMD. These model predictions were compared to the available observations. So, for every individual we calculated the difference between baseline value and the latest available measurement for LS- and TH-BMD and plotted this difference vs. YSM. As shown in Figure 4c,d, we also found the biggest change (i.e., decrease) in BMD at the LS and TH shortly after onset of menopause, and the absolute mean change was ∼25% higher in BMD at the lumbar spine compared to the total hip.

Figure 4

Comparison of zeroth-order and indirect response model for BMD dynamics. (a) Changes in the net bone cell activity (S = z/y) vs. time for the zeroth-order process (ZO, solid line) and the indirect response model (IR, dashed line). The arrows indicate the start and end of the placebo treatment in the two studies. (b). Changes in LS-BMD (gray lines) and TH-BMD (black lines) for the zeroth-order process (solid line) and the indirect response model (dashed line). Scatterplot showing the change between baseline and the latest observation available for LS-BMD (c) and TH-BMD (d) vs. years since menopause.

DISCUSSION

In this study we applied a recently proposed disease systems analysis model for postmenopausal osteoporosis to data from another population that was used for model development and parameter fitting. The model is based on a mechanistic model describing osteoblast and osteoclast activity13 that has been reduced in order to apply it in a population approach. It was shown that the model reduction did not jeopardize the dynamic properties of the model.14 The reduced model was successfully applied to describe responses in treatment with various doses of tibolone and/or calcium in postmenopausal women.15 The current systems pharmacology approach allowed us to test whether the original values of the system-specific parameters allowed for an adequate description of disease progression in a new (external) population. To this aim, we used the placebo arm of the EPIC study.16 This led to a total of 470 women who received placebo treatment for a period of 4 years. Although others have looked at the modeling of placebo response,25,26 to our knowledge this is the first attempt to use a systems-pharmacology population approach for the description of osteoporosis disease progression in response to placebo treatment in an external population. Robust descriptions of the placebo response and disease progression are crucial for the quantification of treatment effects. The placebo effect, PCA (Eq. 3), was updated with an onset and offset parameter and was incorporated into the model in a mechanism-based manner: as an inhibitory factor of the RANK receptor occupancy as calcium promotes its effects on the RANK- RANKL-OPG system through PTH.13,15

It should be noted that all subjects in the tibolone study received 500 mg calcium carbonate once daily. In the EPIC study, dietary calcium intake from food sources was monitored using validated dietary questionnaires. Women with a calcium intake of less than 500 mg per day were advised to increase their intake (either by diet or supplements). In principle, all women in the EPIC study should thus also have received at least 500 mg calcium per day. Since this was monitored at center visits only via questionnaires, more variation in calcium intake could be expected for this study. In order to characterize the true calcium effect, data obtained following the administration of different doses of calcium would be needed.

The disease systems analysis model is a mathematical reduced version14 of the cell interaction model as proposed by Lemaire et al.13 The model reduction did not influence the model dynamics and bears the advantage that all remaining parameter values are identifiable despite the fact that the biomarkers available in this study respond at different timescales upon treatment with placebo. With all system-specific parameter values fixed at their published values, we showed that the model is able to adequately describe the dynamics of the biomarkers within the external population. Furthermore, a sensitivity analysis revealed that, with the exception of kestrogen, all systemic parameters could be reestimated with the current data and their values were in the same range as the previously reported values. Possibly, the heterogeneity in YSM in the population resulted in the large uncertainty of the kestrogen parameter. We anticipate that on the basis of individual serum estrogen levels a better description of the underlying disease dynamics will be feasible.

The most notable difference between the two clinical studies is the YSM at baseline. In the model the YSM are incorporated as the timescale of disease progression (Figure 2). Inclusion of the entire range of YSM in our analysis required the adaptation of the equation for the BMD dynamics. The indirect response model23 used in this study can be seen as a more extensive form of the (originally proposed) zeroth-order process. With this small adaptation we were able to describe the dynamics in all biomarkers for the entire range of YSM with good precision. Furthermore, the fact that the model also accurately describes disease progression far outside the timeframe compared to the timeframe that was used originally in the tibolone study provides evidence for the values of the system-specific parameters. These results provide a strong basis for the potential to use this disease systems analysis model to evaluate the (comparative) effect of other drugs with different mechanisms (e.g., bisphosphonates, unpublished data), or link the mechanism-based core of this model to different biomarkers.

The core model is defined to describe the dynamics as a relative change related to the information in the dynamics of the markers linked to it (based on their relationship to the system). We believe, and this should be tested, that all types of markers can be linked and inform such a system. Whether that holds remains to be seen based on available data. In contrast to the tibolone study, we did not have measurements on the biomarker osteocalcin available. Osteocalcin is a biomarker for the combined osteoblast and osteoclast activity. Osteocalcin is built into bone—when bone is degraded it will become available again in the systemic circulation so the measurement is a composite of building and degradation of bone. We found that removal of osteocalcin from the model did not alter its performance. This also shows the strength of a disease systems analysis model; the core remains the same and the markers supplying the required information on the system can change, which is in contrast to empirical models where (single) marker(s) are directly linked. On the other hand, it remains to be seen how this removal will affect future applications of the model when other drug treatments are considered.

In conclusion, we have shown that a recently proposed disease systems analysis model for osteoporosis can be applied to describe the effect of the placebo treatment on disease progression in a new population. The mechanism-based model allowed for 1) an adequate description of the available biomarkers, 2) allowed for the inclusion of the years since menopause (the characteristic timescale in the model) that was more than five times longer than in the original model. These findings can be considered a first step towards model validation and qualification of this osteoporosis model, which will be further elaborated by also applying the model to other therapeutic interventions.