---
title: "Antithrombin (Moffett 2017)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Antithrombin (Moffett 2017)}
  %\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("Moffett_2017_antithrombin"))$meta
```

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

This vignette validates the packaged `Moffett_2017_antithrombin` model
against Moffett 2017 Table 3 (final-model equations), Table 4 (bootstrap
point estimates), and Table 5 (simulated pharmacokinetic parameters and
predicted concentrations for the median cohort patient at 1, 12, and 24 h
post-infusion).

# Population

Moffett 2017 studied 184 paediatric inpatients (birth to <19 y, median
age 0.35 y, median weight 5.2 kg, 53.3% female) at a single U.S. centre
who received human plasma-derived antithrombin (AT) concentrate
(Thrombate, Grifols) for consumptive AT depletion. The dominant clinical
context was concurrent continuous-infusion unfractionated heparin (UFH)
therapy: 87.5% of patients received concurrent UFH at a mean per-kg dose
of 34.1 +/- 22.7 units/kg/h; 8.2% were treated for coagulopathy and 4.3%
for concurrent enoxaparin therapy. Patients on mechanical circulatory
support during AT administration or monitoring were excluded. Baseline
demographics and laboratory values are in Table 1 and Table 2 of the
source paper. The same information is available programmatically via
`readModelDb("Moffett_2017_antithrombin")()$population`.

# Source trace

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

| Parameter / equation | Value | Source location |
|----------------------|-------|-----------------|
| `lcl` (log CL_ref at WT=70, UFH=0)   | log(0.917)   | Table 3 & Table 4 CL row: 0.917 dL/h (RSE 16.8%, 95% CI 0.615-1.219) |
| `lvc` (log VD_ref at WT=70, BASE=60) | log(67.9)    | Table 3 & Table 4 VD row: 67.9 dL (RSE 4.9%, 95% CI 61.41-74.39) |
| `e_wt_cl` (fixed at 0.75)            | 0.75         | Methods "Model building" and Table 3 (theoretical allometric exponent, fixed) |
| `e_wt_vc` (fixed at 1.00)            | 1.00         | Methods "Model building" and Table 3 (theoretical allometric exponent, fixed) |
| `e_ufh_cl` (linear-additive slope)   | 0.129        | Table 3 & Table 4 UFH row: 0.129 dL/h per (UFH/173) unit (RSE 15.5%) |
| `e_at_bl_vc` (power exponent on VD)  | -0.389       | Table 3 & Table 4 BASE row: -0.389 (RSE 37.5%) |
| `etalcl ~ 0.20353` (IIV variance)    | 47.5% CV     | Table 4 omega1: 47.5% CV (RSE 26.5%); variance = log(0.475^2 + 1) |
| `etalvc ~ 0.07649` (IIV variance)    | 28.2% CV     | Table 4 omega2: 28.2% CV (RSE 31.2%); variance = log(0.282^2 + 1) |
| `expSd <- 0.11056` (log-scale SD)    | 11.1% CV     | Table 4 exponential-error: 11.1% CV (RSE 26.9%); expSd = sqrt(log(0.111^2 + 1)) |
| `tvcl = exp(lcl)*(WT/70)^0.75 + e_ufh_cl*(UFH/173)` | n/a | Table 3 final-model CL equation (additive weight + UFH structure) |
| `tvvc = exp(lvc)*(WT/70)^1*(BASE/60)^-0.389`        | n/a | Table 3 final-model VD equation (multiplicative weight * BASE power) |
| `cl = tvcl * exp(etalcl)`            | n/a | Methods "Model building": exponential IIV on CL |
| `vc = tvvc * exp(etalvc)`            | n/a | Methods "Model building": exponential IIV on VD |
| `d/dt(central) <- -kel * central`    | n/a | Results "Pharmacokinetic modelling": one-compartment model |
| `Cc <- central/vc + AT_BL_UDL`       | n/a | Methods "Data handling": Dansirikul B3 method for endogenous baseline |
| `Cc ~ lnorm(expSd)`                  | n/a | Results "Pharmacokinetic modelling": exponential residual error model |

# Virtual cohort

To reproduce Moffett 2017 Table 5 the vignette simulates the paper's
"median patient" (WT = 5.2 kg, baseline AT activity 60 units/dL) under
two concurrent-UFH conditions: (a) the cohort-mean UFH infusion of 34
units/kg/h absolute (= 176.8 units/h for a 5.2 kg patient) and (b) no
concurrent UFH. Each subject receives a single 60 units/kg IV dose of
AT concentrate (= 312 units) infused over 15 minutes, matching the
Table 5 simulation dose. 200 subjects per arm are simulated with IIV
per Table 4 (47.5% CV on CL, 28.2% CV on VD) plus exponential residual
error.

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

n_per_arm <- 200L
wt_med    <- 5.2   # kg, cohort median (Moffett 2017 Table 1)
base_med  <- 60    # units/dL, reference cohort baseline (Table 4)
dose_amt  <- 60 * wt_med           # 60 units/kg * 5.2 kg = 312 units
inf_dur   <- 0.25                  # 15-min infusion, in hours
dose_rate <- dose_amt / inf_dur    # units/h during the infusion

obs_grid <- sort(unique(c(
  0.0, 0.25,
  seq(0.5,  4.0, by = 0.25),
  seq(4.5, 12.0, by = 0.5),
  seq(13.0, 72.0, by = 1.0)
)))

make_arm <- function(id_offset, ufh_absolute_uh, label) {
  ids <- id_offset + seq_len(n_per_arm)
  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)),
    time = rep(obs_grid, times = n_per_arm),
    evid = 0L,
    cmt  = "central",
    amt  = NA_real_,
    rate = NA_real_
  )
  bind_rows(doses, obs) |>
    mutate(
      arm         = label,
      WT          = wt_med,
      DOSE_UFH_UH = ufh_absolute_uh,
      AT_BL_UDL   = base_med
    ) |>
    arrange(id, time, desc(evid))
}

events <- bind_rows(
  make_arm(id_offset = 0L,           ufh_absolute_uh = 34 * wt_med, label = "UFH 34 units/kg/h"),
  make_arm(id_offset = n_per_arm,    ufh_absolute_uh = 0,           label = "No UFH")
)

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

# Simulation

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

sim <- rxode2::rxSolve(
  object = mod,
  events = events,
  keep   = c("arm")
) |>
  as.data.frame()
```

