---
title: "Abatacept (Lon 2013)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Abatacept (Lon 2013)}
  %\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: `r rxode2::rxode(readModelDb("Lon_2013_abatacept"))$reference`
* Description: `r rxode2::rxode(readModelDb("Lon_2013_abatacept"))$description`
* Article: [J Pharmacokinet Pharmacodyn 40(6):701-712](https://doi.org/10.1007/s10928-013-9341-1)

# Abatacept population PK simulation (Lon 2013)

Simulate abatacept plasma concentration-time profiles in collagen-induced
arthritic (CIA) Lewis rats using the final population PK model from
Lon et al. 2013. Abatacept (CTLA-4Ig Fc-fusion protein; ~92 kDa) blocks
the CD28 co-stimulatory signal required for T-cell activation. Lon 2013
characterized abatacept PK in a rat CIA model using IV and SC single-dose
and SC multiple-dose regimens, fitting all three arms simultaneously with
a two-compartment linear-elimination model and a short-term zero-order SC
absorption pattern (Eq. 3 of the paper).

## Population

From Lon 2013 Methods (Animals and Experimental Design): fifty male
Lewis rats (150-175 g at arrival, 6-9 weeks old) were purchased. CIA was
induced by intradermal injection of bovine type II collagen on day 0 and
boosted on day 7. On day 20, twenty-three rats with >= 50% paw-swelling
response in one or two hind paws were randomized to four treatment arms:

* Vehicle control (single n = 3, multiple n = 3)
* IV 10 mg/kg single dose on day 21 (n = 5)
* SC 20 mg/kg single dose on day 21 (n = 6)
* SC 20 mg/kg on day 21 followed by SC 10 mg/kg on days 23, 25, 27, and
  29 (n = 6)

The population PK fit used the 17 abatacept-treated rats. Plasma
abatacept was quantified by a soluble CTLA-4 ELISA (LLOQ 0.16 ng/mL,
inter-day CV ~15%).

The same information is available programmatically as
`readModelDb("Lon_2013_abatacept")` (the returned function's body holds
the `population` list literal).

## Source trace

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

| Element                         | Source location                                 | Value / form                                        |
| ------------------------------- | ----------------------------------------------- | --------------------------------------------------- |
| Structural model                | Lon 2013 Figure 1 and Results, "Pharmacokinetics" | 2-compartment linear elimination; short-term zero-order SC input |
| SC absorption (Eq. 3)           | Lon 2013 Methods, "Pharmacokinetic Model"       | `dAsc/dt = -kzero`; `kzero = F*Dose/tau` for `t < tau`, else 0 |
| CL  (per kg)                    | Lon 2013 Table 1                                | 21.8 mL/day/kg                                      |
| CLD (per kg)                    | Lon 2013 Table 1                                | 27.5 mL/day/kg                                      |
| V1  (per kg)                    | Lon 2013 Table 1                                | 69.5 mL/kg                                          |
| V2  (per kg)                    | Lon 2013 Table 1                                | 61.9 mL/kg                                          |
| F (SC bioavailability)          | Lon 2013 Table 1                                | 59.2% (0.592)                                       |
| tau (SC input duration)         | Lon 2013 Table 1                                | 2.67 day                                            |
| IIV on CL                       | Lon 2013 Table 1                                | 9.67% CV (omega^2 = log(CV^2 + 1) = 0.009307)       |
| IIV on V1                       | Lon 2013 Table 1                                | 56.2% CV (omega^2 = 0.274478)                       |
| Off-diagonal omega              | Lon 2013 Table 1                                | None reported; IIV treated as diagonal              |
| Proportional residual error     | Lon 2013 Table 1, epsilon_1                     | 16.1% -> propSd = 0.161 (fraction)                  |
| Additive residual error         | Lon 2013 Table 1, epsilon_2                     | 0.0365 umol/L (matches concentration units)         |
| Dose regimens                   | Lon 2013 Experimental Design                    | IV 10 mg/kg; SC 20 mg/kg single; SC 20 + 10x4 mg/kg |
| Body-weight scaling             | Lon 2013 (per-kg PK parameters)                 | CL/CLD (mL/day) = per-kg * WT; V1/V2 (mL) = per-kg * WT |
| Molecular weight                | Lon 2013 Discussion: "IC50 0.731 umol/L, circa 67.3 ug/mL" implies ~92 kDa; matches Orencia label (~92,300 Da) | 92,000 g/mol (fixed, used only for mg/mL -> umol/L conversion) |

## Virtual cohort

Individual body weights were not reported. We build a 17-rat virtual
cohort whose arm sizes match Lon 2013 exactly (IV n = 5, SC single
n = 6, SC multiple n = 6) and whose weights are sampled from a uniform
distribution across the published 150-175 g range at arrival (no weight
gain is modeled over the short post-dosing observation window).

```{r virtual-pop}
set.seed(2013)

cohorts <- tribble(
  ~treatment,             ~n,
  "IV 10 mg/kg",           5,
  "SC 20 mg/kg single",    6,
  "SC multiple",           6
)

pop <- cohorts |>
  group_by(treatment) |>
  summarise(ID = list(seq_len(n[1])), .groups = "drop") |>
  tidyr::unnest(ID) |>
  mutate(
    ID = dplyr::row_number(),
    WT = runif(dplyr::n(), min = 0.150, max = 0.175)
  )

pop
```

## Dosing and event table

Doses are delivered into `depot` for SC arms and into `central` for IV.
The packaged model uses `podo(depot)` and `tad(depot)` to drive a
zero-order release of `F * Dose / tau` into `central` for the first
`tau = 2.67` days after each SC dose; IV doses bypass this mechanism.

```{r events}
obs_times <- sort(unique(c(
  seq(0,   1,   by = 1/12),           # every 2 h on day 1
  seq(1,   3,   by = 1/4),            # every 6 h through day 3
  seq(3,  10,   by = 0.5),            # twice daily through day 10
  seq(10, 30,   by = 1)               # daily through day 30
)))

sc_multiple_times <- c(0, 2, 4, 6, 8)
sc_multiple_amts  <- c(20,  10, 10, 10, 10)  # mg/kg

dose_rows_iv <- pop |>
  filter(treatment == "IV 10 mg/kg") |>
  transmute(
    ID, treatment, WT,
    time = 0,
    amt  = 10 * WT,           # mg/kg * kg = mg
    evid = 1L,
    cmt  = "central",
    dv   = NA_real_
  )

dose_rows_sc_single <- pop |>
  filter(treatment == "SC 20 mg/kg single") |>
  transmute(
    ID, treatment, WT,
    time = 0,
    amt  = 20 * WT,
    evid = 1L,
    cmt  = "depot",
    dv   = NA_real_
  )

dose_rows_sc_multi <- pop |>
  filter(treatment == "SC multiple") |>
  tidyr::crossing(
    dose_idx = seq_along(sc_multiple_times)
  ) |>
  transmute(
    ID, treatment, WT,
    time = sc_multiple_times[dose_idx],
    amt  = sc_multiple_amts[dose_idx] * WT,
    evid = 1L,
    cmt  = "depot",
    dv   = NA_real_
  )

dose_rows <- bind_rows(dose_rows_iv, dose_rows_sc_single, dose_rows_sc_multi)

obs_rows <- pop |>
  select(ID, treatment, WT) |>
  tidyr::crossing(time = obs_times) |>
  mutate(
    amt  = NA_real_,
    evid = 0L,
    cmt  = NA_character_,
    dv   = NA_real_
  )

events <- bind_rows(dose_rows, obs_rows) |>
  arrange(ID, time, desc(evid))
```

## Simulation

Simulate with between-subject variability on CL and V1 so the spread
across the virtual cohort matches the paper's individual variability.

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

events_sim <- events |> rename(id = ID)
# `treatment` is already on every row of `events_sim` (each per-arm dose
# builder carries it through, and the obs rows pull it from `pop`). Carry
# it through rxSolve via `keep = ` rather than a fragile post-hoc
# left_join from pop.
sim <- rxSolve(object = mod, events = events_sim, returnType = "data.frame",
               keep = "treatment") |>
  as_tibble()
```

## Replicate Figure 2: plasma concentration-time profiles by treatment

Figure 2 of Lon 2013 shows abatacept plasma concentration (umol/L) vs.
time following IV 10 mg/kg, SC 20 mg/kg single, and SC multiple-dose
regimens, with 90% prediction intervals from the population model.

```{r figure-2, fig.width = 8, fig.height = 5}
fig2 <- sim |>
  filter(time > 0, !is.na(Cc), Cc > 0) |>
  group_by(treatment, time) |>
  summarise(
    med = median(Cc),
    q05 = quantile(Cc, 0.05),
    q95 = quantile(Cc, 0.95),
    .groups = "drop"
  )

ggplot(fig2, aes(x = time, y = med, colour = treatment, fill = treatment)) +
  geom_ribbon(aes(ymin = q05, ymax = q95), alpha = 0.2, colour = NA) +
  geom_line(linewidth = 0.8) +
  facet_wrap(~treatment, nrow = 1) +
  scale_y_log10() +
  labs(
    x       = "Time since first dose (days)",
    y       = "Abatacept plasma concentration (umol/L, log scale)",
    colour  = "Treatment",
    fill    = "Treatment",
    title   = "Figure 2 replication - simulated abatacept PK by treatment arm",
    caption = "Replicates Figure 2 of Lon et al. 2013 (median and 90% interval across virtual cohort)."
  ) +
  theme_bw() +
  theme(legend.position = "none")
```

# PKNCA validation

Compute non-compartmental CL (IV arm), AUC, Cmax, Tmax, and terminal
half-life per subject per treatment arm. The PKNCA formula groups
concentrations by `treatment + id` so summaries are per-arm. For the SC
single-dose arm we also derive bioavailability F by comparing dose-
normalized AUC(0,inf) against the IV arm.

```{r pknca-prep}
nca_conc <- sim |>
  filter(time >= 0, !is.na(Cc), Cc > 0) |>
  select(id, time, Cc, treatment)

nca_dose <- dose_rows |>
  transmute(id = ID, time, amt, treatment)
```

```{r pknca-run, message = FALSE, warning = FALSE}
conc_obj <- PKNCAconc(nca_conc, Cc ~ time | treatment + id)
dose_obj <- PKNCAdose(nca_dose, amt ~ time | treatment + id)

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

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

sim_nca <- as.data.frame(nca_res$result) |>
  filter(PPTESTCD %in% c("cmax", "tmax", "aucinf.obs", "half.life", "cl.obs")) |>
  group_by(treatment, PPTESTCD) |>
  summarise(
    mean = mean(PPORRES, na.rm = TRUE),
    sd   = sd(PPORRES,   na.rm = TRUE),
    .groups = "drop"
  )

knitr::kable(sim_nca, digits = 3, caption = "Simulated NCA summaries by treatment arm.")
```

## Comparison against Lon 2013 NCA values

Lon 2013 Results (Pharmacokinetics) reports the following NCA-derived
parameters across the abatacept-treated arms:

* CL = 20.8 mL/day/kg
* Vss = 146 mL/kg
* F = 57.7% (SC vs. IV)

The model-estimated values (Table 1) are CL = 21.8 mL/day/kg, V1 + V2 =
69.5 + 61.9 = 131.4 mL/kg, and F = 59.2%. The simulation below converts
the PKNCA output to per-kg clearance (using each rat's WT) and derives F
from dose-normalized AUC(0,inf).

```{r nca-compare}
# Per-kg CL for the IV arm
iv_cl <- as.data.frame(nca_res$result) |>
  filter(treatment == "IV 10 mg/kg", PPTESTCD == "cl.obs") |>
  inner_join(pop |> select(id = ID, WT), by = "id") |>
  mutate(cl_mL_per_day_per_kg = PPORRES / WT * 1000)  # L/day -> mL/day, then / kg

# Dose-normalized AUC for F calculation
per_subject_auc <- as.data.frame(nca_res$result) |>
  filter(PPTESTCD == "aucinf.obs") |>
  inner_join(pop |> select(id = ID, WT), by = "id") |>
  mutate(
    dose_mgkg = case_when(
      treatment == "IV 10 mg/kg"        ~ 10,
      treatment == "SC 20 mg/kg single" ~ 20,
      TRUE                              ~ NA_real_
    ),
    auc_per_dose = PPORRES / dose_mgkg
  ) |>
  filter(!is.na(dose_mgkg))

mean_iv_auc_per_dose <- per_subject_auc |>
  filter(treatment == "IV 10 mg/kg") |>
  summarise(mean_auc_per_dose = mean(auc_per_dose)) |>
  pull(mean_auc_per_dose)

mean_sc_auc_per_dose <- per_subject_auc |>
  filter(treatment == "SC 20 mg/kg single") |>
  summarise(mean_auc_per_dose = mean(auc_per_dose)) |>
  pull(mean_auc_per_dose)

F_sim <- mean_sc_auc_per_dose / mean_iv_auc_per_dose

comparison <- tribble(
  ~parameter,                            ~source,                  ~value,
  "CL (mL/day/kg), NCA IV arm",          "Lon 2013",                20.8,
  "CL (mL/day/kg), NCA IV arm",          "Simulation",              round(mean(iv_cl$cl_mL_per_day_per_kg), 2),
  "F (SC vs IV, fraction), NCA",         "Lon 2013",                0.577,
  "F (SC vs IV, fraction), NCA",         "Simulation",              round(F_sim, 3),
  "CL (mL/day/kg), model theta",         "Lon 2013 Table 1",        21.8,
  "Vss (mL/kg), model V1+V2",            "Lon 2013 Table 1",        131.4,
  "F (fraction), model theta",           "Lon 2013 Table 1",        0.592,
  "Vss (mL/kg), NCA",                    "Lon 2013",                146.0
)

knitr::kable(comparison, digits = 3,
             caption = "Simulated NCA vs. Lon 2013 NCA / model values.")
```

The simulated IV CL and SC F track the published NCA values within the
~20% tolerance described in the verification checklist (differences are
driven by two effects - the structural-model CL (21.8) and NCA CL (20.8)
already differed by ~5% in the source paper, and the finite-sample mean
across 5 simulated rats has its own variability under the 9.67% CV on
CL).

## Assumptions and deviations

* **Body weight distribution.** Lon 2013 reports only that rats weighed
  150-175 g at arrival. We sample `WT` uniformly across that range; no
  weight gain is modeled over the 8-day (SC multiple) or 20-30 day
  observation windows.
* **Sampling schedule.** Individual sampling times were not published
  in full detail. We use a dense early schedule (every 2-6 h for 3
  days, twice daily through day 10, daily thereafter) that supports
  stable NCA on both the early-distribution and late-terminal phases.
* **Assay.** Observations are simulated above the lower limit of
  quantification (0.16 ng/mL) by construction; the filtering step
  `!is.na(Cc) & Cc > 0` removes any post-elimination near-zero
  concentrations that would destabilize the log-linear terminal
  regression in PKNCA.
* **Residual error.** The combined prop + add error model yields an
  additive floor of 0.0365 umol/L (~ 3.4 ng/mL) plus a 16.1%
  proportional component. Simulated concentrations include residual
  variability; the Figure 2 replication plots the median and 90%
  interval across the virtual cohort to match the paper's prediction
  band.
* **PD and disease progression.** Lon 2013 also reports a PD model for
  paw edema (Table 2) built on top of these PK parameters. The PD arm
  is not packaged in this model; `Lon_2013_abatacept()` is PK-only.

## Notes

* **Structural model:** 2-compartment with linear elimination from the
  central compartment. SC doses go to `depot`; a zero-order rate
  `F * Dose / tau` is driven from `depot` into `central` for `t < tau`
  and is zero thereafter, leaving `(1 - F) * Dose` as an unabsorbed
  residual in `depot`.
* **Units.** Dose in mg and volumes in mL give `central / vc` in
  mg/mL; the observation equation divides by MW (92,000 g/mol) and
  multiplies by 1e6 to report concentrations in umol/L, matching the
  units in Lon 2013 Table 1.
* **IIV:** diagonal on CL and V1 only, per Lon 2013 Table 1; no
  IIV was reported on CLD, V2, F, or tau.

## Reference

* `r readModelDb("Lon_2013_abatacept")()$meta$reference`
