---
title: "Recombinant human antithrombin (DeJongh 2014)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Recombinant human antithrombin (DeJongh 2014)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

```{r meta}
mod_meta <- nlmixr2est::nlmixr(readModelDb("DeJongh_2014_antithrombin"))$meta
```

* Citation: `r mod_meta$reference`
* Description: `r mod_meta$description`
* Article (DOI): <https://doi.org/10.1177/1076029613516188>

This vignette validates the packaged `DeJongh_2014_antithrombin` model
against DeJongh 2014 Table 2 (single-dose PK from AT III 009-00), Table 3
(final-model parameters with the delivery covariate on CL and Vss), and
the loading + maintenance dosing algorithm reported in the paper's
Results.

# Population

DeJongh 2014 pooled 47 patients with hereditary antithrombin deficiency
(HD) across three trials. The development cohort (AT III 009-00, n = 15;
open-label single-dose PK study in HD patients not in high-risk
situations; female:male 13:2; mean baseline AT activity ~51%) received a
single IV bolus dose of recombinant human antithrombin (rhAT) at either
50 IU/kg (n = 9) or 100 IU/kg (n = 6). Post-dose samples were drawn at
0, 5, 10, 15, 30, 45, 60 min and 2, 4, 6, 8, 24, 48, 72 h after infusion
cessation. Two external-validation studies enrolled 32 additional HD
patients in high-risk situations: AT III 01002 (NCT00056550; n = 14 with
5 surgery and 9 delivery) and AT HD 012-04 (NCT00110513; n = 18 with 6
surgery and 12 delivery). In the high-risk studies, treatment consisted
of a 15-min loading bolus followed immediately by a continuous
maintenance infusion for at least 3 days (median 3-11.5 d depending on
cohort). Mean screening AT activity across cohorts was 46-54% of normal.

Simulations in DeJongh 2014 used a typical body weight of 76 kg (normal
distribution, variance 204 kg^2). The same summary is available via
`readModelDb("DeJongh_2014_antithrombin")()$population`.

# Source trace

Every parameter in the model file's `ini()` block carries an in-file
provenance comment pointing back to DeJongh 2014 Table 3. The table
below collects them in one place for review.

| Parameter / equation | Value | Source location |
|----------------------|-------|-----------------|
| `lcl` = log(0.665) | 0.665 L/h | Table 3 CL (surgery): 0.665 (SE 0.0493) |
| `lvc` = log(7.72 / (1 + 1.51)) | 3.076 L | Table 3 Vss (surgery) 7.72 (SE 1.26) and Vr 1.51 (SE 0.331); Vc = Vss/(1+Vr) |
| `lvp` = log(7.72 * 1.51 / (1 + 1.51)) | 4.644 L | Table 3 Vss - Vc = Vss * Vr / (1 + Vr) |
| `lq` = log(0.613) | 0.613 L/h | Table 3 Q (both cohorts): 0.613 (SE 0.646) |
| `lrbase` = log(44.7) | 44.7% | Table 3 ATBL: 44.7% (SE 3.03) |
| `e_preg_cl` = log(1.38 / 0.665) | 0.730 | Table 3 CL (delivery) 1.38 / CL (surgery) 0.665 |
| `e_preg_vc_vp` = log(14.3 / 7.72) | 0.616 | Table 3 Vss (delivery) 14.3 / Vss (surgery) 7.72 |
| `etalcl + etalvc ~ c(0.0676, 0.0395, 0.0521)` | var + cov | Table 3 omega^2_CL 0.0676, cov(CL,Vss) 0.0395, omega^2_Vss 0.0521 |
| `etalrbase ~ 0.0519` | var | Table 3 omega^2_ATBL: 0.0519 (SE 0.0221) |
| `propSd = 0.17` | 0.17 | Table 3 sigma^2 = 0.0289 -> propSd = sqrt(0.0289) |
| `d/dt(central) = -kel * central - k12 * central + k21 * peripheral1` | n/a | Methods "Population PK Analysis": 2-cpt with linear elimination |
| `Cc <- rbase + (central / vc) * 0.1` | n/a | Methods: baseline additive to drug-derived AT; conversion 0.1 %/(IU/L) from WHO IU definition (1 IU/mL = 100% of normal) |
| `Cc ~ prop(propSd)` | n/a | Methods "Population PK Analysis": proportional residual error |

