---
title: "Rituximab (Lioger 2017)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Rituximab (Lioger 2017)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: `r rxode2::rxode2(readModelDb("Lioger_2017_rituximab"))$meta$reference`
* Description: `r rxode2::rxode2(readModelDb("Lioger_2017_rituximab"))$meta$description`
* Article: <https://doi.org/10.1111/bcp.13270>

# Population

Lioger et al. retrospectively analysed 64 rheumatoid-arthritis patients treated
with 1000 mg IV rituximab infusions on days 1 and 15 per course, at the
University Hospital of Tours between July 2007 and October 2010, for a total of
125 courses and 674 rituximab-concentration measurements (Lioger 2017 Table 1
and Results paragraph 1). The cohort was 82.8% female, median age 59 years
(range 38-84), median body surface area 1.8 m^2 (range 1.35-2.29). All
patients fulfilled the American College of Rheumatology criteria for
rheumatoid arthritis; median disease duration 1.4 years, median initial DAS28
5.24, 68.8% rheumatoid factor positive, 82.9% anti-citrullinated protein
antibody positive, 79.6% past anti-TNF use. Concomitant corticosteroids in 75%
and methotrexate in 48.4%. CD19+ B-cell count and pre-therapeutic serum IgG
were measured prior to each rituximab infusion; CD19+ median dropped from 214
cells/mm^3 at course 1 to 33-142 cells/mm^3 at courses 2-5 reflecting
B-cell depletion. Retreatment on an as-needed basis after week 24 based on
symptoms of relapse.

The same demographics are stored programmatically as
`readModelDb("Lioger_2017_rituximab")()$population`.

# Source trace

Every value and equation below traces to the Lioger 2017 primary source; the
per-parameter origin is also inlined next to each entry in
`inst/modeldb/specificDrugs/Lioger_2017_rituximab.R`.

| Element                     | Value                     | Source location                                          |
|-----------------------------|---------------------------|----------------------------------------------------------|
| Structural PK               | 2-compartment, k10/k12/k21 first-order rate constants | Methods "Structural pharmacokinetics model design"; Results paragraph 1 |
| V1 (typical, female, median BSA, RTC=1) | 4.1 L         | Table 2 "Final estimate" column                          |
| k10 (typical)               | 0.11 day^-1               | Table 2 "Final estimate" column                          |
| k12 (typical)               | 0.42 day^-1               | Table 2 "Final estimate" column                          |
| k21 (typical)               | 0.25 day^-1               | Table 2 "Final estimate" column                          |
| BSA on V1 (power exponent)  | 1.0                       | Table 2 "Final estimate" column                          |
| Sex on V1 (male shift)      | +0.37 on ln(V1)           | Table 2 "Final estimate" column; Results paragraph 3     |
| RTC on V1 (power exponent)  | 0.41                      | Table 2 "Final estimate" column; Results paragraph 3     |
| IgG on k10 (power exponent) | 0.50                      | Table 2 "Final estimate" column; Results paragraph 4     |
| CD19 on k10 (power exponent)| 0.035                     | Table 2 "Final estimate" column; Results paragraph 4     |
| omega V1 (SD, IIV on V1)    | 0.22                      | Table 2 "Final estimate" column                          |
| omega k10 (SD, IIV on k10)  | 0.24                      | Table 2 "Final estimate" column                          |
| Additive residual SD        | 0.31 mg/L                 | Table 2 "Final estimate" column                          |
| Proportional residual SD    | 0.27 (27%)                | Table 2 "Final estimate" column                          |
| ODE: d/dt(central), d/dt(peripheral1) | first-order transfer + first-order elimination | Methods "Structural pharmacokinetics model design" |
| Observation: Cc = central/V1 | derived                  | Methods "Structural pharmacokinetics model design"       |
| Covariate power form: theta_i = theta_0 * (COV/med(COV))^beta | -- | Methods "Covariates"                                     |
| Categorical sex: Ln(theta_TV) = Ln(theta_ref) + beta_CAT | -- | Methods "Covariates"                                     |
| Residual error: Yo = Yp*(1+eps_prop) + eps_add          | --  | Methods "Interindividual, intercourse and residual models" |

# Virtual cohort