Deterministic typical-value trajectories (no IIV, no residual error) are
used to reproduce Moffett 2017 Table 5's median-patient predicted
concentrations at 1, 12, and 24 h post-infusion.

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

sim_typical <- rxode2::rxSolve(
  object = mod_typical,
  events = events |> filter(id %in% c(1L, n_per_arm + 1L)),  # one subject per arm
  keep   = c("arm")
) |>
  as.data.frame()
```

# Replicate Table 5 -- predicted total AT activity at 1, 12, and 24 h

Moffett 2017 Table 5 reports median-patient predicted post-infusion AT
activity levels (dose-derived AT plus endogenous baseline, per the
Dansirikul B3 handling) for the two UFH conditions. The block below
extracts the packaged model's typical-value predictions at the same
post-infusion time-points.

```{r table5-compare}
# Time reported in Table 5 is "hours post-infusion"; the 15-min infusion
# starts at time 0 and ends at time 0.25 h, so "1 h post-infusion" is
# t = 0.25 + 1 = 1.25 h in the simulation clock. Match the three Table 5
# time-points.
table5_times <- c("1 h" = 0.25 + 1, "12 h" = 0.25 + 12, "24 h" = 0.25 + 24)

sim_at_table5 <- sim_typical |>
  filter(time %in% unname(table5_times)) |>
  mutate(hours_post_infusion = names(table5_times)[match(time, unname(table5_times))]) |>
  select(arm, hours_post_infusion, Cc_sim = Cc)

published <- tibble::tribble(
  ~arm,                     ~hours_post_infusion, ~Cc_paper_mean, ~Cc_paper_sd,
  "UFH 34 units/kg/h",      "1 h",                120,            19,
  "UFH 34 units/kg/h",      "12 h",               92,             14,
  "UFH 34 units/kg/h",      "24 h",               78,             13,
  "No UFH",                 "1 h",                120,            20,
  "No UFH",                 "12 h",               105,            15,
  "No UFH",                 "24 h",               92,             14
)

