Skip to contents

Methoxy polyethylene glycol-epoetin beta / C.E.R.A. (Chanu 2010)

Source: vignettes/articles/Chanu_2010_methoxyPolyethyleneGlycolEpoetinBeta.Rmd
Chanu_2010_methoxyPolyethyleneGlycolEpoetinBeta.Rmd

Model and source

  • Citation: Chanu P, Gieschke R, Charoin JE, Pannier A, Reigner B. Population pharmacokinetic/pharmacodynamic model for C.E.R.A. in both ESA-naive and ESA-treated chronic kidney disease patients with renal anemia. Journal of Clinical Pharmacology. 2010;50(5):507-520. doi:10.1177/0091270009343931
  • Description: One-compartment population PK with first-order subcutaneous absorption for methoxy polyethylene glycol-epoetin beta (C.E.R.A.), coupled to a red-blood-cell life-span hemoglobin PD model (Chanu 2010). Both IV and SC routes are supported. Body weight modifies clearance and volume and age modifies volume through normalized power models. The PD layer is the Krzyzanski life-span indirect-response model implemented as a delay differential equation: hemoglobin production is stimulated by an Emax function of serum C.E.R.A. concentration, and the same production term delayed by one red-blood-cell life span acts as the loss term. An additional SESA term, active only during the first life span after the switch, returns previously ESA-treated patients from their switch hemoglobin to their pre-ESA baseline. C-reactive protein and the prior weekly epoetin dose modify SC50 through normalized power models. Proportional-plus-additive residual error on concentration; additive residual error on hemoglobin.
  • Article: https://doi.org/10.1177/0091270009343931

C.E.R.A. (continuous erythropoietin receptor activator) is the development name for methoxy polyethylene glycol-epoetin beta, a pegylated erythropoiesis- stimulating agent (ESA) given intravenously every 2 weeks or subcutaneously every 4 weeks to correct or maintain haemoglobin in chronic kidney disease.

Chanu 2010 fits one PK model and one PD model sequentially: individual post hoc PK parameters drive the haemoglobin model. Because the PD layer consumes the PK concentration time course, the two are packaged here as a single coupled model with two endpoints (Cc in ng/mL and hb in g/dL) rather than two files.

The PD layer is the Krzyzanski red-blood-cell life-span indirect-response model. Its defining feature is that the loss term is the production term delayed by one red-blood-cell life span, which makes it a delay differential equation rather than an ordinary one:

Hb′(t)=S(t)−S(t−LS)+SESA,S(t)=Hb0LS(1+E(C(t))),E(C)=SmaxCSC50+C \mathrm{Hb}'(t) = S(t) - S(t - \mathrm{LS}) + S_{\mathrm{ESA}},\qquad S(t) = \frac{\mathrm{Hb}_0}{\mathrm{LS}}\bigl(1 + E(C(t))\bigr),\qquad E(C) = \frac{S_{\max} C}{SC_{50} + C}

SESA=Hb0−HbswLSfor 0≤t≤LS,0 thereafter S_{\mathrm{ESA}} = \frac{\mathrm{Hb}_0 - \mathrm{Hb}_{sw}}{\mathrm{LS}} \ \ \text{for } 0 \le t \le \mathrm{LS},\ \ 0 \text{ thereafter}

rxode2’s delay(state, T) supplies C(t−LS)C(t-\mathrm{LS}) from the solver’s dense output; before t=0t = 0 the constant initial condition central(0) = 0 is used, which is exactly the paper’s assumption that baseline production Hb0/LS\mathrm{Hb}_0/\mathrm{LS} “was considered to be constant during one LS before time of the Hb0 assessment”. Delay models are solved on a dense path, so rxSolve() falls back to dop853; the warning that announces this is expected.

Population 

pop <- readModelDb("Chanu_2010_methoxyPolyethyleneGlycolEpoetinBeta")()$population
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_cl_1, etaiov_cl_2, etaiov_cl_3
#> as a work-around try putting the mu-referenced expression on a simple line
str(pop, max.level = 1, give.attr = FALSE)
#> List of 13
#>  $ species       : chr "human"
#>  $ n_subjects    : int 400
#>  $ n_studies     : int 3
#>  $ age_range     : chr "18-89 years (AMICUS 18-89, MAXIMA 23-87, PROTOS 21-85; Chanu 2010 Table I)"
#>  $ age_median    : chr "54 years (AMICUS), 57 years (MAXIMA), 68 years (PROTOS); pooled median not printed"
#>  $ weight_range  : chr "36-131 kg (AMICUS 43-106, MAXIMA 36-131, PROTOS 43-115; Chanu 2010 Table I)"
#>  $ weight_median : chr "66 kg (AMICUS), 70 kg (MAXIMA), 70 kg (PROTOS); pooled median not printed"
#>  $ sex_female_pct: num 40.5
#>  $ race_ethnicity: Named num [1:4] 76.5 17.8 4.5 1.3
#>  $ disease_state : chr "Chronic kidney disease with renal anemia, on dialysis (hemodialysis 386 of 400, peritoneal dialysis 14 of 400)."| __truncated__
#>  $ dose_range    : chr "C.E.R.A. intravenously or subcutaneously, once every 2 weeks (Q2W; n = 263) or once every 4 weeks (Q4W; n = 137"| __truncated__
#>  $ regions       : chr "Not stated; three multicenter phase III studies (AMICUS, MAXIMA, PROTOS)"
#>  $ notes         : chr "PK analysis population 400 patients with 4554 measurable serum concentrations; PD analysis population 400 patie"| __truncated__

The model was built on 400 patients from 3 phase III studies of C.E.R.A. in chronic kidney disease with renal anaemia, all on dialysis (Chanu 2010 Tables I and II):

  • AMICUS (n = 135) enrolled ESA-naive patients and studied correction of anaemia: IV C.E.R.A. Q2W, titrated to Hb >= 11 g/dL with a rise of >= 1 g/dL from baseline. Median age 54 y, median weight 66 kg.
  • MAXIMA (n = 122) and PROTOS (n = 143) enrolled patients already on another ESA and studied maintenance: IV (MAXIMA) or SC (PROTOS) C.E.R.A. Q2W or Q4W, titrated to hold Hb within +/- 1 g/dL of baseline and between 10 and 13 g/dL. Median ages 57 and 68 y, median weight 70 kg in both.

Pooled: 40.5% female; 76.5% White, 17.8% Black, 4.5% Asian, 1.3% other; 96.5% on haemodialysis. 4554 measurable serum concentrations (LOQ 0.15 ng/mL) and 10 089 evaluable haemoglobin assessments.

Source trace

Every ini() entry in inst/modeldb/specificDrugs/Chanu_2010_methoxyPolyethyleneGlycolEpoetinBeta.R carries an in-file comment naming its source row. Collected here for review:

