---
title: "Darbepoetin alfa (Takama 2007)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Darbepoetin alfa (Takama 2007)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5)
library(nlmixr2lib)
library(PKNCA)
library(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
```

# Model and source

* Citation: `r rxode2::rxode(readModelDb("Takama_2007_darbepoetin"))$reference`
* Description: `r rxode2::rxode(readModelDb("Takama_2007_darbepoetin"))$description`
* Article: <https://doi.org/10.1111/j.1365-2125.2006.02756.x>

# Population

131 Japanese adult dialysis patients (63 receiving intermittent haemodialysis,
68 receiving peritoneal dialysis) were pooled from four clinical studies for
the population pharmacokinetic analysis. Baseline median (range) characteristics
from Takama 2007 Table 1: weight 54.7 kg (35.5-132.0); age 60 y (23-84);
serum albumin 3.8 g/dL (2.1-4.8); serum creatinine 10.60 mg/dL (3.68-16.50);
red blood cell counts 317 x 10^4/uL (range 211-422); platelet counts 17.7 x
10^4/uL (5.1-50.1); white blood cell counts 5500/uL (2400-12500). The cohort
is 81 male / 50 female (38.2 percent female). Subjects received 10-90 ug
darbepoetin alfa intravenously, single or multiple administration every 1
or 2 weeks. HD subjects had serial sampling (predose and 0.5, 1, 2, 5, 8,
12, 24, 48, 96, 168 h after dosing, plus 336 h for the 90 ug group); PD
subjects had sparse sampling (predose in weeks 1-4 plus optional 0.5-1 h
post-dose in week 2). Serum darbepoetin alfa was quantified with the
Quantikine in-vitro-diagnostic rHuEPO ELISA kit (R&D Systems) against a
darbepoetin alfa standard curve; the assay LOQ was 0.078 ng/mL and inter-
assay precision was within 20 percent for clinical samples. 917
concentration measurements were available.

The same information is available programmatically via the model's
`population` metadata
(`readModelDb("Takama_2007_darbepoetin")$population`).

# Source trace

The per-parameter origin is recorded as an in-file comment next to each
`ini()` entry in
`inst/modeldb/specificDrugs/Takama_2007_darbepoetin.R`. The table below
collects them in one place for review.

| Equation / parameter                                              | Value                              | Source location |
|-------------------------------------------------------------------|------------------------------------|-----------------|
| `lcl`     (clearance CL at WT = 54 kg)                            | `log(0.0807)` (CL = 0.0807 L/h)    | Table 4 row "theta_CL" (SE 0.00195) |
| `lq`      (intercompartmental clearance Q)                        | `log(0.0616)` (Q = 0.0616 L/h)     | Table 4 row "theta_Q"  (SE 0.00755) |
| `lvc`     (central volume V1 at WT = 54 kg, HD)                   | `log(2.51)`   (V1 = 2.51 L)        | Table 4 row "theta_V1" (SE 0.0629)  |
| `lvp`     (peripheral volume V2)                                  | `log(0.522)`  (V2 = 0.522 L)       | Table 4 row "theta_V2" (SE 0.0407)  |
| `lrbase`  (baseline endogenous erythropoietin level k0)           | `log(0.167)`  (k0 = 0.167 ng/mL)   | Table 4 row "theta_k0" (SE 0.00629) |
| `e_wt_cl` (linear-deviation coefficient of (WT - 54) on CL)        | 0.0195 (1/kg)                      | Table 4 row "theta_CL_WT" (SE 0.00246; 95% CI 0.0147-0.0243) |
| `e_wt_vc` (linear-deviation coefficient of (WT - 54) on V1)        | 0.0163 (1/kg)                      | Table 4 row "theta_V1_WT" (SE 0.00214; 95% CI 0.0121-0.0205) |
| `e_perit_dial_vc` (additive PD-modality fractional increment on V1) | 0.170                            | Table 4 row "theta_V1_DIA" (SE 0.0501; 95% CI 0.0718-0.268) |
| `etalcl`  variance                                                | 0.04153 = log(0.206^2 + 1)         | Table 4 row "omega_CL"  = 20.6 percent CV |
| `etalq`   variance                                                | 0.56449 = log(0.871^2 + 1)         | Table 4 row "omega_Q"   = 87.1 percent CV |
| `etalvc`  variance                                                | 0.04644 = log(0.218^2 + 1)         | Table 4 row "omega_V1"  = 21.8 percent CV |
| `etalvp`  variance                                                | 0.21205 = log(0.486^2 + 1)         | Table 4 row "omega_V2"  = 48.6 percent CV |
| `etalrbase` variance                                              | 0.18741 = log(0.454^2 + 1)         | Table 4 row "omega_k0"  = 45.4 percent CV |
| `addSd`                                                           | 0.0634 ng/mL                       | Table 4 row "sigma_2" (additive) |
| `propSd`                                                          | 0.0653 (= 6.53 percent)            | Table 4 row "sigma_1" (CCV / proportional) |
| `CL = theta_CL * [1 + theta_CL_WT * (WT - 54)]`                   | n/a                                | Methods, "Population model building" |
| `V1 = theta_V1 * [1 + theta_V1_WT * (WT - 54) + theta_V1_DIA * DIA]` | n/a                            | Methods, "Population model building" |
| `d/dt(central) = -kel*central - k12*central + k21*peripheral1`    | n/a                                | Standard two-compartment IV bolus open model; Methods "Step 1: basic pharmacokinetic modelling" |
| `Cc = central / V1 + k0`                                          | n/a                                | Methods, "Cp_ij = Cpm_ij * (1 + eps_1) + eps_2 with k0_j as the j-th individual's baseline" |

# Virtual cohort

Original individual-level data are not publicly available. The figures below
use a virtual population whose covariate distributions approximate the
published trial demographics (Takama 2007 Table 1).

```{r cohort}
set.seed(20060815L)  # paper OnlineEarly publication date

# Per-subject weight draws by dialysis modality. Table 1 reports
# medians and ranges; we draw from truncated normals matched to the
# group-specific median +/- approx half-range / 2 spread, clamped to
# the Table 1 ranges so simulations stay within the cohort.
make_cohort <- function(n, modality, dose_ug, id_offset = 0L) {
  if (modality == "HD") {
    wt_mean <- 53.8; wt_sd <- 7;   wt_lo <- 35.5; wt_hi <- 67.0
    perit_dial <- 0L
  } else {
    wt_mean <- 57.4; wt_sd <- 12;  wt_lo <- 38.8; wt_hi <- 132.0
    perit_dial <- 1L
  }
  tibble::tibble(
    id         = id_offset + seq_len(n),
    WT         = pmin(wt_hi, pmax(wt_lo, round(rnorm(n, mean = wt_mean, sd = wt_sd), 1))),
    PERIT_DIAL = perit_dial,
    modality   = modality,
    treatment  = paste(modality, dose_ug, "ug"),
    amt        = dose_ug
  )
}

cohort <- dplyr::bind_rows(
  make_cohort(n = 30L, modality = "HD", dose_ug = 60, id_offset =   0L),
  make_cohort(n = 30L, modality = "PD", dose_ug = 60, id_offset = 100L)
)

stopifnot(!anyDuplicated(cohort$id))

knitr::kable(
  cohort |>
    dplyr::group_by(modality, treatment) |>
    dplyr::summarise(
      n            = dplyr::n(),
      mean_WT      = round(mean(WT), 1),
      median_WT    = round(median(WT), 1),
      .groups      = "drop"
    ),
  caption = "Per-cohort weight summary (virtual)."
)
```

```{r events}
# Build dose + observation event table. IV bolus into central at t = 0,
# observation grid mirroring the paper's HD serial-sampling design
# (predose, 0.5, 1, 2, 5, 8, 12, 24, 48, 96, 168 h plus a 336 h tail).
obs_times <- c(0, 0.5, 1, 2, 5, 8, 12, 24, 48, 96, 168, 240, 336)

dose_rows <- cohort |>
  dplyr::transmute(
    id = id, time = 0,
    amt = amt, evid = 1L, cmt = "central",
    WT, PERIT_DIAL, modality, treatment
  )

obs_rows <- cohort |>
  tidyr::expand_grid(time = obs_times) |>
  dplyr::transmute(
    id = id, time = time,
    amt = 0, evid = 0L, cmt = "central",
    WT, PERIT_DIAL, modality, treatment
  )

events <- dplyr::bind_rows(dose_rows, obs_rows) |>
  dplyr::arrange(id, time, dplyr::desc(evid))

stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
```

# Simulation

```{r simulate}
mod <- readModelDb("Takama_2007_darbepoetin")

sim <- rxode2::rxSolve(
  mod, events = events,
  keep = c("treatment", "modality", "WT", "PERIT_DIAL")
) |>
  as.data.frame()
```

For deterministic replication (reproducing the typical-subject profile
without between-subject variability), we also simulate the typical
trajectory for the median-weight HD and PD subjects on a fine grid.

```{r simulate-typical}
mod_typical <- mod |> rxode2::zeroRe()

simulate_typical_one <- function(modality, dose_ug, id) {
  perit_dial <- if (modality == "HD") 0L else 1L
  wt <- if (modality == "HD") 53.8 else 57.4   # group-specific median
  dose_row <- tibble::tibble(
    id = id, time = 0, amt = dose_ug, evid = 1L, cmt = "central",
    WT = wt, PERIT_DIAL = perit_dial
  )
  obs_rows <- tibble::tibble(
    id = id, time = seq(0, 336, by = 1),
    amt = 0, evid = 0L, cmt = "central",
    WT = wt, PERIT_DIAL = perit_dial
  )
  ev <- dplyr::bind_rows(dose_row, obs_rows) |>
    dplyr::arrange(time, dplyr::desc(evid))
  rxode2::rxSolve(mod_typical, events = ev) |>
    as.data.frame() |>
    dplyr::mutate(modality = modality,
                  treatment = paste(modality, dose_ug, "ug"))
}

typ_sim <- dplyr::bind_rows(
  simulate_typical_one("HD", 60, 1L),
  simulate_typical_one("PD", 60, 2L)
)
```

# Replicate Figure 1 of the source

Takama 2007 Figure 1 shows observed darbepoetin alfa serum concentrations
versus time, with separate panels for HD and PD patients. The replication
here uses the typical-subject profile for a 60 ug IV bolus, plus
stochastic-sample concentrations from the virtual cohort, on a log y-axis
to match the source.

```{r figure-1, fig.cap = "Replicates Figure 1 of Takama 2007: simulated serum darbepoetin alfa concentration vs. time after a single 60 ug IV bolus, faceted by dialysis modality. Solid line = typical subject (median-weight); points = stochastic-sample subjects."}
ggplot(sim |> dplyr::filter(time > 0),
       aes(time, Cc, group = id, colour = modality)) +
  geom_point(alpha = 0.35, size = 1) +
  geom_line(data = typ_sim |> dplyr::filter(time > 0),
            aes(time, Cc, group = modality), inherit.aes = FALSE,
            linewidth = 1, colour = "black") +
  facet_wrap(~ modality) +
  scale_y_log10() +
  scale_x_continuous(breaks = c(0, 24, 48, 96, 168, 240, 336)) +
  labs(
    x       = "Time after dose (h)",
    y       = "Serum darbepoetin alfa (ng/mL, log scale)",
    colour  = "Modality",
    caption = "Replicates Figure 1 of Takama 2007 (single 60 ug IV bolus per subject)."
  ) +
  theme_minimal() +
  theme(legend.position = "none")
```

# Replicate Figure 3 of the source

Takama 2007 Figure 3 shows the distribution of the initial serum
concentration (C0) and AUC across weight tertiles (panels a, b) and
between dialysis modalities (panels c, d). The paper's narrative
emphasises that "there were no obvious differences in the distributions
and the geometric mean values in the C0 and the AUC between weight
categories and between dialysis techniques", and concludes that no
dosage adjustment is warranted across the covariate ranges studied.

```{r figure-3}
# Per-subject summaries: C0 (post-IV-bolus serum concentration at t = 0)
# and AUC0-336 from trapezoidal integration on the simulation grid.
# rxSolve already returns the post-dose value at time == 0 for an IV
# bolus into central, so no evid filter is needed.
post_dose_t0 <- sim |>
  dplyr::filter(time == 0) |>
  dplyr::transmute(id, modality, WT, C0 = Cc)

per_subject_auc <- sim |>
  dplyr::group_by(id, modality, WT) |>
  dplyr::summarise(
    auc_trap = sum(diff(time) * (head(Cc, -1) + tail(Cc, -1)) / 2),
    .groups  = "drop"
  )

c0_auc <- post_dose_t0 |>
  dplyr::inner_join(per_subject_auc |> dplyr::select(id, auc_trap),
                    by = "id") |>
  dplyr::mutate(
    weight_tertile = dplyr::case_when(
      WT <  quantile(WT, 1/3) ~ "1 (low)",
      WT >= quantile(WT, 2/3) ~ "3 (high)",
      TRUE                    ~ "2 (mid)"
    )
  )
```

```{r figure-3a, fig.cap = "Replicates Figure 3a of Takama 2007: simulated initial serum concentration C0 by weight tertile. Bars are geometric means."}
gm <- function(x) exp(mean(log(x)))

c0_by_tertile <- c0_auc |>
  dplyr::group_by(weight_tertile) |>
  dplyr::summarise(C0_gm = gm(C0), .groups = "drop")

ggplot(c0_auc, aes(weight_tertile, C0)) +
  geom_jitter(width = 0.2, alpha = 0.5) +
  geom_crossbar(data = c0_by_tertile,
                aes(x = weight_tertile, y = C0_gm, ymin = C0_gm, ymax = C0_gm),
                width = 0.5, colour = "red") +
  labs(
    x       = "Weight tertile",
    y       = "Initial serum concentration C0 (ng/mL)",
    caption = "Replicates Figure 3a of Takama 2007."
  ) +
  theme_minimal()
```

```{r figure-3c, fig.cap = "Replicates Figure 3c-d of Takama 2007: simulated initial serum concentration C0 and AUC0-inf by dialysis modality (HD vs PD). The paper concludes there is no marked difference between modalities."}
c0_by_modality <- c0_auc |>
  dplyr::group_by(modality) |>
  dplyr::summarise(
    n       = dplyr::n(),
    C0_gm   = gm(C0),
    AUC_gm  = gm(auc_trap),
    .groups = "drop"
  )

knitr::kable(
  c0_by_modality |>
    dplyr::transmute(
      `Modality`   = modality,
      `n`          = n,
      `GM C0 (ng/mL)` = signif(C0_gm, 3),
      `GM AUC0-336 (ng*h/mL)` = signif(AUC_gm, 4)
    ),
  caption = "Per-modality geometric mean C0 and AUC. Takama 2007 Discussion: differences between HD and PD distributions are not marked."
)

ggplot(c0_auc, aes(modality, C0, fill = modality)) +
  geom_violin(alpha = 0.3, draw_quantiles = c(0.25, 0.5, 0.75)) +
  geom_jitter(width = 0.15, alpha = 0.5) +
  labs(
    x       = "Dialysis modality",
    y       = "Initial serum concentration C0 (ng/mL)",
    caption = "Replicates Figure 3c of Takama 2007."
  ) +
  theme_minimal() +
  theme(legend.position = "none")
```

# PKNCA validation

The paper does not report a tabular NCA (Cmax / Tmax / AUC / half-life)
in the manuscript text or any table; per the Methods, C0 and AUC are
read off the POSTHOC estimates and reported in Figure 3 only as
distributions. We therefore compute NCA on the virtual cohort here to
verify (a) the typical-subject Cmax and AUC are consistent with the
paper's structural-parameter back-calculation
(`Cmax = dose / V1 + k0`; `AUC = dose / CL`) and (b) the apparent
terminal half-life is consistent with the literature value of
approximately 21-25 h for IV darbepoetin alfa (Aranesp prescribing
information).

```{r pknca}
# Time-zero guarantee: IV bolus simulation already produces a t = 0
# row, but defensively bind one in case future grid changes drop it.
# For IV bolus the t = 0 concentration is `dose / V1 + k0`; we let
# PKNCA back-extrapolate via lambda.z when needed, so seed any missing
# t = 0 with Cc = 0 (extravascular convention is harmless here because
# the explicit IV bolus row would dominate when present).
sim_nca <- sim |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, treatment)

