---
title: "Levamisole (Kreeftmeijer-Vegter 2015)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Levamisole (Kreeftmeijer-Vegter 2015)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

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

# Population

Kreeftmeijer-Vegter 2015 reports a population PK analysis of oral
levamisole in 38 children (27 male, 11 female) with frequently relapsing
idiopathic steroid-sensitive nephrotic syndrome (SSNS). The cohort
spanned ages 2.35 to 13.10 years (median 6.28 y) and weights 11 to
68 kg (median 21 kg), recruited from six study centres in India
(47.4%), the Netherlands (15.8%), Belgium (13.2%), France (10.5%),
Poland (10.5%) and Italy (2.6%). Ethnicity was Caucasian 47.4% and
Asian 50.0% per Table 1 of the source. Patients received oral
levamisole 2.5 mg/kg every other day (mean administered dose
2.45 mg/kg, SD 0.24) delivered as 5, 10, 25 or 50 mg film-coated
tablets per a weight-band dosing schedule. The PK dataset comprised
136 plasma concentrations sampled at four routine visits (weeks 8, 12,
20 and 24); samples taken pre-dose or 24 h post-dose were discarded
because 102 of 121 such samples fell below the assay LOQ (5 ng/mL).
The trial was registered as EudraCT 2005-005745-18.

The same population metadata is available programmatically via
`readModelDb("KreeftmeijerVegter_2015_levamisole")$population`.

# Source trace

Per-parameter provenance is recorded as inline `ini()` comments in
`inst/modeldb/specificDrugs/KreeftmeijerVegter_2015_levamisole.R`. The
table below collects them in one place.

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lka`                | log(1.2) (fixed)        | Table 2: ka 1.2 /h fixed in final model |
| `lcl`                | log(44)                 | Table 2: CL/F 44 L/h per 70 kg (RSE 8.5%) |
| `lvc`                | log(236)                | Table 2: V/F 236 L per 70 kg (RSE 13.3%) |
| `allo_cl` (fixed)    | 0.75                    | Results paragraph: "Allometric scaling of CL/F (power value of 0.75) to standard body weight (70 kg)" |
| `allo_vc` (fixed)    | 1                       | Table 3 row 2: "Allometric scaling of CL/F and V/F to BW (70 kg)" (V/F scales linearly) |
| `e_age_cl`           | -0.101                  | Table 2: -10.1% per life year (RSE 25%) |
| `age_ref` (fixed)    | 6.28 y                  | Table 2 caption: population median age 6.28 y (centring) |
| IIV CL/F  (CV 31.6%) | omega^2 = log(1+0.316^2) = 0.0953 | Table 2 |
| IIV V/F   (CV 41.7%) | omega^2 = log(1+0.417^2) = 0.1603 | Table 2 |
| Cov(CL,V) (rho=0.99) | 0.99 * sqrt(0.0953 * 0.1603) = 0.1223 | Methods, IIV paragraph ("V/F correlated 100% with ... CL/F"); rho clipped from 1.0 to 0.99 for PD omega (see Errata) |
| IIV ka    (CV 92.2%) | omega^2 = log(1+0.922^2) = 0.6151 | Table 2 |
| `propSd`             | 0.207                   | Table 2: proportional residual variability 20.7% (RSE 25.4%) |
| Storage-time effect  | 0.05%/day at -20 C      | Methods, PK data analysis paragraph - applied to observations BEFORE fitting, NOT a model term; see Errata |
| 1-compartment ODE     | n/a                    | Results paragraph: "A one compartment model with first order absorption and first order elimination from the central compartment fitted the data best" |
| `Cc <- 1000*central/vc` | n/a                  | Unit conversion: dose mg, V/F L => mg/L; * 1000 for ng/mL |
| `Cc ~ prop(propSd)`  | n/a                     | Methods: "Residual variability was modelled using a proportional error model" |

# Virtual cohort

The original observed data are not publicly available; the cohort
below is a 38-subject virtual reproduction of the published cohort
covariates (Table 1 of the source). Each subject receives a single
oral dose of 2.5 mg/kg levamisole. Body weight is sampled jointly with
age using a simple paediatric WT-for-AGE relationship (linear central
tendency with a small log-normal residual) covering the observed
11-68 kg and 2.35-13.10 y ranges. The simulation grid is denser
(every 5 minutes through 12 h) than the paper's sparse sampling so
that NCA half-life and AUC0-inf can be estimated reliably from the
simulated profile.

```{r cohort}
set.seed(20260618L)