Equation / parameter Value Source location
d/dt(depot), d/dt(central), f(depot) n/a Methods, “Structural PK Model (Stage 1)” (closed-form SC and IV solutions; the equivalent ODE system is used here)
lcl 0.749 L/day Table III, Fixed effects, CL (RSE 2%)
lvc 4.72 L Table III, Fixed effects, V (RSE 2%)
lka 0.825 1/day, fixed Table III, Fixed effects, ka (footnote a “Fixed”)
lfdepot 0.394 Table III, Fixed effects, F (RSE 4%)
e_wt_cl 0.571 Table III, Covariate effects, “Effect of BW on CL” (RSE 13%)
e_wt_vc 0.443 Table III, Covariate effects, “Effect of BW on V” (RSE 17%)
e_age_vc 0.267 Table III, Covariate effects, “Effect of age on V” (RSE 19%)
covariate form TVP = theta_P * (COV/rCOV)^theta_COV, combined multiplicatively n/a Methods, “Covariate Model (Stage 2)”
etalcl 0.0784 (CV% 28) Table III, Random effects IIV, CL (RSE 9%)
etalvc 0.0729 (CV% 27) Table III, Random effects IIV, V (RSE 11%)
etalka 0.6724 (CV% 82), fixed Table III, Random effects IIV, ka (footnote a “Fixed”)
IIV on F CV% 0, fixed (no eta encoded) Table III, Random effects IIV, F (footnote a “Fixed”)
etaiov_cl_1..3 0.0081 (CV% 9) Table III, Random effects IOV, CL (RSE 32%); occasion definition from Methods, “Interoccasion Variability Model”, occasion count from Table II
propSd, addSd from m = 0.150 ng/mL (fixed), sigma1^2 = 0.141, sigma2^2 = 0.691 Table III, Error model; functional form from Results, “Population PK Analysis” (double-exponential error model)
Hb'(t), S(t), S_ESA, E(C(t)) n/a Methods, “Structural PD Model (Stage 1)”
lemax (Smax) 0.425 Table IV, Fixed effects, Smax (RSE 13%)
lec50 (SC50) 0.898 ng/mL Table IV, Fixed effects, SC50 (RSE 34%)
lmtt_rbc (LS) 61.3 day, fixed Table IV, Fixed effects, LS (footnote a “Fixed”)
lrbase (Hb0) 9.30 g/dL, fixed Table IV, Fixed effects, Hb0 (footnote a “Fixed”)
e_crp_ec50 0.319 Table IV, Covariate effects, “Effect of CRP on SC50” (RSE 52%)
e_esad_ec50 0.303 Table IV, Covariate effects, “Effect of DEPO on SC50” (RSE 36%)
etalemax 2.0164 (CV% 142) Table IV, Random effects IIV, Smax (RSE 40%)
etalec50 31.2481 (CV% 559) Table IV, Random effects IIV, SC50 (RSE 34%)
etalmtt_rbc 0.1024 (CV% 32) Table IV, Random effects IIV, LS (RSE 45%)
etalrbase 0.0625 (CV% 25) Table IV, Random effects IIV, Hb0 (RSE 53%)
addSd_hb sqrt(0.357) = 0.5975 g/dL Table IV, Error model, sigma^2 (RSE 4%)
reference WT 70 kg, AGE 60 y, CRP 5 mg/L, ESAD 7000 IU/week n/a not printed; derived from Table I medians – see “Assumptions and deviations”

Event-table helper

The model declares two endpoints, so every observation row must say which endpoint it belongs to. The form that works across multi-endpoint rxUi models is dvid alone on observation rows with cmt = NA_character_, and the real ODE state on dose rows. rxSolve() returns both observables as columns on every row, so a single dvid value is enough to recover the full output frame. useLinCmt = FALSE is mandatory: the default ODE-to-linCmt() auto-conversion corrupts the dvid mapping (and would also discard the delay term).

OCC follows the paper’s design: PK sampling followed the doses given on study days 1, 57 and 127/141 (Table II), so the three IOV occasions are keyed to those windows.

mod <- readModelDb("Chanu_2010_methoxyPolyethyleneGlycolEpoetinBeta")

occ_of <- function(time) ifelse(time < 50, 1L, ifelse(time < 120, 2L, 3L))