sim_nca <- dplyr::bind_rows(
  sim_nca,
  sim_nca |> dplyr::distinct(id, treatment) |>
    dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, treatment, time, .keep_all = TRUE) |>
  dplyr::arrange(id, treatment, time)

conc_obj <- PKNCA::PKNCAconc(
  sim_nca, Cc ~ time | treatment + id,
  concu = "ng/mL",
  timeu = "h"
)

dose_df <- events |>
  dplyr::filter(evid == 1L) |>
  dplyr::transmute(id, time, amt, treatment)

dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | treatment + id,
                             doseu = "ug")

intervals <- data.frame(
  start       = 0,
  end         = Inf,
  cmax        = TRUE,
  tmax        = TRUE,
  aucinf.obs  = TRUE,
  auclast     = TRUE,
  half.life   = TRUE
)

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

## Comparison against published-structural back-calculation

Takama 2007 does not publish an NCA table. The closest deterministic
checks the manuscript supports are (a) the back-calculation
`Cmax = dose / V1 + k0` for an IV bolus and (b) the back-calculation
`AUC0-inf = dose / CL` from the steady-state mass balance. We compare
the simulated typical-subject NCA against these structural targets.

```{r nca_compare}
# Structural back-calculation at the per-modality median weight.
struct_back <- function(modality, wt, dose_ug) {
  cl    <- 0.0807 * (1 + 0.0195 * (wt - 54))
  perit <- if (modality == "HD") 0 else 1
  v1    <- 2.51   * (1 + 0.0163 * (wt - 54) + 0.170 * perit)
  k0    <- 0.167
  list(
    cmax       = dose_ug / v1 + k0,
    aucinf.obs = dose_ug / cl,
    tmax       = 0
  )
}

published <- tibble::tribble(
  ~treatment,    ~cmax,                                  ~tmax, ~aucinf.obs,
  "HD 60 ug",    struct_back("HD", 53.8, 60)$cmax,       0,     struct_back("HD", 53.8, 60)$aucinf.obs,
  "PD 60 ug",    struct_back("PD", 57.4, 60)$cmax,       0,     struct_back("PD", 57.4, 60)$aucinf.obs
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated     = nca_res,
  reference     = published,
  by            = "treatment",
  units         = c(cmax       = "ng/mL",
                    aucinf.obs = "ng*h/mL",
                    tmax       = "h"),
  tolerance_pct = 20
)

knitr::kable(
  cmp,
  caption = "Simulated stochastic-cohort NCA vs. structural back-calculation at the per-modality median weight. * marks rows that differ by >20 percent; these reflect the cohort's between-subject variability in CL and V1 (etalcl = 20.6 percent CV, etalvc = 21.8 percent CV) and are not parameter-tuning targets.",
  align   = c("l", "l", "r", "r", "r")
)
```

