---
title: "Bedaquiline (Svensson 2013)"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Bedaquiline (Svensson 2013)}
  %\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(tibble)
library(ggplot2)
```

# Model and source

* Citation: `r rxode2::rxode(readModelDb("Svensson_2013_bedaquiline"))$reference`
* Description: `r rxode2::rxode(readModelDb("Svensson_2013_bedaquiline"))$description`
* Article: <https://doi.org/10.1128/AAC.00191-13>

# Population

The packaged model was fit to data from ACTG study A5267, a phase I
single-dose drug-drug-interaction (DDI) study designed to characterise the
effect of multiple-dose efavirenz on single-dose bedaquiline pharmacokinetics
in healthy adult volunteers in the United States. Thirty-seven subjects (18
to 65 years, no clinical evidence of TB, HIV-negative, normal QT intervals)
were enrolled at four ACTG sites; 35 completed the first PK sampling period
and were included in the modelling analysis, with 33 completing both
sampling periods. Baseline demographics (Svensson 2013 Table 2): median age
44 years (range 19 to 62), median weight 82.3 kg (range 57.3 to 118.8),
median BMI 25.9 kg/m^2 (range 19.0 to 36.1), 91% male, 69% white
non-Hispanic, 23% Black non-Hispanic, 6% Hispanic (regardless of race), and
3% Asian or Pacific Islander. CYP2B6 composite 516/983 metabolizer genotype
was assigned to 54% extensive, 37% intermediate, and 9% slow metabolizers
and tested as a covariate on the magnitude of the EFV induction effect; it
was not retained in the final model.

Each subject received a single 400 mg oral bedaquiline dose on Day 1 after
a standardised 670 kcal / 33% fat breakfast, with dense PK sampling at
predose and at 1, 2, 3, 4, 5, 6, 8, 12, 24, 48, 72, 120, 168, 216, 264, and
336 hours postdose (up to 14 days). From Day 15 to Day 42, subjects took
600 mg efavirenz fasting each evening. On Day 29, after 14 days of daily
EFV, a second 400 mg bedaquiline dose was administered with the same dense
14-day sampling schedule. Plasma BDQ and M2 concentrations were quantified
in the original study (Dooley et al. 2012); for this analysis, M3 was
additionally measured at a subset of timepoints (predose and 3, 6, 24, 72,
168, 264, and 336 h) yielding 1,152 BDQ + 1,152 M2 + 560 M3 observations
in total.

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

# Source trace

The per-parameter origin is recorded as an in-file comment next to each
`ini()` entry in
`inst/modeldb/specificDrugs/Svensson_2013_bedaquiline.R`. The table below
collects them in one place for review.

| Equation / parameter | Value | Source location |
|----------------------|-------|-----------------|
| `lmtt` (MTT)         | log(1.31) | Table 3 "MTT = 1.31 h" (RSE 12.6%) |
| `lka` (ka)           | log(0.128) | Table 3 "KA = 0.128 1/h" (RSE 8.7%) |
| `lnn` (NN)           | log(5.21) | Table 3 "NN = 5.21" (RSE 20.5%) |
| `lcl` (CL/F BDQ)     | log(2.96) | Table 3 "CL = 2.96 L/h" (RSE 9.5%) |
| `lvc` (V/F BDQ)      | log(17.3) | Table 3 "V = 17.3 L" (RSE 18.7%) |
| `lq` (Q1/F BDQ)      | log(5.01) | Table 3 "Q1 = 5.01 L/h" (RSE 8.3%) |
| `lvp` (VP1/F BDQ)    | log(2870) | Table 3 "VP1 = 2,870 L" (RSE 15.3%) |
| `lq2` (Q2/F BDQ)     | log(4.16) | Table 3 "Q2 = 4.16 L/h" (RSE 10.2%) |
| `lvp2` (VP2/F BDQ)   | log(136)  | Table 3 "VP2 = 136 L" (RSE 9.0%) |
| `lcl_m2` (CL_M2)     | log(12.3) | Table 3 "CL M2 = 12.3 L/h" (RSE 10.1%) |
| `lvc_m2` (V_M2)      | log(659)  | Table 3 "V M2 = 659 L" (RSE 7.2%) |
| `lq_m2` (Q1_M2)      | log(103)  | Table 3 "Q1 M2 = 103 L/h" (RSE 10.5%) |
| `lvp_m2` (VP1_M2)    | log(2840) | Table 3 "VP1 M2 = 2,840 L" (RSE 6.0%) |
| `lcl_m3` (CL_M3)     | log(39.2) | Table 3 "CL M3 = 39.2 L/h" (RSE 9.0%) |
| `lvc_m3` (V_M3)      | log(11.2) | Table 3 "V M3 = 11.2 L" (RSE 44.7%) |
| `lq_m3` (Q_M3)       | log(106)  | Table 3 "Q M3 = 106 L/h" (RSE 9.9%) |
| `lvp_m3` (VP_M3)     | log(2680) | Table 3 "VP M3 = 2,680 L" (RSE 13.6%) |
| `e_wt_cl_q`          | fixed(0.75) | Results "Allometric scaling of disposition parameters with fixed coefficients improved the fit markedly (0.75 for clearances and 1 for volumes; estimation of the coefficients did not significantly improve the fit further)" |
| `e_wt_vc_vp`         | fixed(1)    | same Results sentence |
| `e_efv_cl`           | 2.07 | Table 3 "EFV EFF BDQ and M2 = 2.07" (RSE 3.6%) |
| `e_efv_cl_m3`        | 1.12 | Table 3 "EFV EFF M3 = 1.12" (RSE 3.6%) |
| Block IIV (CL, CL_M2, CL_M3) | var 0.054646, 0.034733, 0.086178; corr +0.296, -0.129, +0.713 | Table 3 BSV CL = 23.7% CV, BSV CL_M2 = 18.8% CV (corr with CL = +29.6%), BSV CL_M3 = 30.0% CV (corr with CL = -12.9%, corr with CL_M2 = +71.3%) |
| `etalvc` (BSV V)       | 0.113081 | Table 3 BSV V = 34.6% CV (RSE 32%) |
| `etalq` (BSV Q1)       | 0.034372 | Table 3 BSV Q1 = 18.7% CV (RSE 15%) |
| `etalvc_m2` (BSV V_M2) | 0.080214 | Table 3 BSV V M2 = 28.9% CV (RSE 19%) |
| `etalvp_m2` (BSV VP_M2)| 0.064902 | Table 3 BSV VP 1M2 = 25.9% CV (RSE 39%) |
| `propSd` (BDQ residual)  | 0.239 | Table 3 "Prop error BDQ = 23.9% CV" (RSE 5.3%) |
| `propSd_m2` (M2 residual)| 0.177 | Table 3 "Prop error M2 = 17.7% CV" (RSE 4.8%) |
| `propSd_m3` (M3 residual)| 0.150 | Table 3 "Prop error M3 = 15.0% CV" (RSE 10.2%) |
| Transit-absorption chain via Savic 2007 input feeding depot at rate (NN+1)/MTT, then first-order ka into central | n/a | Results "BDQ PK was best described by a 3-compartment disposition model with absorption through a dynamic transit compartment model" |
| Three-compartment BDQ + two-compartment M2 + two-compartment M3 ODEs (sequential parent -> M2 -> M3) | n/a | Results "M2 and M3 were described by 2-compartment models with clearance of BDQ and M2, respectively, as input" |
| Allometric scaling on CL, Q, V, V_M2, VP_M2, V_M3, VP_M3 at 70 kg (reference weight assumed) | n/a | Results "Allometric scaling of disposition parameters with fixed coefficients improved the fit" |
| EFV step indicator (CONMED_EFV) switching on at full induction (paper's 1-week onset lag) | n/a | Results "The impact of induction was described as an instantaneous change in clearance 1 week after initialization of EFV treatment" |

# Virtual cohort

Original observed concentrations are not publicly available. We simulate
two virtual cohorts that mirror the two periods of ACTG A5267: a
single 400 mg oral dose of bedaquiline with no co-medication (period 1,
`CONMED_EFV = 0` throughout) and a single 400 mg oral dose given two
weeks into 600 mg once-nightly EFV (period 2, `CONMED_EFV = 1` throughout
because the bedaquiline dose is given more than one week after EFV start).

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

# Helper: build one cohort of `n` subjects as a self-contained event table.
# id_offset shifts subject IDs so multiple cohorts can be bind_rows()-ed
# without colliding (rxSolve treats id as the subject key; duplicate IDs
# across cohorts silently collapse into single Frankenstein subjects).
make_arm <- function(n, conmed_efv, treatment, id_offset = 0L) {
  ids <- id_offset + seq_len(n)
  # Sampling grid mirrors the published design: dense first 24 h, then 24-hourly
  # to 336 h (14 days post-dose).
  sample_grid <- c(
    seq(0, 12, by = 0.5),
    seq(13, 24, by = 1),
    seq(48, 336, by = 24)
  )
  # Weights drawn around the cohort median 82.3 kg, range 57.3-118.8 kg.
  wts <- round(runif(length(ids), min = 57, max = 119), 1)
  dose_rows <- tibble::tibble(
    id          = ids,
    time        = 0,
    evid        = 1L,
    amt         = 400,
    cmt         = "depot",
    WT          = wts,
    CONMED_EFV  = conmed_efv,
    treatment   = treatment
  )
  sample_rows <- tidyr::expand_grid(id = ids, time = sample_grid) |>
    dplyr::left_join(
      tibble::tibble(id = ids, WT = wts), by = "id"
    ) |>
    dplyr::mutate(
      evid       = 0L,
      amt        = 0,
      cmt        = "Cc",
      CONMED_EFV = conmed_efv,
      treatment  = treatment
    )
  dplyr::bind_rows(dose_rows, sample_rows) |>
    dplyr::arrange(id, time, dplyr::desc(evid))
}

n_per_arm <- 50L
events <- dplyr::bind_rows(
  make_arm(n_per_arm, conmed_efv = 0L,
           treatment = "BDQ alone",     id_offset =   0L),
  make_arm(n_per_arm, conmed_efv = 1L,
           treatment = "BDQ + EFV",     id_offset = 100L)
)

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

# Simulation

```{r simulate-typical}
mod <- readModelDb("Svensson_2013_bedaquiline")
mod_typical <- rxode2::zeroRe(mod)