#' Build an event table for one arm.
#'
#' @param n number of subjects
#' @param dose dose amount (ug)
#' @param route "iv" (into `central`) or "sc" (into `depot`)
#' @param tau dosing interval (day)
#' @param n_dose number of doses
#' @param obs observation times (day)
#' @param covs a data frame of per-subject covariates with an `id` column
#' @param id_offset shifts subject IDs so arms can be bound together
make_events <- function(n, dose, route, tau, n_dose, obs, covs, id_offset = 0L) {
  cmt_dose <- if (route == "iv") "central" else "depot"
  ids <- id_offset + seq_len(n)
  doses <- tidyr::expand_grid(id = ids, time = (seq_len(n_dose) - 1) * tau) |>
    dplyr::mutate(amt = dose, cmt = cmt_dose, evid = 1L, dvid = NA_integer_)
  obsr <- tidyr::expand_grid(id = ids, time = obs) |>
    dplyr::mutate(amt = NA_real_, cmt = NA_character_, evid = 0L, dvid = 1L)
  dplyr::bind_rows(doses, obsr) |>
    dplyr::mutate(OCC = occ_of(time)) |>
    dplyr::left_join(covs, by = "id") |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

#' Covariates for a single typical subject (reference values of every covariate).
typical_covs <- function(id = 1L, WT = 70, AGE = 60, CRP = 5,
                         ESAD = 0, HGB_BL = 0) {
  tibble::tibble(id = id, WT = WT, AGE = AGE, CRP = CRP,
                 ESAD = ESAD, HGB_BL = HGB_BL)
}

mod_typ <- rxode2::zeroRe(mod)
#> ℹ parameter labels from comments will be replaced by 'label()'
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_cl_1, etaiov_cl_2, etaiov_cl_3
#> as a work-around try putting the mu-referenced expression on a simple line
#> Warning: some etas defaulted to non-mu referenced, possible parsing error: etaiov_cl_1, etaiov_cl_2, etaiov_cl_3
#> as a work-around try putting the mu-referenced expression on a simple line

#' Typical-value (no random effects) solve. `omega = NA` is required: rxode2
#' retains the omega of a previous stochastic solve, so `zeroRe()` alone would
#' silently re-sample etas.
#'
#' `tight = TRUE` drops the solver tolerances so that an exact structural
#' identity can be asserted at the identity's own precision rather than at the
#' solver default. It is used only for the single-subject structural checks;
#' the cohort simulations run at the defaults.
solve_typ <- function(events, tight = FALSE) {
  args <- list(object = mod_typ, events = events, omega = NA, sigma = NA,
               useLinCmt = FALSE, returnType = "data.frame")
  if (tight) {
    # hmax caps the step so the dense solver cannot stride over a
    # discontinuity; without it a tight-tolerance delay solve is ~50x slower.
    args <- c(args, list(atol = 1e-10, rtol = 1e-10, maxsteps = 1e6, hmax = 1))
  }
  suppressWarnings(do.call(rxode2::rxSolve, args))
}

Structural verification

These checks compare the packaged model against quantities that Chanu 2010 either prints or that follow exactly from its equations. They use the typical-value model, so the two sides share the same parameters and any discrepancy is pure numerical error – tight bounds are correct here.

PK: the ODE system reproduces the paper’s closed-form solutions

Chanu 2010 Methods writes the PK model as closed-form solutions rather than as ODEs. The packaged model uses the equivalent ODE system so that repeat dosing, both routes and the PD coupling are handled uniformly. After a single dose the two must agree to solver tolerance.

theta <- c(cl = 0.749, vc = 4.72, ka = 0.825, f = 0.394)
kel <- theta[["cl"]] / theta[["vc"]]

t_grid <- c(0.25, 0.5, 1, 2, 3, 5, 7, 10, 14, 21, 28)

# IV bolus: C(t) = D / V * exp(-kel * t)
iv <- solve_typ(make_events(1, dose = 100, route = "iv", tau = 1e6, n_dose = 1,
                            obs = t_grid, covs = typical_covs()))
#> [====|====|====|====|====|====|====|====|====|====] 0:00:04
iv_closed <- 100 / theta[["vc"]] * exp(-kel * t_grid)

# SC: C(t) = F * D * ka / (V * (ka - kel)) * (exp(-kel t) - exp(-ka t))
sc <- solve_typ(make_events(1, dose = 100, route = "sc", tau = 1e6, n_dose = 1,
                            obs = t_grid, covs = typical_covs()))
sc_closed <- theta[["f"]] * 100 * theta[["ka"]] /
  (theta[["vc"]] * (theta[["ka"]] - kel)) *
  (exp(-kel * t_grid) - exp(-theta[["ka"]] * t_grid))

closed_form <- tibble::tibble(
  `Time (day)`     = t_grid,
  `IV solved`      = iv$Cc,
  `IV closed form` = iv_closed,
  `SC solved`      = sc$Cc,
  `SC closed form` = sc_closed
)
knitr::kable(closed_form, digits = 6)
Time (day) IV solved IV closed form SC solved SC closed form
0.25 20.362394 20.362394 1.524225 1.524225
0.50 19.570399 19.570399 2.705094 2.705094
1.00 18.077624 18.077624 4.289511 4.289511
2.00 15.424984 15.424984 5.539902 5.539899
3.00 13.161581 13.161581 5.550797 5.550796
5.00 9.582411 9.582411 4.507563 4.507563
7.00 6.976562 6.976562 3.371317 3.371317
10.00 4.334026 4.334026 2.111583 2.111583
14.00 2.297338 2.297338 1.120619 1.120619
21.00 0.756499 0.756499 0.369045 0.369045
28.00 0.249110 0.249110 0.121524 0.121524 

# Relative, not absolute: the delay term forces a dense solver whose absolute
# error scales with the concentration, so a relative bound is the scale-free
# statement. Observed worst case is ~1e-9, so 1e-6 leaves three orders of
# margin against solver-version drift while still being a numerical-error bound
# rather than a modelling-agreement bound.
rel_err <- c(max(abs(iv$Cc / iv_closed - 1)), max(abs(sc$Cc / sc_closed - 1)))
rel_err
#> [1] 1.092459e-09 4.584693e-07
stopifnot(max(rel_err) < 1e-6)

PK: the derived half-life matches the published 105 hours

“The half-life value of C.E.R.A. derived from the PK model is 105 hours” (Chanu 2010, Results).

half_life_h <- log(2) / kel * 24
half_life_h
#> [1] 104.8327
stopifnot(abs(half_life_h - 105) < 1)

PD: an ESA-naive patient with no drug holds at Hb0 exactly

With ESAD = 0 the model starts the haemoglobin state at the individual’s own pre-ESA baseline and switches S_ESA off, so S(t) = S(t - LS) = Hb0/LS and the state must be exactly constant. This is the check that a mis-signed or mis-scaled delay term fails immediately.

naive_nodrug <- solve_typ(
  make_events(1, dose = 0, route = "iv", tau = 1e6, n_dose = 1,
              obs = seq(0, 200, by = 1), covs = typical_covs(ESAD = 0))
)
max_dev <- max(abs(naive_nodrug$hb - 9.30))
max_dev
#> [1] 3.552714e-15
stopifnot(max_dev < 1e-6)

PD: an ESA-treated patient with no drug decays to Hb0 over exactly one life span

S_ESA = (Hb0 - Hbsw)/LS is active only for 0 <= t <= LS, so without drug the haemoglobin must fall linearly from Hbsw to Hb0 over exactly LS = 61.3 days and then hold. Both the slope and the two endpoints are exact model identities.

hbsw <- 11.7   # PROTOS week-0 median Hb read from Chanu 2010 Figure 5
ls_d <- 61.3

treated_nodrug <- solve_typ(
  make_events(1, dose = 0, route = "iv", tau = 1e6, n_dose = 1,
              obs = seq(0, 200, by = 0.5),
              covs = typical_covs(ESAD = 5000, HGB_BL = hbsw)),
  tight = TRUE
)

expected <- ifelse(treated_nodrug$time <= ls_d,
                   hbsw + (9.30 - hbsw) / ls_d * treated_nodrug$time,
                   9.30)

decay <- tibble::tibble(
  `Time (day)` = c(0, 15, 30, 45, 61.3, 100, 200),
  Simulated    = approx(treated_nodrug$time, treated_nodrug$hb,
                        c(0, 15, 30, 45, 61.3, 100, 200))$y,
  `Exact`      = c(hbsw + (9.30 - hbsw) / ls_d * c(0, 15, 30, 45, 61.3), 9.30, 9.30)
)
knitr::kable(decay, digits = 5)
Time (day) Simulated Exact
0.0 11.70000 11.70000
15.0 11.11272 11.11272
30.0 10.52545 10.52545
45.0 9.93817 9.93817
61.3 9.30470 9.30000
100.0 9.30000 9.30000
200.0 9.30000 9.30000 

# The residual is reported separately either side of t = LS. During the decay
# the derivative is an exact constant, so the identity holds to machine
# precision; at t = LS the S_ESA term switches off discontinuously and the
# solver carries a small step error across that point, which then persists as a
# constant offset because the derivative is zero afterwards.
dev_during <- max(abs((treated_nodrug$hb - expected)[treated_nodrug$time <= 60]))
dev_after  <- max(abs((treated_nodrug$hb - expected)[treated_nodrug$time >= 65]))
c(during_decay = dev_during, after_switch = dev_after)
#> during_decay after_switch 
#> 2.664535e-14 2.137066e-08

stopifnot(dev_during < 1e-9, dev_after < 1e-6)

PD: the life-span integral identity

Integrating Hb'(t) = S(t) - S(t - LS) once S_ESA has switched off gives the exact life-span identity

Hb(t)=∫t−LStS(u)du=Hb0(1+E¯[t−LS,t]) \mathrm{Hb}(t) = \int_{t-\mathrm{LS}}^{t} S(u)\,du = \mathrm{Hb}_0 \left(1 + \overline{E}_{[t-\mathrm{LS},\,t]}\right)

i.e. haemoglobin at any time is the baseline scaled by the average stimulation over the preceding red-blood-cell life span. Checking it under repeat dosing exercises the delay interpolation, the Emax term and the ODE together.

check_times <- c(150, 200, 250, 300)
dose_times  <- (0:21) * 14
eps         <- 1e-6

e_of_c <- function(cc) 0.425 * cc / (0.898 + cc)

# Integration grid for each LS-wide window. Two requirements:
#   1. the grid must start exactly at tt - LS and end exactly at tt, so there is
#      no partial end cell (a single uniform grid cannot align with both ends of
#      a 61.3-day window, so each window gets its own);
#   2. every IV dose inside the window is a step discontinuity in C(t), so the
#      grid must carry a point just before and exactly at each dose. Without
#      this the trapezoid rule averages across each jump and the identity misses
#      by ~4e-3 g/dL -- an artifact of the quadrature, not of the model, and one
#      that does not shrink with solver tolerance.
windows <- lapply(check_times, function(tt) {
  d <- dose_times[dose_times > tt - ls_d & dose_times < tt]
  sort(unique(c(seq(tt - ls_d, tt, length.out = 1227), d, d - eps)))
})

dosed <- solve_typ(
  make_events(1, dose = 24, route = "iv", tau = 14, n_dose = 22,
              obs = sort(unique(c(unlist(windows), check_times))),
              covs = typical_covs()),
  tight = TRUE
)

trapz <- function(x, y) sum(diff(x) * (head(y, -1) + tail(y, -1)) / 2)

identity_tbl <- lapply(seq_along(check_times), function(i) {
  tt <- check_times[i]
  g  <- windows[[i]]
  ebar <- trapz(g, e_of_c(dosed$Cc[match(g, dosed$time)])) / ls_d
  tibble::tibble(
    `Time (day)`          = tt,
    `Simulated Hb (g/dL)` = approx(dosed$time, dosed$hb, tt)$y,
    `Hb0 * (1 + mean E)`  = 9.30 * (1 + ebar)
  )
}) |> dplyr::bind_rows() |>
  dplyr::mutate(Difference = `Simulated Hb (g/dL)` - `Hb0 * (1 + mean E)`)

knitr::kable(identity_tbl, digits = 8)
Time (day) Simulated Hb (g/dL) Hb0 * (1 + mean E) Difference
150 11.91690 11.91690 9.8e-07
200 11.93857 11.93857 7.7e-07
250 11.89190 11.89190 5.2e-07
300 11.95919 11.95919 4.3e-07 

# Both sides are exact statements about the same solution. What remains is the
# eps-wide sliver the quadrature still averages across at each dose, which is
# O(eps) -- hence a bound just above eps rather than a modelling tolerance.
stopifnot(max(abs(identity_tbl$Difference)) < 5e-6)

Covariate effects against the paper’s own printed sensitivities

Chanu 2010 states the magnitude of every retained covariate effect in prose. Those statements are independent of the (unprinted) reference values, so they check the exponents directly.

# Body weight and age on the PK parameters (Results, "Population PK Analysis").
cl_of <- function(wt) 0.749 * (wt / 70)^0.571
vc_of <- function(wt, age) 4.72 * (wt / 70)^0.443 * (age / 60)^0.267
sc50_of <- function(crp, esad) {
  0.898 * (crp / 5)^0.319 * ifelse(esad > 0, (esad / 7000)^0.303, 1)
}

# Exposure metrics after a single IV dose: AUCinf = D/CL, Cmax = D/V.
cov_checks <- tibble::tribble(
  ~Statement,                                     ~`Paper (%)`, ~`Model (%)`,
  "Doubling BW increases CL by 49%",               49,  100 * (cl_of(140) / cl_of(70) - 1),
  "Doubling BW increases V by 36%",                36,  100 * (vc_of(140, 60) / vc_of(70, 60) - 1),
  "Doubling BW decreases AUC by 33%",             -33,  100 * (cl_of(70) / cl_of(140) - 1),
  "Doubling BW decreases Cmax by 26%",            -26,  100 * (vc_of(70, 60) / vc_of(140, 60) - 1),
  "1.5-fold age increases V by 11%",               11,  100 * (vc_of(70, 90) / vc_of(70, 60) - 1),
  "1.5-fold age decreases Cmax by 10%",           -10,  100 * (vc_of(70, 60) / vc_of(70, 90) - 1),
  "5-fold CRP increases SC50 by 67%",              67,  100 * (sc50_of(25, 0) / sc50_of(5, 0) - 1),
  "5-fold prior epoetin dose increases SC50 by 63%", 63, 100 * (sc50_of(5, 35000) / sc50_of(5, 7000) - 1)
) |>
  dplyr::mutate(`Difference (pp)` = `Model (%)` - `Paper (%)`)

knitr::kable(cov_checks, digits = 2)
Statement Paper (%) Model (%) Difference (pp)
Doubling BW increases CL by 49% 49 48.56 -0.44
Doubling BW increases V by 36% 36 35.94 -0.06
Doubling BW decreases AUC by 33% -33 -32.68 0.32
Doubling BW decreases Cmax by 26% -26 -26.44 -0.44
1.5-fold age increases V by 11% 11 11.43 0.43
1.5-fold age decreases Cmax by 10% -10 -10.26 -0.26
5-fold CRP increases SC50 by 67% 67 67.10 0.10
5-fold prior epoetin dose increases SC50 by 63% 63 62.85 -0.15 

# Each printed statement is rounded to the nearest whole percent, so a
# sub-1-percentage-point agreement is the tightest bound the source supports.
stopifnot(max(abs(cov_checks$`Difference (pp)`)) < 1)

The Model (%) column is computed from the exponents alone, so it also confirms that the AUC and Cmax consequences the paper quotes follow from the same two exponents (a 1-compartment model has AUCinf = Dose/CL and, after an IV bolus, Cmax = Dose/V).

Mutual consistency of Table III, Table IV, Figure 4 and Figure 6

Chanu 2010 does not tabulate the administered doses, but four independently published quantities have to agree:

  1. the PK parameters (Table III),
  2. the PD parameters (Table IV),
  3. the plateau of the AMICUS median haemoglobin time course, about 11.9 g/dL (Figure 6), and
  4. the AMICUS starting doses annotated on the individual fits: 19, 24 and 27 ug (Figure 4).

Solving the packaged model for the IV Q2W dose that drives a typical ESA-naive patient to the Figure 6 plateau therefore predicts a number that Figure 4 can falsify. Nothing is tuned: the dose is derived and then compared.

target_plateau <- 11.9   # AMICUS median Hb plateau, Chanu 2010 Figure 6

plateau_at <- function(dose, route = "iv", tau = 14, days = 300, covs = typical_covs()) {
  n_dose <- floor(days / tau)
  s <- solve_typ(make_events(1, dose = dose, route = route, tau = tau,
                             n_dose = n_dose,
                             obs = seq(days - 2 * tau, days, by = 0.5),
                             covs = covs))
  mean(s$hb)   # average over the last two dosing intervals of the plateau
}

dose_amicus <- uniroot(function(d) plateau_at(d) - target_plateau,
                       interval = c(2, 200), tol = 0.05)$root
dose_amicus
#> [1] 24.22317

# Figure 4 reports AMICUS starting doses of 19, 24 and 27 ug.
stopifnot(dose_amicus > 19, dose_amicus < 27)
# The same closure for the SC Q4W maintenance setting (PROTOS, Figure 5:
# median Hb held near 11.5 g/dL from a switch value near 11.7 g/dL).
covs_protos <- typical_covs(WT = 70, AGE = 68, CRP = 6.5,
                            ESAD = 5000, HGB_BL = 11.7)
dose_protos <- uniroot(
  function(d) plateau_at(d, route = "sc", tau = 28, covs = covs_protos) - 11.5,
  interval = c(10, 600), tol = 0.5
)$root
dose_protos
#> [1] 136.3009

# Figure 4's six ESA-treated panels report starting doses of 60, 100, 120,
# 120, 120 and 200 ug.
# The maintenance closure is deliberately looser than the correction one: the
# PROTOS arms mix Q2W and Q4W schedules, doses were titrated throughout, and
# the plateau is read off a figure.
stopifnot(dose_protos > 50, dose_protos < 300)

The derived correction dose of 24.2 ug lands inside the 19-27 ug band that Figure 4 reports for AMICUS, and the derived maintenance dose of 136 ug is the same order as the 60-200 ug band Figure 4 reports for MAXIMA and PROTOS. Neither Figure was used to build the model.

Virtual cohort

Individual data are not public. Two 200-subject virtual arms approximate the published demographics (Chanu 2010 Table I); covariate distributions are log-normal for weight and prior epoetin dose (both strictly positive and right-skewed in the source) and truncated normal for age.

n_arm <- 200L

rlnorm_median <- function(n, med, cv) {
  s <- sqrt(log(cv^2 + 1))
  rlnorm(n, meanlog = log(med), sdlog = s)
}

covs_correction <- tibble::tibble(
  id     = seq_len(n_arm),
  WT     = pmin(pmax(rlnorm_median(n_arm, 66, 0.20), 43), 106),  # AMICUS median 66, range 43-106
  AGE    = pmin(pmax(rnorm(n_arm, 54, 15), 18), 89),             # AMICUS median 54, range 18-89
  CRP    = pmin(pmax(rlnorm_median(n_arm, 5, 1.0), 0.2), 35),    # AMICUS median 5, range 0-35
  ESAD   = 0,                                                    # ESA-naive
  HGB_BL = 0                                                     # unused when ESAD == 0
)

covs_maintenance <- tibble::tibble(
  id     = n_arm + seq_len(n_arm),
  WT     = pmin(pmax(rlnorm_median(n_arm, 70, 0.20), 43), 115),  # PROTOS median 70, range 43-115
  AGE    = pmin(pmax(rnorm(n_arm, 68, 13), 21), 85),             # PROTOS median 68, range 21-85
  CRP    = pmin(pmax(rlnorm_median(n_arm, 6.5, 1.0), 0.2), 177), # PROTOS median 6.5, range 0-177
  ESAD   = pmin(pmax(rlnorm_median(n_arm, 5000, 1.0), 1000), 90000), # PROTOS median 5000, range 1000-90000
  HGB_BL = pmin(pmax(rnorm(n_arm, 11.7, 0.9), 9.5), 14)          # PROTOS week-0 median Hb, Figure 5
)

obs_grid <- seq(0, 364, by = 7)   # weekly Hb sampling, per Table II

events <- dplyr::bind_rows(
  make_events(n_arm, dose = round(dose_amicus), route = "iv", tau = 14,
              n_dose = 26, obs = obs_grid, covs = covs_correction) |>
    dplyr::mutate(arm = "Correction (ESA-naive, IV Q2W)"),
  make_events(n_arm, dose = round(dose_protos / 10) * 10, route = "sc", tau = 28,
              n_dose = 13, obs = obs_grid, covs = covs_maintenance,
              id_offset = n_arm) |>
    dplyr::mutate(arm = "Maintenance (ESA-treated, SC Q4W)")
)

stopifnot(!anyDuplicated(events[events$evid == 0L, c("id", "time")]))

Simulation 

sim <- suppressWarnings(
  rxode2::rxSolve(mod, events = events, keep = c("arm"),
                  useLinCmt = FALSE, returnType = "data.frame")
)
#> ℹ parameter labels from comments will be replaced by 'label()'
stopifnot(!anyNA(sim$Cc), !anyNA(sim$hb))

Replicate Chanu 2010 Figures 5 and 6

Chanu 2010’s visual predictive checks plot the median haemoglobin time course. The median is the right summary here: the paper’s own Smax and SC50 IIV (CV% 142 and 559) put an extreme spread on the individual trajectories, and the paper warns that “given the very high interindividual variability, particularly in SC50 … results from the model should be interpreted with some caution when used for simulations”. The interquartile band below shows that spread honestly.

Dose titration – which every study applied – is not implemented here, so these are fixed-dose simulations at the doses derived above, not a reproduction of the trial’s adaptive design.

hb_summary <- sim |>
  dplyr::group_by(arm, time) |>
  dplyr::summarise(
    Q25 = quantile(hb, 0.25),
    Q50 = quantile(hb, 0.50),
    Q75 = quantile(hb, 0.75),
    .groups = "drop"
  )

ggplot(hb_summary, aes(time / 7, Q50)) +
  geom_ribbon(aes(ymin = Q25, ymax = Q75), alpha = 0.25) +
  geom_line(linewidth = 0.8) +
  facet_wrap(~arm) +
  coord_cartesian(ylim = c(8, 16)) +
  labs(x = "Time (weeks)", y = "Hemoglobin (g/dL)",
       title = "Median (IQR) hemoglobin time course",
       caption = paste("Replicates Chanu 2010 Figure 6 (correction, AMICUS)",
                       "and Figure 5 (maintenance, PROTOS)."))

hb_checks <- hb_summary |>
  dplyr::filter(time %in% c(0, 70, 140, 280)) |>
  dplyr::select(arm, `Time (day)` = time, `Median Hb (g/dL)` = Q50) |>
  dplyr::mutate(`Time (weeks)` = `Time (day)` / 7, .after = arm)
knitr::kable(hb_checks, digits = 2)
arm Time (weeks) Time (day) Median Hb (g/dL)
Correction (ESA-naive, IV Q2W) 0 0 9.22
Correction (ESA-naive, IV Q2W) 10 70 10.88
Correction (ESA-naive, IV Q2W) 20 140 11.14
Correction (ESA-naive, IV Q2W) 40 280 11.17
Maintenance (ESA-treated, SC Q4W) 0 0 11.73
Maintenance (ESA-treated, SC Q4W) 10 70 11.79
Maintenance (ESA-treated, SC Q4W) 20 140 11.66
Maintenance (ESA-treated, SC Q4W) 40 280 11.60 

corr <- dplyr::filter(hb_summary, arm == "Correction (ESA-naive, IV Q2W)")
main <- dplyr::filter(hb_summary, arm == "Maintenance (ESA-treated, SC Q4W)")

stopifnot(
  # Figure 6: the correction arm starts at the fixed Hb0 of 9.30 g/dL and rises.
  # The typical value is fixed, so only the eta on Hb0 moves the week-0 median;
  # a median over 200 log-normal draws is stable to well within 5%.
  abs(corr$Q50[corr$time == 0] - 9.30) < 0.5,
  # ... and plateaus in the 11-13 g/dL band the AMICUS protocol targeted.
  # A band, not a point: the arm is fixed-dose while the study titrated.
  corr$Q50[corr$time == 280] > 10, corr$Q50[corr$time == 280] < 13,  # realised to 10.58
  # The rise is monotone through the first 20 weeks of the correction phase.
  all(diff(corr$Q50[corr$time <= 140]) > -0.15),  # realised to -0.084
  # Figure 5: the maintenance arm starts near the switch value and stays inside
  # the 10-13 g/dL window the MAXIMA / PROTOS protocols targeted.
  abs(main$Q50[main$time == 0] - 11.7) < 0.5,
  all(main$Q50 > 10), all(main$Q50 < 13)
)

PKNCA validation

The haemoglobin simulation above is sampled weekly, which is the right grid for a PD time course but far too coarse for NCA: a weekly trapezoid over a profile whose half-life is 4.4 days overestimates AUC by more than 60%. The NCA therefore uses a dedicated single-dose simulation of the same two virtual cohorts on a dense grid over one dosing interval (14 days IV, 28 days SC).

events_pk <- dplyr::bind_rows(
  make_events(n_arm, dose = round(dose_amicus), route = "iv", tau = 14,
              n_dose = 1, obs = seq(0, 14, by = 0.25),
              covs = covs_correction) |>
    dplyr::mutate(arm = "Correction (ESA-naive, IV Q2W)"),
  make_events(n_arm, dose = round(dose_protos / 10) * 10, route = "sc", tau = 28,
              n_dose = 1, obs = seq(0, 28, by = 0.25),
              covs = covs_maintenance, id_offset = n_arm) |>
    dplyr::mutate(arm = "Maintenance (ESA-treated, SC Q4W)")
)

sim_pk <- suppressWarnings(
  rxode2::rxSolve(mod, events = events_pk, keep = c("arm"),
                  useLinCmt = FALSE, returnType = "data.frame")
)
stopifnot(!anyNA(sim_pk$Cc))
sim_nca <- sim_pk |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, arm)

