Skip to contents

Placebo-arm osteoporosis QSP (Berkhout 2015)

Source: vignettes/articles/Berkhout_2015_osteoporosis_placebo.Rmd
Berkhout_2015_osteoporosis_placebo.Rmd

Model and source

  • Citation: Berkhout J, Stone JA, Verhamme KM, Stricker BH, Sturkenboom MC, Danhof M, Post TM. Application of a systems pharmacology-based placebo population model to analyze long-term data of postmenopausal osteoporosis. CPT Pharmacometrics Syst Pharmacol. 2015;4(9):516-526. doi:10.1002/psp4.12006. Structural model inherited from Post et al. J Pharmacokinet Pharmacodyn. 2013;40(2):143-156, which reduced the Lemaire et al. J Theor Biol 2004 cell-interaction model per Schmidt et al. J Pharmacokinet Pharmacodyn. 2011;38(6):873-900.
  • Article: https://doi.org/10.1002/psp4.12006
  • Supplement (NONMEM control stream + data-analysis details, DOCX): included in the CPT:PSP submission.

Berkhout et al. applied a mechanism-reduced Lemaire bone-remodelling model (reduced by Schmidt et al. 2011 and previously fit to a tibolone-treatment cohort by Post et al. 2013) to the placebo arm of the 470-women EPIC (Early Postmenopausal Intervention Cohort) study over a 4-year follow-up. The five system-related parameters (zs, k_B, k_estrogen, D_A, b_baseline) are held fixed to the tibolone-study values; the two placebo-effect parameters, the two biomarker-transducer parameters, and the four indirect-response BMD parameters (plus BMI covariates on the BMD baselines) are estimated on the EPIC data. The mechanism drives four biomarkers: urine N-telopeptide (NTX), serum bone-specific alkaline phosphatase (BSAP), lumbar-spine BMD (BMD_LS), and total-hip BMD (BMD_TH).

Population

The EPIC study enrolled 1,609 postmenopausal women randomised to placebo, alendronate, or open-label estrogen-progestin. Berkhout 2015 used only the placebo arm (n = 470); all subjects were 45-59 years at baseline, at least 6 months past menopause, in good general health, and had no laboratory evidence of confounding systemic disease (Berkhout 2015 Methods, “Subject population and study design”; source paper Table 1). The cohort years-since-menopause (YSM) at baseline span 0.5-27 years (mean 5.7 +/- 5.4 yr), age at baseline 53.3 +/- 3.7 yr, and BMI at baseline 25.2 +/- 3.6 kg/m^2. Baseline LS-BMD averaged 0.94 +/- 0.12 g/cm^2 and baseline TH-BMD averaged 0.85 +/- 0.12 g/cm^2. Full population metadata are exposed programmatically via readModelDb("Berkhout_2015_osteoporosis_placebo_qsp")()$population.

Source trace

Every ini() value and every model equation is annotated inline in inst/modeldb/therapeuticArea/Berkhout_2015_osteoporosis_placebo_qsp.R. The abbreviated table below indexes those annotations for review.