table5_compare <- sim_at_table5 |>
  left_join(published, by = c("arm", "hours_post_infusion")) |>
  mutate(
    hours_post_infusion = factor(hours_post_infusion, levels = c("1 h", "12 h", "24 h")),
    rel_diff_pct = round(100 * (Cc_sim - Cc_paper_mean) / Cc_paper_mean, 1)
  ) |>
  arrange(arm, hours_post_infusion) |>
  rename(
    "Arm"                             = arm,
    "Time post-infusion"              = hours_post_infusion,
    "Simulated typical-value (units/dL)" = Cc_sim,
    "Table 5 mean (units/dL)"         = Cc_paper_mean,
    "Table 5 SD (units/dL)"           = Cc_paper_sd,
    "Relative difference (%)"         = rel_diff_pct
  )

knitr::kable(
  table5_compare,
  digits  = c(NA, NA, 1, 0, 0, 1),
  caption = paste(
    "Table 5 replication: packaged typical-value predictions vs Moffett 2017",
    "Table 5 published means for the median patient (WT = 5.2 kg,",
    "baseline AT = 60 units/dL) after a 60 units/kg 15-min AT infusion."
  )
)
```

Every simulated typical-value point sits inside the 1-SD band of the
published Table 5 mean. The +2 units/dL residual on the "1 h" row (both
arms) is the difference between the paper's residual-error-laden
simulation mean and the deterministic typical-value point; the paper
reports Cc for a 60 units/kg dose that the paper writes as producing
"an hour after initiation of infusion" 120 units/dL (Results paragraph
5, quoted verbatim), which the typical-value model reproduces to
~1 unit/dL.

# Trajectories

```{r vpc}
sim |>
  filter(time > 0) |>
  group_by(arm, time) |>
  summarise(
    Q05 = quantile(Cc, 0.05, na.rm = TRUE),
    Q50 = quantile(Cc, 0.50, na.rm = TRUE),
    Q95 = quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  ggplot(aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.20) +
  geom_line() +
  facet_wrap(~ arm) +
  labs(
    x        = "Time (h post-dose)",
    y        = "Total AT activity Cc (units/dL)",
    title    = "Simulated total AT activity by concurrent-UFH arm",
    subtitle = paste(
      "Median (line) and 5-95% interval (ribbon); n = 200 subjects per arm,",
      "WT = 5.2 kg, baseline AT = 60 units/dL. Includes IIV + exponential residual error."
    )
  ) +
  theme_minimal()
```

# PKNCA validation of the drug-derived component

Moffett 2017 Table 5 reports the pharmacokinetic half-life, CL, and VD
computed for the median patient in a 10 000-replication simulation with
IIV. The published `Cc` values include the endogenous baseline
(Dansirikul B3 handling), so a naive PKNCA half-life on the total
observed AT would asymptote to the baseline and give a nonsensical
lambda-z. To validate against Table 5 we extract the drug-derived
concentration `Cc_drug = Cc - AT_BL_UDL` (the exogenous AT contribution
alone) and run PKNCA on that.

```{r pknca}
sim_for_nca <- sim |>
  filter(!is.na(Cc)) |>
  mutate(Cc_drug = Cc - AT_BL_UDL) |>
  select(id, time, Cc = Cc_drug, arm) |>
  as.data.frame()

# Guarantee a time = 0 pre-dose row for every (id, arm) so PKNCA anchors
# AUC0-Inf at the true baseline. Pre-dose Cc_drug = 0 (no exogenous AT).
sim_for_nca <- bind_rows(
  sim_for_nca,
  sim_for_nca |>
    distinct(id, arm) |>
    mutate(time = 0, Cc = 0)
) |>
  distinct(id, arm, time, .keep_all = TRUE) |>
  arrange(id, arm, time) |>
  as.data.frame()

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

conc_obj <- PKNCA::PKNCAconc(
  data    = sim_for_nca,
  formula = Cc ~ time | arm + id,
  concu   = "units/dL",
  timeu   = "hr"
)
dose_obj <- PKNCA::PKNCAdose(
  data    = dose_df,
  formula = amt ~ time | arm + id,
  doseu   = "units"
)

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 5 published values

```{r nca_compare}
# Table 5 reports CL and VD in mL/kg/h and mL/kg. For the median patient
# (WT = 5.2 kg) convert to the model's absolute units (dL/h and dL) to
# compare directly: absolute_CL_dLh = CL_ml_kg_h * WT / 100 (1 dL = 100 mL).
# But the direct NCA outputs are Cmax (units/dL) and t-half (h), which
# are cohort-mean-with-IIV values in Table 5 -- so compare those directly.

published_nca <- tibble::tribble(
  ~arm,                  ~half.life,
  "UFH 34 units/kg/h",   15.5,
  "No UFH",              31.6
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated     = nca_res,
  reference     = published_nca,
  by            = "arm",
  params        = "half.life",
  units         = c(half.life = "h"),
  tolerance_pct = 20
)

knitr::kable(
  cmp,
  caption = paste(
    "Table 5 half-life comparison: simulated median across 200 subjects with",
    "IIV vs Moffett 2017 Table 5 half-life. * flags rows differing by >20%."
  ),
  align   = c("l", "l", "r", "r", "r")
)
cat(attr(cmp, "footnote"))
```

The simulated half-life sits well within the 20% tolerance for both
UFH conditions and matches the published paper's simulation
(Moffett 2017 Table 5 reports 15.5 +/- 9.3 h with UFH and 31.6 +/- 19.7
h without). The mean-with-IIV half-life for the UFH arm is
approximately ln(2) * E[VD/CL] under log-normal disposition, which under
the packaged model's IIV variances evaluates to ~15.3 h, matching the
paper's reported 15.5 h to ~1%.

Because Table 5 does not tabulate Cmax / AUCinf per arm as free
comparators, only the half-life comparison is included here. The
qualitative CL and VD comparison (paper Table 5: 5.7 +/- 2.8 mL/kg/h
with UFH vs 2.8 +/- 1.4 mL/kg/h without; VD ~101 mL/kg in both arms)
is exercised implicitly by the half-life match and the Table 5
concentration-point replication above.

# Assumptions and deviations

- **Endogenous-baseline handling.** The observed AT activity is the sum
  of dose-derived AT and the endogenous baseline
  (`Cc = central / vc + AT_BL_UDL`) per the Dansirikul et al. 2008 B3
  method described in Moffett 2017 Methods "Data handling". The
  packaged model applies the baseline as an additive constant in the
  observation equation; sensitivity to the assumption that endogenous
  AT is time-invariant over the observation window (a simplification
  the source paper also makes) has not been re-tested here.

- **Allometric exponents are fixed at theoretical values (0.75 on CL,
  1.00 on VD).** Moffett 2017 Methods "Model building" states these
  exponents were the theoretical values used for the primary model; the
  paper does not report per-parameter uncertainty on the exponents,
  consistent with the fixed status. The `ini()` block wraps both
  exponents in `fixed()` to preserve this provenance.

- **UFH covariate is time-varying at the record level in production
  data; held constant per subject here.** The source paper's dataset
  carried UFH as a per-observation infusion-rate covariate; the
  vignette's simulation uses a single-arm-constant UFH rate (34
  units/kg/h absolute for the UFH arm, 0 for the no-UFH arm) to match
  Moffett 2017 Table 5's median-patient simulation, which the paper
  also treats as a constant per arm.

