---
title: "Ethambutol (Jonsson 2011)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Ethambutol (Jonsson 2011)}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Model and source

* Citation: `r rxode2::rxode2(readModelDb("Jonsson_2011_ethambutol"))$reference`
* Description: `r rxode2::rxode2(readModelDb("Jonsson_2011_ethambutol"))$description`
* Article: <https://doi.org/10.1128/AAC.00274-11>
* DDMoRE-bundle replicate of the same fit: `modellib("Jonsson_2011_ethambutol_ddmore")` (DDMODEL00000220; see `vignettes/articles/Jonsson_2011_ethambutol_ddmore.Rmd` for the .lst / .mod-sourced encoding).

This model is the publication-Table-2 extraction of Jonsson et al.'s
two-compartment ethambutol popPK model for South African pulmonary
tuberculosis patients. Parameter values come directly from the paper's
Table 2 ("Final parameter estimates by NONMEM, together with bootstrap
estimates presented as medians", NONMEM Estimate column). The companion
`Jonsson_2011_ethambutol_ddmore` entry encodes the same fit from the
upstream NONMEM control stream + listing (DDMODEL00000220) and therefore
carries slightly more decimal-precision on the IIV / IOV variances
(the paper rounds those to integer % CV). The two entries are
intentionally redundant under the standing convention that DDMoRE-source
and paper-source extractions are tracked separately because the
extraction-route fingerprints differ.

# Population

Jonsson et al. analyzed 1,869 plasma ethambutol concentrations from 189
adults with pulmonary tuberculosis pooled across two South African
centers (DP Marais SANTA Centre, n = 60; Brewelskloof Hospital,
n = 129). Demographics (Table 1):

- Sex: 54% male, 46% female.
- Age: median 36 years, range 16-72.
- Body weight: median 47 kg, range 29-86 (cohort-mean reference for the
  50 kg allometric anchor).
- BMI: median 18 kg/m^2, range 12-36; 18% had BMI < 16 kg/m^2 (severe
  malnutrition), 47% had BMI >= 18.5 kg/m^2 (not malnourished).
- Race: 16% Black, 83% Coloured, 1% White.
- HIV positive: 13% (24 of 189).
- Estimated baseline creatinine clearance: 79 mL/min, range 23-150
  (Cockcroft-Gault, truncated at 150).
- Daily ethambutol dose: 800 mg (29%), 1000 mg (17%), 1200 mg (53%),
  1500 mg (1%); combined with isoniazid, pyrazinamide, rifampin,
  and (sometimes) streptomycin.

The DP Marais SANTA Centre arm was sampled on four occasions per
subject (three random plasma samples per occasion within 0-12 h
postdose), while the Brewelskloof Hospital arm contributed a single
rich-sampling occasion per subject (predose plus 0.5, 1, 1.5, 2, 2.5,
3, 4, 6, 8 h postdose). This sampling design is what supports the
four-occasion inter-occasion variability (IOV) term on apparent oral
clearance.

The same information is available programmatically via the model's
`population` metadata
(`readModelDb("Jonsson_2011_ethambutol")$population` after the model is
loaded).

# Source trace

Per-parameter origin (also recorded as in-file comments next to each
`ini()` entry of `inst/modeldb/specificDrugs/Jonsson_2011_ethambutol.R`):

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lcl`                | log(39.9)  | Table 2 row 'CL/F (liter/h)', NONMEM Estimate column. |
| `lvc`                | log(82.4)  | Table 2 row 'V1/F (liters)', NONMEM Estimate column. |
| `lka`                | log(0.474) | Table 2 row 'ka (h^-1)', NONMEM Estimate column. |
| `lvp`                | log(623)   | Table 2 row 'V2/F (liters)', NONMEM Estimate column. |
| `lq`                 | log(34.3)  | Table 2 row 'Q (liters/h)', NONMEM Estimate column. |
| `lmtt`               | log(0.789) | Table 2 row 'MTT (h)', NONMEM Estimate column. |
| `e_wt_cl_q = 0.75`   | fixed      | Methods 'Pharmacokinetic analysis' paragraph: 'allometric scaling by body weight was introduced on all clearance and volume terms with powers of 3/4 and 1, respectively'. |
| `e_wt_vc_vp = 1.0`   | fixed      | Methods 'Pharmacokinetic analysis' paragraph, same sentence. |
| `e_hiv_pos_f`        | -0.154     | Table 2 row 'Effect of presence of HIV on F (fractional change)', NONMEM Estimate column. Table 2 footnote a closed form: F = 1 - 0.154 x HIV. |
| `etalka`             | omega^2 = 0.1521 | Table 2 row 'IIV on ka (% CV)' = 39; paper convention CV% ~= sqrt(omega^2) x 100, so omega^2 = 0.39^2 = 0.1521. Eta shrinkage 40% (Table 2 footnote d). |
| `etalmtt`            | omega^2 = 0.8649 | Table 2 row 'IIV on MTT (% CV)' = 93; omega^2 = 0.93^2 = 0.8649. Eta shrinkage 22% (Table 2 footnote e). |
| `etalcl`             | omega^2 = 0.0400 | Table 2 row 'IIV on CL/F (% CV)' = 20; omega^2 = 0.20^2 = 0.0400. Eta shrinkage 49% (Table 2 footnote f). |
| `etaiov_cl_1..4`     | omega^2 = 0.1296 (each) | Table 2 row 'IOV on CL/F (% CV)' = 36; omega^2 = 0.36^2 = 0.1296. Single value applied to all four occasions. |
| `propSd`             | 0.318      | Table 2 row 'Proportional residual error (% CV)' Estimate = 31.8; reported on the linear scale (proportion). |
| `addSd`              | 0.107      | Table 2 row 'Additive residual error (mg/liter)' Estimate = 0.107; reported on the linear scale (mg/L). |
| `d/dt(transit1)` -> `d/dt(depot)` -> `d/dt(central)` -> `d/dt(peripheral1)` | n/a | Results 'Population pharmacokinetics' paragraph: 'a two-compartment model with one transit compartment prior to first-order absorption'; Methods 'Pharmacokinetic analysis' paragraph: 'Single-compartment and multicompartment models with first-order absorption and elimination were fitted'. |
| `f(transit1) = 1 + e_hiv_pos_f * HIV_POS` | n/a | Table 2 footnote a closed form: F = 1 - 0.154 x HIV, applied to the dosing compartment. |
| `Cc = central / vc`  | n/a   | Dose units mg, vc units L -> Cc in mg/L (equal to ug/mL). Concentration units stated in Methods 'Drug quantification' paragraph. |
| `Cc ~ add(addSd) + prop(propSd)` | n/a | Methods 'Pharmacokinetic analysis' paragraph: 'Residual variability was described with both additive and proportional error terms', back-transformed from log-transformed observation space. |
| `(WT / 50)^0.75` on `cl`, `q` | n/a | Table 2 footnote a closed forms for CL/F and Q. |
| `(WT / 50)^1` on `vc`, `vp` | n/a | Table 2 footnote a closed forms for V1/F and V2/F. |

# Virtual cohort

Original observed data are not publicly available. The figures and
NCA summary below use a virtual cohort whose covariate distribution
approximates the Table 1 demographics (weight: log-normal around the
47 kg median; sex 46% female; 13% HIV-positive; dose 1000 mg QD as a
representative single-dose-level scenario), simulated to steady state
under the four-occasion IOV design that drove the original IOV term
on apparent oral clearance.

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

n_subjects     <- 60L                 # condensed from the 189-subject paper cohort
                                      # to keep the vignette wall-clock under the
                                      # 5-minute pkgdown gate while preserving
                                      # 5-95 percentile envelope shape
dose_amt_mg    <- 1000                # representative middle dose level (17% of cohort)
dose_interval  <- 24                  # hours between QD doses
n_doses        <- 21                  # 21 daily doses to approach steady state
sample_hours   <- c(0, 0.25, 0.5, 0.75, 1, 1.5, 2, 3, 4, 6, 8, 12, 18, 24)
ss_dose_index  <- n_doses - 1
ss_clock_start <- ss_dose_index * dose_interval

# Build per-subject covariate draws approximating Table 1.
subjects <- tibble::tibble(
  id      = seq_len(n_subjects),
  WT      = round(stats::rlnorm(n_subjects,
                                meanlog = log(47),
                                sdlog   = 0.20)),  # 47 kg median, modest spread
  HIV_POS = stats::rbinom(n_subjects, size = 1L, prob = 0.13),
  OCC     = 1L                                    # primary occasion for the QD cohort
)

dose_rows <- subjects |>
  dplyr::select(id, WT, HIV_POS, OCC) |>
  tidyr::expand_grid(
    time = seq.int(0L, by = dose_interval, length.out = n_doses)
  ) |>
  dplyr::mutate(
    amt  = dose_amt_mg,
    evid = 1L,
    cmt  = 1L                       # transit1
  )

obs_rows <- subjects |>
  dplyr::select(id, WT, HIV_POS, OCC) |>
  tidyr::expand_grid(
    time = ss_clock_start + sample_hours
  ) |>
  dplyr::mutate(
    amt  = 0,
    evid = 0L,
    cmt  = NA_integer_
  )

events <- dplyr::bind_rows(dose_rows, obs_rows) |>
  dplyr::arrange(id, time, dplyr::desc(evid))

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

# Simulation

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

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

For the typical-value (Figure-like) replication used in the figures
below, zero out the random effects so the prediction is deterministic.

```{r simulate-typical}
mod_typical <- mod |> rxode2::zeroRe()

# Typical reference subjects: 50 kg HIV-negative and 50 kg HIV-positive
# (so the body-weight allometric multiplier is 1.0 and the HIV
# bioavailability effect is the only difference between the two
# trajectories).
typical_events <- dplyr::bind_rows(
  events |> dplyr::filter(id == 1L) |>
    dplyr::mutate(WT = 50, HIV_POS = 0L, id = 1L),
  events |> dplyr::filter(id == 1L) |>
    dplyr::mutate(WT = 50, HIV_POS = 1L, id = 2L)
) |>
  dplyr::arrange(id, time, dplyr::desc(evid))

sim_typical <- rxode2::rxSolve(
  mod_typical,
  events = typical_events,
  keep   = c("WT", "HIV_POS", "OCC")
) |>
  as.data.frame()
```

# Replicate published figures

```{r figure-3-stoch, fig.width = 7, fig.height = 4.5}
# Approximate replication of the steady-state shape in Figure 3 (paper's
# VPC, n = 1000): the per-time 5-50-95 percentile envelope after a typical
# dose. Direct overlay against Figure 3 itself is out of scope -- the paper
# shows VPCs stratified by study location with location-specific dosing
# schedules; the figure here uses the QD vignette cohort and is intended
# as a shape-check rather than a panel-for-panel replication.
stoch_quantiles <- sim |>
  dplyr::filter(time >= ss_clock_start, time <= ss_clock_start + 24) |>
  dplyr::mutate(tad = time - ss_clock_start) |>
  dplyr::group_by(tad) |>
  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(stoch_quantiles, aes(tad, 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 (mg/L)",
    title = paste0("Steady-state Cc envelope (5/50/95th percentile), ",
                   n_subjects, " virtual subjects, 1000 mg QD"),
    caption = "Shape-check approximating the steady-state structure shown in Figure 3 of Jonsson 2011."
  )
```

```{r hiv-effect, fig.width = 7, fig.height = 4}
# HIV effect on bioavailability (paper Discussion paragraph 4, Table 2
# footnote a closed form): F = 1 - 0.154 * HIV_POS gives a 15.4%
# reduction in steady-state exposure for HIV-positive subjects.
hiv_compare <- sim_typical |>
  dplyr::filter(time >= ss_clock_start, time <= ss_clock_start + 24) |>
  dplyr::mutate(
    tad = time - ss_clock_start,
    hiv = ifelse(HIV_POS == 1L, "HIV-positive", "HIV-negative (reference)")
  ) |>
  dplyr::distinct(hiv, tad, Cc)

ggplot(hiv_compare, aes(tad, Cc, colour = hiv)) +
  geom_line(linewidth = 0.8) +
  scale_y_log10() +
  labs(
    x = "Time after dose (h)",
    y = "Cc (mg/L)",
    colour = NULL,
    title = "Typical-value steady-state profile by HIV status (1000 mg QD, WT = 50 kg)"
  ) +
  theme(legend.position = "bottom")
```

# PKNCA validation

Steady-state NCA on the simulated cohort: Cmax, Tmax, and AUC over
the 24-hour dose interval (PKNCA `auclast` between `start = 0` and
`end = 24` post-dose).

```{r pknca}
# Keep pre-dose Cc = 0 / Cmin: filter only on missing Cc. (The `Cc > 0`
# filter that used to live here was log-scale-plotting leftover; for
# PKNCA it drops the time = 0 anchor and triggers the AUC-range warning.)
pkn_in <- sim |>
  dplyr::filter(time >= ss_clock_start, time <= ss_clock_start + 24) |>
  dplyr::mutate(
    tad       = time - ss_clock_start,
    treatment = "1000 mg QD"
  ) |>
  dplyr::filter(!is.na(Cc))

# Defensive: ensure tad = 0 exists per (id, treatment).
pkn_in <- dplyr::bind_rows(
  pkn_in,
  pkn_in |> dplyr::distinct(id, treatment) |>
    dplyr::mutate(tad = 0, Cc = 0)
) |>
  dplyr::distinct(id, treatment, tad, .keep_all = TRUE) |>
  dplyr::arrange(id, treatment, tad)

dose_pkn <- events |>
  dplyr::filter(evid == 1L, time == ss_clock_start) |>
  dplyr::mutate(treatment = "1000 mg QD")

conc_obj <- PKNCA::PKNCAconc(pkn_in, Cc ~ tad | treatment + id)
dose_obj <- PKNCA::PKNCAdose(dose_pkn, amt ~ time | treatment + id,
                             route = "extravascular")

intervals <- data.frame(
  start    = 0,
  end      = 24,
  cmax     = TRUE,
  tmax     = TRUE,
  auclast  = TRUE
)

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

nca_res$result |>
  dplyr::filter(PPTESTCD %in% c("cmax", "tmax", "auclast")) |>
  dplyr::group_by(treatment, PPTESTCD) |>
  dplyr::summarise(
    median = stats::median(PPORRES, na.rm = TRUE),
    p05    = stats::quantile(PPORRES, 0.05, na.rm = TRUE),
    p95    = stats::quantile(PPORRES, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  dplyr::mutate(`NCA parameter` = nlmixr2lib::ncaParamLabel(PPTESTCD),
                .keep = "unused", .before = "treatment") |>
  knitr::kable(
    caption = paste0("Simulated steady-state NCA parameters (1000 mg QD, ",
                     n_subjects, " virtual subjects, mixed HIV status).")
  )
```

## Comparison against the paper

The paper reports the typical CL/F at 50 kg as 39.9 L/h (Table 2).
At steady state with `Dose = 1000 mg`, `tau = 24 h`, and the
HIV-negative typical bioavailability `F = 1`, the expected AUCss(0-tau)
for a 50 kg typical subject is

`AUCss = Dose * F / CL = 1000 / 39.9 ~= 25.1 mg*h/L`.

For a typical 47 kg subject (cohort median),
`CL_47 = 39.9 * (47/50)^0.75 = 38.13 L/h` and
`AUCss = 1000 / 38.13 ~= 26.2 mg*h/L`. The simulated `auclast`
median in the cohort table above should land near 25-27 mg*h/L given
the lognormal weight spread. Tmax should be ~2-3 h (consistent with
the paper's Figure 1 plasma profiles peaking in that window).

The paper does not report an NCA summary table directly comparable to
this PKNCA output -- the only typical-exposure benchmark Jonsson et
al. give is the closed-form typical CL/F itself plus the Figure 3 VPC.
Differences greater than ~20% from the closed-form expectation should
trigger investigation of the model encoding rather than tuning.

# Inter-occasion variability (optional check)

The IOV term contributes additional stochastic variation in `cl_i`
across occasions. Setting `OCC = 2`, `3`, or `4` swaps which of the
four `etaiov_cl_<k>` slots multiplies into `cl`; with all four slots
sharing a common variance (0.1296), the marginal IOV-only spread is
the same across occasions, consistent with the paper's single IOV
estimate of 36% CV.

```{r iov-spread, fig.width = 6, fig.height = 4}
iov_events <- dplyr::bind_rows(
  events |> dplyr::mutate(OCC = 1L),
  events |> dplyr::mutate(OCC = 2L)
)

set.seed(20260609L + 7L)
sim_iov <- rxode2::rxSolve(
  mod, events = iov_events,
  keep = c("OCC")
) |>
  as.data.frame() |>
  dplyr::filter(time >= ss_clock_start, time <= ss_clock_start + 24) |>
  dplyr::mutate(tad = time - ss_clock_start)

iov_summary <- sim_iov |>
  dplyr::group_by(OCC, tad) |>
  dplyr::summarise(
    Q50 = stats::quantile(Cc, 0.50, na.rm = TRUE),
    Q05 = stats::quantile(Cc, 0.05, na.rm = TRUE),
    Q95 = stats::quantile(Cc, 0.95, na.rm = TRUE),
    .groups = "drop"
  )

ggplot(iov_summary, aes(tad, Q50, colour = factor(OCC))) +
  geom_line(linewidth = 0.7) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95, fill = factor(OCC)),
              alpha = 0.15, colour = NA) +
  scale_y_log10() +
  labs(
    x = "Time after dose (h)",
    y = "Cc (mg/L)",
    colour = "OCC", fill = "OCC",
    title = paste0("Steady-state Cc by occasion (1000 mg QD, ",
                   n_subjects, " per OCC, 90% spread)")
  ) +
  theme(legend.position = "bottom")
```

# Assumptions and deviations

- **No observed data on disk.** Direct overlay against the paper's
  Figure 1 (raw concentrations) or Figure 3 (VPC) is not part of
  this vignette -- the paper's observed dataset is not redistributable.
  Validation is by closed-form CL/F vs. AUCss check + PKNCA roundtrip
  + paper-stated benchmark comparison.
- **IIV / IOV variance scale.** Table 2 reports IIV / IOV as integer
  % CV without disambiguating whether the underlying value is the
  log-scale omega-squared, the linear-scale variance, or sqrt(omega^2).
  The encoding here uses the paper's stated convention
  `CV% ~= sqrt(omega^2) * 100`, i.e. `omega^2 = (CV/100)^2`, which
  matches the rounding shown in Table 2 against the DDMORE bundle's
  byte-level OMEGA values (e.g. paper IIV CL = 20%, bundle
  OMEGA(1,1) = 0.0381 -- sqrt(0.0381) = 0.195 ~= 20%). For higher
  decimal precision, use `Jonsson_2011_ethambutol_ddmore` which
  reads the values directly from the .lst.
- **Allometric scaling is fixed, not estimated.** Methods
  'Pharmacokinetic analysis' paragraph hard-codes the theory-based
  exponents (0.75 on CL/F and Q/F, 1.0 on V1/F and V2/F) on a 50 kg
  reference weight. The packaged model wraps these in `fixed()`.
- **No IIV on V1, V2, or Q.** Results 'Population pharmacokinetics'
  paragraph: 'IIV terms were included on the absorption rate
  coefficient, CL/F, and the mean transit time, and an IOV term was
  included on CL/F.' The packaged model omits the corresponding
  `etalvc`, `etalvp`, `etalq` and the vignette simulations therefore
  have no between-subject spread on those parameters.
- **IOV `$OMEGA BLOCK(1) SAME` unrolled into four independent
  fixed-equal etas.** nlmixr2 has no NONMEM `SAME` shortcut, so the
  four occasion slots `etaiov_cl_1..4` each get their own declaration
  with the first estimated and the next three hard-fixed at the
  shared value (0.1296) to preserve the source's IOV parameterization.
  Matches the convention used by `Jonsson_2011_ethambutol_ddmore.R`,
  `Aregbe_2012_alvespimycin.R`, and `Xie_2019_agomelatine.R`.
- **Single-occasion cohorts pass OCC = 1.** For the Brewelskloof
  Hospital arm (and any downstream consumer simulating a
  single-occasion design), set `OCC = 1` so the first IOV eta is the
  one that applies; the other three slots are mutually exclusive
  zero-weighted.
- **Reference weight 50 kg, not the more common 70 kg adult.** Jonsson
  2011's South African TB cohort had a median WT of 47 kg; the model
  uses 50 kg as the allometric reference (paper Table 2 footnote a's
  closed-form covariate equations).
- **Combined-error parameterisation is `combined2` (default
  Pythagorean SD).** The paper used log-transformed observations
  with combined additive + proportional error terms; on the
  back-transformed linear scale this is the nlmixr2 default
  Pythagorean-SD form, so `Cc ~ add(addSd) + prop(propSd)` preserves
  it without an explicit `combined1()` / `combined2()` qualifier.
- **Cohort condensed from 189 to 60.** The full 189-subject paper
  cohort would double the vignette wall-clock under `rxode2::rxSolve`
  with IIV / IOV / EPS layered on; the condensed 60-subject cohort
  keeps the vignette under the 5-minute pkgdown gate without
  materially changing the per-time-point median / 5-95% spread.
  Consumers who need the full 189-subject reproduction can re-run
  with `n_subjects <- 189L` locally.
- **DDMoRE replicate carries higher decimal precision on IIV / IOV.**
  See `Jonsson_2011_ethambutol_ddmore` for the .lst-sourced encoding
  (`OMEGA(1,1) = 0.0381`, `OMEGA(3,3) = 0.153`, `OMEGA(6,6) = 0.862`,
  `OMEGA(7,7) = 0.127`) when the paper-rounded values used here are
  not sufficient.