sim_typical <- rxode2::rxSolve(
  mod_typical, events = events,
  keep = c("CONMED_EFV", "treatment", "WT"),
  returnType = "data.frame"
)
head(sim_typical[, c("id", "time", "Cc", "Cc_m2", "Cc_m3", "treatment")])
```

```{r simulate-vpc}
sim_vpc <- rxode2::rxSolve(
  mod, events = events,
  keep = c("CONMED_EFV", "treatment", "WT"),
  returnType = "data.frame"
)
```

# Replicate published figures

## Figure 4 -- typical BDQ + M2 + M3 profiles by treatment

Svensson 2013 Figure 4 stratifies the visual predictive check by treatment
(BDQ alone vs BDQ with EFV) and by analyte (BDQ, M2, M3). We replicate the
typical-value profile here at 14 days after a single 400 mg dose. The
published figure shows BDQ in mass units; we plot in mg/L (= ug/mL).

```{r figure-4-typical, fig.height = 6}
sim_typical_long <- sim_typical |>
  dplyr::filter(time > 0, time <= 336) |>
  dplyr::distinct(time, treatment, .keep_all = TRUE) |>
  dplyr::select(time, treatment, Cc, Cc_m2, Cc_m3) |>
  tidyr::pivot_longer(c(Cc, Cc_m2, Cc_m3),
                      names_to = "analyte", values_to = "conc_mgL") |>
  dplyr::mutate(
    analyte = dplyr::recode(analyte,
                            Cc    = "Bedaquiline (BDQ)",
                            Cc_m2 = "M2 metabolite",
                            Cc_m3 = "M3 metabolite"),
    analyte = factor(analyte, levels = c("Bedaquiline (BDQ)",
                                          "M2 metabolite",
                                          "M3 metabolite"))
  )