- **No correlation between IIV(CL) and IIV(VD).** Moffett 2017 Table 4
  reports two separate omega estimates (omega1 for CL, omega2 for VD)
  with no cross-term; the packaged model preserves the diagonal
  structure (`etalcl` and `etalvc` as independent one-parameter
  blocks). The source paper's Discussion notes that including
  baseline AT effects on both CL and VD "resulted in unacceptable
  estimates of IIV", implying that a fully correlated OMEGA was
  explored during model development but the final published model is
  diagonal; the packaged model faithfully reproduces the published
  diagonal form.

- **Race and ethnicity not modelled.** Moffett 2017 did not report the
  cohort's racial or ethnic distribution and did not test race as a
  covariate; `population$race_ethnicity` is `"not reported by the
  paper"` and no race indicator appears in `covariateData`.

- **No maturation function.** Age, postconceptional age, and serum
  creatinine were all tested and rejected during backward elimination
  (Appendix); the final model has no maturation term on either CL or
  VD despite the paediatric range of the cohort (birth to <19 y). The
  packaged model reflects the paper's final structure, so simulations
  in very young infants (< 30 days) rely on the fixed allometric
  exponents rather than an explicit maturation function.

- **Erratum search.** A search of the Wiley British Journal of
  Clinical Pharmacology landing page for doi:10.1111/bcp.13359 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.