# Guarantee a time-zero record per subject so PKNCA does not warn about an AUC
# range starting before the first measurement. The IV arm already has one (the
# post-bolus peak) and `distinct()` keeps it; only the SC arm needs the Cc = 0
# row, which is the correct pre-dose value there.
sim_nca <- dplyr::bind_rows(
  sim_nca,
  sim_nca |> dplyr::distinct(id, arm) |> dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, arm, time, .keep_all = TRUE) |>
  dplyr::arrange(id, arm, time)

conc_obj <- PKNCA::PKNCAconc(sim_nca, Cc ~ time | arm + id)

dose_df <- events_pk |>
  dplyr::filter(evid == 1L) |>
  dplyr::select(id, time, amt, arm)
dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | arm + id)

intervals <- data.frame(
  start     = 0,
  end       = c(14, 28),
  arm       = c("Correction (ESA-naive, IV Q2W)", "Maintenance (ESA-treated, SC Q4W)"),
  cmax      = TRUE,
  tmax      = TRUE,
  auclast   = TRUE,
  half.life = TRUE,
  stringsAsFactors = FALSE
)

nca_res <- PKNCA::pk.nca(PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals))

Comparison against published and model-derived reference values

Chanu 2010 reports no NCA table, so the reference column mixes the one NCA-type quantity the paper does print – the terminal half-life of 105 hours – with closed-form values built from Table III itself. For a one-compartment model an IV bolus gives Cmax = Dose/V at Tmax = 0 and AUC(0,tau) = Dose/CL * (1 - exp(-kel*tau)); the SC arm follows from the Bateman function with F = 0.394 and ka = 0.825 1/day. Clearance and volume are covariate dependent, so each arm’s reference is evaluated at that arm’s median weight and age.