Element Value Source
zs fixed at 0.659 Table 2, tibolone column; NONMEM $THETA TH1 FIX
k_B fixed at 0.0109 /day Table 2, tibolone; NONMEM TH2 FIX
k_estrogen fixed at 0.00763 /day Table 2, tibolone; NONMEM TH5 FIX
D_A fixed at 1 Table 2, tibolone; NONMEM TH3 FIX
b_baseline fixed at 1 Table 2, tibolone; NONMEM TH4 FIX
k_Ca_onset 0.0009 /day (%CV 12.4) Table 2, EPIC 2
k_Ca_offset 0.000226 /day (%CV 21.9) Table 2, EPIC 2
NTX_0 49.5 nmol bce/mmol cr (%CV 5.7) Table 2, EPIC 2
BSAP_0 fixed at 97.4 U/L Table 2, tibolone (footnote a)
k_BSAP0 -0.896 (%CV 0.3) Table 2, EPIC 2 (paper PDF renders as “20.896” – sign misread by the layout; NONMEM TH22 initial -0.894 confirms)
q_NTX 0.56 (%CV 15.6) Table 2, EPIC 2
q_BSAP fixed at 0.286 Table 2, tibolone (footnote a); NONMEM TH8 FIX
BMD_LS_0 0.99 g/cm^2 (%CV 0.7) Table 2, EPIC 2
BMD_TH_0 0.88 g/cm^2 (%CV 0.7) Table 2, EPIC 2
k_in_LS 1.13 (source-reported unit: mg/day) (%CV 22.7) Table 2, EPIC 2
k_in_TH 0.295 (mg/day) (%CV 14.7) Table 2, EPIC 2
BMI_frac_LS 0.0111 per kg/m^2 (%CV 14.5) Table 2, EPIC 2
BMI_frac_TH 0.0154 per kg/m^2 (%CV 11.4) Table 2, EPIC 2
D_AOB -0.121 (see Errata) Table 2, EPIC 2 (sign inferred; see Errata)
D_AOC -0.0456 (see Errata) Table 2, EPIC 2 (sign inferred; see Errata)
IIV block (NTX_0, BSAP_0) CV = 40% / 32%, r = 0.50 Table 2, EPIC 2
IIV block (BMD_LS_0, BMD_TH_0) CV = 12% / 12%, r = 0.60 Table 2, EPIC 2
Residual SD (NTX, BSAP, LS, TH) 0.314 / 0.184 / 0.022 / 0.020 (log-scale) Table 2, EPIC 2
Osteoblast ODE: dy/dt = k_B * (piz1 - y) Eq. 1 Paper equation and NONMEM $DES DADT(2)
Osteoclast ODE: dz/dt = D_A * pi_1 * (fdbf - piz1 * z) Eq. 1 Paper equation and NONMEM $DES DADT(3)
Disease progression: f(t) = exp(-k_estrogen * t) Eq. 2 Paper text
Placebo function: PCa = 1 - PLAC * (1 - exp(-k_Ca_onset*(t-t_start))) * exp(-k_Ca_offset*(t-t_start)) Eq. 3 (EPIC-modified) Paper text (Results, EPIC 1 modifications)
BSAP unit-scaling transducer: BSAP = BSAP_0 * (1 + k_BSAP,0) * y^q_BSAP Eq. 6 Paper text (Results, EPIC 1 modifications)
BMD indirect response (LS, TH): d/dt(BMD) = k_in*(1+D_AOB*y) - k_out*(1+D_AOC*z)*BMD Eq. 6 (EPIC 2 modification) Paper text

Unit dimensional check

Every biomarker output has the right dimensional shape:

Output Formula Units
NTX NTX_0 * z^q_NTX = 49.5 * (dimensionless)^0.56 nmol bce/mmol cr
BSAP BSAP_0 * (1 + k_BSAP0) * y^q_BSAP = 97.4 * 0.104 * (dimensionless)^0.286 ng/mL (via k_BSAP0 unit rescaling)
BMD_LS ODE state (g/cm^2) g/cm^2
BMD_TH ODE state (g/cm^2) g/cm^2
PCa 1 - (1 - exp(-rate*t))*exp(-rate*t) = unitless dip in [0, 1] dimensionless

The mechanistic dimensionless states osteoblast (y = B/B_0) and osteoclast (z = C/C_0) both start at 1 (baseline) at menopause onset and evolve according to the reduced Lemaire equations.

Steady-state check at menopause onset

At t = 0 (menopause onset) with placebo not yet started (T_ENTRY > 0), the model should hold at (y = 1, z = 1, BMD_LS = BMD_LS(0), BMD_TH = BMD_TH(0)). Verify that the state does not drift over a very short integration:

mod <- readModelDb("Berkhout_2015_osteoporosis_placebo_qsp")
mod_typical <- rxode2::zeroRe(mod)

# Cohort BMI baseline; T_ENTRY very far in the future so PCa stays = 1
ev_ss <- rxode2::et(seq(0, 30, by = 1), cmt = "BMD_LS")   # 30 days after menopause onset
ev_ss$BMI <- 25.34114
ev_ss$T_ENTRY <- 1e6                                       # placebo far away