The geometric-mean simulated Cmax and AUC inevitably differ from the
median-weight back-calculation because the cohort sample drawn from the
truncated-normal weight distribution does not have median = 53.8 / 57.4
kg exactly, and because IIV on CL and V1 spreads the per-subject NCA
values. The paper's only quantitative claims are that (a) the structural
parameters listed in Table 4 are stable under 1000-replicate bootstrap
resampling (Table 4 right panel; the bootstrap mean of theta_CL was
0.0807 L/h, identical to the original-dataset estimate) and (b)
distributions of C0 and AUC are not markedly different between
weight tertiles or between HD and PD modalities (Figure 3 and
Discussion). The simulations above are consistent with both claims.

## Terminal half-life

Simulated half-life across the cohort:

```{r half-life}
hl_summary <- as.data.frame(nca_res$result) |>
  dplyr::filter(PPTESTCD == "half.life") |>
  dplyr::group_by(treatment) |>
  dplyr::summarise(
    n          = dplyr::n(),
    median_hl  = signif(median(PPORRES, na.rm = TRUE), 3),
    q25_hl     = signif(quantile(PPORRES, 0.25, na.rm = TRUE), 3),
    q75_hl     = signif(quantile(PPORRES, 0.75, na.rm = TRUE), 3),
    .groups    = "drop"
  )

knitr::kable(
  hl_summary |>
    dplyr::transmute(
      `Group`             = treatment,
      `n`                 = n,
      `Median t1/2 (h)`   = median_hl,
      `Q25 t1/2 (h)`      = q25_hl,
      `Q75 t1/2 (h)`      = q75_hl
    ),
  caption = "Simulated terminal half-life by group. Literature value for IV darbepoetin alfa in dialysis patients is approximately 21-25 h (Aranesp prescribing information, derived from Macdougall 1999); the simulated medians fall within or near this range."
)
```