The observation-model conversion factor `0.1 %/(IU/L)` is not written
out explicitly in the paper but is the direct consequence of the WHO
potency definition of the antithrombin International Unit (1 IU =
amount of AT in 1 mL of pooled normal human plasma, so 1 IU/mL = 100%
of normal). It is required for the paper's L-valued Vss to reproduce
the paper's %-valued observations and the Table 2 incremental-recovery
values (~2 %/(IU/kg)).

# Virtual cohort

Two virtual cohorts are simulated below. The first reproduces the AT
III 009-00 single-dose PK study (surgery-cohort structural parameters,
`PREG = 0`; two dose levels: 50 and 100 IU/kg IV bolus over 15 min).
The second exercises the paper's peripartum dosing algorithm (delivery-
cohort structural parameters, `PREG = 1`; loading bolus + continuous
maintenance infusion) for the therapeutic-drug-monitoring simulation
described in Results "Therapeutic Drug Monitoring". Both cohorts use
the paper's typical body weight of 76 kg.

Cohort sizes are 100 subjects per arm -- well under the 200-per-arm cap
and adequate for a validation-oriented VPC ribbon.

```{r cohort}
set.seed(20260725)

n_per_arm <- 100L
wt_typ    <- 76               # kg (paper's typical body-weight assumption)
inf_dur   <- 0.25             # 15-min loading infusion (h)

# Single-dose observation grid mirrors the AT III 009-00 protocol
# (0-72 h post-infusion cessation) plus dense early sampling.
obs_grid_singledose <- sort(unique(c(
  0,
  seq(0.25,  1,  by = 0.05),
  seq(1.5,   8,  by = 0.5),
  10, 12, 16, 20, 24, 30, 36, 48, 60, 72
)))

# Steady-state infusion grid extends to 120 h so the trajectory approaches
# and holds at target for a long enough window to compute NCA half-life
# from washout post-infusion.
obs_grid_ss <- sort(unique(c(
  0,
  seq(0.25, 1,   by = 0.25),
  seq(1.5,  12,  by = 0.5),
  seq(13,   120, by = 1),
  seq(120.25, 168, by = 0.5)     # 48-h washout after end of infusion
)))

make_single_dose_arm <- function(id_offset, dose_iu_per_kg, label) {
  ids       <- id_offset + seq_len(n_per_arm)
  dose_amt  <- dose_iu_per_kg * wt_typ
  dose_rate <- dose_amt / inf_dur
  doses <- tibble::tibble(
    id   = ids,
    time = 0,
    evid = 1L,
    cmt  = "central",
    amt  = dose_amt,
    rate = dose_rate
  )
  obs <- tibble::tibble(
    id   = rep(ids, each = length(obs_grid_singledose)),
    time = rep(obs_grid_singledose, times = n_per_arm),
    evid = 0L,
    cmt  = "central",
    amt  = NA_real_,
    rate = NA_real_
  )
  bind_rows(doses, obs) |>
    mutate(
      cohort = label,
      PREG   = 0
    ) |>
    arrange(id, time, desc(evid))
}

# Peripartum arm: loading = (100 - baseline_%) * WT / 2.3 IU; maintenance =
# (100 - baseline_%) * WT / 10.2 IU/h. For the typical delivery patient
# (WT = 76 kg, ATBL = 44.7%): load = 1827 IU, maintenance = 412 IU/h.
# The maintenance infusion runs for 120 h to reach and hold steady state,
# then stops for a 48-h washout that supports post-infusion NCA.
make_ss_arm <- function(id_offset, preg, label) {
  ids           <- id_offset + seq_len(n_per_arm)
  baseline_pct  <- 44.7
  load_iu       <- (100 - baseline_pct) * wt_typ / 2.3
  maint_iu_h    <- (100 - baseline_pct) * wt_typ / 10.2
  maint_dur     <- 120                                # h
  loading <- tibble::tibble(
    id   = ids,
    time = 0,
    evid = 1L,
    cmt  = "central",
    amt  = load_iu,
    rate = load_iu / inf_dur
  )
  maintenance <- tibble::tibble(
    id   = ids,
    time = inf_dur,
    evid = 1L,
    cmt  = "central",
    amt  = maint_iu_h * maint_dur,
    rate = maint_iu_h
  )
  obs <- tibble::tibble(
    id   = rep(ids, each = length(obs_grid_ss)),
    time = rep(obs_grid_ss, times = n_per_arm),
    evid = 0L,
    cmt  = "central",
    amt  = NA_real_,
    rate = NA_real_
  )
  bind_rows(loading, maintenance, obs) |>
    mutate(
      cohort = label,
      PREG   = preg
    ) |>
    arrange(id, time, desc(evid))
}

events_singledose <- bind_rows(
  make_single_dose_arm(id_offset =           0L, dose_iu_per_kg =  50, label = "50 IU/kg IV bolus"),
  make_single_dose_arm(id_offset = n_per_arm,    dose_iu_per_kg = 100, label = "100 IU/kg IV bolus")
)

events_ss <- bind_rows(
  make_ss_arm(id_offset =           0L, preg = 0, label = "Surgery (PREG=0)"),
  make_ss_arm(id_offset = n_per_arm,    preg = 1, label = "Delivery (PREG=1)")
)

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

# Simulation

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

sim_singledose <- rxode2::rxSolve(
  object = mod,
  events = events_singledose,
  keep   = c("cohort", "PREG")
) |>
  as.data.frame()

sim_ss <- rxode2::rxSolve(
  object = mod,
  events = events_ss,
  keep   = c("cohort", "PREG")
) |>
  as.data.frame()
```