n_subjects <- 38L
age_min    <- 2.35
age_max    <- 13.10
dose_mg_per_kg <- 2.5
ref_age <- 6.28
ref_wt  <- 70
sample_hours <- seq(0, 12, by = 1 / 12)

# Sample AGE from a log-normal centred at the paper's median (6.28 y)
# so the simulated cohort's median age matches Table 1, not the
# uniform-range midpoint (7.73 y). Truncate to the published 2.35-13.10
# y range, oversampling so 38 in-range values are available.
ages_raw <- exp(rnorm(n_subjects * 4, log(ref_age), 0.40))
ages     <- ages_raw[ages_raw >= age_min & ages_raw <= age_max][seq_len(n_subjects)]

# Approximate paediatric WT-for-AGE central tendency: WT ~ 2.5*AGE + 7 kg
# spans ~13 kg at 2.4 y to ~40 kg at 13 y, similar to CDC 50th percentile.
# Add a log-normal residual (sigma ~ 0.20) so the spread covers the
# published 11-68 kg range.
wt_cent <- 2.5 * ages + 7
weights <- pmax(11, pmin(68, wt_cent * exp(rnorm(n_subjects, 0, 0.20))))

subjects <- tibble::tibble(
  id   = seq_len(n_subjects),
  AGE  = ages,
  WT   = round(weights, 1),
  dose = round(dose_mg_per_kg * weights, 1)
)

dose_rows <- subjects |>
  dplyr::transmute(id, time = 0, amt = dose, cmt = "depot",
                   evid = 1L, AGE, WT)

obs_rows <- subjects |>
  tidyr::crossing(time = sample_hours) |>
  dplyr::transmute(id, time, amt = 0, cmt = NA_character_,
                   evid = 0L, AGE, WT)

events <- dplyr::bind_rows(dose_rows, obs_rows) |>
  dplyr::arrange(id, time, dplyr::desc(evid)) |>
  dplyr::mutate(treatment = "single 2.5 mg/kg oral")

stopifnot(!anyDuplicated(unique(events[, c("id", "time", "evid")])))
```

# Simulation

```{r simulate}
mod <- rxode2::rxode2(readModelDb("KreeftmeijerVegter_2015_levamisole"))

sim <- rxode2::rxSolve(
  mod,
  events = events,
  keep   = c("WT", "AGE", "treatment")
) |>
  as.data.frame()
```

For a deterministic typical-value replication (no IIV / residual
scatter), zero out the random effects:

```{r simulate-typical}
mod_typical <- mod |> rxode2::zeroRe()
sim_typical <- rxode2::rxSolve(
  mod_typical,
  events = events,
  keep   = c("WT", "AGE", "treatment")
) |>
  as.data.frame()