# Assumptions and deviations

- **Weight distribution.** Original per-subject covariates are not
  publicly available. The virtual cohort draws WT from a truncated
  normal matched to the per-modality median and the Table 1 range; the
  cohort mean approximates but does not exactly equal the published
  median.
- **Sex, age, baseline lab covariates.** The paper screened sex,
  creatinine, RBC, WBC, platelets, and albumin via GAM and reported
  large LLD changes for WT and DIA only; sex, age, creatinine, RBC,
  WBC, platelets, and albumin were either not selected by GAM or were
  not retained at p < 0.005 in the backward-elimination step (Table 3).
  None of these enter the final model, so they are not represented in
  the cohort beyond their summary in `population$disease_state`.
- **Multiple-dose data.** PD subjects received darbepoetin alfa every
  1 or 2 weeks; the paper treats interoccasion variability as not
  required because the pharmacokinetics did not depend on duration of
  administration (citing reference [7]). The vignette simulation uses
  a single IV bolus to match the cleanest replicate of Figure 1; the
  packaged model itself accommodates multiple-dose events without
  modification.
- **Endogenous EPO encoding.** Takama 2007 Methods describe the
  baseline as `k0_j` for the j-th individual entering the predicted
  concentration `Cpm_ij`. The packaged model encodes this as an
  additive offset on the observation (`Cc = central / V1 + rbase`)
  rather than as a continuous endogenous-production rate input to the
  central compartment. The two formulations are mathematically
  equivalent for the linear two-compartment structure and the observed
  Cc(t) trajectory; the additive-offset form keeps the IV-bolus
  initial condition simple (`central(0) = 0`).
- **No tabulated NCA in the source.** Takama 2007 does not publish a
  Cmax / AUC / half-life table; per the Methods, C0 and AUC are read
  off POSTHOC estimates and reported as distributions in Figure 3
  only. The PKNCA validation block here compares simulated NCA to a
  structural back-calculation (`Cmax = dose / V1 + k0`,
  `AUC = dose / CL`) rather than to a published table.
- **Linear-deviation parameterisation.** The paper's
  `CL = theta_CL * [1 + theta_CL_WT * (WT - 54)]` form is preserved
  exactly. The form is positive across the cohort weight range
  (35.5-132 kg) but can in principle go negative for WT far below the
  cohort; simulation weights should therefore stay within the study
  range.
- **Diagonal IIV.** The paper reports each omega individually with no
  covariance matrix; the packaged model uses a diagonal $OMEGA
  (independent etas) consistent with the standard NONMEM convention.