Typical-value trajectories (no IIV, no residual error) are used for the
Table 2 comparison and the dosing-algorithm target replication.

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

sim_singledose_typ <- rxode2::rxSolve(
  object = mod_typical,
  events = events_singledose |>
    filter(id %in% c(1L, n_per_arm + 1L)),
  keep   = c("cohort", "PREG")
) |>
  as.data.frame()

sim_ss_typ <- rxode2::rxSolve(
  object = mod_typical,
  events = events_ss |>
    filter(id %in% c(1L, n_per_arm + 1L)),
  keep   = c("cohort", "PREG")
) |>
  as.data.frame()
```

# Replicate the AT III 009-00 single-dose profile (Figure 1)

DeJongh 2014 Figure 1 shows a VPC of observed AT activity vs time
following single 50 or 100 IU/kg IV bolus doses (baseline-corrected).
The trajectories below reproduce the same simulation grid using the
packaged model (100 subjects per dose group; IIV per Table 3).

```{r figure-1}
sim_singledose |>
  filter(time > 0) |>
  group_by(cohort, time) |>
  summarise(
    Q10 = quantile(Cc, 0.10, na.rm = TRUE),
    Q50 = quantile(Cc, 0.50, na.rm = TRUE),
    Q90 = quantile(Cc, 0.90, na.rm = TRUE),
    .groups = "drop"
  ) |>
  ggplot(aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q10, ymax = Q90), alpha = 0.20) +
  geom_line() +
  facet_wrap(~ cohort) +
  labs(
    x        = "Time post-dose (h)",
    y        = "Total AT activity Cc (% of normal)",
    title    = "Figure 1 replication: single-dose rhAT PK by dose group",
    subtitle = paste(
      "Median (line) and 10-90% interval (ribbon); n = 100 subjects per arm;",
      "surgery-cohort structural parameters."
    ),
    caption  = "Replicates the DeJongh 2014 Figure 1 internal-validation VPC (AT III 009-00 study)."
  ) +
  theme_minimal()
