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

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

# Model and source

* Citation: `r rxode2::rxode2(readModelDb("Gandhi_2021_abatacept"))$reference`
* Description: `r rxode2::rxode2(readModelDb("Gandhi_2021_abatacept"))$description`
* Article: [J Clin Pharmacol. 2021;61(5):651-661](https://doi.org/10.1002/jcph.1781)
  (open access via [PMC8359474](https://pmc.ncbi.nlm.nih.gov/articles/PMC8359474/))

# Population

Gandhi 2021 pooled 13 phase 2/3 abatacept clinical studies into a final
population PK dataset of **2 616 patients (12 759 abatacept serum
concentrations)**: 2 213 adults with rheumatoid arthritis (RA) contributing
9 420 samples, and 403 patients aged 2-17 years with polyarticular juvenile
idiopathic arthritis (pJIA) contributing 3 339 samples. Concentrations
were quantified by validated ELISA (LLOQ 1.0 ng/mL); 8.7% of samples below
LLOQ were excluded.

Gandhi 2021 used this pooled population PK model to select a weight-tiered
weekly subcutaneous (SC) abatacept regimen for pJIA: 50 mg for patients
< 25 kg, 87.5 mg for 25 to < 50 kg, and 125 mg for >= 50 kg, designed to
deliver a steady-state trough concentration (C~minss~) >= 10 µg/mL in
> 90 % of patients across the 2-17 year age range. The selected regimen
was subsequently confirmed in a phase 3 pJIA trial (130 of 131 evaluable
patients met the C~minss~ target).

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

# Source trace

Every structural parameter, covariate effect, IIV element, and
residual-error term below is taken from Gandhi 2021 Table 2 (final
population PK model fitted to the expanded RA + pJIA dataset). Reference
covariate values come from the Figure 1 caption: 49-year-old reference
patient, baseline body weight 68 kg, albumin 4.1 g/dL, calculated GFR
99.18 mL/min/1.73 m², swollen joint count 15, not on NSAIDs. Male is the
reference category for sex per the Methods text "for sex, male was used as
the reference"; the reference patient's female designation in Figure 1 is
the visualization baseline used for the covariate-effect forest plot.

| Equation / parameter                                    | Value                                | Source location |
|---------------------------------------------------------|--------------------------------------|-----------------|
| `lcl` (CL)                                              | `log(0.0179 * 24)` L/day             | Table 2: CL = 0.0179 L/h |
| `lvc` (VC)                                              | `log(3.29)` L                        | Table 2: VC = 3.29 L |
| `lq`  (Q)                                               | `log(0.0231 * 24)` L/day             | Table 2: Q = 0.0231 L/h |
| `lvp` (VP)                                              | `log(3.67)` L                        | Table 2: VP = 3.67 L |
| `lka` (KA~TV~, relative absorption)                     | `log(0.00521 * 24)` 1/day            | Table 2: KA = 0.00521 1/h (paper "L/h" is a units typo) |
| `logitfdepot` (logit F~TV,ref~)                         | `1.20831`                            | Table 2: F = 0.770; logit(0.770) = 1.20831 |
| `e_wt_cl`  ((WT/68)^exp^ on CL)                         | `0.706`                              | Table 2: Power of body weight on CL |
| `e_wt_vc`  ((WT/68)^exp^ on VC)                         | `0.603`                              | Table 2: Power of body weight on VC |
| `e_wt_vp`  ((WT/68)^exp^ on VP)                         | `0.575`                              | Table 2: Power of body weight on VP |
| `e_wt_f`   (slope of log(WT/68) on logit-F)             | `-0.506`                             | Table 2: Power of body weight on bioavailability |
| `e_age_vc` ((AGE/49)^exp^ on VC)                        | `0.114`                              | Table 2: Power of age on VC |
| `e_age_f`  (slope of log(AGE/49) on logit-F)            | `0.487`                              | Table 2: Power of age on bioavailability |
| `e_alb_cl` ((ALB/4.1)^exp^ on CL)                       | `-0.722`                             | Table 2: Power of albumin on CL |
| `e_crcl_cl` ((CRCL/99.18)^exp^ on CL)                   | `0.259`                              | Table 2: Power of GFR on CL |
| `e_swol_cl` (((SWOL+1)/16)^exp^ on CL)                  | `0.0742`                             | Table 2: Power of SJC on CL |
| `e_sexf_cl` (exp(SEXF·coef) on CL)                      | `0.0674`                             | Table 2: Exponent of male sex on CL |
| `e_nsaid_cl` (exp(CONMED_NSAID·coef) on CL)             | `0.102`                              | Table 2: Exponent of NSAID on CL |
| `e_jia_f` (additive on logit-F for DIS_PJIA=1)          | `3.08`                               | Table 2: Exponent of JIA on bioavailability |
| `var(etalka)`                                           | `1.11`                               | Table 2 IIV/residual column, KA row |
| `var(etalvc)`                                           | `0.0464`                             | Table 2 IIV/residual column, VC row |
| `var(etalcl)`                                           | `0.0637`                             | Table 2 IIV/residual column, CL row |
| `var(etalvp)`                                           | `0.154`                              | Table 2 IIV/residual column, VP row |
| `var(etalogitfdepot)`                                   | `0.516`                              | Table 2 IIV/residual column, F row |
| `propSd` (= sqrt(SIGMA~PROP~))                          | `sqrt(0.0615)` ≈ 0.248               | Table 2 IIV/residual column: Proportional residual error = 0.0615 (variance) |
| `addSd`  (= sqrt(SIGMA~ADD~))                           | `sqrt(0.00134)` ≈ 0.0366 mg/L        | Table 2 IIV/residual column: Additive residual error = 0.00134 (mg²/L² variance) |
| Structure (2-cmt + first-order SC / instantaneous IV input + logit-F + KA > k~el~ constraint + combined residual error) | n/a | Methods, Population PK Analysis section, paragraphs on equations S1 and S2 |

## Parameterization notes

- **Time-unit conversion.** Gandhi 2021 reports `CL` and `Q` in L/h and
  `KA` in 1/h (the published Table 2 row header "KA (L/h)" is a units typo
  for the rate constant). The nlmixr2lib convention is time in days, so
  each of these values is multiplied by 24 inside `log(...)` in `ini()`.
- **Logit-F parameterisation.** Gandhi 2021 constrains absolute
  bioavailability to (0, 1) via an inverse-logit link (Methods Eq. S1):
  `F_abs = 1 / (1 + exp(-F_TV))` with `F_TV = logitfdepot +
  etalogitfdepot + DIS_PJIA * e_jia_f + e_wt_f * log(WT/68) + e_age_f *
  log(AGE/49)`. The reference value `logitfdepot = 1.20831` gives `F_abs
  ≈ 0.770` for the adult RA reference patient (68 kg, 49 yr). For a
  typical pJIA patient the +3.08 disease shift drives logit-F to ≈ 4.288,
  i.e. `F_abs ≈ 0.987`.
- **KA > k~el~ constraint.** Gandhi 2021 Methods Eq. S2 reparameterises
  absorption as `KA_i = KA_TV * exp(etaKA) + k_el,i` with `k_el,i =
  CL_i / VC_i`, to prevent flip-flop of parameter estimates during
  fitting. The model file implements this verbatim. With typical-value
  parameters (`KA_TV * 24 = 0.12504 1/day`, `k_el = 0.4296 / 3.29 ≈
  0.131 1/day`), the effective absorption rate constant is ≈ 0.256 1/day
  (half-life ≈ 2.7 days), giving an SC `Tmax` near 5 days that matches
  abatacept SC label values.
- **Independent IIVs.** The paper does not report a full IIV
  variance-covariance block (unlike Li 2019, which does). The IIVs are
  treated as independent here, matching the Table 2 layout where each ETA
  carries its own variance and shrinkage value.
- **Residual error reported as variance.** Gandhi 2021 Table 2 reports
  the SIGMA matrix variances in the same "Estimate" column as the IIV
  variances (see Errata below for the rationale). The proportional and
  additive variances are converted to standard deviations by `sqrt()` for
  nlmixr2's `add()` / `prop()` conventions: `propSd = sqrt(0.0615) ≈
  0.248` and `addSd = sqrt(0.00134) ≈ 0.0366 mg/L`.

# Virtual cohort

The simulations below use a virtual cohort whose covariate distributions
approximate the Gandhi 2021 study populations (Table S1 listing the 13
phase 2/3 contributing studies is not embedded in the PMC full text).
Subject-level observed data were not released with the paper.

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

# Adult RA cohort
n_ra <- 100  # downsampled from 300 for vignette build budget; VPC band shape preserved
ra <- tibble::tibble(
  id              = seq_len(n_ra),
  WT              = pmin(pmax(rnorm(n_ra, mean = 70, sd = 18),  40, 160)),
  AGE             = pmin(pmax(rnorm(n_ra, mean = 49, sd = 13),  18,  90)),
  ALB             = pmin(pmax(rnorm(n_ra, mean = 41, sd =    3.5), 25, 50)),
  CRCL            = pmin(pmax(rnorm(n_ra, mean = 99,  sd = 25), 30, 180)),
  SWOL_28JOINT    = pmin(pmax(round(rnorm(n_ra, mean = 15, sd = 6)), 0, 28)),
  SEXF            = rbinom(n_ra, 1, 0.78),       # RA cohorts ~75-80% female
  CONMED_NSAID    = rbinom(n_ra, 1, 0.55),
  DIS_PJIA        = 0L,
  cohort          = "RA"
)

# Pediatric pJIA cohort spanning 2-17 years and the 3 weight-tier groups.
n_pjia <- 100  # downsampled from 200 for vignette build budget; trough distribution shape preserved
pjia <- tibble::tibble(
  id              = seq.int(from = n_ra + 1, length.out = n_pjia),
  AGE             = runif(n_pjia, 2, 17),
  WT              = pmax(pmin(8 + (AGE - 2) * 4 + rnorm(n_pjia, 0, 6), 100), 8),
  ALB             = pmin(pmax(rnorm(n_pjia, mean = 41, sd = 4), 25, 50)),
  CRCL            = pmin(pmax(rnorm(n_pjia, mean = 110, sd = 25), 60, 200)),
  SWOL_28JOINT    = pmin(pmax(round(rnorm(n_pjia, mean = 6, sd = 4)), 0, 28)),
  SEXF            = rbinom(n_pjia, 1, 0.70),     # pJIA cohorts skew female
  CONMED_NSAID    = rbinom(n_pjia, 1, 0.40),
  DIS_PJIA        = 1L,
  cohort          = "pJIA"
)

pjia <- pjia |>
  dplyr::mutate(
    pjia_dose = dplyr::case_when(
      WT <  25 ~ 50,
      WT <  50 ~ 87.5,
      TRUE     ~ 125
    )
  )
```

Three regimens are simulated:

* **Adult RA, 125 mg SC QW** — the approved adult SC regimen (DIS_PJIA = 0).
* **pJIA weight-tiered SC QW** — 50 mg for < 25 kg, 87.5 mg for 25 to
  < 50 kg, 125 mg for >= 50 kg (Gandhi 2021 Tables; DIS_PJIA = 1).
* **Adult RA weight-tiered ~10 mg/kg IV Q4W** — 500 mg if < 60 kg,
  750 mg if 60-100 kg, 1000 mg if > 100 kg (the approved IV dosing
  recommendation).

```{r events}
tau_sc <- 7    # SC QW
tau_iv <- 28   # IV Q4W
n_sc   <- 26   # 26 weekly doses -> 182 days, deeply into SS
n_iv   <- 7    # 7 q4w doses -> 196 days

dose_days_sc <- seq(0, tau_sc * (n_sc - 1), by = tau_sc)
dose_days_iv <- seq(0, tau_iv * (n_iv - 1), by = tau_iv)

ra_iv <- ra |>
  dplyr::mutate(
    amt_iv = dplyr::case_when(
      WT <  60  ~ 500,
      WT <= 100 ~ 750,
      TRUE      ~ 1000
    )
  )

build_sc_events <- function(cohort_df, dose_days, treatment, dose_col = NULL,
                            fixed_amt = NULL) {
  ev_dose <- cohort_df |>
    tidyr::crossing(time = dose_days) |>
    dplyr::mutate(
      amt = if (!is.null(fixed_amt)) fixed_amt else .data[[dose_col]],
      cmt = "depot", evid = 1L, treatment = treatment
    )
  # downsampled from by = 1 for vignette build budget; PKNCA needs the
  # dose +1/+3 extras for the steady-state interval, but the smooth-line VPC
  # is visually identical at by = 2
  obs_days <- sort(unique(c(seq(0, max(dose_days) + tau_sc, by = 2),
                            dose_days + 1, dose_days + 3)))
  ev_obs <- cohort_df |>
    tidyr::crossing(time = obs_days) |>
    dplyr::mutate(amt = 0, cmt = NA_character_, evid = 0L, treatment = treatment)
  dplyr::bind_rows(ev_dose, ev_obs) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

build_iv_events <- function(cohort_iv, dose_days, treatment) {
  ev_dose <- cohort_iv |>
    tidyr::crossing(time = dose_days) |>
    dplyr::mutate(amt = amt_iv, cmt = "central", evid = 1L, treatment = treatment)
  # downsampled from by = 1 for vignette build budget; IV profile is smooth
  obs_days <- sort(unique(c(seq(0, max(dose_days) + tau_iv, by = 2),
                            dose_days + 1, dose_days + 7)))
  ev_obs <- cohort_iv |>
    tidyr::crossing(time = obs_days) |>
    dplyr::mutate(amt = 0, cmt = NA_character_, evid = 0L, treatment = treatment)
  dplyr::bind_rows(ev_dose, ev_obs) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

events_ra_sc   <- build_sc_events(ra,   dose_days_sc, "RA_SC_125_QW",
                                  fixed_amt = 125)
events_pjia_sc <- build_sc_events(pjia, dose_days_sc, "pJIA_SC_weight_tiered",
                                  dose_col = "pjia_dose")
events_ra_iv   <- build_iv_events(ra_iv, dose_days_iv, "RA_IV_weight_tiered_Q4W")
```

# Simulation

```{r simulate}
mod <- rxode2::rxode2(readModelDb("Gandhi_2021_abatacept"))
keep_cols <- c("WT", "AGE", "ALB", "CRCL", "SWOL_28JOINT",
               "SEXF", "CONMED_NSAID", "DIS_PJIA", "cohort", "treatment")

sim_ra_sc   <- as.data.frame(rxode2::rxSolve(mod, events = events_ra_sc,
                                             keep = keep_cols))
sim_pjia_sc <- as.data.frame(rxode2::rxSolve(mod, events = events_pjia_sc,
                                             keep = keep_cols))
sim_ra_iv   <- as.data.frame(rxode2::rxSolve(mod, events = events_ra_iv,
                                             keep = keep_cols))
sim         <- dplyr::bind_rows(sim_ra_sc, sim_pjia_sc, sim_ra_iv)
```

# Replicate published figures

## SC absorption Tmax in a typical adult RA patient

Gandhi 2021 Methods Eq. S2 reparameterises KA to ensure KA > k~el~. The
effective absorption rate constant for a typical adult RA patient is
`KA_TV + k_el = 0.00521 * 24 + 0.0179 * 24 / 3.29 ≈ 0.256 1/day`,
giving an SC absorption half-life of ≈ 2.7 days and a single-dose Tmax
near 5 days — consistent with the abatacept SC label.

```{r tmax-single-sc}
mod_typ <- mod |> rxode2::zeroRe()
typ_cov <- tibble::tibble(
  id = 1L, WT = 68, AGE = 49, ALB = 41, CRCL = 99.18,
  SWOL_28JOINT = 15, SEXF = 0L, CONMED_NSAID = 0L, DIS_PJIA = 0L
)
ev_single_sc <- typ_cov |>
  tidyr::crossing(time = c(0, seq(0.25, 28, by = 0.25))) |>
  dplyr::mutate(
    amt  = ifelse(time == 0, 125, 0),
    cmt  = ifelse(time == 0, "depot", NA_character_),
    evid = ifelse(time == 0, 1L, 0L)
  ) |>
  dplyr::arrange(id, time, dplyr::desc(evid))
sim_single_sc <- as.data.frame(rxode2::rxSolve(mod_typ, events = ev_single_sc))

tmax_day  <- sim_single_sc$time[which.max(sim_single_sc$Cc)]
cmax_mgL  <- max(sim_single_sc$Cc)

ggplot(sim_single_sc, aes(time, Cc)) +
  geom_line(linewidth = 0.8) +
  geom_vline(xintercept = tmax_day, linetype = "dashed", colour = "steelblue") +
  labs(
    x = "Time after single SC dose (days)",
    y = "Abatacept Cc (mg/L = µg/mL)",
    title = sprintf("Typical adult RA single 125 mg SC: Tmax = %.1f days, Cmax = %.2f mg/L",
                    tmax_day, cmax_mgL)
  ) +
  theme_minimal()
```

## pJIA C~minss~ by age category (replicates Figure 3)

Gandhi 2021 Figure 3 plots the predicted distribution of steady-state
trough concentration (C~minss~) following SC and IV abatacept in pJIA
patients, stratified by age category (2-<6, 6-<12, 12-17 years). The
paper states the SC C~minss~ in pJIA is approximately 3-fold higher than
the IV C~minss~, and the C~minss~ exceeds the 10 µg/mL near-maximal
efficacy target.

```{r fig3-pjia-cminss, fig.width = 7.5, fig.height = 4.5}
ss_start_sc <- tau_sc * (n_sc - 1)
ss_end_sc   <- ss_start_sc + tau_sc

trough_pjia <- sim_pjia_sc |>
  dplyr::filter(time == ss_start_sc + tau_sc) |>
  dplyr::mutate(age_grp = cut(AGE, breaks = c(2, 6, 12, 17.001),
                              right = FALSE,
                              labels = c("2 to <6", "6 to <12", "12-17"),
                              include.lowest = TRUE))

ggplot(trough_pjia, aes(age_grp, Cc)) +
  geom_boxplot(fill = "steelblue", alpha = 0.5, outlier.size = 0.6) +
  geom_hline(yintercept = 10, linetype = "dashed", colour = "red") +
  scale_y_log10() +
  labs(
    x = "Age category (years)",
    y = "Steady-state trough Cc (mg/L, log scale)",
    title = "Replicates Gandhi 2021 Figure 3 (pJIA SC weight-tiered Cminss)",
    caption = "Red dashed line = 10 µg/mL near-maximal-efficacy target."
  ) +
  theme_minimal()
```

## SC vs IV concentration-time profiles in adult RA

```{r ra-sc-iv, fig.width = 8, fig.height = 4.5}
vpc_ra <- dplyr::bind_rows(
  sim_ra_sc |> dplyr::mutate(regimen = "SC 125 mg QW"),
  sim_ra_iv |> dplyr::mutate(regimen = "IV weight-tiered Q4W")
) |>
  dplyr::filter(!is.na(Cc), time > 0, time <= 84) |>
  dplyr::group_by(regimen, time) |>
  dplyr::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(vpc_ra, aes(time, Q50, colour = regimen, fill = regimen)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.2, colour = NA) +
  geom_line(linewidth = 0.8) +
  labs(
    x = "Time (days)",
    y = "Abatacept Cc (mg/L)",
    title = "Simulated 5-50-95 percentile profiles in adult RA: SC vs IV",
    caption = "Virtual RA cohort (N = 300); first 12 weeks."
  ) +
  theme_minimal()
```

# PKNCA validation

Non-compartmental analysis of the steady-state SC and IV dosing intervals.
PKNCA computes per-subject C~max~, C~min~, C~avg~, and AUC~tau~. Two
PKNCA blocks are run (one per regimen-cohort grouping) so that each
formula carries `id/treatment` per the skill's PKNCA recipe.

```{r pknca-pjia-sc}
nca_conc_pjia <- sim_pjia_sc |>
  dplyr::filter(time >= ss_start_sc, time <= ss_end_sc, !is.na(Cc)) |>
  dplyr::mutate(time_nom = time - ss_start_sc) |>
  dplyr::select(id, time = time_nom, Cc, treatment)

nca_dose_pjia <- pjia |>
  dplyr::mutate(time = 0, amt = pjia_dose, treatment = "pJIA_SC_weight_tiered") |>
  dplyr::select(id, time, amt, treatment)

conc_obj_pjia <- PKNCA::PKNCAconc(nca_conc_pjia, Cc ~ time | treatment + id)
dose_obj_pjia <- PKNCA::PKNCAdose(nca_dose_pjia, amt ~ time | treatment + id)

intervals_sc <- data.frame(start = 0, end = tau_sc,
                           cmax = TRUE, cmin = TRUE, tmax = TRUE,
                           auclast = TRUE, cav = TRUE)
nca_pjia <- PKNCA::pk.nca(PKNCA::PKNCAdata(conc_obj_pjia, dose_obj_pjia,
                                           intervals = intervals_sc))
summary(nca_pjia)
```

```{r pknca-ra-sc}
nca_conc_ra <- sim_ra_sc |>
  dplyr::filter(time >= ss_start_sc, time <= ss_end_sc, !is.na(Cc)) |>
  dplyr::mutate(time_nom = time - ss_start_sc) |>
  dplyr::select(id, time = time_nom, Cc, treatment)

nca_dose_ra <- ra |>
  dplyr::mutate(time = 0, amt = 125, treatment = "RA_SC_125_QW") |>
  dplyr::select(id, time, amt, treatment)

conc_obj_ra <- PKNCA::PKNCAconc(nca_conc_ra, Cc ~ time | treatment + id)
dose_obj_ra <- PKNCA::PKNCAdose(nca_dose_ra, amt ~ time | treatment + id)

nca_ra <- PKNCA::pk.nca(PKNCA::PKNCAdata(conc_obj_ra, dose_obj_ra,
                                         intervals = intervals_sc))
summary(nca_ra)
```

## Comparison against the Gandhi 2021 C~minss~ >= 10 µg/mL benchmark

Gandhi 2021 reports that the selected weight-tiered subcutaneous regimen
delivered C~minss~ >= 10 µg/mL in **130 of 131** evaluable pJIA patients
across the 2-17 year age range — i.e. > 99 %. The cohort-fraction check
below is a direct numerical replicate of that claim against the virtual
pJIA cohort.

```{r benchmark}
pct_above_pjia <- trough_pjia |>
  dplyr::group_by(age_grp) |>
  dplyr::summarise(
    N             = dplyr::n(),
    pct_ge_10_mgL = 100 * mean(Cc >= 10, na.rm = TRUE),
    median_Cmin   = median(Cc),
    .groups = "drop"
  )

knitr::kable(pct_above_pjia, digits = 1,
  caption = "Fraction of virtual pJIA subjects achieving steady-state Cmin >= 10 mg/L on the weight-tiered SC regimen (Gandhi 2021: 130/131 = 99.2%).")

trough_ra <- sim_ra_sc |>
  dplyr::filter(time == ss_start_sc + tau_sc) |>
  dplyr::mutate(wt_grp = cut(WT, breaks = c(-Inf, 60, 100, Inf),
                             labels = c("<60 kg", "60-100 kg", ">100 kg")))

pct_above_ra <- trough_ra |>
  dplyr::group_by(wt_grp) |>
  dplyr::summarise(
    N             = dplyr::n(),
    pct_ge_10_mgL = 100 * mean(Cc >= 10, na.rm = TRUE),
    median_Cmin   = median(Cc),
    .groups = "drop"
  )

knitr::kable(pct_above_ra, digits = 1,
  caption = "Fraction of virtual adult RA subjects achieving steady-state Cmin >= 10 mg/L on 125 mg SC QW (Gandhi 2021: comparable to IV ~10 mg/kg Q4W).")
```

# Assumptions and deviations

- **Time-unit rescaling.** Gandhi 2021 reports CL, Q, and KA in per-hour
  units; the model file converts each to per-day by the factor of 24 to
  match the nlmixr2lib `units$time = "day"` convention. VC, VP, F, IIV
  variances, and residual-error magnitudes carry through unchanged.
- **KA "L/h" units typo.** Gandhi 2021 Table 2 lists KA's units as "L/h",
  but KA is a first-order absorption rate constant with units of 1/time.
  The model file treats KA as 1/h (then converts to 1/day) consistent
  with the rest of Methods (Eq. S2 expresses KA as a sum of rate
  constants).
- **SJC 28-joint scale assumed.** Gandhi 2021 reports a swollen joint
  count covariate (SJC) with reference value 15 but does not explicitly
  identify the joint-count scale. Per operator decision, the canonical
  `SWOL_28JOINT` is used because (1) the same author group used the
  28-joint count in the prior Li 2019 RA-only analysis, (2) the reference
  value 15 is consistent with the 28-joint scale's 0-28 range, and (3)
  pediatric pJIA patients are a minority of the pooled cohort
  (403/2616 ≈ 15 %).
- **Sex covariate direction.** Gandhi 2021 Methods state "for sex, male
  was used as the reference"; the model codes `SEXF = 1` for female and
  applies the published Table 2 coefficient `+0.0674` as
  `exp(SEXF * 0.0674)`, so females have ≈ 7 % higher CL than males. The
  sign is opposite to Li 2019 (`-0.0722` for SEX in the RA-only cohort);
  the dataset expansion to pJIA in Gandhi 2021 produced the sign flip.
  Convention follows the existing Li 2019 implementation in nlmixr2lib.
- **Independent IIVs (no full block).** Gandhi 2021 Table 2 reports each
  ETA's variance and shrinkage independently, with no
  variance-covariance block reported in the paper text. The model treats
  the IIVs as independent. Li 2019 (the same group's prior model) did
  report a 4×4 full block on `etalcl + etalvc + etalq + etalvp`, so users
  who want a closer match to the Li 2019 covariance structure can
  override the IIV block downstream.
- **pJIA-on-CL not retained.** Gandhi 2021 explicitly states that
  "Disease (pJIA vs RA) did not have a clinically relevant effect on
  abatacept CL" and the final-model Table 2 does not include a JIA
  effect on CL. The pJIA covariate enters the bioavailability submodel
  only.
- **Virtual-cohort covariate distributions.** Table S1 (the per-study
  demographic summary) is referenced in the paper but is not embedded in
  the PMC full text, so the virtual-cohort distributions are approximate.
  Adult RA: WT ~ N(70, 18), AGE ~ N(49, 13), ALB ~ N(4.1, 0.35),
  CRCL ~ N(99, 25), SWOL_28JOINT ~ rounded N(15, 6), SEXF = 78 % female,
  CONMED_NSAID = 55 %. pJIA: AGE uniform on [2, 17], WT linearly scaled
  with age plus noise (8 + 4·(AGE-2) ± 6 kg), SEXF = 70 % female. The
  reference covariate values match the Gandhi 2021 Figure 1 caption; the
  dispersions approximate a typical phase 3 RA / pJIA cohort.
- **IV as instantaneous bolus.** Gandhi 2021 administers IV abatacept as
  a 30-minute infusion. This simulation treats IV doses as an
  instantaneous bolus to `central`, which slightly overstates the early
  post-dose Cmax but has no effect on AUC, trough, or half-life
  estimates. Users wanting exact infusion kinetics should set `rate` (or
  `dur`) in the rxode2 event table.
- **No subject-level observed data.** Gandhi 2021 does not release
  individual subject data; the figures reported here are a pure forward
  simulation of the published final-model parameters against the virtual
  cohort.

# Errata

The published Gandhi 2021 paper contains one ambiguity that a reviewer
re-reading the source alongside this model should be aware of:

- **Residual error scale (Table 2).** Gandhi 2021 Table 2 reports the
  proportional and additive residual errors (0.0615 and 0.00134) in the
  same "Estimate" column that carries the IIV variances for each ETA.
  The IIV column entries are unambiguously variances (`var(etalka) =
  1.11` is too large to be a SD), so by symmetry the residual-error
  entries are likely SIGMA (variance) values too. Per operator decision,
  this model treats the residual errors as variances and converts to SDs
  via `sqrt()`: `propSd = sqrt(0.0615) ≈ 0.248` (24.8 % proportional
  error) and `addSd = sqrt(0.00134) ≈ 0.0366 mg/L` (≈ 37 ng/mL
  additive). The alternative interpretation — SDs directly, matching the
  Li 2019 nlmixr2lib convention from the same author group — would give
  `propSd = 0.0615` (6.15 %) and `addSd = 0.00134 mg/L` (1.34 ng/mL),
  which is unusually tight for a biologic but plausible given Gandhi
  2021's larger and more richly-covariated dataset. The paper does not
  explicitly label the residual-error column as variance vs SD; without
  the supplement (control stream) on disk this cannot be disambiguated
  beyond the column-symmetry argument used here. If a future reader
  obtains the supplement and the residual-error parameterisation
  conflicts with this interpretation, the `propSd` / `addSd` lines in
  `ini()` should be updated.
- **No formal erratum.** A search of PubMed and the J Clin Pharmacol
  corrections feed on 2026-04-25 returned no published corrections.