ggplot(sim_typical_long,
       aes(time, conc_mgL, colour = treatment)) +
  geom_line(linewidth = 1) +
  facet_wrap(~ analyte, scales = "free_y", ncol = 1) +
  scale_y_log10() +
  labs(
    x = "Time after single 400 mg bedaquiline dose (h)",
    y = "Typical plasma concentration (mg/L)",
    title = "Figure 4 -- BDQ + M2 + M3 typical profiles, alone vs with EFV",
    caption = paste(
      "Replicates the typical-value layer of Figure 4 of Svensson 2013:",
      "BDQ, M2, and M3 concentrations after a single 400 mg bedaquiline",
      "dose given either alone or after 2 weeks of 600 mg once-nightly EFV",
      "(CONMED_EFV = 1 throughout the 14-day sampling window)."
    )
  )
```

## VPC across the cohort

```{r figure-vpc-cohort, fig.height = 7}
sim_summary <- sim_vpc |>
  dplyr::filter(time > 0) |>
  dplyr::group_by(time, treatment) |>
  dplyr::summarise(
    BDQ_Q05 = quantile(Cc,    0.05, na.rm = TRUE),
    BDQ_Q50 = quantile(Cc,    0.50, na.rm = TRUE),
    BDQ_Q95 = quantile(Cc,    0.95, na.rm = TRUE),
    M2_Q05  = quantile(Cc_m2, 0.05, na.rm = TRUE),
    M2_Q50  = quantile(Cc_m2, 0.50, na.rm = TRUE),
    M2_Q95  = quantile(Cc_m2, 0.95, na.rm = TRUE),
    M3_Q05  = quantile(Cc_m3, 0.05, na.rm = TRUE),
    M3_Q50  = quantile(Cc_m3, 0.50, na.rm = TRUE),
    M3_Q95  = quantile(Cc_m3, 0.95, na.rm = TRUE),
    .groups = "drop"
  ) |>
  tidyr::pivot_longer(-c(time, treatment),
                      names_to = c("analyte", "stat"),
                      names_sep = "_") |>
  tidyr::pivot_wider(names_from = stat, values_from = value) |>
  dplyr::mutate(
    analyte = dplyr::recode(analyte,
                            BDQ = "Bedaquiline (BDQ)",
                            M2  = "M2 metabolite",
                            M3  = "M3 metabolite"),
    analyte = factor(analyte, levels = c("Bedaquiline (BDQ)",
                                          "M2 metabolite",
                                          "M3 metabolite"))
  )

