---
title: "Infliximab (Hanzel 2021)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Infliximab (Hanzel 2021)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

```{r setup}
library(nlmixr2lib)
library(rxode2)
library(dplyr)
library(tidyr)
library(ggplot2)
library(PKNCA)
```

# Model and source

* Citation: Hanzel J, Bukkems LH, Gecse KB, D'Haens GR, Mathot RAA.
  Population pharmacokinetics of subcutaneous infliximab CT-P13 in
  Crohn's disease and ulcerative colitis. Aliment Pharmacol Ther.
  2021;54(10):1309-1319.
* Article: [doi:10.1111/apt.16609](https://doi.org/10.1111/apt.16609)
* Open access full text:
  <https://onlinelibrary.wiley.com/doi/10.1111/apt.16609>

CT-P13 is an infliximab biosimilar (chimeric anti-TNF-alpha IgG1
monoclonal antibody). Hanzel et al. characterised the pharmacokinetics
of the subcutaneous formulation by jointly modelling intravenous
(5 mg/kg) induction and variable-dose subcutaneous maintenance phases
of the CT-P13 1.6 phase I study (NCT02883452) in 175 adults with active
Crohn's disease or ulcerative colitis. The final model is a
two-compartment structure with first-order subcutaneous absorption,
linear elimination, and bioavailability < 1, with body weight on CL,
Vc, Vp, and Q; serum albumin on CL; and antidrug-antibody (ATI)
positivity on CL.

# Population

Hanzel 2021 Table 2 summarises the 175-subject pooled cohort:

* **Sex:** 54% female, 46% male.
* **Age:** median 36 years (IQR 29-48, range 18-70).
* **Disease:** Crohn's disease 55% (n = 97), ulcerative colitis 45%
  (n = 78). Moderate-to-severe activity at baseline.
* **Body weight:** median 69 kg (IQR 60-80, range 43-118).
* **Albumin:** median 44 g/L (IQR 40-46, range 28-54). Note SI units
  (g/L), distinct from the g/dL convention used in some other
  infliximab popPK papers (e.g., Fasanmade 2009).
* **C-reactive protein:** median 3.2 mg/L (IQR 1.2-7.9, range 0.2-89.4).
* **Concomitant immunomodulators:** 84% pre-treatment, 46% during
  treatment.
* **Concomitant corticosteroids:** 61% pre-treatment, 41% during
  treatment.
* **Treatment-emergent antibodies:** 33% developed neutralising
  antibodies during treatment.

Per-cohort design (Hanzel 2021 Methods, "Clinical study design and
patients"):

| Part | n | Induction | Maintenance |
|------|---|-----------|-------------|
| 1 cohort 1 | 13 | 5 mg/kg IV at wk 0, 2 | 5 mg/kg IV q8w to wk 54 |
| 1 cohort 2 | 11 | 5 mg/kg IV at wk 0, 2 | 120 mg SC q2w to wk 54 |
| 1 cohort 3 | 12 | 5 mg/kg IV at wk 0, 2 | 180 mg SC q2w to wk 54 |
| 1 cohort 4 |  8 | 5 mg/kg IV at wk 0, 2 | 240 mg SC q2w to wk 54 |
| 2 IV arm   | 65 | 5 mg/kg IV wk 0, 2, 6, 14, 22 | weight-based SC q2w from wk 30 |
| 2 SC arm   | 66 | 5 mg/kg IV at wk 0, 2 | weight-based SC q2w from wk 6 |

Weight-based SC dose: 120 mg if <= 80 kg, 240 mg if > 80 kg.

The same information is available programmatically via
`readModelDb("Hanzel_2021_infliximab")$population`.

# Source trace

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

| Element | Source | Value / form |
| ------- | ------ | ------------ |
| CL      | Hanzel 2021 Table 3 (Final model) | 0.355 L/d (70 kg, ALB 44 g/L, ADA-) |
| Vc      | Table 3 | 3.10 L (70 kg) |
| Vp      | Table 3 | 1.93 L (70 kg) |
| Q       | Table 3 | 0.598 L/d (70 kg) |
| Ka      | Table 3 | 0.273 /d |
| F1      | Table 3 | 0.791 (subcutaneous) |
| WT on CL  | Table 3 | Power, `(WT/70)^0.666` |
| WT on Vc  | Table 3 | Power, `(WT/70)^0.385` |
| WT on Vp  | Table 3 | Power, `(WT/70)^1.08` |
| WT on Q   | Table 3 | Power, `(WT/70)^1.26` |
| ALB on CL | Table 3 | Power, `(ALB/44)^(-0.826)` |
| ATI on CL | Table 3 | Power-of-coefficient, `1.39^ADA_POS` (Methods: categorical "on-off" power switch) |
| IIV CL    | Table 3 | CV 27.7% (omega^2 = log(1 + 0.277^2) = 0.07393) |
| IIV F1    | Table 3 | CV 16.4% (omega^2 = 0.02654) |
| IIV Vc    | Table 3 | CV 21.4% (omega^2 = 0.04478) |
| IIV Ka    | Table 3 | CV 48.5% (omega^2 = 0.21125) |
| Corr(CL, F1) | Table 3 | -0.013 |
| Corr(CL, Vc) | Table 3 |  0.028 |
| Corr(CL, Ka) | Table 3 | -0.046 |
| Corr(F1, Vc) | Table 3 |  0.00008 |
| Corr(F1, Ka) | Table 3 |  0.003 |
| Corr(Vc, Ka) | Table 3 | -0.069 |
| Residual additive    | Table 3 | 1.66 mg/L |
| Residual proportional | Table 3 | 0.102 (10.2%) |
| ODE structure: 2-cmt + depot | Methods, "Population PK model development" | first-order absorption, linear elimination, F applied at depot |

# Covariate column naming

| Source column | Canonical column used here |
| ------------- | -------------------------- |
| `WT`  | `WT`  |
| `ALB` (g/L)   | `ALB` |
| `ATI` (binary anti-drug antibody indicator) | `ADA_POS` |

# Virtual population

Original per-subject data are not publicly available. The virtual
population below is constructed to match the paper's Monte Carlo
simulation setup (Methods, "Monte Carlo simulations"): three weight
strata (50, 70, 120 kg) of 1000 patients each, with albumin sampled
from a distribution roughly matching Table 2 (median 44 g/L,
IQR 40-46, range 28-54), and all subjects modelled as ADA-negative.

```{r virtual-pop}
set.seed(2021)
n_per_stratum <- 50L  # smaller than the paper's 1000 for vignette runtime (downsampled from 200 for vignette build budget)

make_pop <- function(weight, n, id_offset) {
  tibble(
    ID      = id_offset + seq_len(n),
    WT      = weight,
    ALB     = pmax(28, pmin(rnorm(n, 44, 4), 54)),  # match Table 2 distribution
    ADA_POS = 0L,
    weight_kg = paste0(weight, " kg")
  )
}

pop <- bind_rows(
  make_pop(50,  n_per_stratum, id_offset =     0L),
  make_pop(70,  n_per_stratum, id_offset =  1000L),
  make_pop(120, n_per_stratum, id_offset =  2000L)
)
```

# Replicating Figure 3: IV induction + SC maintenance

Hanzel 2021 Figure 3 simulates IV 5 mg/kg at days 0 and 14 followed by
subcutaneous 120 mg every 2 weeks starting day 42, for 50, 70, and
120 kg virtual patients. We reproduce the same regimen.

```{r f3-events}
maintenance_days <- seq(42, 42 + 14 * 26, by = 14)   # ~52 weeks of q2w SC dosing

dose_iv <- pop %>%
  crossing(TIME = c(0, 14)) %>%
  mutate(
    AMT  = 5 * WT,           # 5 mg/kg
    EVID = 1,
    CMT  = "central",
    DV   = NA_real_
  ) %>%
  select(ID, TIME, AMT, EVID, CMT, DV, WT, ALB, ADA_POS, weight_kg)

dose_sc <- pop %>%
  crossing(TIME = maintenance_days) %>%
  mutate(
    AMT  = 120,              # 120 mg SC fixed dose
    EVID = 1,
    CMT  = "depot",
    DV   = NA_real_
  ) %>%
  select(ID, TIME, AMT, EVID, CMT, DV, WT, ALB, ADA_POS, weight_kg)

obs_times <- sort(unique(c(
  seq(0, 14, by = 1),
  seq(14, 42, by = 2),
  # every-2-days during maintenance is smooth-profile-equivalent
  # (trough_days every 14d are a subset since 14 is divisible by 2)
  seq(42, max(maintenance_days) + 14, by = 2)
)))

obs <- pop %>%
  crossing(TIME = obs_times) %>%
  mutate(
    AMT  = NA_real_,
    EVID = 0,
    CMT  = "central",
    DV   = NA_real_
  ) %>%
  select(ID, TIME, AMT, EVID, CMT, DV, WT, ALB, ADA_POS, weight_kg)

events_f3 <- bind_rows(dose_iv, dose_sc, obs) %>%
  arrange(ID, TIME, desc(EVID))

stopifnot(!anyDuplicated(events_f3[, c("ID", "TIME", "EVID")]))
```

```{r f3-simulate}
mod <- readModelDb("Hanzel_2021_infliximab")
sim_f3 <- rxSolve(mod, events_f3, keep = "weight_kg", returnType = "data.frame")
```

```{r f3-plot, fig.width = 7, fig.height = 5}
sim_f3 %>%
  filter(time > 0) %>%
  group_by(time, weight_kg) %>%
  summarise(
    median = median(Cc, na.rm = TRUE),
    lo     = quantile(Cc, 0.05, na.rm = TRUE),
    hi     = quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  ggplot(aes(x = time, color = weight_kg, fill = weight_kg)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.15, linetype = 0) +
  geom_line(aes(y = median), linewidth = 0.8) +
  geom_hline(yintercept = c(3, 7), linetype = "dotted", colour = "grey40") +
  scale_y_log10(limits = c(0.5, 200)) +
  labs(
    x = "Time (days)",
    y = "CT-P13 concentration (mg/L)",
    color = "Body weight",
    fill = "Body weight",
    title = "Replicates Hanzel 2021 Figure 3",
    subtitle = "5 mg/kg IV at days 0 and 14, then 120 mg SC q2w from day 42",
    caption = "Median and 90% prediction interval; dotted lines at 3 and 7 mg/L."
  ) +
  theme_bw()
```

The simulated profiles reproduce the qualitative pattern of Figure 3:
all three weight strata reach a 5-7 mg/L trough plateau by ~12 weeks,
with heavier patients sitting below the lighter ones at steady state
because the per-kg subcutaneous maintenance dose is lower for heavier
patients on the fixed 120 mg regimen.

# Replicating Figure 2 (VPC) for the IV 5 mg/kg arm

Figure 2 of Hanzel 2021 shows a prediction-corrected VPC by route /
dose. We reproduce the IV 5 mg/kg arm here. Subjects receive IV doses
on days 0, 14, 42 (week 6), 98 (week 14), and 154 (week 22), as in the
part-2 IV arm of the trial.

```{r vpc-iv-events}
n_vpc <- 100L  # downsampled from 300 for vignette build budget

set.seed(202101)
pop_iv <- tibble(
  ID      = 5000L + seq_len(n_vpc),
  WT      = pmax(43, pmin(rlnorm(n_vpc, log(69), 0.22), 118)),
  ALB     = pmax(28, pmin(rnorm(n_vpc, 44, 4), 54)),
  ADA_POS = rbinom(n_vpc, 1, 0.10),
  weight_kg = "VPC"
)

iv_days <- c(0, 14, 42, 98, 154)

dose_vpc <- pop_iv %>%
  crossing(TIME = iv_days) %>%
  mutate(
    AMT  = 5 * WT,
    EVID = 1,
    CMT  = "central",
    DV   = NA_real_
  ) %>%
  select(ID, TIME, AMT, EVID, CMT, DV, WT, ALB, ADA_POS)

obs_vpc_times <- sort(unique(c(
  seq(0,   14, by = 1),
  seq(14,  42, by = 2),
  seq(42, 220, by = 4)
)))

obs_vpc <- pop_iv %>%
  crossing(TIME = obs_vpc_times) %>%
  mutate(
    AMT  = NA_real_,
    EVID = 0,
    CMT  = "central",
    DV   = NA_real_
  ) %>%
  select(ID, TIME, AMT, EVID, CMT, DV, WT, ALB, ADA_POS)

events_vpc <- bind_rows(dose_vpc, obs_vpc) %>%
  arrange(ID, TIME, desc(EVID))

sim_vpc <- rxSolve(mod, events_vpc, returnType = "data.frame")
```

```{r vpc-plot, fig.width = 7, fig.height = 5}
sim_vpc %>%
  filter(time > 0) %>%
  group_by(time) %>%
  summarise(
    median = median(Cc, na.rm = TRUE),
    lo     = quantile(Cc, 0.05, na.rm = TRUE),
    hi     = quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  ggplot(aes(x = time)) +
  geom_ribbon(aes(ymin = lo, ymax = hi), alpha = 0.20, fill = "steelblue") +
  geom_line(aes(y = median), linewidth = 0.8, colour = "steelblue") +
  scale_y_log10() +
  labs(
    x = "Time (days)",
    y = "CT-P13 concentration (mg/L)",
    title = "Simulated PI for IV 5 mg/kg arm of CT-P13 1.6 part 2",
    subtitle = "5 mg/kg IV on days 0, 14, 42, 98, 154; replicates Hanzel 2021 Fig. 2 (upper-left panel) qualitatively",
    caption = "Median and 90% prediction interval."
  ) +
  theme_bw()
```

# PKNCA validation

Run PKNCA for Cmax, Tmax, AUClast, and half-life on the IV induction
phase (single 5 mg/kg dose, day 0 to day 14) of the Figure 3 cohort.
Stratification is by the `weight_kg` treatment grouping so per-weight
results can be compared.

```{r pknca}
single_dose_window <- 14
sim_nca <- sim_f3 %>%
  filter(time <= single_dose_window) %>%
  transmute(ID = id, time, Cc, treatment = weight_kg)

dose_nca <- pop %>%
  transmute(ID, time = 0, amt = 5 * WT, treatment = weight_kg)

conc_obj <- PKNCAconc(sim_nca,  Cc  ~ time | treatment + ID,
                      concu = "mg/L", timeu = "day")
dose_obj <- PKNCAdose(dose_nca, amt ~ time | treatment + ID,
                      doseu = "mg")

intervals <- data.frame(
  start     = 0,
  end       = single_dose_window,
  cmax      = TRUE,
  tmax      = TRUE,
  auclast   = TRUE,
  half.life = TRUE
)
nca_data    <- PKNCAdata(conc_obj, dose_obj, intervals = intervals)
nca_results <- pk.nca(nca_data)
nca_summary <- summary(nca_results)
knitr::kable(
  nca_summary,
  digits  = 3,
  caption = "PKNCA summary after 5 mg/kg IV induction (day 0 to day 14), by weight stratum."
)
```

## Comparison against the published typical-value half-life

Hanzel 2021 Abstract: "The half-life in a typical 70 kg patient with
serum albumin of 44 g/L was 10.8 days." This is the typical-value
terminal half-life from the model parameters CL, Vc, Vp, and Q. We
compute it analytically from the packaged model:

```{r half-life-typ}
pars <- list(CL = 0.355, Vc = 3.10, Vp = 1.93, Q = 0.598)
k10 <- pars$CL / pars$Vc
k12 <- pars$Q  / pars$Vc
k21 <- pars$Q  / pars$Vp
disc <- (k10 + k12 + k21)^2 - 4 * k10 * k21
lambda_z <- ((k10 + k12 + k21) - sqrt(disc)) / 2

tibble(
  source        = c("Hanzel 2021 Abstract (typical 70 kg, ALB 44 g/L)", "Computed from model parameters"),
  half_life_d   = c(10.8, log(2) / lambda_z)
) %>%
  knitr::kable(
    digits  = 2,
    caption = "Typical-value terminal half-life (days)."
  )
```

The analytic value is identical to the paper's reported 10.8 days,
confirming the structural-parameter set was transcribed correctly.

## Comparison against the published trough-coverage Monte Carlo

Hanzel 2021 reports Monte Carlo results for steady-state SC trough
coverage (Results, "Monte Carlo simulations"): for the 50 / 70 / 120 kg
strata, the percentage of trough concentrations exceeding 7 mg/L during
maintenance dosing was 96.4% / 90.3% / 71.9%, respectively. We
reproduce the Monte Carlo end-of-interval troughs from the Figure 3
simulation (every-2-weeks SC 120 mg phase, after week 30 to capture
near-steady-state behaviour):

```{r trough-coverage}
trough_days <- seq(42 + 14 * 14, max(maintenance_days), by = 14)  # day 238 onward
troughs <- sim_f3 %>%
  filter(time %in% trough_days) %>%
  group_by(weight_kg) %>%
  summarise(
    n     = n(),
    pct_above_3  = mean(Cc > 3) * 100,
    pct_above_5  = mean(Cc > 5) * 100,
    pct_above_7  = mean(Cc > 7) * 100,
    .groups = "drop"
  )

published <- tribble(
  ~weight_kg, ~paper_pct_above_3, ~paper_pct_above_5, ~paper_pct_above_7,
  "50 kg",   99.9, 99.0, 96.4,
  "70 kg",   99.9, 98.3, 90.3,
  "120 kg",  99.5, 90.6, 71.9
)

troughs %>%
  left_join(published, by = "weight_kg") %>%
  knitr::kable(
    digits  = 1,
    caption = "Trough coverage during 120 mg SC q2w maintenance, simulation vs Hanzel 2021."
  )
```

Differences from the published values are driven by sample size (the
vignette uses fewer subjects per stratum than the paper's 1000) and by
the simplified albumin distribution; the qualitative ordering
(50 kg > 70 kg > 120 kg) is reproduced.

# Assumptions and deviations

- **Inter-occasion variability (IOV) on CL not implemented.** Hanzel
  2021 Table 3 reports IOV on CL of 17.5% with each dosing event
  defined as a new occasion. Library models conventionally carry only
  subject-level IIV; the IOV term is documented here and dropped for
  simulation. Subject-level variability in this model is therefore
  somewhat lower than the published combined IIV+IOV total.
- **ATI is treated as binary, time-varying.** The source Methods Table 1
  lists ATI as a time-dependent on/off covariate; we map it to the
  canonical `ADA_POS` column and document the rename in
  `covariateData[[ADA_POS]]$source_name`.
- **Albumin distribution.** Per-subject albumin trajectories are not
  publicly available; we sample ALB ~ N(44, 4) g/L clipped to the
  [28, 54] g/L range reported in Table 2. Time-variation is not
  modelled.
- **Race / ethnicity not reported.** The source paper does not provide
  race / ethnicity demographics, so the `population` metadata records
  this as "Not reported."
- **PKNCA half-life vs typical-value half-life.** Per the vignette
  template guidance, the PKNCA half-life from the 14-day IV phase
  reflects whatever sampling window is given; the analytic
  typical-value half-life of 10.8 days from CL/Vc/Vp/Q is what the
  model parameters are meant to reproduce, which they do exactly.
- **Off-diagonal IIV interpretation.** Table 3 reports the off-diagonal
  IIV elements with the column header "Correlation between ..."; we
  take these at face value as correlation coefficients in [-1, 1] and
  convert to covariances via cov_ij = corr_ij * sqrt(var_i * var_j) for
  the rxode2 block-eta specification.

# Reference

* `r rxode2::rxode(readModelDb("Hanzel_2021_infliximab"))$reference`