```

# Replicate Figure 1: prediction-corrected VPC shape

Figure 1 of Kreeftmeijer-Vegter 2015 is a prediction-corrected visual
predictive check; the model's stochastic simulation reproduces the
overall shape (rapid absorption peak near 1.5-2 h, exponential decline
with apparent terminal half-life ~2.6 h).

```{r figure-1, fig.width = 7, fig.height = 4.5}
sim_quantiles <- sim |>
  dplyr::filter(time > 0, !is.na(Cc)) |>
  dplyr::group_by(time) |>
  dplyr::summarise(
    Q05 = stats::quantile(Cc, 0.05, na.rm = TRUE),
    Q50 = stats::quantile(Cc, 0.50, na.rm = TRUE),
    Q95 = stats::quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(sim_quantiles, aes(time, Q50)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.20, fill = "steelblue") +
  geom_line(linewidth = 0.7) +
  scale_y_log10() +
  labs(
    x = "Time after dose (h)",
    y = "Cc (ng/mL)",
    title = "Replicates Kreeftmeijer-Vegter 2015 Figure 1 (VPC shape)",
    caption = "Stochastic VPC, median + 5-95% interval, n = 38, single 2.5 mg/kg oral dose."
  )
```

# Age effect on CL/F

The model encodes a linear proportional age effect on CL/F (slope
-10.1% per additional life year above the population median of
6.28 y). The plot below isolates that effect at the reference 70 kg
body weight. Negative AGE values are unphysical; the function is
plotted over the in-cohort 2.35-13.10 y range.

```{r age-effect, fig.width = 7, fig.height = 4}
e_age_cl <- -0.101
age_grid <- seq(age_min, age_max, by = 0.1)
cl_ratio <- 1 + e_age_cl * (age_grid - ref_age)

ggplot(data.frame(age = age_grid, cl = cl_ratio), aes(age, cl)) +
  geom_line(linewidth = 0.8, colour = "steelblue") +
  geom_hline(yintercept = 1, linetype = "dashed", colour = "grey50") +
  geom_vline(xintercept = ref_age, linetype = "dashed", colour = "grey50") +
  labs(
    x = "Age (years)",
    y = "CL/F multiplier (relative to age 6.28 y)",
    title = "Linear age effect on CL/F (Kreeftmeijer-Vegter 2015)",
    caption = "CL/F * (1 - 0.101 * (AGE - 6.28)); centred at population median age 6.28 y."
  )
```

# PKNCA validation

NCA on the simulated stochastic cohort over the single-dose 0-12 h
window. The simulation grid is dense enough to estimate `lambda.z` for
each subject, so AUC0-inf, half-life and Cmax / Tmax are reported.

```{r pknca}
sim_nca <- sim |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, treatment)

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

dose_pkn <- events |>
  dplyr::filter(evid == 1L) |>
  dplyr::select(id, time, amt, treatment)

conc_obj <- PKNCA::PKNCAconc(sim_nca, Cc ~ time | treatment + id,
                             concu = "ng/mL", timeu = "h")
dose_obj <- PKNCA::PKNCAdose(dose_pkn, amt ~ time | treatment + id,
                             doseu = "mg", route = "extravascular")

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

nca_res <- PKNCA::pk.nca(PKNCA::PKNCAdata(conc_obj, dose_obj, intervals = intervals))
```

## Comparison against published NCA

Kreeftmeijer-Vegter 2015 Table 4 reports the median (IQR) secondary PK
parameters derived from the individual empirical Bayes estimates of
each subject:

```{r nca_compare}
published <- tibble::tribble(
  ~treatment,               ~cmax,   ~tmax, ~aucinf.obs, ~half.life,
  "single 2.5 mg/kg oral", 438.3,    1.65,  2847,         2.60
)

cmp <- nlmixr2lib::ncaComparisonTable(
  simulated     = nca_res,
  reference     = published,
  by            = "treatment",
  units         = c(cmax = "ng/mL", aucinf.obs = "ng*h/mL",
                    tmax = "h",      half.life  = "h"),
  tolerance_pct = 20
)