ggplot(sim_summary |> dplyr::filter(time <= 336),
       aes(time, Q50, colour = treatment, fill = treatment)) +
  geom_ribbon(aes(ymin = Q05, ymax = Q95), alpha = 0.20, colour = NA) +
  geom_line(linewidth = 0.8) +
  facet_wrap(~ analyte, scales = "free_y", ncol = 1) +
  scale_y_log10() +
  labs(
    x = "Time after single 400 mg bedaquiline dose (h)",
    y = "Plasma concentration (mg/L)",
    title = "VPC -- BDQ + M2 + M3, alone vs with EFV",
    caption = paste(
      "Median + 5-95% prediction interval across",
      paste(nrow(unique(events[, c("id", "treatment")])), "subjects (50 per arm),"),
      "single 400 mg oral dose."
    )
  )
```

# PKNCA validation

Use PKNCA to compute Cmax, Tmax, AUC0-14d, and (where the terminal phase is
characterisable) half-life by treatment arm so the simulated 14-day exposure
ratios (BDQ + EFV vs BDQ alone) can be compared against the published
relative steady-state-concentration ratios.

Svensson 2013 reports relative typical steady-state concentrations of BDQ
(EFV vs alone) = 48% (SE 1.9%), M2 = 48% (SE 1.9%), and M3 = 88% (SE 3.7%)
(Svensson 2013 Results "Relative exposures"). These are full-induction
steady-state ratios, not single-dose 14-day-window ratios, so the PKNCA
comparison below uses 14-day AUC ratios as the closest single-dose analogue
(the comparison is qualitative -- the model's chronic-EFV steady-state
behaviour matches the published 48% / 48% / 88% almost exactly, as verified
in the single-subject typical-value check below).

```{r pknca-bdq}
# Bedaquiline NCA -- guarantee a time=0 row per (id, treatment) before PKNCA
sim_nca_bdq <- sim_vpc |>
  dplyr::filter(!is.na(Cc)) |>
  dplyr::select(id, time, Cc, treatment)