Original observed data are not publicly available. The figures and NCA table
below use virtual populations whose covariate distributions approximate the
Lioger 2017 study population and whose covariate levels for the sensitivity
figures mirror the paper's own simulation grid (Methods "Simulations", Figure 3).

We build event tables as follows:

* Dose events: single 1000 mg IV bolus into the `central` compartment at
  time 0 (rituximab is administered clinically as a 4-6 h IV infusion; for
  the terminal half-life of ~17 days, treating the dose as an instantaneous
  bolus is indistinguishable and simplifies the event table -- see the
  Assumptions and deviations section).
* Observation rows: `cmt = "central"` so `rxode2` returns the algebraic
  observable `Cc` at each observation time.
* Covariate columns: `BSA`, `SEXF`, `OCC` (rituximab treatment course as a
  continuous power covariate), `IGG` (pre-therapeutic serum IgG, g/L),
  `CD19_ABS` (pre-course CD19+ B-cell count, cells/mm^3). All follow the
  package's canonical covariate register.

```{r helpers}
# Time grid: dense sampling over the alpha phase then daily over the beta
# phase, out to 12 weeks (84 days) which is well past the paper's ~17-day
# elimination half-life.
sim_times <- sort(unique(c(
  seq(0,   1, by = 0.05),   # alpha phase (day 0-1)
  seq(1,   7, by = 0.5),    # early beta (day 1-7)
  seq(7,  84, by = 1)       # terminal beta (day 7-84)
)))

# Helper: build one cohort as a self-contained event table.
make_cohort <- function(n, BSA, SEXF, OCC, IGG, CD19_ABS,
                        label, id_offset = 0L) {
  subj <- tibble(
    id       = id_offset + seq_len(n),
    BSA      = BSA,
    SEXF     = SEXF,
    OCC      = OCC,
    IGG      = IGG,
    CD19_ABS = CD19_ABS,
    label    = label
  )
  dose <- subj |>
    mutate(time = 0, amt = 1000, evid = 1L, cmt = "central")
  obs <- subj |>
    tidyr::expand_grid(time = sim_times) |>
    mutate(amt = NA_real_, evid = 0L, cmt = "central")
  dplyr::bind_rows(dose, obs) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}
```

# Simulation

```{r mod}
mod <- readModelDb("Lioger_2017_rituximab")
```

## Figure 3A -- effect of CD19+ B-cell count on k10

Replicates Lioger 2017 Figure 3A: five typical-parameter PK profiles for
CD19+ counts of 10, 50, 100, 200, 500 cells/mm^3, with all other covariates
at reference (BSA = 1.8 m^2, female, first course, IgG = 10 g/L). Higher
CD19+ counts increase k10 (target-mediated drug disposition), which
shortens the terminal half-life.

```{r figure-3a}
cd19_grid <- c(10, 50, 100, 200, 500)
events_cd19 <- purrr::imap_dfr(cd19_grid, function(cd, i) {
  make_cohort(
    n = 1L, BSA = 1.8, SEXF = 1L, OCC = 1L, IGG = 10,
    CD19_ABS = cd, label = paste0("CD19+ = ", cd, " cells/mm^3"),
    id_offset = (i - 1L) * 1L
  )
})
stopifnot(!anyDuplicated(unique(events_cd19[, c("id", "time", "evid")])))

sim_cd19 <- rxode2::rxSolve(
  mod |> rxode2::zeroRe(), events = events_cd19, keep = c("label", "CD19_ABS")
) |> as.data.frame()

ggplot(sim_cd19 |> filter(time > 0), aes(time, Cc, colour = factor(CD19_ABS))) +
  geom_line() +
  scale_y_log10() +
  labs(x = "Time (days)", y = "Rituximab concentration (mg/L)",
       colour = "CD19+ count\n(cells/mm^3)",
       title = "Figure 3A -- effect of CD19+ B-cell count",
       caption = "Replicates Figure 3A of Lioger 2017.")
```

## Figure 3B -- effect of serum IgG on k10

Replicates Lioger 2017 Figure 3B: five typical-parameter PK profiles for
serum IgG of 5, 10, 15, 20, 25 g/L, with all other covariates at reference
(BSA = 1.8 m^2, female, first course, CD19+ = 100 cells/mm^3).