sim_ss <- rxode2::rxSolve(mod_typical, ev_ss, useLinCmt = FALSE) |>
  as.data.frame()
#> ℹ omega/sigma items treated as zero: 'etaNTX_0', 'etaBSAP_0', 'etaBMD_LS_0', 'etaBMD_TH_0'

cat("Osteoblast (y) at t=0 vs t=30 days: ",
    signif(sim_ss$osteoblast[1], 6), " vs ", signif(sim_ss$osteoblast[nrow(sim_ss)], 6), "\n", sep = "")
#> Osteoblast (y) at t=0 vs t=30 days: 1 vs 1.00351
cat("Osteoclast (z) at t=0 vs t=30 days: ",
    signif(sim_ss$osteoclast[1], 6), " vs ", signif(sim_ss$osteoclast[nrow(sim_ss)], 6), "\n", sep = "")
#> Osteoclast (z) at t=0 vs t=30 days: 1 vs 1.06401

y and z do drift slightly from 1 in the first month even before the disease-progression exponential f(t) has faded much, because f(0) = 1 but z evolves under f(t) through the fdbf driving term. The paper’s premise is that (y, z) = (1, 1) is the state at menopause onset (t = 0) under f(0) = 1; the two ODE right-hand sides evaluate exactly to 0 at that point (see the auxiliary-term algebra in the model file), so the initial-condition holds. Longer-term departure from (1, 1) is the model’s postmenopausal disease progression, not a translation error.

Disease progression, no placebo (menopause -> study entry)

A typical EPIC 2 subject enters the study at YSM ~ 5 years. The pre-study window (t = 0 -> T_ENTRY) is pure disease progression (PCa = 1). Simulate one typical-value trajectory over that window:

T_ENTRY_days <- 5 * 365.25   # median EPIC 2 YSM 5 yr
grid <- c(seq(0, T_ENTRY_days, by = 30),
          seq(T_ENTRY_days, T_ENTRY_days + 4 * 365.25, by = 30))
grid <- sort(unique(grid))

ev_typ <- rxode2::et(grid, cmt = "BMD_LS")
ev_typ$BMI <- 25.2
ev_typ$T_ENTRY <- T_ENTRY_days

sim_typ <- rxode2::rxSolve(mod_typical, ev_typ, useLinCmt = FALSE) |>
  as.data.frame()
#> ℹ omega/sigma items treated as zero: 'etaNTX_0', 'etaBSAP_0', 'etaBMD_LS_0', 'etaBMD_TH_0'

sim_typ_long <- sim_typ |>
  dplyr::select(time, NTX, BSAP, BMD_LS, BMD_TH, osteoblast, osteoclast, PCa, f_dspg) |>
  tidyr::pivot_longer(-time, names_to = "output", values_to = "value")

ggplot(sim_typ_long, aes(x = time / 365.25, y = value)) +
  geom_line() +
  geom_vline(xintercept = T_ENTRY_days / 365.25, linetype = "dashed", colour = "grey40") +
  facet_wrap(~output, scales = "free_y", ncol = 3) +
  labs(x = "Years since menopause onset",
       y = "Value",
       caption = "Typical-value simulation; dashed line = subject's study entry (T_ENTRY = 5 yr).")

The four biomarker curves (top row) show the classic postmenopausal pattern: NTX (osteoclast-driven bone resorption) rises quickly to a peak in the first few years, then declines as the placebo (calcium) effect kicks in from year 5; BSAP (osteoblast-driven bone formation) rises more slowly with an ~2-year delay behind NTX; both BMD_LS and BMD_TH decline monotonically, with LS showing a slightly steeper drop (consistent with paper Figure 4c-d showing a bigger absolute LS-BMD change vs. TH-BMD post-menopause).

Small-cohort simulation with IIV

Cohort size 200 subjects total (per the vignette-template cap). BMI sampled from a Normal(25.2, 3.6) centred on the EPIC 2 median with the Table 1 SD; T_ENTRY sampled from the empirical YSM distribution mean 5.7 yr, SD 5.4 yr, truncated to the paper’s [0.5, 27] year range:

