Skip to contents

Rituximab (Lioger 2017)

Source: vignettes/articles/Lioger_2017_rituximab.Rmd
Lioger_2017_rituximab.Rmd

Model and source

  • Citation: Lioger B, Edupuganti SR, Mulleman D, Passot C, Desvignes C, Bejan-Angoulvant T, Thibault G, Gouilleux-Gruart V, Melet J, Paintaud G, Ternant D. Antigenic burden and serum IgG concentrations influence rituximab pharmacokinetics in rheumatoid arthritis patients. Br J Clin Pharmacol. 2017 Sep;83(9):1773-1781. doi:10.1111/bcp.13270. PMID:28230269.
  • Description: Two-compartment population PK model of rituximab in rheumatoid arthritis patients, parameterised with first-order distribution (k12, k21) and elimination (k10) rate constants (per Lioger 2017 Methods ‘Structural pharmacokinetics model design’: rate-constant parameterisation gave lower AIC, shrinkages, and Vc/elimination correlation than CL/Vc). Final covariate model: BSA, sex, and rituximab treatment course on V1; time-varying CD19+ B-cell count and IgG serum concentration on k10 (Lioger 2017 Table 2).
  • 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.
# 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 

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.

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()
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalvc', 'etalk10'
#> Warning: multi-subject simulation without without 'omega'

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).

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()
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalvc', 'etalk10'
#> Warning: multi-subject simulation without without 'omega'

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.

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()
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalvc', 'etalk10'
#> Warning: multi-subject simulation without without 'omega'

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.

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()
#> ℹ parameter labels from comments will be replaced by 'label()'

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).

# 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()
#> ℹ parameter labels from comments will be replaced by 'label()'
#> ℹ omega/sigma items treated as zero: 'etalvc', 'etalk10'
#> Warning: multi-subject simulation without without 'omega'

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))
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")
  )
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.
Scenario Cmax (mg/L) AUCinf (mg*day/L) T1/2beta simulated (days) T1/2beta Lioger 2017 (days) % difference
CD19+=10 cells/mm^3 243.902 2403.8 20.1 NA NA
CD19+=500 cells/mm^3 243.902 2096.3 17.7 NA NA
Course 5 126.078 1146.4 18.7 23.0 -18.7
IgG=25 g/L 243.902 1403.1 12.5 11.0 13.6
IgG=5 g/L 243.902 3135.7 25.6 28.0 -8.6
Reference 243.902 2217.8 18.7 17.3 8.1

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.