knitr::kable(
  cmp,
  caption = "Simulated vs. published NCA (Kreeftmeijer-Vegter 2015 Table 4). * differs from reference by >20%."
)
```

The published secondary PK parameters in Table 4 are derived from
individual Empirical Bayes estimates of the 38 actual subjects; the
simulated medians above come from a 38-subject virtual cohort whose
weight-age covariates approximate Table 1 but are independently
sampled. Small discrepancies (typically within +/- 20%) reflect the
virtual-vs-actual covariate-sampling difference; large discrepancies
(starred rows) would indicate a structural extraction error and
warrant investigation rather than parameter tuning.

# Assumptions and deviations

- **Storage-time correction is NOT a model term.** Methods of the
  source state that "storage duration was included in the model as a
  fixed effect (0.0005 x days in storage at -20 C) directly on the
  observed concentrations, and the model was evaluated using the
  storage time-corrected observed concentrations." This is a
  pre-analytical correction applied to the observations before
  NONMEM fitting (not an in-model covariate effect on Cc). The
  packaged model file accordingly does not encode storage time;
  the typical-value predictions correspond to storage-time-corrected
  (i.e. assay-degradation-free) concentrations. Consumers fitting
  fresh samples need no adjustment; consumers using long-stored
  samples should apply the same multiplicative correction at
  data-assembly time.
- **IIV on V/F encoded as a near-perfect-correlation omega block
  (rho = 0.99).** The source's Methods paragraph states "the
  interindividual variability of V/F correlated 100% with the
  interindividual variability of CL/F (resulting in
  overparameterization of our model), the interindividual variability
  for CL/F was estimated and used to estimate V/F using an additional
  scaling parameter." Two encodings are structurally equivalent to
  that prose: (a) a single eta on CL/F multiplied by a fixed scale on
  V/F, and (b) a full 2x2 omega block with covariance ~= sqrt(var_CL
  * var_V). This model file uses (b) -- the canonical
  `etalcl + etalvc ~ c(...)` block form `checkModelConventions()`
  recognises -- with rho clipped to 0.99 (cov = 0.1223) rather than
  exactly 1.0 (cov = 0.1235). The clip is required for rxode2 to
  Cholesky-decompose the omega during stochastic simulation; with rho
  = 1 the matrix is singular and `rxSolve` errors with
  "chol(): decomposition failed". The 1% relaxation is unobservable
  in any practical simulation summary (Cmax / AUC / Tmax / half-life
  medians and IQRs are indistinguishable between rho = 0.99 and
  rho = 1) but should be noted by consumers who plan to refit this
  model to data.
- **Apparent half-life in the simulated cohort is approximate.** The
  paper's Table 4 median t1/2 (2.60 h) comes from empirical-Bayes
  individual estimates over 38 actual children; the simulated median
  in the table above uses PKNCA's lambda-z fit over a dense
  observation grid. With Vss-only kinetics (1-compartment, no
  terminal-phase complication), the two estimates should agree
  closely.
- **Cohort weight-age coupling is synthetic.** The published trial
  recruited heterogeneous Indian and European children; the
  simulated cohort here samples WT-for-AGE from a simple linear
  central-tendency formula with a log-normal residual. The intent is
  to span the published 11-68 kg / 2.35-13.10 y ranges with a
  plausible joint distribution, not to reproduce the specific
  Table 1 individuals. Consumers who need to reproduce the
  source's Table 4 secondary parameters exactly should set
  `WT` and `AGE` to the actual per-subject values from the trial
  data.
- **`linCmt()` not used.** The model is written with explicit
  `d/dt(depot)` and `d/dt(central)` ODEs so the age and weight
  effects on CL/F are maximally visible alongside the ODE structure.
  A `linCmt()` parameterisation would be equally correct.
- **No lag time, no second compartment.** Per the source's Results
  paragraph 3, a fixed or estimated absorption lag did not improve
  the fit and was therefore excluded ("insufficient data ... within
  1 h of dosing ... mechanistically highly implausible"). Multi-
  compartmental structural models were evaluated but the 1-compartment
  form fit best (Results paragraph 1).
- **Between-occasion variability (BOV) on bioavailability not
  encoded.** The source paragraph after Table 2 reports that BOV on F
  could not be reliably estimated due to data sparseness and was
  therefore omitted; the packaged model file follows the published
  final model and does not encode BOV.
- **Storage-time stability footnote: ng/mL not mg/L.** The source's
  bioanalytical assay was calibrated 5-2000 ng/mL; the model file's
  units are reported in ng/mL to match. The internal arithmetic
  multiplies `central / vc` (mg / L) by 1000 to convert to ng/mL for
  the `Cc` observation.