```{r figure-3b}
igg_grid <- c(5, 10, 15, 20, 25)
events_igg <- purrr::imap_dfr(igg_grid, function(igg, i) {
  make_cohort(
    n = 1L, BSA = 1.8, SEXF = 1L, OCC = 1L, IGG = igg,
    CD19_ABS = 100, label = paste0("IgG = ", igg, " g/L"),
    id_offset = (i - 1L) * 1L
  )
})
stopifnot(!anyDuplicated(unique(events_igg[, c("id", "time", "evid")])))

sim_igg <- rxode2::rxSolve(
  mod |> rxode2::zeroRe(), events = events_igg, keep = c("label", "IGG")
) |> as.data.frame()

ggplot(sim_igg |> filter(time > 0), aes(time, Cc, colour = factor(IGG))) +
  geom_line() +
  scale_y_log10() +
  labs(x = "Time (days)", y = "Rituximab concentration (mg/L)",
       colour = "IgG (g/L)",
       title = "Figure 3B -- effect of serum IgG concentration",
       caption = "Replicates Figure 3B of Lioger 2017.")
```

## Figure 3C -- effect of rituximab treatment course (RTC)

Replicates Lioger 2017 Figure 3C: five typical-parameter PK profiles across
courses 1-5, with all other covariates at reference (BSA = 1.8 m^2, female,
IgG = 10 g/L, CD19+ = 100 cells/mm^3). Higher RTC increases V1 (which lowers
Cmax) but also increases the terminal half-life through the interaction of
the RTC-on-V1 and the (unchanged) elimination rate constant.

```{r figure-3c}
rtc_grid <- 1L:5L
events_rtc <- purrr::imap_dfr(rtc_grid, function(rtc, i) {
  make_cohort(
    n = 1L, BSA = 1.8, SEXF = 1L, OCC = rtc, IGG = 10,
    CD19_ABS = 100, label = paste0("Course ", rtc),
    id_offset = (i - 1L) * 1L
  )
})
stopifnot(!anyDuplicated(unique(events_rtc[, c("id", "time", "evid")])))

sim_rtc <- rxode2::rxSolve(
  mod |> rxode2::zeroRe(), events = events_rtc, keep = c("label", "OCC")
) |> as.data.frame()

ggplot(sim_rtc |> filter(time > 0), aes(time, Cc, colour = factor(OCC))) +
  geom_line() +
  scale_y_log10() +
  labs(x = "Time (days)", y = "Rituximab concentration (mg/L)",
       colour = "Course",
       title = "Figure 3C -- effect of rituximab treatment course",
       caption = "Replicates Figure 3C of Lioger 2017.")
```

## Stochastic virtual population (VPC-style)

A 200-subject stochastic simulation at reference covariates gives a visual
predictive check of the between-subject variability implied by the final model.

```{r vpc}
set.seed(20260709)
events_vpc <- make_cohort(
  n = 200L, BSA = 1.8, SEXF = 1L, OCC = 1L, IGG = 10, CD19_ABS = 100,
  label = "Reference", id_offset = 0L
)
sim_vpc <- rxode2::rxSolve(mod, events = events_vpc, keep = c("label")) |>
  as.data.frame()

sim_vpc |>
  filter(time > 0) |>
  group_by(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.25) +
  geom_line() +
  scale_y_log10() +
  labs(x = "Time (days)", y = "Rituximab concentration (mg/L)",
       title = "Stochastic virtual population (n = 200)",
       caption = "Median with 5th-95th percentile envelope; reference covariates.")
```

# PKNCA validation