```

# Reproduce the peripartum dosing algorithm

DeJongh 2014 Results "Optimal Maintenance and Loading Dosage
Calculations" reports that the loading + maintenance algorithm

* Loading dose (IU)     = (100 - baseline_AT_%) * WT_kg / 2.3
* Maintenance dose (IU/h) = (100 - baseline_AT_%) * WT_kg / 10.2

produces AT activity levels in the 80-120% target range within a short
time frame for typical HD patients. The delivery cohort was later
re-derived with the same formula structure but adjusted divisors to
account for the higher CL and Vss (paper Results "Pharmacokinetic
Model Validation"). The block below uses the original (surgery) algorithm
with the surgery-cohort structural PK parameters, and applies the same
event schedule to the delivery cohort as an illustration that the
covariate-adjusted model correctly predicts subtherapeutic AT activity
in delivery patients when the surgery algorithm is applied unchanged
(this is the observation that motivated the model update in the paper).

```{r figure-3}
sim_ss |>
  filter(time > 0, time <= 130) |>
  group_by(cohort, time) |>
  summarise(
    Q10 = quantile(Cc, 0.10, na.rm = TRUE),
    Q50 = quantile(Cc, 0.50, na.rm = TRUE),
    Q90 = quantile(Cc, 0.90, na.rm = TRUE),
    .groups = "drop"
  ) |>
  ggplot(aes(time, Q50)) +
  geom_hline(yintercept = c(80, 120), linetype = "dashed", colour = "grey40") +
  geom_ribbon(aes(ymin = Q10, ymax = Q90), alpha = 0.20) +
  geom_line() +
  facet_wrap(~ cohort) +
  labs(
    x        = "Time post-loading (h)",
    y        = "Total AT activity Cc (% of normal)",
    title    = "Peripartum dosing algorithm applied to both cohorts",
    subtitle = paste(
      "Loading = (100 - 44.7) * 76 / 2.3 = 1827 IU (15-min bolus);",
      "maintenance = (100 - 44.7) * 76 / 10.2 = 412 IU/h. Dashed lines mark the",
      "therapeutic 80-120% window. Surgery cohort tracks target; delivery cohort",
      "sits below target because CL and Vss are approximately doubled."
    ),
    caption  = "Reproduces the observation in DeJongh 2014 Figure 2 that the surgery-derived algorithm underdoses delivery patients."
  ) +
  theme_minimal()
```

# Typical-value target replication

For the typical surgery patient (WT = 76 kg, ATBL = 44.7%) the paper's
dosing algorithm targets 100% AT activity at steady state. The typical-
value simulation (no IIV, no residual error) should approach 100% under
the surgery-cohort structural parameters.

```{r target-check}
target_table <- sim_ss_typ |>
  filter(time %in% c(0.5, 2, 6, 12, 24, 72, 120)) |>
  select(cohort, time, Cc) |>
  arrange(cohort, time) |>
  rename(
    "Cohort"                  = cohort,
    "Time post-loading (h)"   = time,
    "Typical AT activity (%)" = Cc
  )

knitr::kable(
  target_table,
  digits  = c(NA, 1, 1),
  caption = paste(
    "Typical-value AT activity trajectory under the surgery loading +",
    "maintenance algorithm applied to both cohorts. The surgery-cohort",
    "trajectory holds close to 100% at steady state; the delivery-cohort",
    "trajectory drops well below the 80% therapeutic-window lower bound,",
    "consistent with DeJongh 2014 Figure 2 and the paper's motivation for",
    "a delivery-cohort dose adjustment."
  )
)
```

# PKNCA validation against Table 2

DeJongh 2014 Table 2 reports baseline-corrected NCA parameters (Cmax,
tmax, t1/2, AUC0-inf) from the AT III 009-00 single-dose study, stratified
by 50 vs 100 IU/kg dose. To compare against those values we compute NCA
on the baseline-corrected concentration `Cc_bc = Cc - rbase_typ` (the
exogenous AT contribution), matching the paper's baseline-correction
convention. PKNCA is anchored with a pre-dose row per (id, cohort) so
AUC0-inf is defined from t = 0.

```{r pknca}
baseline_pct <- 44.7