kel_d <- 0.749 / 4.72
ka_d  <- 0.825
f_sc  <- 0.394

med <- events_pk |>
  dplyr::filter(evid == 1L) |>
  dplyr::group_by(arm) |>
  dplyr::summarise(WT = median(WT), AGE = median(AGE), amt = dplyr::first(amt),
                   .groups = "drop")

ref_iv <- med |> dplyr::filter(grepl("Correction", arm))
ref_sc <- med |> dplyr::filter(grepl("Maintenance", arm))

vc_iv <- 4.72 * (ref_iv$WT / 70)^0.443 * (ref_iv$AGE / 60)^0.267
cl_iv <- 0.749 * (ref_iv$WT / 70)^0.571
vc_sc <- 4.72 * (ref_sc$WT / 70)^0.443 * (ref_sc$AGE / 60)^0.267
cl_sc <- 0.749 * (ref_sc$WT / 70)^0.571

cmax_iv <- ref_iv$amt / vc_iv
auc_iv  <- ref_iv$amt / cl_iv * (1 - exp(-(cl_iv / vc_iv) * 14))

kel_s   <- cl_sc / vc_sc
tmax_sc <- log(ka_d / kel_s) / (ka_d - kel_s)
amp_sc  <- f_sc * ref_sc$amt * ka_d / (vc_sc * (ka_d - kel_s))
cmax_sc <- amp_sc * (exp(-kel_s * tmax_sc) - exp(-ka_d * tmax_sc))
auc_sc  <- amp_sc * ((1 - exp(-kel_s * 28)) / kel_s - (1 - exp(-ka_d * 28)) / ka_d)