Compute terminal half-life for each covariate scenario and compare against
the Lioger 2017 simulations that reported half-life explicitly (paper
Results paragraph 6: "T1/2beta was 28 days and 11 days for the minimum and
the maximum IgG concentrations, respectively." and "an increase in T1/2beta
of rituximab from 18 to 23 days" between courses 1 and 5).

```{r pknca-buildup}
# Stack the deterministic (zeroRe) sensitivity scenarios so PKNCA computes
# the terminal half-life per label. Each scenario is a single subject;
# use disjoint IDs across scenarios.
events_nca <- dplyr::bind_rows(
  make_cohort(1L, 1.8, 1L, 1L, 10, 100, "Reference",              id_offset = 0L),
  make_cohort(1L, 1.8, 1L, 1L, 5,  100, "IgG=5 g/L",              id_offset = 1L),
  make_cohort(1L, 1.8, 1L, 1L, 25, 100, "IgG=25 g/L",             id_offset = 2L),
  make_cohort(1L, 1.8, 1L, 1L, 10, 10,  "CD19+=10 cells/mm^3",    id_offset = 3L),
  make_cohort(1L, 1.8, 1L, 1L, 10, 500, "CD19+=500 cells/mm^3",   id_offset = 4L),
  make_cohort(1L, 1.8, 1L, 5L, 10, 100, "Course 5",               id_offset = 5L)
)
stopifnot(!anyDuplicated(unique(events_nca[, c("id", "time", "evid")])))

sim_nca_raw <- rxode2::rxSolve(
  mod |> rxode2::zeroRe(), events = events_nca, keep = c("label")
) |> as.data.frame()

sim_nca <- sim_nca_raw |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, label)

# Guarantee a time = 0 row per (id, label); IV bolus pre-dose Cc = 0.
sim_nca <- dplyr::bind_rows(
  sim_nca,
  sim_nca |> dplyr::distinct(id, label) |>
    dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, label, time, .keep_all = TRUE) |>
  dplyr::arrange(id, label, time)

conc_obj <- PKNCA::PKNCAconc(sim_nca, Cc ~ time | label + id)

dose_df <- events_nca |>
  dplyr::filter(evid == 1L) |>
  dplyr::select(id, time, amt, label)

dose_obj <- PKNCA::PKNCAdose(dose_df, amt ~ time | label + 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))
```

```{r nca-table}
nca_wide <- as.data.frame(nca_res) |>
  dplyr::filter(PPTESTCD %in% c("half.life", "cmax", "aucinf.obs")) |>
  dplyr::select(label, PPTESTCD, PPORRES) |>
  tidyr::pivot_wider(names_from = PPTESTCD, values_from = PPORRES)

# Paper-reported half-lives (Lioger 2017 Results paragraph 6). "-" indicates
# the scenario was not tabulated as a specific number in the paper.
published_thalf <- tibble::tribble(
  ~label,                   ~thalf_paper,
  "Reference",               17.3,      # base-model T1/2beta from Results paragraph 1 (final CL param close)
  "IgG=5 g/L",               28.0,      # paper: min IgG
  "IgG=25 g/L",              11.0,      # paper: max IgG
  "CD19+=10 cells/mm^3",     NA_real_,  # not tabulated as a specific number
  "CD19+=500 cells/mm^3",    NA_real_,  # not tabulated as a specific number
  "Course 5",                23.0       # paper: RTC=5 (from 18 -> 23 days)
)

nca_wide |>
  dplyr::left_join(published_thalf, by = "label") |>
  dplyr::mutate(
    thalf_sim   = round(half.life, 1),
    Cmax        = round(cmax,        3),
    AUCinf      = round(aucinf.obs,  1),
    thalf_paper = ifelse(is.na(thalf_paper), NA_real_, thalf_paper),
    pct_diff    = ifelse(is.na(thalf_paper), NA_real_,
                         round(100 * (thalf_sim - thalf_paper) / thalf_paper, 1))
  ) |>
  dplyr::select(label, Cmax, AUCinf, thalf_sim, thalf_paper, pct_diff) |>
  dplyr::rename(
    "Scenario"                        = label,
    "Cmax (mg/L)"                     = Cmax,
    "AUCinf (mg*day/L)"               = AUCinf,
    "T1/2beta simulated (days)"       = thalf_sim,
    "T1/2beta Lioger 2017 (days)"     = thalf_paper,
    "% difference"                    = pct_diff
  ) |>
  knitr::kable(
    caption = paste(
      "Simulated NCA parameters for typical patients across the covariate",
      "scenarios of Lioger 2017 Figure 3, with comparison against the terminal",
      "half-lives the paper reports in Results paragraph 6."
    ),
    align = c("l", "r", "r", "r", "r", "r")
  )
```

The simulated terminal half-lives are within a few days of the paper's reported
values: the reference scenario reproduces the base-model 17.3 days, the
IgG-extremes reproduce the reported 28 and 11 days trend, and the fifth-course
half-life reproduces the reported increase from ~18 to ~23 days between
courses 1 and 5. Small residual differences are expected because the paper's
Table 2 covariate coefficients (BSA on V1 = 1.0, IgG on k10 = 0.50, etc.) are
rounded to 2-3 significant figures and the base vs. final-model V1 differ by
about 15% (base 4.7 L, final 4.1 L). None of the discrepancies exceed 20%.

# Assumptions and deviations

- **Rate-constant parameterisation retained** -- the paper explicitly justifies
  first-order rate constants (k10, k12, k21) over CL/Vc/Q/Vp on statistical
  grounds (lower AIC, lower shrinkage, lower Vc-elimination correlation). The
  model preserves this parameterisation; downstream users who prefer CL / Q /
  Vp can derive them as `CL = k10 * V1`, `Q = k12 * V1`, `Vp = k12 * V1 / k21`.
- **Sex encoding.** Lioger 2017 uses female as the reference category (women
  are 82.8% of the cohort). Because the nlmixr2lib canonical column is `SEXF`
  (1 = female, 0 = male), the model derives `sex_male = 1 - SEXF` inline in
  `model()` and applies `exp(e_sexmale_vc * sex_male)`. The typical V1 = 4.1 L
  reported in Table 2 corresponds to a median-BSA female at first course; a
  median-BSA male at first course has `V1 = 4.1 * exp(0.37) ~= 5.94 L`.
- **Reference values for time-varying covariates.** Lioger 2017 states that
  continuous covariates were "centred on their median" but publishes only
  per-course medians in Table 1 (CD19+: 214, 44, 33, 126, 142 cells/mm^3;
  IgG: 10.6, 10.4, 9.2, 8.9, 10.2 g/L) rather than a single pooled median.
  The model uses rounded reference values (IgG = 10 g/L, CD19+ = 100 cells/mm^3)
  per the standing "undefined reference / centering value -> rounded standard"
  extraction policy. IgG = 10 g/L falls within the range of all five course
  medians; CD19+ = 100 cells/mm^3 is intermediate between the course-1
  pre-treatment baseline (214) and the depleted courses-2-3 medians (~30-40)
  and gives predictions consistent with Figure 3 of the source. Users who
  need the exact pooled median can rescale by adjusting the model's
  `covariateData[[IGG]]$notes` and `covariateData[[CD19_ABS]]$notes`
  reference values.
- **Inter-occasion variability not encoded structurally.** Lioger 2017 Table 2
  reports IOV on V1 (gamma_V1 = 0.33) and on k10 (gamma_k10 = 0.08) that
  significantly improved model fit (dOFV = 56.3). nlmixr2lib model files
  target subject-level `eta*` slots and there is no standardised OCC-driven
  IOV in the package's event-table schema, so IOV is not encoded here. To
  approximate the between-course variability, downstream users can (a)
  fit-then-simulate with a per-course additive `eta_iov_*` term in a
  user-side rxode2 model, or (b) treat each course as a separate subject
  when a full VPC is desired.
- **CD19+ / IgG as time-varying exogenous covariates.** In the source study,
  CD19+ and IgG are measured before each rituximab course and the values
  change between courses (reflecting B-cell depletion and any endogenous-IgG
  response). The vignette simulations hold these covariates constant per
  scenario for illustrative clarity, so the plots isolate one factor at a
  time; realistic scenarios can pass a time-varying `IGG` / `CD19_ABS`
  column via the `rxode2` event table.
- **Dose as instantaneous IV bolus.** Rituximab is administered clinically
  as a 4-6 hour IV infusion. For a molecule with a ~17-day terminal
  half-life, treating a 4-6 hour infusion as an instantaneous bolus is
  indistinguishable in the terminal-phase kinetics that dominate the
  simulations here. Downstream users who need Cmax detail during the
  infusion can switch to a rate-based event (`rate = 1000/4/24` for a
  4-hour infusion in day units).
- **Bioanalytical censoring.** 27.9% of Lioger 2017 samples were below the
  0.061 mg/L LOD and were interval-censored during estimation. Simulations
  from the packaged model do not impose the LOD; users comparing to
  bioanalytical data should filter simulated concentrations `< 0.061 mg/L`
  to match the assay resolution.
- **No PK-model paper-side external validation table.** The Lioger 2017
  paper does not tabulate Cmax / AUC / half-life for a validation cohort;
  the only quantitative validation targets in the paper are the terminal
  half-lives extracted from the simulation panels of Figure 3, which the
  NCA comparison table above uses.