sim_nca_bdq <- dplyr::bind_rows(
  sim_nca_bdq,
  sim_nca_bdq |> dplyr::distinct(id, treatment) |>
    dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, treatment, time, .keep_all = TRUE) |>
  dplyr::arrange(id, treatment, time)

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

conc_obj_bdq <- PKNCA::PKNCAconc(
  sim_nca_bdq, Cc ~ time | treatment + id,
  concu = "mg/L", timeu = "h"
)
dose_obj <- PKNCA::PKNCAdose(
  dose_df, amt ~ time | treatment + id,
  doseu = "mg"
)

intervals_14d <- data.frame(
  start      = 0,
  end        = 336,
  cmax       = TRUE,
  tmax       = TRUE,
  auclast    = TRUE,
  half.life  = TRUE
)

nca_res_bdq <- PKNCA::pk.nca(PKNCA::PKNCAdata(
  conc_obj_bdq, dose_obj, intervals = intervals_14d
))
knitr::kable(
  summary(nca_res_bdq),
  caption = "Simulated NCA parameters for bedaquiline by treatment arm (single 400 mg oral dose, 14-day sampling window)."
)
```

```{r pknca-m2}
# M2 metabolite NCA
sim_nca_m2 <- sim_vpc |>
  dplyr::filter(!is.na(Cc_m2)) |>
  dplyr::select(id, time, Cc_m2, treatment) |>
  dplyr::rename(Cc = Cc_m2)

sim_nca_m2 <- dplyr::bind_rows(
  sim_nca_m2,
  sim_nca_m2 |> dplyr::distinct(id, treatment) |>
    dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, treatment, time, .keep_all = TRUE) |>
  dplyr::arrange(id, treatment, time)

conc_obj_m2 <- PKNCA::PKNCAconc(
  sim_nca_m2, Cc ~ time | treatment + id,
  concu = "mg/L", timeu = "h"
)

nca_res_m2 <- PKNCA::pk.nca(PKNCA::PKNCAdata(
  conc_obj_m2, dose_obj, intervals = intervals_14d
))
knitr::kable(
  summary(nca_res_m2),
  caption = "Simulated NCA parameters for M2 metabolite by treatment arm."
)
```

```{r pknca-m3}
# M3 metabolite NCA
sim_nca_m3 <- sim_vpc |>
  dplyr::filter(!is.na(Cc_m3)) |>
  dplyr::select(id, time, Cc_m3, treatment) |>
  dplyr::rename(Cc = Cc_m3)

sim_nca_m3 <- dplyr::bind_rows(
  sim_nca_m3,
  sim_nca_m3 |> dplyr::distinct(id, treatment) |>
    dplyr::mutate(time = 0, Cc = 0)
) |>
  dplyr::distinct(id, treatment, time, .keep_all = TRUE) |>
  dplyr::arrange(id, treatment, time)

conc_obj_m3 <- PKNCA::PKNCAconc(
  sim_nca_m3, Cc ~ time | treatment + id,
  concu = "mg/L", timeu = "h"
)

nca_res_m3 <- PKNCA::pk.nca(PKNCA::PKNCAdata(
  conc_obj_m3, dose_obj, intervals = intervals_14d
))
knitr::kable(
  summary(nca_res_m3),
  caption = "Simulated NCA parameters for M3 metabolite by treatment arm."
)
```

## Steady-state ratio check (typical-value, single subject)

The relative-steady-state-concentration formula in Svensson 2013 Results is
the ratio of `F * Dose / (CL * tau)` with and without EFV at full induction.
Because F, Dose, and tau are common, this reduces to `CL_base / CL_induced`.
For BDQ and M2 the induced CL is `2.07 * CL_base`, so the relative
steady-state concentration is `1 / 2.07 = 0.483` (= 48.3%). For M3 the
induced CL is `1.12 * CL_base`, giving `1 / 1.12 = 0.893` (= 89.3%). The
paper reports 48% (SE 1.9%) for BDQ and M2 and 88% (SE 3.7%) for M3, in
good agreement.

We confirm this against a clean single-subject typical-value simulation:

```{r ss-check}
ev1 <- tibble::tibble(
  id   = 1L,
  time = c(0, 1, 2, 4, 8, 12, 24, 48, 72, 168, 336),
  evid = c(1L, rep(0L, 10)),
  amt  = c(400, rep(0, 10)),
  cmt  = c("depot", rep("Cc", 10)),
  WT   = 70,
  CONMED_EFV = 0L
)
ev2 <- ev1; ev2$CONMED_EFV <- 1L