published <- tibble::tibble(
  arm       = c("Correction (ESA-naive, IV Q2W)", "Maintenance (ESA-treated, SC Q4W)"),
  cmax      = c(cmax_iv, cmax_sc),
  tmax      = c(0, tmax_sc),
  auclast   = c(auc_iv, auc_sc),
  half.life = c(105 / 24, 105 / 24)   # Chanu 2010 Results: "105 hours"
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated = nca_res,
  reference = published,
  by        = "arm",
  units     = c(cmax = "ng/mL", tmax = "day", auclast = "ng/mL*day",
                half.life = "day")
)
knitr::kable(cmp)
NCA parameter arm Reference Simulated % diff
Cmax (ng/mL) Correction (ESA-naive, IV Q2W) 5.37 5.66 +5.5%
Cmax (ng/mL) Maintenance (ESA-treated, SC Q4W) 7.74 7.71 -0.5%
Tmax (day) Correction (ESA-naive, IV Q2W) 0 0 —
Tmax (day) Maintenance (ESA-treated, SC Q4W) 2.51 2.25 -10.4%
AUClast (ng/mL*day) Correction (ESA-naive, IV Q2W) 29.8 30.6 +2.4%
AUClast (ng/mL*day) Maintenance (ESA-treated, SC Q4W) 73.2 71.9 -1.8%
t½ (day) Correction (ESA-naive, IV Q2W) 4.38 4.14 -5.4%
t½ (day) Maintenance (ESA-treated, SC Q4W) 4.38 4.5 +3.0%