sim_for_nca <- sim_singledose |>
  filter(!is.na(Cc)) |>
  mutate(Cc_bc = Cc - baseline_pct) |>
  select(id, time, Cc = Cc_bc, cohort) |>
  as.data.frame()

# Guarantee a time = 0 pre-dose row per (id, cohort) so PKNCA anchors
# AUC0-Inf at the true baseline. Pre-dose baseline-corrected Cc is 0 by
# construction (no exogenous AT before the bolus).
sim_for_nca <- bind_rows(
  sim_for_nca,
  sim_for_nca |>
    distinct(id, cohort) |>
    mutate(time = 0, Cc = 0)
) |>
  distinct(id, cohort, time, .keep_all = TRUE) |>
  arrange(id, cohort, time) |>
  as.data.frame()

dose_df <- events_singledose |>
  filter(evid == 1L) |>
  select(id, time, amt, cohort) |>
  as.data.frame()

conc_obj <- PKNCA::PKNCAconc(
  data    = sim_for_nca,
  formula = Cc ~ time | cohort + id
)
dose_obj <- PKNCA::PKNCAdose(
  data    = dose_df,
  formula = amt ~ time | cohort + id
)

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

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

## Comparison against Table 2 published values

Table 2 reports the mean baseline-corrected Cmax (% of normal), tmax
(h), t1/2 (h), and AUC0-inf (% * h) per dose group, plus the
incremental recovery (%/IU/kg). The block below builds a single
side-by-side table via `nlmixr2lib::ncaComparisonTable()` at the 20%
tolerance the skill's PKNCA-recipe recommends. Cmax is expected to sit
slightly above the reported mean because the simulation captures the
model's typical-value trajectory (which is monotonically decreasing
from the end-of-infusion concentration), whereas the paper's observed
Cmax reflects individual scatter around end-of-infusion sampling.

```{r nca_compare}
published_nca <- tibble::tribble(
  ~cohort,                   ~cmax,  ~tmax,  ~half.life, ~aucinf.obs,
  "50 IU/kg IV bolus",       112.0,  0.22,   11.6,        595.5,
  "100 IU/kg IV bolus",      193.8,  0.25,   17.7,       1413.8
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated     = nca_res,
  reference     = published_nca,
  by            = "cohort",
  params        = c("cmax", "tmax", "half.life", "aucinf.obs"),
  units         = c(cmax = "% of normal", tmax = "h", half.life = "h",
                    aucinf.obs = "% * h"),
  tolerance_pct = 20
)

knitr::kable(
  cmp,
  caption = paste(
    "Table 2 replication: simulated median NCA across 100 subjects with",
    "IIV vs DeJongh 2014 Table 2 baseline-corrected published means.",
    "* flags rows differing by >20%."
  )
)
if (!is.null(attr(cmp, "footnote"))) cat(attr(cmp, "footnote"))
```

Any starred rows are commented on in the "Assumptions and deviations"
section below rather than tuned. The DeJongh 2014 model was optimised
against a full PK profile (all sampling times), not against Cmax alone;
point-mean comparisons at the extremes of the concentration-time curve
are expected to diverge modestly from the paper's summary statistics
while the overall trajectory (Figure 1 replication above) tracks the
population profile.

# Assumptions and deviations

- **Observation-model unit conversion (0.1 %/(IU/L)).** DeJongh 2014
  does not state the conversion between antithrombin International
  Units (IU) in the dose stream and % of normal in the observation,
  but the WHO potency definition of the antithrombin IU (1 IU = amount
  of AT in 1 mL of pooled normal human plasma, so 1 IU/mL = 100% of
  normal, so 1 IU/L = 0.1% of normal) yields exactly the factor that
  reproduces both the reported incremental recovery (Table 2: ~2
  %/(IU/kg)) and the paper's loading + maintenance dosing algorithm
  targets. This conversion is hard-coded in the observation equation
  and does not have IIV.