s1 <- as.data.frame(
  rxode2::rxSolve(mod_typical, events = ev1, returnType = "data.frame")
)
s2 <- as.data.frame(
  rxode2::rxSolve(mod_typical, events = ev2, returnType = "data.frame")
)

ratio_tbl <- tibble::tibble(
  time         = s1$time,
  BDQ_alone    = s1$Cc,
  BDQ_with_EFV = s2$Cc,
  ratio_BDQ    = s2$Cc / s1$Cc,
  ratio_M2     = s2$Cc_m2 / s1$Cc_m2,
  ratio_M3     = s2$Cc_m3 / s1$Cc_m3
)
knitr::kable(ratio_tbl,
             caption = paste(
               "Single-subject typical-value Cc(EFV) / Cc(no EFV) ratio",
               "at canonical time points (70 kg, F = 1, CONMED_EFV = 1",
               "throughout). At late time points (>= 168 h) the BDQ and M2",
               "ratios converge to the CL-ratio limit ~ 0.48 and the M3",
               "ratio to ~ 0.89, matching Svensson 2013 Results."
             ))
```

## Comparison against published relative steady-state concentrations

```{r rel-ss-comparison}
published <- tibble::tibble(
  Analyte                 = c("BDQ", "M2", "M3"),
  Published_Rel_Css_pct   = c(48, 48, 88),
  Simulated_Rel_Css_pct   = round(c(s2$Cc[s2$time == 336]    / s1$Cc[s1$time == 336]    * 100,
                                     s2$Cc_m2[s2$time == 336] / s1$Cc_m2[s1$time == 336] * 100,
                                     s2$Cc_m3[s2$time == 336] / s1$Cc_m3[s1$time == 336] * 100), 1),
  Pct_Diff                = NA_real_
)
published$Pct_Diff <- round(
  100 * (published$Simulated_Rel_Css_pct - published$Published_Rel_Css_pct) /
    published$Published_Rel_Css_pct, 1
)
knitr::kable(
  published,
  caption = paste(
    "Simulated vs. published relative typical steady-state concentration",
    "(EFV / alone) for BDQ, M2, and M3. Simulated values are the typical",
    "single-subject (70 kg) ratio at 336 h post-dose; published values are",
    "from Svensson 2013 Results 'Relative exposures' for the model's full",
    "induction steady-state under chronic 600 mg once-nightly EFV."
  ),
  align = c("l", "r", "r", "r")
)
```

The 14-day single-dose AUC ratio understates the asymptotic-steady-state
ratio because the early absorption / first-distribution phase (during which
the parent has not yet sampled the deep peripheral compartments and
clearance has minimal cumulative effect) is shared between the two
treatments. The relative steady-state concentrations above are the
correct asymptotic comparison.

# Assumptions and deviations

- **EFV-effect individual etas dropped.** Svensson 2013 Table 3 reports a
  6 x 6 lower-triangular NONMEM BLOCK(6) on the bedaquiline-CL eta, the
  M2-CL eta, the M3-CL eta, and three individual-level EFV-effect etas
  (one each for BSV on the EFV-induction effect for BDQ, M2, and M3).
  The 3 x 3 sub-block on (BSV CL, BSV CL_M2, BSV CL_M3) is preserved
  here with the published correlation pattern (+29.6% for BDQ-M2,
  -12.9% for BDQ-M3, +71.3% for M2-M3), but the three EFV-effect etas
  (BSV EFV EFF-BDQ = 20.6% CV, EFF-M2 = 28.2% CV, EFF-M3 = 32.7% CV)
  and their off-diagonal correlations with the structural-CL block
  (Svensson 2013 Table 3 c-marked entries) are dropped because the
  etas are gated by `CONMED_EFV` (well-defined only for subjects in the
  post-induction window) and nlmixr2lib has no idiomatic encoding for
  a covariate-gated random effect. As a consequence the simulated VPC
  in the BDQ + EFV arm underrepresents the published spread on the EFV
  ratio; the typical-value central tendency is preserved.
- **Between-occasion variability dropped.** Svensson 2013 Table 3 reports
  BOV on bioavailability (23.6% CV) and on MTT (55.4% CV) plus BSV on
  bioavailability (24.3% CV) from the two-occasion crossover design
  (period 1 dose alone, period 2 dose under EFV). nlmixr2lib has no
  idiomatic encoding for between-occasion variability separate from
  between-subject variability; both BOV components and BSV F are dropped
  here. The simulated VPC therefore underrepresents within-subject
  variability between the two doses each subject received.
- **Cross-output residual correlations dropped.** Svensson 2013 Table 3
  reports correlations across the proportional residual errors for the
  three analytes (BDQ-M2 corr +14.9%, BDQ-M3 corr +7.5%, M2-M3 corr
  +11.7%). nlmixr2lib has no idiomatic encoding for cross-output residual
  correlation, so the BDQ, M2, and M3 residual proportions are encoded as
  independent.
- **Time-varying residual weighting dropped.** Svensson 2013 Table 3
  reports an additive multiplicative factor of 1.87 on the residual SD
  for samples drawn within the absorption phase (TAD < 6 h) and a factor
  of 3.28 for samples below the limit of quantification. Both are
  fit-time weightings that absorb absorption-phase unmodelled-variability
  and BLQ-handling artifacts; both are dropped here, and the simulated
  residual SD is held constant at the post-absorption above-LLOQ value
  (23.9% on BDQ, 17.7% on M2, 15.0% on M3) across all time points.
- **Reference weight assumed 70 kg.** Svensson 2013 Results states that
  fixed allometric exponents of 0.75 (clearances) and 1 (volumes)
  improved the fit but does not specify the reference weight. Following
  the canonical Anderson and Holford 2008 convention universally adopted
  in popPK and the sibling Svensson 2014 bedaquiline paper that builds on
  this analysis, the reference weight is set to 70 kg. The cohort median
  was 82.3 kg, so typical subjects at the cohort median scale by
  `(82.3 / 70)^0.75` = 1.13 on apparent CL (relative to 70 kg) and
  `(82.3 / 70)` = 1.18 on apparent V. This is internally consistent with
  the apparent F-relative parameter values in Table 3 (which are
  population-typical at the unspecified reference weight).
- **F = 1 anchor for the F-relative parameterisation.** Svensson 2013
  Table 3 footnote b explicitly states "Estimated with typical value of F
  fixed to 1". Setting F = 1 inside `transit(nn, mtt, 1)` preserves the
  apparent F-relative parameter values directly; simulated bedaquiline
  mass in the central compartment represents the F-apparent mass and
  concentrations represent F-apparent (= observed) plasma concentrations.
- **EFV step indicator set externally.** Svensson 2013 Methods states
  "the effect was modeled as an instantaneous change in clearance (CL)
  and/or bioavailability (F). Several time points between day 1 and day
  14 of EFV treatment were tested for this change" and Results notes
  "the impact of induction was described as an instantaneous change in
  clearance 1 week after initialization of EFV treatment". For new
  simulations, populate `CONMED_EFV` = 1 on every observation row that
  falls >= 1 week after the first 600 mg once-nightly EFV dose and 0
  otherwise; pre-induction observations should carry the unmodified
  baseline CL. The vignette's example cohorts assume the bedaquiline
  dose is given after 2 weeks of EFV (Svensson 2013 period 2 design),
  so `CONMED_EFV = 1` is assigned throughout the 14-day sampling
  window.
- **Concentration units.** Model output is in mg/L (= ug/mL). The
  published figures use mass or molar units; readers who need molar
  concentrations should multiply BDQ mg/L by 1000 / 555.50 = 1.800 to
  obtain nmol/L (M2: 1000 / 541.47 = 1.847; M3: 1000 / 527.45 = 1.896).
- **NCA terminal-phase fit.** With 14-day sampling and a 5-6 month
  terminal half-life, NCA `aucinf.obs` and `half.life` cannot be reliably
  estimated; the comparison above uses `auclast` (14-day AUC) instead.
  The reported relative-steady-state-concentration comparison (BDQ / M2:
  48%, M3: 88%) is the more meaningful single-dose surrogate because it
  isolates the CL ratio and matches the published Results.