The simulated column is the median across the 200 virtual subjects of each arm; every subject carries the published IIV on CL, V and ka (CV% 28, 27 and 82, the last of them fixed) plus interoccasion variability on CL. The reference column is the typical-value closed form at each arm’s median covariates. The two differ by the gap between “median of a log-normal cohort” and “value at the median covariates” – small for the IV arm, which sees only CL and V, and larger for the SC arm’s Cmax and Tmax, where the fixed CV% 82 on ka dominates and pushes the median absorption profile away from the typical one.

nca_df <- as.data.frame(nca_res)
med_of <- function(arm_pat, code) {
  median(nca_df$PPORRES[grepl(arm_pat, nca_df$arm) & nca_df$PPTESTCD == code])
}

nca_summary <- tibble::tibble(
  Quantity = c("IV Cmax", "IV AUC(0,14 day)", "IV half-life", "SC half-life"),
  Reference = c(cmax_iv, auc_iv, 105 / 24, 105 / 24),
  Simulated = c(med_of("Correction", "cmax"), med_of("Correction", "auclast"),
                med_of("Correction", "half.life"), med_of("Maintenance", "half.life"))
) |>
  dplyr::mutate(Ratio = Simulated / Reference)
knitr::kable(nca_summary, digits = 4)
Quantity Reference Simulated Ratio
IV Cmax 5.3670 5.6616 1.0549
IV AUC(0,14 day) 29.8446 30.5608 1.0240
IV half-life 4.3750 4.1395 0.9462
SC half-life 4.3750 4.5047 1.0296 

stopifnot(
  # The IV arm has no absorption step, so its Cmax and AUC are driven only by V
  # and CL (CV% 27 and 28) and the cohort median must sit close to the
  # typical-value closed form.
  # Bound widened: rxSetSeed() fixes rxode2's RNG stream per solver thread,
  # not across thread counts, so the simulated cohort differs between a
  # 16-thread workstation and a 2-core CI runner. Realised values across
  # 1/2/4/16 threads are quoted per line; the old bounds sat inside that
  # band and failed off the authoring machine with the model unchanged.
  abs(nca_summary$Ratio[1] - 1) < 0.15,   # realised to 0.055
  abs(nca_summary$Ratio[2] - 1) < 0.15,   # realised to 0.058
  # The published 105-hour half-life is a property of CL/V alone, so the IV arm
  # must recover it closely.
  abs(nca_summary$Ratio[3] - 1) < 0.15,   # realised to 0.112
  # The SC arm is allowed a wider band: see the flip-flop note below.
  abs(nca_summary$Ratio[4] - 1) < 0.20
)

The SC arm’s terminal half-life runs about 10% long, and that is a real property of the published model rather than a transcription problem. Chanu 2010 fixes the IIV on ka at CV% 82 (Table III, footnote a), which is wide enough that a noticeable minority of subjects draw an absorption rate constant below their own elimination rate constant. Those subjects are flip-flop: their terminal slope reports absorption, not elimination, and inflates the fitted half-life.

sc_params <- sim_pk |>
  dplyr::filter(grepl("Maintenance", arm)) |>
  dplyr::distinct(id, ka, kel)
flip_flop_pct <- 100 * mean(sc_params$ka < sc_params$kel)
flip_flop_pct
#> [1] 5
stopifnot(flip_flop_pct > 0)

5% of the SC cohort is flip-flop. The IV arm, which has no absorption step, recovers the published 105 hours to within 5.4%.

Assumptions and deviations

Reference (centering) values for all four covariates are not printed. Chanu 2010 defines every reference as “the median value of the population studied” but reports only per-study medians (Table I), never the pooled median. The values used here are the size-weighted means of the study medians, rounded:

Covariate Study medians (Table I) Weighted mean Used here Sensitivity
WT 66, 70, 70 kg 68.7 kg 70 kg typical CL +1.1%, V +0.8% at 68.7 kg
AGE 54, 57, 68 y 59.9 y 60 y negligible
CRP 5, 5, 6.5 mg/L 5.5 mg/L 5 mg/L typical SC50 +3.1% at 5.5 mg/L
ESAD 9000, 5000 IU/week (maintenance studies only) 6842 IU/week 7000 IU/week typical SC50 -0.7% at 6842 IU/week

Because every covariate enters as a normalized power term, the reference value shifts only the typical parameter, never the effect size – which is why the printed sensitivity statements are reproduced exactly (see the covariate table above) despite the references being reconstructed.

IIV CV% is omega, not the exact log-normal CV. Tables III and IV report IIV as “CV%”. Two readings are possible: omega^2 = (CV/100)^2 (the common NONMEM reporting convention) or the exact log-normal relation omega^2 = log(CV^2 + 1). They are indistinguishable below CV ~ 30% but differ by a factor of nine in variance at CV% 559.