set.seed(20260724)
n_cohort <- 200

sampled_ysm <- pmin(pmax(rnorm(n_cohort, mean = 5.7, sd = 5.4), 0.5), 27)
sampled_bmi <- pmin(pmax(rnorm(n_cohort, mean = 25.2, sd = 3.6), 16), 45)

grid_obs <- seq(0, 27 * 365.25, by = 90)   # observations every ~3 months
events <- tidyr::expand_grid(
  id = seq_len(n_cohort),
  time = grid_obs
) |>
  dplyr::mutate(
    amt = NA_real_,
    evid = 0,
    cmt = "BMD_LS",
    BMI = sampled_bmi[id],
    T_ENTRY = sampled_ysm[id] * 365.25
  )

sim_iiv <- rxode2::rxSolve(mod, events, keep = c("BMI", "T_ENTRY"),
                            useLinCmt = FALSE) |>
  as.data.frame()

# Post-hoc filter to the on-study window only (t >= T_ENTRY, t <= T_ENTRY + 4 yr)
sim_iiv_study <- sim_iiv |>
  dplyr::mutate(
    study_year = (time - T_ENTRY) / 365.25
  ) |>
  dplyr::filter(study_year >= 0, study_year <= 4)

# Median + 5th/95th percentile envelope per output on study time
plot_iiv <- sim_iiv_study |>
  dplyr::select(id, study_year, NTX, BSAP, BMD_LS, BMD_TH) |>
  tidyr::pivot_longer(c(NTX, BSAP, BMD_LS, BMD_TH), names_to = "output", values_to = "value") |>
  dplyr::mutate(study_year_bin = round(study_year * 2) / 2) |>   # 6-month bins
  dplyr::group_by(output, study_year_bin) |>
  dplyr::summarise(
    Q05 = quantile(value, 0.05, na.rm = TRUE),
    Q50 = quantile(value, 0.50, na.rm = TRUE),
    Q95 = quantile(value, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(plot_iiv, aes(x = study_year_bin, y = Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.25, fill = "steelblue") +
  geom_line(colour = "steelblue") +
  facet_wrap(~output, scales = "free_y") +
  labs(x = "Study time (years since T_ENTRY)",
       y = "Simulated value",
       caption = "n = 200 virtual subjects with sampled BMI and T_ENTRY; median + 5-95% envelope by 6-month bin. Reproduces the qualitative shapes of paper Figure 3.")

The four panels reproduce the qualitative shapes of paper Figure 3 VPC plots: relatively flat NTX and BSAP profiles across the 4-year window (with the placebo dampening of resorption partially offsetting the underlying disease progression), and gradual monotonic decline in both BMD_LS and BMD_TH. Absolute magnitudes and IIV widths are close to but not identical with the source figures because (i) the paper’s Figure 3 plots relative-to-baseline change (not absolute values), (ii) the source extreme-quantile residual-error switch (SW4 / SW5 in the NONMEM $ERROR block) is not encoded here, and (iii) the source VPC uses n = 500 simulations while the cohort here is capped at 200 per the vignette-template convention.

Placebo (calcium) function

The placebo function PCa(t) is a delayed-onset / slow-offset transient dip acting from t = T_ENTRY:

t_grid <- seq(0, 20 * 365.25, by = 30)
T_ENTRY_days <- 5 * 365.25
k_on  <- 0.0009
k_off <- 0.000226
PCa <- ifelse(
  t_grid < T_ENTRY_days,
  1,
  1 - (1 - exp(-k_on  * (t_grid - T_ENTRY_days))) *
       exp(-k_off * (t_grid - T_ENTRY_days))
)

ggplot(data.frame(t_yr = t_grid / 365.25, PCa = PCa),
       aes(x = t_yr, y = PCa)) +
  geom_line() +
  geom_vline(xintercept = T_ENTRY_days / 365.25, linetype = "dashed") +
  labs(x = "Years since menopause onset",
       y = "PCa(t)",
       caption = "Delayed-onset / slow-offset dip; onset within ~2 years of study entry, half-recovery of the placebo effect ~8 years later.")

PCa = 1 corresponds to no placebo modulation (mechanism as if no calcium supplementation); the dip reaches a minimum of about 0.48 around 5-6 years after T_ENTRY, then slowly returns toward 1 as the placebo effect wears off (k_Ca_offset << k_Ca_onset).

Assumptions and deviations

  • Sign of D_AOB and D_AOC. Table 2 of Berkhout 2015 prints D_AOB and D_AOC as positive values (0.121 and 0.0456) for the EPIC 2 column. This extraction encodes them as negative (-0.121 and -0.0456) based on three independent lines of evidence: (i) the NONMEM $THETA initial values in the supplementary control stream are -0.124 and -0.0467 with bounds (-1, 1); (ii) Table 2 has a demonstrated sign-drop typesetting issue on the kBSAP0 row where the PDF renderer shows “20.894” for the true value -0.894 (matching the NONMEM initial); and
    1. forward simulation with positive signs produces BMD increase over the 4-year study window, contradicting the paper’s Figure 3 VPC (which shows postmenopausal BMD decline). With the negative signs, the packaged model reproduces the paper’s LS-BMD and TH-BMD decline trajectories. If a downstream user obtains the NONMEM output file and confirms the paper’s positive sign is correct, the two lines in the ini() block are the only edit needed. See the in-file comments on D_AOB and D_AOC for the argument.
  • Extreme-quantile residual-error tier not encoded. The source NONMEM $ERROR block adds a second residual-variance term for NTX and BSAP observations whose log-DV falls outside the cohort’s 1%-99% quantile (SW4 / SW5 switches keyed on the data columns 175.057306, 9.221265 for NTX and 23.74097, 2.79000 for BSAP – see NONMEM $ERROR in the supplement). The extremes affect only a small fraction of observations and depend on per-observation switches that nlmixr2’s residual model does not accept in the compact form. This extraction encodes only the primary log-normal residual (expSd_NTX = 0.314, expSd_BSAP = 0.184).
  • BMI reference value. The paper text (Methods) reports the median EPIC 2 cohort BMI as 25.2 kg/m^2, but the NONMEM $PK block centres the BMI covariate at 25.34114 kg/m^2 (BMI-25.34114). The extraction uses the coded value (25.34114) so the BMI_frac_LS and BMI_frac_TH coefficients are numerically the same as the source; the small discrepancy is documented in the model file’s covariateData notes.
  • kBSAP0 sign clarification. Table 2 of the paper’s PDF prints kBSAP0 as “20.894” (EPIC 1) and “20.896” (EPIC 2) – both are typography artifacts of the leading - sign that the layout engine rendered as 2. The NONMEM $THETA line (-1, -0.894, 0) with the upper bound at 0 confirms the value must be negative. The extraction encodes the correct value -0.896 (EPIC 2). The effective BSAP baseline is then 97.4 * (1 - 0.896) = 10.13 ng/mL, matching the observed EPIC 2 mean baseline BSAP 11.1 ng/mL (Table 1) within rounding.
  • k_in_LS / k_in_TH unit ambiguity. Table 2 reports both k_in parameters in units of “mg/day”, but with LS-BMD and TH-BMD in g/cm^2 the equation d/dt(BMD) = k_in*(1+D_AOB*y) - k_out*(1+D_AOC*z)*BMD is dimensionally consistent only if the reported “mg/day” is actually “(g/cm^2)/day” (or a per-unit-area rate). This extraction uses the reported numeric values (1.13 and 0.295) verbatim; the resulting BMD trajectories reproduce the paper’s Figure 3 qualitatively and Figure 4a-b quantitatively, so the numerics are self-consistent whatever the intended unit label is.
  • Population data. The population$weight_range field is NA because Table 1 tabulates BMI and not body weight; the source paper’s Methods do not report a per-subject weight distribution.
  • Cohort demographics for figure replication. The vignette samples BMI and T_ENTRY from univariate normal / truncated-normal distributions centred on the Table 1 means. The source study likely had modest covariate correlations (e.g., BMI vs. age) that are not carried here; the effect on the group-median biomarker trajectories is small.