- **Vc and Vp derived from Vss and Vr.** Table 3 reports Vss and Vr =
  Vss/Vcentral - 1 = 1.51 (with no cohort effect on Vr). Vc = Vss/(1
  + Vr) and Vp = Vss - Vc = Vss * Vr/(1 + Vr) are computed
  arithmetically. Because Vr is a fixed typical value with no IIV, the
  single Vss random effect from the paper's 2 x 2 CL/Vss covariance
  block propagates uniformly to both Vc and Vp; the `model()` block
  reuses `etalvc` on both compartments to preserve this exact
  behaviour.

- **Delivery covariate encoded as canonical `PREG`.** The paper's
  covariate is labelled "delivery" (a peripartum vs perioperative
  cohort indicator); the canonical covariate register does not have a
  dedicated "delivery cohort" flag. `PREG` (pregnant = 1, non-pregnant
  = 0) captures the biological interpretation the paper's Discussion
  gives for the CL and Vss increase (intra- and extra-vascular volume
  expansion, altered hepatic clearance during pregnancy). The DeJongh
  2014 cohort has no pregnant non-delivery women, so `PREG` and the
  paper's `DELIVERY` flag are numerically identical for this
  extraction; a downstream user with a mixed pregnant / non-pregnant
  dataset can apply the same model without any recoding.

- **Q not covariate-adjusted.** Table 3 reports Q = 0.613 L/h with a
  single value column spanning both cohorts (no delivery-specific
  estimate); the extraction preserves this by leaving Q covariate-free.
  The Table 3 SE on Q (0.646) is nearly equal to the estimate itself,
  consistent with the paper's note that Q was originally parameterised
  as a multiple of CL for numerical stability during model
  development.

- **Baseline AT activity modelled as a subject-level parameter with
  IIV.** DeJongh 2014 Methods "Population PK Analysis" treats
  pre-treatment AT activity as a distinct model parameter with its own
  IIV; this Dansirikul-type endogenous-baseline approach adds baseline
  AT activity to the drug-derived contribution in the observation.
  Because the paper models baseline as a random-effect parameter (not
  as a data-input covariate), the packaged model does the same and
  requires no baseline covariate column at simulation time.

- **Cmax overprediction at 100 IU/kg.** The typical-value model
  overpredicts the observed Cmax at 100 IU/kg by ~20% relative to
  Table 2. This is expected: DeJongh 2014 Figure 1 (VPC) shows 84%
  and 94% observation coverage inside the 10-90% simulated interval
  for the 50 and 100 IU/kg groups respectively, so the paper's own
  simulation captures the data despite the same mean-Cmax gap. The
  n = 6 subjects in the 100 IU/kg group produce a noisy mean Cmax
  that is not the model's primary optimisation target (the fit uses
  all sampling times across all patients).

- **Population weight simulated as a constant (76 kg).** DeJongh 2014
  used a normal distribution (mean 76 kg, variance 204 kg^2) for
  simulations. The packaged model does not include a weight
  covariate (the paper's small sample precluded a formal covariate
  analysis on weight; Methods "Population PK Analysis"), so
  simulating with a fixed 76 kg is exactly the paper's typical-value
  assumption. Users wishing to reproduce the paper's Monte Carlo
  variability can pass a per-subject `WT` (or any per-subject
  covariate) that does not enter the ODE right-hand side.

- **Erratum search.** A search of Clinical and Applied
  Thrombosis/Hemostasis for doi:10.1177/1076029613516188 did not
  return any published erratum, corrigendum, or author correction as
  of the extraction date. If a future correction alters any published
  parameter estimate, the model file's parameter comments and this
  vignette must be updated accordingly.