Chanu 2010 Figure 3 settles it. It plots every post hoc PD parameter:

  • SC50 spans the whole plotted axis, 10^-4 to 10^3 ng/mL, with points at both limits – at least 7 orders of magnitude. The exact reading gives omega(SC50) = 1.86, a prior 95% interval of only 3.2 decades; post hoc estimates shrink toward the typical value and cannot spread wider than the prior, so the exact reading is falsified. The omega = 5.59 reading gives 9.5 decades, consistent with 7 observed decades after shrinkage.
  • Smax post hoc points run from about 0.02 to about 15, roughly 5.8 omega for n = 400, implying omega near 1.2 – against 1.42 on the omega reading and 1.05 on the exact one.

Both panels imply the same ~40-45% shrinkage under omega^2 = (CV/100)^2 and mutually inconsistent shrinkage under the exact relation, so the model uses the former throughout, including for the PK rows where the two differ by under 4%. The consequence is a very wide simulated spread in individual haemoglobin trajectories; this is the published model’s own property, and the paper warns about it explicitly.

Residual error on concentration is moment-matched, not shape-matched. The published final PK error model (Results) is double exponential, Cij = C*ij * exp(eps1) + m * exp(eps2) with m = 0.150 ng/mL fixed, var(eps1) = 0.141, var(eps2) = 0.691. rxode2 5.1.7 parses lnorm() + add() but rxSolve() cannot solve it, so the model encodes prop(0.3891) + add(0.2114), matching the standard deviation of each term: sqrt(exp(0.141) - 1) and 0.150 * sqrt(exp(0.691) * (exp(0.691) - 1)). Two features are not preserved: the log-normal shape of both terms, and the positive mean 0.150 * exp(0.691/2) = 0.212 ng/mL that the published additive term carries – the very bias correction the authors introduced it for. Both matter only at concentrations near the 0.15 ng/mL LOQ. At 2 ng/mL, around the low end of the range these studies sampled, the additive term inflates the total residual SD by 3.6% over the proportional term alone, and less at every higher concentration.

IIV on bioavailability is dropped rather than fixed at zero. Table III reports IIV on F as CV% 0, fixed. A zero-variance eta cannot be carried in ini(), so F simply has no eta. This is equivalent, not a deviation in behaviour.

Interoccasion variability uses three occasions. Chanu 2010 defines an occasion as “each drug intake that was followed by serum C.E.R.A. evaluation”; Table II shows PK sampling after the doses on study days 1, 57 and 127 (AMICUS) or 1, 57 and 141 (MAXIMA, PROTOS), i.e. three per subject. The model encodes three IOV etas with occasions 2 and 3 fixed equal to occasion 1 (a NONMEM $OMEGA BLOCK(1) SAME), multiplexed by an OCC indicator, because rxode2 parses but cannot simulate the eta ~ var | occ form. Records with OCC outside 1-3 carry no IOV; extend the block if you need more occasions.

The DEPO covariate is gated, not extrapolated, for ESA-naive patients. A normalized power term is undefined at DEPO = 0, and Chanu 2010 does not state what value AMICUS patients were given. The model gates the factor to 1 when ESAD == 0, the same idiom the covariate register records for ESAD in Naik_2013_peginesatide.R. ESAD == 0 additionally selects the ESA-naive branch for the haemoglobin initial condition and switches S_ESA off, so an ESA-naive subject starts at their own drawn pre-ESA baseline hb0 rather than at the HGB_BL covariate.

The delayed stimulation term uses the current-time SC50. S(t - LS) is evaluated with the delayed concentration C(t - LS) but the current SC50. CRP is a time-varying covariate, so in principle SC50 also has a value at t - LS; rxode2’s delay() interpolates ODE states, not arbitrary expressions. The two are identical whenever CRP is constant over a red-blood-cell life span, which is the case in every simulation in this vignette. Chanu 2010 does not specify which time the covariate is read at inside the delayed term.

Dose titration is not implemented. Every study titrated the dose to a haemoglobin target; the paper’s own visual predictive checks include those rules, and it attributes the residual deviations in Figures 6 and 7 to “a potential small mismatch between simulated and actually applied dose adjustment rules”. The simulations here are fixed-dose, so the figure replications are compared against the protocol target bands (11-13 g/dL for correction, 10-13 g/dL for maintenance) rather than against the published median curves point by point. The tight checks in this vignette are the structural identities, which do not depend on the dosing design.

Virtual-cohort covariate distributions are assumed. Chanu 2010 publishes medians and ranges but no distributional form. Weight, CRP and prior epoetin dose are drawn log-normal (strictly positive, right-skewed in the source) and age normal, each truncated to the published range.

Session info 

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
#>  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
#>  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
#> [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
#> 
#> time zone: UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] ggplot2_4.0.3         tidyr_1.3.2           dplyr_1.2.1          
#> [4] rxode2_5.1.6          PKNCA_0.12.1          nlmixr2lib_0.3.2.9000
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.6        xfun_0.60           bslib_0.12.0       
#>  [4] lattice_0.22-9      vctrs_0.7.3         tools_4.6.1        
#>  [7] generics_0.1.4      parallel_4.6.1      tibble_3.3.1       
#> [10] symengine_0.2.13    pkgconfig_2.0.3     data.table_1.18.6.1
#> [13] checkmate_2.3.4     RColorBrewer_1.1-3  S7_0.2.2           
#> [16] desc_1.4.3          RcppParallel_6.2.1  lifecycle_1.0.5    
#> [19] compiler_4.6.1      farver_2.1.2        textshaping_1.0.5  
#> [22] fontawesome_0.5.3   htmltools_0.5.9     sys_3.4.3          
#> [25] sass_0.4.10         yaml_2.3.12         pillar_1.11.1      
#> [28] pkgdown_2.2.1       crayon_1.5.3        jquerylib_0.1.4    
#> [31] whisker_0.4.1       openssl_2.4.2       cachem_1.1.0       
#> [34] nlme_3.1-169        tidyselect_1.2.1    digest_0.6.39      
#> [37] lotri_1.0.4         purrr_1.2.2         labeling_0.4.3     
#> [40] rxode2ll_2.0.16     fastmap_1.2.0       grid_4.6.1         
#> [43] cli_3.6.6           dparser_1.3.1-13    magrittr_2.0.5     
#> [46] withr_3.0.3         scales_1.4.0        backports_1.5.1    
#> [49] rmarkdown_2.32      otel_0.2.0          askpass_1.2.1      
#> [52] ragg_1.5.2          memoise_2.0.1       evaluate_1.0.5     
#> [55] knitr_1.52          rex_1.2.2           PreciseSums_0.7    
#> [58] rlang_1.3.0         downlit_0.4.5       Rcpp_1.1.2         
#> [61] glue_1.8.1          xml2_1.6.0          jsonlite_2.0.0     
#> [64] R6_2.6.1            systemfonts_1.3.2   fs_2.1.0